1 //#include "UserTreeAnalysis.H" // remove if copied to user working directory
6 #include "HEPParticle.H"
11 #define PI 3.141592653
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
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]){
25 TH1D *h=(TH1D*)(Setup::user_histograms->FindObject(name));
27 h=new TH1D(name,name,Setup::nbins[0][0],min,max);
29 Setup::user_histograms->Add(h);
30 // printf("user histogram created %s\n", name);
35 double angle(double X, double Y){
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;
43 if(TMath::Abs(X)/R<0.8)
46 if(Y<0 && an>0) an=-an;
47 if(Y>0 && an<0) an=-an;
52 if(X<0 && an>=0.) an=PI-an;
53 else if(X<0.) an=-PI-an;
61 int ZtautauAnalysis(HEPParticle *mother,HEPParticleList *stableDaughters, int nparams, double *params)
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.
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.
85 assert(stableDaughters!=0);
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];
96 HEPParticleList *lostDaughters=new HEPParticleList;
97 double pt_limit=threshold*mother->GetM();
99 HEPParticleListIterator daughters (*stableDaughters);
101 double ephot=pow(10,22);
103 for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
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());
111 double p_t=d4.GetX0(); // default
112 switch ((int) selector)
114 case 1: p_t=part->GetE(); break;
115 case 2: p_t=d4.Xt(); break;
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++);
123 lostDaughters->push_back(part);
124 stableDaughters->remove(part);
125 part=daughters.first();
129 if( ifSofty(part->GetPDGId(),nparams,params) && ephot>p_t) ephot=p_t;
131 while(nphot>(int) leftmax)
133 double ephot1=pow(10,22);
135 for (HEPParticle *part=daughters.first(); part!=0; (redo)? part=part:part=daughters.next() ) {
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());
142 double p_t=d4.GetX0(); // default
143 switch ((int) selector)
145 case 1: p_t=part->GetE(); break;
146 case 2: p_t=d4.Xt(); break;
150 if ( ifSofty(part->GetPDGId(),nparams,params) && p_t == ephot) {
153 lostDaughters->push_back(part);
154 stableDaughters->remove(part);
156 part=daughters.first();
160 if( ifSofty(part->GetPDGId(),nparams,params) && ephot1>p_t) ephot1=p_t;
166 // ##############################################################
167 // here functionality of removig some daughters is finished
168 // now we reconsider what to do with them
169 // ##############################################################
170 // delete lostDaughters;
174 // ##############################################################
175 // Now: What to do with lost daughters?
176 // ##############################################################
177 // lostDaughters->ls();
178 int version=(int) actLost;
182 case 1: // add lost to charged
184 HEPParticleListIterator lost (*lostDaughters);
185 for (HEPParticle *lostpart=lost.first(); lostpart!=0; lostpart=lost.next() ) {
187 double searchvirt=pow(10.0,22);
189 for( HEPParticle *part=daughters.first(); part!=0; part=daughters.next() ){
190 if(part->GetCharge()==0) continue;
191 VV=lostpart->GetP4()+part->GetP4();
193 if (VV.GetM()<searchvirt) {searchvirt=VV.GetM(); Best=part;}
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 ());
204 default: break; // do nothing
207 delete lostDaughters;
213 double px=0,py=0,pz=0,e=0;
214 double px1=0,py1=0,pz1=0,e1=0;
216 //loop over all daughters
217 for(HEPParticle *part=daughters.first();part!=0;part=daughters.next())
219 if(abs(part->GetPDGId())==211)
226 if(abs(part->GetPDGId())==211) continue;
227 //Sum e+ e- 4-vectors
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);
236 char mm[] = "pi-energy";
237 fillUserHisto(mm,mass,1.0,0.0,1.1);
239 //For ZtautauAnalysis - other histograms are not needed
240 bool activateUserHist=false;
241 if(!activateUserHist) return 1;
243 // segmet of user defined histograms
245 double X=mother->GetPx(), Y=mother->GetPy(), Z=mother->GetPz();
246 // double E=mother->GetE(), MM=mother->GetM();
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);