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