ROOT logo
/*
   UserTreeAnalysis implementation

   Details regarding the working of the analysis itself are given below,
   just after the function header.
*/

//#include "UserTreeAnalysis.H"   // remove if copied to user working directory
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "MC4Vector.H"
#include "HEPParticle.H"
#include "TH1.h"
#include "Setup.H"
#include "TObjArray.h"
#include "TMath.h"
#define PI 3.141592653

inline bool ifSofty(int Id,int nparams, double *params){
  if (nparams<5 && Id==22) return true;   // to remove photons only
  for (int i=nparams-1; i>3; i--)
    if (Id==params[i]) return true;       // to remove all what is in params from nparams down to 4
  return false;
}

// very similar to  MC_FillUserHistogram from Generate.cxx
inline void fillUserHisto(char *name,double val, double weight=1.0, 
                          double min=Setup::bin_min[0][0], 
                          double max=Setup::bin_max[0][0]){

    TH1D *h=(TH1D*)(Setup::user_histograms->FindObject(name));
    if(!h){
      h=new TH1D(name,name,Setup::nbins[0][0],min,max);
      if(!h) return;
      Setup::user_histograms->Add(h);
      //      printf("user histogram created %s\n", name);
}
    h->Fill(val,weight);

}
double angle(double X, double Y){

 
    double  an=0.0; 
    double R=sqrt(X*X+Y*Y); 
    //    if(R<pow(10,-20)) printf(" angle this time X %f\n", 10000*X);
    if(R<pow(10,-20)) return an;

    if(TMath::Abs(X)/R<0.8) 
    { 
        an=acos(X/R);
	if(Y<0 && an>0) an=-an;
        if(Y>0 && an<0) an=-an; 
    }
    else 
    { 
      an=asin(Y/R);
      if(X<0 && an>=0.)  an=PI-an; 
      else if(X<0.)      an=-PI-an; 
       
    } 
  return an;
}

