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