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