In AOD Event Reader: TOF mismatch with new method + macros for proton femtoscopy...
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoESDTrackCut.cxx
1 /*
2 ***************************************************************************
3 *
4 * $Id$
5 *
6 *
7 ***************************************************************************
8 *
9 *
10 *
11 *
12 ***************************************************************************
13 *
14 * $Log$
15 * Revision 1.3  2007/05/22 09:01:42  akisiel
16 * Add the possibiloity to save cut settings in the ROOT file
17 *
18 * Revision 1.2  2007/05/21 10:38:25  akisiel
19 * More coding rule conformance
20 *
21 * Revision 1.1  2007/05/16 10:25:06  akisiel
22 * Making the directory structure of AliFemtoUser flat. All files go into one common directory
23 *
24 * Revision 1.4  2007/05/03 09:46:10  akisiel
25 * Fixing Effective C++ warnings
26 *
27 * Revision 1.3  2007/04/27 07:25:59  akisiel
28 * Make revisions needed for compilation from the main AliRoot tree
29 *
30 * Revision 1.1.1.1  2007/04/25 15:38:41  panos
31 * Importing the HBT code dir
32 *
33 * Revision 1.4  2007-04-03 16:00:08  mchojnacki
34 * Changes to iprove memory managing
35 *
36 * Revision 1.3  2007/03/13 15:30:03  mchojnacki
37 * adding reader for simulated data
38 *
39 * Revision 1.2  2007/03/08 14:58:03  mchojnacki
40 * adding some alice stuff
41 *
42 * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
43 * First version on CVS
44 *
45 **************************************************************************/
46
47 #include "AliFemtoESDTrackCut.h"
48 #include <cstdio>
49
50 #ifdef __ROOT__
51 ClassImp(AliFemtoESDTrackCut)
52 #endif
53
54
55 // electron
56 // 0.13 - 1.8
57 // 0       7.594129e-02    8.256141e-03
58 // 1       -5.535827e-01   8.170825e-02
59 // 2       1.728591e+00    3.104210e-01
60 // 3       -2.827893e+00   5.827802e-01
61 // 4       2.503553e+00    5.736207e-01
62 // 5       -1.125965e+00   2.821170e-01
63 // 6       2.009036e-01    5.438876e-02
64
65 // pion
66 // 0.13 - 2.0
67 // 0       1.063457e+00    8.872043e-03
68 // 1       -4.222208e-01   2.534402e-02
69 // 2       1.042004e-01    1.503945e-02
70
71 // kaon
72 // 0.18 - 2.0
73 // 0       -7.289406e-02   1.686074e-03
74 // 1       4.415666e-01    1.143939e-02
75 // 2       -2.996790e-01   1.840964e-02
76 // 3       6.704652e-02    7.783990e-03
77
78 // proton
79 // 0.26 - 2.0
80 // 0       -3.730200e-02   2.347311e-03
81 // 1       1.163684e-01    1.319316e-02
82 // 2       8.354116e-02    1.997948e-02
83 // 3       -4.608098e-02   8.336400e-03
84
85
86   AliFemtoESDTrackCut::AliFemtoESDTrackCut() :
87     fCharge(0),
88     fLabel(0),
89     fStatus(0),
90     fPIDMethod(knSigma),
91   fNsigmaTPCTOF(kFALSE),
92   fNsigmaTPConly(kFALSE),
93   fNsigma(3.),
94     fminTPCclsF(0),
95     fminTPCncls(0),
96     fminITScls(0),
97     fMaxITSchiNdof(1000.0),
98     fMaxTPCchiNdof(1000.0),
99     fMaxSigmaToVertex(1000.0),
100     fNTracksPassed(0),
101     fNTracksFailed(0),
102     fRemoveKinks(kFALSE),
103     fRemoveITSFake(kFALSE),
104     fMostProbable(0),
105     fMaxImpactXY(1000.0),
106   fMinImpactXY(-1000.0),
107     fMaxImpactZ(1000.0),
108     fMaxImpactXYPtOff(1000.0),
109     fMaxImpactXYPtNrm(1000.0),
110     fMaxImpactXYPtPow(1000.0),
111     fMinPforTOFpid(0.0),
112     fMaxPforTOFpid(10000.0),
113     fMinPforTPCpid(0.0),
114     fMaxPforTPCpid(10000.0),
115     fMinPforITSpid(0.0),
116     fMaxPforITSpid(10000.0)
117 {
118   // Default constructor
119   fNTracksPassed = fNTracksFailed = 0;
120   fCharge = 0;  // takes both charges 0
121   fPt[0]=0.0;              fPt[1] = 100.0;//100
122   fRapidity[0]=-2;       fRapidity[1]=2;//-2 2
123   fEta[0]=-2;       fEta[1]=2;//-2 2
124   fPidProbElectron[0]=-1;fPidProbElectron[1]=2;
125   fPidProbPion[0]=-1;    fPidProbPion[1]=2;
126   fPidProbKaon[0]=-1;fPidProbKaon[1]=2;
127   fPidProbProton[0]=-1;fPidProbProton[1]=2;
128   fPidProbMuon[0]=-1;fPidProbMuon[1]=2;
129   for (Int_t i = 0; i < 3; i++)
130     fCutClusterRequirementITS[i] = AliESDtrackCuts::kOff;
131   fLabel=false;
132   fStatus=0;
133   fminTPCclsF=0;
134   fminITScls=0;
135   fPIDMethod=knSigma;
136   fNsigmaTPCTOF=kFALSE;
137   fNsigmaTPConly=kFALSE;
138   fNsigma=3.;
139 }
140 //------------------------------
141 AliFemtoESDTrackCut::~AliFemtoESDTrackCut(){
142   /* noop */
143 }
144 //------------------------------
145 bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
146 {
147   //cout<<"AliFemtoESDTrackCut::Pass"<<endl;
148
149   // test the particle and return
150   // true if it meets all the criteria
151   // false if it doesn't meet at least one of the criteria
152   float tMost[5];
153
154   //cout<<"AliFemtoESD  cut"<<endl;
155   //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
156   if (fStatus!=0)
157     {
158       //cout<<" status "<<track->Label()<<" "<<track->Flags()<<" "<<track->TPCnclsF()<<" "<<track->ITSncls()<<endl;
159       if ((track->Flags()&fStatus)!=fStatus)
160         {
161           //      cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
162           return false;
163         }
164
165     }
166   if (fRemoveKinks) {
167     if ((track->KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2)))
168       return false;
169   }
170   if (fRemoveITSFake) {
171     if (track->ITSncls() < 0)
172       return false;
173   }
174   if (fminTPCclsF>track->TPCnclsF())
175     {
176       //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
177       return false;
178     }
179   if (fminTPCncls>track->TPCncls())
180     {
181       //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
182       return false;
183     }
184   if (fminITScls>track->ITSncls())
185     {
186       //cout<<" No go because ITS Number of Cls"<<fminITScls<< " "<<track->ITSncls()<<endl;
187       return false;
188     }
189
190   if (fMaxImpactXY < TMath::Abs(track->ImpactD()))
191     return false;
192
193   if (fMinImpactXY > TMath::Abs(track->ImpactD()))
194     return false;
195
196   if (fMaxImpactZ < TMath::Abs(track->ImpactZ()))
197     return false;
198
199   if (fMaxSigmaToVertex < track->SigmaToVertex()) {
200     return false;
201   }
202
203   if (track->ITSncls() > 0)
204     if ((track->ITSchi2()/track->ITSncls()) > fMaxITSchiNdof) {
205       return false;
206     }
207
208   if (track->TPCncls() > 0)
209     if ((track->TPCchi2()/track->TPCncls()) > fMaxTPCchiNdof) {
210       return false;
211     }
212   //ITS cluster requirenments
213   for (Int_t i = 0; i < 3; i++)
214     if(!CheckITSClusterRequirement(fCutClusterRequirementITS[i], track->HasPointOnITSLayer(i*2), track->HasPointOnITSLayer(i*2+1)))
215       return false;
216
217   if (fLabel)
218     {
219       //cout<<"labels"<<endl;
220       if(track->Label()<0)
221         {
222           fNTracksFailed++;
223           //   cout<<"No Go Through the cut"<<endl;
224           //  cout<<fLabel<<" Label="<<track->Label()<<endl;
225           return false;
226         }
227     }
228   if (fCharge!=0)
229     {
230       //cout<<"AliFemtoESD  cut ch "<<endl;
231       //cout<<fCharge<<" Charge="<<track->Charge()<<endl;
232       if (track->Charge()!= fCharge)
233         {
234           fNTracksFailed++;
235           //  cout<<"No Go Through the cut"<<endl;
236           // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
237           return false;
238         }
239     }
240
241
242
243
244   Bool_t tTPCPidIn = (track->Flags()&AliFemtoTrack::kTPCpid)>0;
245   Bool_t tITSPidIn = (track->Flags()&AliFemtoTrack::kITSpid)>0;
246   Bool_t tTOFPidIn = (track->Flags()&AliFemtoTrack::kTOFpid)>0;
247
248   if(fMinPforTOFpid > 0 && track->P().Mag() > fMinPforTOFpid &&
249      track->P().Mag() < fMaxPforTOFpid && !tTOFPidIn)
250     {
251       fNTracksFailed++;
252       return false;
253     }
254
255   if(fMinPforTPCpid > 0 && track->P().Mag() > fMinPforTPCpid &&
256      track->P().Mag() < fMaxPforTPCpid && !tTPCPidIn)
257     {
258       fNTracksFailed++;
259       return false;
260     }
261
262   if(fMinPforITSpid > 0 && track->P().Mag() > fMinPforITSpid &&
263      track->P().Mag() < fMaxPforITSpid && !tITSPidIn)
264     {
265       fNTracksFailed++;
266       return false;
267     }
268
269
270   float tEnergy = ::sqrt(track->P().Mag2()+fMass*fMass);
271   float tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
272   float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
273   float tEta = track->P().PseudoRapidity();
274
275   if (fMaxImpactXYPtOff < 999.0) {
276     if ((fMaxImpactXYPtOff + fMaxImpactXYPtNrm*TMath::Power(tPt, fMaxImpactXYPtPow)) < TMath::Abs(track->ImpactD())) {
277       fNTracksFailed++;
278       return false;
279     }
280   }
281
282   if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
283     {
284       fNTracksFailed++;
285       //cout<<"No Go Through the cut"<<endl;
286       //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
287       return false;
288     }
289   if ((tEta<fEta[0])||(tEta>fEta[1]))
290     {
291       fNTracksFailed++;
292       //cout<<"No Go Through the cut"<<endl;
293       //cout<<fEta[0]<<" < Eta ="<<tEta<<" <"<<fEta[1]<<endl;
294       return false;
295     }
296   if ((tPt<fPt[0])||(tPt>fPt[1]))
297     {
298       fNTracksFailed++;
299       //cout<<"No Go Through the cut"<<endl;
300       //cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
301       return false;
302     }
303
304
305
306
307   //   cout << "Track has pids: "
308   //        << track->PidProbElectron() << " "
309   //        << track->PidProbMuon() << " "
310   //        << track->PidProbPion() << " "
311   //        << track->PidProbKaon() << " "
312   //        << track->PidProbProton() << " "
313   //        << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
314
315
316   if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
317     {
318       fNTracksFailed++;
319       //cout<<"No Go Through the cut"<<endl;
320       //cout<<fPidProbElectron[0]<<" < e ="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
321       return false;
322     }
323   if ((track->PidProbPion()<fPidProbPion[0])||(track->PidProbPion()>fPidProbPion[1]))
324     {
325       fNTracksFailed++;
326       //cout<<"No Go Through the cut"<<endl;
327       //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
328       return false;
329     }
330   if ((track->PidProbKaon()<fPidProbKaon[0])||(track->PidProbKaon()>fPidProbKaon[1]))
331     {
332       fNTracksFailed++;
333       //cout<<"No Go Through the cut"<<endl;
334       //cout<<fPidProbKaon[0]<<" < k ="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
335       return false;
336     }
337   if ((track->PidProbProton()<fPidProbProton[0])||(track->PidProbProton()>fPidProbProton[1]))
338     {
339       fNTracksFailed++;
340       //cout<<"No Go Through the cut"<<endl;
341       //cout<<fPidProbProton[0]<<" < p  ="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
342       return false;
343     }
344   if ((track->PidProbMuon()<fPidProbMuon[0])||(track->PidProbMuon()>fPidProbMuon[1]))
345     {
346       fNTracksFailed++;
347       //cout<<"No Go Through the cut"<<endl;
348       //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
349       return false;
350     }
351
352   if (fMostProbable) {
353
354     int imost=0;
355     tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().Mag());
356     tMost[1] = 0.0;
357     tMost[2] = track->PidProbPion()*PidFractionPion(track->P().Mag());
358     tMost[3] = track->PidProbKaon()*PidFractionKaon(track->P().Mag());
359     tMost[4] = track->PidProbProton()*PidFractionProton(track->P().Mag());
360     float ipidmax = 0.0;
361
362
363     //****N Sigma Method****
364         if(fPIDMethod==0){
365           // Looking for pions
366           if (fMostProbable == 2) {
367             if (IsPionNSigma(track->P().Mag(), track->NSigmaTPCPi(), track->NSigmaTOFPi()))
368               imost = 2;
369
370           }
371           else if (fMostProbable == 3) {
372
373
374           if (IsKaonNSigma(track->P().Mag(), track->NSigmaTPCK(), track->NSigmaTOFK())){
375
376               imost = 3;
377             }
378
379           }
380       else if (fMostProbable == 4) { // proton nsigma-PID required contour adjusting (in LHC10h)
381         if ( IsProtonNSigma(track->P().Mag(), track->NSigmaTPCP(), track->NSigmaTOFP()) && (TMath::Abs(track->NSigmaTPCP()) < TMath::Abs(track->NSigmaTPCPi())) && (TMath::Abs(track->NSigmaTPCP()) < TMath::Abs(track->NSigmaTPCK())) && (TMath::Abs(track->NSigmaTOFP()) < TMath::Abs(track->NSigmaTOFPi())) && (TMath::Abs(track->NSigmaTOFP()) < TMath::Abs(track->NSigmaTOFK()))
382              // && IsProtonTPCdEdx(track->P().Mag(), track->TPCsignal())
383             )
384               imost = 4;
385           }
386
387         }
388
389
390
391     //****Contour Method****
392         if(fPIDMethod==1){
393           for (int ip=0; ip<5; ip++)
394             if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; };
395
396           // Looking for pions
397           if (fMostProbable == 2) {
398             if (imost == 2) {
399               // Using the TPC to reject non-pions
400               if (!(IsPionTPCdEdx(track->P().Mag(), track->TPCsignal())))
401                 imost = 0;
402               if (0) {
403                 // Using the TOF to reject non-pions
404                 if (track->P().Mag() < 0.6) {
405                   if (tTOFPidIn)
406                     if (!IsPionTOFTime(track->P().Mag(), track->TOFpionTime()))
407                       imost = 0;
408                 }
409                 else {
410                   if (tTOFPidIn) {
411                     if (!IsPionTOFTime(track->P().Mag(), track->TOFpionTime()))
412                       imost = 0;
413                   }
414                   else {
415                     imost = 0;
416                   }
417                 }
418               }
419             }
420           }
421
422           // Looking for kaons
423           else if (fMostProbable == 3) {
424             //       if (imost == 3) {
425             // Using the TPC to reject non-kaons
426             if (track->P().Mag() < 0.6) {
427               if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
428                 imost = 0;
429               else imost = 3;
430               if (1) {
431                 // Using the TOF to reject non-kaons
432                 if (tTOFPidIn)
433                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
434                     imost = 0;
435               }
436             }
437             else {
438               if (1) {
439                 if (tTOFPidIn) {
440                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
441                     imost = 0;
442                   else
443                     imost = 3;
444                 }
445                 else {
446                   if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
447                     imost = 0;
448                   else
449                     imost = 3;
450                 }
451               }
452             }
453             //       }
454           }
455
456           // Looking for protons
457           else if (fMostProbable == 4) {
458             //       if (imost == 3) {
459             // Using the TPC to reject non-kaons
460             if (track->P().Mag() < 0.8) {
461               if (!(IsProtonTPCdEdx(track->P().Mag(), track->TPCsignal())))
462                 imost = 0;
463               else imost = 4;
464               if (0) {
465                 // Using the TOF to reject non-kaons
466                 if (tTOFPidIn)
467                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
468                     imost = 0;
469               }
470             }
471             else {
472               if (0) {
473                 if (tTOFPidIn) {
474                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
475                     imost = 0;
476                   else
477                     imost = 3;
478                 }
479                 else {
480                   if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
481                     imost = 0;
482                   else
483                     imost = 3;
484                 }
485               }
486             }
487             //       }
488           }
489         }
490     if (imost != fMostProbable) return false;
491   }
492
493   //fan
494   //cout<<"****** Go Through the cut ******"<<endl;
495   // cout<<fLabel<<" Label="<<track->Label()<<endl;
496   // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
497   // cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
498   //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
499   //cout<<fPidProbElectron[0]<<" <  e="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
500   //cout<<fPidProbPion[0]<<" <  pi="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
501   //cout<<fPidProbKaon[0]<<" <  k="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
502   //cout<<fPidProbProton[0]<<" <  p="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
503   //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
504   fNTracksPassed++ ;
505   return true;
506
507
508 }
509 //------------------------------
510 AliFemtoString AliFemtoESDTrackCut::Report()
511 {
512   // Prepare report from the execution
513   string tStemp;
514   char tCtemp[100];
515   snprintf(tCtemp , 100, "Particle mass:\t%E\n",this->Mass());
516   tStemp=tCtemp;
517   snprintf(tCtemp , 100, "Particle charge:\t%d\n",fCharge);
518   tStemp+=tCtemp;
519   snprintf(tCtemp , 100, "Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
520   tStemp+=tCtemp;
521   snprintf(tCtemp , 100, "Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
522   tStemp+=tCtemp;
523   snprintf(tCtemp , 100, "Particle eta:\t%E - %E\n",fEta[0],fEta[1]);
524   tStemp+=tCtemp;
525   snprintf(tCtemp , 100, "Number of tracks which passed:\t%ld  Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
526   tStemp += tCtemp;
527   AliFemtoString returnThis = tStemp;
528   return returnThis;
529 }
530 TList *AliFemtoESDTrackCut::ListSettings()
531 {
532   // return a list of settings in a writable form
533   TList *tListSetttings = new TList();
534   char buf[200];
535   snprintf(buf, 200, "AliFemtoESDTrackCut.mass=%f", this->Mass());
536   tListSetttings->AddLast(new TObjString(buf));
537
538   snprintf(buf, 200, "AliFemtoESDTrackCut.charge=%i", fCharge);
539   tListSetttings->AddLast(new TObjString(buf));
540   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.minimum=%f", fPidProbPion[0]);
541   tListSetttings->AddLast(new TObjString(buf));
542   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.maximum=%f", fPidProbPion[1]);
543   tListSetttings->AddLast(new TObjString(buf));
544   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.minimum=%f", fPidProbKaon[0]);
545   tListSetttings->AddLast(new TObjString(buf));
546   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.maximum=%f", fPidProbKaon[1]);
547   tListSetttings->AddLast(new TObjString(buf));
548   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.minimum=%f", fPidProbProton[0]);
549   tListSetttings->AddLast(new TObjString(buf));
550   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.maximum=%f", fPidProbProton[1]);
551   tListSetttings->AddLast(new TObjString(buf));
552   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.minimum=%f", fPidProbElectron[0]);
553   tListSetttings->AddLast(new TObjString(buf));
554   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.maximum=%f", fPidProbElectron[1]);
555   tListSetttings->AddLast(new TObjString(buf));
556   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.minimum=%f", fPidProbMuon[0]);
557   tListSetttings->AddLast(new TObjString(buf));
558   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.maximum=%f", fPidProbMuon[1]);
559   tListSetttings->AddLast(new TObjString(buf));
560   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumtpcclusters=%i", fminTPCclsF);
561   tListSetttings->AddLast(new TObjString(buf));
562   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumitsclusters=%i", fminTPCclsF);
563   tListSetttings->AddLast(new TObjString(buf));
564   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.minimum=%f", fPt[0]);
565   tListSetttings->AddLast(new TObjString(buf));
566   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.maximum=%f", fPt[1]);
567   tListSetttings->AddLast(new TObjString(buf));
568   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.minimum=%f", fRapidity[0]);
569   tListSetttings->AddLast(new TObjString(buf));
570   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.maximum=%f", fRapidity[1]);
571   tListSetttings->AddLast(new TObjString(buf));
572   snprintf(buf, 200, "AliFemtoESDTrackCut.removekinks=%i", fRemoveKinks);
573   tListSetttings->AddLast(new TObjString(buf));
574   snprintf(buf, 200, "AliFemtoESDTrackCut.maxitschindof=%f", fMaxITSchiNdof);
575   tListSetttings->AddLast(new TObjString(buf));
576   snprintf(buf, 200, "AliFemtoESDTrackCut.maxtpcchindof=%f", fMaxTPCchiNdof);
577   tListSetttings->AddLast(new TObjString(buf));
578   snprintf(buf, 200, "AliFemtoESDTrackCut.maxsigmatovertex=%f", fMaxSigmaToVertex);
579   tListSetttings->AddLast(new TObjString(buf));
580   snprintf(buf, 200, "AliFemtoESDTrackCut.maximpactxy=%f", fMaxImpactXY);
581   tListSetttings->AddLast(new TObjString(buf));
582   snprintf(buf, 200, "AliFemtoESDTrackCut.maximpactz=%f", fMaxImpactZ);
583   tListSetttings->AddLast(new TObjString(buf));
584   if (fMostProbable) {
585     if (fMostProbable == 2)
586       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Pion");
587     if (fMostProbable == 3)
588       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Kaon");
589     if (fMostProbable == 4)
590       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Proton");
591     tListSetttings->AddLast(new TObjString(buf));
592   }
593   return tListSetttings;
594 }
595 void AliFemtoESDTrackCut::SetRemoveKinks(const bool& flag)
596 {
597   fRemoveKinks = flag;
598 }
599
600 void AliFemtoESDTrackCut::SetRemoveITSFake(const bool& flag)
601 {
602   fRemoveITSFake = flag;
603 }
604
605 // electron
606 // 0.13 - 1.8
607 // 0       7.594129e-02    8.256141e-03
608 // 1       -5.535827e-01   8.170825e-02
609 // 2       1.728591e+00    3.104210e-01
610 // 3       -2.827893e+00   5.827802e-01
611 // 4       2.503553e+00    5.736207e-01
612 // 5       -1.125965e+00   2.821170e-01
613 // 6       2.009036e-01    5.438876e-02
614 float AliFemtoESDTrackCut::PidFractionElectron(float mom) const
615 {
616   // Provide a parameterized fraction of electrons dependent on momentum
617   if (mom<0.13)
618     return (7.594129e-02
619             -5.535827e-01*0.13
620             +1.728591e+00*0.13*0.13
621             -2.827893e+00*0.13*0.13*0.13
622             +2.503553e+00*0.13*0.13*0.13*0.13
623             -1.125965e+00*0.13*0.13*0.13*0.13*0.13
624             +2.009036e-01*0.13*0.13*0.13*0.13*0.13*0.13);
625
626   if (mom>1.8)
627     return (7.594129e-02
628             -5.535827e-01*1.8
629             +1.728591e+00*1.8*1.8
630             -2.827893e+00*1.8*1.8*1.8
631             +2.503553e+00*1.8*1.8*1.8*1.8
632             -1.125965e+00*1.8*1.8*1.8*1.8*1.8
633             +2.009036e-01*1.8*1.8*1.8*1.8*1.8*1.8);
634   return (7.594129e-02
635           -5.535827e-01*mom
636           +1.728591e+00*mom*mom
637           -2.827893e+00*mom*mom*mom
638           +2.503553e+00*mom*mom*mom*mom
639           -1.125965e+00*mom*mom*mom*mom*mom
640           +2.009036e-01*mom*mom*mom*mom*mom*mom);
641 }
642
643 // pion
644 // 0.13 - 2.0
645 // 0       1.063457e+00    8.872043e-03
646 // 1       -4.222208e-01   2.534402e-02
647 // 2       1.042004e-01    1.503945e-02
648 float AliFemtoESDTrackCut::PidFractionPion(float mom) const
649 {
650   // Provide a parameterized fraction of pions dependent on momentum
651   if (mom<0.13)
652     return ( 1.063457e+00
653              -4.222208e-01*0.13
654              +1.042004e-01*0.0169);
655   if (mom>2.0)
656     return ( 1.063457e+00
657              -4.222208e-01*2.0
658              +1.042004e-01*4.0);
659   return ( 1.063457e+00
660            -4.222208e-01*mom
661            +1.042004e-01*mom*mom);
662 }
663
664 // kaon
665 // 0.18 - 2.0
666 // 0       -7.289406e-02   1.686074e-03
667 // 1       4.415666e-01    1.143939e-02
668 // 2       -2.996790e-01   1.840964e-02
669 // 3       6.704652e-02    7.783990e-03
670 float AliFemtoESDTrackCut::PidFractionKaon(float mom) const
671 {
672   // Provide a parameterized fraction of kaons dependent on momentum
673   if (mom<0.18)
674     return (-7.289406e-02
675             +4.415666e-01*0.18
676             -2.996790e-01*0.18*0.18
677             +6.704652e-02*0.18*0.18*0.18);
678   if (mom>2.0)
679     return (-7.289406e-02
680             +4.415666e-01*2.0
681             -2.996790e-01*2.0*2.0
682             +6.704652e-02*2.0*2.0*2.0);
683   return (-7.289406e-02
684           +4.415666e-01*mom
685           -2.996790e-01*mom*mom
686           +6.704652e-02*mom*mom*mom);
687 }
688
689 // proton
690 // 0.26 - 2.0
691 // 0       -3.730200e-02   2.347311e-03
692 // 1       1.163684e-01    1.319316e-02
693 // 2       8.354116e-02    1.997948e-02
694 // 3       -4.608098e-02   8.336400e-03
695 float AliFemtoESDTrackCut::PidFractionProton(float mom) const
696 {
697   // Provide a parameterized fraction of protons dependent on momentum
698   if (mom<0.26) return  0.0;
699   if (mom>2.0)
700     return (-3.730200e-02
701             +1.163684e-01*2.0
702             +8.354116e-02*2.0*2.0
703             -4.608098e-02*2.0*2.0*2.0);
704   return (-3.730200e-02
705           +1.163684e-01*mom
706           +8.354116e-02*mom*mom
707           -4.608098e-02*mom*mom*mom);
708 }
709
710 void AliFemtoESDTrackCut::SetMomRangeTOFpidIs(const float& minp, const float& maxp)
711 {
712   fMinPforTOFpid = minp;
713   fMaxPforTOFpid = maxp;
714 }
715
716 void AliFemtoESDTrackCut::SetMomRangeTPCpidIs(const float& minp, const float& maxp)
717 {
718   fMinPforTPCpid = minp;
719   fMaxPforTPCpid = maxp;
720 }
721
722 void AliFemtoESDTrackCut::SetMomRangeITSpidIs(const float& minp, const float& maxp)
723 {
724   fMinPforITSpid = minp;
725   fMaxPforITSpid = maxp;
726 }
727
728 bool AliFemtoESDTrackCut::IsPionTPCdEdx(float mom, float dEdx)
729 {
730   //   double a1 = -95.4545, b1 = 86.5455;
731   //   double a2 = 0.0,      b2 = 56.0;
732   double a1 = -343.75,  b1 = 168.125;
733   double a2 = 0.0,      b2 = 65.0;
734
735   if (mom < 0.32) {
736     if (dEdx < a1*mom+b1) return true;
737   }
738   if (dEdx < a2*mom+b2) return true;
739
740   return false;
741 }
742
743 bool AliFemtoESDTrackCut::IsKaonTPCdEdx(float mom, float dEdx)
744 {
745
746 //   double a1 = -547.0; double b1 =  297.0;
747 //   double a2 = -125.0; double b2 =  145.0;
748 //   double a3 = -420.0; double b3 =  357.0;
749 //   double a4 = -110.0; double b4 =  171.0;
750 //   double b5 =   72.0;
751
752 //   if (mom<0.2) return false;
753
754 //   if (mom<0.36) {
755 //     if (dEdx < a1*mom+b1) return false;
756 //     if (dEdx > a3*mom+b3) return false;
757 //   }
758 //   else if (mom<0.6) {
759 //     if (dEdx < a2*mom+b2) return false;
760 //     if (dEdx > a3*mom+b3) return false;
761 //   }
762 //   else if (mom<0.9) {
763 //     if (dEdx > a4*mom+b4) return false;
764 //     if (dEdx <        b5) return false;
765 //   }
766 //   else
767 //     return false;
768 //   //   else {
769 //   //     if (dEdx > b5) return false;
770 //   //   }
771
772 //   return true;
773
774   double a1 = -268.896; double b1 =  198.669;
775   double a2 = -49.0012;  double b2 =  88.7214;
776
777   if (mom<0.2) return false;
778
779   if (mom>0.3 && mom<0.5) {
780     if (dEdx < a1*mom+b1) return false;
781   }
782   else  if (mom<1.2) {
783     if (dEdx < a2*mom+b2) return false;
784   }
785
786   return true;
787
788 }
789
790 bool AliFemtoESDTrackCut::IsProtonTPCdEdx(float mom, float dEdx)
791 {
792   double a1 = -1800.0; double b1 =  940.0;
793   double a2 = -500.0;  double b2 =  420.0;
794   double a3 = -216.7;  double b3 =  250.0;
795
796   if (mom<0.2) return false;
797
798   if (mom>0.3 && mom<0.4) {
799     if (dEdx < a1*mom+b1) return false;
800   }
801   else  if (mom<0.6) {
802     if (dEdx < a2*mom+b2) return false;
803   }
804   else  if (mom<0.9) {
805     if (dEdx < a3*mom+b3) return false;
806   }
807
808   return true;
809
810 }
811
812 bool AliFemtoESDTrackCut::IsPionTOFTime(float mom, float ttof)
813 {
814   double a1 = -427.0; double b1 =  916.0;
815   double a2 =  327.0; double b2 = -888.0;
816   if (mom<0.3) return kFALSE;
817   if (mom>2.0) return kFALSE;
818   if (ttof > a1*mom+b1) return kFALSE;
819   if (ttof < a2*mom+b2) return kFALSE;
820
821   return kTRUE;
822 }
823
824 bool AliFemtoESDTrackCut::IsKaonTOFTime(float mom, float ttof)
825 {
826   double a1 =   000.0; double b1 =  -500.0;
827   double a2 =   000.0; double b2 =   500.0;
828   double a3 =   850.0; double b3 = -1503.0;
829   double a4 = -1637.0; double b4 =  3621.0;
830
831   if (mom<0.3) return kFALSE;
832   if (mom>2.06) return kFALSE;
833   if (mom<1.2) {
834     if (ttof > a2*mom+b2) return kFALSE;
835     if (ttof < a1*mom+b1) return kFALSE;
836   }
837   if (mom<1.9) {
838     if (ttof > a2*mom+b2) return kFALSE;
839     if (ttof < a3*mom+b3) return kFALSE;
840   }
841   if (mom<2.06) {
842     if (ttof > a4*mom+b4) return kFALSE;
843     if (ttof < a3*mom+b3) return kFALSE;
844   }
845   return kTRUE;
846 }
847
848 bool AliFemtoESDTrackCut::IsProtonTOFTime(float mom, float ttof)
849 {
850   double a1 =   000.0; double b1 =  -915.0;
851   double a2 =   000.0; double b2 =   600.0;
852   double a3 =   572.0; double b3 = -1715.0;
853
854   if (mom<0.3) return kFALSE;
855   if (mom>3.0) return kFALSE;
856   if (mom<1.4) {
857     if (ttof > a2*mom+b2) return kFALSE;
858     if (ttof < a1*mom+b1) return kFALSE;
859   }
860   if (mom<3.0) {
861     if (ttof > a2*mom+b2) return kFALSE;
862     if (ttof < a3*mom+b3) return kFALSE;
863   }
864   return kTRUE;
865 }
866
867
868
869
870 bool AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma(float mom, float nsigmaK)
871 {
872 //  cout<<" AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma "<<mom<<" "<<nsigmaK<<endl;
873
874
875   if(mom<0.35 && TMath::Abs(nsigmaK)<5.0)return true;
876   if(mom>=0.35 && mom<0.5 && TMath::Abs(nsigmaK)<3.0)return true;
877   if(mom>=0.5 && mom<0.7 && TMath::Abs(nsigmaK)<2.0)return true;
878
879   return false;
880 }
881
882 bool AliFemtoESDTrackCut::IsKaonTOFNSigma(float mom, float nsigmaK)
883 {
884 //  cout<<" AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma "<<mom<<" "<<nsigmaK<<endl;
885   //fan
886   //  if(mom<1.5 && TMath::Abs(nsigmaK)<3.0)return true;
887   if(mom>=1.5 && TMath::Abs(nsigmaK)<2.0)return true;
888   return false;
889 }
890
891 /*
892 bool AliFemtoESDTrackCut::IsKaonNSigma(float mom, float nsigmaTPCK, float nsigmaTOFK)
893 {
894
895
896   if(mom<0.5)
897     {
898           if(TMath::Abs(nsigmaTPCK)<2.0)
899            {
900            return true;
901            }
902            else
903            {
904            return false;
905            }
906     }
907
908
909    if(mom>=0.5)
910     {
911          if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0)
912          {
913          return true;
914          }
915          else
916          {
917          return false;
918          }
919     }
920
921 //   if(mom>1.5 || mom<0.15)return false;
922
923
924 }
925
926 */
927
928 //old
929 bool AliFemtoESDTrackCut::IsKaonNSigma(float mom, float nsigmaTPCK, float nsigmaTOFK)
930 {
931
932   if(mom<0.4)
933     {
934       if(nsigmaTOFK<-999.)
935         {
936           if(TMath::Abs(nsigmaTPCK)<2.0) return true;
937         }
938       else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
939     }
940   else if(mom>=0.4 && mom<=0.6)
941     {
942       if(nsigmaTOFK<-999.)
943         {
944           if(TMath::Abs(nsigmaTPCK)<2.0) return true;
945         }
946       else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
947     }
948   else if(nsigmaTOFK<-999.)
949     {
950       return false;
951     }
952   else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
953
954   return false;
955 }
956
957
958
959 bool AliFemtoESDTrackCut::IsPionNSigma(float mom, float nsigmaTPCPi, float nsigmaTOFPi)
960 {
961   if(mom<0.65)
962     {
963       if(nsigmaTOFPi<-999.)
964         {
965           if(mom<0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
966           else if(mom<0.5 && mom>=0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
967           else if(mom>=0.5 && TMath::Abs(nsigmaTPCPi)<2.0) return true;
968           else return false;
969         }
970       else if(TMath::Abs(nsigmaTOFPi)<3.0 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
971     }
972   else if(nsigmaTOFPi<-999.)
973     {
974       return false;
975     }
976   else if(mom<1.5 && TMath::Abs(nsigmaTOFPi)<3.0 && TMath::Abs(nsigmaTPCPi)<5.0) return true;
977   else if(mom>=1.5 && TMath::Abs(nsigmaTOFPi)<2.0 && TMath::Abs(nsigmaTPCPi)<5.0) return true;
978
979   return false;
980 }
981
982
983 bool AliFemtoESDTrackCut::IsProtonNSigma(float mom, float nsigmaTPCP, float nsigmaTOFP)
984 {
985
986   if (fNsigmaTPCTOF) {
987     if (mom > 0.8) {
988 //        if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP )/TMath::Sqrt(2) < 3.0)
989         if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP ) < fNsigma)
990             return true;
991         }
992     else {
993         if (TMath::Abs(nsigmaTPCP) < fNsigma)
994             return true;
995     }
996   }
997   else if (fNsigmaTPConly) {
998     if (TMath::Abs(nsigmaTPCP) < fNsigma)
999       return true;
1000   }
1001   else {
1002     if (mom > 0.8 && mom < 2.5) {
1003       if ( TMath::Abs(nsigmaTPCP) < 3.0 && TMath::Abs(nsigmaTOFP) < 3.0)
1004         return true;
1005     }
1006     else if (mom > 2.5) {
1007       if ( TMath::Abs(nsigmaTPCP) < 3.0 && TMath::Abs(nsigmaTOFP) < 2.0)
1008         return true;
1009     }
1010     else {
1011       if (TMath::Abs(nsigmaTPCP) < 3.0)
1012         return true;
1013     }
1014   }
1015
1016   return false;
1017 }
1018
1019
1020 void AliFemtoESDTrackCut::SetPIDMethod(ReadPIDMethodType newMethod)
1021 {
1022   fPIDMethod = newMethod;
1023 }
1024
1025
1026 void AliFemtoESDTrackCut::SetNsigmaTPCTOF(Bool_t nsigma)
1027 {
1028   fNsigmaTPCTOF = nsigma;
1029 }
1030
1031 void AliFemtoESDTrackCut::SetNsigmaTPConly(Bool_t nsigma)
1032 {
1033   fNsigmaTPConly = nsigma;
1034 }
1035
1036 void AliFemtoESDTrackCut::SetNsigma(Double_t nsigma)
1037 {
1038   fNsigma = nsigma;
1039 }
1040
1041 void AliFemtoESDTrackCut::SetClusterRequirementITS(AliESDtrackCuts::Detector det, AliESDtrackCuts::ITSClusterRequirement req)
1042 {
1043   fCutClusterRequirementITS[det] = req;
1044 }
1045
1046 Bool_t AliFemtoESDTrackCut::CheckITSClusterRequirement(AliESDtrackCuts::ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1047 {
1048   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1049
1050   switch (req)
1051     {
1052     case AliESDtrackCuts::kOff:        return kTRUE;
1053     case AliESDtrackCuts::kNone:       return !clusterL1 && !clusterL2;
1054     case AliESDtrackCuts::kAny:        return clusterL1 || clusterL2;
1055     case AliESDtrackCuts::kFirst:      return clusterL1;
1056     case AliESDtrackCuts::kOnlyFirst:  return clusterL1 && !clusterL2;
1057     case AliESDtrackCuts::kSecond:     return clusterL2;
1058     case AliESDtrackCuts::kOnlySecond: return clusterL2 && !clusterL1;
1059     case AliESDtrackCuts::kBoth:       return clusterL1 && clusterL2;
1060   }
1061
1062   return kFALSE;
1063 }