TEvtGen/Photos/examples/testing/ZmumuNLO/ZmumuAnalysis.C

   1 /*
   2    UserTreeAnalysis implementation
   3
   4    Details regarding the working of the analysis itself are given below,
   5    just after the function header.
   6 */
   7
   8 //#include "UserTreeAnalysis.H"   // remove if copied to user working directory
   9 #include <stdio.h>
  10 #include <assert.h>
  11 #include <math.h>
  12 #include "MC4Vector.H"
  13 #include "HEPParticle.H"
  14 #include "TH1.h"
  15 #include "Setup.H"
  16 #include "TObjArray.h"
  17 #include "TMath.h"
  18 #define PI 3.141592653
  19
  20 inline bool ifSofty(int Id,int nparams, double *params){
  21   if (nparams<5 && Id==22) return true;   // to remove photons only
  22   for (int i=nparams-1; i>3; i--)
  23     if (Id==params[i]) return true;       // to remove all what is in params from nparams down to 4
  24   return false;
  25 }
  26
  27 // very similar to  MC_FillUserHistogram from Generate.cxx
  28 inline void fillUserHisto(char *name,double val, double weight=1.0,
  29                           double min=Setup::bin_min[0][0],
  30                           double max=Setup::bin_max[0][0]){
  31
  32     TH1D *h=(TH1D*)(Setup::user_histograms->FindObject(name));
  33     if(!h){
  34       h=new TH1D(name,name,Setup::nbins[0][0],min,max);
  35       if(!h) return;
  36       Setup::user_histograms->Add(h);
  37       //      printf("user histogram created %s\n", name);
  38 }
  39     h->Fill(val,weight);
  40
  41 }
  42 double angle(double X, double Y){
  43
  44
  45     double  an=0.0;
  46     double R=sqrt(X*X+Y*Y);
  47     //    if(R<pow(10,-20)) printf(" angle this time X %f\n", 10000*X);
  48     if(R<pow(10,-20)) return an;
  49
  50     if(TMath::Abs(X)/R<0.8)
  51     {
  52         an=acos(X/R);
  53         if(Y<0 && an>0) an=-an;
  54         if(Y>0 && an<0) an=-an;
  55     }
  56     else
  57     {
  58       an=asin(Y/R);
  59       if(X<0 && an>=0.)  an=PI-an;
  60       else if(X<0.)      an=-PI-an;
  61
  62     }
  63   return an;
  64 }
  65
  66 int ZmumuAnalysis(HEPParticle *mother,HEPParticleList *stableDaughters, int nparams, double *params)
  67 {
  68   // THIS IS EXAMPLE of userTreeAnalysis. It acts on list of particle X (under our MC-test) final
  69   // daughters. Of course in real life many options may need to be introduced by the user.
  70     // PARAMETERS:
  71     // params[0] threshold on Energy (or p_T) of particle expressed as a fraction of mother's
  72     //           Energy (or p_T) in mother's frame. If not specified - default is 0.05
  73     //           Note that setting it to zero is possible.
  74     // params[1] maximum number of left soft-suspected particles.
  75     //           0: (default) all listed particles are removed, even if hard
  76     //           1: 2:  one two etc  removable particles are kept (at most)
  77     //          -1: this option is off, all hard particles stay.
  78     // params[2] control tag on discriminating particles,
  79     //           0: (default) Energy in decaying particle rest frame
  80     //           1: Energy in lab frame
  81     //           2: p_T with respect to beam direction in lab frame.
  82     // params[3] type of action to be applied on removed daughters
  83     //           0: (default) nothing is done, they are lost
  84     //           1:  algorithm as in studies on PHOTOS is used, see papers P. Golonka and Z. Was.
  85     // params[4] from this paramter on PDG id-s of particles to be removed can be added.
  86     //           Default is that only photons are removed.
  87     // To get to origin of our particle (mother) one need to go after UserEventAnalysis,
  88     // instructive example will be given later.
  89     assert(mother!=0);
  90     assert(stableDaughters!=0);
  91
  92
  93     double threshold=0.05, leftmax=0.0, selector=0.0, actLost=0.0;
  94     if (nparams>0 && params==0)  return -1;
  95     if (nparams>0) threshold=params[0];
  96     if (nparams>1) leftmax  =params[1];
  97     if (nparams>2) selector =params[2];
  98     if (nparams>3) actLost  =params[3];
  99
 100
 101     HEPParticleList *lostDaughters=new HEPParticleList;
 102     double pt_limit=threshold*mother->GetM();
 103
 104     HEPParticleListIterator daughters (*stableDaughters);
 105     int nphot=0;
 106     double ephot=pow(10,22);
 107     bool redo=false;
 108     for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
 109       if(redo) redo=false;
 110         MC4Vector d4(part->GetE(),part->GetPx(),part->GetPy(),part->GetPz(),part->GetM());
 111         // boost to mother's frame:
 112         d4.Boost(mother->GetPx(),mother->GetPy(),mother->GetPz(),mother->GetE(),mother->GetM());
 113
 114
 115
 116         double p_t=d4.GetX0();  // default
 117         switch ((int) selector)
 118         {
 119         case 1: p_t=part->GetE(); break;
 120         case 2: p_t=d4.Xt(); break;
 121         }
 122
 123
 124         if(ifSofty(part->GetPDGId(),nparams,params)) nphot++;
 125         if ( ifSofty(part->GetPDGId(),nparams,params) && p_t < pt_limit) {
 126           //      printf("Odrzucamy! %i\n",count++);
 127           nphot=0;
 128             lostDaughters->push_back(part);
 129             stableDaughters->remove(part);
 130             part=daughters.first();
 131             redo=true;
 132             continue;
 133         }
 134         if( ifSofty(part->GetPDGId(),nparams,params) && ephot>p_t) ephot=p_t;
 135     }
 136     while(nphot>(int) leftmax)
 137     {
 138       double ephot1=pow(10,22);
 139       redo=false;
 140       for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
 141         if(redo) redo=false;
 142         MC4Vector d4(part->GetE(),part->GetPx(),part->GetPy(),part->GetPz(),part->GetM());
 143         // boost to mother's frame:
 144         d4.Boost(mother->GetPx(),mother->GetPy(),mother->GetPz(),mother->GetE(),mother->GetM());
 145
 146
 147         double p_t=d4.GetX0();  // default
 148         switch ((int) selector)
 149         {
 150         case 1: p_t=part->GetE(); break;
 151         case 2: p_t=d4.Xt(); break;
 152         }
 153
 154
 155         if ( ifSofty(part->GetPDGId(),nparams,params) && p_t == ephot) {
 156
 157           nphot--;
 158             lostDaughters->push_back(part);
 159             stableDaughters->remove(part);
 160             ephot1=pow(10,22);
 161              part=daughters.first();
 162              redo=true;
 163              continue;
 164         }
 165         if( ifSofty(part->GetPDGId(),nparams,params) && ephot1>p_t) ephot1=p_t;
 166       }
 167       ephot=ephot1;
 168
 169     }
 170
 171     // ##############################################################
 172     // SPECIAL CASE - ordered or symmetrical photons
 173     //                This code changes order of photons
 174     // ##############################################################
 175     HEPParticleListIterator symmetry (*stableDaughters);
 176     HEPParticle *phot1 = NULL;
 177     HEPParticle *phot2 = NULL;
 178     for (HEPParticle *part=symmetry.first(); part!=0; part=symmetry.next() )
 179     {
 180       if(part->GetPDGId()==22)
 181       {
 182         if(!phot1) phot1=part;
 183         else
 184         {
 185           phot2=part;
 186           break;
 187         }
 188       }
 189     }
 190     if(phot1 && phot2)
 191     {
 192       if(phot1->GetE()<phot2->GetE()) // for ordered photons
 193       //if( rand()*1.0/RAND_MAX > 0.5) // for photon symmetrization
 194       {
 195         double buf_x = phot1->GetPx();
 196         double buf_y = phot1->GetPy();
 197         double buf_z = phot1->GetPz();
 198         double buf_e = phot1->GetE();
 199
 200         phot1->SetPx(phot2->GetPx());
 201         phot1->SetPy(phot2->GetPy());
 202         phot1->SetPz(phot2->GetPz());
 203         phot1->SetE (phot2->GetE() );
 204
 205         phot2->SetPx(buf_x);
 206         phot2->SetPy(buf_y);
 207         phot2->SetPz(buf_z);
 208         phot2->SetE (buf_e);
 209       }
 210     }
 211
 212     // ##############################################################
 213     // here functionality of removig some daughters is finished
 214     // now we reconsider what to do with them
 215     // ##############################################################
 216     // delete lostDaughters;
 217     // return 1;
 218
 219
 220     // ##############################################################
 221     // Now: What to do with lost daughters?
 222     // ##############################################################
 223     //    lostDaughters->ls();
 224     int version=(int) actLost;
 225
 226     switch(version)
 227       {
 228       case 1:  // add lost to charged
 229         {
 230           HEPParticleListIterator lost (*lostDaughters);
 231           for (HEPParticle *lostpart=lost.first(); lostpart!=0; lostpart=lost.next() ) {
 232             HEPParticle* Best=0;
 233             double searchvirt=pow(10.0,22);
 234             MC4Vector VV;
 235             for( HEPParticle *part=daughters.first(); part!=0; part=daughters.next() ){
 236               if(part->GetCharge()==0) continue;
 237               VV=lostpart->GetP4()+part->GetP4();
 238               VV.AdjustM();
 239               if (VV.GetM()<searchvirt) {searchvirt=VV.GetM(); Best=part;}
 240             }
 241             if(Best) {
 242               Best->SetPx(Best->GetPx()+lostpart->GetPx());
 243               Best->SetPy(Best->GetPy()+lostpart->GetPy());
 244               Best->SetPz(Best->GetPz()+lostpart->GetPz());
 245               Best->SetE (Best->GetE ()+lostpart->GetE ());
 246             }
 247           }
 248         break;
 249         }
 250       default: break; // do nothing
 251       }
 252
 253     delete lostDaughters;
 254
 255
 256     bool activateUserHist=true;
 257     if(!activateUserHist) return 1;
 258
 259     // segmet of user defined histograms
 260
 261     double X=mother->GetPx(), Y=mother->GetPy(), Z=mother->GetPz();
 262     // double E=mother->GetE(),  MM=mother->GetM();
 263
 264
 265     double pt=sqrt(X*X+Y*Y);
 266     double eta=log((sqrt(pt*pt+Z*Z)+TMath::Abs(Z))/pt);
 267     if(Z<0 && eta>0) eta=-eta;
 268     if(Z>0 && eta<0) eta=-eta;
 269     double phi=angle(X,Y);
 270     char hist1[] = "mother-PT";
 271     char hist2[] = "mother-eta";
 272     char hist3[] = "mother-phi";
 273     fillUserHisto(hist1,pt,1.0,0.0,100.0);
 274     fillUserHisto(hist2,eta,1.0,-8.0,8.0);
 275     fillUserHisto(hist3,phi,1.0,-PI,PI);
 276     return 1;
 277 }
 278