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