some coverty fixes
[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   switch (fParamType)
660   {
661     case kGlobal:
662       fTrack = track;
663       break;
664     case kESD_TPConly:
665       if (!track->FillTPCOnlyTrack(fTPCtrack))
666       {
667         fTrack=NULL;
668         fMCparticle=NULL;
669         fTrackLabel=-1;
670         return;
671       }
672       fTrack = &fTPCtrack;
673       //recalculate the label and mc particle, they may differ as TPClabel != global label
674       fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
675       if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
676       else fMCparticle=NULL;
677       break;
678     default:
679       fTrack = track;
680       break;
681   }
682 }
683
684 //-----------------------------------------------------------------------
685 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts()
686 {
687   //get standard cuts
688   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly cuts");
689   delete cuts->fAliESDtrackCuts;
690   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
691   cuts->SetParamType(kESD_TPConly);
692   return cuts;
693 }
694
695 //-----------------------------------------------------------------------
696 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
697 {
698   //get standard cuts
699   AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
700   delete cuts->fAliESDtrackCuts;
701   cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
702   cuts->SetParamType(kGlobal);
703   return cuts;
704 }
705
706 //-----------------------------------------------------------------------
707 AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack(int index) const
708 {
709   //get a flow track constructed from whatever we applied cuts on
710   //caller is resposible for deletion
711   //if construction fails return NULL
712   AliFlowTrack* flowtrack=NULL;
713   TParticle *tmpTParticle=NULL;
714   AliMCParticle* tmpAliMCParticle=NULL;
715   if (fParamType==kESD_SPDtracklet)
716   {
717     switch (fParamMix)
718     {
719       case kPure:
720         flowtrack = new AliFlowTrack();
721         flowtrack->SetPhi(fTrackPhi);
722         flowtrack->SetEta(fTrackEta);
723         break;
724       case kTrackWithMCkine:
725         if (!fMCparticle) return NULL;
726         flowtrack = new AliFlowTrack();
727         flowtrack->SetPhi( fMCparticle->Phi() );
728         flowtrack->SetEta( fMCparticle->Eta() );
729         flowtrack->SetPt( fMCparticle->Pt() );
730         break;
731       case kTrackWithMCpt:
732         if (!fMCparticle) return NULL;
733         flowtrack = new AliFlowTrack();
734         flowtrack->SetPhi(fTrackPhi);
735         flowtrack->SetEta(fTrackEta);
736         flowtrack->SetPt(fMCparticle->Pt());
737         break;
738       case kTrackWithPtFromFirstMother:
739         if (!fMCparticle) return NULL;
740         flowtrack = new AliFlowTrack();
741         flowtrack->SetPhi(fTrackPhi);
742         flowtrack->SetEta(fTrackEta);
743         tmpTParticle = fMCparticle->Particle();
744         tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
745         flowtrack->SetPt(tmpAliMCParticle->Pt());
746         break;
747       default:
748         flowtrack = new AliFlowTrack();
749         flowtrack->SetPhi(fTrackPhi);
750         flowtrack->SetEta(fTrackEta);
751         break;
752     }
753     flowtrack->SetSource(AliFlowTrack::kFromTracklet);
754   }
755   else
756   {
757     if (!fTrack) return NULL;
758     switch(fParamMix)
759     {
760       case kPure:
761         flowtrack = new AliFlowTrack(fTrack);
762         break;
763       case kTrackWithMCkine:
764         flowtrack = new AliFlowTrack(fMCparticle);
765         break;
766       case kTrackWithMCPID:
767         flowtrack = new AliFlowTrack(fTrack);
768         //flowtrack->setPID(...) from mc, when implemented
769         break;
770       case kTrackWithMCpt:
771         if (!fMCparticle) return NULL;
772         flowtrack = new AliFlowTrack(fTrack);
773         flowtrack->SetPt(fMCparticle->Pt());
774         break;
775       case kTrackWithPtFromFirstMother:
776         if (!fMCparticle) return NULL;
777         flowtrack = new AliFlowTrack(fTrack);
778         tmpTParticle = fMCparticle->Particle();
779         tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
780         flowtrack->SetPt(tmpAliMCParticle->Pt());
781         break;
782       default:
783         flowtrack = new AliFlowTrack(fTrack);
784         break;
785     }
786     if (fParamType==kMC) flowtrack->SetSource(AliFlowTrack::kFromMC);
787     else if (dynamic_cast<AliESDtrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromESD);
788     else if (dynamic_cast<AliAODTrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromAOD);
789     else if (dynamic_cast<AliMCParticle*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromMC);
790   }
791
792   if(flowtrack)
793     flowtrack->SetIndexOnCollection(index);
794   return flowtrack;
795 }
796
797 //-----------------------------------------------------------------------
798 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
799 {
800   //check if current particle is a physical primary
801   if (!fMCevent) return kFALSE;
802   if (fTrackLabel<0) return kFALSE;
803   return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
804 }
805
806 //-----------------------------------------------------------------------
807 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
808 {
809   //check if current particle is a physical primary
810   Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
811   AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
812   if (!track) return kFALSE;
813   TParticle* particle = track->Particle();
814   Bool_t transported = particle->TestBit(kTransportBit);
815   //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
816         //(physprim && (transported || !requiretransported))?"YES":"NO"  );
817   return (physprim && (transported || !requiretransported));
818 }
819
820 //-----------------------------------------------------------------------
821 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type) 
822 {
823   //return the name of the selected parameter type
824   switch (type)
825   {
826     case kMC:
827       return "MC";
828     case kGlobal:
829       return "ESD global";
830     case kESD_TPConly:
831       return "TPC only";
832     case kESD_SPDtracklet:
833       return "SPD tracklet";
834     default:
835       return "unknown";
836   }
837 }
838
839 //-----------------------------------------------------------------------
840 void AliFlowTrackCuts::DefineHistograms()
841 {
842   //define qa histograms
843 }
844
845 //-----------------------------------------------------------------------
846 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
847 {
848   //get the number of tracks in the input event according source
849   //selection (ESD tracks, tracklets, MC particles etc.)
850   AliESDEvent* esd=NULL;
851   switch (fParamType)
852   {
853     case kESD_SPDtracklet:
854       esd = dynamic_cast<AliESDEvent*>(fEvent);
855       if (!esd) return 0;
856       return esd->GetMultiplicity()->GetNumberOfTracklets();
857     case kMC:
858       if (!fMCevent) return 0;
859       return fMCevent->GetNumberOfTracks();
860     default:
861       if (!fEvent) return 0;
862       return fEvent->GetNumberOfTracks();
863   }
864   return 0;
865 }
866
867 //-----------------------------------------------------------------------
868 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
869 {
870   //get the input object according the data source selection:
871   //(esd tracks, traclets, mc particles,etc...)
872   AliESDEvent* esd=NULL;
873   switch (fParamType)
874   {
875     case kESD_SPDtracklet:
876       esd = dynamic_cast<AliESDEvent*>(fEvent);
877       if (!esd) return NULL;
878       return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
879     case kMC:
880       if (!fMCevent) return NULL;
881       return fMCevent->GetTrack(i);
882     default:
883       if (!fEvent) return NULL;
884       return fEvent->GetTrack(i);
885   }
886 }
887
888 //-----------------------------------------------------------------------
889 void AliFlowTrackCuts::Clear(Option_t*)
890 {
891   //clean up
892   fTrack=NULL;
893   fMCevent=NULL;
894   fMCparticle=NULL;
895   fTrackLabel=0;
896   fTrackWeight=0.0;
897   fTrackEta=0.0;
898   fTrackPhi=0.0;
899 }
900
901 //-----------------------------------------------------------------------
902 Bool_t AliFlowTrackCuts::PassesTOFpidCut(AliESDtrack* t )
903 {
904   //check if passes PID cut using timing in TOF
905   Bool_t goodtrack = (t) && 
906                      (t->GetStatus() & AliESDtrack::kTOFpid) && 
907                      (t->GetTOFsignal() > 12000) && 
908                      (t->GetTOFsignal() < 100000) && 
909                      (t->GetIntegratedLength() > 365) && 
910                     !(t->GetStatus() & AliESDtrack::kTOFmismatch);
911
912   if (!goodtrack) return kFALSE;
913   
914   Float_t pt = t->Pt();
915   Float_t p = t->GetP();
916   Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
917   Float_t timeTOF = t->GetTOFsignal()- trackT0; 
918   //2=pion 3=kaon 4=protons
919   Double_t inttimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
920   t->GetIntegratedTimes(inttimes);
921   //construct the pid index because it's screwed up in TOF
922   Int_t pid = 0;
923   switch (fAliPID)
924   {
925     case AliPID::kPion:
926       pid=2;
927       break;
928     case AliPID::kKaon:
929       pid=3;
930       break;
931     case AliPID::kProton:
932       pid=4;
933       break;
934     default:
935       return kFALSE;
936   }
937   Float_t s = timeTOF-inttimes[pid];
938
939   Float_t* arr = fTOFpidCuts->GetMatrixArray();
940   Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
941   if (col<0) return kFALSE;
942   Float_t min = (*fTOFpidCuts)(1,col);
943   Float_t max = (*fTOFpidCuts)(2,col);
944
945   //printf("--------------TOF pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
946   return (s>min && s<max);
947 }
948
949 //-----------------------------------------------------------------------
950 Bool_t AliFlowTrackCuts::PassesTPCpidCut(AliESDtrack* track)
951 {
952   //check if passes PID cut using dedx signal in the TPC
953   if (!fTPCpidCuts)
954   {
955     printf("no TPCpidCuts\n");
956     return kFALSE;
957   }
958
959   const AliExternalTrackParam* tpcparam = track->GetInnerParam();
960   if (!tpcparam) return kFALSE;
961   Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(tpcparam->GetP(), fAliPID);
962   Float_t sigTPC = track->GetTPCsignal();
963   Float_t s = (sigTPC-sigExp)/sigExp;
964   Double_t pt = track->Pt();
965
966   Float_t* arr = fTPCpidCuts->GetMatrixArray();
967   Int_t col = TMath::BinarySearch(fTPCpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
968   if (col<0) return kFALSE;
969   Float_t min = (*fTPCpidCuts)(1,col);
970   Float_t max = (*fTPCpidCuts)(2,col);
971
972   //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
973   return (s>min && s<max);
974 }
975
976 //-----------------------------------------------------------------------
977 void AliFlowTrackCuts::InitPIDcuts()
978 {
979   //init matrices with PID cuts
980   TMatrixF* t = NULL;
981   if (!fTPCpidCuts)
982   {
983     if (fAliPID==AliPID::kPion)
984     {
985       t = new TMatrixF(3,10);
986       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.4;  (*t)(2,0)  =   0.2;
987       (*t)(0,1)  = 0.25;  (*t)(1,1)  = -0.4;  (*t)(2,1)  =   0.2;
988       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.4;  (*t)(2,2)  =  0.25;
989       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.4;  (*t)(2,3)  =  0.25;
990       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.4;  (*t)(2,4)  =   0.3;
991       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.4;  (*t)(2,5)  =   0.3;
992       (*t)(0,6)  = 0.50;  (*t)(1,6)  = -0.4;  (*t)(2,6)  =  0.25;
993       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.4;  (*t)(2,7)  =  0.15;
994       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.4;  (*t)(2,8)  =   0.1;
995       (*t)(0,9)  = 0.65;  (*t)(1,9)  =    0;  (*t)(2,9)  =     0;
996     }
997     else
998     if (fAliPID==AliPID::kKaon)
999     {
1000       t = new TMatrixF(3,7);
1001       (*t)(0,0)  = 0.20;  (*t)(1,0)  = -0.2;  (*t)(2,0)  = 0.4; 
1002       (*t)(0,1)  = 0.25;  (*t)(1,1)  =-0.15;  (*t)(2,1)  = 0.4;
1003       (*t)(0,2)  = 0.30;  (*t)(1,2)  = -0.1;  (*t)(2,2)  = 0.4;
1004       (*t)(0,3)  = 0.35;  (*t)(1,3)  = -0.1;  (*t)(2,3)  = 0.4;
1005       (*t)(0,4)  = 0.40;  (*t)(1,4)  = -0.1;  (*t)(2,4)  = 0.6;
1006       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.1;  (*t)(2,5)  = 0.6;
1007       (*t)(0,6)  = 0.50;  (*t)(1,6)  =    0;  (*t)(2,6)  =   0;
1008     }
1009     else
1010     if (fAliPID==AliPID::kProton)
1011     {
1012       t = new TMatrixF(3,16);
1013       (*t)(0,0)  = 0.20;  (*t)(1,0)  =     0;  (*t)(2,0)  =    0; 
1014       (*t)(0,1)  = 0.25;  (*t)(1,1)  =  -0.2;  (*t)(2,1)  =  0.3; 
1015       (*t)(0,2)  = 0.30;  (*t)(1,2)  =  -0.2;  (*t)(2,2)  =  0.6; 
1016       (*t)(0,3)  = 0.35;  (*t)(1,3)  =  -0.2;  (*t)(2,3)  =  0.6; 
1017       (*t)(0,4)  = 0.40;  (*t)(1,4)  =  -0.2;  (*t)(2,4)  =  0.6; 
1018       (*t)(0,5)  = 0.45;  (*t)(1,5)  = -0.15;  (*t)(2,5)  =  0.6; 
1019       (*t)(0,6)  = 0.50;  (*t)(1,6)  =  -0.1;  (*t)(2,6)  =  0.6; 
1020       (*t)(0,7)  = 0.55;  (*t)(1,7)  = -0.05;  (*t)(2,7)  =  0.6; 
1021       (*t)(0,8)  = 0.60;  (*t)(1,8)  = -0.05;  (*t)(2,8)  = 0.45; 
1022       (*t)(0,9)  = 0.65;  (*t)(1,9)  = -0.05;  (*t)(2,9)  = 0.45; 
1023       (*t)(0,10) = 0.70;  (*t)(1,10) = -0.05;  (*t)(2,10) = 0.45; 
1024       (*t)(0,11) = 0.75;  (*t)(1,11) = -0.05;  (*t)(2,11) = 0.45; 
1025       (*t)(0,12) = 0.80;  (*t)(1,12) =     0;  (*t)(2,12) = 0.45; 
1026       (*t)(0,13) = 0.85;  (*t)(1,13) =     0;  (*t)(2,13) = 0.45; 
1027       (*t)(0,14) = 0.90;  (*t)(1,14) =     0;  (*t)(2,14) = 0.45;
1028       (*t)(0,15) = 0.95;  (*t)(1,15) =     0;  (*t)(2,15) =    0;
1029     }
1030     fTPCpidCuts=t;
1031   }
1032   t = NULL;
1033   if (!fTOFpidCuts)
1034   {
1035     if (fAliPID==AliPID::kPion)
1036     {
1037       t = new TMatrixF(3,31);
1038       (*t)(0,0)  = 0.3;   (*t)(1,0)  = -700;  (*t)(2,0)  = 700;
1039       (*t)(0,1)  = 0.35;  (*t)(1,1)  = -800;  (*t)(2,1)  = 800;
1040       (*t)(0,2)  = 0.40;  (*t)(1,2)  = -600;  (*t)(2,2)  = 800;
1041       (*t)(0,3)  = 0.45;  (*t)(1,3)  = -500;  (*t)(2,3)  = 700;
1042       (*t)(0,4)  = 0.50;  (*t)(1,4)  = -400;  (*t)(2,4)  = 700;
1043       (*t)(0,5)  = 0.55;  (*t)(1,5)  = -400;  (*t)(2,5)  = 700;
1044       (*t)(0,6)  = 0.60;  (*t)(1,6)  = -400;  (*t)(2,6)  = 700;
1045       (*t)(0,7)  = 0.65;  (*t)(1,7)  = -400;  (*t)(2,7)  = 700;
1046       (*t)(0,8)  = 0.70;  (*t)(1,8)  = -400;  (*t)(2,8)  = 700;
1047       (*t)(0,9)  = 0.75;  (*t)(1,9)  = -400;  (*t)(2,9)  = 700;
1048       (*t)(0,10) = 0.80;  (*t)(1,10) = -400;  (*t)(2,10) = 600;
1049       (*t)(0,11) = 0.85;  (*t)(1,11) = -400;  (*t)(2,11) = 600;
1050       (*t)(0,12) = 0.90;  (*t)(1,12) = -400;  (*t)(2,12) = 600;
1051       (*t)(0,13) = 0.95;  (*t)(1,13) = -400;  (*t)(2,13) = 600;
1052       (*t)(0,14) = 1.00;  (*t)(1,14) = -400;  (*t)(2,14) = 550;
1053       (*t)(0,15) = 1.10;  (*t)(1,15) = -400;  (*t)(2,15) = 450;
1054       (*t)(0,16) = 1.20;  (*t)(1,16) = -400;  (*t)(2,16) = 400;
1055       (*t)(0,17) = 1.30;  (*t)(1,17) = -400;  (*t)(2,17) = 300;
1056       (*t)(0,18) = 1.40;  (*t)(1,18) = -400;  (*t)(2,18) = 300;
1057       (*t)(0,19) = 1.50;  (*t)(1,19) = -400;  (*t)(2,19) = 250;
1058       (*t)(0,20) = 1.60;  (*t)(1,20) = -400;  (*t)(2,20) = 200;
1059       (*t)(0,21) = 1.70;  (*t)(1,21) = -400;  (*t)(2,21) = 150;
1060       (*t)(0,22) = 1.80;  (*t)(1,22) = -400;  (*t)(2,22) = 100;
1061       (*t)(0,23) = 1.90;  (*t)(1,23) = -400;  (*t)(2,23) =  70;
1062       (*t)(0,24) = 2.00;  (*t)(1,24) = -400;  (*t)(2,24) =  50;
1063       (*t)(0,25) = 2.10;  (*t)(1,25) = -400;  (*t)(2,25) =   0;
1064       (*t)(0,26) = 2.20;  (*t)(1,26) = -400;  (*t)(2,26) =   0;
1065       (*t)(0,27) = 2.30;  (*t)(1,26) = -400;  (*t)(2,26) =   0;
1066       (*t)(0,28) = 2.40;  (*t)(1,26) = -400;  (*t)(2,26) = -50;
1067       (*t)(0,29) = 2.50;  (*t)(1,26) = -400;  (*t)(2,26) = -50;
1068       (*t)(0,30) = 2.60;  (*t)(1,26) =    0;  (*t)(2,26) =   0;
1069     }
1070     else
1071     if (fAliPID==AliPID::kProton)
1072     {
1073       t = new TMatrixF(3,39);
1074       (*t)(0,0)  = 0.3;  (*t)(1,0)   = 0;     (*t)(2,0) = 0;
1075       (*t)(0,1)  = 0.35;  (*t)(1,1)  = 0;     (*t)(2,1) = 0;
1076       (*t)(0,2)  = 0.40;  (*t)(1,2)  = 0;     (*t)(2,2) = 0;
1077       (*t)(0,3)  = 0.45;  (*t)(1,3)  = 0;     (*t)(2,3) = 0;
1078       (*t)(0,4)  = 0.50;  (*t)(1,4)  = 0;     (*t)(2,4) = 0;
1079       (*t)(0,5)  = 0.55;  (*t)(1,5)  = -900;  (*t)(2,5)  = 600;
1080       (*t)(0,6)  = 0.60;  (*t)(1,6)  = -800;  (*t)(2,6)  = 600;
1081       (*t)(0,7)  = 0.65;  (*t)(1,7)  = -800;  (*t)(2,7)  = 600;
1082       (*t)(0,8)  = 0.70;  (*t)(1,8)  = -800;  (*t)(2,8)  = 600;
1083       (*t)(0,9)  = 0.75;  (*t)(1,9)  = -700;  (*t)(2,9)  = 500;
1084       (*t)(0,10) = 0.80;  (*t)(1,10) = -700;  (*t)(2,10) = 500;
1085       (*t)(0,11) = 0.85;  (*t)(1,11) = -700;  (*t)(2,11) = 500;
1086       (*t)(0,12) = 0.90;  (*t)(1,12) = -600;  (*t)(2,12) = 500;
1087       (*t)(0,13) = 0.95;  (*t)(1,13) = -600;  (*t)(2,13) = 500;
1088       (*t)(0,14) = 1.00;  (*t)(1,14) = -600;  (*t)(2,14) = 500;
1089       (*t)(0,15) = 1.10;  (*t)(1,15) = -600;  (*t)(2,15) = 500;
1090       (*t)(0,16) = 1.20;  (*t)(1,16) = -500;  (*t)(2,16) = 500;
1091       (*t)(0,17) = 1.30;  (*t)(1,17) = -500;  (*t)(2,17) = 500;
1092       (*t)(0,18) = 1.40;  (*t)(1,18) = -500;  (*t)(2,18) = 500;
1093       (*t)(0,19) = 1.50;  (*t)(1,19) = -500;  (*t)(2,19) = 500;
1094       (*t)(0,20) = 1.60;  (*t)(1,20) = -400;  (*t)(2,20) = 500;
1095       (*t)(0,21) = 1.70;  (*t)(1,21) = -400;  (*t)(2,21) = 500;
1096       (*t)(0,22) = 1.80;  (*t)(1,22) = -400;  (*t)(2,22) = 500;
1097       (*t)(0,23) = 1.90;  (*t)(1,23) = -400;  (*t)(2,23) = 500;
1098       (*t)(0,24) = 2.00;  (*t)(1,24) = -400;  (*t)(2,24) = 500;
1099       (*t)(0,25) = 2.10;  (*t)(1,25) = -350;  (*t)(2,25) = 500;
1100       (*t)(0,26) = 2.20;  (*t)(1,26) = -350;  (*t)(2,26) = 500;
1101       (*t)(0,27) = 2.30;  (*t)(1,27) = -300;  (*t)(2,27) = 500;
1102       (*t)(0,28) = 2.40;  (*t)(1,28) = -300;  (*t)(2,28) = 500;
1103       (*t)(0,29) = 2.50;  (*t)(1,29) = -300;  (*t)(2,29) = 500;
1104       (*t)(0,30) = 2.60;  (*t)(1,30) = -250;  (*t)(2,30) = 500;
1105       (*t)(0,31) = 2.70;  (*t)(1,31) = -200;  (*t)(2,31) = 500;
1106       (*t)(0,32) = 2.80;  (*t)(1,32) = -150;  (*t)(2,32) = 500;
1107       (*t)(0,33) = 2.90;  (*t)(1,33) = -150;  (*t)(2,33) = 500;
1108       (*t)(0,34) = 3.00;  (*t)(1,34) = -100;  (*t)(2,34) = 400;
1109       (*t)(0,35) = 3.10;  (*t)(1,35) = -100;  (*t)(2,35) = 400;
1110       (*t)(0,36) = 3.50;  (*t)(1,36) = -100;  (*t)(2,36) = 400;
1111       (*t)(0,37) = 3.60;  (*t)(1,37) =    0;  (*t)(2,37) = 0;
1112       (*t)(0,38) = 3.70;  (*t)(1,38) =    0;  (*t)(2,38) = 0;
1113     }
1114     else
1115     if (fAliPID==AliPID::kKaon)
1116     {
1117       t = new TMatrixF(3,26);
1118       (*t)(0,0)  = 0.3;   (*t)(1,0)  =    0;  (*t)(2,0)  =    0;
1119       (*t)(0,1)  = 0.35;  (*t)(1,1)  =    0;  (*t)(2,1)  =    0;
1120       (*t)(0,2)  = 0.40;  (*t)(1,2)  = -800;  (*t)(2,2)  =  600;
1121       (*t)(0,3)  = 0.45;  (*t)(1,3)  = -800;  (*t)(2,3)  =  600;
1122       (*t)(0,4)  = 0.50;  (*t)(1,4)  = -800;  (*t)(2,4)  =  600;
1123       (*t)(0,5)  = 0.55;  (*t)(1,5)  = -800;  (*t)(2,5)  =  600;
1124       (*t)(0,6)  = 0.60;  (*t)(1,6)  = -800;  (*t)(2,6)  =  600;
1125       (*t)(0,7)  = 0.65;  (*t)(1,7)  = -700;  (*t)(2,7)  =  600;
1126       (*t)(0,8)  = 0.70;  (*t)(1,8)  = -600;  (*t)(2,8)  =  600;
1127       (*t)(0,9)  = 0.75;  (*t)(1,9)  = -600;  (*t)(2,9)  =  500;
1128       (*t)(0,10) = 0.80;  (*t)(1,10) = -500;  (*t)(2,10) =  500;
1129       (*t)(0,11) = 0.85;  (*t)(1,11) = -500;  (*t)(2,11) =  500;
1130       (*t)(0,12) = 0.90;  (*t)(1,12) = -400;  (*t)(2,12) =  500;
1131       (*t)(0,13) = 0.95;  (*t)(1,13) = -400;  (*t)(2,13) =  500;
1132       (*t)(0,14) = 1.00;  (*t)(1,14) = -400;  (*t)(2,14) =  500;
1133       (*t)(0,15) = 1.10;  (*t)(1,15) = -350;  (*t)(2,15) =  450;
1134       (*t)(0,16) = 1.20;  (*t)(1,16) = -300;  (*t)(2,16) =  400;
1135       (*t)(0,17) = 1.30;  (*t)(1,17) = -300;  (*t)(2,17) =  400;
1136       (*t)(0,18) = 1.40;  (*t)(1,18) = -250;  (*t)(2,18) =  400;
1137       (*t)(0,19) = 1.50;  (*t)(1,19) = -200;  (*t)(2,19) =  400;
1138       (*t)(0,20) = 1.60;  (*t)(1,20) = -150;  (*t)(2,20) =  400;
1139       (*t)(0,21) = 1.70;  (*t)(1,21) = -100;  (*t)(2,21) =  400;
1140       (*t)(0,22) = 1.80;  (*t)(1,22) =  -50;  (*t)(2,22) =  400;
1141       (*t)(0,23) = 1.90;  (*t)(1,22) =    0;  (*t)(2,22) =  400;
1142       (*t)(0,24) = 2.00;  (*t)(1,22) =   50;  (*t)(2,22) =  400;
1143       (*t)(0,25) = 2.10;  (*t)(1,22) =    0;  (*t)(2,22) =    0;
1144     }
1145     fTOFpidCuts=t;
1146   }
1147 }
1148
1149 //-----------------------------------------------------------------------
1150 // part added by F. Noferini (some methods)
1151 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(AliESDtrack* track){
1152
1153   Bool_t goodtrack = track && (track->GetStatus() & AliESDtrack::kTOFpid) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000) && (track->GetIntegratedLength() > 365) && !(track->GetStatus() & AliESDtrack::kTOFmismatch);
1154
1155   if (! goodtrack)
1156        return kFALSE;
1157
1158   Int_t pdg = GetESDPdg(track,"bayesianALL");
1159   //  printf("pdg set to %i\n",pdg);
1160
1161   Int_t pid = 0;
1162   Float_t prob = 0;
1163   switch (fAliPID)
1164   {
1165     case AliPID::kPion:
1166       pid=211;
1167       prob = fProbBayes[2];
1168       break;
1169     case AliPID::kKaon:
1170       pid=321;
1171       prob = fProbBayes[3];
1172      break;
1173     case AliPID::kProton:
1174       pid=2212;
1175       prob = fProbBayes[4];
1176       break;
1177     case AliPID::kElectron:
1178       pid=-11;
1179        prob = fProbBayes[0];
1180      break;
1181     default:
1182       return kFALSE;
1183   }
1184
1185   //  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);
1186   if(TMath::Abs(pdg) == TMath::Abs(pid) && prob > 0.8){
1187     if(!fCutCharge)
1188       return kTRUE;
1189     else if (fCutCharge && fCharge * track->GetSign() > 0)
1190       return kTRUE;
1191   }
1192   return kFALSE;
1193 }
1194 //-----------------------------------------------------------------------
1195 Int_t AliFlowTrackCuts::GetESDPdg(AliESDtrack *track,Option_t *option,Int_t ipart,Float_t cPi,Float_t cKa,Float_t cPr){
1196   Int_t pdg = 0;
1197   Int_t pdgvalues[5] = {-11,-13,211,321,2212};
1198   Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
1199
1200   if(strstr(option,"bayesianTOF")){ // Bayesian TOF PID
1201     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1202     Double_t rcc=0.;
1203     
1204     Float_t pt = track->Pt();
1205     
1206     Int_t iptesd = 0;
1207     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1208   
1209     if(cPi < 0){
1210       c[0] = fC[iptesd][0];
1211       c[1] = fC[iptesd][1];
1212       c[2] = fC[iptesd][2];
1213       c[3] = fC[iptesd][3];
1214       c[4] = fC[iptesd][4];
1215     }
1216     else{
1217       c[0] = 0.0;
1218       c[1] = 0.0;
1219       c[2] = cPi;
1220       c[3] = cKa;
1221       c[4] = cPr;      
1222     }
1223
1224     Double_t r1[10]; track->GetTOFpid(r1);
1225     
1226     Int_t i;
1227     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1228     
1229     Double_t w[10];
1230     for (i=0; i<5; i++){
1231         w[i]=c[i]*r1[i]/rcc;
1232         fProbBayes[i] = w[i];
1233     }
1234     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1235       pdg = 211*Int_t(track->GetSign());
1236     }
1237     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1238       pdg = 2212*Int_t(track->GetSign());
1239     }
1240     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1241       pdg = 321*Int_t(track->GetSign());
1242     }
1243     else if (w[0]>=w[1]) { //electrons
1244       pdg = -11*Int_t(track->GetSign());
1245     }
1246     else{ // muon
1247       pdg = -13*Int_t(track->GetSign());
1248     }
1249   }
1250
1251   else if(strstr(option,"bayesianTPC")){ // Bayesian TPC PID
1252     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1253     Double_t rcc=0.;
1254     
1255     Float_t pt = track->Pt();
1256     
1257     Int_t iptesd = 0;
1258     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1259   
1260     if(cPi < 0){
1261       c[0] = fC[iptesd][0];
1262       c[1] = fC[iptesd][1];
1263       c[2] = fC[iptesd][2];
1264       c[3] = fC[iptesd][3];
1265       c[4] = fC[iptesd][4];
1266     }
1267     else{
1268       c[0] = 0.0;
1269       c[1] = 0.0;
1270       c[2] = cPi;
1271       c[3] = cKa;
1272       c[4] = cPr;      
1273     }
1274
1275     Double_t r1[10]; track->GetTPCpid(r1);
1276     
1277     Int_t i;
1278     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1279     
1280     Double_t w[10];
1281     for (i=0; i<5; i++){
1282         w[i]=c[i]*r1[i]/rcc;
1283         fProbBayes[i] = w[i];
1284     }
1285     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1286       pdg = 211*Int_t(track->GetSign());
1287     }
1288     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1289       pdg = 2212*Int_t(track->GetSign());
1290     }
1291     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1292       pdg = 321*Int_t(track->GetSign());
1293     }
1294     else if (w[0]>=w[1]) { //electrons
1295       pdg = -11*Int_t(track->GetSign());
1296     }
1297     else{ // muon
1298       pdg = -13*Int_t(track->GetSign());
1299     }
1300   }
1301   
1302   else if(strstr(option,"bayesianALL")){
1303     Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1304     Double_t rcc=0.;
1305     
1306     Float_t pt = track->Pt();
1307     
1308     Int_t iptesd = 0;
1309     while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1310
1311     if(cPi < 0){
1312       c[0] = fC[iptesd][0];
1313       c[1] = fC[iptesd][1];
1314       c[2] = fC[iptesd][2];
1315       c[3] = fC[iptesd][3];
1316       c[4] = fC[iptesd][4];
1317     }
1318     else{
1319       c[0] = 0.0;
1320       c[1] = 0.0;
1321       c[2] = cPi;
1322       c[3] = cKa;
1323       c[4] = cPr;      
1324     }
1325
1326     Double_t r1[10]; track->GetTOFpid(r1);
1327     Double_t r2[10]; track->GetTPCpid(r2);
1328
1329     Int_t i;
1330     for (i=0; i<5; i++) rcc+=(c[i]*r1[i]*r2[i]);
1331     
1332
1333     Double_t w[10];
1334     for (i=0; i<5; i++){
1335         w[i]=c[i]*r1[i]*r2[i]/rcc;
1336         fProbBayes[i] = w[i];
1337     }
1338
1339     if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1340       pdg = 211*Int_t(track->GetSign());
1341     }
1342     else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1343       pdg = 2212*Int_t(track->GetSign());
1344     }
1345     else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1346       pdg = 321*Int_t(track->GetSign());
1347     }
1348     else if (w[0]>=w[1]) { //electrons
1349       pdg = -11*Int_t(track->GetSign());
1350     }
1351     else{ // muon
1352       pdg = -13*Int_t(track->GetSign());
1353     }
1354   }
1355
1356   else if(strstr(option,"sigmacutTOF")){
1357     printf("PID not implemented yet: %s\nNO PID!!!!\n",option);
1358     Float_t p = track->P();
1359
1360     // Take expected times
1361     Double_t exptimes[5];
1362     track->GetIntegratedTimes(exptimes);
1363
1364     // Take resolution for TOF response
1365     // like fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1366     Float_t resolution = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1367
1368     if(TMath::Abs(exptimes[ipart] - track->GetTOFsignal()) < 3 * resolution){
1369       pdg = pdgvalues[ipart] * Int_t(track->GetSign());
1370     }
1371   }
1372
1373   else{
1374     printf("Invalid PID option: %s\nNO PID!!!!\n",option);
1375   }
1376
1377   return pdg;
1378 }
1379 //-----------------------------------------------------------------------
1380 void AliFlowTrackCuts::SetPriors(){
1381   // set abbundancies
1382     fBinLimitPID[0] = 0.30;
1383     fC[0][0] = 0.015;
1384     fC[0][1] = 0.015;
1385     fC[0][2] = 1;
1386     fC[0][3] = 0.0025;
1387     fC[0][4] = 0.000015;
1388     fBinLimitPID[1] = 0.35;
1389     fC[1][0] = 0.015;
1390     fC[1][1] = 0.015;
1391     fC[1][2] = 1;
1392     fC[1][3] = 0.01;
1393     fC[1][4] = 0.001;
1394     fBinLimitPID[2] = 0.40;
1395     fC[2][0] = 0.015;
1396     fC[2][1] = 0.015;
1397     fC[2][2] = 1;
1398     fC[2][3] = 0.026;
1399     fC[2][4] = 0.004;
1400     fBinLimitPID[3] = 0.45;
1401     fC[3][0] = 0.015;
1402     fC[3][1] = 0.015;
1403     fC[3][2] = 1;
1404     fC[3][3] = 0.026;
1405     fC[3][4] = 0.004;
1406     fBinLimitPID[4] = 0.50;
1407     fC[4][0] = 0.015;
1408     fC[4][1] = 0.015;
1409     fC[4][2] = 1.000000;
1410     fC[4][3] = 0.05;
1411     fC[4][4] = 0.01;
1412     fBinLimitPID[5] = 0.60;
1413     fC[5][0] = 0.012;
1414     fC[5][1] = 0.012;
1415     fC[5][2] = 1;
1416     fC[5][3] = 0.085;
1417     fC[5][4] = 0.022;
1418     fBinLimitPID[6] = 0.70;
1419     fC[6][0] = 0.01;
1420     fC[6][1] = 0.01;
1421     fC[6][2] = 1;
1422     fC[6][3] = 0.12;
1423     fC[6][4] = 0.036;
1424     fBinLimitPID[7] = 0.80;
1425     fC[7][0] = 0.0095;
1426     fC[7][1] = 0.0095;
1427     fC[7][2] = 1;
1428     fC[7][3] = 0.15;
1429     fC[7][4] = 0.05;
1430     fBinLimitPID[8] = 0.90;
1431     fC[8][0] = 0.0085;
1432     fC[8][1] = 0.0085;
1433     fC[8][2] = 1;
1434     fC[8][3] = 0.18;
1435     fC[8][4] = 0.074;
1436     fBinLimitPID[9] = 1;
1437     fC[9][0] = 0.008;
1438     fC[9][1] = 0.008;
1439     fC[9][2] = 1;
1440     fC[9][3] = 0.22;
1441     fC[9][4] = 0.1;
1442     fBinLimitPID[10] = 1.20;
1443     fC[10][0] = 0.007;
1444     fC[10][1] = 0.007;
1445     fC[10][2] = 1;
1446     fC[10][3] = 0.28;
1447     fC[10][4] = 0.16;
1448     fBinLimitPID[11] = 1.40;
1449     fC[11][0] = 0.0066;
1450     fC[11][1] = 0.0066;
1451     fC[11][2] = 1;
1452     fC[11][3] = 0.35;
1453     fC[11][4] = 0.23;
1454     fBinLimitPID[12] = 1.60;
1455     fC[12][0] = 0.0075;
1456     fC[12][1] = 0.0075;
1457     fC[12][2] = 1;
1458     fC[12][3] = 0.40;
1459     fC[12][4] = 0.31;
1460     fBinLimitPID[13] = 1.80;
1461     fC[13][0] = 0.0062;
1462     fC[13][1] = 0.0062;
1463     fC[13][2] = 1;
1464     fC[13][3] = 0.45;
1465     fC[13][4] = 0.39;
1466     fBinLimitPID[14] = 2.00;
1467     fC[14][0] = 0.005;
1468     fC[14][1] = 0.005;
1469     fC[14][2] = 1;
1470     fC[14][3] = 0.46;
1471     fC[14][4] = 0.47;
1472     fBinLimitPID[15] = 2.20;
1473     fC[15][0] = 0.0042;
1474     fC[15][1] = 0.0042;
1475     fC[15][2] = 1;
1476     fC[15][3] = 0.5;
1477     fC[15][4] = 0.55;
1478     fBinLimitPID[16] = 2.40;
1479     fC[16][0] = 0.007;
1480     fC[16][1] = 0.007;
1481     fC[16][2] = 1;
1482     fC[16][3] = 0.5;
1483     fC[16][4] = 0.6;
1484     
1485     for(Int_t i=17;i<fnPIDptBin;i++){
1486         fBinLimitPID[i] = 2.0 + 0.2 * (i-14);
1487         fC[i][0] = fC[13][0];
1488         fC[i][1] = fC[13][1];
1489         fC[i][2] = fC[13][2];
1490         fC[i][3] = fC[13][3];
1491         fC[i][4] = fC[13][4];
1492     }  
1493 }
1494 // end part added by F. Noferini
1495 //-----------------------------------------------------------------------