]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG2/FLOW/AliFlowTasks/AliFlowTrackCuts.cxx
hooks for PMD flow analysis
[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 #include "AliESDPmdTrack.h"
55
56 ClassImp(AliFlowTrackCuts)
57
58 //-----------------------------------------------------------------------
59 AliFlowTrackCuts::AliFlowTrackCuts():
60   AliFlowTrackSimpleCuts(),
61   fAliESDtrackCuts(NULL),
62   fQA(NULL),
63   fCutMC(kFALSE),
64   fCutMCprocessType(kFALSE),
65   fMCprocessType(kPNoProcess),
66   fCutMCPID(kFALSE),
67   fMCPID(0),
68   fIgnoreSignInPID(kFALSE),
69   fCutMCisPrimary(kFALSE),
70   fRequireTransportBitForPrimaries(kTRUE),
71   fMCisPrimary(kFALSE),
72   fRequireCharge(kFALSE),
73   fFakesAreOK(kTRUE),
74   fCutSPDtrackletDeltaPhi(kFALSE),
75   fSPDtrackletDeltaPhiMax(FLT_MAX),
76   fSPDtrackletDeltaPhiMin(-FLT_MAX),
77   fIgnoreTPCzRange(kFALSE),
78   fIgnoreTPCzRangeMax(FLT_MAX),
79   fIgnoreTPCzRangeMin(-FLT_MAX),
80   fCutChi2PerClusterTPC(kFALSE),
81   fMaxChi2PerClusterTPC(FLT_MAX),
82   fMinChi2PerClusterTPC(-FLT_MAX),
83   fCutNClustersTPC(kFALSE),
84   fNClustersTPCMax(INT_MAX),
85   fNClustersTPCMin(INT_MIN),  
86   fCutNClustersITS(kFALSE),
87   fNClustersITSMax(INT_MAX),
88   fNClustersITSMin(INT_MIN),  
89   fUseAODFilterBit(kFALSE),
90   fAODFilterBit(0),
91   fCutDCAToVertexXY(kFALSE),
92   fCutDCAToVertexZ(kFALSE),
93   fCutMinimalTPCdedx(kFALSE),
94   fMinimalTPCdedx(0.),
95   fParamType(kGlobal),
96   fParamMix(kPure),
97   fTrack(NULL),
98   fTrackPhi(0.),
99   fTrackEta(0.),
100   fTrackWeight(0.),
101   fTrackLabel(INT_MIN),
102   fMCevent(NULL),
103   fMCparticle(NULL),
104   fEvent(NULL),
105   fTPCtrack(),
106   fESDpid(),
107   fPIDsource(kTOFpid),
108   fTPCpidCuts(NULL),
109   fTOFpidCuts(NULL),
110   fParticleID(AliPID::kPion),
111   fParticleProbability(.9)
112 {
113   //io constructor 
114   for ( Int_t i=0; i<5; i++ ) { fProbBayes[i]=0.0; }
115   SetPriors(); //init arrays
116 }
117
118 //-----------------------------------------------------------------------
119 AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
120   AliFlowTrackSimpleCuts(),
121   fAliESDtrackCuts(new AliESDtrackCuts()),
122   fQA(NULL),
123   fCutMC(kFALSE),
124   fCutMCprocessType(kFALSE),
125   fMCprocessType(kPNoProcess),
126   fCutMCPID(kFALSE),
127   fMCPID(0),
128   fIgnoreSignInPID(kFALSE),
129   fCutMCisPrimary(kFALSE),
130   fRequireTransportBitForPrimaries(kTRUE),
131   fMCisPrimary(kFALSE),
132   fRequireCharge(kFALSE),
133   fFakesAreOK(kTRUE),
134   fCutSPDtrackletDeltaPhi(kFALSE),
135   fSPDtrackletDeltaPhiMax(FLT_MAX),
136   fSPDtrackletDeltaPhiMin(-FLT_MAX),
137   fIgnoreTPCzRange(kFALSE),
138   fIgnoreTPCzRangeMax(FLT_MAX),
139   fIgnoreTPCzRangeMin(-FLT_MAX),
140   fCutChi2PerClusterTPC(kFALSE),
141   fMaxChi2PerClusterTPC(FLT_MAX),
142   fMinChi2PerClusterTPC(-FLT_MAX),
143   fCutNClustersTPC(kFALSE),
144   fNClustersTPCMax(INT_MAX),
145   fNClustersTPCMin(INT_MIN),  
146   fCutNClustersITS(kFALSE),
147   fNClustersITSMax(INT_MAX),
148   fNClustersITSMin(INT_MIN),
149   fUseAODFilterBit(kFALSE),
150   fAODFilterBit(0),
151   fCutDCAToVertexXY(kFALSE),
152   fCutDCAToVertexZ(kFALSE),
153   fCutMinimalTPCdedx(kFALSE),
154   fMinimalTPCdedx(0.),
155   fParamType(kGlobal),
156   fParamMix(kPure),
157   fTrack(NULL),
158   fTrackPhi(0.),
159   fTrackEta(0.),
160   fTrackWeight(0.),
161   fTrackLabel(INT_MIN),
162   fMCevent(NULL),
163   fMCparticle(NULL),
164   fEvent(NULL),
165   fTPCtrack(),
166   fESDpid(),
167   fPIDsource(kTOFpid),
168   fTPCpidCuts(NULL),
169   fTOFpidCuts(NULL),
170   fParticleID(AliPID::kPion),
171   fParticleProbability(.9)
172 {
173   //constructor 
174   SetName(name);
175   SetTitle("AliFlowTrackCuts");
176   fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
177                                                     2.63394e+01,
178                                                     5.04114e-11,
179                                                     2.12543e+00,
180                                                     4.88663e+00 );
181   for ( Int_t i=0; i<5; i++ ) { fProbBayes[i]=0.0; }
182   SetPriors(); //init arrays
183 }
184
185 //-----------------------------------------------------------------------
186 AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
187   AliFlowTrackSimpleCuts(that),
188   fAliESDtrackCuts(NULL),
189   fQA(NULL),
190   fCutMC(that.fCutMC),
191   fCutMCprocessType(that.fCutMCprocessType),
192   fMCprocessType(that.fMCprocessType),
193   fCutMCPID(that.fCutMCPID),
194   fMCPID(that.fMCPID),
195   fIgnoreSignInPID(that.fIgnoreSignInPID),
196   fCutMCisPrimary(that.fCutMCisPrimary),
197   fRequireTransportBitForPrimaries(that.fRequireTransportBitForPrimaries),
198   fMCisPrimary(that.fMCisPrimary),
199   fRequireCharge(that.fRequireCharge),
200   fFakesAreOK(that.fFakesAreOK),
201   fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
202   fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
203   fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
204   fIgnoreTPCzRange(that.fIgnoreTPCzRange),
205   fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
206   fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
207   fCutChi2PerClusterTPC(that.fCutChi2PerClusterTPC),
208   fMaxChi2PerClusterTPC(that.fMaxChi2PerClusterTPC),
209   fMinChi2PerClusterTPC(that.fMinChi2PerClusterTPC),
210   fCutNClustersTPC(that.fCutNClustersTPC),
211   fNClustersTPCMax(that.fNClustersTPCMax),
212   fNClustersTPCMin(that.fNClustersTPCMin),
213   fCutNClustersITS(that.fCutNClustersITS),
214   fNClustersITSMax(that.fNClustersITSMax),
215   fNClustersITSMin(that.fNClustersITSMin),
216   fUseAODFilterBit(that.fUseAODFilterBit),
217   fAODFilterBit(that.fAODFilterBit),
218   fCutDCAToVertexXY(that.fCutDCAToVertexXY),
219   fCutDCAToVertexZ(that.fCutDCAToVertexZ),
220   fCutMinimalTPCdedx(that.fCutMinimalTPCdedx),
221   fMinimalTPCdedx(that.fMinimalTPCdedx),
222   fParamType(that.fParamType),
223   fParamMix(that.fParamMix),
224   fTrack(NULL),
225   fTrackPhi(0.),
226   fTrackEta(0.),
227   fTrackWeight(0.),
228   fTrackLabel(INT_MIN),
229   fMCevent(NULL),
230   fMCparticle(NULL),
231   fEvent(NULL),
232   fTPCtrack(),
233   fESDpid(that.fESDpid),
234   fPIDsource(that.fPIDsource),
235   fTPCpidCuts(NULL),
236   fTOFpidCuts(NULL),
237   fParticleID(that.fParticleID),
238   fParticleProbability(that.fParticleProbability)
239 {
240   //copy constructor
241   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
242   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
243   if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
244   memcpy(fProbBayes,that.fProbBayes,sizeof(fProbBayes));
245   SetPriors(); //init arrays
246 }
247
248 //-----------------------------------------------------------------------
249 AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
250 {
251   //assignment
252   if (this==&that) return *this;
253
254   AliFlowTrackSimpleCuts::operator=(that);
255   if (that.fAliESDtrackCuts) *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
256   fQA=NULL;
257   fCutMC=that.fCutMC;
258   fCutMCprocessType=that.fCutMCprocessType;
259   fMCprocessType=that.fMCprocessType;
260   fCutMCPID=that.fCutMCPID;
261   fMCPID=that.fMCPID;
262   fIgnoreSignInPID=that.fIgnoreSignInPID,
263   fCutMCisPrimary=that.fCutMCisPrimary;
264   fRequireTransportBitForPrimaries=that.fRequireTransportBitForPrimaries;
265   fMCisPrimary=that.fMCisPrimary;
266   fRequireCharge=that.fRequireCharge;
267   fFakesAreOK=that.fFakesAreOK;
268   fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
269   fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
270   fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
271   fIgnoreTPCzRange=that.fIgnoreTPCzRange;
272   fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
273   fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
274   fCutChi2PerClusterTPC=that.fCutChi2PerClusterTPC;
275   fMaxChi2PerClusterTPC=that.fMaxChi2PerClusterTPC;
276   fMinChi2PerClusterTPC=that.fMinChi2PerClusterTPC;
277   fCutNClustersTPC=that.fCutNClustersTPC;
278   fNClustersTPCMax=that.fNClustersTPCMax;
279   fNClustersTPCMin=that.fNClustersTPCMin;  
280   fCutNClustersITS=that.fCutNClustersITS;
281   fNClustersITSMax=that.fNClustersITSMax;
282   fNClustersITSMin=that.fNClustersITSMin;  
283   fUseAODFilterBit=that.fUseAODFilterBit;
284   fAODFilterBit=that.fAODFilterBit;
285   fCutDCAToVertexXY=that.fCutDCAToVertexXY;
286   fCutDCAToVertexZ=that.fCutDCAToVertexZ;
287   fCutMinimalTPCdedx=that.fCutMinimalTPCdedx;
288   fMinimalTPCdedx=that.fMinimalTPCdedx;
289   fParamType=that.fParamType;
290   fParamMix=that.fParamMix;
291
292   fTrack=NULL;
293   fTrackPhi=0.;
294   fTrackPhi=0.;
295   fTrackWeight=0.;
296   fTrackLabel=INT_MIN;
297   fMCevent=NULL;
298   fMCparticle=NULL;
299   fEvent=NULL;
300
301   fESDpid = that.fESDpid;
302   fPIDsource = that.fPIDsource;
303
304   delete fTPCpidCuts;
305   delete fTOFpidCuts;
306   if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
307   if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
308
309   fParticleID=that.fParticleID;
310   fParticleProbability=that.fParticleProbability;
311   memcpy(fProbBayes,that.fProbBayes,sizeof(fProbBayes));
312
313   return *this;
314 }
315
316 //-----------------------------------------------------------------------
317 AliFlowTrackCuts::~AliFlowTrackCuts()
318 {
319   //dtor
320   delete fAliESDtrackCuts;
321   delete fTPCpidCuts;
322   delete fTOFpidCuts;
323 }
324
325 //-----------------------------------------------------------------------
326 void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
327 {
328   //set the event
329   Clear();
330   fEvent=event;
331   fMCevent=mcEvent;
332
333   //do the magic for ESD
334   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(event);
335   if (fCutPID && myESD)
336   {
337     //TODO: maybe call it only for the TOF options?
338     // Added by F. Noferini for TOF PID
339     fESDpid.SetTOFResponse(myESD,AliESDpid::kTOF_T0);
340     fESDpid.MakePID(myESD,kFALSE);
341     // End F. Noferini added part
342   }
343
344   //TODO: AOD
345 }
346
347 //-----------------------------------------------------------------------
348 void AliFlowTrackCuts::SetCutMC( Bool_t b )
349 {
350   //will we be cutting on MC information?
351   fCutMC=b;
352
353   //if we cut on MC info then also the Bethe Bloch should be the one tuned for MC
354   if (fCutMC)
355   {
356     fESDpid.GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
357                                                        1.75295e+01,
358                                                        3.40030e-09,
359                                                        1.96178e+00,
360                                                        3.91720e+00);
361   }
362 }
363
364 //-----------------------------------------------------------------------
365 Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
366 {
367   //check cuts
368   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
369   if (vparticle) return PassesCuts(vparticle);
370   AliFlowTrackSimple* flowtrack = dynamic_cast<AliFlowTrackSimple*>(obj);
371   if (flowtrack) return PassesCuts(flowtrack);
372   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
373   if (tracklets) return PassesCuts(tracklets,id);
374   AliESDPmdTrack* pmdtrack = dynamic_cast<AliESDPmdTrack*>(obj);
375   if (pmdtrack) return PassesPMDcuts(pmdtrack);
376   return kFALSE;  //default when passed wrong type of object
377 }
378
379 //-----------------------------------------------------------------------
380 Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
381 {
382   //check cuts
383   AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
384   if (vparticle) 
385   {
386     return PassesMCcuts(fMCevent,vparticle->GetLabel());
387   }
388   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
389   if (tracklets)
390   {
391     Int_t label0 = tracklets->GetLabel(id,0);
392     Int_t label1 = tracklets->GetLabel(id,1);
393     Int_t label = (label0==label1)?tracklets->GetLabel(id,1):-666;
394     return PassesMCcuts(fMCevent,label);
395   }
396   return kFALSE;  //default when passed wrong type of object
397 }
398
399 //-----------------------------------------------------------------------
400 Bool_t AliFlowTrackCuts::PassesCuts(AliFlowTrackSimple* track)
401 {
402   //check cuts on a flowtracksimple
403
404   //clean up from last iteration
405   fTrack = NULL;
406   return AliFlowTrackSimpleCuts::PassesCuts(track);
407 }
408
409 //-----------------------------------------------------------------------
410 Bool_t AliFlowTrackCuts::PassesCuts(AliMultiplicity* tracklet, Int_t id)
411 {
412   //check cuts on a tracklets
413
414   //clean up from last iteration, and init label
415   fTrack = NULL;
416   fMCparticle=NULL;
417   fTrackLabel=-1;
418
419   fTrackPhi = tracklet->GetPhi(id);
420   fTrackEta = tracklet->GetEta(id);
421   fTrackWeight = 1.0;
422   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
423   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
424
425   //check MC info if available
426   //if the 2 clusters have different label track cannot be good
427   //and should therefore not pass the mc cuts
428   Int_t label0 = tracklet->GetLabel(id,0);
429   Int_t label1 = tracklet->GetLabel(id,1);
430   //if possible get label and mcparticle
431   fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
432   if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
433   if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
434   //check MC cuts
435   if (fCutMC && !PassesMCcuts()) return kFALSE;
436   return kTRUE;
437 }
438
439 //-----------------------------------------------------------------------
440 Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
441 {
442   //check the MC info
443   if (!mcEvent) return kFALSE;
444   if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
445   AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
446   if (!mcparticle) {AliError("no MC track"); return kFALSE;}
447
448   if (fCutMCisPrimary)
449   {
450     if (IsPhysicalPrimary(mcEvent,label,fRequireTransportBitForPrimaries) != fMCisPrimary) return kFALSE;
451   }
452   if (fCutMCPID)
453   {
454     Int_t pdgCode = mcparticle->PdgCode();
455     if (fIgnoreSignInPID) 
456     {
457       if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
458     }
459     else 
460     {
461       if (fMCPID != pdgCode) return kFALSE;
462     }
463   }
464   if ( fCutMCprocessType )
465   {
466     TParticle* particle = mcparticle->Particle();
467     Int_t processID = particle->GetUniqueID();
468     if (processID != fMCprocessType ) return kFALSE;
469   }
470   return kTRUE;
471 }
472 //-----------------------------------------------------------------------
473 Bool_t AliFlowTrackCuts::PassesMCcuts()
474 {
475   if (!fMCevent) return kFALSE;
476   if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
477   fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
478   return PassesMCcuts(fMCevent,fTrackLabel);
479 }
480
481 //-----------------------------------------------------------------------
482 Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
483 {
484   //check cuts for an ESD vparticle
485
486   ////////////////////////////////////////////////////////////////
487   //  start by preparing the track parameters to cut on //////////
488   ////////////////////////////////////////////////////////////////
489   //clean up from last iteration
490   fTrack=NULL; 
491
492   //get the label and the mc particle
493   fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
494   if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
495   else fMCparticle=NULL;
496
497   Bool_t isMCparticle = kFALSE; //some things are different for MC particles, check!
498   AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*>(vparticle);
499   AliAODTrack* aodTrack = NULL;
500   if (esdTrack)
501     //for an ESD track we do some magic sometimes like constructing TPC only parameters
502     //or doing some hybrid, handle that here
503     HandleESDtrack(esdTrack);
504   else
505   {
506     HandleVParticle(vparticle);
507     //now check if produced particle is MC
508     isMCparticle = (dynamic_cast<AliMCParticle*>(fTrack))!=NULL;
509     aodTrack = dynamic_cast<AliAODTrack*>(vparticle); //keep the additional dynamic cast out of the way for ESDs
510   }
511   ////////////////////////////////////////////////////////////////
512   ////////////////////////////////////////////////////////////////
513
514   if (!fTrack) return kFALSE;
515   //because it may be different from global, not needed for aodTrack because we dont do anything funky there
516   if (esdTrack) esdTrack = static_cast<AliESDtrack*>(fTrack);
517   
518   Bool_t pass=kTRUE;
519   //check the common cuts for the current particle fTrack (MC,AOD,ESD)
520   Double_t pt = fTrack->Pt();
521   if (!fFakesAreOK) {if (fTrackLabel<0) pass=kFALSE;}
522   if (fCutPt) {if (pt < fPtMin || pt >= fPtMax ) pass=kFALSE;}
523   if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
524   if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
525   if (fRequireCharge) {if (fTrack->Charge() == 0) pass=kFALSE;}
526   if (fCutCharge && !isMCparticle) {if (fTrack->Charge() != fCharge) pass=kFALSE;}
527   if (fCutCharge && isMCparticle)
528   { 
529     //in case of an MC particle the charge is stored in units of 1/3|e| 
530     Int_t charge = TMath::Nint(fTrack->Charge()/3.0); //mc particles have charge in units of 1/3e
531     if (charge!=fCharge) pass=kFALSE;
532   }
533   //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
534
535   //when additionally MC info is required
536   if (fCutMC && !PassesMCcuts()) pass=kFALSE;
537
538   //the case of ESD or AOD
539   if (esdTrack) { if (!PassesESDcuts(esdTrack)) { pass=kFALSE; } }
540   if (aodTrack) { if (!PassesAODcuts(aodTrack)) { pass=kFALSE; } }
541
542   //true by default, if we didn't set any cuts
543   return pass;
544 }
545
546 //_______________________________________________________________________
547 Bool_t AliFlowTrackCuts::PassesAODcuts(AliAODTrack* track)
548 {
549   Bool_t pass = kTRUE;
550
551   if (fCutNClustersTPC)
552   {
553     Int_t ntpccls = track->GetTPCNcls();
554     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
555   }
556
557   if (fCutNClustersITS)
558   {
559     Int_t nitscls = track->GetITSNcls();
560     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
561   }
562   
563    if (fCutChi2PerClusterTPC)
564   {
565     Double_t chi2tpc = track->Chi2perNDF();
566     if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
567   }
568   
569   if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
570   if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
571   
572   if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
573   
574   if (fCutDCAToVertexXY && track->DCA()>GetMaxDCAToVertexXY()) pass=kFALSE;
575     
576
577   return pass;
578 }
579
580 //_______________________________________________________________________
581 Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
582 {
583   Bool_t pass=kTRUE;
584   if (fIgnoreTPCzRange)
585   {
586     const AliExternalTrackParam* pin = track->GetOuterParam();
587     const AliExternalTrackParam* pout = track->GetInnerParam();
588     if (pin&&pout)
589     {
590       Double_t zin = pin->GetZ();
591       Double_t zout = pout->GetZ();
592       if (zin*zout<0) pass=kFALSE;   //reject if cross the membrane
593       if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
594       if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
595     }
596   }
597  
598   Int_t ntpccls = ( fParamType==kESD_TPConly )?
599                     track->GetTPCNclsIter1():track->GetTPCNcls();    
600   if (fCutChi2PerClusterTPC)
601   {
602     Float_t tpcchi2 = (fParamType==kESD_TPConly)?
603                        track->GetTPCchi2Iter1():track->GetTPCchi2();
604     tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
605     if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
606       pass=kFALSE;
607   }
608
609   if (fCutMinimalTPCdedx) 
610   {
611     if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
612   }
613
614   if (fCutNClustersTPC)
615   {
616     if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
617   }
618
619   Int_t nitscls = track->GetNcls(0);
620   if (fCutNClustersITS)
621   {
622     if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
623   }
624
625   if (fCutPID)
626   {
627     switch (fPIDsource)    
628     {
629       case kTPCpid:
630         if (!PassesTPCpidCut(track)) pass=kFALSE;
631         break;
632       case kTPCdedx:
633         if (!PassesTPCdedxCut(track)) pass=kFALSE;
634         break;
635       case kTOFpid:
636         if (!PassesTOFpidCut(track)) pass=kFALSE;
637         break;
638       case kTOFbeta:
639         if (!PassesTOFbetaCut(track)) pass=kFALSE;
640         break;
641             // part added by F. Noferini
642       case kTOFbayesian:
643               if (!PassesTOFbayesianCut(track)) pass=kFALSE;
644               break;
645             // end part added by F. Noferini
646       default:
647         printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
648         pass=kFALSE;
649         break;
650     }
651   }    
652
653   //some stuff is still handled by AliESDtrackCuts class - delegate
654   if (fAliESDtrackCuts)
655   {
656     if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
657   }
658  
659   return pass;
660 }
661
662 //-----------------------------------------------------------------------
663 void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
664 {
665   //handle the general case
666   switch (fParamType)
667   {
668     default:
669       fTrack = track;
670       break;
671   }
672 }
673
674 //-----------------------------------------------------------------------
675 void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
676 {
677   //handle esd track
678   switch (fParamType)
679   {
680     case kGlobal:
681       fTrack = track;
682       break;
683     case kESD_TPConly:
684       if (!track->FillTPCOnlyTrack(fTPCtrack)) 
685       {
686         fTrack=NULL;
687         fMCparticle=NULL;
688         fTrackLabel=-1;
689         return;
690       }
691       fTrack = &fTPCtrack;
692       //recalculate the label and mc particle, they may differ as TPClabel != global label
693       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
694       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
695       else fMCparticle=NULL;
696       break;
697     default:
698       fTrack = track;
699       break;
700   }
701 }
702
703 //-----------------------------------------------------------------------
704 Int_t AliFlowTrackCuts::Count(AliVEvent* event)
705 {
706   //calculate the number of track in given event.
707   //if argument is NULL(default) take the event attached
708   //by SetEvent()
709   Int_t multiplicity = 0;
710   if (!event)
711   {
712     for (Int_t i=0; i<GetNumberOfInputObjects(); i++)
713     {
714       if (IsSelected(GetInputObject(i))) multiplicity++;
715     }
716   }
717   else
718   {
719     for (Int_t i=0; i<event->GetNumberOfTracks(); i++)
720     {
721       if (IsSelected(event->GetTrack(i))) multiplicity++;
722     }
723   }
724   return multiplicity;
725 }
726
727 //-----------------------------------------------------------------------
728 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()
729 {
730   //get standard cuts
731   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard Global flow cuts");
732   cuts->SetParamType(kGlobal);
733   cuts->SetPtRange(0.2,5.);
734   cuts->SetEtaRange(-0.8,0.8);
735   cuts->SetMinNClustersTPC(70);
736   cuts->SetMinChi2PerClusterTPC(0.1);
737   cuts->SetMaxChi2PerClusterTPC(4.0);
738   cuts->SetMinNClustersITS(2);
739   cuts->SetRequireITSRefit(kTRUE);
740   cuts->SetRequireTPCRefit(kTRUE);
741   cuts->SetMaxDCAToVertexXY(0.3);
742   cuts->SetMaxDCAToVertexZ(0.3);
743   cuts->SetAcceptKinkDaughters(kFALSE);
744   cuts->SetMinimalTPCdedx(10.);
745
746   return cuts;
747 }
748
749 //-----------------------------------------------------------------------
750 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts2010()
751 {
752   //get standard cuts
753   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly flow cuts");
754   cuts->SetParamType(kESD_TPConly);
755   cuts->SetPtRange(0.2,5.);
756   cuts->SetEtaRange(-0.8,0.8);
757   cuts->SetMinNClustersTPC(70);
758   cuts->SetMinChi2PerClusterTPC(0.2);
759   cuts->SetMaxChi2PerClusterTPC(4.0);
760   cuts->SetMaxDCAToVertexXY(3.0);
761   cuts->SetMaxDCAToVertexZ(3.0);
762   cuts->SetDCAToVertex2D(kTRUE);
763   cuts->SetAcceptKinkDaughters(kFALSE);
764   cuts->SetMinimalTPCdedx(10.);
765
766   return cuts;
767 }
768
769 //-----------------------------------------------------------------------
770 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts()
771 {
772   //get standard cuts
773   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly flow cuts");
774   cuts->SetParamType(kESD_TPConly);
775   cuts->SetPtRange(0.2,5.);
776   cuts->SetEtaRange(-0.8,0.8);
777   cuts->SetMinNClustersTPC(70);
778   cuts->SetMinChi2PerClusterTPC(0.2);
779   cuts->SetMaxChi2PerClusterTPC(4.0);
780   cuts->SetMaxDCAToVertexXY(3.0);
781   cuts->SetMaxDCAToVertexZ(3.0);
782   cuts->SetDCAToVertex2D(kTRUE);
783   cuts->SetAcceptKinkDaughters(kFALSE);
784   cuts->SetMinimalTPCdedx(10.);
785
786   return cuts;
787 }
788
789 //-----------------------------------------------------------------------
790 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
791 {
792   //get standard cuts
793   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
794   delete cuts->fAliESDtrackCuts;
795   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
796   cuts->SetParamType(kGlobal);
797   return cuts;
798 }
799
800 //-----------------------------------------------------------------------
801 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackSPDtracklet() const
802 {
803   //make a flow track from tracklet
804   AliFlowTrack* flowtrack=NULL;
805   TParticle *tmpTParticle=NULL;
806   AliMCParticle* tmpAliMCParticle=NULL;
807   switch (fParamMix)
808   {
809     case kPure:
810       flowtrack = new AliFlowTrack();
811       flowtrack->SetPhi(fTrackPhi);
812       flowtrack->SetEta(fTrackEta);
813       break;
814     case kTrackWithMCkine:
815       if (!fMCparticle) return NULL;
816       flowtrack = new AliFlowTrack();
817       flowtrack->SetPhi( fMCparticle->Phi() );
818       flowtrack->SetEta( fMCparticle->Eta() );
819       flowtrack->SetPt( fMCparticle->Pt() );
820       break;
821     case kTrackWithMCpt:
822       if (!fMCparticle) return NULL;
823       flowtrack = new AliFlowTrack();
824       flowtrack->SetPhi(fTrackPhi);
825       flowtrack->SetEta(fTrackEta);
826       flowtrack->SetPt(fMCparticle->Pt());
827       break;
828     case kTrackWithPtFromFirstMother:
829       if (!fMCparticle) return NULL;
830       flowtrack = new AliFlowTrack();
831       flowtrack->SetPhi(fTrackPhi);
832       flowtrack->SetEta(fTrackEta);
833       tmpTParticle = fMCparticle->Particle();
834       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
835       flowtrack->SetPt(tmpAliMCParticle->Pt());
836       break;
837     default:
838       flowtrack = new AliFlowTrack();
839       flowtrack->SetPhi(fTrackPhi);
840       flowtrack->SetEta(fTrackEta);
841       break;
842   }
843   flowtrack->SetSource(AliFlowTrack::kFromTracklet);
844   return flowtrack;
845 }
846
847 //-----------------------------------------------------------------------
848 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVParticle() const
849 {
850   //make flow track from AliVParticle (ESD,AOD,MC)
851   if (!fTrack) return NULL;
852   AliFlowTrack* flowtrack=NULL;
853   TParticle *tmpTParticle=NULL;
854   AliMCParticle* tmpAliMCParticle=NULL;
855   switch(fParamMix)
856   {
857     case kPure:
858       flowtrack = new AliFlowTrack(fTrack);
859       break;
860     case kTrackWithMCkine:
861       flowtrack = new AliFlowTrack(fMCparticle);
862       break;
863     case kTrackWithMCPID:
864       flowtrack = new AliFlowTrack(fTrack);
865       //flowtrack->setPID(...) from mc, when implemented
866       break;
867     case kTrackWithMCpt:
868       if (!fMCparticle) return NULL;
869       flowtrack = new AliFlowTrack(fTrack);
870       flowtrack->SetPt(fMCparticle->Pt());
871       break;
872     case kTrackWithPtFromFirstMother:
873       if (!fMCparticle) return NULL;
874       flowtrack = new AliFlowTrack(fTrack);
875       tmpTParticle = fMCparticle->Particle();
876       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
877       flowtrack->SetPt(tmpAliMCParticle->Pt());
878       break;
879     default:
880       flowtrack = new AliFlowTrack(fTrack);
881       break;
882   }
883   if (fParamType==kMC) flowtrack->SetSource(AliFlowTrack::kFromMC);
884   else if (dynamic_cast<AliESDtrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromESD);
885   else if (dynamic_cast<AliAODTrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromAOD);
886   else if (dynamic_cast<AliMCParticle*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromMC);
887   return flowtrack;
888 }
889
890 //-----------------------------------------------------------------------
891 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackPMDtrack() const
892 {
893   //make a flow track from PMD track
894   AliFlowTrack* flowtrack=NULL;
895   TParticle *tmpTParticle=NULL;
896   AliMCParticle* tmpAliMCParticle=NULL;
897   switch (fParamMix)
898   {
899     case kPure:
900       flowtrack = new AliFlowTrack();
901       flowtrack->SetPhi(fTrackPhi);
902       flowtrack->SetEta(fTrackEta);
903       flowtrack->SetWeight(fTrackWeight);
904       break;
905     case kTrackWithMCkine:
906       if (!fMCparticle) return NULL;
907       flowtrack = new AliFlowTrack();
908       flowtrack->SetPhi( fMCparticle->Phi() );
909       flowtrack->SetEta( fMCparticle->Eta() );
910       flowtrack->SetWeight(fTrackWeight);
911       flowtrack->SetPt( fMCparticle->Pt() );
912       break;
913     case kTrackWithMCpt:
914       if (!fMCparticle) return NULL;
915       flowtrack = new AliFlowTrack();
916       flowtrack->SetPhi(fTrackPhi);
917       flowtrack->SetEta(fTrackEta);
918       flowtrack->SetWeight(fTrackWeight);
919       flowtrack->SetPt(fMCparticle->Pt());
920       break;
921     case kTrackWithPtFromFirstMother:
922       if (!fMCparticle) return NULL;
923       flowtrack = new AliFlowTrack();
924       flowtrack->SetPhi(fTrackPhi);
925       flowtrack->SetEta(fTrackEta);
926       flowtrack->SetWeight(fTrackWeight);
927       tmpTParticle = fMCparticle->Particle();
928       tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
929       flowtrack->SetPt(tmpAliMCParticle->Pt());
930       break;
931     default:
932       flowtrack = new AliFlowTrack();
933       flowtrack->SetPhi(fTrackPhi);
934       flowtrack->SetEta(fTrackEta);
935       flowtrack->SetWeight(fTrackWeight);
936       break;
937   }
938
939   flowtrack->SetSource(AliFlowTrack::kFromPMD);
940   return flowtrack;
941 }
942
943 //-----------------------------------------------------------------------
944 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
945 {
946   //get a flow track constructed from whatever we applied cuts on
947   //caller is resposible for deletion
948   //if construction fails return NULL
949   switch (fParamType)
950   {
951     case kESD_SPDtracklet:
952       return MakeFlowTrackSPDtracklet();
953     case kPMD:
954       return MakeFlowTrackPMDtrack();
955     default:
956       return MakeFlowTrackVParticle();
957   }
958 }
959
960 //-----------------------------------------------------------------------
961 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
962 {
963   //check if current particle is a physical primary
964   if (!fMCevent) return kFALSE;
965   if (fTrackLabel<0) return kFALSE;
966   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
967 }
968
969 //-----------------------------------------------------------------------
970 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
971 {
972   //check if current particle is a physical primary
973   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
974   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
975   if (!track) return kFALSE;
976   TParticle* particle = track->Particle();
977   Bool_t transported = particle->TestBit(kTransportBit);
978   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
979         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
980   return (physprim && (transported || !requiretransported));
981 }
982
983 //-----------------------------------------------------------------------
984 void AliFlowTrackCuts::DefineHistograms()
985 {
986   //define qa histograms
987 }
988
989 //-----------------------------------------------------------------------
990 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
991 {
992   //get the number of tracks in the input event according source
993   //selection (ESD tracks, tracklets, MC particles etc.)
994   AliESDEvent* esd=NULL;
995   switch (fParamType)
996   {
997     case kESD_SPDtracklet:
998       esd = dynamic_cast<AliESDEvent*>(fEvent);
999       if (!esd) return 0;
1000       return esd->GetMultiplicity()->GetNumberOfTracklets();
1001     case kMC:
1002       if (!fMCevent) return 0;
1003       return fMCevent->GetNumberOfTracks();
1004     case kPMD:
1005       esd = dynamic_cast<AliESDEvent*>(fEvent);
1006       if (!esd) return 0;
1007       return esd->GetNumberOfPmdTracks();
1008     default:
1009       if (!fEvent) return 0;
1010       return fEvent->GetNumberOfTracks();
1011   }
1012   return 0;
1013 }
1014
1015 //-----------------------------------------------------------------------
1016 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
1017 {
1018   //get the input object according the data source selection:
1019   //(esd tracks, traclets, mc particles,etc...)
1020   AliESDEvent* esd=NULL;
1021   switch (fParamType)
1022   {
1023     case kESD_SPDtracklet:
1024       esd = dynamic_cast<AliESDEvent*>(fEvent);
1025       if (!esd) return NULL;
1026       return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
1027     case kMC:
1028       if (!fMCevent) return NULL;
1029       return fMCevent->GetTrack(i);
1030     case kPMD:
1031       esd = dynamic_cast<AliESDEvent*>(fEvent);
1032       if (!esd) return NULL;
1033       return esd->GetPmdTrack(i);
1034     default:
1035       if (!fEvent) return NULL;
1036       return fEvent->GetTrack(i);
1037   }
1038 }
1039
1040 //-----------------------------------------------------------------------
1041 void AliFlowTrackCuts::Clear(Option_t*)
1042 {
1043   //clean up
1044   fTrack=NULL;
1045   fMCevent=NULL;
1046   fMCparticle=NULL;
1047   fTrackLabel=0;
1048   fTrackWeight=0.0;
1049   fTrackEta=0.0;
1050   fTrackPhi=0.0;
1051 }
1052
1053 //-----------------------------------------------------------------------
1054 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(AliESDtrack* track )
1055 {
1056   //check if passes PID cut using timing in TOF
1057   Bool_t goodtrack = (track) && 
1058                      (track->GetStatus() & AliESDtrack::kTOFpid) && 
1059                      (track->GetTOFsignal() > 12000) && 
1060                      (track->GetTOFsignal() < 100000) && 
1061                      (track->GetIntegratedLength() > 365) && 
1062                     !(track->GetStatus() & AliESDtrack::kTOFmismatch);
1063
1064   if (!goodtrack) return kFALSE;
1065   
1066   const Float_t c = 2.99792457999999984e-02;  
1067   Float_t p = track->GetP();
1068   Float_t L = track->GetIntegratedLength();  
1069   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
1070   Float_t timeTOF = track->GetTOFsignal()- trackT0; 
1071   Float_t beta = L/timeTOF/c;
1072   Double_t integratedTimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
1073   track->GetIntegratedTimes(integratedTimes);
1074
1075   //construct the pid index because it's not AliPID::EParticleType
1076   Int_t pid = 0;
1077   switch (fParticleID)
1078   {
1079     case AliPID::kPion:
1080       pid=2;
1081       break;
1082     case AliPID::kKaon:
1083       pid=3;
1084       break;
1085     case AliPID::kProton:
1086       pid=4;
1087       break;
1088     default:
1089       return kFALSE;
1090   }
1091
1092   //signal to cut on
1093   Float_t s = beta-L/integratedTimes[pid]/c;
1094
1095   Float_t* arr = fTOFpidCuts->GetMatrixArray();
1096   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
1097   if (col<0) return kFALSE;
1098   Float_t min = (*fTOFpidCuts)(1,col);
1099   Float_t max = (*fTOFpidCuts)(2,col);
1100
1101   //printf("--------------TOF beta cut %s\n",(s>min && s<max)?"PASS":"FAIL");
1102   return (s>min && s<max);
1103 }
1104
1105 //-----------------------------------------------------------------------
1106 Bool_t AliFlowTrackCuts::PassesTOFpidCut(AliESDtrack* track)
1107 {
1108   //check if passes PID cut using default TOF pid
1109   Double_t pidTOF[AliPID::kSPECIES];
1110   track->GetTOFpid(pidTOF);
1111   if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
1112   return kFALSE;
1113 }
1114
1115 //-----------------------------------------------------------------------
1116 Bool_t AliFlowTrackCuts::PassesTPCpidCut(AliESDtrack* track)
1117 {
1118   //check if passes PID cut using default TPC pid
1119   Double_t pidTPC[AliPID::kSPECIES];
1120   track->GetTPCpid(pidTPC);
1121   Double_t probablity = 0.;
1122   switch (fParticleID)
1123   {
1124     case AliPID::kPion:
1125       probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
1126       break;
1127     default:
1128       probablity = pidTPC[fParticleID];
1129   }
1130   if (probablity >= fParticleProbability) return kTRUE;
1131   return kFALSE;
1132 }
1133
1134 //-----------------------------------------------------------------------
1135 Bool_t AliFlowTrackCuts::PassesTPCdedxCut(AliESDtrack* track)
1136 {
1137   //check if passes PID cut using dedx signal in the TPC
1138   if (!fTPCpidCuts)
1139   {
1140     printf("no TPCpidCuts\n");
1141     return kFALSE;
1142   }
1143
1144   const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
1145   if (!tpcparam) return kFALSE;
1146   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(tpcparam->GetP(), fParticleID);
1147   Float_t sigTPC = track->GetTPCsignal();
1148   Float_t s = (sigTPC-sigExp)/sigExp;
1149   Double_t pt = track->Pt();
1150
1151   Float_t* arr = fTPCpidCuts->GetMatrixArray();
1152   Int_t col = TMath::BinarySearch(fTPCpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
1153   if (col<0) return kFALSE;
1154   Float_t min = (*fTPCpidCuts)(1,col);
1155   Float_t max = (*fTPCpidCuts)(2,col);
1156
1157   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
1158   return (s>min && s<max);
1159 }
1160
1161 //-----------------------------------------------------------------------
1162 void AliFlowTrackCuts::InitPIDcuts()
1163 {
1164   //init matrices with PID cuts
1165   TMatrixF* t = NULL;
1166   if (!fTPCpidCuts)
1167   {
1168     if (fParticleID==AliPID::kPion)
1169     {
1170       t = new TMatrixF(3,10);
1171       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.2;
1172       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.2;
1173       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.25;
1174       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.25;
1175       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
1176       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
1177       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
1178       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
1179       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
1180       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  =     0;
1181     }
1182     else
1183     if (fParticleID==AliPID::kKaon)
1184     {
1185       t = new TMatrixF(3,7);
1186       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.4; 
1187       (*t)(0,1)  = 0.25;  (*t)(1,1)  =-0.15;  (*t)(2,1)  = 0.4;
1188       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.1;  (*t)(2,2)  = 0.4;
1189       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.1;  (*t)(2,3)  = 0.4;
1190       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.6;
1191       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.6;
1192       (*t)(0,6)  = 0.50;  (*t)(1,6)  =    0;  (*t)(2,6)  =   0;
1193     }
1194     else
1195     if (fParticleID==AliPID::kProton)
1196     {
1197       t = new TMatrixF(3,16);
1198       (*t)(0,0)  = 0.20;  (*t)(1,0)  =     0;  (*t)(2,0)  =    0; 
1199       (*t)(0,1)  = 0.25;  (*t)(1,1)  =  -0.2;  (*t)(2,1)  =  0.3; 
1200       (*t)(0,2)  = 0.30;  (*t)(1,2)  =  -0.2;  (*t)(2,2)  =  0.6; 
1201       (*t)(0,3)  = 0.35;  (*t)(1,3)  =  -0.2;  (*t)(2,3)  =  0.6; 
1202       (*t)(0,4)  = 0.40;  (*t)(1,4)  =  -0.2;  (*t)(2,4)  =  0.6; 
1203       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.15;  (*t)(2,5)  =  0.6; 
1204       (*t)(0,6)  = 0.50;  (*t)(1,6)  =  -0.1;  (*t)(2,6)  =  0.6; 
1205       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =  0.6; 
1206       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.05;  (*t)(2,8)  = 0.45; 
1207       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.05;  (*t)(2,9)  = 0.45; 
1208       (*t)(0,10) = 0.70;  (*t)(1,10) = -0.05;  (*t)(2,10) = 0.45; 
1209       (*t)(0,11) = 0.75;  (*t)(1,11) = -0.05;  (*t)(2,11) = 0.45; 
1210       (*t)(0,12) = 0.80;  (*t)(1,12) =     0;  (*t)(2,12) = 0.45; 
1211       (*t)(0,13) = 0.85;  (*t)(1,13) =     0;  (*t)(2,13) = 0.45; 
1212       (*t)(0,14) = 0.90;  (*t)(1,14) =     0;  (*t)(2,14) = 0.45;
1213       (*t)(0,15) = 0.95;  (*t)(1,15) =     0;  (*t)(2,15) =    0;
1214     }
1215     fTPCpidCuts=t;
1216   }
1217   t = NULL;
1218   if (!fTOFpidCuts)
1219   {
1220     if (fParticleID==AliPID::kPion)
1221     {
1222       //TOF pions, 0.9 purity
1223       t = new TMatrixF(3,61);
1224       (*t)(0,0)  = 0.000;  (*t)(2,0)  = 0.000;  (*t)(2,0)  =   0.000;
1225       (*t)(0,1)  = 0.050;  (*t)(2,1)  = 0.000;  (*t)(2,1)  =   0.000;
1226       (*t)(0,2)  = 0.100;  (*t)(2,2)  = 0.000;  (*t)(2,2)  =   0.000;
1227       (*t)(0,3)  = 0.150;  (*t)(2,3)  = 0.000;  (*t)(2,3)  =   0.000;
1228       (*t)(0,4)  = 0.200;  (*t)(2,4)  = 0.000;  (*t)(2,4)  =   0.000;
1229       (*t)(0,5)  = 0.250;  (*t)(2,5)  = -0.046;  (*t)(2,5)  =   0.046;
1230       (*t)(0,6)  = 0.300;  (*t)(2,6)  = -0.038;  (*t)(2,6)  =   0.038;
1231       (*t)(0,7)  = 0.350;  (*t)(2,7)  = -0.034;  (*t)(2,7)  =   0.034;
1232       (*t)(0,8)  = 0.400;  (*t)(2,8)  = -0.032;  (*t)(2,8)  =   0.032;
1233       (*t)(0,9)  = 0.450;  (*t)(2,9)  = -0.030;  (*t)(2,9)  =   0.030;
1234       (*t)(0,10)  = 0.500;  (*t)(2,10)  = -0.030;  (*t)(2,10)  =   0.030;
1235       (*t)(0,11)  = 0.550;  (*t)(2,11)  = -0.030;  (*t)(2,11)  =   0.030;
1236       (*t)(0,12)  = 0.600;  (*t)(2,12)  = -0.030;  (*t)(2,12)  =   0.030;
1237       (*t)(0,13)  = 0.650;  (*t)(2,13)  = -0.030;  (*t)(2,13)  =   0.030;
1238       (*t)(0,14)  = 0.700;  (*t)(2,14)  = -0.030;  (*t)(2,14)  =   0.030;
1239       (*t)(0,15)  = 0.750;  (*t)(2,15)  = -0.030;  (*t)(2,15)  =   0.030;
1240       (*t)(0,16)  = 0.800;  (*t)(2,16)  = -0.030;  (*t)(2,16)  =   0.030;
1241       (*t)(0,17)  = 0.850;  (*t)(2,17)  = -0.030;  (*t)(2,17)  =   0.030;
1242       (*t)(0,18)  = 0.900;  (*t)(2,18)  = -0.030;  (*t)(2,18)  =   0.030;
1243       (*t)(0,19)  = 0.950;  (*t)(2,19)  = -0.028;  (*t)(2,19)  =   0.028;
1244       (*t)(0,20)  = 1.000;  (*t)(2,20)  = -0.028;  (*t)(2,20)  =   0.028;
1245       (*t)(0,21)  = 1.100;  (*t)(2,21)  = -0.028;  (*t)(2,21)  =   0.028;
1246       (*t)(0,22)  = 1.200;  (*t)(2,22)  = -0.026;  (*t)(2,22)  =   0.028;
1247       (*t)(0,23)  = 1.300;  (*t)(2,23)  = -0.024;  (*t)(2,23)  =   0.028;
1248       (*t)(0,24)  = 1.400;  (*t)(2,24)  = -0.020;  (*t)(2,24)  =   0.028;
1249       (*t)(0,25)  = 1.500;  (*t)(2,25)  = -0.018;  (*t)(2,25)  =   0.028;
1250       (*t)(0,26)  = 1.600;  (*t)(2,26)  = -0.016;  (*t)(2,26)  =   0.028;
1251       (*t)(0,27)  = 1.700;  (*t)(2,27)  = -0.014;  (*t)(2,27)  =   0.028;
1252       (*t)(0,28)  = 1.800;  (*t)(2,28)  = -0.012;  (*t)(2,28)  =   0.026;
1253       (*t)(0,29)  = 1.900;  (*t)(2,29)  = -0.010;  (*t)(2,29)  =   0.026;
1254       (*t)(0,30)  = 2.000;  (*t)(2,30)  = -0.008;  (*t)(2,30)  =   0.026;
1255       (*t)(0,31)  = 2.100;  (*t)(2,31)  = -0.008;  (*t)(2,31)  =   0.024;
1256       (*t)(0,32)  = 2.200;  (*t)(2,32)  = -0.006;  (*t)(2,32)  =   0.024;
1257       (*t)(0,33)  = 2.300;  (*t)(2,33)  = -0.004;  (*t)(2,33)  =   0.024;
1258       (*t)(0,34)  = 2.400;  (*t)(2,34)  = -0.004;  (*t)(2,34)  =   0.024;
1259       (*t)(0,35)  = 2.500;  (*t)(2,35)  = -0.002;  (*t)(2,35)  =   0.024;
1260       (*t)(0,36)  = 2.600;  (*t)(2,36)  = -0.002;  (*t)(2,36)  =   0.024;
1261       (*t)(0,37)  = 2.700;  (*t)(2,37)  = 0.000;  (*t)(2,37)  =   0.024;
1262       (*t)(0,38)  = 2.800;  (*t)(2,38)  = 0.000;  (*t)(2,38)  =   0.026;
1263       (*t)(0,39)  = 2.900;  (*t)(2,39)  = 0.000;  (*t)(2,39)  =   0.024;
1264       (*t)(0,40)  = 3.000;  (*t)(2,40)  = 0.002;  (*t)(2,40)  =   0.026;
1265       (*t)(0,41)  = 3.100;  (*t)(2,41)  = 0.002;  (*t)(2,41)  =   0.026;
1266       (*t)(0,42)  = 3.200;  (*t)(2,42)  = 0.002;  (*t)(2,42)  =   0.026;
1267       (*t)(0,43)  = 3.300;  (*t)(2,43)  = 0.002;  (*t)(2,43)  =   0.026;
1268       (*t)(0,44)  = 3.400;  (*t)(2,44)  = 0.002;  (*t)(2,44)  =   0.026;
1269       (*t)(0,45)  = 3.500;  (*t)(2,45)  = 0.002;  (*t)(2,45)  =   0.026;
1270       (*t)(0,46)  = 3.600;  (*t)(2,46)  = 0.002;  (*t)(2,46)  =   0.026;
1271       (*t)(0,47)  = 3.700;  (*t)(2,47)  = 0.002;  (*t)(2,47)  =   0.026;
1272       (*t)(0,48)  = 3.800;  (*t)(2,48)  = 0.002;  (*t)(2,48)  =   0.026;
1273       (*t)(0,49)  = 3.900;  (*t)(2,49)  = 0.004;  (*t)(2,49)  =   0.024;
1274       (*t)(0,50)  = 4.000;  (*t)(2,50)  = 0.004;  (*t)(2,50)  =   0.026;
1275       (*t)(0,51)  = 4.100;  (*t)(2,51)  = 0.004;  (*t)(2,51)  =   0.026;
1276       (*t)(0,52)  = 4.200;  (*t)(2,52)  = 0.004;  (*t)(2,52)  =   0.024;
1277       (*t)(0,53)  = 4.300;  (*t)(2,53)  = 0.006;  (*t)(2,53)  =   0.024;
1278       (*t)(0,54)  = 4.400;  (*t)(2,54)  = 0.000;  (*t)(2,54)  =   0.000;
1279       (*t)(0,55)  = 4.500;  (*t)(2,55)  = 0.000;  (*t)(2,55)  =   0.000;
1280       (*t)(0,56)  = 4.600;  (*t)(2,56)  = 0.000;  (*t)(2,56)  =   0.000;
1281       (*t)(0,57)  = 4.700;  (*t)(2,57)  = 0.000;  (*t)(2,57)  =   0.000;
1282       (*t)(0,58)  = 4.800;  (*t)(2,58)  = 0.000;  (*t)(2,58)  =   0.000;
1283       (*t)(0,59)  = 4.900;  (*t)(2,59)  = 0.000;  (*t)(2,59)  =   0.000;
1284       (*t)(0,60)  = 5.900;  (*t)(2,60)  = 0.000;  (*t)(2,60)  =   0.000;
1285     }
1286     else
1287     if (fParticleID==AliPID::kProton)
1288     {
1289       //TOF protons, 0.9 purity
1290       t = new TMatrixF(3,61);
1291       (*t)(0,0)  = 0.000;  (*t)(2,0)  = 0.000;  (*t)(2,0)  =   0.000;
1292       (*t)(0,1)  = 0.050;  (*t)(2,1)  = 0.000;  (*t)(2,1)  =   0.000;
1293       (*t)(0,2)  = 0.100;  (*t)(2,2)  = 0.000;  (*t)(2,2)  =   0.000;
1294       (*t)(0,3)  = 0.150;  (*t)(2,3)  = 0.000;  (*t)(2,3)  =   0.000;
1295       (*t)(0,4)  = 0.200;  (*t)(2,4)  = 0.000;  (*t)(2,4)  =   0.000;
1296       (*t)(0,5)  = 0.250;  (*t)(2,5)  = 0.000;  (*t)(2,5)  =   0.000;
1297       (*t)(0,6)  = 0.300;  (*t)(2,6)  = 0.000;  (*t)(2,6)  =   0.000;
1298       (*t)(0,7)  = 0.350;  (*t)(2,7)  = 0.000;  (*t)(2,7)  =   0.000;
1299       (*t)(0,8)  = 0.400;  (*t)(2,8)  = 0.000;  (*t)(2,8)  =   0.000;
1300       (*t)(0,9)  = 0.450;  (*t)(2,9)  = 0.000;  (*t)(2,9)  =   0.000;
1301       (*t)(0,10)  = 0.500;  (*t)(2,10)  = 0.000;  (*t)(2,10)  =   0.000;
1302       (*t)(0,11)  = 0.550;  (*t)(2,11)  = 0.000;  (*t)(2,11)  =   0.000;
1303       (*t)(0,12)  = 0.600;  (*t)(2,12)  = 0.000;  (*t)(2,12)  =   0.000;
1304       (*t)(0,13)  = 0.650;  (*t)(2,13)  = 0.000;  (*t)(2,13)  =   0.000;
1305       (*t)(0,14)  = 0.700;  (*t)(2,14)  = 0.000;  (*t)(2,14)  =   0.000;
1306       (*t)(0,15)  = 0.750;  (*t)(2,15)  = 0.000;  (*t)(2,15)  =   0.000;
1307       (*t)(0,16)  = 0.800;  (*t)(2,16)  = 0.000;  (*t)(2,16)  =   0.000;
1308       (*t)(0,17)  = 0.850;  (*t)(2,17)  = -0.070;  (*t)(2,17)  =   0.070;
1309       (*t)(0,18)  = 0.900;  (*t)(2,18)  = -0.072;  (*t)(2,18)  =   0.072;
1310       (*t)(0,19)  = 0.950;  (*t)(2,19)  = -0.072;  (*t)(2,19)  =   0.072;
1311       (*t)(0,20)  = 1.000;  (*t)(2,20)  = -0.074;  (*t)(2,20)  =   0.074;
1312       (*t)(0,21)  = 1.100;  (*t)(2,21)  = -0.032;  (*t)(2,21)  =   0.032;
1313       (*t)(0,22)  = 1.200;  (*t)(2,22)  = -0.026;  (*t)(2,22)  =   0.026;
1314       (*t)(0,23)  = 1.300;  (*t)(2,23)  = -0.026;  (*t)(2,23)  =   0.026;
1315       (*t)(0,24)  = 1.400;  (*t)(2,24)  = -0.024;  (*t)(2,24)  =   0.024;
1316       (*t)(0,25)  = 1.500;  (*t)(2,25)  = -0.024;  (*t)(2,25)  =   0.024;
1317       (*t)(0,26)  = 1.600;  (*t)(2,26)  = -0.026;  (*t)(2,26)  =   0.026;
1318       (*t)(0,27)  = 1.700;  (*t)(2,27)  = -0.026;  (*t)(2,27)  =   0.026;
1319       (*t)(0,28)  = 1.800;  (*t)(2,28)  = -0.026;  (*t)(2,28)  =   0.026;
1320       (*t)(0,29)  = 1.900;  (*t)(2,29)  = -0.026;  (*t)(2,29)  =   0.026;
1321       (*t)(0,30)  = 2.000;  (*t)(2,30)  = -0.026;  (*t)(2,30)  =   0.026;
1322       (*t)(0,31)  = 2.100;  (*t)(2,31)  = -0.026;  (*t)(2,31)  =   0.026;
1323       (*t)(0,32)  = 2.200;  (*t)(2,32)  = -0.026;  (*t)(2,32)  =   0.024;
1324       (*t)(0,33)  = 2.300;  (*t)(2,33)  = -0.028;  (*t)(2,33)  =   0.022;
1325       (*t)(0,34)  = 2.400;  (*t)(2,34)  = -0.028;  (*t)(2,34)  =   0.020;
1326       (*t)(0,35)  = 2.500;  (*t)(2,35)  = -0.028;  (*t)(2,35)  =   0.018;
1327       (*t)(0,36)  = 2.600;  (*t)(2,36)  = -0.028;  (*t)(2,36)  =   0.016;
1328       (*t)(0,37)  = 2.700;  (*t)(2,37)  = -0.028;  (*t)(2,37)  =   0.016;
1329       (*t)(0,38)  = 2.800;  (*t)(2,38)  = -0.030;  (*t)(2,38)  =   0.014;
1330       (*t)(0,39)  = 2.900;  (*t)(2,39)  = -0.030;  (*t)(2,39)  =   0.012;
1331       (*t)(0,40)  = 3.000;  (*t)(2,40)  = -0.030;  (*t)(2,40)  =   0.012;
1332       (*t)(0,41)  = 3.100;  (*t)(2,41)  = -0.030;  (*t)(2,41)  =   0.010;
1333       (*t)(0,42)  = 3.200;  (*t)(2,42)  = -0.030;  (*t)(2,42)  =   0.010;
1334       (*t)(0,43)  = 3.300;  (*t)(2,43)  = -0.030;  (*t)(2,43)  =   0.010;
1335       (*t)(0,44)  = 3.400;  (*t)(2,44)  = -0.030;  (*t)(2,44)  =   0.008;
1336       (*t)(0,45)  = 3.500;  (*t)(2,45)  = -0.030;  (*t)(2,45)  =   0.008;
1337       (*t)(0,46)  = 3.600;  (*t)(2,46)  = -0.030;  (*t)(2,46)  =   0.008;
1338       (*t)(0,47)  = 3.700;  (*t)(2,47)  = -0.030;  (*t)(2,47)  =   0.006;
1339       (*t)(0,48)  = 3.800;  (*t)(2,48)  = -0.030;  (*t)(2,48)  =   0.006;
1340       (*t)(0,49)  = 3.900;  (*t)(2,49)  = -0.030;  (*t)(2,49)  =   0.006;
1341       (*t)(0,50)  = 4.000;  (*t)(2,50)  = -0.028;  (*t)(2,50)  =   0.004;
1342       (*t)(0,51)  = 4.100;  (*t)(2,51)  = -0.030;  (*t)(2,51)  =   0.004;
1343       (*t)(0,52)  = 4.200;  (*t)(2,52)  = -0.030;  (*t)(2,52)  =   0.004;
1344       (*t)(0,53)  = 4.300;  (*t)(2,53)  = -0.028;  (*t)(2,53)  =   0.002;
1345       (*t)(0,54)  = 4.400;  (*t)(2,54)  = -0.030;  (*t)(2,54)  =   0.002;
1346       (*t)(0,55)  = 4.500;  (*t)(2,55)  = -0.028;  (*t)(2,55)  =   0.002;
1347       (*t)(0,56)  = 4.600;  (*t)(2,56)  = -0.028;  (*t)(2,56)  =   0.002;
1348       (*t)(0,57)  = 4.700;  (*t)(2,57)  = -0.028;  (*t)(2,57)  =   0.000;
1349       (*t)(0,58)  = 4.800;  (*t)(2,58)  = -0.028;  (*t)(2,58)  =   0.002;
1350       (*t)(0,59)  = 4.900;  (*t)(2,59)  = 0.000;  (*t)(2,59)  =   0.000;
1351       (*t)(0,60)  = 5.900;  (*t)(2,60)  = 0.000;  (*t)(2,60)  =   0.000; 
1352     }
1353     else
1354     if (fParticleID==AliPID::kKaon)
1355     {
1356       //TOF kaons, 0.9 purity
1357       t = new TMatrixF(3,61);
1358       (*t)(0,0)  = 0.000;  (*t)(2,0)  = 0.000;  (*t)(2,0)  =   0.000;
1359       (*t)(0,1)  = 0.050;  (*t)(2,1)  = 0.000;  (*t)(2,1)  =   0.000;
1360       (*t)(0,2)  = 0.100;  (*t)(2,2)  = 0.000;  (*t)(2,2)  =   0.000;
1361       (*t)(0,3)  = 0.150;  (*t)(2,3)  = 0.000;  (*t)(2,3)  =   0.000;
1362       (*t)(0,4)  = 0.200;  (*t)(2,4)  = 0.000;  (*t)(2,4)  =   0.000;
1363       (*t)(0,5)  = 0.250;  (*t)(2,5)  = 0.000;  (*t)(2,5)  =   0.000;
1364       (*t)(0,6)  = 0.300;  (*t)(2,6)  = 0.000;  (*t)(2,6)  =   0.000;
1365       (*t)(0,7)  = 0.350;  (*t)(2,7)  = 0.000;  (*t)(2,7)  =   0.000;
1366       (*t)(0,8)  = 0.400;  (*t)(2,8)  = 0.000;  (*t)(2,8)  =   0.000;
1367       (*t)(0,9)  = 0.450;  (*t)(2,9)  = 0.000;  (*t)(2,9)  =   0.000;
1368       (*t)(0,10)  = 0.500;  (*t)(2,10)  = 0.000;  (*t)(2,10)  =   0.000;
1369       (*t)(0,11)  = 0.550;  (*t)(2,11)  = -0.026;  (*t)(2,11)  =   0.026;
1370       (*t)(0,12)  = 0.600;  (*t)(2,12)  = -0.026;  (*t)(2,12)  =   0.026;
1371       (*t)(0,13)  = 0.650;  (*t)(2,13)  = -0.026;  (*t)(2,13)  =   0.026;
1372       (*t)(0,14)  = 0.700;  (*t)(2,14)  = -0.026;  (*t)(2,14)  =   0.026;
1373       (*t)(0,15)  = 0.750;  (*t)(2,15)  = -0.026;  (*t)(2,15)  =   0.026;
1374       (*t)(0,16)  = 0.800;  (*t)(2,16)  = -0.026;  (*t)(2,16)  =   0.026;
1375       (*t)(0,17)  = 0.850;  (*t)(2,17)  = -0.024;  (*t)(2,17)  =   0.024;
1376       (*t)(0,18)  = 0.900;  (*t)(2,18)  = -0.024;  (*t)(2,18)  =   0.024;
1377       (*t)(0,19)  = 0.950;  (*t)(2,19)  = -0.024;  (*t)(2,19)  =   0.024;
1378       (*t)(0,20)  = 1.000;  (*t)(2,20)  = -0.024;  (*t)(2,20)  =   0.024;
1379       (*t)(0,21)  = 1.100;  (*t)(2,21)  = -0.024;  (*t)(2,21)  =   0.024;
1380       (*t)(0,22)  = 1.200;  (*t)(2,22)  = -0.024;  (*t)(2,22)  =   0.022;
1381       (*t)(0,23)  = 1.300;  (*t)(2,23)  = -0.024;  (*t)(2,23)  =   0.020;
1382       (*t)(0,24)  = 1.400;  (*t)(2,24)  = -0.026;  (*t)(2,24)  =   0.016;
1383       (*t)(0,25)  = 1.500;  (*t)(2,25)  = -0.028;  (*t)(2,25)  =   0.014;
1384       (*t)(0,26)  = 1.600;  (*t)(2,26)  = -0.028;  (*t)(2,26)  =   0.012;
1385       (*t)(0,27)  = 1.700;  (*t)(2,27)  = -0.028;  (*t)(2,27)  =   0.010;
1386       (*t)(0,28)  = 1.800;  (*t)(2,28)  = -0.028;  (*t)(2,28)  =   0.010;
1387       (*t)(0,29)  = 1.900;  (*t)(2,29)  = -0.028;  (*t)(2,29)  =   0.008;
1388       (*t)(0,30)  = 2.000;  (*t)(2,30)  = -0.028;  (*t)(2,30)  =   0.006;
1389       (*t)(0,31)  = 2.100;  (*t)(2,31)  = -0.026;  (*t)(2,31)  =   0.006;
1390       (*t)(0,32)  = 2.200;  (*t)(2,32)  = -0.024;  (*t)(2,32)  =   0.004;
1391       (*t)(0,33)  = 2.300;  (*t)(2,33)  = -0.020;  (*t)(2,33)  =   0.002;
1392       (*t)(0,34)  = 2.400;  (*t)(2,34)  = -0.020;  (*t)(2,34)  =   0.002;
1393       (*t)(0,35)  = 2.500;  (*t)(2,35)  = -0.018;  (*t)(2,35)  =   0.000;
1394       (*t)(0,36)  = 2.600;  (*t)(2,36)  = -0.016;  (*t)(2,36)  =   0.000;
1395       (*t)(0,37)  = 2.700;  (*t)(2,37)  = -0.014;  (*t)(2,37)  =   -0.002;
1396       (*t)(0,38)  = 2.800;  (*t)(2,38)  = -0.014;  (*t)(2,38)  =   -0.004;
1397       (*t)(0,39)  = 2.900;  (*t)(2,39)  = -0.012;  (*t)(2,39)  =   -0.004;
1398       (*t)(0,40)  = 3.000;  (*t)(2,40)  = -0.010;  (*t)(2,40)  =   -0.006;
1399       (*t)(0,41)  = 3.100;  (*t)(2,41)  = 0.000;  (*t)(2,41)  =   0.000;
1400       (*t)(0,42)  = 3.200;  (*t)(2,42)  = 0.000;  (*t)(2,42)  =   0.000;
1401       (*t)(0,43)  = 3.300;  (*t)(2,43)  = 0.000;  (*t)(2,43)  =   0.000;
1402       (*t)(0,44)  = 3.400;  (*t)(2,44)  = 0.000;  (*t)(2,44)  =   0.000;
1403       (*t)(0,45)  = 3.500;  (*t)(2,45)  = 0.000;  (*t)(2,45)  =   0.000;
1404       (*t)(0,46)  = 3.600;  (*t)(2,46)  = 0.000;  (*t)(2,46)  =   0.000;
1405       (*t)(0,47)  = 3.700;  (*t)(2,47)  = 0.000;  (*t)(2,47)  =   0.000;
1406       (*t)(0,48)  = 3.800;  (*t)(2,48)  = 0.000;  (*t)(2,48)  =   0.000;
1407       (*t)(0,49)  = 3.900;  (*t)(2,49)  = 0.000;  (*t)(2,49)  =   0.000;
1408       (*t)(0,50)  = 4.000;  (*t)(2,50)  = 0.000;  (*t)(2,50)  =   0.000;
1409       (*t)(0,51)  = 4.100;  (*t)(2,51)  = 0.000;  (*t)(2,51)  =   0.000;
1410       (*t)(0,52)  = 4.200;  (*t)(2,52)  = 0.000;  (*t)(2,52)  =   0.000;
1411       (*t)(0,53)  = 4.300;  (*t)(2,53)  = 0.000;  (*t)(2,53)  =   0.000;
1412       (*t)(0,54)  = 4.400;  (*t)(2,54)  = 0.000;  (*t)(2,54)  =   0.000;
1413       (*t)(0,55)  = 4.500;  (*t)(2,55)  = 0.000;  (*t)(2,55)  =   0.000;
1414       (*t)(0,56)  = 4.600;  (*t)(2,56)  = 0.000;  (*t)(2,56)  =   0.000;
1415       (*t)(0,57)  = 4.700;  (*t)(2,57)  = 0.000;  (*t)(2,57)  =   0.000;
1416       (*t)(0,58)  = 4.800;  (*t)(2,58)  = 0.000;  (*t)(2,58)  =   0.000;
1417       (*t)(0,59)  = 4.900;  (*t)(2,59)  = 0.000;  (*t)(2,59)  =   0.000;
1418       (*t)(0,60)  = 5.900;  (*t)(2,60)  = 0.000;  (*t)(2,60)  =   0.000;
1419     }
1420     fTOFpidCuts=t;
1421   }
1422 }
1423
1424 //-----------------------------------------------------------------------
1425 // part added by F. Noferini (some methods)
1426 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(AliESDtrack* track){
1427
1428   Bool_t goodtrack = track && (track->GetStatus() & AliESDtrack::kTOFpid) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000) && (track->GetIntegratedLength() > 365) && !(track->GetStatus() & AliESDtrack::kTOFmismatch);
1429
1430   if (! goodtrack)
1431        return kFALSE;
1432
1433   Int_t pdg = GetESDPdg(track,"bayesianALL");
1434   //  printf("pdg set to %i\n",pdg);
1435
1436   Int_t pid = 0;
1437   Float_t prob = 0;
1438   switch (fParticleID)
1439   {
1440     case AliPID::kPion:
1441       pid=211;
1442       prob = fProbBayes[2];
1443       break;
1444     case AliPID::kKaon:
1445       pid=321;
1446       prob = fProbBayes[3];
1447      break;
1448     case AliPID::kProton:
1449       pid=2212;
1450       prob = fProbBayes[4];
1451       break;
1452     case AliPID::kElectron:
1453       pid=-11;
1454        prob = fProbBayes[0];
1455      break;
1456     default:
1457       return kFALSE;
1458   }
1459
1460   //  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);
1461   if(TMath::Abs(pdg) == TMath::Abs(pid) && prob > 0.8){
1462     if(!fCutCharge)
1463       return kTRUE;
1464     else if (fCutCharge && fCharge * track->GetSign() > 0)
1465       return kTRUE;
1466   }
1467   return kFALSE;
1468 }
1469 //-----------------------------------------------------------------------
1470 Int_t AliFlowTrackCuts::GetESDPdg(AliESDtrack *track,Option_t *option,Int_t ipart,Float_t cPi,Float_t cKa,Float_t cPr){
1471   Int_t pdg = 0;
1472   Int_t pdgvalues[5] = {-11,-13,211,321,2212};
1473   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
1474
1475   if(strstr(option,"bayesianTOF")){ // Bayesian TOF PID
1476     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1477     Double_t rcc=0.;
1478     
1479     Float_t pt = track->Pt();
1480     
1481     Int_t iptesd = 0;
1482     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1483   
1484     if(cPi < 0){
1485       c[0] = fC[iptesd][0];
1486       c[1] = fC[iptesd][1];
1487       c[2] = fC[iptesd][2];
1488       c[3] = fC[iptesd][3];
1489       c[4] = fC[iptesd][4];
1490     }
1491     else{
1492       c[0] = 0.0;
1493       c[1] = 0.0;
1494       c[2] = cPi;
1495       c[3] = cKa;
1496       c[4] = cPr;      
1497     }
1498
1499     Double_t r1[10]; track->GetTOFpid(r1);
1500     
1501     Int_t i;
1502     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1503     
1504     Double_t w[10];
1505     for (i=0; i<5; i++){
1506         w[i]=c[i]*r1[i]/rcc;
1507         fProbBayes[i] = w[i];
1508     }
1509     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1510       pdg = 211*Int_t(track->GetSign());
1511     }
1512     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1513       pdg = 2212*Int_t(track->GetSign());
1514     }
1515     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1516       pdg = 321*Int_t(track->GetSign());
1517     }
1518     else if (w[0]>=w[1]) { //electrons
1519       pdg = -11*Int_t(track->GetSign());
1520     }
1521     else{ // muon
1522       pdg = -13*Int_t(track->GetSign());
1523     }
1524   }
1525
1526   else if(strstr(option,"bayesianTPC")){ // Bayesian TPC PID
1527     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1528     Double_t rcc=0.;
1529     
1530     Float_t pt = track->Pt();
1531     
1532     Int_t iptesd = 0;
1533     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1534   
1535     if(cPi < 0){
1536       c[0] = fC[iptesd][0];
1537       c[1] = fC[iptesd][1];
1538       c[2] = fC[iptesd][2];
1539       c[3] = fC[iptesd][3];
1540       c[4] = fC[iptesd][4];
1541     }
1542     else{
1543       c[0] = 0.0;
1544       c[1] = 0.0;
1545       c[2] = cPi;
1546       c[3] = cKa;
1547       c[4] = cPr;      
1548     }
1549
1550     Double_t r1[10]; track->GetTPCpid(r1);
1551     
1552     Int_t i;
1553     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1554     
1555     Double_t w[10];
1556     for (i=0; i<5; i++){
1557         w[i]=c[i]*r1[i]/rcc;
1558         fProbBayes[i] = w[i];
1559     }
1560     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1561       pdg = 211*Int_t(track->GetSign());
1562     }
1563     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1564       pdg = 2212*Int_t(track->GetSign());
1565     }
1566     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1567       pdg = 321*Int_t(track->GetSign());
1568     }
1569     else if (w[0]>=w[1]) { //electrons
1570       pdg = -11*Int_t(track->GetSign());
1571     }
1572     else{ // muon
1573       pdg = -13*Int_t(track->GetSign());
1574     }
1575   }
1576   
1577   else if(strstr(option,"bayesianALL")){
1578     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1579     Double_t rcc=0.;
1580     
1581     Float_t pt = track->Pt();
1582     
1583     Int_t iptesd = 0;
1584     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1585
1586     if(cPi < 0){
1587       c[0] = fC[iptesd][0];
1588       c[1] = fC[iptesd][1];
1589       c[2] = fC[iptesd][2];
1590       c[3] = fC[iptesd][3];
1591       c[4] = fC[iptesd][4];
1592     }
1593     else{
1594       c[0] = 0.0;
1595       c[1] = 0.0;
1596       c[2] = cPi;
1597       c[3] = cKa;
1598       c[4] = cPr;      
1599     }
1600
1601     Double_t r1[10]; track->GetTOFpid(r1);
1602     Double_t r2[10]; track->GetTPCpid(r2);
1603
1604     Int_t i;
1605     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]*r2[i]);
1606     
1607
1608     Double_t w[10];
1609     for (i=0; i<5; i++){
1610         w[i]=c[i]*r1[i]*r2[i]/rcc;
1611         fProbBayes[i] = w[i];
1612     }
1613
1614     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1615       pdg = 211*Int_t(track->GetSign());
1616     }
1617     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1618       pdg = 2212*Int_t(track->GetSign());
1619     }
1620     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1621       pdg = 321*Int_t(track->GetSign());
1622     }
1623     else if (w[0]>=w[1]) { //electrons
1624       pdg = -11*Int_t(track->GetSign());
1625     }
1626     else{ // muon
1627       pdg = -13*Int_t(track->GetSign());
1628     }
1629   }
1630
1631   else if(strstr(option,"sigmacutTOF")){
1632     printf("PID not implemented yet: %s\nNO PID!!!!\n",option);
1633     Float_t p = track->P();
1634
1635     // Take expected times
1636     Double_t exptimes[5];
1637     track->GetIntegratedTimes(exptimes);
1638
1639     // Take resolution for TOF response
1640     // like fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1641     Float_t resolution = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1642
1643     if(TMath::Abs(exptimes[ipart] - track->GetTOFsignal()) < 3 * resolution){
1644       pdg = pdgvalues[ipart] * Int_t(track->GetSign());
1645     }
1646   }
1647
1648   else{
1649     printf("Invalid PID option: %s\nNO PID!!!!\n",option);
1650   }
1651
1652   return pdg;
1653 }
1654 //-----------------------------------------------------------------------
1655 void AliFlowTrackCuts::SetPriors(){
1656   // set abbundancies
1657     fBinLimitPID[0] = 0.30;
1658     fC[0][0] = 0.015;
1659     fC[0][1] = 0.015;
1660     fC[0][2] = 1;
1661     fC[0][3] = 0.0025;
1662     fC[0][4] = 0.000015;
1663     fBinLimitPID[1] = 0.35;
1664     fC[1][0] = 0.015;
1665     fC[1][1] = 0.015;
1666     fC[1][2] = 1;
1667     fC[1][3] = 0.01;
1668     fC[1][4] = 0.001;
1669     fBinLimitPID[2] = 0.40;
1670     fC[2][0] = 0.015;
1671     fC[2][1] = 0.015;
1672     fC[2][2] = 1;
1673     fC[2][3] = 0.026;
1674     fC[2][4] = 0.004;
1675     fBinLimitPID[3] = 0.45;
1676     fC[3][0] = 0.015;
1677     fC[3][1] = 0.015;
1678     fC[3][2] = 1;
1679     fC[3][3] = 0.026;
1680     fC[3][4] = 0.004;
1681     fBinLimitPID[4] = 0.50;
1682     fC[4][0] = 0.015;
1683     fC[4][1] = 0.015;
1684     fC[4][2] = 1.000000;
1685     fC[4][3] = 0.05;
1686     fC[4][4] = 0.01;
1687     fBinLimitPID[5] = 0.60;
1688     fC[5][0] = 0.012;
1689     fC[5][1] = 0.012;
1690     fC[5][2] = 1;
1691     fC[5][3] = 0.085;
1692     fC[5][4] = 0.022;
1693     fBinLimitPID[6] = 0.70;
1694     fC[6][0] = 0.01;
1695     fC[6][1] = 0.01;
1696     fC[6][2] = 1;
1697     fC[6][3] = 0.12;
1698     fC[6][4] = 0.036;
1699     fBinLimitPID[7] = 0.80;
1700     fC[7][0] = 0.0095;
1701     fC[7][1] = 0.0095;
1702     fC[7][2] = 1;
1703     fC[7][3] = 0.15;
1704     fC[7][4] = 0.05;
1705     fBinLimitPID[8] = 0.90;
1706     fC[8][0] = 0.0085;
1707     fC[8][1] = 0.0085;
1708     fC[8][2] = 1;
1709     fC[8][3] = 0.18;
1710     fC[8][4] = 0.074;
1711     fBinLimitPID[9] = 1;
1712     fC[9][0] = 0.008;
1713     fC[9][1] = 0.008;
1714     fC[9][2] = 1;
1715     fC[9][3] = 0.22;
1716     fC[9][4] = 0.1;
1717     fBinLimitPID[10] = 1.20;
1718     fC[10][0] = 0.007;
1719     fC[10][1] = 0.007;
1720     fC[10][2] = 1;
1721     fC[10][3] = 0.28;
1722     fC[10][4] = 0.16;
1723     fBinLimitPID[11] = 1.40;
1724     fC[11][0] = 0.0066;
1725     fC[11][1] = 0.0066;
1726     fC[11][2] = 1;
1727     fC[11][3] = 0.35;
1728     fC[11][4] = 0.23;
1729     fBinLimitPID[12] = 1.60;
1730     fC[12][0] = 0.0075;
1731     fC[12][1] = 0.0075;
1732     fC[12][2] = 1;
1733     fC[12][3] = 0.40;
1734     fC[12][4] = 0.31;
1735     fBinLimitPID[13] = 1.80;
1736     fC[13][0] = 0.0062;
1737     fC[13][1] = 0.0062;
1738     fC[13][2] = 1;
1739     fC[13][3] = 0.45;
1740     fC[13][4] = 0.39;
1741     fBinLimitPID[14] = 2.00;
1742     fC[14][0] = 0.005;
1743     fC[14][1] = 0.005;
1744     fC[14][2] = 1;
1745     fC[14][3] = 0.46;
1746     fC[14][4] = 0.47;
1747     fBinLimitPID[15] = 2.20;
1748     fC[15][0] = 0.0042;
1749     fC[15][1] = 0.0042;
1750     fC[15][2] = 1;
1751     fC[15][3] = 0.5;
1752     fC[15][4] = 0.55;
1753     fBinLimitPID[16] = 2.40;
1754     fC[16][0] = 0.007;
1755     fC[16][1] = 0.007;
1756     fC[16][2] = 1;
1757     fC[16][3] = 0.5;
1758     fC[16][4] = 0.6;
1759     
1760     for(Int_t i=17;i<fnPIDptBin;i++){
1761         fBinLimitPID[i] = 2.0 + 0.2 * (i-14);
1762         fC[i][0] = fC[13][0];
1763         fC[i][1] = fC[13][1];
1764         fC[i][2] = fC[13][2];
1765         fC[i][3] = fC[13][3];
1766         fC[i][4] = fC[13][4];
1767     }  
1768 }
1769 // end part added by F. Noferini
1770 //-----------------------------------------------------------------------
1771
1772
1773 //-----------------------------------------------------------------------
1774 const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
1775 {
1776   //get the name of the particle id source
1777   switch (s)
1778   {
1779     case kTPCdedx:
1780       return "TPCdedx";
1781     case kTOFbeta:
1782       return "TOFbeta";
1783     case kTPCpid:
1784       return "TPCpid";
1785     case kTOFpid:
1786       return "TOFpid";
1787     case kTOFbayesian:
1788       return "TOFbayesianPID";
1789     default:
1790       return "NOPID";
1791   }
1792 }
1793
1794 //-----------------------------------------------------------------------
1795 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type) 
1796 {
1797   //return the name of the selected parameter type
1798   switch (type)
1799   {
1800     case kMC:
1801       return "MC";
1802     case kGlobal:
1803       return "ESD global";
1804     case kESD_TPConly:
1805       return "TPC only";
1806     case kESD_SPDtracklet:
1807       return "SPD tracklet";
1808     case kPMD:
1809       return "PMD";
1810     default:
1811       return "unknown";
1812   }
1813 }
1814
1815 //-----------------------------------------------------------------------
1816 Bool_t AliFlowTrackCuts::PassesPMDcuts(AliESDPmdTrack* /*track*/ )
1817 {
1818   //check PMD specific cuts
1819   //clean up from last iteration, and init label
1820
1821   fTrack = NULL;
1822   fMCparticle=NULL;
1823   fTrackLabel=-1;
1824
1825   fTrackPhi = 0.;
1826   fTrackEta = 0.;
1827   fTrackWeight = 1.0;
1828
1829   Bool_t pass=kTRUE;
1830   if (fCutEta) {if (  fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
1831   if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
1832
1833   return kFALSE;
1834 }