Port of changes from v4-07-Release and additional rule conformance
[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     fNTracksPassed(0),
92     fNTracksFailed(0),
93     fRemoveKinks(kFALSE),
94     fMostProbable(0)
95 {
96   // Default constructor
97   fNTracksPassed = fNTracksFailed = 0;
98   fCharge = 0;  // takes both charges 0
99   fPt[0]=0.0;              fPt[1] = 100.0;//100
100   fRapidity[0]=-2;       fRapidity[1]=2;//-2 2
101   fPidProbElectron[0]=-1;fPidProbElectron[1]=2;
102   fPidProbPion[0]=-1;    fPidProbPion[1]=2;
103   fPidProbKaon[0]=-1;fPidProbKaon[1]=2;
104   fPidProbProton[0]=-1;fPidProbProton[1]=2;
105   fPidProbMuon[0]=-1;fPidProbMuon[1]=2;
106   fLabel=false;
107   fStatus=0;
108   fminTPCclsF=0;
109   fminITScls=0;
110 }
111 //------------------------------
112 AliFemtoESDTrackCut::~AliFemtoESDTrackCut(){
113   /* noop */
114 }
115 //------------------------------
116 bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
117 {
118   // test the particle and return 
119   // true if it meets all the criteria
120   // false if it doesn't meet at least one of the criteria
121   float tMost[5];
122   
123   //cout<<"AliFemtoESD  cut"<<endl;
124   //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
125   if (fStatus!=0)
126     {
127       //cout<<" status "<<track->Label()<<" "<<track->Flags()<<" "<<track->TPCnclsF()<<" "<<track->ITSncls()<<endl;
128       if ((track->Flags()&fStatus)!=fStatus)
129         {
130           //      cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
131           return false;
132         }
133         
134     }
135   if (fminTPCclsF>track->TPCnclsF())
136     {
137       //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
138       return false;
139     }
140   if (fminITScls>track->ITSncls())
141     {
142       //cout<<" No go because ITS Number of Cls"<<fminITScls<< " "<<track->ITSncls()<<endl;
143       return false;
144     }
145         
146   if (fLabel)
147     {
148       //cout<<"labels"<<endl;
149       if(track->Label()<0)
150         {
151           fNTracksFailed++;
152           //   cout<<"No Go Through the cut"<<endl;
153           //  cout<<fLabel<<" Label="<<track->Label()<<endl;
154           return false;
155         }    
156     }
157   if (fCharge!=0)
158     {              
159       //cout<<"AliFemtoESD  cut ch "<<endl;
160       //cout<<fCharge<<" Charge="<<track->Charge()<<endl;
161       if (track->Charge()!= fCharge)    
162         {
163           fNTracksFailed++;
164           //  cout<<"No Go Through the cut"<<endl;
165           // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
166           return false;
167         }
168     }
169   float tEnergy = ::sqrt(track->P().mag2()+fMass*fMass);
170   float tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
171   float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
172   if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
173     {
174       fNTracksFailed++;
175       //cout<<"No Go Through the cut"<<endl;   
176       //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
177       return false;
178     }
179   if ((tPt<fPt[0])||(tPt>fPt[1]))
180     {
181       fNTracksFailed++;
182       //cout<<"No Go Through the cut"<<endl;
183       //cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
184       return false;
185     }
186 //   cout << "Track has pids: " 
187 //        << track->PidProbElectron() << " " 
188 //        << track->PidProbMuon() << " " 
189 //        << track->PidProbPion() << " " 
190 //        << track->PidProbKaon() << " " 
191 //        << track->PidProbProton() << " " 
192 //        << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
193
194     
195   if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
196     {
197       fNTracksFailed++;
198       //cout<<"No Go Through the cut"<<endl;
199       //cout<<fPidProbElectron[0]<<" < e ="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
200       return false;
201     }
202   if ((track->PidProbPion()<fPidProbPion[0])||(track->PidProbPion()>fPidProbPion[1]))
203     {
204       fNTracksFailed++;
205       //cout<<"No Go Through the cut"<<endl;
206       //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
207       return false;
208     }
209   if ((track->PidProbKaon()<fPidProbKaon[0])||(track->PidProbKaon()>fPidProbKaon[1]))
210     {
211       fNTracksFailed++;
212       //cout<<"No Go Through the cut"<<endl;
213       //cout<<fPidProbKaon[0]<<" < k ="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
214       return false;
215     }
216   if ((track->PidProbProton()<fPidProbProton[0])||(track->PidProbProton()>fPidProbProton[1]))
217     {
218       fNTracksFailed++;
219       //cout<<"No Go Through the cut"<<endl;
220       //cout<<fPidProbProton[0]<<" < p  ="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
221       return false;
222     }
223   if ((track->PidProbMuon()<fPidProbMuon[0])||(track->PidProbMuon()>fPidProbMuon[1]))
224     {
225       fNTracksFailed++;
226       //cout<<"No Go Through the cut"<<endl;
227       //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
228       return false;
229     }
230   if (fRemoveKinks) {
231     if ((track->KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2)))
232       return false;
233   }
234
235   if (fMostProbable) {
236     tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().mag());
237     tMost[1] = 0.0;
238     tMost[2] = track->PidProbPion()*PidFractionPion(track->P().mag());
239     tMost[3] = track->PidProbKaon()*PidFractionKaon(track->P().mag());
240     tMost[4] = track->PidProbProton()*PidFractionProton(track->P().mag());
241     int imost=0;
242     float ipidmax = 0.0;
243     for (int ip=0; ip<5; ip++)
244       if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; };
245     if (imost != fMostProbable) return false;
246   }
247   
248   // cout<<"Go Through the cut"<<endl;
249   // cout<<fLabel<<" Label="<<track->Label()<<endl;
250   // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
251   // cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
252   //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
253   //cout<<fPidProbElectron[0]<<" <  e="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
254   //cout<<fPidProbPion[0]<<" <  pi="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
255   //cout<<fPidProbKaon[0]<<" <  k="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
256   //cout<<fPidProbProton[0]<<" <  p="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
257   //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
258   fNTracksPassed++ ;
259   return true;
260     
261     
262 }
263 //------------------------------
264 AliFemtoString AliFemtoESDTrackCut::Report()
265 {
266   // Prepare report from the execution
267   string tStemp;
268   char tCtemp[100];
269   sprintf(tCtemp,"Particle mass:\t%E\n",this->Mass());
270   tStemp=tCtemp;
271   sprintf(tCtemp,"Particle charge:\t%d\n",fCharge);
272   tStemp+=tCtemp;
273   sprintf(tCtemp,"Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
274   tStemp+=tCtemp;
275   sprintf(tCtemp,"Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
276   tStemp+=tCtemp;
277   sprintf(tCtemp,"Number of tracks which passed:\t%ld  Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
278   tStemp += tCtemp;
279   AliFemtoString returnThis = tStemp;
280   return returnThis;
281 }
282 TList *AliFemtoESDTrackCut::ListSettings()
283 {
284   // return a list of settings in a writable form
285   TList *tListSetttings = new TList();
286   char buf[200];
287   snprintf(buf, 200, "AliFemtoESDTrackCut.mass=%lf", this->Mass());
288   tListSetttings->AddLast(new TObjString(buf));
289
290   snprintf(buf, 200, "AliFemtoESDTrackCut.charge=%i", fCharge);
291   tListSetttings->AddLast(new TObjString(buf));
292   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.minimum=%lf", fPidProbPion[0]);
293   tListSetttings->AddLast(new TObjString(buf));
294   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.maximum=%lf", fPidProbPion[1]);
295   tListSetttings->AddLast(new TObjString(buf));
296   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.minimum=%lf", fPidProbKaon[0]);
297   tListSetttings->AddLast(new TObjString(buf));
298   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.maximum=%lf", fPidProbKaon[1]);
299   tListSetttings->AddLast(new TObjString(buf));
300   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.minimum=%lf", fPidProbProton[0]);
301   tListSetttings->AddLast(new TObjString(buf));
302   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.maximum=%lf", fPidProbProton[1]);
303   tListSetttings->AddLast(new TObjString(buf));
304   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.minimum=%lf", fPidProbElectron[0]);
305   tListSetttings->AddLast(new TObjString(buf));
306   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.maximum=%lf", fPidProbElectron[1]);
307   tListSetttings->AddLast(new TObjString(buf));
308   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.minimum=%lf", fPidProbMuon[0]);
309   tListSetttings->AddLast(new TObjString(buf));
310   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.maximum=%lf", fPidProbMuon[1]);
311   tListSetttings->AddLast(new TObjString(buf));
312   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumtpcclusters=%i", fminTPCclsF);
313   tListSetttings->AddLast(new TObjString(buf));
314   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumitsclusters=%i", fminTPCclsF);
315   tListSetttings->AddLast(new TObjString(buf));
316   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.minimum=%lf", fPt[0]);
317   tListSetttings->AddLast(new TObjString(buf));
318   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.maximum=%lf", fPt[1]);
319   tListSetttings->AddLast(new TObjString(buf));
320   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.minimum=%lf", fRapidity[0]);
321   tListSetttings->AddLast(new TObjString(buf));
322   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.maximum=%lf", fRapidity[1]);
323   tListSetttings->AddLast(new TObjString(buf));
324   snprintf(buf, 200, "AliFemtoESDTrackCut.removekinks=%i", fRemoveKinks);
325   tListSetttings->AddLast(new TObjString(buf));
326   if (fMostProbable) {
327     if (fMostProbable == 2)
328       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Pion");
329     if (fMostProbable == 3)
330       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Kaon");
331     if (fMostProbable == 4)
332       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Proton");
333     tListSetttings->AddLast(new TObjString(buf));
334   }
335   return tListSetttings;
336 }
337 void AliFemtoESDTrackCut::SetRemoveKinks(const bool& flag)
338 {
339   fRemoveKinks = flag;
340 }
341                             
342                             // electron
343 // 0.13 - 1.8
344 // 0       7.594129e-02    8.256141e-03
345 // 1       -5.535827e-01   8.170825e-02
346 // 2       1.728591e+00    3.104210e-01
347 // 3       -2.827893e+00   5.827802e-01
348 // 4       2.503553e+00    5.736207e-01
349 // 5       -1.125965e+00   2.821170e-01
350 // 6       2.009036e-01    5.438876e-02
351 float AliFemtoESDTrackCut::PidFractionElectron(float mom) const
352 {
353   // Provide a parameterized fraction of electrons dependent on momentum
354   if (mom<0.13) return 0.0;
355   if (mom>1.8) return 0.0;
356   return (7.594129e-02 
357           -5.535827e-01*mom        
358           +1.728591e+00*mom*mom    
359           -2.827893e+00*mom*mom*mom 
360           +2.503553e+00*mom*mom*mom*mom    
361           -1.125965e+00*mom*mom*mom*mom*mom      
362           +2.009036e-01*mom*mom*mom*mom*mom*mom);   
363 }
364
365 // pion
366 // 0.13 - 2.0
367 // 0       1.063457e+00    8.872043e-03
368 // 1       -4.222208e-01   2.534402e-02
369 // 2       1.042004e-01    1.503945e-02
370 float AliFemtoESDTrackCut::PidFractionPion(float mom) const
371 {
372   // Provide a parameterized fraction of pions dependent on momentum
373   if (mom<0.13) return 0.0;
374   if (mom>2.0) return 0.0;
375   return ( 1.063457e+00
376            -4.222208e-01*mom
377            +1.042004e-01*mom*mom);
378 }
379
380 // kaon
381 // 0.18 - 2.0
382 // 0       -7.289406e-02   1.686074e-03
383 // 1       4.415666e-01    1.143939e-02
384 // 2       -2.996790e-01   1.840964e-02
385 // 3       6.704652e-02    7.783990e-03
386 float AliFemtoESDTrackCut::PidFractionKaon(float mom) const
387 {
388   // Provide a parameterized fraction of kaons dependent on momentum
389   if (mom<0.18) return 0.0;
390   if (mom>2.0) return 0.0;
391   return (-7.289406e-02
392           +4.415666e-01*mom        
393           -2.996790e-01*mom*mom    
394           +6.704652e-02*mom*mom*mom);
395 }
396
397 // proton
398 // 0.26 - 2.0
399 // 0       -3.730200e-02   2.347311e-03
400 // 1       1.163684e-01    1.319316e-02
401 // 2       8.354116e-02    1.997948e-02
402 // 3       -4.608098e-02   8.336400e-03
403 float AliFemtoESDTrackCut::PidFractionProton(float mom) const
404 {
405   // Provide a parameterized fraction of protons dependent on momentum
406   if (mom<0.26) return  0.0;
407   if (mom>2.0) return 0.0;
408   return (-3.730200e-02  
409           +1.163684e-01*mom           
410           +8.354116e-02*mom*mom       
411           -4.608098e-02*mom*mom*mom);  
412 }