Bring AliFemto up to date with latest code developements
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoESDTrackCut.cxx
1 /***************************************************************************
2  *
3  * $Id$ 
4  *
5  * 
6  ***************************************************************************
7  *
8  * 
9  *              
10  *
11  ***************************************************************************
12  *
13  * $Log$
14  * Revision 1.3  2007/05/22 09:01:42  akisiel
15  * Add the possibiloity to save cut settings in the ROOT file
16  *
17  * Revision 1.2  2007/05/21 10:38:25  akisiel
18  * More coding rule conformance
19  *
20  * Revision 1.1  2007/05/16 10:25:06  akisiel
21  * Making the directory structure of AliFemtoUser flat. All files go into one common directory
22  *
23  * Revision 1.4  2007/05/03 09:46:10  akisiel
24  * Fixing Effective C++ warnings
25  *
26  * Revision 1.3  2007/04/27 07:25:59  akisiel
27  * Make revisions needed for compilation from the main AliRoot tree
28  *
29  * Revision 1.1.1.1  2007/04/25 15:38:41  panos
30  * Importing the HBT code dir
31  *
32  * Revision 1.4  2007-04-03 16:00:08  mchojnacki
33  * Changes to iprove memory managing
34  *
35  * Revision 1.3  2007/03/13 15:30:03  mchojnacki
36  * adding reader for simulated data
37  *
38  * Revision 1.2  2007/03/08 14:58:03  mchojnacki
39  * adding some alice stuff
40  *
41  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
42  * First version on CVS
43  *
44  **************************************************************************/
45
46 #include "AliFemtoESDTrackCut.h"
47 #include <cstdio>
48
49 #ifdef __ROOT__ 
50 ClassImp(AliFemtoESDTrackCut)
51 #endif
52
53
54 // electron
55 // 0.13 - 1.8
56 // 0       7.594129e-02    8.256141e-03
57 // 1       -5.535827e-01   8.170825e-02
58 // 2       1.728591e+00    3.104210e-01
59 // 3       -2.827893e+00   5.827802e-01
60 // 4       2.503553e+00    5.736207e-01
61 // 5       -1.125965e+00   2.821170e-01
62 // 6       2.009036e-01    5.438876e-02
63
64 // pion
65 // 0.13 - 2.0
66 // 0       1.063457e+00    8.872043e-03
67 // 1       -4.222208e-01   2.534402e-02
68 // 2       1.042004e-01    1.503945e-02
69
70 // kaon
71 // 0.18 - 2.0
72 // 0       -7.289406e-02   1.686074e-03
73 // 1       4.415666e-01    1.143939e-02
74 // 2       -2.996790e-01   1.840964e-02
75 // 3       6.704652e-02    7.783990e-03
76
77 // proton
78 // 0.26 - 2.0
79 // 0       -3.730200e-02   2.347311e-03
80 // 1       1.163684e-01    1.319316e-02
81 // 2       8.354116e-02    1.997948e-02
82 // 3       -4.608098e-02   8.336400e-03
83
84
85 AliFemtoESDTrackCut::AliFemtoESDTrackCut() :
86     fCharge(0),
87     fLabel(0),
88     fStatus(0),
89     fminTPCclsF(0),
90     fminITScls(0),
91     fMaxITSchiNdof(1000.0),
92     fMaxTPCchiNdof(1000.0),
93     fMaxSigmaToVertex(1000.0),
94     fNTracksPassed(0),
95     fNTracksFailed(0),
96     fRemoveKinks(kFALSE),
97     fMostProbable(0)
98 {
99   // Default constructor
100   fNTracksPassed = fNTracksFailed = 0;
101   fCharge = 0;  // takes both charges 0
102   fPt[0]=0.0;              fPt[1] = 100.0;//100
103   fRapidity[0]=-2;       fRapidity[1]=2;//-2 2
104   fPidProbElectron[0]=-1;fPidProbElectron[1]=2;
105   fPidProbPion[0]=-1;    fPidProbPion[1]=2;
106   fPidProbKaon[0]=-1;fPidProbKaon[1]=2;
107   fPidProbProton[0]=-1;fPidProbProton[1]=2;
108   fPidProbMuon[0]=-1;fPidProbMuon[1]=2;
109   fLabel=false;
110   fStatus=0;
111   fminTPCclsF=0;
112   fminITScls=0;
113 }
114 //------------------------------
115 AliFemtoESDTrackCut::~AliFemtoESDTrackCut(){
116   /* noop */
117 }
118 //------------------------------
119 bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
120 {
121   // test the particle and return 
122   // true if it meets all the criteria
123   // false if it doesn't meet at least one of the criteria
124   float tMost[5];
125   
126   //cout<<"AliFemtoESD  cut"<<endl;
127   //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
128   if (fStatus!=0)
129     {
130       //cout<<" status "<<track->Label()<<" "<<track->Flags()<<" "<<track->TPCnclsF()<<" "<<track->ITSncls()<<endl;
131       if ((track->Flags()&fStatus)!=fStatus)
132         {
133           //      cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
134           return false;
135         }
136         
137     }
138   if (fRemoveKinks) {
139     if ((track->KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2)))
140       return false;
141   }
142   if (fminTPCclsF>track->TPCnclsF())
143     {
144       //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
145       return false;
146     }
147   if (fminTPCncls>track->TPCncls())
148     {
149       //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
150       return false;
151     }
152   if (fminITScls>track->ITSncls())
153     {
154       //cout<<" No go because ITS Number of Cls"<<fminITScls<< " "<<track->ITSncls()<<endl;
155       return false;
156     }
157         
158   if (fMaxSigmaToVertex < track->SigmaToVertex()) {
159     return false;
160   }
161   
162   if (track->ITSncls() > 0) 
163     if ((track->ITSchi2()/track->ITSncls()) > fMaxITSchiNdof) {
164       return false;
165     }
166
167   if (track->TPCncls() > 0)
168     if ((track->TPCchi2()/track->TPCncls()) > fMaxTPCchiNdof) {
169       return false;
170     }
171
172   if (fLabel)
173     {
174       //cout<<"labels"<<endl;
175       if(track->Label()<0)
176         {
177           fNTracksFailed++;
178           //   cout<<"No Go Through the cut"<<endl;
179           //  cout<<fLabel<<" Label="<<track->Label()<<endl;
180           return false;
181         }    
182     }
183   if (fCharge!=0)
184     {              
185       //cout<<"AliFemtoESD  cut ch "<<endl;
186       //cout<<fCharge<<" Charge="<<track->Charge()<<endl;
187       if (track->Charge()!= fCharge)    
188         {
189           fNTracksFailed++;
190           //  cout<<"No Go Through the cut"<<endl;
191           // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
192           return false;
193         }
194     }
195   float tEnergy = ::sqrt(track->P().mag2()+fMass*fMass);
196   float tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
197   float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
198   if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
199     {
200       fNTracksFailed++;
201       //cout<<"No Go Through the cut"<<endl;   
202       //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
203       return false;
204     }
205   if ((tPt<fPt[0])||(tPt>fPt[1]))
206     {
207       fNTracksFailed++;
208       //cout<<"No Go Through the cut"<<endl;
209       //cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
210       return false;
211     }
212 //   cout << "Track has pids: " 
213 //        << track->PidProbElectron() << " " 
214 //        << track->PidProbMuon() << " " 
215 //        << track->PidProbPion() << " " 
216 //        << track->PidProbKaon() << " " 
217 //        << track->PidProbProton() << " " 
218 //        << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
219
220     
221   if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
222     {
223       fNTracksFailed++;
224       //cout<<"No Go Through the cut"<<endl;
225       //cout<<fPidProbElectron[0]<<" < e ="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
226       return false;
227     }
228   if ((track->PidProbPion()<fPidProbPion[0])||(track->PidProbPion()>fPidProbPion[1]))
229     {
230       fNTracksFailed++;
231       //cout<<"No Go Through the cut"<<endl;
232       //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
233       return false;
234     }
235   if ((track->PidProbKaon()<fPidProbKaon[0])||(track->PidProbKaon()>fPidProbKaon[1]))
236     {
237       fNTracksFailed++;
238       //cout<<"No Go Through the cut"<<endl;
239       //cout<<fPidProbKaon[0]<<" < k ="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
240       return false;
241     }
242   if ((track->PidProbProton()<fPidProbProton[0])||(track->PidProbProton()>fPidProbProton[1]))
243     {
244       fNTracksFailed++;
245       //cout<<"No Go Through the cut"<<endl;
246       //cout<<fPidProbProton[0]<<" < p  ="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
247       return false;
248     }
249   if ((track->PidProbMuon()<fPidProbMuon[0])||(track->PidProbMuon()>fPidProbMuon[1]))
250     {
251       fNTracksFailed++;
252       //cout<<"No Go Through the cut"<<endl;
253       //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
254       return false;
255     }
256
257   if (fMostProbable) {
258     tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().mag());
259     tMost[1] = 0.0;
260     tMost[2] = track->PidProbPion()*PidFractionPion(track->P().mag());
261     tMost[3] = track->PidProbKaon()*PidFractionKaon(track->P().mag());
262     tMost[4] = track->PidProbProton()*PidFractionProton(track->P().mag());
263     int imost=0;
264     float ipidmax = 0.0;
265     for (int ip=0; ip<5; ip++)
266       if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; };
267     if (imost != fMostProbable) return false;
268   }
269   
270   // cout<<"Go Through the cut"<<endl;
271   // cout<<fLabel<<" Label="<<track->Label()<<endl;
272   // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
273   // cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
274   //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
275   //cout<<fPidProbElectron[0]<<" <  e="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
276   //cout<<fPidProbPion[0]<<" <  pi="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
277   //cout<<fPidProbKaon[0]<<" <  k="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
278   //cout<<fPidProbProton[0]<<" <  p="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
279   //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
280   fNTracksPassed++ ;
281   return true;
282     
283     
284 }
285 //------------------------------
286 AliFemtoString AliFemtoESDTrackCut::Report()
287 {
288   // Prepare report from the execution
289   string tStemp;
290   char tCtemp[100];
291   sprintf(tCtemp,"Particle mass:\t%E\n",this->Mass());
292   tStemp=tCtemp;
293   sprintf(tCtemp,"Particle charge:\t%d\n",fCharge);
294   tStemp+=tCtemp;
295   sprintf(tCtemp,"Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
296   tStemp+=tCtemp;
297   sprintf(tCtemp,"Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
298   tStemp+=tCtemp;
299   sprintf(tCtemp,"Number of tracks which passed:\t%ld  Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
300   tStemp += tCtemp;
301   AliFemtoString returnThis = tStemp;
302   return returnThis;
303 }
304 TList *AliFemtoESDTrackCut::ListSettings()
305 {
306   // return a list of settings in a writable form
307   TList *tListSetttings = new TList();
308   char buf[200];
309   snprintf(buf, 200, "AliFemtoESDTrackCut.mass=%lf", this->Mass());
310   tListSetttings->AddLast(new TObjString(buf));
311
312   snprintf(buf, 200, "AliFemtoESDTrackCut.charge=%i", fCharge);
313   tListSetttings->AddLast(new TObjString(buf));
314   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.minimum=%lf", fPidProbPion[0]);
315   tListSetttings->AddLast(new TObjString(buf));
316   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.maximum=%lf", fPidProbPion[1]);
317   tListSetttings->AddLast(new TObjString(buf));
318   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.minimum=%lf", fPidProbKaon[0]);
319   tListSetttings->AddLast(new TObjString(buf));
320   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.maximum=%lf", fPidProbKaon[1]);
321   tListSetttings->AddLast(new TObjString(buf));
322   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.minimum=%lf", fPidProbProton[0]);
323   tListSetttings->AddLast(new TObjString(buf));
324   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.maximum=%lf", fPidProbProton[1]);
325   tListSetttings->AddLast(new TObjString(buf));
326   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.minimum=%lf", fPidProbElectron[0]);
327   tListSetttings->AddLast(new TObjString(buf));
328   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.maximum=%lf", fPidProbElectron[1]);
329   tListSetttings->AddLast(new TObjString(buf));
330   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.minimum=%lf", fPidProbMuon[0]);
331   tListSetttings->AddLast(new TObjString(buf));
332   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.maximum=%lf", fPidProbMuon[1]);
333   tListSetttings->AddLast(new TObjString(buf));
334   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumtpcclusters=%i", fminTPCclsF);
335   tListSetttings->AddLast(new TObjString(buf));
336   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumitsclusters=%i", fminTPCclsF);
337   tListSetttings->AddLast(new TObjString(buf));
338   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.minimum=%lf", fPt[0]);
339   tListSetttings->AddLast(new TObjString(buf));
340   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.maximum=%lf", fPt[1]);
341   tListSetttings->AddLast(new TObjString(buf));
342   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.minimum=%lf", fRapidity[0]);
343   tListSetttings->AddLast(new TObjString(buf));
344   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.maximum=%lf", fRapidity[1]);
345   tListSetttings->AddLast(new TObjString(buf));
346   snprintf(buf, 200, "AliFemtoESDTrackCut.removekinks=%i", fRemoveKinks);
347   tListSetttings->AddLast(new TObjString(buf));
348   snprintf(buf, 200, "AliFemtoESDTrackCut.maxitschindof=%lf", fMaxITSchiNdof);
349   tListSetttings->AddLast(new TObjString(buf));
350   snprintf(buf, 200, "AliFemtoESDTrackCut.maxtpcchindof=%lf", fMaxTPCchiNdof);
351   tListSetttings->AddLast(new TObjString(buf));
352   snprintf(buf, 200, "AliFemtoESDTrackCut.maxsigmatovertex=%lf", fMaxSigmaToVertex);
353   tListSetttings->AddLast(new TObjString(buf));
354   if (fMostProbable) {
355     if (fMostProbable == 2)
356       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Pion");
357     if (fMostProbable == 3)
358       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Kaon");
359     if (fMostProbable == 4)
360       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Proton");
361     tListSetttings->AddLast(new TObjString(buf));
362   }
363   return tListSetttings;
364 }
365 void AliFemtoESDTrackCut::SetRemoveKinks(const bool& flag)
366 {
367   fRemoveKinks = flag;
368 }
369                             
370                             // electron
371 // 0.13 - 1.8
372 // 0       7.594129e-02    8.256141e-03
373 // 1       -5.535827e-01   8.170825e-02
374 // 2       1.728591e+00    3.104210e-01
375 // 3       -2.827893e+00   5.827802e-01
376 // 4       2.503553e+00    5.736207e-01
377 // 5       -1.125965e+00   2.821170e-01
378 // 6       2.009036e-01    5.438876e-02
379 float AliFemtoESDTrackCut::PidFractionElectron(float mom) const
380 {
381   // Provide a parameterized fraction of electrons dependent on momentum
382   if (mom<0.13) return 0.0;
383   if (mom>1.8) return 0.0;
384   return (7.594129e-02 
385           -5.535827e-01*mom        
386           +1.728591e+00*mom*mom    
387           -2.827893e+00*mom*mom*mom 
388           +2.503553e+00*mom*mom*mom*mom    
389           -1.125965e+00*mom*mom*mom*mom*mom      
390           +2.009036e-01*mom*mom*mom*mom*mom*mom);   
391 }
392
393 // pion
394 // 0.13 - 2.0
395 // 0       1.063457e+00    8.872043e-03
396 // 1       -4.222208e-01   2.534402e-02
397 // 2       1.042004e-01    1.503945e-02
398 float AliFemtoESDTrackCut::PidFractionPion(float mom) const
399 {
400   // Provide a parameterized fraction of pions dependent on momentum
401   if (mom<0.13) return 0.0;
402   if (mom>2.0) return 0.0;
403   return ( 1.063457e+00
404            -4.222208e-01*mom
405            +1.042004e-01*mom*mom);
406 }
407
408 // kaon
409 // 0.18 - 2.0
410 // 0       -7.289406e-02   1.686074e-03
411 // 1       4.415666e-01    1.143939e-02
412 // 2       -2.996790e-01   1.840964e-02
413 // 3       6.704652e-02    7.783990e-03
414 float AliFemtoESDTrackCut::PidFractionKaon(float mom) const
415 {
416   // Provide a parameterized fraction of kaons dependent on momentum
417   if (mom<0.18) return 0.0;
418   if (mom>2.0) return 0.0;
419   return (-7.289406e-02
420           +4.415666e-01*mom        
421           -2.996790e-01*mom*mom    
422           +6.704652e-02*mom*mom*mom);
423 }
424
425 // proton
426 // 0.26 - 2.0
427 // 0       -3.730200e-02   2.347311e-03
428 // 1       1.163684e-01    1.319316e-02
429 // 2       8.354116e-02    1.997948e-02
430 // 3       -4.608098e-02   8.336400e-03
431 float AliFemtoESDTrackCut::PidFractionProton(float mom) const
432 {
433   // Provide a parameterized fraction of protons dependent on momentum
434   if (mom<0.26) return  0.0;
435   if (mom>2.0) return 0.0;
436   return (-3.730200e-02  
437           +1.163684e-01*mom           
438           +8.354116e-02*mom*mom       
439           -4.608098e-02*mom*mom*mom);  
440 }