changes in the PID of protons + new config for the proton femtoscopy train
[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     else if (fMostProbable == 5) { // no-protons
387       if ( !IsProtonNSigma(track->P().Mag(), track->NSigmaTPCP(), track->NSigmaTOFP()) )
388               imost = 5;
389           }
390
391         }
392
393
394
395     //****Contour Method****
396         if(fPIDMethod==1){
397           for (int ip=0; ip<5; ip++)
398             if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; };
399
400           // Looking for pions
401           if (fMostProbable == 2) {
402             if (imost == 2) {
403               // Using the TPC to reject non-pions
404               if (!(IsPionTPCdEdx(track->P().Mag(), track->TPCsignal())))
405                 imost = 0;
406               if (0) {
407                 // Using the TOF to reject non-pions
408                 if (track->P().Mag() < 0.6) {
409                   if (tTOFPidIn)
410                     if (!IsPionTOFTime(track->P().Mag(), track->TOFpionTime()))
411                       imost = 0;
412                 }
413                 else {
414                   if (tTOFPidIn) {
415                     if (!IsPionTOFTime(track->P().Mag(), track->TOFpionTime()))
416                       imost = 0;
417                   }
418                   else {
419                     imost = 0;
420                   }
421                 }
422               }
423             }
424           }
425
426           // Looking for kaons
427           else if (fMostProbable == 3) {
428             //       if (imost == 3) {
429             // Using the TPC to reject non-kaons
430             if (track->P().Mag() < 0.6) {
431               if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
432                 imost = 0;
433               else imost = 3;
434               if (1) {
435                 // Using the TOF to reject non-kaons
436                 if (tTOFPidIn)
437                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
438                     imost = 0;
439               }
440             }
441             else {
442               if (1) {
443                 if (tTOFPidIn) {
444                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
445                     imost = 0;
446                   else
447                     imost = 3;
448                 }
449                 else {
450                   if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
451                     imost = 0;
452                   else
453                     imost = 3;
454                 }
455               }
456             }
457             //       }
458           }
459
460           // Looking for protons
461           else if (fMostProbable == 4) {
462             //       if (imost == 3) {
463             // Using the TPC to reject non-kaons
464             if (track->P().Mag() < 0.8) {
465               if (!(IsProtonTPCdEdx(track->P().Mag(), track->TPCsignal())))
466                 imost = 0;
467               else imost = 4;
468               if (0) {
469                 // Using the TOF to reject non-kaons
470                 if (tTOFPidIn)
471                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
472                     imost = 0;
473               }
474             }
475             else {
476               if (0) {
477                 if (tTOFPidIn) {
478                   if (!IsKaonTOFTime(track->P().Mag(), track->TOFkaonTime()))
479                     imost = 0;
480                   else
481                     imost = 3;
482                 }
483                 else {
484                   if (!(IsKaonTPCdEdx(track->P().Mag(), track->TPCsignal())))
485                     imost = 0;
486                   else
487                     imost = 3;
488                 }
489               }
490             }
491             //       }
492           }
493         }
494     if (imost != fMostProbable) return false;
495   }
496
497   //fan
498   //cout<<"****** Go Through the cut ******"<<endl;
499   // cout<<fLabel<<" Label="<<track->Label()<<endl;
500   // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
501   // cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
502   //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
503   //cout<<fPidProbElectron[0]<<" <  e="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
504   //cout<<fPidProbPion[0]<<" <  pi="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
505   //cout<<fPidProbKaon[0]<<" <  k="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
506   //cout<<fPidProbProton[0]<<" <  p="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
507   //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
508   fNTracksPassed++ ;
509   return true;
510
511
512 }
513 //------------------------------
514 AliFemtoString AliFemtoESDTrackCut::Report()
515 {
516   // Prepare report from the execution
517   string tStemp;
518   char tCtemp[100];
519   snprintf(tCtemp , 100, "Particle mass:\t%E\n",this->Mass());
520   tStemp=tCtemp;
521   snprintf(tCtemp , 100, "Particle charge:\t%d\n",fCharge);
522   tStemp+=tCtemp;
523   snprintf(tCtemp , 100, "Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
524   tStemp+=tCtemp;
525   snprintf(tCtemp , 100, "Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
526   tStemp+=tCtemp;
527   snprintf(tCtemp , 100, "Particle eta:\t%E - %E\n",fEta[0],fEta[1]);
528   tStemp+=tCtemp;
529   snprintf(tCtemp , 100, "Number of tracks which passed:\t%ld  Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
530   tStemp += tCtemp;
531   AliFemtoString returnThis = tStemp;
532   return returnThis;
533 }
534 TList *AliFemtoESDTrackCut::ListSettings()
535 {
536   // return a list of settings in a writable form
537   TList *tListSetttings = new TList();
538   char buf[200];
539   snprintf(buf, 200, "AliFemtoESDTrackCut.mass=%f", this->Mass());
540   tListSetttings->AddLast(new TObjString(buf));
541
542   snprintf(buf, 200, "AliFemtoESDTrackCut.charge=%i", fCharge);
543   tListSetttings->AddLast(new TObjString(buf));
544   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.minimum=%f", fPidProbPion[0]);
545   tListSetttings->AddLast(new TObjString(buf));
546   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.maximum=%f", fPidProbPion[1]);
547   tListSetttings->AddLast(new TObjString(buf));
548   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.minimum=%f", fPidProbKaon[0]);
549   tListSetttings->AddLast(new TObjString(buf));
550   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.maximum=%f", fPidProbKaon[1]);
551   tListSetttings->AddLast(new TObjString(buf));
552   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.minimum=%f", fPidProbProton[0]);
553   tListSetttings->AddLast(new TObjString(buf));
554   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.maximum=%f", fPidProbProton[1]);
555   tListSetttings->AddLast(new TObjString(buf));
556   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.minimum=%f", fPidProbElectron[0]);
557   tListSetttings->AddLast(new TObjString(buf));
558   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.maximum=%f", fPidProbElectron[1]);
559   tListSetttings->AddLast(new TObjString(buf));
560   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.minimum=%f", fPidProbMuon[0]);
561   tListSetttings->AddLast(new TObjString(buf));
562   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.maximum=%f", fPidProbMuon[1]);
563   tListSetttings->AddLast(new TObjString(buf));
564   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumtpcclusters=%i", fminTPCclsF);
565   tListSetttings->AddLast(new TObjString(buf));
566   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumitsclusters=%i", fminTPCclsF);
567   tListSetttings->AddLast(new TObjString(buf));
568   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.minimum=%f", fPt[0]);
569   tListSetttings->AddLast(new TObjString(buf));
570   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.maximum=%f", fPt[1]);
571   tListSetttings->AddLast(new TObjString(buf));
572   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.minimum=%f", fRapidity[0]);
573   tListSetttings->AddLast(new TObjString(buf));
574   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.maximum=%f", fRapidity[1]);
575   tListSetttings->AddLast(new TObjString(buf));
576   snprintf(buf, 200, "AliFemtoESDTrackCut.removekinks=%i", fRemoveKinks);
577   tListSetttings->AddLast(new TObjString(buf));
578   snprintf(buf, 200, "AliFemtoESDTrackCut.maxitschindof=%f", fMaxITSchiNdof);
579   tListSetttings->AddLast(new TObjString(buf));
580   snprintf(buf, 200, "AliFemtoESDTrackCut.maxtpcchindof=%f", fMaxTPCchiNdof);
581   tListSetttings->AddLast(new TObjString(buf));
582   snprintf(buf, 200, "AliFemtoESDTrackCut.maxsigmatovertex=%f", fMaxSigmaToVertex);
583   tListSetttings->AddLast(new TObjString(buf));
584   snprintf(buf, 200, "AliFemtoESDTrackCut.maximpactxy=%f", fMaxImpactXY);
585   tListSetttings->AddLast(new TObjString(buf));
586   snprintf(buf, 200, "AliFemtoESDTrackCut.maximpactz=%f", fMaxImpactZ);
587   tListSetttings->AddLast(new TObjString(buf));
588   if (fMostProbable) {
589     if (fMostProbable == 2)
590       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Pion");
591     if (fMostProbable == 3)
592       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Kaon");
593     if (fMostProbable == 4)
594       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Proton");
595     tListSetttings->AddLast(new TObjString(buf));
596   }
597   return tListSetttings;
598 }
599 void AliFemtoESDTrackCut::SetRemoveKinks(const bool& flag)
600 {
601   fRemoveKinks = flag;
602 }
603
604 void AliFemtoESDTrackCut::SetRemoveITSFake(const bool& flag)
605 {
606   fRemoveITSFake = flag;
607 }
608
609 // electron
610 // 0.13 - 1.8
611 // 0       7.594129e-02    8.256141e-03
612 // 1       -5.535827e-01   8.170825e-02
613 // 2       1.728591e+00    3.104210e-01
614 // 3       -2.827893e+00   5.827802e-01
615 // 4       2.503553e+00    5.736207e-01
616 // 5       -1.125965e+00   2.821170e-01
617 // 6       2.009036e-01    5.438876e-02
618 float AliFemtoESDTrackCut::PidFractionElectron(float mom) const
619 {
620   // Provide a parameterized fraction of electrons dependent on momentum
621   if (mom<0.13)
622     return (7.594129e-02
623             -5.535827e-01*0.13
624             +1.728591e+00*0.13*0.13
625             -2.827893e+00*0.13*0.13*0.13
626             +2.503553e+00*0.13*0.13*0.13*0.13
627             -1.125965e+00*0.13*0.13*0.13*0.13*0.13
628             +2.009036e-01*0.13*0.13*0.13*0.13*0.13*0.13);
629
630   if (mom>1.8)
631     return (7.594129e-02
632             -5.535827e-01*1.8
633             +1.728591e+00*1.8*1.8
634             -2.827893e+00*1.8*1.8*1.8
635             +2.503553e+00*1.8*1.8*1.8*1.8
636             -1.125965e+00*1.8*1.8*1.8*1.8*1.8
637             +2.009036e-01*1.8*1.8*1.8*1.8*1.8*1.8);
638   return (7.594129e-02
639           -5.535827e-01*mom
640           +1.728591e+00*mom*mom
641           -2.827893e+00*mom*mom*mom
642           +2.503553e+00*mom*mom*mom*mom
643           -1.125965e+00*mom*mom*mom*mom*mom
644           +2.009036e-01*mom*mom*mom*mom*mom*mom);
645 }
646
647 // pion
648 // 0.13 - 2.0
649 // 0       1.063457e+00    8.872043e-03
650 // 1       -4.222208e-01   2.534402e-02
651 // 2       1.042004e-01    1.503945e-02
652 float AliFemtoESDTrackCut::PidFractionPion(float mom) const
653 {
654   // Provide a parameterized fraction of pions dependent on momentum
655   if (mom<0.13)
656     return ( 1.063457e+00
657              -4.222208e-01*0.13
658              +1.042004e-01*0.0169);
659   if (mom>2.0)
660     return ( 1.063457e+00
661              -4.222208e-01*2.0
662              +1.042004e-01*4.0);
663   return ( 1.063457e+00
664            -4.222208e-01*mom
665            +1.042004e-01*mom*mom);
666 }
667
668 // kaon
669 // 0.18 - 2.0
670 // 0       -7.289406e-02   1.686074e-03
671 // 1       4.415666e-01    1.143939e-02
672 // 2       -2.996790e-01   1.840964e-02
673 // 3       6.704652e-02    7.783990e-03
674 float AliFemtoESDTrackCut::PidFractionKaon(float mom) const
675 {
676   // Provide a parameterized fraction of kaons dependent on momentum
677   if (mom<0.18)
678     return (-7.289406e-02
679             +4.415666e-01*0.18
680             -2.996790e-01*0.18*0.18
681             +6.704652e-02*0.18*0.18*0.18);
682   if (mom>2.0)
683     return (-7.289406e-02
684             +4.415666e-01*2.0
685             -2.996790e-01*2.0*2.0
686             +6.704652e-02*2.0*2.0*2.0);
687   return (-7.289406e-02
688           +4.415666e-01*mom
689           -2.996790e-01*mom*mom
690           +6.704652e-02*mom*mom*mom);
691 }
692
693 // proton
694 // 0.26 - 2.0
695 // 0       -3.730200e-02   2.347311e-03
696 // 1       1.163684e-01    1.319316e-02
697 // 2       8.354116e-02    1.997948e-02
698 // 3       -4.608098e-02   8.336400e-03
699 float AliFemtoESDTrackCut::PidFractionProton(float mom) const
700 {
701   // Provide a parameterized fraction of protons dependent on momentum
702   if (mom<0.26) return  0.0;
703   if (mom>2.0)
704     return (-3.730200e-02
705             +1.163684e-01*2.0
706             +8.354116e-02*2.0*2.0
707             -4.608098e-02*2.0*2.0*2.0);
708   return (-3.730200e-02
709           +1.163684e-01*mom
710           +8.354116e-02*mom*mom
711           -4.608098e-02*mom*mom*mom);
712 }
713
714 void AliFemtoESDTrackCut::SetMomRangeTOFpidIs(const float& minp, const float& maxp)
715 {
716   fMinPforTOFpid = minp;
717   fMaxPforTOFpid = maxp;
718 }
719
720 void AliFemtoESDTrackCut::SetMomRangeTPCpidIs(const float& minp, const float& maxp)
721 {
722   fMinPforTPCpid = minp;
723   fMaxPforTPCpid = maxp;
724 }
725
726 void AliFemtoESDTrackCut::SetMomRangeITSpidIs(const float& minp, const float& maxp)
727 {
728   fMinPforITSpid = minp;
729   fMaxPforITSpid = maxp;
730 }
731
732 bool AliFemtoESDTrackCut::IsPionTPCdEdx(float mom, float dEdx)
733 {
734   //   double a1 = -95.4545, b1 = 86.5455;
735   //   double a2 = 0.0,      b2 = 56.0;
736   double a1 = -343.75,  b1 = 168.125;
737   double a2 = 0.0,      b2 = 65.0;
738
739   if (mom < 0.32) {
740     if (dEdx < a1*mom+b1) return true;
741   }
742   if (dEdx < a2*mom+b2) return true;
743
744   return false;
745 }
746
747 bool AliFemtoESDTrackCut::IsKaonTPCdEdx(float mom, float dEdx)
748 {
749
750 //   double a1 = -547.0; double b1 =  297.0;
751 //   double a2 = -125.0; double b2 =  145.0;
752 //   double a3 = -420.0; double b3 =  357.0;
753 //   double a4 = -110.0; double b4 =  171.0;
754 //   double b5 =   72.0;
755
756 //   if (mom<0.2) return false;
757
758 //   if (mom<0.36) {
759 //     if (dEdx < a1*mom+b1) return false;
760 //     if (dEdx > a3*mom+b3) return false;
761 //   }
762 //   else if (mom<0.6) {
763 //     if (dEdx < a2*mom+b2) return false;
764 //     if (dEdx > a3*mom+b3) return false;
765 //   }
766 //   else if (mom<0.9) {
767 //     if (dEdx > a4*mom+b4) return false;
768 //     if (dEdx <        b5) return false;
769 //   }
770 //   else
771 //     return false;
772 //   //   else {
773 //   //     if (dEdx > b5) return false;
774 //   //   }
775
776 //   return true;
777
778   double a1 = -268.896; double b1 =  198.669;
779   double a2 = -49.0012;  double b2 =  88.7214;
780
781   if (mom<0.2) return false;
782
783   if (mom>0.3 && mom<0.5) {
784     if (dEdx < a1*mom+b1) return false;
785   }
786   else  if (mom<1.2) {
787     if (dEdx < a2*mom+b2) return false;
788   }
789
790   return true;
791
792 }
793
794 bool AliFemtoESDTrackCut::IsProtonTPCdEdx(float mom, float dEdx)
795 {
796   double a1 = -1800.0; double b1 =  940.0;
797   double a2 = -500.0;  double b2 =  420.0;
798   double a3 = -216.7;  double b3 =  250.0;
799
800   if (mom<0.2) return false;
801
802   if (mom>0.3 && mom<0.4) {
803     if (dEdx < a1*mom+b1) return false;
804   }
805   else  if (mom<0.6) {
806     if (dEdx < a2*mom+b2) return false;
807   }
808   else  if (mom<0.9) {
809     if (dEdx < a3*mom+b3) return false;
810   }
811
812   return true;
813
814 }
815
816 bool AliFemtoESDTrackCut::IsPionTOFTime(float mom, float ttof)
817 {
818   double a1 = -427.0; double b1 =  916.0;
819   double a2 =  327.0; double b2 = -888.0;
820   if (mom<0.3) return kFALSE;
821   if (mom>2.0) return kFALSE;
822   if (ttof > a1*mom+b1) return kFALSE;
823   if (ttof < a2*mom+b2) return kFALSE;
824
825   return kTRUE;
826 }
827
828 bool AliFemtoESDTrackCut::IsKaonTOFTime(float mom, float ttof)
829 {
830   double a1 =   000.0; double b1 =  -500.0;
831   double a2 =   000.0; double b2 =   500.0;
832   double a3 =   850.0; double b3 = -1503.0;
833   double a4 = -1637.0; double b4 =  3621.0;
834
835   if (mom<0.3) return kFALSE;
836   if (mom>2.06) return kFALSE;
837   if (mom<1.2) {
838     if (ttof > a2*mom+b2) return kFALSE;
839     if (ttof < a1*mom+b1) return kFALSE;
840   }
841   if (mom<1.9) {
842     if (ttof > a2*mom+b2) return kFALSE;
843     if (ttof < a3*mom+b3) return kFALSE;
844   }
845   if (mom<2.06) {
846     if (ttof > a4*mom+b4) return kFALSE;
847     if (ttof < a3*mom+b3) return kFALSE;
848   }
849   return kTRUE;
850 }
851
852 bool AliFemtoESDTrackCut::IsProtonTOFTime(float mom, float ttof)
853 {
854   double a1 =   000.0; double b1 =  -915.0;
855   double a2 =   000.0; double b2 =   600.0;
856   double a3 =   572.0; double b3 = -1715.0;
857
858   if (mom<0.3) return kFALSE;
859   if (mom>3.0) return kFALSE;
860   if (mom<1.4) {
861     if (ttof > a2*mom+b2) return kFALSE;
862     if (ttof < a1*mom+b1) return kFALSE;
863   }
864   if (mom<3.0) {
865     if (ttof > a2*mom+b2) return kFALSE;
866     if (ttof < a3*mom+b3) return kFALSE;
867   }
868   return kTRUE;
869 }
870
871
872
873
874 bool AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma(float mom, float nsigmaK)
875 {
876 //  cout<<" AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma "<<mom<<" "<<nsigmaK<<endl;
877
878
879   if(mom<0.35 && TMath::Abs(nsigmaK)<5.0)return true;
880   if(mom>=0.35 && mom<0.5 && TMath::Abs(nsigmaK)<3.0)return true;
881   if(mom>=0.5 && mom<0.7 && TMath::Abs(nsigmaK)<2.0)return true;
882
883   return false;
884 }
885
886 bool AliFemtoESDTrackCut::IsKaonTOFNSigma(float mom, float nsigmaK)
887 {
888 //  cout<<" AliFemtoESDTrackCut::IsKaonTPCdEdxNSigma "<<mom<<" "<<nsigmaK<<endl;
889   //fan
890   //  if(mom<1.5 && TMath::Abs(nsigmaK)<3.0)return true;
891   if(mom>=1.5 && TMath::Abs(nsigmaK)<2.0)return true;
892   return false;
893 }
894
895 /*
896 bool AliFemtoESDTrackCut::IsKaonNSigma(float mom, float nsigmaTPCK, float nsigmaTOFK)
897 {
898
899
900   if(mom<0.5)
901     {
902           if(TMath::Abs(nsigmaTPCK)<2.0)
903            {
904            return true;
905            }
906            else
907            {
908            return false;
909            }
910     }
911
912
913    if(mom>=0.5)
914     {
915          if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0)
916          {
917          return true;
918          }
919          else
920          {
921          return false;
922          }
923     }
924
925 //   if(mom>1.5 || mom<0.15)return false;
926
927
928 }
929
930 */
931
932 //old
933 bool AliFemtoESDTrackCut::IsKaonNSigma(float mom, float nsigmaTPCK, float nsigmaTOFK)
934 {
935
936   if(mom<0.4)
937     {
938       if(nsigmaTOFK<-999.)
939         {
940           if(TMath::Abs(nsigmaTPCK)<2.0) return true;
941         }
942       else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
943     }
944   else if(mom>=0.4 && mom<=0.6)
945     {
946       if(nsigmaTOFK<-999.)
947         {
948           if(TMath::Abs(nsigmaTPCK)<2.0) return true;
949         }
950       else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
951     }
952   else if(nsigmaTOFK<-999.)
953     {
954       return false;
955     }
956   else if(TMath::Abs(nsigmaTOFK)<3.0 && TMath::Abs(nsigmaTPCK)<3.0) return true;
957
958   return false;
959 }
960
961
962
963 bool AliFemtoESDTrackCut::IsPionNSigma(float mom, float nsigmaTPCPi, float nsigmaTOFPi)
964 {
965   if(mom<0.65)
966     {
967       if(nsigmaTOFPi<-999.)
968         {
969           if(mom<0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
970           else if(mom<0.5 && mom>=0.35 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
971           else if(mom>=0.5 && TMath::Abs(nsigmaTPCPi)<2.0) return true;
972           else return false;
973         }
974       else if(TMath::Abs(nsigmaTOFPi)<3.0 && TMath::Abs(nsigmaTPCPi)<3.0) return true;
975     }
976   else if(nsigmaTOFPi<-999.)
977     {
978       return false;
979     }
980   else if(mom<1.5 && TMath::Abs(nsigmaTOFPi)<3.0 && TMath::Abs(nsigmaTPCPi)<5.0) return true;
981   else if(mom>=1.5 && TMath::Abs(nsigmaTOFPi)<2.0 && TMath::Abs(nsigmaTPCPi)<5.0) return true;
982
983   return false;
984 }
985
986
987 bool AliFemtoESDTrackCut::IsProtonNSigma(float mom, float nsigmaTPCP, float nsigmaTOFP)
988 {
989
990   if (fNsigmaTPCTOF) {
991     if (mom > 0.8) {
992 //        if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP )/TMath::Sqrt(2) < 3.0)
993         if (TMath::Hypot( nsigmaTOFP, nsigmaTPCP ) < fNsigma)
994             return true;
995         }
996     else {
997         if (TMath::Abs(nsigmaTPCP) < fNsigma)
998             return true;
999     }
1000   }
1001   else if (fNsigmaTPConly) {
1002     if (TMath::Abs(nsigmaTPCP) < fNsigma)
1003       return true;
1004   }
1005   else {
1006     if (mom > 0.8 && mom < 2.5) {
1007       if ( TMath::Abs(nsigmaTPCP) < 3.0 && TMath::Abs(nsigmaTOFP) < 3.0)
1008         return true;
1009     }
1010     else if (mom > 2.5) {
1011       if ( TMath::Abs(nsigmaTPCP) < 3.0 && TMath::Abs(nsigmaTOFP) < 2.0)
1012         return true;
1013     }
1014     else {
1015       if (TMath::Abs(nsigmaTPCP) < 3.0)
1016         return true;
1017     }
1018   }
1019
1020   return false;
1021 }
1022
1023
1024 void AliFemtoESDTrackCut::SetPIDMethod(ReadPIDMethodType newMethod)
1025 {
1026   fPIDMethod = newMethod;
1027 }
1028
1029
1030 void AliFemtoESDTrackCut::SetNsigmaTPCTOF(Bool_t nsigma)
1031 {
1032   fNsigmaTPCTOF = nsigma;
1033 }
1034
1035 void AliFemtoESDTrackCut::SetNsigmaTPConly(Bool_t nsigma)
1036 {
1037   fNsigmaTPConly = nsigma;
1038 }
1039
1040 void AliFemtoESDTrackCut::SetNsigma(Double_t nsigma)
1041 {
1042   fNsigma = nsigma;
1043 }
1044
1045 void AliFemtoESDTrackCut::SetClusterRequirementITS(AliESDtrackCuts::Detector det, AliESDtrackCuts::ITSClusterRequirement req)
1046 {
1047   fCutClusterRequirementITS[det] = req;
1048 }
1049
1050 Bool_t AliFemtoESDTrackCut::CheckITSClusterRequirement(AliESDtrackCuts::ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
1051 {
1052   // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
1053
1054   switch (req)
1055     {
1056     case AliESDtrackCuts::kOff:        return kTRUE;
1057     case AliESDtrackCuts::kNone:       return !clusterL1 && !clusterL2;
1058     case AliESDtrackCuts::kAny:        return clusterL1 || clusterL2;
1059     case AliESDtrackCuts::kFirst:      return clusterL1;
1060     case AliESDtrackCuts::kOnlyFirst:  return clusterL1 && !clusterL2;
1061     case AliESDtrackCuts::kSecond:     return clusterL2;
1062     case AliESDtrackCuts::kOnlySecond: return clusterL2 && !clusterL1;
1063     case AliESDtrackCuts::kBoth:       return clusterL1 && clusterL2;
1064   }
1065
1066   return kFALSE;
1067 }