d2a222c87d9107fa0a7d5ee86272c468d7b8dad5
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliFlowTrackCuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */ 
17
18 // AliFlowTrackCuts:
19 // ESD track cuts for flow framework 
20 //
21 // origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
22 //
23 // This class gurantees consistency of cut methods, trackparameter
24 // selection (global tracks, TPC only, etc..) and parameter mixing
25 // in the flow framework. Transparently handles different input types:
26 // ESD, MC, AOD.
27 // This class works in 2 steps: first the requested track parameters are
28 // constructed (to be set by SetParamType() ), then cuts are applied.
29 // the constructed track can be requested AFTER checking the cuts by
30 // calling GetTrack(), in this case the cut object stays in control,
31 // caller does not have to delete the track.
32 // Additionally caller can request an AliFlowTrack object to be constructed
33 // according the parameter mixing scenario requested by SetParamMix().
34 // AliFlowTrack is made using MakeFlowTrack() method, its an 'object factory'
35 // so caller needs to take care of the freshly created object.
36
37 #include <limits.h>
38 #include <float.h>
39 #include <TMatrix.h>
40 #include "TParticle.h"
41 #include "TObjArray.h"
42 #include "AliStack.h"
43 #include "AliMCEvent.h"
44 #include "AliESDEvent.h"
45 #include "AliVParticle.h"
46 #include "AliMCParticle.h"
47 #include "AliESDtrack.h"
48 #include "AliMultiplicity.h"
49 #include "AliAODTrack.h"
50 #include "AliFlowTrack.h"
51 #include "AliFlowTrackCuts.h"
52 #include "AliLog.h"
53 #include "AliESDpid.h"
54
55 ClassImp(AliFlowTrackCuts)
56
57 //-----------------------------------------------------------------------
58 AliFlowTrackCuts::AliFlowTrackCuts():
59   AliFlowTrackSimpleCuts(),
60   fAliESDtrackCuts(NULL),
61   fQA(NULL),
62   fCutMC(kFALSE),
63   fCutMCprocessType(kFALSE),
64   fMCprocessType(kPNoProcess),
65   fCutMCPID(kFALSE),
66   fMCPID(0),
67   fIgnoreSignInPID(kFALSE),
68   fCutMCisPrimary(kFALSE),
69   fRequireTransportBitForPrimaries(kTRUE),
70   fMCisPrimary(kFALSE),
71   fRequireCharge(kFALSE),
72   fFakesAreOK(kTRUE),
73   fCutSPDtrackletDeltaPhi(kFALSE),
74   fSPDtrackletDeltaPhiMax(FLT_MAX),
75   fSPDtrackletDeltaPhiMin(-FLT_MAX),
76   fIgnoreTPCzRange(kFALSE),
77   fIgnoreTPCzRangeMax(FLT_MAX),
78   fIgnoreTPCzRangeMin(-FLT_MAX),
79   fCutChi2PerClusterTPC(kFALSE),
80   fMaxChi2PerClusterTPC(FLT_MAX),
81   fMinChi2PerClusterTPC(-FLT_MAX),
82   fCutNClustersTPC(kFALSE),
83   fNClustersTPCMax(INT_MAX),
84   fNClustersTPCMin(INT_MIN),  
85   fCutNClustersITS(kFALSE),
86   fNClustersITSMax(INT_MAX),
87   fNClustersITSMin(INT_MIN),  
88   fUseAODFilterBit(kFALSE),
89   fAODFilterBit(0),
90   fCutDCAToVertexXY(kFALSE),
91   fCutDCAToVertexZ(kFALSE),
92   fParamType(kGlobal),
93   fParamMix(kPure),
94   fTrack(NULL),
95   fTrackPhi(0.),
96   fTrackEta(0.),
97   fTrackWeight(0.),
98   fTrackLabel(INT_MIN),
99   fMCevent(NULL),
100   fMCparticle(NULL),
101   fEvent(NULL),
102   fTPCtrack(),
103   fESDpid(),
104   fPIDsource(kTPCTOFpid),
105   fTPCpidCuts(NULL),
106   fTOFpidCuts(NULL),
107   fTPCTOFpidCrossOverPt(0.4),
108   fAliPID(AliPID::kPion)
109 {
110   //io constructor 
111 }
112
113 //-----------------------------------------------------------------------
114 AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
115   AliFlowTrackSimpleCuts(),
116   fAliESDtrackCuts(new AliESDtrackCuts()),
117   fQA(NULL),
118   fCutMC(kFALSE),
119   fCutMCprocessType(kFALSE),
120   fMCprocessType(kPNoProcess),
121   fCutMCPID(kFALSE),
122   fMCPID(0),
123   fIgnoreSignInPID(kFALSE),
124   fCutMCisPrimary(kFALSE),
125   fRequireTransportBitForPrimaries(kTRUE),
126   fMCisPrimary(kFALSE),
127   fRequireCharge(kFALSE),
128   fFakesAreOK(kTRUE),
129   fCutSPDtrackletDeltaPhi(kFALSE),
130   fSPDtrackletDeltaPhiMax(FLT_MAX),
131   fSPDtrackletDeltaPhiMin(-FLT_MAX),
132   fIgnoreTPCzRange(kFALSE),
133   fIgnoreTPCzRangeMax(FLT_MAX),
134   fIgnoreTPCzRangeMin(-FLT_MAX),
135   fCutChi2PerClusterTPC(kFALSE),
136   fMaxChi2PerClusterTPC(FLT_MAX),
137   fMinChi2PerClusterTPC(-FLT_MAX),
138   fCutNClustersTPC(kFALSE),
139   fNClustersTPCMax(INT_MAX),
140   fNClustersTPCMin(INT_MIN),  
141   fCutNClustersITS(kFALSE),
142   fNClustersITSMax(INT_MAX),
143   fNClustersITSMin(INT_MIN),
144   fUseAODFilterBit(kFALSE),
145   fAODFilterBit(0),
146   fCutDCAToVertexXY(kFALSE),
147   fCutDCAToVertexZ(kFALSE),
148   fParamType(kGlobal),
149   fParamMix(kPure),
150   fTrack(NULL),
151   fTrackPhi(0.),
152   fTrackEta(0.),
153   fTrackWeight(0.),
154   fTrackLabel(INT_MIN),
155   fMCevent(NULL),
156   fMCparticle(NULL),
157   fEvent(NULL),
158   fTPCtrack(),
159   fESDpid(),
160   fPIDsource(kTPCTOFpid),
161   fTPCpidCuts(NULL),
162   fTOFpidCuts(NULL),
163   fTPCTOFpidCrossOverPt(0.4),
164   fAliPID(AliPID::kPion)
165 {
166   //constructor 
167   SetName(name);
168   SetTitle("AliFlowTrackCuts");
169   fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
170                                                     2.63394e+01,
171                                                     5.04114e-11,
172                                                     2.12543e+00,
173                                                     4.88663e+00 );
174
175 }
176
177 //-----------------------------------------------------------------------
178 AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
179   AliFlowTrackSimpleCuts(that),
180   fAliESDtrackCuts(NULL),
181   fQA(NULL),
182   fCutMC(that.fCutMC),
183   fCutMCprocessType(that.fCutMCprocessType),
184   fMCprocessType(that.fMCprocessType),
185   fCutMCPID(that.fCutMCPID),
186   fMCPID(that.fMCPID),
187   fIgnoreSignInPID(that.fIgnoreSignInPID),
188   fCutMCisPrimary(that.fCutMCisPrimary),
189   fRequireTransportBitForPrimaries(that.fRequireTransportBitForPrimaries),
190   fMCisPrimary(that.fMCisPrimary),
191   fRequireCharge(that.fRequireCharge),
192   fFakesAreOK(that.fFakesAreOK),
193   fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
194   fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
195   fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
196   fIgnoreTPCzRange(that.fIgnoreTPCzRange),
197   fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
198   fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
199   fCutChi2PerClusterTPC(that.fCutChi2PerClusterTPC),
200   fMaxChi2PerClusterTPC(that.fMaxChi2PerClusterTPC),
201   fMinChi2PerClusterTPC(that.fMinChi2PerClusterTPC),
202   fCutNClustersTPC(that.fCutNClustersTPC),
203   fNClustersTPCMax(that.fNClustersTPCMax),
204   fNClustersTPCMin(that.fNClustersTPCMin),
205   fCutNClustersITS(that.fCutNClustersITS),
206   fNClustersITSMax(that.fNClustersITSMax),
207   fNClustersITSMin(that.fNClustersITSMin),
208   fUseAODFilterBit(that.fUseAODFilterBit),
209   fAODFilterBit(that.fAODFilterBit),
210   fCutDCAToVertexXY(that.fCutDCAToVertexXY),
211   fCutDCAToVertexZ(that.fCutDCAToVertexZ),
212   fParamType(that.fParamType),
213   fParamMix(that.fParamMix),
214   fTrack(NULL),
215   fTrackPhi(0.),
216   fTrackEta(0.),
217   fTrackWeight(0.),
218   fTrackLabel(INT_MIN),
219   fMCevent(NULL),
220   fMCparticle(NULL),
221   fEvent(NULL),
222   fTPCtrack(),
223   fESDpid(that.fESDpid),
224   fPIDsource(that.fPIDsource),
225   fTPCpidCuts(NULL),
226   fTOFpidCuts(NULL),
227   fTPCTOFpidCrossOverPt(that.fTPCTOFpidCrossOverPt),
228   fAliPID(that.fAliPID)
229 {
230   //copy constructor
231   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
232   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
233   if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
234 }
235
236 //-----------------------------------------------------------------------
237 AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
238 {
239   //assignment
240   AliFlowTrackSimpleCuts::operator=(that);
241   if (that.fAliESDtrackCuts) *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
242   fQA=NULL;
243   fCutMC=that.fCutMC;
244   fCutMCprocessType=that.fCutMCprocessType;
245   fMCprocessType=that.fMCprocessType;
246   fCutMCPID=that.fCutMCPID;
247   fMCPID=that.fMCPID;
248   fIgnoreSignInPID=that.fIgnoreSignInPID,
249   fCutMCisPrimary=that.fCutMCisPrimary;
250   fRequireTransportBitForPrimaries=that.fRequireTransportBitForPrimaries;
251   fMCisPrimary=that.fMCisPrimary;
252   fRequireCharge=that.fRequireCharge;
253   fFakesAreOK=that.fFakesAreOK;
254   fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
255   fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
256   fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
257   fIgnoreTPCzRange=that.fIgnoreTPCzRange;
258   fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
259   fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
260   fCutChi2PerClusterTPC=that.fCutChi2PerClusterTPC;
261   fMaxChi2PerClusterTPC=that.fMaxChi2PerClusterTPC;
262   fMinChi2PerClusterTPC=that.fMinChi2PerClusterTPC;
263   fCutNClustersTPC=that.fCutNClustersTPC;
264   fNClustersTPCMax=that.fNClustersTPCMax;
265   fNClustersTPCMin=that.fNClustersTPCMin;  
266   fCutNClustersITS=that.fCutNClustersITS;
267   fNClustersITSMax=that.fNClustersITSMax;
268   fNClustersITSMin=that.fNClustersITSMin;  
269   fUseAODFilterBit=that.fUseAODFilterBit;
270   fAODFilterBit=that.fAODFilterBit;
271   fCutDCAToVertexXY=that.fCutDCAToVertexXY;
272   fCutDCAToVertexZ=that.fCutDCAToVertexZ;
273   fParamType=that.fParamType;
274   fParamMix=that.fParamMix;
275
276   fTrack=NULL;
277   fTrackPhi=0.;
278   fTrackPhi=0.;
279   fTrackWeight=0.;
280   fTrackLabel=INT_MIN;
281   fMCevent=NULL;
282   fMCparticle=NULL;
283   fEvent=NULL;
284
285   fESDpid = that.fESDpid;
286   fPIDsource = that.fPIDsource;
287
288   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
289   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
290   fTPCTOFpidCrossOverPt=that.fTPCTOFpidCrossOverPt;
291
292   fAliPID=that.fAliPID;
293
294   return *this;
295 }
296
297 //-----------------------------------------------------------------------
298 AliFlowTrackCuts::~AliFlowTrackCuts()
299 {
300   //dtor
301   delete fAliESDtrackCuts;
302   delete fTPCpidCuts;
303   delete fTOFpidCuts;
304 }
305
306 //-----------------------------------------------------------------------
307 void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
308 {
309   //set the event
310   Clear();
311   fEvent=event;
312   fMCevent=mcEvent;
313
314   //do the magic for ESD
315   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(event);
316   if (fCutPID && myESD)
317   {
318     //TODO: maybe call it only for the TOF options?
319     // Added by F. Noferini for TOF PID
320     fESDpid.SetTOFResponse(myESD,AliESDpid::kTOF_T0);
321     fESDpid.MakePID(myESD,kFALSE);
322     // End F. Noferini added part
323   }
324
325   //TODO: AOD
326 }
327
328 //-----------------------------------------------------------------------
329 void AliFlowTrackCuts::SetCutMC( Bool_t b )
330 {
331   //will we be cutting on MC information?
332   fCutMC=b;
333
334   //if we cut on MC info then also the Bethe Bloch should be the one tuned for MC
335   if (fCutMC)
336   {
337     fESDpid.GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
338                                                        1.75295e+01,
339                                                        3.40030e-09,
340                                                        1.96178e+00,
341                                                        3.91720e+00);
342   }
343 }
344
345 //-----------------------------------------------------------------------
346 Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
347 {
348   //check cuts
349   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
350   if (vparticle) return PassesCuts(vparticle);
351   AliFlowTrackSimple* flowtrack = dynamic_cast<AliFlowTrackSimple*>(obj);
352   if (flowtrack) return PassesCuts(flowtrack);
353   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
354   if (tracklets) return PassesCuts(tracklets,id);
355   return kFALSE;  //default when passed wrong type of object
356 }
357
358 //-----------------------------------------------------------------------
359 Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
360 {
361   //check cuts
362   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
363   if (vparticle) 
364   {
365     return PassesMCcuts(fMCevent,vparticle->GetLabel());
366   }
367   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
368   if (tracklets)
369   {
370     Int_t label0 = tracklets->GetLabel(id,0);
371     Int_t label1 = tracklets->GetLabel(id,1);
372     Int_t label = (label0==label1)?tracklets->GetLabel(id,1):-666;
373     return PassesMCcuts(fMCevent,label);
374   }
375   return kFALSE;  //default when passed wrong type of object
376 }
377
378 //-----------------------------------------------------------------------
379 Bool_t AliFlowTrackCuts::PassesCuts(AliFlowTrackSimple* track)
380 {
381   //check cuts on a flowtracksimple
382
383   //clean up from last iteration
384   fTrack = NULL;
385   return AliFlowTrackSimpleCuts::PassesCuts(track);
386 }
387
388 //-----------------------------------------------------------------------
389 Bool_t AliFlowTrackCuts::PassesCuts(AliMultiplicity* tracklet, Int_t id)
390 {
391   //check cuts on a tracklets
392
393   //clean up from last iteration, and init label
394   fTrack = NULL;
395   fMCparticle=NULL;
396   fTrackLabel=-1;
397
398   fTrackPhi = tracklet->GetPhi(id);
399   fTrackEta = tracklet->GetEta(id);
400   fTrackWeight = 1.0;
401   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
402   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
403
404   //check MC info if available
405   //if the 2 clusters have different label track cannot be good
406   //and should therefore not pass the mc cuts
407   Int_t label0 = tracklet->GetLabel(id,0);
408   Int_t label1 = tracklet->GetLabel(id,1);
409   //if possible get label and mcparticle
410   fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
411   if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
412   if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
413   //check MC cuts
414   if (fCutMC && !PassesMCcuts()) return kFALSE;
415   return kTRUE;
416 }
417
418 //-----------------------------------------------------------------------
419 Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
420 {
421   //check the MC info
422   if (!mcEvent) return kFALSE;
423   if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
424   AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
425   if (!mcparticle) {AliError("no MC track"); return kFALSE;}
426
427   if (fCutMCisPrimary)
428   {
429     if (IsPhysicalPrimary(mcEvent,label,fRequireTransportBitForPrimaries) != fMCisPrimary) return kFALSE;
430   }
431   if (fCutMCPID)
432   {
433     Int_t pdgCode = mcparticle->PdgCode();
434     if (fIgnoreSignInPID) 
435     {
436       if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
437     }
438     else 
439     {
440       if (fMCPID != pdgCode) return kFALSE;
441     }
442   }
443   if ( fCutMCprocessType )
444   {
445     TParticle* particle = mcparticle->Particle();
446     Int_t processID = particle->GetUniqueID();
447     if (processID != fMCprocessType ) return kFALSE;
448   }
449   return kTRUE;
450 }
451 //-----------------------------------------------------------------------
452 Bool_t AliFlowTrackCuts::PassesMCcuts()
453 {
454   if (!fMCevent) return kFALSE;
455   if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
456   fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
457   return PassesMCcuts(fMCevent,fTrackLabel);
458 }
459
460 //-----------------------------------------------------------------------
461 Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
462 {
463   //check cuts for an ESD vparticle
464
465   ////////////////////////////////////////////////////////////////
466   //  start by preparing the track parameters to cut on //////////
467   ////////////////////////////////////////////////////////////////
468   //clean up from last iteration
469   fTrack=NULL; 
470
471   //get the label and the mc particle
472   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
473   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
474   else fMCparticle=NULL;
475
476   Bool_t isMCparticle = kFALSE; //some things are different for MC particles, check!
477   AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*>(vparticle);
478   AliAODTrack* aodTrack = NULL;
479   if (esdTrack)
480     //for an ESD track we do some magic sometimes like constructing TPC only parameters
481     //or doing some hybrid, handle that here
482     HandleESDtrack(esdTrack);
483   else
484   {
485     HandleVParticle(vparticle);
486     //now check if produced particle is MC
487     isMCparticle = (dynamic_cast<AliMCParticle*>(fTrack))!=NULL;
488     aodTrack = dynamic_cast<AliAODTrack*>(vparticle); //keep the additional dynamic cast out of the way for ESDs
489   }
490   ////////////////////////////////////////////////////////////////
491   ////////////////////////////////////////////////////////////////
492
493   if (!fTrack) return kFALSE;
494   //because it may be different from global, not needed for aodTrack because we dont do anything funky there
495   if (esdTrack) esdTrack = static_cast<AliESDtrack*>(fTrack);
496   
497   Bool_t pass=kTRUE;
498   //check the common cuts for the current particle fTrack (MC,AOD,ESD)
499   Double_t pt = fTrack->Pt();
500   if (!fFakesAreOK) {if (fTrackLabel<0) pass=kFALSE;}
501   if (fCutPt) {if (pt < fPtMin || pt >= fPtMax ) pass=kFALSE;}
502   if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
503   if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
504   if (fRequireCharge) {if (fTrack->Charge() == 0) pass=kFALSE;}
505   if (fCutCharge && !isMCparticle) {if (fTrack->Charge() != fCharge) pass=kFALSE;}
506   if (fCutCharge && isMCparticle)
507   { 
508     //in case of an MC particle the charge is stored in units of 1/3|e| 
509     Int_t charge = TMath::Nint(fTrack->Charge()/3.0); //mc particles have charge in units of 1/3e
510     if (charge!=fCharge) pass=kFALSE;
511   }
512   //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
513
514   //when additionally MC info is required
515   if (fCutMC && !PassesMCcuts()) pass=kFALSE;
516
517   //the case of ESD or AOD
518   if (esdTrack) PassesESDcuts(esdTrack);
519   if (aodTrack) PassesAODcuts(aodTrack);
520
521   //true by default, if we didn't set any cuts
522   return pass;
523 }
524
525 //_______________________________________________________________________
526 Bool_t AliFlowTrackCuts::PassesAODcuts(AliAODTrack* track)
527 {
528   Bool_t pass = kTRUE;
529
530   if (fCutNClustersTPC)
531   {
532     Int_t ntpccls = track->GetTPCNcls();
533     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
534   }
535
536   if (fCutNClustersITS)
537   {
538     Int_t nitscls = track->GetITSNcls();
539     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
540   }
541   
542    if (fCutChi2PerClusterTPC)
543   {
544     Double_t chi2tpc = track->Chi2perNDF();
545     if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
546   }
547   
548   if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
549   if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
550   
551   if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
552   
553   if (fCutDCAToVertexXY && track->DCA()>GetMaxDCAToVertexXY()) pass=kFALSE;
554     
555
556   return pass;
557 }
558
559 //_______________________________________________________________________
560 Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
561 {
562   Bool_t pass=kTRUE;
563   if (fIgnoreTPCzRange)
564   {
565     const AliExternalTrackParam* pin = track->GetOuterParam();
566     const AliExternalTrackParam* pout = track->GetInnerParam();
567     if (pin&&pout)
568     {
569       Double_t zin = pin->GetZ();
570       Double_t zout = pout->GetZ();
571       if (zin*zout<0) pass=kFALSE;   //reject if cross the membrane
572       if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
573       if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
574     }
575   }
576  
577   Int_t ntpccls = ( fParamType==kESD_TPConly )?
578                     track->GetTPCNclsIter1():track->GetTPCNcls();    
579   if (fCutChi2PerClusterTPC)
580   {
581     Float_t tpcchi2 = (fParamType==kESD_TPConly)?
582                        track->GetTPCchi2Iter1():track->GetTPCchi2();
583     tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
584     if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
585       pass=kFALSE;
586   }
587
588   if (fCutNClustersTPC)
589   {
590     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
591   }
592
593   Int_t nitscls = track->GetNcls(0);
594   if (fCutNClustersITS)
595   {
596     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
597   }
598
599   Double_t pt = fTrack->Pt();
600   if (fCutPID)
601   {
602     switch (fPIDsource)    
603     {
604       case kTPCpid:
605         if (!PassesTPCpidCut(track)) pass=kFALSE;
606         break;
607       case kTOFpid:
608         if (!PassesTOFpidCut(track)) pass=kFALSE;
609         break;
610       case kTPCTOFpid:
611         if (pt< fTPCTOFpidCrossOverPt)
612         {
613           if (!PassesTPCpidCut(track)) pass=kFALSE;
614         }
615         else //if (pt>=fTPCTOFpidCrossOverPt)
616         {
617           if (!PassesTOFpidCut(track)) pass=kFALSE;
618         }
619         break;
620             // part added by F. Noferini
621       case kTOFbayesian:
622               if (!PassesTOFbayesianCut(track)) pass=kFALSE;
623               break;
624             // end part added by F. Noferini
625       default:
626         printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
627         pass=kFALSE;
628         break;
629     }
630   }    
631
632   //some stuff is still handled by AliESDtrackCuts class - delegate
633   if (fAliESDtrackCuts)
634   {
635     if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
636   }
637  
638   return pass;
639 }
640
641 //-----------------------------------------------------------------------
642 void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
643 {
644   //handle the general case
645   switch (fParamType)
646   {
647     default:
648       fTrack = track;
649       break;
650   }
651 }
652
653 //-----------------------------------------------------------------------
654 void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
655 {
656   //handle esd track
657   switch (fParamType)
658   {
659     case kGlobal:
660       fTrack = track;
661       break;
662     case kESD_TPConly:
663       if (!track->FillTPCOnlyTrack(fTPCtrack))
664       {
665         fTrack=NULL;
666         fMCparticle=NULL;
667         fTrackLabel=-1;
668         return;
669       }
670       fTrack = &fTPCtrack;
671       //recalculate the label and mc particle, they may differ as TPClabel != global label
672       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
673       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
674       else fMCparticle=NULL;
675       break;
676     default:
677       fTrack = track;
678       break;
679   }
680 }
681
682 //-----------------------------------------------------------------------
683 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts()
684 {
685   //get standard cuts
686   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly cuts");
687   delete cuts->fAliESDtrackCuts;
688   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
689   cuts->SetParamType(kESD_TPConly);
690   return cuts;
691 }
692
693 //-----------------------------------------------------------------------
694 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
695 {
696   //get standard cuts
697   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
698   delete cuts->fAliESDtrackCuts;
699   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
700   cuts->SetParamType(kGlobal);
701   return cuts;
702 }
703
704 //-----------------------------------------------------------------------
705 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack(int index) const
706 {
707   //get a flow track constructed from whatever we applied cuts on
708   //caller is resposible for deletion
709   //if construction fails return NULL
710   AliFlowTrack* flowtrack=NULL;
711   TParticle *tmpTParticle=NULL;
712   AliMCParticle* tmpAliMCParticle=NULL;
713   if (fParamType==kESD_SPDtracklet)
714   {
715     switch (fParamMix)
716     {
717       case kPure:
718         flowtrack = new AliFlowTrack();
719         flowtrack->SetPhi(fTrackPhi);
720         flowtrack->SetEta(fTrackEta);
721         break;
722       case kTrackWithMCkine:
723         if (!fMCparticle) return NULL;
724         flowtrack = new AliFlowTrack();
725         flowtrack->SetPhi( fMCparticle->Phi() );
726         flowtrack->SetEta( fMCparticle->Eta() );
727         flowtrack->SetPt( fMCparticle->Pt() );
728         break;
729       case kTrackWithMCpt:
730         if (!fMCparticle) return NULL;
731         flowtrack = new AliFlowTrack();
732         flowtrack->SetPhi(fTrackPhi);
733         flowtrack->SetEta(fTrackEta);
734         flowtrack->SetPt(fMCparticle->Pt());
735         break;
736       case kTrackWithPtFromFirstMother:
737         if (!fMCparticle) return NULL;
738         flowtrack = new AliFlowTrack();
739         flowtrack->SetPhi(fTrackPhi);
740         flowtrack->SetEta(fTrackEta);
741         tmpTParticle = fMCparticle->Particle();
742         tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
743         flowtrack->SetPt(tmpAliMCParticle->Pt());
744         break;
745       default:
746         flowtrack = new AliFlowTrack();
747         flowtrack->SetPhi(fTrackPhi);
748         flowtrack->SetEta(fTrackEta);
749         break;
750     }
751     flowtrack->SetSource(AliFlowTrack::kFromTracklet);
752   }
753   else
754   {
755     if (!fTrack) return NULL;
756     switch(fParamMix)
757     {
758       case kPure:
759         flowtrack = new AliFlowTrack(fTrack);
760         break;
761       case kTrackWithMCkine:
762         flowtrack = new AliFlowTrack(fMCparticle);
763         break;
764       case kTrackWithMCPID:
765         flowtrack = new AliFlowTrack(fTrack);
766         //flowtrack->setPID(...) from mc, when implemented
767         break;
768       case kTrackWithMCpt:
769         if (!fMCparticle) return NULL;
770         flowtrack = new AliFlowTrack(fTrack);
771         flowtrack->SetPt(fMCparticle->Pt());
772         break;
773       case kTrackWithPtFromFirstMother:
774         if (!fMCparticle) return NULL;
775         flowtrack = new AliFlowTrack(fTrack);
776         tmpTParticle = fMCparticle->Particle();
777         tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
778         flowtrack->SetPt(tmpAliMCParticle->Pt());
779         break;
780       default:
781         flowtrack = new AliFlowTrack(fTrack);
782         break;
783     }
784     if (fParamType==kMC) flowtrack->SetSource(AliFlowTrack::kFromMC);
785     else if (dynamic_cast<AliESDtrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromESD);
786     else if (dynamic_cast<AliAODTrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromAOD);
787     else if (dynamic_cast<AliMCParticle*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromMC);
788   }
789
790   if(flowtrack)
791     flowtrack->SetIndexOnCollection(index);
792   return flowtrack;
793 }
794
795 //-----------------------------------------------------------------------
796 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
797 {
798   //check if current particle is a physical primary
799   if (!fMCevent) return kFALSE;
800   if (fTrackLabel<0) return kFALSE;
801   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
802 }
803
804 //-----------------------------------------------------------------------
805 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
806 {
807   //check if current particle is a physical primary
808   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
809   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
810   if (!track) return kFALSE;
811   TParticle* particle = track->Particle();
812   Bool_t transported = particle->TestBit(kTransportBit);
813   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
814         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
815   return (physprim && (transported || !requiretransported));
816 }
817
818 //-----------------------------------------------------------------------
819 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type) 
820 {
821   //return the name of the selected parameter type
822   switch (type)
823   {
824     case kMC:
825       return "MC";
826     case kGlobal:
827       return "ESD global";
828     case kESD_TPConly:
829       return "TPC only";
830     case kESD_SPDtracklet:
831       return "SPD tracklet";
832     default:
833       return "unknown";
834   }
835 }
836
837 //-----------------------------------------------------------------------
838 void AliFlowTrackCuts::DefineHistograms()
839 {
840   //define qa histograms
841 }
842
843 //-----------------------------------------------------------------------
844 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
845 {
846   //get the number of tracks in the input event according source
847   //selection (ESD tracks, tracklets, MC particles etc.)
848   AliESDEvent* esd=NULL;
849   switch (fParamType)
850   {
851     case kESD_SPDtracklet:
852       esd = dynamic_cast<AliESDEvent*>(fEvent);
853       if (!esd) return 0;
854       return esd->GetMultiplicity()->GetNumberOfTracklets();
855     case kMC:
856       if (!fMCevent) return 0;
857       return fMCevent->GetNumberOfTracks();
858     default:
859       if (!fEvent) return 0;
860       return fEvent->GetNumberOfTracks();
861   }
862   return 0;
863 }
864
865 //-----------------------------------------------------------------------
866 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
867 {
868   //get the input object according the data source selection:
869   //(esd tracks, traclets, mc particles,etc...)
870   AliESDEvent* esd=NULL;
871   switch (fParamType)
872   {
873     case kESD_SPDtracklet:
874       esd = dynamic_cast<AliESDEvent*>(fEvent);
875       if (!esd) return NULL;
876       return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
877     case kMC:
878       if (!fMCevent) return NULL;
879       return fMCevent->GetTrack(i);
880     default:
881       if (!fEvent) return NULL;
882       return fEvent->GetTrack(i);
883   }
884 }
885
886 //-----------------------------------------------------------------------
887 void AliFlowTrackCuts::Clear(Option_t*)
888 {
889   //clean up
890   fTrack=NULL;
891   fMCevent=NULL;
892   fMCparticle=NULL;
893   fTrackLabel=0;
894   fTrackWeight=0.0;
895   fTrackEta=0.0;
896   fTrackPhi=0.0;
897 }
898
899 //-----------------------------------------------------------------------
900 Bool_t AliFlowTrackCuts::PassesTOFpidCut(AliESDtrack* t )
901 {
902   //check if passes PID cut using timing in TOF
903   Bool_t goodtrack = (t) && 
904                      (t->GetStatus() & AliESDtrack::kTOFpid) && 
905                      (t->GetTOFsignal() > 12000) && 
906                      (t->GetTOFsignal() < 100000) && 
907                      (t->GetIntegratedLength() > 365) && 
908                     !(t->GetStatus() & AliESDtrack::kTOFmismatch);
909
910   if (!goodtrack) return kFALSE;
911   
912   Float_t pt = t->Pt();
913   Float_t p = t->GetP();
914   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
915   Float_t timeTOF = t->GetTOFsignal()- trackT0; 
916   //2=pion 3=kaon 4=protons
917   Double_t inttimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
918   t->GetIntegratedTimes(inttimes);
919   //construct the pid index because it's screwed up in TOF
920   Int_t pid = 0;
921   switch (fAliPID)
922   {
923     case AliPID::kPion:
924       pid=2;
925       break;
926     case AliPID::kKaon:
927       pid=3;
928       break;
929     case AliPID::kProton:
930       pid=4;
931       break;
932     default:
933       return kFALSE;
934   }
935   Float_t s = timeTOF-inttimes[pid];
936
937   Float_t* arr = fTOFpidCuts->GetMatrixArray();
938   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
939   if (col<0) return kFALSE;
940   Float_t min = (*fTOFpidCuts)(1,col);
941   Float_t max = (*fTOFpidCuts)(2,col);
942
943   //printf("--------------TOF pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
944   return (s>min && s<max);
945 }
946
947 //-----------------------------------------------------------------------
948 Bool_t AliFlowTrackCuts::PassesTPCpidCut(AliESDtrack* track)
949 {
950   //check if passes PID cut using dedx signal in the TPC
951   if (!fTPCpidCuts)
952   {
953     printf("no TPCpidCuts\n");
954     return kFALSE;
955   }
956
957   const AliExternalTrackParam* tpcparam = track->GetInnerParam();
958   if (!tpcparam) return kFALSE;
959   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(tpcparam->GetP(), fAliPID);
960   Float_t sigTPC = track->GetTPCsignal();
961   Float_t s = (sigTPC-sigExp)/sigExp;
962   Double_t pt = track->Pt();
963
964   Float_t* arr = fTPCpidCuts->GetMatrixArray();
965   Int_t col = TMath::BinarySearch(fTPCpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
966   if (col<0) return kFALSE;
967   Float_t min = (*fTPCpidCuts)(1,col);
968   Float_t max = (*fTPCpidCuts)(2,col);
969
970   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
971   return (s>min && s<max);
972 }
973
974 //-----------------------------------------------------------------------
975 void AliFlowTrackCuts::InitPIDcuts()
976 {
977   //init matrices with PID cuts
978   TMatrixF* t = NULL;
979   if (!fTPCpidCuts)
980   {
981     if (fAliPID==AliPID::kPion)
982     {
983       t = new TMatrixF(3,10);
984       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.2;
985       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.2;
986       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.25;
987       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.25;
988       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
989       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
990       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
991       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
992       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
993       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  =     0;
994     }
995     else
996     if (fAliPID==AliPID::kKaon)
997     {
998       t = new TMatrixF(3,7);
999       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.4; 
1000       (*t)(0,1)  = 0.25;  (*t)(1,1)  =-0.15;  (*t)(2,1)  = 0.4;
1001       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.1;  (*t)(2,2)  = 0.4;
1002       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.1;  (*t)(2,3)  = 0.4;
1003       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.6;
1004       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.6;
1005       (*t)(0,6)  = 0.50;  (*t)(1,6)  =    0;  (*t)(2,6)  =   0;
1006     }
1007     else
1008     if (fAliPID==AliPID::kProton)
1009     {
1010       t = new TMatrixF(3,16);
1011       (*t)(0,0)  = 0.20;  (*t)(1,0)  =     0;  (*t)(2,0)  =    0; 
1012       (*t)(0,1)  = 0.25;  (*t)(1,1)  =  -0.2;  (*t)(2,1)  =  0.3; 
1013       (*t)(0,2)  = 0.30;  (*t)(1,2)  =  -0.2;  (*t)(2,2)  =  0.6; 
1014       (*t)(0,3)  = 0.35;  (*t)(1,3)  =  -0.2;  (*t)(2,3)  =  0.6; 
1015       (*t)(0,4)  = 0.40;  (*t)(1,4)  =  -0.2;  (*t)(2,4)  =  0.6; 
1016       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.15;  (*t)(2,5)  =  0.6; 
1017       (*t)(0,6)  = 0.50;  (*t)(1,6)  =  -0.1;  (*t)(2,6)  =  0.6; 
1018       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =  0.6; 
1019       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.05;  (*t)(2,8)  = 0.45; 
1020       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.05;  (*t)(2,9)  = 0.45; 
1021       (*t)(0,10) = 0.70;  (*t)(1,10) = -0.05;  (*t)(2,10) = 0.45; 
1022       (*t)(0,11) = 0.75;  (*t)(1,11) = -0.05;  (*t)(2,11) = 0.45; 
1023       (*t)(0,12) = 0.80;  (*t)(1,12) =     0;  (*t)(2,12) = 0.45; 
1024       (*t)(0,13) = 0.85;  (*t)(1,13) =     0;  (*t)(2,13) = 0.45; 
1025       (*t)(0,14) = 0.90;  (*t)(1,14) =     0;  (*t)(2,14) = 0.45;
1026       (*t)(0,15) = 0.95;  (*t)(1,15) =     0;  (*t)(2,15) =    0;
1027     }
1028     fTPCpidCuts=t;
1029   }
1030   t = NULL;
1031   if (!fTOFpidCuts)
1032   {
1033     if (fAliPID==AliPID::kPion)
1034     {
1035       t = new TMatrixF(3,31);
1036       (*t)(0,0)  = 0.3;   (*t)(1,0)  = -700;  (*t)(2,0)  = 700;
1037       (*t)(0,1)  = 0.35;  (*t)(1,1)  = -800;  (*t)(2,1)  = 800;
1038       (*t)(0,2)  = 0.40;  (*t)(1,2)  = -600;  (*t)(2,2)  = 800;
1039       (*t)(0,3)  = 0.45;  (*t)(1,3)  = -500;  (*t)(2,3)  = 700;
1040       (*t)(0,4)  = 0.50;  (*t)(1,4)  = -400;  (*t)(2,4)  = 700;
1041       (*t)(0,5)  = 0.55;  (*t)(1,5)  = -400;  (*t)(2,5)  = 700;
1042       (*t)(0,6)  = 0.60;  (*t)(1,6)  = -400;  (*t)(2,6)  = 700;
1043       (*t)(0,7)  = 0.65;  (*t)(1,7)  = -400;  (*t)(2,7)  = 700;
1044       (*t)(0,8)  = 0.70;  (*t)(1,8)  = -400;  (*t)(2,8)  = 700;
1045       (*t)(0,9)  = 0.75;  (*t)(1,9)  = -400;  (*t)(2,9)  = 700;
1046       (*t)(0,10) = 0.80;  (*t)(1,10) = -400;  (*t)(2,10) = 600;
1047       (*t)(0,11) = 0.85;  (*t)(1,11) = -400;  (*t)(2,11) = 600;
1048       (*t)(0,12) = 0.90;  (*t)(1,12) = -400;  (*t)(2,12) = 600;
1049       (*t)(0,13) = 0.95;  (*t)(1,13) = -400;  (*t)(2,13) = 600;
1050       (*t)(0,14) = 1.00;  (*t)(1,14) = -400;  (*t)(2,14) = 550;
1051       (*t)(0,15) = 1.10;  (*t)(1,15) = -400;  (*t)(2,15) = 450;
1052       (*t)(0,16) = 1.20;  (*t)(1,16) = -400;  (*t)(2,16) = 400;
1053       (*t)(0,17) = 1.30;  (*t)(1,17) = -400;  (*t)(2,17) = 300;
1054       (*t)(0,18) = 1.40;  (*t)(1,18) = -400;  (*t)(2,18) = 300;
1055       (*t)(0,19) = 1.50;  (*t)(1,19) = -400;  (*t)(2,19) = 250;
1056       (*t)(0,20) = 1.60;  (*t)(1,20) = -400;  (*t)(2,20) = 200;
1057       (*t)(0,21) = 1.70;  (*t)(1,21) = -400;  (*t)(2,21) = 150;
1058       (*t)(0,22) = 1.80;  (*t)(1,22) = -400;  (*t)(2,22) = 100;
1059       (*t)(0,23) = 1.90;  (*t)(1,23) = -400;  (*t)(2,23) =  70;
1060       (*t)(0,24) = 2.00;  (*t)(1,24) = -400;  (*t)(2,24) =  50;
1061       (*t)(0,25) = 2.10;  (*t)(1,25) = -400;  (*t)(2,25) =   0;
1062       (*t)(0,26) = 2.20;  (*t)(1,26) = -400;  (*t)(2,26) =   0;
1063       (*t)(0,27) = 2.30;  (*t)(1,26) = -400;  (*t)(2,26) =   0;
1064       (*t)(0,28) = 2.40;  (*t)(1,26) = -400;  (*t)(2,26) = -50;
1065       (*t)(0,29) = 2.50;  (*t)(1,26) = -400;  (*t)(2,26) = -50;
1066       (*t)(0,30) = 2.60;  (*t)(1,26) =    0;  (*t)(2,26) =   0;
1067     }
1068     else
1069     if (fAliPID==AliPID::kProton)
1070     {
1071       t = new TMatrixF(3,39);
1072       (*t)(0,0)  = 0.3;  (*t)(1,0)   = 0;     (*t)(2,0) = 0;
1073       (*t)(0,1)  = 0.35;  (*t)(1,1)  = 0;     (*t)(2,1) = 0;
1074       (*t)(0,2)  = 0.40;  (*t)(1,2)  = 0;     (*t)(2,2) = 0;
1075       (*t)(0,3)  = 0.45;  (*t)(1,3)  = 0;     (*t)(2,3) = 0;
1076       (*t)(0,4)  = 0.50;  (*t)(1,4)  = 0;     (*t)(2,4) = 0;
1077       (*t)(0,5)  = 0.55;  (*t)(1,5)  = -900;  (*t)(2,5)  = 600;
1078       (*t)(0,6)  = 0.60;  (*t)(1,6)  = -800;  (*t)(2,6)  = 600;
1079       (*t)(0,7)  = 0.65;  (*t)(1,7)  = -800;  (*t)(2,7)  = 600;
1080       (*t)(0,8)  = 0.70;  (*t)(1,8)  = -800;  (*t)(2,8)  = 600;
1081       (*t)(0,9)  = 0.75;  (*t)(1,9)  = -700;  (*t)(2,9)  = 500;
1082       (*t)(0,10) = 0.80;  (*t)(1,10) = -700;  (*t)(2,10) = 500;
1083       (*t)(0,11) = 0.85;  (*t)(1,11) = -700;  (*t)(2,11) = 500;
1084       (*t)(0,12) = 0.90;  (*t)(1,12) = -600;  (*t)(2,12) = 500;
1085       (*t)(0,13) = 0.95;  (*t)(1,13) = -600;  (*t)(2,13) = 500;
1086       (*t)(0,14) = 1.00;  (*t)(1,14) = -600;  (*t)(2,14) = 500;
1087       (*t)(0,15) = 1.10;  (*t)(1,15) = -600;  (*t)(2,15) = 500;
1088       (*t)(0,16) = 1.20;  (*t)(1,16) = -500;  (*t)(2,16) = 500;
1089       (*t)(0,17) = 1.30;  (*t)(1,17) = -500;  (*t)(2,17) = 500;
1090       (*t)(0,18) = 1.40;  (*t)(1,18) = -500;  (*t)(2,18) = 500;
1091       (*t)(0,19) = 1.50;  (*t)(1,19) = -500;  (*t)(2,19) = 500;
1092       (*t)(0,20) = 1.60;  (*t)(1,20) = -400;  (*t)(2,20) = 500;
1093       (*t)(0,21) = 1.70;  (*t)(1,21) = -400;  (*t)(2,21) = 500;
1094       (*t)(0,22) = 1.80;  (*t)(1,22) = -400;  (*t)(2,22) = 500;
1095       (*t)(0,23) = 1.90;  (*t)(1,23) = -400;  (*t)(2,23) = 500;
1096       (*t)(0,24) = 2.00;  (*t)(1,24) = -400;  (*t)(2,24) = 500;
1097       (*t)(0,25) = 2.10;  (*t)(1,25) = -350;  (*t)(2,25) = 500;
1098       (*t)(0,26) = 2.20;  (*t)(1,26) = -350;  (*t)(2,26) = 500;
1099       (*t)(0,27) = 2.30;  (*t)(1,27) = -300;  (*t)(2,27) = 500;
1100       (*t)(0,28) = 2.40;  (*t)(1,28) = -300;  (*t)(2,28) = 500;
1101       (*t)(0,29) = 2.50;  (*t)(1,29) = -300;  (*t)(2,29) = 500;
1102       (*t)(0,30) = 2.60;  (*t)(1,30) = -250;  (*t)(2,30) = 500;
1103       (*t)(0,31) = 2.70;  (*t)(1,31) = -200;  (*t)(2,31) = 500;
1104       (*t)(0,32) = 2.80;  (*t)(1,32) = -150;  (*t)(2,32) = 500;
1105       (*t)(0,33) = 2.90;  (*t)(1,33) = -150;  (*t)(2,33) = 500;
1106       (*t)(0,34) = 3.00;  (*t)(1,34) = -100;  (*t)(2,34) = 400;
1107       (*t)(0,35) = 3.10;  (*t)(1,35) = -100;  (*t)(2,35) = 400;
1108       (*t)(0,36) = 3.50;  (*t)(1,36) = -100;  (*t)(2,36) = 400;
1109       (*t)(0,37) = 3.60;  (*t)(1,37) =    0;  (*t)(2,37) = 0;
1110       (*t)(0,38) = 3.70;  (*t)(1,38) =    0;  (*t)(2,38) = 0;
1111     }
1112     else
1113     if (fAliPID==AliPID::kKaon)
1114     {
1115       t = new TMatrixF(3,26);
1116       (*t)(0,0)  = 0.3;   (*t)(1,0)  =    0;  (*t)(2,0)  =    0;
1117       (*t)(0,1)  = 0.35;  (*t)(1,1)  =    0;  (*t)(2,1)  =    0;
1118       (*t)(0,2)  = 0.40;  (*t)(1,2)  = -800;  (*t)(2,2)  =  600;
1119       (*t)(0,3)  = 0.45;  (*t)(1,3)  = -800;  (*t)(2,3)  =  600;
1120       (*t)(0,4)  = 0.50;  (*t)(1,4)  = -800;  (*t)(2,4)  =  600;
1121       (*t)(0,5)  = 0.55;  (*t)(1,5)  = -800;  (*t)(2,5)  =  600;
1122       (*t)(0,6)  = 0.60;  (*t)(1,6)  = -800;  (*t)(2,6)  =  600;
1123       (*t)(0,7)  = 0.65;  (*t)(1,7)  = -700;  (*t)(2,7)  =  600;
1124       (*t)(0,8)  = 0.70;  (*t)(1,8)  = -600;  (*t)(2,8)  =  600;
1125       (*t)(0,9)  = 0.75;  (*t)(1,9)  = -600;  (*t)(2,9)  =  500;
1126       (*t)(0,10) = 0.80;  (*t)(1,10) = -500;  (*t)(2,10) =  500;
1127       (*t)(0,11) = 0.85;  (*t)(1,11) = -500;  (*t)(2,11) =  500;
1128       (*t)(0,12) = 0.90;  (*t)(1,12) = -400;  (*t)(2,12) =  500;
1129       (*t)(0,13) = 0.95;  (*t)(1,13) = -400;  (*t)(2,13) =  500;
1130       (*t)(0,14) = 1.00;  (*t)(1,14) = -400;  (*t)(2,14) =  500;
1131       (*t)(0,15) = 1.10;  (*t)(1,15) = -350;  (*t)(2,15) =  450;
1132       (*t)(0,16) = 1.20;  (*t)(1,16) = -300;  (*t)(2,16) =  400;
1133       (*t)(0,17) = 1.30;  (*t)(1,17) = -300;  (*t)(2,17) =  400;
1134       (*t)(0,18) = 1.40;  (*t)(1,18) = -250;  (*t)(2,18) =  400;
1135       (*t)(0,19) = 1.50;  (*t)(1,19) = -200;  (*t)(2,19) =  400;
1136       (*t)(0,20) = 1.60;  (*t)(1,20) = -150;  (*t)(2,20) =  400;
1137       (*t)(0,21) = 1.70;  (*t)(1,21) = -100;  (*t)(2,21) =  400;
1138       (*t)(0,22) = 1.80;  (*t)(1,22) =  -50;  (*t)(2,22) =  400;
1139       (*t)(0,23) = 1.90;  (*t)(1,22) =    0;  (*t)(2,22) =  400;
1140       (*t)(0,24) = 2.00;  (*t)(1,22) =   50;  (*t)(2,22) =  400;
1141       (*t)(0,25) = 2.10;  (*t)(1,22) =    0;  (*t)(2,22) =    0;
1142     }
1143     fTOFpidCuts=t;
1144   }
1145 }
1146
1147 //-----------------------------------------------------------------------
1148 // part added by F. Noferini (some methods)
1149 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(AliESDtrack* track){
1150   SetPriors();
1151
1152   Bool_t goodtrack = track && (track->GetStatus() & AliESDtrack::kTOFpid) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000) && (track->GetIntegratedLength() > 365) && !(track->GetStatus() & AliESDtrack::kTOFmismatch);
1153
1154   if (! goodtrack)
1155        return kFALSE;
1156
1157   Int_t pdg = GetESDPdg(track,"bayesianALL");
1158   //  printf("pdg set to %i\n",pdg);
1159
1160   Int_t pid = 0;
1161   Float_t prob = 0;
1162   switch (fAliPID)
1163   {
1164     case AliPID::kPion:
1165       pid=211;
1166       prob = fProbBayes[2];
1167       break;
1168     case AliPID::kKaon:
1169       pid=321;
1170       prob = fProbBayes[3];
1171      break;
1172     case AliPID::kProton:
1173       pid=2212;
1174       prob = fProbBayes[4];
1175       break;
1176     case AliPID::kElectron:
1177       pid=-11;
1178        prob = fProbBayes[0];
1179      break;
1180     default:
1181       return kFALSE;
1182   }
1183
1184   //  printf("pt = %f -- all prob = [%4.2f,%4.2f,%4.2f,%4.2f,%4.2f] -- prob = %f\n",track->Pt(),fProbBayes[0],fProbBayes[1],fProbBayes[2],fProbBayes[3],fProbBayes[4],prob);
1185   if(TMath::Abs(pdg) == TMath::Abs(pid) && prob > 0.8){
1186     if(!fCutCharge)
1187       return kTRUE;
1188     else if (fCutCharge && fCharge * track->GetSign() > 0)
1189       return kTRUE;
1190   }
1191   return kFALSE;
1192 }
1193 //-----------------------------------------------------------------------
1194 Int_t AliFlowTrackCuts::GetESDPdg(AliESDtrack *track,Option_t *option,Int_t ipart,Float_t cPi,Float_t cKa,Float_t cPr){
1195   Int_t pdg = 0;
1196   Int_t pdgvalues[5] = {-11,-13,211,321,2212};
1197   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
1198
1199   if(strstr(option,"bayesianTOF")){ // Bayesian TOF PID
1200     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1201     Double_t rcc=0.;
1202     
1203     Float_t pt = track->Pt();
1204     
1205     Int_t iptesd = 0;
1206     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1207   
1208     if(cPi < 0){
1209       c[0] = fC[iptesd][0];
1210       c[1] = fC[iptesd][1];
1211       c[2] = fC[iptesd][2];
1212       c[3] = fC[iptesd][3];
1213       c[4] = fC[iptesd][4];
1214     }
1215     else{
1216       c[0] = 0.0;
1217       c[1] = 0.0;
1218       c[2] = cPi;
1219       c[3] = cKa;
1220       c[4] = cPr;      
1221     }
1222
1223     Double_t r1[10]; track->GetTOFpid(r1);
1224     
1225     Int_t i;
1226     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1227     
1228     Double_t w[10];
1229     for (i=0; i<5; i++){
1230         w[i]=c[i]*r1[i]/rcc;
1231         fProbBayes[i] = w[i];
1232     }
1233     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1234       pdg = 211*Int_t(track->GetSign());
1235     }
1236     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1237       pdg = 2212*Int_t(track->GetSign());
1238     }
1239     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1240       pdg = 321*Int_t(track->GetSign());
1241     }
1242     else if (w[0]>=w[1]) { //electrons
1243       pdg = -11*Int_t(track->GetSign());
1244     }
1245     else{ // muon
1246       pdg = -13*Int_t(track->GetSign());
1247     }
1248   }
1249
1250   else if(strstr(option,"bayesianTPC")){ // Bayesian TPC PID
1251     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1252     Double_t rcc=0.;
1253     
1254     Float_t pt = track->Pt();
1255     
1256     Int_t iptesd = 0;
1257     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1258   
1259     if(cPi < 0){
1260       c[0] = fC[iptesd][0];
1261       c[1] = fC[iptesd][1];
1262       c[2] = fC[iptesd][2];
1263       c[3] = fC[iptesd][3];
1264       c[4] = fC[iptesd][4];
1265     }
1266     else{
1267       c[0] = 0.0;
1268       c[1] = 0.0;
1269       c[2] = cPi;
1270       c[3] = cKa;
1271       c[4] = cPr;      
1272     }
1273
1274     Double_t r1[10]; track->GetTPCpid(r1);
1275     
1276     Int_t i;
1277     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1278     
1279     Double_t w[10];
1280     for (i=0; i<5; i++){
1281         w[i]=c[i]*r1[i]/rcc;
1282         fProbBayes[i] = w[i];
1283     }
1284     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1285       pdg = 211*Int_t(track->GetSign());
1286     }
1287     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1288       pdg = 2212*Int_t(track->GetSign());
1289     }
1290     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1291       pdg = 321*Int_t(track->GetSign());
1292     }
1293     else if (w[0]>=w[1]) { //electrons
1294       pdg = -11*Int_t(track->GetSign());
1295     }
1296     else{ // muon
1297       pdg = -13*Int_t(track->GetSign());
1298     }
1299   }
1300   
1301   else if(strstr(option,"bayesianALL")){
1302     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1303     Double_t rcc=0.;
1304     
1305     Float_t pt = track->Pt();
1306     
1307     Int_t iptesd = 0;
1308     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1309
1310     if(cPi < 0){
1311       c[0] = fC[iptesd][0];
1312       c[1] = fC[iptesd][1];
1313       c[2] = fC[iptesd][2];
1314       c[3] = fC[iptesd][3];
1315       c[4] = fC[iptesd][4];
1316     }
1317     else{
1318       c[0] = 0.0;
1319       c[1] = 0.0;
1320       c[2] = cPi;
1321       c[3] = cKa;
1322       c[4] = cPr;      
1323     }
1324
1325     Double_t r1[10]; track->GetTOFpid(r1);
1326     Double_t r2[10]; track->GetTPCpid(r2);
1327
1328     Int_t i;
1329     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]*r2[i]);
1330     
1331
1332     Double_t w[10];
1333     for (i=0; i<5; i++){
1334         w[i]=c[i]*r1[i]*r2[i]/rcc;
1335         fProbBayes[i] = w[i];
1336     }
1337
1338     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1339       pdg = 211*Int_t(track->GetSign());
1340     }
1341     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1342       pdg = 2212*Int_t(track->GetSign());
1343     }
1344     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1345       pdg = 321*Int_t(track->GetSign());
1346     }
1347     else if (w[0]>=w[1]) { //electrons
1348       pdg = -11*Int_t(track->GetSign());
1349     }
1350     else{ // muon
1351       pdg = -13*Int_t(track->GetSign());
1352     }
1353   }
1354
1355   else if(strstr(option,"sigmacutTOF")){
1356     printf("PID not implemented yet: %s\nNO PID!!!!\n",option);
1357     Float_t p = track->P();
1358
1359     // Take expected times
1360     Double_t exptimes[5];
1361     track->GetIntegratedTimes(exptimes);
1362
1363     // Take resolution for TOF response
1364     // like fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1365     Float_t resolution = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1366
1367     if(TMath::Abs(exptimes[ipart] - track->GetTOFsignal()) < 3 * resolution){
1368       pdg = pdgvalues[ipart] * Int_t(track->GetSign());
1369     }
1370   }
1371
1372   else{
1373     printf("Invalid PID option: %s\nNO PID!!!!\n",option);
1374   }
1375
1376   return pdg;
1377 }
1378 //-----------------------------------------------------------------------
1379 void AliFlowTrackCuts::SetPriors(){
1380   // set abbundancies
1381     fBinLimitPID[0] = 0.30;
1382     fC[0][0] = 0.015;
1383     fC[0][1] = 0.015;
1384     fC[0][2] = 1;
1385     fC[0][3] = 0.0025;
1386     fC[0][4] = 0.000015;
1387     fBinLimitPID[1] = 0.35;
1388     fC[1][0] = 0.015;
1389     fC[1][1] = 0.015;
1390     fC[1][2] = 1;
1391     fC[1][3] = 0.01;
1392     fC[1][4] = 0.001;
1393     fBinLimitPID[2] = 0.40;
1394     fC[2][0] = 0.015;
1395     fC[2][1] = 0.015;
1396     fC[2][2] = 1;
1397     fC[2][3] = 0.026;
1398     fC[2][4] = 0.004;
1399     fBinLimitPID[3] = 0.45;
1400     fC[3][0] = 0.015;
1401     fC[3][1] = 0.015;
1402     fC[3][2] = 1;
1403     fC[3][3] = 0.026;
1404     fC[3][4] = 0.004;
1405     fBinLimitPID[4] = 0.50;
1406     fC[4][0] = 0.015;
1407     fC[4][1] = 0.015;
1408     fC[4][2] = 1.000000;
1409     fC[4][3] = 0.05;
1410     fC[4][4] = 0.01;
1411     fBinLimitPID[5] = 0.60;
1412     fC[5][0] = 0.012;
1413     fC[5][1] = 0.012;
1414     fC[5][2] = 1;
1415     fC[5][3] = 0.085;
1416     fC[5][4] = 0.022;
1417     fBinLimitPID[6] = 0.70;
1418     fC[6][0] = 0.01;
1419     fC[6][1] = 0.01;
1420     fC[6][2] = 1;
1421     fC[6][3] = 0.12;
1422     fC[6][4] = 0.036;
1423     fBinLimitPID[7] = 0.80;
1424     fC[7][0] = 0.0095;
1425     fC[7][1] = 0.0095;
1426     fC[7][2] = 1;
1427     fC[7][3] = 0.15;
1428     fC[7][4] = 0.05;
1429     fBinLimitPID[8] = 0.90;
1430     fC[8][0] = 0.0085;
1431     fC[8][1] = 0.0085;
1432     fC[8][2] = 1;
1433     fC[8][3] = 0.18;
1434     fC[8][4] = 0.074;
1435     fBinLimitPID[9] = 1;
1436     fC[9][0] = 0.008;
1437     fC[9][1] = 0.008;
1438     fC[9][2] = 1;
1439     fC[9][3] = 0.22;
1440     fC[9][4] = 0.1;
1441     fBinLimitPID[10] = 1.20;
1442     fC[10][0] = 0.007;
1443     fC[10][1] = 0.007;
1444     fC[10][2] = 1;
1445     fC[10][3] = 0.28;
1446     fC[10][4] = 0.16;
1447     fBinLimitPID[11] = 1.40;
1448     fC[11][0] = 0.0066;
1449     fC[11][1] = 0.0066;
1450     fC[11][2] = 1;
1451     fC[11][3] = 0.35;
1452     fC[11][4] = 0.23;
1453     fBinLimitPID[12] = 1.60;
1454     fC[12][0] = 0.0075;
1455     fC[12][1] = 0.0075;
1456     fC[12][2] = 1;
1457     fC[12][3] = 0.40;
1458     fC[12][4] = 0.31;
1459     fBinLimitPID[13] = 1.80;
1460     fC[13][0] = 0.0062;
1461     fC[13][1] = 0.0062;
1462     fC[13][2] = 1;
1463     fC[13][3] = 0.45;
1464     fC[13][4] = 0.39;
1465     fBinLimitPID[14] = 2.00;
1466     fC[14][0] = 0.005;
1467     fC[14][1] = 0.005;
1468     fC[14][2] = 1;
1469     fC[14][3] = 0.46;
1470     fC[14][4] = 0.47;
1471     fBinLimitPID[15] = 2.20;
1472     fC[15][0] = 0.0042;
1473     fC[15][1] = 0.0042;
1474     fC[15][2] = 1;
1475     fC[15][3] = 0.5;
1476     fC[15][4] = 0.55;
1477     fBinLimitPID[16] = 2.40;
1478     fC[16][0] = 0.007;
1479     fC[16][1] = 0.007;
1480     fC[16][2] = 1;
1481     fC[16][3] = 0.5;
1482     fC[16][4] = 0.6;
1483     
1484     for(Int_t i=17;i<fnPIDptBin;i++){
1485         fBinLimitPID[i] = 2.0 + 0.2 * (i-14);
1486         fC[i][0] = fC[13][0];
1487         fC[i][1] = fC[13][1];
1488         fC[i][2] = fC[13][2];
1489         fC[i][3] = fC[13][3];
1490         fC[i][4] = fC[13][4];
1491     }  
1492 }
1493 // end part added by F. Noferini
1494 //-----------------------------------------------------------------------