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