int ZmumuAnalysis(HEPParticle *mother,HEPParticleList *stableDaughters, int nparams, double *params)
{
  // THIS IS EXAMPLE of userTreeAnalysis. It acts on list of particle X (under our MC-test) final
  // daughters. Of course in real life many options may need to be introduced by the user.
    // PARAMETERS:
    // params[0] threshold on Energy (or p_T) of particle expressed as a fraction of mother's 
    //		 Energy (or p_T) in mother's frame. If not specified - default is 0.05 
    //           Note that setting it to zero is possible.
    // params[1] maximum number of left soft-suspected particles. 
    //           0: (default) all listed particles are removed, even if hard
    //           1: 2:  one two etc  removable particles are kept (at most)
    //          -1: this option is off, all hard particles stay.
    // params[2] control tag on discriminating particles, 
    //           0: (default) Energy in decaying particle rest frame
    //           1: Energy in lab frame
    //           2: p_T with respect to beam direction in lab frame.
    // params[3] type of action to be applied on removed daughters 
    //           0: (default) nothing is done, they are lost
    //           1:  algorithm as in studies on PHOTOS is used, see papers P. Golonka and Z. Was.
    // params[4] from this paramter on PDG id-s of particles to be removed can be added. 
    //           Default is that only photons are removed.
    // To get to origin of our particle (mother) one need to go after UserEventAnalysis,
    // instructive example will be given later. 
    assert(mother!=0);
    assert(stableDaughters!=0);


    double threshold=0.05, leftmax=0.0, selector=0.0, actLost=0.0; 
    if (nparams>0 && params==0)  return -1;
    if (nparams>0) threshold=params[0]; 
    if (nparams>1) leftmax  =params[1]; 
    if (nparams>2) selector =params[2]; 
    if (nparams>3) actLost  =params[3]; 

    
    HEPParticleList *lostDaughters=new HEPParticleList;    
    double pt_limit=threshold*mother->GetM();

    HEPParticleListIterator daughters (*stableDaughters);
    int nphot=0;
    double ephot=pow(10,22);
    bool redo=false;
    for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
      if(redo) redo=false;
	MC4Vector d4(part->GetE(),part->GetPx(),part->GetPy(),part->GetPz(),part->GetM());
	// boost to mother's frame:
	d4.Boost(mother->GetPx(),mother->GetPy(),mother->GetPz(),mother->GetE(),mother->GetM());



	double p_t=d4.GetX0();  // default
	switch ((int) selector)
        {
	case 1: p_t=part->GetE(); break;
	case 2: p_t=d4.Xt(); break;
	}


        if(ifSofty(part->GetPDGId(),nparams,params)) nphot++;
	if ( ifSofty(part->GetPDGId(),nparams,params) && p_t < pt_limit) {
	  //	  printf("Odrzucamy! %i\n",count++);
          nphot=0;
            lostDaughters->push_back(part);
	    stableDaughters->remove(part);
	    part=daughters.first();
            redo=true;
            continue;
	}
        if( ifSofty(part->GetPDGId(),nparams,params) && ephot>p_t) ephot=p_t;
    }
    while(nphot>(int) leftmax)
    {
      double ephot1=pow(10,22);
      redo=false;
      for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
        if(redo) redo=false;
	MC4Vector d4(part->GetE(),part->GetPx(),part->GetPy(),part->GetPz(),part->GetM());
	// boost to mother's frame:
	d4.Boost(mother->GetPx(),mother->GetPy(),mother->GetPz(),mother->GetE(),mother->GetM());


	double p_t=d4.GetX0();  // default
	switch ((int) selector)
        {
	case 1: p_t=part->GetE(); break;
	case 2: p_t=d4.Xt(); break;
	}


	if ( ifSofty(part->GetPDGId(),nparams,params) && p_t == ephot) {

	  nphot--;
            lostDaughters->push_back(part);
	    stableDaughters->remove(part);
            ephot1=pow(10,22);
	     part=daughters.first();
             redo=true;
             continue;
	}
        if( ifSofty(part->GetPDGId(),nparams,params) && ephot1>p_t) ephot1=p_t;
      }
      ephot=ephot1;

    }

    // ##############################################################
    // SPECIAL CASE - ordered or symmetrical photons
    //                This code changes order of photons
    // ##############################################################
    HEPParticleListIterator symmetry (*stableDaughters);
    HEPParticle *phot1 = NULL;
    HEPParticle *phot2 = NULL;
    for (HEPParticle *part=symmetry.first(); part!=0; part=symmetry.next() )
    {
      if(part->GetPDGId()==22)
      {
        if(!phot1) phot1=part;
        else
        {
          phot2=part;
          break;
        }
      }
    }
    if(phot1 && phot2)
    {
      if(phot1->GetE()<phot2->GetE()) // for ordered photons
      //if( rand()*1.0/RAND_MAX > 0.5) // for photon symmetrization
      {
        double buf_x = phot1->GetPx();
        double buf_y = phot1->GetPy();
        double buf_z = phot1->GetPz();
        double buf_e = phot1->GetE();

        phot1->SetPx(phot2->GetPx());
        phot1->SetPy(phot2->GetPy());
        phot1->SetPz(phot2->GetPz());
        phot1->SetE (phot2->GetE() );

        phot2->SetPx(buf_x);
        phot2->SetPy(buf_y);
        phot2->SetPz(buf_z);
        phot2->SetE (buf_e);
      }
    }
 
    // ##############################################################
    // here functionality of removig some daughters is finished
    // now we reconsider what to do with them
    // ##############################################################
    // delete lostDaughters;
    // return 1; 


    // ##############################################################
    // Now: What to do with lost daughters?
    // ##############################################################
    //    lostDaughters->ls();
    int version=(int) actLost;
        
    switch(version)
      {
      case 1:  // add lost to charged
	{
          HEPParticleListIterator lost (*lostDaughters);
	  for (HEPParticle *lostpart=lost.first(); lostpart!=0; lostpart=lost.next() ) {
            HEPParticle* Best=0;
	    double searchvirt=pow(10.0,22);
            MC4Vector VV;
	    for( HEPParticle *part=daughters.first(); part!=0; part=daughters.next() ){
	      if(part->GetCharge()==0) continue;
              VV=lostpart->GetP4()+part->GetP4();
	      VV.AdjustM();
	      if (VV.GetM()<searchvirt) {searchvirt=VV.GetM(); Best=part;}
	    }
	    if(Best) {
              Best->SetPx(Best->GetPx()+lostpart->GetPx());
              Best->SetPy(Best->GetPy()+lostpart->GetPy());
              Best->SetPz(Best->GetPz()+lostpart->GetPz());
              Best->SetE (Best->GetE ()+lostpart->GetE ());
	    }
	  }
	break;
	}
      default: break; // do nothing
      }

    delete lostDaughters;

    
    bool activateUserHist=true;
    if(!activateUserHist) return 1;

    // segmet of user defined histograms

    double X=mother->GetPx(), Y=mother->GetPy(), Z=mother->GetPz(); 
    // double E=mother->GetE(),  MM=mother->GetM();


    double pt=sqrt(X*X+Y*Y);
    double eta=log((sqrt(pt*pt+Z*Z)+TMath::Abs(Z))/pt);
    if(Z<0 && eta>0) eta=-eta;
    if(Z>0 && eta<0) eta=-eta;
    double phi=angle(X,Y);
    char hist1[] = "mother-PT";
    char hist2[] = "mother-eta";
    char hist3[] = "mother-phi";
    fillUserHisto(hist1,pt,1.0,0.0,100.0);
    fillUserHisto(hist2,eta,1.0,-8.0,8.0);
    fillUserHisto(hist3,phi,1.0,-PI,PI);
    return 1;
}

 ZmumuAnalysis.C:1
 ZmumuAnalysis.C:2
 ZmumuAnalysis.C:3
 ZmumuAnalysis.C:4
 ZmumuAnalysis.C:5
 ZmumuAnalysis.C:6
 ZmumuAnalysis.C:7
 ZmumuAnalysis.C:8
 ZmumuAnalysis.C:9
 ZmumuAnalysis.C:10
 ZmumuAnalysis.C:11
 ZmumuAnalysis.C:12
 ZmumuAnalysis.C:13
 ZmumuAnalysis.C:14
 ZmumuAnalysis.C:15
 ZmumuAnalysis.C:16
 ZmumuAnalysis.C:17
 ZmumuAnalysis.C:18
 ZmumuAnalysis.C:19
 ZmumuAnalysis.C:20
 ZmumuAnalysis.C:21
 ZmumuAnalysis.C:22
 ZmumuAnalysis.C:23
 ZmumuAnalysis.C:24
 ZmumuAnalysis.C:25
 ZmumuAnalysis.C:26
 ZmumuAnalysis.C:27
 ZmumuAnalysis.C:28
 ZmumuAnalysis.C:29
 ZmumuAnalysis.C:30
 ZmumuAnalysis.C:31
 ZmumuAnalysis.C:32
 ZmumuAnalysis.C:33
 ZmumuAnalysis.C:34
 ZmumuAnalysis.C:35
 ZmumuAnalysis.C:36
 ZmumuAnalysis.C:37
 ZmumuAnalysis.C:38
 ZmumuAnalysis.C:39
 ZmumuAnalysis.C:40
 ZmumuAnalysis.C:41
 ZmumuAnalysis.C:42
 ZmumuAnalysis.C:43
 ZmumuAnalysis.C:44
 ZmumuAnalysis.C:45
 ZmumuAnalysis.C:46
 ZmumuAnalysis.C:47
 ZmumuAnalysis.C:48
 ZmumuAnalysis.C:49
 ZmumuAnalysis.C:50
 ZmumuAnalysis.C:51
 ZmumuAnalysis.C:52
 ZmumuAnalysis.C:53
 ZmumuAnalysis.C:54
 ZmumuAnalysis.C:55
 ZmumuAnalysis.C:56
 ZmumuAnalysis.C:57
 ZmumuAnalysis.C:58
 ZmumuAnalysis.C:59
 ZmumuAnalysis.C:60
 ZmumuAnalysis.C:61
 ZmumuAnalysis.C:62
 ZmumuAnalysis.C:63
 ZmumuAnalysis.C:64
 ZmumuAnalysis.C:65
 ZmumuAnalysis.C:66
 ZmumuAnalysis.C:67
 ZmumuAnalysis.C:68
 ZmumuAnalysis.C:69
 ZmumuAnalysis.C:70
 ZmumuAnalysis.C:71
 ZmumuAnalysis.C:72
 ZmumuAnalysis.C:73
 ZmumuAnalysis.C:74
 ZmumuAnalysis.C:75
 ZmumuAnalysis.C:76
 ZmumuAnalysis.C:77
 ZmumuAnalysis.C:78
 ZmumuAnalysis.C:79
 ZmumuAnalysis.C:80
 ZmumuAnalysis.C:81
 ZmumuAnalysis.C:82
 ZmumuAnalysis.C:83
 ZmumuAnalysis.C:84
 ZmumuAnalysis.C:85
 ZmumuAnalysis.C:86
 ZmumuAnalysis.C:87
 ZmumuAnalysis.C:88
 ZmumuAnalysis.C:89
 ZmumuAnalysis.C:90
 ZmumuAnalysis.C:91
 ZmumuAnalysis.C:92
 ZmumuAnalysis.C:93
 ZmumuAnalysis.C:94
 ZmumuAnalysis.C:95
 ZmumuAnalysis.C:96
 ZmumuAnalysis.C:97
 ZmumuAnalysis.C:98
 ZmumuAnalysis.C:99
 ZmumuAnalysis.C:100
 ZmumuAnalysis.C:101
 ZmumuAnalysis.C:102
 ZmumuAnalysis.C:103
 ZmumuAnalysis.C:104
 ZmumuAnalysis.C:105
 ZmumuAnalysis.C:106
 ZmumuAnalysis.C:107
 ZmumuAnalysis.C:108
 ZmumuAnalysis.C:109
 ZmumuAnalysis.C:110
 ZmumuAnalysis.C:111
 ZmumuAnalysis.C:112
 ZmumuAnalysis.C:113
 ZmumuAnalysis.C:114
 ZmumuAnalysis.C:115
 ZmumuAnalysis.C:116
 ZmumuAnalysis.C:117
 ZmumuAnalysis.C:118
 ZmumuAnalysis.C:119
 ZmumuAnalysis.C:120
 ZmumuAnalysis.C:121
 ZmumuAnalysis.C:122
 ZmumuAnalysis.C:123
 ZmumuAnalysis.C:124
 ZmumuAnalysis.C:125
 ZmumuAnalysis.C:126
 ZmumuAnalysis.C:127
 ZmumuAnalysis.C:128
 ZmumuAnalysis.C:129
 ZmumuAnalysis.C:130
 ZmumuAnalysis.C:131
 ZmumuAnalysis.C:132
 ZmumuAnalysis.C:133
 ZmumuAnalysis.C:134
 ZmumuAnalysis.C:135
 ZmumuAnalysis.C:136
 ZmumuAnalysis.C:137
 ZmumuAnalysis.C:138
 ZmumuAnalysis.C:139
 ZmumuAnalysis.C:140
 ZmumuAnalysis.C:141
 ZmumuAnalysis.C:142
 ZmumuAnalysis.C:143
 ZmumuAnalysis.C:144
 ZmumuAnalysis.C:145
 ZmumuAnalysis.C:146
 ZmumuAnalysis.C:147
 ZmumuAnalysis.C:148
 ZmumuAnalysis.C:149
 ZmumuAnalysis.C:150
 ZmumuAnalysis.C:151
 ZmumuAnalysis.C:152
 ZmumuAnalysis.C:153
 ZmumuAnalysis.C:154
 ZmumuAnalysis.C:155
 ZmumuAnalysis.C:156
 ZmumuAnalysis.C:157
 ZmumuAnalysis.C:158
 ZmumuAnalysis.C:159
 ZmumuAnalysis.C:160
 ZmumuAnalysis.C:161
 ZmumuAnalysis.C:162
 ZmumuAnalysis.C:163
 ZmumuAnalysis.C:164
 ZmumuAnalysis.C:165
 ZmumuAnalysis.C:166
 ZmumuAnalysis.C:167
 ZmumuAnalysis.C:168
 ZmumuAnalysis.C:169
 ZmumuAnalysis.C:170
 ZmumuAnalysis.C:171
 ZmumuAnalysis.C:172
 ZmumuAnalysis.C:173
 ZmumuAnalysis.C:174
 ZmumuAnalysis.C:175
 ZmumuAnalysis.C:176
 ZmumuAnalysis.C:177
 ZmumuAnalysis.C:178
 ZmumuAnalysis.C:179
 ZmumuAnalysis.C:180
 ZmumuAnalysis.C:181
 ZmumuAnalysis.C:182
 ZmumuAnalysis.C:183
 ZmumuAnalysis.C:184
 ZmumuAnalysis.C:185
 ZmumuAnalysis.C:186
 ZmumuAnalysis.C:187
 ZmumuAnalysis.C:188
 ZmumuAnalysis.C:189
 ZmumuAnalysis.C:190
 ZmumuAnalysis.C:191
 ZmumuAnalysis.C:192
 ZmumuAnalysis.C:193
 ZmumuAnalysis.C:194
 ZmumuAnalysis.C:195
 ZmumuAnalysis.C:196
 ZmumuAnalysis.C:197
 ZmumuAnalysis.C:198
 ZmumuAnalysis.C:199
 ZmumuAnalysis.C:200
 ZmumuAnalysis.C:201
 ZmumuAnalysis.C:202
 ZmumuAnalysis.C:203
 ZmumuAnalysis.C:204
 ZmumuAnalysis.C:205
 ZmumuAnalysis.C:206
 ZmumuAnalysis.C:207
 ZmumuAnalysis.C:208
 ZmumuAnalysis.C:209
 ZmumuAnalysis.C:210
 ZmumuAnalysis.C:211
 ZmumuAnalysis.C:212
 ZmumuAnalysis.C:213
 ZmumuAnalysis.C:214
 ZmumuAnalysis.C:215
 ZmumuAnalysis.C:216
 ZmumuAnalysis.C:217
 ZmumuAnalysis.C:218
 ZmumuAnalysis.C:219
 ZmumuAnalysis.C:220
 ZmumuAnalysis.C:221
 ZmumuAnalysis.C:222
 ZmumuAnalysis.C:223
 ZmumuAnalysis.C:224
 ZmumuAnalysis.C:225
 ZmumuAnalysis.C:226
 ZmumuAnalysis.C:227
 ZmumuAnalysis.C:228
 ZmumuAnalysis.C:229
 ZmumuAnalysis.C:230
 ZmumuAnalysis.C:231
 ZmumuAnalysis.C:232
 ZmumuAnalysis.C:233
 ZmumuAnalysis.C:234
 ZmumuAnalysis.C:235
 ZmumuAnalysis.C:236
 ZmumuAnalysis.C:237
 ZmumuAnalysis.C:238
 ZmumuAnalysis.C:239
 ZmumuAnalysis.C:240
 ZmumuAnalysis.C:241
 ZmumuAnalysis.C:242
 ZmumuAnalysis.C:243
 ZmumuAnalysis.C:244
 ZmumuAnalysis.C:245
 ZmumuAnalysis.C:246
 ZmumuAnalysis.C:247
 ZmumuAnalysis.C:248
 ZmumuAnalysis.C:249
 ZmumuAnalysis.C:250
 ZmumuAnalysis.C:251
 ZmumuAnalysis.C:252
 ZmumuAnalysis.C:253
 ZmumuAnalysis.C:254
 ZmumuAnalysis.C:255
 ZmumuAnalysis.C:256
 ZmumuAnalysis.C:257
 ZmumuAnalysis.C:258
 ZmumuAnalysis.C:259
 ZmumuAnalysis.C:260
 ZmumuAnalysis.C:261
 ZmumuAnalysis.C:262
 ZmumuAnalysis.C:263
 ZmumuAnalysis.C:264
 ZmumuAnalysis.C:265
 ZmumuAnalysis.C:266
 ZmumuAnalysis.C:267
 ZmumuAnalysis.C:268
 ZmumuAnalysis.C:269
 ZmumuAnalysis.C:270
 ZmumuAnalysis.C:271
 ZmumuAnalysis.C:272
 ZmumuAnalysis.C:273
 ZmumuAnalysis.C:274
 ZmumuAnalysis.C:275
 ZmumuAnalysis.C:276
 ZmumuAnalysis.C:277
 ZmumuAnalysis.C:278
 ZmumuAnalysis.C:279