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