]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliFlowTrackCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / 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. *
86a97121 14 **********************************************************i****************/
daf66719 15
16/* $Id$ */
17
18// AliFlowTrackCuts:
3c67c769 19// Data selection for flow framework
daf66719 20//
21// origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
86a97121 22// mods: Redmer A. Bertens (rbertens@cern.ch)
5559ce24 23//
24// This class gurantees consistency of cut methods, trackparameter
25// selection (global tracks, TPC only, etc..) and parameter mixing
26// in the flow framework. Transparently handles different input types:
27// ESD, MC, AOD.
28// This class works in 2 steps: first the requested track parameters are
29// constructed (to be set by SetParamType() ), then cuts are applied.
30// the constructed track can be requested AFTER checking the cuts by
31// calling GetTrack(), in this case the cut object stays in control,
32// caller does not have to delete the track.
33// Additionally caller can request an AliFlowTrack object to be constructed
34// according the parameter mixing scenario requested by SetParamMix().
35// AliFlowTrack is made using MakeFlowTrack() method, its an 'object factory'
36// so caller needs to take care of the freshly created object.
daf66719 37
67006c99 38#include "TFile.h"
daf66719 39#include <limits.h>
40#include <float.h>
32b846cd 41#include <TMatrix.h>
875ee12a 42#include "TMCProcess.h"
9a0783cc 43#include "TParticle.h"
1a80f9f6 44#include "TH2F.h"
9a0783cc 45#include "AliStack.h"
499fe731 46#include "TBrowser.h"
84b44ed4 47#include "TArrayD.h"
daf66719 48#include "AliMCEvent.h"
9a0783cc 49#include "AliESDEvent.h"
8eca5d19 50#include "AliAODEvent.h"
daf66719 51#include "AliVParticle.h"
86a97121 52#include "AliVVZERO.h"
daf66719 53#include "AliMCParticle.h"
dd294bb6
MK		 54#include "AliESDkink.h"
55#include "AliESDv0.h"
daf66719 56#include "AliESDtrack.h"
2e5052c5 57#include "AliESDMuonTrack.h" // XZhang 20120604
12b2b8bc 58#include "AliMultiplicity.h"
daf66719 59#include "AliAODTrack.h"
2e5052c5 60#include "AliAODTracklets.h" // XZhang 20120615
8eca5d19 61#include "AliFlowTrackSimple.h"
daf66719 62#include "AliFlowTrack.h"
63#include "AliFlowTrackCuts.h"
64#include "AliLog.h"
32b846cd 65#include "AliESDpid.h"
d148af7e 66#include "AliESDPmdTrack.h"
92b507bb 67#include "AliESDUtils.h" //TODO
cc0afcfc 68#include "AliFlowBayesianPID.h"
dd294bb6
MK		 69#include "AliFlowCandidateTrack.h"
70#include "AliKFParticle.h"
71#include "AliESDVZERO.h"
72#include "AliFlowCommonConstants.h"
043e9010 73#include "AliAnalysisManager.h"
74#include "AliPIDResponse.h"
daf66719 75
dc3e3e91 76using std::cout;
77using std::endl;
78
daf66719 79ClassImp(AliFlowTrackCuts)
80
81//-----------------------------------------------------------------------
82AliFlowTrackCuts::AliFlowTrackCuts():
441ea1cf 83 AliFlowTrackSimpleCuts(),
84 fAliESDtrackCuts(NULL),
2e5052c5 85 fMuonTrackCuts(NULL), // XZhang 20120604
441ea1cf 86 fQA(NULL),
87 fCutMC(kFALSE),
f21bdf48 88 fCutMChasTrackReferences(kFALSE),
441ea1cf 89 fCutMCprocessType(kFALSE),
90 fMCprocessType(kPNoProcess),
91 fCutMCPID(kFALSE),
92 fMCPID(0),
2ed97fda 93 fCutMCfirstMotherPID(kFALSE),
94 fMCfirstMotherPID(0),
a14b8f3c 95 fIgnoreSignInMCPID(kFALSE),
441ea1cf 96 fCutMCisPrimary(kFALSE),
97 fRequireTransportBitForPrimaries(kTRUE),
98 fMCisPrimary(kFALSE),
99 fRequireCharge(kFALSE),
100 fFakesAreOK(kTRUE),
101 fCutSPDtrackletDeltaPhi(kFALSE),
102 fSPDtrackletDeltaPhiMax(FLT_MAX),
103 fSPDtrackletDeltaPhiMin(-FLT_MAX),
104 fIgnoreTPCzRange(kFALSE),
105 fIgnoreTPCzRangeMax(FLT_MAX),
106 fIgnoreTPCzRangeMin(-FLT_MAX),
107 fCutChi2PerClusterTPC(kFALSE),
108 fMaxChi2PerClusterTPC(FLT_MAX),
109 fMinChi2PerClusterTPC(-FLT_MAX),
110 fCutNClustersTPC(kFALSE),
111 fNClustersTPCMax(INT_MAX),
112 fNClustersTPCMin(INT_MIN),
42655d16 113 fCutNClustersITS(kFALSE),
114 fNClustersITSMax(INT_MAX),
115 fNClustersITSMin(INT_MIN),
549b5f40
RAB		 116 fUseAODFilterBit(kTRUE),
117 fAODFilterBit(1),
a63303bd 118 fCutDCAToVertexXY(kFALSE),
119 fCutDCAToVertexZ(kFALSE),
2b1eaa10 120 fCutMinimalTPCdedx(kFALSE),
121 fMinimalTPCdedx(0.),
92b507bb 122 fLinearizeVZEROresponse(kFALSE),
67006c99 123 fCentralityPercentileMin(0.),
124 fCentralityPercentileMax(5.),
1a80f9f6 125 fCutPmdDet(kFALSE),
126 fPmdDet(0),
127 fCutPmdAdc(kFALSE),
128 fPmdAdc(0.),
129 fCutPmdNcell(kFALSE),
130 fPmdNcell(0.),
dd294bb6
MK		 131 fMinKinkAngle(TMath::DegToRad()*2.),
132 fMinKinkRadius(130.),
133 fMaxKinkRadius(200.),
134 fMinKinkQt(.05),
135 fMaxKinkQt(.5),
136 fMinKinkInvMassKmu(0.),
137 fMaxKinkInvMassKmu(0.6),
138 fForceTPCstandalone(kFALSE),
139 fRequireKinkDaughters(kFALSE),
441ea1cf 140 fParamType(kGlobal),
141 fParamMix(kPure),
dd294bb6
MK		 142 fKink(NULL),
143 fV0(NULL),
441ea1cf 144 fTrack(NULL),
dd294bb6
MK		 145 fTrackMass(0.),
146 fTrackPt(0.),
441ea1cf 147 fTrackPhi(0.),
148 fTrackEta(0.),
2281e7c5 149 fTrackWeight(1.),
441ea1cf 150 fTrackLabel(INT_MIN),
151 fMCevent(NULL),
152 fMCparticle(NULL),
153 fEvent(NULL),
154 fTPCtrack(),
aab6527a 155 fESDpid(),
9126fd0c 156 fBayesianResponse(NULL),
2b1eaa10 157 fPIDsource(kTOFpid),
441ea1cf 158 fTPCpidCuts(NULL),
159 fTOFpidCuts(NULL),
a14b8f3c 160 fParticleID(AliPID::kUnknown),
499fe731 161 fParticleProbability(.9),
875ee12a 162 fAllowTOFmismatchFlag(kFALSE),
0b14d0f4 163 fRequireStrictTOFTPCagreement(kFALSE),
9126fd0c 164 fCutRejectElectronsWithTPCpid(kFALSE),
67006c99 165 fUseTPCTOFNsigmaCutContours(kFALSE),
8b649d8c 166 fProbBayes(0.0),
86a97121 167 fCurrCentr(0.0),
c27333c2 168 fVZEROgainEqualization(NULL),
86a97121 169 fApplyRecentering(kFALSE),
043e9010 170 fVZEROgainEqualizationPerRing(kFALSE),
84b44ed4 171 fChi2A(0x0),
172 fChi2C(0x0),
173 fChi3A(0x0),
174 fChi3C(0x0),
043e9010 175 fPIDResponse(NULL),
67006c99 176 fNsigmaCut2(9),
177 fContoursFile(0),
178 fCutContourList(0),
179 fMaxITSclusterShared(0),
180 fMaxITSChi2(0)
441ea1cf 181{
182 //io constructor
422c61c7 183 SetPriors(); //init arrays
b147dd2f 184 // New PID procedure (Bayesian Combined PID)
185 // allocating here is necessary because we don't
186 // stream this member
187 // TODO: fix streaming problems AliFlowBayesianPID
188 fBayesianResponse = new AliFlowBayesianPID();
189 fBayesianResponse->SetNewTrackParam();
86a97121 190 for(Int_t i(0); i < 4; i++) {
c27333c2
MK		 191 fVZEROApol[i] = 0;
192 fVZEROCpol[i] = 0;
86a97121 193 }
c2a07d6a 194 for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = kTRUE;
67006c99 195
196 for(int i=0;i<50;i++){
197 fCutContour[i]= NULL;
198 fCutGraph[i]=NULL;
199 }
b147dd2f 200
441ea1cf 201}
202
203//-----------------------------------------------------------------------
204AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
cc0afcfc 205 AliFlowTrackSimpleCuts(name),
1a80f9f6 206 fAliESDtrackCuts(NULL),
2e5052c5 207 fMuonTrackCuts(NULL), // XZhang 20120604
a1c43d26 208 fQA(NULL),
1ff4bda1 209 fCutMC(kFALSE),
f21bdf48 210 fCutMChasTrackReferences(kFALSE),
daf66719 211 fCutMCprocessType(kFALSE),
212 fMCprocessType(kPNoProcess),
213 fCutMCPID(kFALSE),
214 fMCPID(0),
2ed97fda 215 fCutMCfirstMotherPID(kFALSE),
216 fMCfirstMotherPID(0),
a14b8f3c 217 fIgnoreSignInMCPID(kFALSE),
daf66719 218 fCutMCisPrimary(kFALSE),
441ea1cf 219 fRequireTransportBitForPrimaries(kTRUE),
daf66719 220 fMCisPrimary(kFALSE),
957517fa 221 fRequireCharge(kFALSE),
127a5825 222 fFakesAreOK(kTRUE),
9a0783cc 223 fCutSPDtrackletDeltaPhi(kFALSE),
224 fSPDtrackletDeltaPhiMax(FLT_MAX),
225 fSPDtrackletDeltaPhiMin(-FLT_MAX),
1ff4bda1 226 fIgnoreTPCzRange(kFALSE),
227 fIgnoreTPCzRangeMax(FLT_MAX),
228 fIgnoreTPCzRangeMin(-FLT_MAX),
2948ac5a 229 fCutChi2PerClusterTPC(kFALSE),
230 fMaxChi2PerClusterTPC(FLT_MAX),
231 fMinChi2PerClusterTPC(-FLT_MAX),
32b846cd 232 fCutNClustersTPC(kFALSE),
233 fNClustersTPCMax(INT_MAX),
234 fNClustersTPCMin(INT_MIN),
42655d16 235 fCutNClustersITS(kFALSE),
236 fNClustersITSMax(INT_MAX),
5ba02b32 237 fNClustersITSMin(INT_MIN),
549b5f40
RAB		 238 fUseAODFilterBit(kTRUE),
239 fAODFilterBit(1),
a63303bd 240 fCutDCAToVertexXY(kFALSE),
241 fCutDCAToVertexZ(kFALSE),
2b1eaa10 242 fCutMinimalTPCdedx(kFALSE),
243 fMinimalTPCdedx(0.),
92b507bb 244 fLinearizeVZEROresponse(kFALSE),
67006c99 245 fCentralityPercentileMin(0.),
246 fCentralityPercentileMax(5.),
1a80f9f6 247 fCutPmdDet(kFALSE),
248 fPmdDet(0),
249 fCutPmdAdc(kFALSE),
250 fPmdAdc(0.),
251 fCutPmdNcell(kFALSE),
252 fPmdNcell(0.),
dd294bb6
MK		 253 fMinKinkAngle(TMath::DegToRad()*2.),
254 fMinKinkRadius(130.),
255 fMaxKinkRadius(200.),
256 fMinKinkQt(0.05),
257 fMaxKinkQt(0.5),
258 fMinKinkInvMassKmu(0.0),
259 fMaxKinkInvMassKmu(0.6),
260 fForceTPCstandalone(kFALSE),
261 fRequireKinkDaughters(kFALSE),
12b2b8bc 262 fParamType(kGlobal),
daf66719 263 fParamMix(kPure),
dd294bb6
MK		 264 fKink(NULL),
265 fV0(NULL),
daf66719 266 fTrack(NULL),
dd294bb6
MK		 267 fTrackMass(0.),
268 fTrackPt(0.),
12b2b8bc 269 fTrackPhi(0.),
270 fTrackEta(0.),
2281e7c5 271 fTrackWeight(1.),
127a5825 272 fTrackLabel(INT_MIN),
957517fa 273 fMCevent(NULL),
9a0783cc 274 fMCparticle(NULL),
1ff4bda1 275 fEvent(NULL),
32b846cd 276 fTPCtrack(),
aab6527a 277 fESDpid(),
9126fd0c 278 fBayesianResponse(NULL),
2b1eaa10 279 fPIDsource(kTOFpid),
32b846cd 280 fTPCpidCuts(NULL),
281 fTOFpidCuts(NULL),
a14b8f3c 282 fParticleID(AliPID::kUnknown),
499fe731 283 fParticleProbability(.9),
875ee12a 284 fAllowTOFmismatchFlag(kFALSE),
0b14d0f4 285 fRequireStrictTOFTPCagreement(kFALSE),
9126fd0c 286 fCutRejectElectronsWithTPCpid(kFALSE),
67006c99 287 fUseTPCTOFNsigmaCutContours(kFALSE),
8b649d8c 288 fProbBayes(0.0),
86a97121 289 fCurrCentr(0.0),
c27333c2 290 fVZEROgainEqualization(NULL),
86a97121 291 fApplyRecentering(kFALSE),
043e9010 292 fVZEROgainEqualizationPerRing(kFALSE),
84b44ed4 293 fChi2A(0x0),
294 fChi2C(0x0),
295 fChi3A(0x0),
296 fChi3C(0x0),
043e9010 297 fPIDResponse(NULL),
67006c99 298 fNsigmaCut2(9),
299 fContoursFile(0),
300 fCutContourList(0),
301 fMaxITSclusterShared(0),
302 fMaxITSChi2(0)
daf66719 303{
67006c99 304 //constructor
441ea1cf 305 SetTitle("AliFlowTrackCuts");
aab6527a 306 fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
307 2.63394e+01,
308 5.04114e-11,
309 2.12543e+00,
310 4.88663e+00 );
422c61c7 311 SetPriors(); //init arrays
9126fd0c 312 // New PID procedure (Bayesian Combined PID)
62f2ea04 313 fBayesianResponse = new AliFlowBayesianPID();
314 fBayesianResponse->SetNewTrackParam();
86a97121 315 for(Int_t i(0); i < 4; i++) {
c27333c2
MK		 316 fVZEROApol[i] = 0;
317 fVZEROCpol[i] = 0;
86a97121 318 }
319 for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = kTRUE;
c2a07d6a 320
321 for(int i=0;i<50;i++){
322 fCutContour[i]= NULL;
323 fCutGraph[i]=NULL;
67006c99 324 }
c2a07d6a 325}
daf66719 326
327//-----------------------------------------------------------------------
ee242db3 328AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
329 AliFlowTrackSimpleCuts(that),
441ea1cf 330 fAliESDtrackCuts(NULL),
2e5052c5 331 fMuonTrackCuts(NULL), // XZhang 20120604
a1c43d26 332 fQA(NULL),
1ff4bda1 333 fCutMC(that.fCutMC),
f21bdf48 334 fCutMChasTrackReferences(that.fCutMChasTrackReferences),
ee242db3 335 fCutMCprocessType(that.fCutMCprocessType),
336 fMCprocessType(that.fMCprocessType),
337 fCutMCPID(that.fCutMCPID),
338 fMCPID(that.fMCPID),
2ed97fda 339 fCutMCfirstMotherPID(that.fCutMCfirstMotherPID),
340 fMCfirstMotherPID(that.fMCfirstMotherPID),
a14b8f3c 341 fIgnoreSignInMCPID(that.fIgnoreSignInMCPID),
ee242db3 342 fCutMCisPrimary(that.fCutMCisPrimary),
441ea1cf 343 fRequireTransportBitForPrimaries(that.fRequireTransportBitForPrimaries),
ee242db3 344 fMCisPrimary(that.fMCisPrimary),
345 fRequireCharge(that.fRequireCharge),
346 fFakesAreOK(that.fFakesAreOK),
347 fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
348 fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
349 fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
1ff4bda1 350 fIgnoreTPCzRange(that.fIgnoreTPCzRange),
351 fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
352 fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
2948ac5a 353 fCutChi2PerClusterTPC(that.fCutChi2PerClusterTPC),
354 fMaxChi2PerClusterTPC(that.fMaxChi2PerClusterTPC),
355 fMinChi2PerClusterTPC(that.fMinChi2PerClusterTPC),
32b846cd 356 fCutNClustersTPC(that.fCutNClustersTPC),
357 fNClustersTPCMax(that.fNClustersTPCMax),
358 fNClustersTPCMin(that.fNClustersTPCMin),
42655d16 359 fCutNClustersITS(that.fCutNClustersITS),
360 fNClustersITSMax(that.fNClustersITSMax),
361 fNClustersITSMin(that.fNClustersITSMin),
5ba02b32 362 fUseAODFilterBit(that.fUseAODFilterBit),
363 fAODFilterBit(that.fAODFilterBit),
a63303bd 364 fCutDCAToVertexXY(that.fCutDCAToVertexXY),
365 fCutDCAToVertexZ(that.fCutDCAToVertexZ),
2b1eaa10 366 fCutMinimalTPCdedx(that.fCutMinimalTPCdedx),
367 fMinimalTPCdedx(that.fMinimalTPCdedx),
92b507bb 368 fLinearizeVZEROresponse(that.fLinearizeVZEROresponse),
67006c99 369 fCentralityPercentileMin(that.fCentralityPercentileMin),
370 fCentralityPercentileMax(that.fCentralityPercentileMax),
1a80f9f6 371 fCutPmdDet(that.fCutPmdDet),
372 fPmdDet(that.fPmdDet),
373 fCutPmdAdc(that.fCutPmdAdc),
374 fPmdAdc(that.fPmdAdc),
375 fCutPmdNcell(that.fCutPmdNcell),
376 fPmdNcell(that.fPmdNcell),
dd294bb6
MK		 377 fMinKinkAngle(that.fMinKinkAngle),
378 fMinKinkRadius(that.fMinKinkRadius),
379 fMaxKinkRadius(that.fMaxKinkRadius),
380 fMinKinkQt(that.fMinKinkQt),
381 fMaxKinkQt(that.fMaxKinkQt),
382 fMinKinkInvMassKmu(that.fMinKinkInvMassKmu),
383 fMaxKinkInvMassKmu(that.fMaxKinkInvMassKmu),
384 fForceTPCstandalone(that.fForceTPCstandalone),
385 fRequireKinkDaughters(that.fRequireKinkDaughters),
ee242db3 386 fParamType(that.fParamType),
387 fParamMix(that.fParamMix),
dd294bb6
MK		 388 fKink(NULL),
389 fV0(NULL),
daf66719 390 fTrack(NULL),
dd294bb6
MK		 391 fTrackMass(0.),
392 fTrackPt(0.),
ee242db3 393 fTrackPhi(0.),
394 fTrackEta(0.),
2281e7c5 395 fTrackWeight(1.),
127a5825 396 fTrackLabel(INT_MIN),
957517fa 397 fMCevent(NULL),
9a0783cc 398 fMCparticle(NULL),
1ff4bda1 399 fEvent(NULL),
32b846cd 400 fTPCtrack(),
401 fESDpid(that.fESDpid),
9126fd0c 402 fBayesianResponse(NULL),
32b846cd 403 fPIDsource(that.fPIDsource),
404 fTPCpidCuts(NULL),
405 fTOFpidCuts(NULL),
2b1eaa10 406 fParticleID(that.fParticleID),
499fe731 407 fParticleProbability(that.fParticleProbability),
875ee12a 408 fAllowTOFmismatchFlag(that.fAllowTOFmismatchFlag),
0b14d0f4 409 fRequireStrictTOFTPCagreement(that.fRequireStrictTOFTPCagreement),
9126fd0c 410 fCutRejectElectronsWithTPCpid(that.fCutRejectElectronsWithTPCpid),
67006c99 411 fUseTPCTOFNsigmaCutContours(that.fUseTPCTOFNsigmaCutContours),
8b649d8c 412 fProbBayes(0.0),
86a97121 413 fCurrCentr(0.0),
c27333c2 414 fVZEROgainEqualization(NULL),
86a97121 415 fApplyRecentering(that.fApplyRecentering),
043e9010 416 fVZEROgainEqualizationPerRing(that.fVZEROgainEqualizationPerRing),
84b44ed4 417 fChi2A(0x0),
418 fChi2C(0x0),
419 fChi3A(0x0),
420 fChi3C(0x0),
043e9010 421 fPIDResponse(that.fPIDResponse),
67006c99 422 fNsigmaCut2(that.fNsigmaCut2),
423 fContoursFile(0),
424 fCutContourList(0),
425 fMaxITSclusterShared(0),
426 fMaxITSChi2(0)
daf66719 427{
428 //copy constructor
32b846cd 429 if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
430 if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
441ea1cf 431 if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
2e5052c5 432 if (that.fMuonTrackCuts) fMuonTrackCuts = new AliMuonTrackCuts(*(that.fMuonTrackCuts)); // XZhang 20120604
422c61c7 433 SetPriors(); //init arrays
67006c99 434 if(fUseTPCTOFNsigmaCutContours) GetTPCTOFPIDContours();
1a80f9f6 435 if (that.fQA) DefineHistograms();
84b44ed4 436 // would be neater via copy ctor of TArrayD
437 fChi2A = new TArrayD(that.fChi2A->GetSize(), that.fChi2A->GetArray());
438 fChi2C = new TArrayD(that.fChi2C->GetSize(), that.fChi2C->GetArray());
439 fChi3A = new TArrayD(that.fChi3A->GetSize(), that.fChi3A->GetArray());
440 fChi3C = new TArrayD(that.fChi3C->GetSize(), that.fChi3C->GetArray());
9126fd0c 441 // New PID procedure (Bayesian Combined PID)
62f2ea04 442 fBayesianResponse = new AliFlowBayesianPID();
443 fBayesianResponse->SetNewTrackParam();
67006c99 444
c27333c2
MK		 445 // VZERO gain calibration
446 // no reason to init fVZEROgainEqualizationPerRing, will be initialized on node if necessary
86a97121 447 // pointer is set to NULL in initialization list of this constructor
c27333c2 448// if (that.fVZEROgainEqualization) fVZEROgainEqualization = new TH1(*(that.fVZEROgainEqualization));
86a97121 449 for(Int_t i(0); i < 4; i++) { // no use to copy these guys since they're only initialized on worked node
c27333c2
MK		 450 fVZEROApol[i] = 0.;
451 fVZEROCpol[i] = 0.;
86a97121 452 }
453 for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = that.fUseVZERORing[i];
67006c99 454
455 for(int i=0;i<50;i++){
456 fCutContour[i]= that.fCutContour[i];
457 fCutGraph[i]=that.fCutGraph[i];
458 }
459
86a97121 460
daf66719 461}
462
463//-----------------------------------------------------------------------
ee242db3 464AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
daf66719 465{
466 //assignment
cdc20344 467 if (this==&that) return *this;
468
ee242db3 469 AliFlowTrackSimpleCuts::operator=(that);
1a80f9f6 470 //the following may seem excessive but if AliESDtrackCuts properly does copy and clone
471 //this approach is better memory-fragmentation-wise in some cases
472 if (that.fAliESDtrackCuts && fAliESDtrackCuts) *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
473 if (that.fAliESDtrackCuts && !fAliESDtrackCuts) fAliESDtrackCuts=new AliESDtrackCuts(*(that.fAliESDtrackCuts));
474 if (!that.fAliESDtrackCuts) delete fAliESDtrackCuts; fAliESDtrackCuts=NULL;
2e5052c5 475
476 if ( that.fMuonTrackCuts && fMuonTrackCuts) *fMuonTrackCuts = *(that.fMuonTrackCuts); // XZhang 20120604
52990716 477 if ( that.fMuonTrackCuts && !fMuonTrackCuts) fMuonTrackCuts = new AliMuonTrackCuts(*(that.fMuonTrackCuts)); // XZhang 20120604
2e5052c5 478 if (!that.fMuonTrackCuts) delete fMuonTrackCuts; fMuonTrackCuts = NULL; // XZhang 20120604
c27333c2 479 if (!that.fVZEROgainEqualization) delete fVZEROgainEqualization; fVZEROgainEqualization = NULL;
1a80f9f6 480 //these guys we don't need to copy, just reinit
481 if (that.fQA) {fQA->Delete(); delete fQA; fQA=NULL; DefineHistograms();}
1ff4bda1 482 fCutMC=that.fCutMC;
f21bdf48 483 fCutMChasTrackReferences=that.fCutMChasTrackReferences;
ee242db3 484 fCutMCprocessType=that.fCutMCprocessType;
485 fMCprocessType=that.fMCprocessType;
486 fCutMCPID=that.fCutMCPID;
487 fMCPID=that.fMCPID;
2ed97fda 488 fCutMCfirstMotherPID=that.fCutMCfirstMotherPID;
489 fMCfirstMotherPID=that.fMCfirstMotherPID;
a14b8f3c 490 fIgnoreSignInMCPID=that.fIgnoreSignInMCPID,
ee242db3 491 fCutMCisPrimary=that.fCutMCisPrimary;
441ea1cf 492 fRequireTransportBitForPrimaries=that.fRequireTransportBitForPrimaries;
ee242db3 493 fMCisPrimary=that.fMCisPrimary;
494 fRequireCharge=that.fRequireCharge;
495 fFakesAreOK=that.fFakesAreOK;
496 fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
497 fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
498 fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
1ff4bda1 499 fIgnoreTPCzRange=that.fIgnoreTPCzRange;
500 fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
501 fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
2948ac5a 502 fCutChi2PerClusterTPC=that.fCutChi2PerClusterTPC;
503 fMaxChi2PerClusterTPC=that.fMaxChi2PerClusterTPC;
504 fMinChi2PerClusterTPC=that.fMinChi2PerClusterTPC;
32b846cd 505 fCutNClustersTPC=that.fCutNClustersTPC;
506 fNClustersTPCMax=that.fNClustersTPCMax;
507 fNClustersTPCMin=that.fNClustersTPCMin;
42655d16 508 fCutNClustersITS=that.fCutNClustersITS;
509 fNClustersITSMax=that.fNClustersITSMax;
510 fNClustersITSMin=that.fNClustersITSMin;
a63303bd 511 fUseAODFilterBit=that.fUseAODFilterBit;
512 fAODFilterBit=that.fAODFilterBit;
513 fCutDCAToVertexXY=that.fCutDCAToVertexXY;
514 fCutDCAToVertexZ=that.fCutDCAToVertexZ;
2b1eaa10 515 fCutMinimalTPCdedx=that.fCutMinimalTPCdedx;
516 fMinimalTPCdedx=that.fMinimalTPCdedx;
92b507bb 517 fLinearizeVZEROresponse=that.fLinearizeVZEROresponse;
67006c99 518 fCentralityPercentileMin=that.fCentralityPercentileMin;
519 fCentralityPercentileMax=that.fCentralityPercentileMax;
1a80f9f6 520 fCutPmdDet=that.fCutPmdDet;
521 fPmdDet=that.fPmdDet;
522 fCutPmdAdc=that.fCutPmdAdc;
523 fPmdAdc=that.fPmdAdc;
524 fCutPmdNcell=that.fCutPmdNcell;
525 fPmdNcell=that.fPmdNcell;
dd294bb6
MK		 526 fMinKinkAngle=that.fMinKinkAngle;
527 fMinKinkRadius=that.fMinKinkRadius;
528 fMaxKinkRadius=that.fMaxKinkRadius;
529 fMinKinkQt=that.fMinKinkQt;
530 fMaxKinkQt=that.fMaxKinkQt;
531 fMinKinkInvMassKmu=that.fMinKinkInvMassKmu;
532 fMaxKinkInvMassKmu=that.fMaxKinkInvMassKmu;
533 fForceTPCstandalone=that.fForceTPCstandalone;
534 fRequireKinkDaughters=that.fRequireKinkDaughters;
1a80f9f6 535
ee242db3 536 fParamType=that.fParamType;
537 fParamMix=that.fParamMix;
daf66719 538
dd294bb6
MK		 539 fKink=NULL;
540 fV0=NULL;
daf66719 541 fTrack=NULL;
dd294bb6
MK		 542 fTrackMass=0.;
543 fTrackPt=0.;
1a80f9f6 544 fTrackEta=0.;
ee242db3 545 fTrackPhi=0.;
2281e7c5 546 fTrackWeight=1.;
127a5825 547 fTrackLabel=INT_MIN;
957517fa 548 fMCevent=NULL;
daf66719 549 fMCparticle=NULL;
9a0783cc 550 fEvent=NULL;
daf66719 551
32b846cd 552 fESDpid = that.fESDpid;
553 fPIDsource = that.fPIDsource;
554
cdc20344 555 delete fTPCpidCuts;
556 delete fTOFpidCuts;
32b846cd 557 if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
558 if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
32b846cd 559
2b1eaa10 560 fParticleID=that.fParticleID;
561 fParticleProbability=that.fParticleProbability;
875ee12a 562 fAllowTOFmismatchFlag=that.fAllowTOFmismatchFlag;
563 fRequireStrictTOFTPCagreement=that.fRequireStrictTOFTPCagreement;
0b14d0f4 564 fCutRejectElectronsWithTPCpid=that.fCutRejectElectronsWithTPCpid;
67006c99 565 fUseTPCTOFNsigmaCutContours=that.fUseTPCTOFNsigmaCutContours;
8b649d8c 566 fProbBayes = that.fProbBayes;
9126fd0c 567 fCurrCentr = that.fCurrCentr;
67006c99 568
86a97121 569 fApplyRecentering = that.fApplyRecentering;
c27333c2 570 fVZEROgainEqualizationPerRing = that.fVZEROgainEqualizationPerRing;
93ec5f2f 571#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
c27333c2 572 if (that.fVZEROgainEqualization) fVZEROgainEqualization = new TH1(*(that.fVZEROgainEqualization));
93ec5f2f 573#else
574 //PH Lets try Clone, however the result might be wrong
c27333c2 575 if (that.fVZEROgainEqualization) fVZEROgainEqualization = (TH1*)that.fVZEROgainEqualization->Clone();
93ec5f2f 576#endif
67006c99 577
86a97121 578 for(Int_t i(0); i < 4; i++) { // no use to copy these guys since they're only initialized on worked node
c27333c2
MK		 579 fVZEROApol[i] = that.fVZEROApol[i];
580 fVZEROCpol[i] = that.fVZEROCpol[i];
86a97121 581 }
582 for(Int_t i(0); i < 8; i++) fUseVZERORing[i] = that.fUseVZERORing[i];
84b44ed4 583 // would be neater via copy ctor of TArrayD
584 fChi2A = new TArrayD(that.fChi2A->GetSize(), that.fChi2A->GetArray());
585 fChi2C = new TArrayD(that.fChi2C->GetSize(), that.fChi2C->GetArray());
586 fChi3A = new TArrayD(that.fChi3A->GetSize(), that.fChi3A->GetArray());
587 fChi3C = new TArrayD(that.fChi3C->GetSize(), that.fChi3C->GetArray());
9126fd0c 588
043e9010 589 fPIDResponse = that.fPIDResponse;
590 fNsigmaCut2 = that.fNsigmaCut2;
591
daf66719 592 return *this;
593}
594
595//-----------------------------------------------------------------------
596AliFlowTrackCuts::~AliFlowTrackCuts()
597{
598 //dtor
daf66719 599 delete fAliESDtrackCuts;
32b846cd 600 delete fTPCpidCuts;
601 delete fTOFpidCuts;
2e5052c5 602 if (fMuonTrackCuts) delete fMuonTrackCuts; // XZhang 20120604
499fe731 603 if (fQA) { fQA->SetOwner(); fQA->Delete(); delete fQA; }
84b44ed4 604 if (fVZEROgainEqualization) {
605 delete fVZEROgainEqualization;
606 fVZEROgainEqualization = 0x0;
607 }
608 if(fChi2A) {
609 delete fChi2A;
610 fChi2A = 0x0;
611 }
612 if(fChi3A) {
613 delete fChi3A;
614 fChi3A = 0x0;
615 }
616 if(fChi2C) {
617 delete fChi2C;
618 fChi2C = 0x0;
619 }
620 if(fChi3C) {
621 delete fChi3C;
622 fChi3C = 0x0;
623 }
67006c99 624 fContoursFile->Close();
625 for(int i=0;i<50;i++){
626 delete fCutContour[i];
627 delete fCutGraph[i];
628 }
daf66719 629}
630
aab6527a 631//-----------------------------------------------------------------------
632void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
633{
634 //set the event
635 Clear();
636 fEvent=event;
637 fMCevent=mcEvent;
638
043e9010 639 // Get PID response
640 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
641 if(man){
642 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
643 if(inputHandler) fPIDResponse=inputHandler->GetPIDResponse();
644 }
645
aab6527a 646 //do the magic for ESD
647 AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(event);
0d16f553 648 AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(event);
aab6527a 649 if (fCutPID && myESD)
650 {
651 //TODO: maybe call it only for the TOF options?
652 // Added by F. Noferini for TOF PID
9126fd0c 653 // old procedure now implemented inside fBayesianResponse
654 // fESDpid.MakePID(myESD,kFALSE);
655 // new procedure
656 fBayesianResponse->SetDetResponse(myESD, fCurrCentr,AliESDpid::kTOF_T0,kTRUE); // centrality = PbPb centrality class (0-100%) or -1 for pp collisions
aab6527a 657 fESDpid.SetTOFResponse(myESD,AliESDpid::kTOF_T0);
aab6527a 658 // End F. Noferini added part
659 }
0d16f553 660 if (fCutPID && myAOD){
661 fBayesianResponse->SetDetResponse(myAOD, fCurrCentr,AliESDpid::kTOF_T0); // centrality = PbPb centrality class (0-100%) or -1 for pp collisions
662 if(myAOD->GetTOFHeader()){
663 fESDpid.SetTOFResponse(myAOD,AliESDpid::kTOF_T0);
664 }
665 else{ // corrected on the fly track by track if tof header is not present (old AODs)
666 }
667 // End F. Noferini added part
668 }
669
670 if(fPIDsource==kTOFbayesian) fBayesianResponse->SetDetAND(1);
671 else if(fPIDsource==kTPCbayesian) fBayesianResponse->ResetDetOR(1);
67006c99 672
673 if(fUseTPCTOFNsigmaCutContours) GetTPCTOFPIDContours();
aab6527a 674}
675
676//-----------------------------------------------------------------------
677void AliFlowTrackCuts::SetCutMC( Bool_t b )
678{
679 //will we be cutting on MC information?
680 fCutMC=b;
681
682 //if we cut on MC info then also the Bethe Bloch should be the one tuned for MC
683 if (fCutMC)
684 {
685 fESDpid.GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
686 1.75295e+01,
687 3.40030e-09,
688 1.96178e+00,
689 3.91720e+00);
690 }
691}
692
daf66719 693//-----------------------------------------------------------------------
12b2b8bc 694Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
daf66719 695{
cac040c5 696 //check cuts
daf66719 697 AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
2e5052c5 698//if (vparticle) return PassesCuts(vparticle); // XZhang 20120604
699 if (vparticle) { // XZhang 20120604
700 if (fParamType==kMUON) return PassesMuonCuts(vparticle); // XZhang 20120604
701 return PassesCuts(vparticle); // XZhang 20120604
702 } // XZhang 20120604
4aae2a93 703
daf66719 704 AliFlowTrackSimple* flowtrack = dynamic_cast<AliFlowTrackSimple*>(obj);
705 if (flowtrack) return PassesCuts(flowtrack);
12b2b8bc 706 AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
707 if (tracklets) return PassesCuts(tracklets,id);
2e5052c5 708 AliAODTracklets* trkletAOD = dynamic_cast<AliAODTracklets*>(obj); // XZhang 20120615
709 if (trkletAOD) return PassesCuts(trkletAOD,id); // XZhang 20120615
d148af7e 710 AliESDPmdTrack* pmdtrack = dynamic_cast<AliESDPmdTrack*>(obj);
711 if (pmdtrack) return PassesPMDcuts(pmdtrack);
3cca9bf4 712 AliVVZERO* vvzero = dynamic_cast<AliVVZERO*>(obj); // downcast to base class
c27333c2 713 if (vvzero) return PassesVZEROcuts(id);
dd294bb6
MK		 714 AliESDkink* kink = dynamic_cast<AliESDkink*>(obj);
715 if (kink) return PassesCuts(kink);
716 //AliESDv0* v0 = dynamic_cast<AliESDv0*>(obj);
717 //if (v0) return PassesCuts(v0);
3cca9bf4 718
daf66719 719 return kFALSE; //default when passed wrong type of object
720}
721
1ff4bda1 722//-----------------------------------------------------------------------
723Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
724{
725 //check cuts
726 AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
727 if (vparticle)
728 {
2281e7c5 729 return PassesMCcuts(fMCevent,(fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel());
1ff4bda1 730 }
731 AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
732 if (tracklets)
733 {
734 Int_t label0 = tracklets->GetLabel(id,0);
735 Int_t label1 = tracklets->GetLabel(id,1);
2281e7c5 736 Int_t label = (fFakesAreOK)?TMath::Abs(label0):label0;
737 if (label0!=label1) label = -666;
1ff4bda1 738 return PassesMCcuts(fMCevent,label);
739 }
740 return kFALSE; //default when passed wrong type of object
741}
742
daf66719 743//-----------------------------------------------------------------------
1a80f9f6 744Bool_t AliFlowTrackCuts::PassesCuts(const AliFlowTrackSimple* track)
daf66719 745{
746 //check cuts on a flowtracksimple
5559ce24 747
748 //clean up from last iteration
dd294bb6 749 ClearTrack();
daf66719 750 return AliFlowTrackSimpleCuts::PassesCuts(track);
751}
752
12b2b8bc 753//-----------------------------------------------------------------------
1a80f9f6 754Bool_t AliFlowTrackCuts::PassesCuts(const AliMultiplicity* tracklet, Int_t id)
12b2b8bc 755{
756 //check cuts on a tracklets
757
b35580e6 758 if (id<0) return kFALSE;
759
9a0783cc 760 //clean up from last iteration, and init label
dd294bb6 761 ClearTrack();
12b2b8bc 762
763 fTrackPhi = tracklet->GetPhi(id);
764 fTrackEta = tracklet->GetEta(id);
765 fTrackWeight = 1.0;
766 if (fCutEta) {if ( fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
767 if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
768
769 //check MC info if available
9a0783cc 770 //if the 2 clusters have different label track cannot be good
771 //and should therefore not pass the mc cuts
772 Int_t label0 = tracklet->GetLabel(id,0);
773 Int_t label1 = tracklet->GetLabel(id,1);
d0471ea0 774 //if possible get label and mcparticle
9a0783cc 775 fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
7afa829c 776 if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
d0471ea0 777 if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
778 //check MC cuts
1ff4bda1 779 if (fCutMC && !PassesMCcuts()) return kFALSE;
12b2b8bc 780 return kTRUE;
781}
782
2e5052c5 783//-----------------------------------------------------------------------
784Bool_t AliFlowTrackCuts::PassesCuts(const AliAODTracklets* tracklet, Int_t id)
785{
786 // XZhang 20120615
787 //check cuts on a tracklets
788
789 if (id<0) return kFALSE;
790
791 //clean up from last iteration, and init label
dd294bb6 792 ClearTrack();
2e5052c5 793
794 fTrackPhi = tracklet->GetPhi(id);
795//fTrackEta = tracklet->GetEta(id);
796 fTrackEta = -1.*TMath::Log(TMath::Tan(tracklet->GetTheta(id)/2.));
797 fTrackWeight = 1.0;
798 if (fCutEta) {if ( fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
799 if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
800
801 //check MC info if available
802 //if the 2 clusters have different label track cannot be good
803 //and should therefore not pass the mc cuts
804 Int_t label0 = tracklet->GetLabel(id,0);
805 Int_t label1 = tracklet->GetLabel(id,1);
806 //if possible get label and mcparticle
807 fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
808 if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
809 if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
810 //check MC cuts
811 if (fCutMC && !PassesMCcuts()) return kFALSE;
812 return kTRUE;
813}
814
12b2b8bc 815//-----------------------------------------------------------------------
1ff4bda1 816Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
12b2b8bc 817{
818 //check the MC info
1ff4bda1 819 if (!mcEvent) return kFALSE;
820 if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
821 AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
822 if (!mcparticle) {AliError("no MC track"); return kFALSE;}
12b2b8bc 823
824 if (fCutMCisPrimary)
825 {
441ea1cf 826 if (IsPhysicalPrimary(mcEvent,label,fRequireTransportBitForPrimaries) != fMCisPrimary) return kFALSE;
12b2b8bc 827 }
828 if (fCutMCPID)
829 {
1ff4bda1 830 Int_t pdgCode = mcparticle->PdgCode();
a14b8f3c 831 if (fIgnoreSignInMCPID)
4cbcbead 832 {
833 if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
834 }
835 else
836 {
837 if (fMCPID != pdgCode) return kFALSE;
838 }
12b2b8bc 839 }
2ed97fda 840 if (fCutMCfirstMotherPID)
841 {
842
843 TParticle* tparticle=mcparticle->Particle();
844 Int_t firstMotherLabel = 0;
845 if (tparticle) { firstMotherLabel = tparticle->GetFirstMother(); }
846 AliVParticle* firstMotherParticle = mcEvent->GetTrack(firstMotherLabel);
847 Int_t pdgcodeFirstMother = 0;
848 if (firstMotherParticle) { pdgcodeFirstMother = firstMotherParticle->PdgCode(); }
849 if (pdgcodeFirstMother != fMCfirstMotherPID) return kFALSE;
850 }
12b2b8bc 851 if ( fCutMCprocessType )
852 {
1ff4bda1 853 TParticle* particle = mcparticle->Particle();
12b2b8bc 854 Int_t processID = particle->GetUniqueID();
855 if (processID != fMCprocessType ) return kFALSE;
856 }
f21bdf48 857 if (fCutMChasTrackReferences)
858 {
859 if (mcparticle->GetNumberOfTrackReferences()<1) return kFALSE;
860 }
12b2b8bc 861 return kTRUE;
862}
1a80f9f6 863
1ff4bda1 864//-----------------------------------------------------------------------
865Bool_t AliFlowTrackCuts::PassesMCcuts()
866{
1a80f9f6 867 //check MC info
1ff4bda1 868 if (!fMCevent) return kFALSE;
869 if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
870 fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
871 return PassesMCcuts(fMCevent,fTrackLabel);
872}
12b2b8bc 873
dd294bb6
MK		 874//-----------------------------------------------------------------------
875Bool_t AliFlowTrackCuts::PassesCuts(const AliESDv0* v0)
876{
877 //check if the v0 passes cuts
878
879 //clean up from last iteration
880 ClearTrack();
881
882 fV0 = const_cast<AliESDv0*>(v0);
883
884 Bool_t pass = kTRUE;
885
886 //first, extract/prepare the v0
887 if (!v0->GetOnFlyStatus()) return kFALSE; //skip offline V0 finder
888 const AliExternalTrackParam *negHelix=v0->GetParamN();
889 const AliExternalTrackParam *posHelix=v0->GetParamP();
890 AliVParticle *v0tracks[2];
891 v0tracks[0] = fEvent->GetTrack(v0->GetNindex());
892 v0tracks[1] = fEvent->GetTrack(v0->GetPindex());
893 if( v0tracks[1]->Charge() < 0)
894 {
895 v0tracks[1] = fEvent->GetTrack(v0->GetNindex());
896 v0tracks[0] = fEvent->GetTrack(v0->GetPindex());
897 negHelix=v0->GetParamP();
898 posHelix=v0->GetParamN();
899 }
900
901 int KalmanPidPairs[4][2] =
902 {
903 {-11,11}, // e-,e+ (gamma)
904 {-211,2212}, // pi- p (lambda)
905 {-2212,211}, // anti-p pi+ (anti-lambda)
906 {-211,211} // pi- pi+ (Kshort)
907 // {-321,321} // K- K+ (phi)
908 };
909
910 Int_t id = 3;
911 //refit using a mass hypothesis
912 AliKFParticle v0trackKFneg(*(negHelix),KalmanPidPairs[id][0]);
913 AliKFParticle v0trackKFpos(*(posHelix),KalmanPidPairs[id][1]);
914 AliKFParticle v0particleRefit;
915 v0particleRefit += v0trackKFneg;
916 v0particleRefit += v0trackKFpos;
917 Double_t invMassErr= -999;
918 v0particleRefit.GetMass(fTrackMass,invMassErr);
919 //Double_t openAngle = v0trackKFneg.GetAngle(v0trackKFpos);
920 fTrackEta = v0particleRefit.GetEta();
921 fTrackPt = v0particleRefit.GetPt();
922 fTrackPhi = TMath::Pi()+v0particleRefit.GetPhi();
923 ////find out which mass bin and put the number in fPOItype
924 //Int_t massBins = AliFlowCommonConstants::GetMaster()->GetNbinsMass();
925 //Double_t massMin = AliFlowCommonConstants::GetMaster()->GetMassMin();
926 //Double_t massMax = AliFlowCommonConstants::GetMaster()->GetMassMax();
927 //fPOItype = 1 + int(massBins*(fTrackMass-massMin)/(massMax-massMin));
928
929 /////////////////////////////////////////////////////////////////////////////
930 //apply cuts
931 if ( v0tracks[0]->Charge() == v0tracks[1]->Charge() ) pass=kFALSE;
932 if ( v0tracks[0]->Pt()<0.15 || v0tracks[1]->Pt()<0.15 ) pass=kFALSE;
933
934 return pass;
935}
936
937//-----------------------------------------------------------------------
938Bool_t AliFlowTrackCuts::PassesCuts(const AliESDkink* kink)
939{
940 //check if the kink passes cuts
941
942 //clean up from last iteration
943 ClearTrack();
944 fKink=const_cast<AliESDkink*>(kink);
945
946 Bool_t pass = kTRUE;
947
948 Float_t kinkAngle = kink->GetAngle(2);
949 if (kinkAngle<fMinKinkAngle) pass = kFALSE;
950 Double_t kinkRadius = kink->GetR();
951 if (kinkRadius<fMinKinkRadius || kinkRadius>fMaxKinkRadius) pass = kFALSE;
952
953 //Qt
954 const TVector3 motherMfromKink(kink->GetMotherP());
955 const TVector3 daughterMfromKink(kink->GetDaughterP());
956 Float_t qt=kink->GetQt();
957 if ( qt < fMinKinkQt || qt > fMaxKinkQt) pass = kFALSE;
958
959 //invariant mass
960 Float_t energyDaughterMu = TMath::Sqrt( daughterMfromKink.Mag()*daughterMfromKink.Mag()+
961 0.105658*0.105658 );
962 Float_t p1XM = motherMfromKink.Px();
963 Float_t p1YM = motherMfromKink.Py();
964 Float_t p1ZM = motherMfromKink.Pz();
965 Float_t p2XM = daughterMfromKink.Px();
966 Float_t p2YM = daughterMfromKink.Py();
967 Float_t p2ZM = daughterMfromKink.Pz();
968 Float_t p3Daughter = TMath::Sqrt( ((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+
969 ((p1ZM-p2ZM)*(p1ZM-p2ZM)) );
970 Double_t invariantMassKmu = TMath::Sqrt( (energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-
971 motherMfromKink.Mag()*motherMfromKink.Mag() );
972 if (invariantMassKmu > fMaxKinkInvMassKmu) pass=kFALSE;
973 if (invariantMassKmu < fMinKinkInvMassKmu) pass=kFALSE;
974 fTrackMass=invariantMassKmu;
975
976 if (fQA)
977 {
978 QAbefore(13)->Fill(qt);
979 if (pass) QAafter(13)->Fill(qt);
980 QAbefore(14)->Fill(invariantMassKmu);
981 if (pass) QAafter(14)->Fill(invariantMassKmu);
982 const Double_t* kinkPosition = kink->GetPosition();
983 QAbefore(15)->Fill(kinkPosition[0],kinkPosition[1]);
984 if (pass) QAafter(15)->Fill(kinkPosition[0],kinkPosition[1]);
985 QAbefore(16)->Fill(motherMfromKink.Mag(),kinkAngle*TMath::RadToDeg());
986 if (pass) QAafter(16)->Fill(motherMfromKink.Mag(),kinkAngle*TMath::RadToDeg());
987 }
988
989 //mother track cuts
990 Int_t indexKinkMother = kink->GetIndex(0);
991 AliESDtrack* motherTrack = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(indexKinkMother));
992 if (!motherTrack) return kFALSE;
993 if (!PassesCuts(motherTrack)) pass = kFALSE;
994
995 return pass;
996}
997
daf66719 998//-----------------------------------------------------------------------
999Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
1000{
1001 //check cuts for an ESD vparticle
1002
127a5825 1003 ////////////////////////////////////////////////////////////////
1004 // start by preparing the track parameters to cut on //////////
1005 ////////////////////////////////////////////////////////////////
5559ce24 1006 //clean up from last iteration
dd294bb6
MK		 1007 ClearTrack();
1008 Bool_t pass=kTRUE;
5559ce24 1009
957517fa 1010 //get the label and the mc particle
127a5825 1011 fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
1012 if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
daf66719 1013 else fMCparticle=NULL;
1014
957517fa 1015 Bool_t isMCparticle = kFALSE; //some things are different for MC particles, check!
daf66719 1016 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*>(vparticle);
42655d16 1017 AliAODTrack* aodTrack = NULL;
daf66719 1018 if (esdTrack)
7d27a354 1019 {
42655d16 1020 //for an ESD track we do some magic sometimes like constructing TPC only parameters
1021 //or doing some hybrid, handle that here
daf66719 1022 HandleESDtrack(esdTrack);
7d27a354 1023 }
daf66719 1024 else
957517fa 1025 {
daf66719 1026 HandleVParticle(vparticle);
957517fa 1027 //now check if produced particle is MC
1028 isMCparticle = (dynamic_cast<AliMCParticle*>(fTrack))!=NULL;
42655d16 1029 aodTrack = dynamic_cast<AliAODTrack*>(vparticle); //keep the additional dynamic cast out of the way for ESDs
957517fa 1030 }
127a5825 1031 ////////////////////////////////////////////////////////////////
1032 ////////////////////////////////////////////////////////////////
1033
1ff4bda1 1034 if (!fTrack) return kFALSE;
42655d16 1035 //because it may be different from global, not needed for aodTrack because we dont do anything funky there
1036 if (esdTrack) esdTrack = static_cast<AliESDtrack*>(fTrack);
2948ac5a 1037
9a0783cc 1038 //check the common cuts for the current particle fTrack (MC,AOD,ESD)
32b846cd 1039 Double_t pt = fTrack->Pt();
f21bdf48 1040 Double_t p = fTrack->P();
7afa829c 1041 if (!fFakesAreOK) {if (fTrackLabel<0) pass=kFALSE;}
32b846cd 1042 if (fCutPt) {if (pt < fPtMin || pt >= fPtMax ) pass=kFALSE;}
924b02b0 1043 if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
1044 if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
1045 if (fRequireCharge) {if (fTrack->Charge() == 0) pass=kFALSE;}
1046 if (fCutCharge && !isMCparticle) {if (fTrack->Charge() != fCharge) pass=kFALSE;}
957517fa 1047 if (fCutCharge && isMCparticle)
1048 {
1049 //in case of an MC particle the charge is stored in units of 1/3|e|
1050 Int_t charge = TMath::Nint(fTrack->Charge()/3.0); //mc particles have charge in units of 1/3e
32b846cd 1051 if (charge!=fCharge) pass=kFALSE;
957517fa 1052 }
924b02b0 1053 //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
daf66719 1054
957517fa 1055 //when additionally MC info is required
1ff4bda1 1056 if (fCutMC && !PassesMCcuts()) pass=kFALSE;
daf66719 1057
42655d16 1058 //the case of ESD or AOD
bb6ca457 1059 if (esdTrack) { if (!PassesESDcuts(esdTrack)) { pass=kFALSE; } }
60875c3c 1060 if (aodTrack) { if (!PassesAODcuts(aodTrack,pass)) { pass=kFALSE; } }
42655d16 1061
f21bdf48 1062 if (fQA)
1063 {
c3eb0b6c 1064 if (fMCparticle)
f21bdf48 1065 {
875ee12a 1066 TParticle* tparticle=fMCparticle->Particle();
1067 Int_t processID = tparticle->GetUniqueID();
8eca5d19 1068 Int_t firstMotherLabel = tparticle->GetFirstMother();
3c67c769 1069 Bool_t primary = IsPhysicalPrimary();
875ee12a 1070 //TLorentzVector v;
1071 //mcparticle->Particle()->ProductionVertex(v);
1072 //Double_t prodvtxX = v.X();
1073 //Double_t prodvtxY = v.Y();
1074
c3eb0b6c 1075 Int_t pdgcode = fMCparticle->PdgCode();
483a17da 1076 AliVParticle* firstMotherParticle = fMCevent->GetTrack(firstMotherLabel);
1077 Int_t pdgcodeFirstMother = 0;
1078 if (firstMotherParticle) {pdgcodeFirstMother = firstMotherParticle->PdgCode();}
8eca5d19 1079
483a17da 1080 Float_t pdg = 0;
c3eb0b6c 1081 switch (TMath::Abs(pdgcode))
1082 {
1083 case 11:
1084 pdg = AliPID::kElectron + 0.5; break;
1085 case 13:
1086 pdg = AliPID::kMuon + 0.5; break;
1087 case 211:
1088 pdg = AliPID::kPion + 0.5; break;
1089 case 321:
1090 pdg = AliPID::kKaon + 0.5; break;
1091 case 2212:
1092 pdg = AliPID::kProton + 0.5; break;
1093 default:
1094 pdg = AliPID::kUnknown + 0.5; break;
1095 }
1096 pdg = TMath::Sign(pdg,static_cast<Float_t>(pdgcode));
8eca5d19 1097
2281e7c5 1098 Float_t geantCode = 0;
1099 switch (pdgcodeFirstMother)
8eca5d19 1100 {
2281e7c5 1101 case 22: //#gamma
1102 geantCode=1;
0284e4ef 1103 break;
2281e7c5 1104 case -11: //e^{+}
1105 geantCode=2;
0284e4ef 1106 break;
2281e7c5 1107 case 11: //e^{-}
1108 geantCode=3;
0284e4ef 1109 break;
2281e7c5 1110 case 12: case 14: case 16: //#nu
1111 geantCode=4;
1112 break;
1113 case -13:
1114 geantCode=5; //#mu^{+}
1115 break;
1116 case 13:
1117 geantCode=6; //#mu^{-}
1118 break;
1119 case 111: //#pi^{0}
1120 geantCode=7;
1121 break;
1122 case 211: //#pi^{+}
1123 geantCode=8;
1124 break;
1125 case -211: //#pi^{-}
1126 geantCode=9;
1127 break;
1128 case 130: //K^{0}_{L}
1129 geantCode=10;
1130 break;
1131 case 321: //K^{+}
1132 geantCode=11;
1133 break;
1134 case -321: //K^{-}
1135 geantCode=12;
1136 break;
1137 case 2112:
1138 geantCode = 13; //n
1139 break;
1140 case 2212:
1141 geantCode = 14; //p
1142 break;
1143 case -2212:
1144 geantCode=15; //#bar{p}
1145 break;
1146 case 310:
1147 geantCode=16; //K^{0}_{S}
1148 break;
1149 case 221:
1150 geantCode=17; //#eta
1151 break;
1152 case 3122:
1153 geantCode=18; //#Lambda
1154 break;
1155 case 3222:
1156 geantCode=19; //#Sigma^{+}
1157 break;
1158 case 3212:
1159 geantCode=20; //#Sigma^{0}
1160 break;
1161 case 3112:
1162 geantCode=21; //#Sigma^{-}
1163 break;
1164 case 3322:
1165 geantCode=22; //#Theta^{0}
1166 break;
1167 case 3312:
1168 geantCode=23; //#Theta^{-}
1169 break;
1170 case 3332:
1171 geantCode=24; //#Omega^{-}
1172 break;
1173 case -2112:
1174 geantCode=25; //#bar{n}
1175 break;
1176 case -3122:
1177 geantCode=26; //#bar{#Lambda}
1178 break;
1179 case -3112:
1180 geantCode=27; //#bar{#Sigma}^{-}
1181 break;
1182 case -3212:
1183 geantCode=28; //#bar{#Sigma}^{0}
1184 break;
1185 case -3222:
1186 geantCode=29; //#bar{#Sigma}^{+}
1187 break;
1188 case -3322:
1189 geantCode=30; //#bar{#Theta}^{0}
1190 break;
1191 case -3312:
1192 geantCode=31; //#bar{#Theta}^{+}
1193 break;
1194 case -3332:
1195 geantCode=32; //#bar{#Omega}^{+}
1196 break;
1197 case -15:
1198 geantCode=33; //#tau^{+}
1199 break;
1200 case 15:
1201 geantCode=34; //#tau^{-}
1202 break;
1203 case 411:
1204 geantCode=35; //D^{+}
1205 break;
1206 case -411:
1207 geantCode=36; //D^{-}
1208 break;
1209 case 421:
1210 geantCode=37; //D^{0}
1211 break;
1212 case -421:
1213 geantCode=38; //#bar{D}^{0}
1214 break;
1215 case 431:
1216 geantCode=39; //D_{s}^{+}
1217 break;
1218 case -431:
1219 geantCode=40; //#bar{D_{s}}^{-}
1220 break;
1221 case 4122:
1222 geantCode=41; //#Lambda_{c}^{+}
1223 break;
1224 case 24:
1225 geantCode=42; //W^{+}
1226 break;
1227 case -24:
1228 geantCode=43; //W^{-}
1229 break;
1230 case 23:
1231 geantCode=44; //Z^{0}
0284e4ef 1232 break;
8eca5d19 1233 default:
2281e7c5 1234 geantCode = 100;
0284e4ef 1235 break;
8eca5d19 1236 }
8eca5d19 1237
c3eb0b6c 1238 QAbefore(2)->Fill(p,pdg);
3c67c769 1239 QAbefore(3)->Fill(p,primary?0.5:-0.5);
875ee12a 1240 QAbefore(4)->Fill(p,static_cast<Float_t>(processID));
2281e7c5 1241 QAbefore(7)->Fill(p,geantCode+0.5);
c3eb0b6c 1242 if (pass) QAafter(2)->Fill(p,pdg);
3c67c769 1243 if (pass) QAafter(3)->Fill(p,primary?0.5:-0.5);
875ee12a 1244 if (pass) QAafter(4)->Fill(p,static_cast<Float_t>(processID));
2281e7c5 1245 if (pass) QAafter(7)->Fill(p,geantCode);
f21bdf48 1246 }
f21bdf48 1247 }
1248
42655d16 1249 //true by default, if we didn't set any cuts
1250 return pass;
1251}
1252
1253//_______________________________________________________________________
60875c3c 1254Bool_t AliFlowTrackCuts::PassesAODcuts(const AliAODTrack* track, Bool_t passedFid)
42655d16 1255{
1a80f9f6 1256 //check cuts for AOD
60875c3c 1257 Bool_t pass = passedFid;
42655d16 1258
1259 if (fCutNClustersTPC)
1ff4bda1 1260 {
42655d16 1261 Int_t ntpccls = track->GetTPCNcls();
1262 if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
1263 }
1264
1265 if (fCutNClustersITS)
1266 {
1267 Int_t nitscls = track->GetITSNcls();
1268 if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
1269 }
5ba02b32 1270
1271 if (fCutChi2PerClusterTPC)
1272 {
1273 Double_t chi2tpc = track->Chi2perNDF();
1274 if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
1275 }
1276
1277 if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
1278 if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
1279
1280 if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
1281
a63303bd 1282 if (fCutDCAToVertexXY && track->DCA()>GetMaxDCAToVertexXY()) pass=kFALSE;
42655d16 1283
8eca5d19 1284 if (fCutDCAToVertexZ && track->ZAtDCA()>GetMaxDCAToVertexZ()) pass=kFALSE;
1285
60875c3c 1286 Double_t dedx = track->GetTPCsignal();
1287 if (dedx < fMinimalTPCdedx) pass=kFALSE;
80efa607 1288 Double_t time[9];
60875c3c 1289 track->GetIntegratedTimes(time);
d0ce6147 1290 if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
1291 {
1292 if (!PassesAODpidCut(track)) pass=kFALSE;
1293 }
67006c99 1294/*
1295 if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid)) ) {
1296 Double_t c = TMath::C()*1.E-9;
1297 Double_t length = -999., beta =-999, tofTime = -999., tof = -999.;
1298 tofTime = track->GetTOFsignal();//in ps
1299 length = track->GetIntegratedLength();
1300
1301 tof = tofTime*1E-3; // ns
1302 if (tof > 0 && length > 0){
1303 length = length*0.01; // in meters
1304 tof = tof*c;
1305 beta = length/tof;
1306 if (fQA){
1307 QAbefore(18)->Fill(track->P()*track->Charge(),beta);
1308 if(pass) QAafter(18)->Fill(track->P()*track->Charge(),beta);
1309 }
1310 }
1311 }*/
1312
60875c3c 1313 if (fQA) {
d0ce6147 1314 // changed 04062014 used to be filled before possible PID cut
60875c3c 1315 Double_t momTPC = track->GetTPCmomentum();
a22c1eb1 1316 QAbefore( 0)->Fill(momTPC,GetBeta(track, kTRUE));
1317 if(pass) QAafter( 0)->Fill(momTPC, GetBeta(track, kTRUE));
60875c3c 1318 QAbefore( 1)->Fill(momTPC,dedx);
1319 QAbefore( 5)->Fill(track->Pt(),track->DCA());
1320 QAbefore( 6)->Fill(track->Pt(),track->ZAtDCA());
1321 if (pass) QAafter( 1)->Fill(momTPC,dedx);
1322 if (pass) QAafter( 5)->Fill(track->Pt(),track->DCA());
1323 if (pass) QAafter( 6)->Fill(track->Pt(),track->ZAtDCA());
1324 QAbefore( 8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1325 if (pass) QAafter( 8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1326 QAbefore( 9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1327 if (pass) QAafter( 9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1328 QAbefore(10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1329 if (pass) QAafter( 10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1330 QAbefore(11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1331 if (pass) QAafter( 11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1332 QAbefore(12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
1333 if (pass) QAafter( 12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
67006c99 1334
1335 //QAbefore(19)->Fill(track->P()*track->Charge(),track->GetDetPid()->GetTPCsignal());
1336 //if(pass) QAafter(19)->Fill(track->P()*track->Charge(),track->GetDetPid()->GetTPCsignal());
1337
60875c3c 1338 }
8eca5d19 1339
d0ce6147 1340
42655d16 1341 return pass;
1342}
1343
1344//_______________________________________________________________________
1345Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
1346{
1a80f9f6 1347 //check cuts on ESD tracks
42655d16 1348 Bool_t pass=kTRUE;
cd62a2a8 1349 Float_t dcaxy=0.0;
1350 Float_t dcaz=0.0;
1351 track->GetImpactParameters(dcaxy,dcaz);
499fe731 1352 const AliExternalTrackParam* pout = track->GetOuterParam();
1353 const AliExternalTrackParam* pin = track->GetInnerParam();
67006c99 1354
42655d16 1355 if (fIgnoreTPCzRange)
1356 {
42655d16 1357 if (pin&&pout)
441ea1cf 1358 {
42655d16 1359 Double_t zin = pin->GetZ();
1360 Double_t zout = pout->GetZ();
1361 if (zin*zout<0) pass=kFALSE; //reject if cross the membrane
1362 if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
1363 if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
441ea1cf 1364 }
42655d16 1365 }
32b846cd 1366
1a80f9f6 1367 Int_t ntpccls = ( fParamType==kTPCstandalone )?
42655d16 1368 track->GetTPCNclsIter1():track->GetTPCNcls();
1369 if (fCutChi2PerClusterTPC)
1370 {
1a80f9f6 1371 Float_t tpcchi2 = (fParamType==kTPCstandalone)?
42655d16 1372 track->GetTPCchi2Iter1():track->GetTPCchi2();
1373 tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
1374 if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
1375 pass=kFALSE;
1376 }
1377
2b1eaa10 1378 if (fCutMinimalTPCdedx)
1379 {
1380 if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
1381 }
1382
42655d16 1383 if (fCutNClustersTPC)
1384 {
1385 if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
1386 }
1387
1388 Int_t nitscls = track->GetNcls(0);
1389 if (fCutNClustersITS)
1390 {
1391 if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
1392 }
32b846cd 1393
1a80f9f6 1394 //some stuff is still handled by AliESDtrackCuts class - delegate
1395 if (fAliESDtrackCuts)
1396 {
1397 if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
1398 }
1399
cd62a2a8 1400 //PID part with pid QA
499fe731 1401 Double_t beta = GetBeta(track);
1402 Double_t dedx = Getdedx(track);
1403 if (fQA)
1404 {
1405 if (pass) QAbefore(0)->Fill(track->GetP(),beta);
1406 if (pass) QAbefore(1)->Fill(pin->GetP(),dedx);
1407 }
a14b8f3c 1408 if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
42655d16 1409 {
3c67c769 1410 if (!PassesESDpidCut(track)) pass=kFALSE;
1411 }
0b14d0f4 1412 if (fCutRejectElectronsWithTPCpid)
1413 {
1414 //reject electrons using standard TPC pid
3c67c769 1415 //TODO this should be rethought....
0b14d0f4 1416 Double_t pidTPC[AliPID::kSPECIES];
1417 track->GetTPCpid(pidTPC);
1418 if (pidTPC[AliPID::kElectron]<fParticleProbability) pass=kFALSE;
1419 }
67006c99 1420
1421 Int_t counterForSharedCluster=MaxSharedITSClusterCuts(fMaxITSclusterShared, track);
1422 if(counterForSharedCluster >= fMaxITSclusterShared) pass=kFALSE;
1423
1424 Double_t chi2perClusterITS = MaxChi2perITSClusterCuts(fMaxITSChi2, track);
1425 if(chi2perClusterITS >= fMaxITSChi2) pass=kFALSE;
1426
499fe731 1427 if (fQA)
1428 {
1429 if (pass) QAafter(0)->Fill(track->GetP(),beta);
1430 if (pass) QAafter(1)->Fill(pin->GetP(),dedx);
1431 }
cd62a2a8 1432 //end pid part with qa
3c67c769 1433
cd62a2a8 1434 if (fQA)
1435 {
1436 Double_t pt = track->Pt();
1437 QAbefore(5)->Fill(pt,dcaxy);
1438 QAbefore(6)->Fill(pt,dcaz);
1439 if (pass) QAafter(5)->Fill(pt,dcaxy);
1440 if (pass) QAafter(6)->Fill(pt,dcaz);
dd294bb6
MK		 1441 QAbefore(17)->Fill(Float_t(track->GetKinkIndex(0)));
1442 if (pass) QAafter(17)->Fill(Float_t(track->GetKinkIndex(0)));
cd62a2a8 1443 }
32b846cd 1444
42655d16 1445 return pass;
daf66719 1446}
1447
1448//-----------------------------------------------------------------------
1449void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
1450{
1451 //handle the general case
daf66719 1452 switch (fParamType)
1453 {
daf66719 1454 default:
daf66719 1455 fTrack = track;
32b846cd 1456 break;
daf66719 1457 }
1458}
1459
1460//-----------------------------------------------------------------------
1461void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
1462{
1463 //handle esd track
daf66719 1464 switch (fParamType)
1465 {
12b2b8bc 1466 case kGlobal:
daf66719 1467 fTrack = track;
daf66719 1468 break;
1a80f9f6 1469 case kTPCstandalone:
2b1eaa10 1470 if (!track->FillTPCOnlyTrack(fTPCtrack))
1ff4bda1 1471 {
1472 fTrack=NULL;
1473 fMCparticle=NULL;
2281e7c5 1474 fTrackLabel=-998;
1ff4bda1 1475 return;
1476 }
1477 fTrack = &fTPCtrack;
957517fa 1478 //recalculate the label and mc particle, they may differ as TPClabel != global label
127a5825 1479 fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1480 if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
957517fa 1481 else fMCparticle=NULL;
daf66719 1482 break;
daf66719 1483 default:
dd294bb6
MK		 1484 if (fForceTPCstandalone)
1485 {
1486 if (!track->FillTPCOnlyTrack(fTPCtrack))
1487 {
1488 fTrack=NULL;
1489 fMCparticle=NULL;
1490 fTrackLabel=-998;
1491 return;
1492 }
1493 fTrack = &fTPCtrack;
1494 //recalculate the label and mc particle, they may differ as TPClabel != global label
1495 fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1496 if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1497 else fMCparticle=NULL;
1498 }
1499 else
1500 fTrack = track;
32b846cd 1501 break;
daf66719 1502 }
1503}
1504
a0241c3a 1505//-----------------------------------------------------------------------
1506Int_t AliFlowTrackCuts::Count(AliVEvent* event)
1507{
1508 //calculate the number of track in given event.
1509 //if argument is NULL(default) take the event attached
1510 //by SetEvent()
67006c99 1511
a0241c3a 1512 Int_t multiplicity = 0;
1513 if (!event)
1514 {
1515 for (Int_t i=0; i<GetNumberOfInputObjects(); i++)
1516 {
1517 if (IsSelected(GetInputObject(i))) multiplicity++;
1518 }
1519 }
1520 else
1521 {
1522 for (Int_t i=0; i<event->GetNumberOfTracks(); i++)
1523 {
1524 if (IsSelected(event->GetTrack(i))) multiplicity++;
1525 }
1526 }
1527 return multiplicity;
1528}
1529
902ad6f1 1530//-----------------------------------------------------------------------
6c9a4feb 1531AliFlowTrackCuts* AliFlowTrackCuts::GetAODTrackCutsForFilterBit(UInt_t bit, TString suffix)
902ad6f1 1532{
1533 // object which in its default form only cuts on filterbit (for AOD analysis)
819dadc0 1534 // NOTE : the cut object is defined in such a way that ONLY filterbit is tested,
1535 // all other paramters are NOT checked
1536 // the exception is TPCdecx which is always cut for reasons of backward compatibility
1537 // but by setting the threshold to larger than -99999999 effectively there is no check
6c9a4feb 1538 AliFlowTrackCuts* cuts = new AliFlowTrackCuts(Form("AOD fitlerbit %i, %s", (int)bit, suffix.Data()));
902ad6f1 1539 cuts->SetMinimalTPCdedx(-999999999);
1540 cuts->SetAODfilterBit(bit);
7d7bf930 1541 cuts->SetParamType(AliFlowTrackCuts::kAODFilterBit);
902ad6f1 1542 return cuts;
1543}
1a80f9f6 1544//-----------------------------------------------------------------------
1545AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts()
54dd223b 1546{
1547 //returns the lhc10h vzero track cuts, this function
1548 //is left here for backward compatibility
1549 //if a run is recognized as 11h, the calibration method will
1550 //switch to 11h calbiration, which means that the cut
1551 //object is updated but not replaced.
1552 //calibratin is only available for PbPb runs
1553 return GetStandardVZEROOnlyTrackCuts2010();
1554}
1555//-----------------------------------------------------------------------
1556AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2010()
1a80f9f6 1557{
c27333c2 1558 //get standard VZERO cuts
07a3fb14 1559 //DISCLAIMER: LHC10h VZERO calibration consists (by default) of two steps
1560 //1) re-weigting of signal
1561 //2) re-centering of q-vectors
1562 //step 2 is available only for n==2 and n==3, for the higher harmonics the user
1563 //is repsonsible for making sure the q-sub distributions are (sufficiently) flat
1564 //or a sensible NUA procedure is applied !
6c9a4feb 1565 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts");
07a3fb14 1566 cuts->SetParamType(AliFlowTrackCuts::kVZERO);
1a80f9f6 1567 cuts->SetEtaRange( -10, +10 );
07a3fb14 1568 cuts->SetEtaGap(-1., 1.);
1a80f9f6 1569 cuts->SetPhiMin( 0 );
1570 cuts->SetPhiMax( TMath::TwoPi() );
86a97121 1571 // options for the reweighting
c27333c2 1572 cuts->SetVZEROgainEqualizationPerRing(kFALSE);
54dd223b 1573 cuts->SetApplyRecentering(kTRUE);
07a3fb14 1574 // to exclude a ring , do e.g.
54dd223b 1575 // cuts->SetUseVZERORing(7, kFALSE);
07a3fb14 1576 // excluding a ring will break the re-centering as re-centering relies on a
1577 // database file which tuned to receiving info from all rings
54dd223b 1578 return cuts;
1579}
1580//-----------------------------------------------------------------------
1581AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2011()
1582{
c27333c2 1583 //get standard VZERO cuts for 2011 data
54dd223b 1584 //in this case, the vzero segments will be weighted by
1585 //VZEROEqMultiplicity,
1586 //if recentering is enableded, the sub-q vectors
1587 //will be taken from the event header, so make sure to run
1588 //the VZERO event plane selection task before this task !
07a3fb14 1589 //DISCLAIMER: recentering is only available for n==2
1590 //for the higher harmonics the user
1591 //is repsonsible for making sure the q-sub distributions are (sufficiently) flat
1592 //or a sensible NUA procedure is applied !
54dd223b 1593 //recentering replaces the already evaluated q-vectors, so
1594 //when chosen, additional settings (e.g. excluding rings)
1595 //have no effect. recentering is true by default
1596 //
1597 //NOTE user is responsible for running the vzero event plane
1598 //selection task in advance, e.g. add to your launcher macro
1599 //
1600 // gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C");
1601 // AddTaskVZEROEPSelection();
1602 //
1603 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts 2011");
c27333c2 1604 cuts->SetParamType(kVZERO);
54dd223b 1605 cuts->SetEtaRange( -10, +10 );
07a3fb14 1606 cuts->SetEtaGap(-1., 1.);
54dd223b 1607 cuts->SetPhiMin( 0 );
1608 cuts->SetPhiMax( TMath::TwoPi() );
1609 cuts->SetApplyRecentering(kTRUE);
d0ce6147 1610 cuts->SetVZEROgainEqualizationPerRing(kFALSE);
1611 return cuts;
1a80f9f6 1612}
a0241c3a 1613//-----------------------------------------------------------------------
1614AliFlowTrackCuts* AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()
1615{
1616 //get standard cuts
1a80f9f6 1617 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard Global tracks");
a0241c3a 1618 cuts->SetParamType(kGlobal);
1619 cuts->SetPtRange(0.2,5.);
1620 cuts->SetEtaRange(-0.8,0.8);
1621 cuts->SetMinNClustersTPC(70);
1622 cuts->SetMinChi2PerClusterTPC(0.1);
1623 cuts->SetMaxChi2PerClusterTPC(4.0);
1624 cuts->SetMinNClustersITS(2);
1625 cuts->SetRequireITSRefit(kTRUE);
1626 cuts->SetRequireTPCRefit(kTRUE);
1627 cuts->SetMaxDCAToVertexXY(0.3);
1628 cuts->SetMaxDCAToVertexZ(0.3);
1629 cuts->SetAcceptKinkDaughters(kFALSE);
2b1eaa10 1630 cuts->SetMinimalTPCdedx(10.);
a0241c3a 1631
1632 return cuts;
1633}
1634
1635//-----------------------------------------------------------------------
1a80f9f6 1636AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010()
a0241c3a 1637{
1638 //get standard cuts
1a80f9f6 1639 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone 2010");
1640 cuts->SetParamType(kTPCstandalone);
a0241c3a 1641 cuts->SetPtRange(0.2,5.);
1642 cuts->SetEtaRange(-0.8,0.8);
1643 cuts->SetMinNClustersTPC(70);
1644 cuts->SetMinChi2PerClusterTPC(0.2);
1645 cuts->SetMaxChi2PerClusterTPC(4.0);
1646 cuts->SetMaxDCAToVertexXY(3.0);
1647 cuts->SetMaxDCAToVertexZ(3.0);
1648 cuts->SetDCAToVertex2D(kTRUE);
1649 cuts->SetAcceptKinkDaughters(kFALSE);
2b1eaa10 1650 cuts->SetMinimalTPCdedx(10.);
a0241c3a 1651
1652 return cuts;
1653}
1654
daf66719 1655//-----------------------------------------------------------------------
1a80f9f6 1656AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts()
daf66719 1657{
1658 //get standard cuts
1a80f9f6 1659 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone");
1660 cuts->SetParamType(kTPCstandalone);
a0241c3a 1661 cuts->SetPtRange(0.2,5.);
1662 cuts->SetEtaRange(-0.8,0.8);
1663 cuts->SetMinNClustersTPC(70);
1664 cuts->SetMinChi2PerClusterTPC(0.2);
1665 cuts->SetMaxChi2PerClusterTPC(4.0);
1666 cuts->SetMaxDCAToVertexXY(3.0);
1667 cuts->SetMaxDCAToVertexZ(3.0);
1668 cuts->SetDCAToVertex2D(kTRUE);
1669 cuts->SetAcceptKinkDaughters(kFALSE);
2b1eaa10 1670 cuts->SetMinimalTPCdedx(10.);
a0241c3a 1671
daf66719 1672 return cuts;
1673}
1674
1675//-----------------------------------------------------------------------
1676AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
1677{
1678 //get standard cuts
441ea1cf 1679 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
daf66719 1680 delete cuts->fAliESDtrackCuts;
1681 cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
12b2b8bc 1682 cuts->SetParamType(kGlobal);
daf66719 1683 return cuts;
1684}
1685
2e5052c5 1686//-----------------------------------------------------------------------------
84714904 1687AliFlowTrackCuts* AliFlowTrackCuts::GetStandardMuonTrackCuts(Bool_t isMC, Int_t passN)
2e5052c5 1688{
1689// XZhang 20120604
1690 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard muon track cuts");
1691 cuts->SetParamType(kMUON);
1692 cuts->SetStandardMuonTrackCuts();
1693 cuts->SetIsMuonMC(isMC);
84714904 1694 cuts->SetMuonPassNumber(passN);
2e5052c5 1695 return cuts;
1696}
1697
dd294bb6
MK		 1698//-----------------------------------------------------------------------
1699//AliFlowTrack* AliFlowTrackCuts::FillFlowTrackV0(TObjArray* trackCollection, Int_t trackIndex) const
1700//{
1701// //fill flow track from a reconstructed V0 (topological)
1702// AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1703// if (flowtrack)
1704// {
1705// flowtrack->Clear();
1706// }
1707// else
1708// {
1709// flowtrack = new AliFlowCandidateTrack();
1710// trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1711// }
1712//
1713// TParticle *tmpTParticle=NULL;
1714// AliMCParticle* tmpAliMCParticle=NULL;
1715// AliExternalTrackParam* externalParams=NULL;
1716// AliESDtrack* esdtrack=NULL;
1717// switch(fParamMix)
1718// {
1719// case kPure:
1720// flowtrack->Set(fTrack);
1721// break;
1722// case kTrackWithMCkine:
1723// flowtrack->Set(fMCparticle);
1724// break;
1725// case kTrackWithMCPID:
1726// flowtrack->Set(fTrack);
1727// //flowtrack->setPID(...) from mc, when implemented
1728// break;
1729// case kTrackWithMCpt:
1730// if (!fMCparticle) return NULL;
1731// flowtrack->Set(fTrack);
1732// flowtrack->SetPt(fMCparticle->Pt());
1733// break;
1734// case kTrackWithPtFromFirstMother:
1735// if (!fMCparticle) return NULL;
1736// flowtrack->Set(fTrack);
1737// tmpTParticle = fMCparticle->Particle();
1738// tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1739// flowtrack->SetPt(tmpAliMCParticle->Pt());
1740// break;
1741// case kTrackWithTPCInnerParams:
1742// esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1743// if (!esdtrack) return NULL;
1744// externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1745// if (!externalParams) return NULL;
1746// flowtrack->Set(externalParams);
1747// break;
1748// default:
1749// flowtrack->Set(fTrack);
1750// break;
1751// }
1752// if (fParamType==kMC)
1753// {
1754// flowtrack->SetSource(AliFlowTrack::kFromMC);
1755// flowtrack->SetID(fTrackLabel);
1756// }
1757// else if (dynamic_cast<AliAODTrack*>(fTrack))
1758// {
1759// if (fParamType==kMUON) // XZhang 20120604
1760// flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
1761// else // XZhang 20120604
1762// flowtrack->SetSource(AliFlowTrack::kFromAOD); // XZhang 20120604
1763// flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1764// }
1765// else if (dynamic_cast<AliMCParticle*>(fTrack))
1766// {
1767// flowtrack->SetSource(AliFlowTrack::kFromMC);
1768// flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1769// }
1770//
1771// if (fV0)
1772// {
1773// //add daughter indices
1774// }
1775//
1776// return flowtrack;
1777//}
1778
1779//-----------------------------------------------------------------------
1780AliFlowTrack* AliFlowTrackCuts::FillFlowTrackKink(TObjArray* trackCollection, Int_t trackIndex) const
1781{
1782 //fill flow track from AliVParticle (ESD,AOD,MC)
1783 AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1784 if (flowtrack)
1785 {
1786 flowtrack->Clear();
1787 }
1788 else
1789 {
1790 flowtrack = new AliFlowCandidateTrack();
1791 trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1792 }
1793
1794 TParticle *tmpTParticle=NULL;
1795 AliMCParticle* tmpAliMCParticle=NULL;
1796 AliExternalTrackParam* externalParams=NULL;
1797 AliESDtrack* esdtrack=NULL;
1798 switch(fParamMix)
1799 {
1800 case kPure:
1801 flowtrack->Set(fTrack);
1802 break;
1803 case kTrackWithMCkine:
1804 flowtrack->Set(fMCparticle);
1805 break;
1806 case kTrackWithMCPID:
1807 flowtrack->Set(fTrack);
1808 //flowtrack->setPID(...) from mc, when implemented
1809 break;
1810 case kTrackWithMCpt:
1811 if (!fMCparticle) return NULL;
1812 flowtrack->Set(fTrack);
1813 flowtrack->SetPt(fMCparticle->Pt());
1814 break;
1815 case kTrackWithPtFromFirstMother:
1816 if (!fMCparticle) return NULL;
1817 flowtrack->Set(fTrack);
1818 tmpTParticle = fMCparticle->Particle();
1819 tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1820 flowtrack->SetPt(tmpAliMCParticle->Pt());
1821 break;
1822 case kTrackWithTPCInnerParams:
1823 esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
1824 if (!esdtrack) return NULL;
1825 externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
1826 if (!externalParams) return NULL;
1827 flowtrack->Set(externalParams);
1828 break;
1829 default:
1830 flowtrack->Set(fTrack);
1831 break;
1832 }
1833 if (fParamType==kMC)
1834 {
1835 flowtrack->SetSource(AliFlowTrack::kFromMC);
1836 flowtrack->SetID(fTrackLabel);
1837 }
1838 else if (dynamic_cast<AliESDtrack*>(fTrack))
1839 {
1840 flowtrack->SetSource(AliFlowTrack::kFromESD);
1841 flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1842 }
1843 else if (dynamic_cast<AliESDMuonTrack*>(fTrack)) // XZhang 20120604
1844 { // XZhang 20120604
1845 flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
1846 flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID()); // XZhang 20120604
1847 } // XZhang 20120604
1848 else if (dynamic_cast<AliAODTrack*>(fTrack))
1849 {
1850 if (fParamType==kMUON) // XZhang 20120604
1851 flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
1852 else // XZhang 20120604
1853 flowtrack->SetSource(AliFlowTrack::kFromAOD); // XZhang 20120604
1854 flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1855 }
1856 else if (dynamic_cast<AliMCParticle*>(fTrack))
1857 {
1858 flowtrack->SetSource(AliFlowTrack::kFromMC);
1859 flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
1860 }
1861
1862 if (fKink)
1863 {
1864 Int_t indexMother = fKink->GetIndex(0);
1865 Int_t indexDaughter = fKink->GetIndex(1);
1866 flowtrack->SetID(indexMother);
1867 flowtrack->AddDaughter(indexDaughter);
1868 flowtrack->SetMass(1.);
1869 flowtrack->SetSource(AliFlowTrack::kFromKink);
1870 }
1871
1872 flowtrack->SetMass(fTrackMass);
1873
1874 return flowtrack;
1875}
1876
1877//-----------------------------------------------------------------------
1878AliFlowTrack* AliFlowTrackCuts::FillFlowTrackGeneric(TObjArray* trackCollection, Int_t trackIndex) const
1879{
1880 //fill a flow track from tracklet,vzero,pmd,...
1881 AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1882 if (flowtrack)
1883 {
1884 flowtrack->Clear();
1885 }
1886 else
1887 {
1888 flowtrack = new AliFlowTrack();
1889 trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1890 }
1891
1892 if (FillFlowTrackGeneric(flowtrack)) return flowtrack;
1893 else
1894 {
1895 trackCollection->RemoveAt(trackIndex);
1896 delete flowtrack;
1897 return NULL;
1898 }
1899}
1900
7d27a354 1901//-----------------------------------------------------------------------
1902Bool_t AliFlowTrackCuts::FillFlowTrackGeneric(AliFlowTrack* flowtrack) const
1903{
1904 //fill a flow track from tracklet,vzero,pmd,...
1905 TParticle *tmpTParticle=NULL;
1906 AliMCParticle* tmpAliMCParticle=NULL;
1907 switch (fParamMix)
1908 {
1909 case kPure:
1910 flowtrack->SetPhi(fTrackPhi);
1911 flowtrack->SetEta(fTrackEta);
2281e7c5 1912 flowtrack->SetWeight(fTrackWeight);
dd294bb6
MK		 1913 flowtrack->SetPt(fTrackPt);
1914 flowtrack->SetMass(fTrackMass);
7d27a354 1915 break;
1916 case kTrackWithMCkine:
1917 if (!fMCparticle) return kFALSE;
1918 flowtrack->SetPhi( fMCparticle->Phi() );
1919 flowtrack->SetEta( fMCparticle->Eta() );
1920 flowtrack->SetPt( fMCparticle->Pt() );
dd294bb6 1921 flowtrack->SetMass(fTrackMass);
7d27a354 1922 break;
1923 case kTrackWithMCpt:
1924 if (!fMCparticle) return kFALSE;
1925 flowtrack->SetPhi(fTrackPhi);
1926 flowtrack->SetEta(fTrackEta);
2281e7c5 1927 flowtrack->SetWeight(fTrackWeight);
7d27a354 1928 flowtrack->SetPt(fMCparticle->Pt());
dd294bb6 1929 flowtrack->SetMass(fTrackMass);
7d27a354 1930 break;
1931 case kTrackWithPtFromFirstMother:
1932 if (!fMCparticle) return kFALSE;
1933 flowtrack->SetPhi(fTrackPhi);
1934 flowtrack->SetEta(fTrackEta);
2281e7c5 1935 flowtrack->SetWeight(fTrackWeight);
7d27a354 1936 tmpTParticle = fMCparticle->Particle();
1937 tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
1938 flowtrack->SetPt(tmpAliMCParticle->Pt());
dd294bb6 1939 flowtrack->SetMass(fTrackMass);
7d27a354 1940 break;
1941 default:
1942 flowtrack->SetPhi(fTrackPhi);
1943 flowtrack->SetEta(fTrackEta);
2281e7c5 1944 flowtrack->SetWeight(fTrackWeight);
dd294bb6
MK		 1945 flowtrack->SetPt(fTrackPt);
1946 flowtrack->SetMass(fTrackMass);
7d27a354 1947 break;
1948 }
1949 flowtrack->SetSource(AliFlowTrack::kFromTracklet);
1950 return kTRUE;
1951}
1952
dd294bb6
MK		 1953//-----------------------------------------------------------------------
1954AliFlowTrack* AliFlowTrackCuts::FillFlowTrackVParticle(TObjArray* trackCollection, Int_t trackIndex) const
1955{
1956 //fill flow track from AliVParticle (ESD,AOD,MC)
1957
1958 AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1959 if (flowtrack)
1960 {
1961 flowtrack->Clear();
1962 }
1963 else
1964 {
1965 flowtrack = new AliFlowTrack();
1966 trackCollection->AddAtAndExpand(flowtrack,trackIndex);
1967 }
1968
1969 if (FillFlowTrackVParticle(flowtrack)) return flowtrack;
1970 else
1971 {
1972 trackCollection->RemoveAt(trackIndex);
1973 delete flowtrack;
1974 return NULL;
1975 }
1976}
1977
7d27a354 1978//-----------------------------------------------------------------------
1979Bool_t AliFlowTrackCuts::FillFlowTrackVParticle(AliFlowTrack* flowtrack) const
1980{
1981 //fill flow track from AliVParticle (ESD,AOD,MC)
1982 if (!fTrack) return kFALSE;
1983 TParticle *tmpTParticle=NULL;
1984 AliMCParticle* tmpAliMCParticle=NULL;
1985 AliExternalTrackParam* externalParams=NULL;
1986 AliESDtrack* esdtrack=NULL;
1987 switch(fParamMix)
1988 {
1989 case kPure:
1990 flowtrack->Set(fTrack);
1991 break;
1992 case kTrackWithMCkine:
1993 flowtrack->Set(fMCparticle);
1994 break;
1995 case kTrackWithMCPID:
1996 flowtrack->Set(fTrack);
1997 //flowtrack->setPID(...) from mc, when implemented
1998 break;
1999 case kTrackWithMCpt:
2000 if (!fMCparticle) return kFALSE;
2001 flowtrack->Set(fTrack);
2002 flowtrack->SetPt(fMCparticle->Pt());
2003 break;
2004 case kTrackWithPtFromFirstMother:
2005 if (!fMCparticle) return kFALSE;
2006 flowtrack->Set(fTrack);
2007 tmpTParticle = fMCparticle->Particle();
2008 tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2009 flowtrack->SetPt(tmpAliMCParticle->Pt());
2010 break;
2011 case kTrackWithTPCInnerParams:
2012 esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2013 if (!esdtrack) return kFALSE;
2014 externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2015 if (!externalParams) return kFALSE;
2016 flowtrack->Set(externalParams);
2017 break;
2018 default:
2019 flowtrack->Set(fTrack);
2020 break;
2021 }
499fe731 2022 if (fParamType==kMC)
2023 {
2024 flowtrack->SetSource(AliFlowTrack::kFromMC);
2281e7c5 2025 flowtrack->SetID(fTrackLabel);
499fe731 2026 }
2027 else if (dynamic_cast<AliESDtrack*>(fTrack))
2028 {
2029 flowtrack->SetSource(AliFlowTrack::kFromESD);
2030 flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2031 }
2e5052c5 2032 else if (dynamic_cast<AliESDMuonTrack*>(fTrack)) // XZhang 20120604
2033 { // XZhang 20120604
2034 flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
2035 flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID()); // XZhang 20120604
2036 } // XZhang 20120604
2037 else if (dynamic_cast<AliAODTrack*>(fTrack))
2038 {
2039 if (fParamType==kMUON) // XZhang 20120604
2040 flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
2041 else // XZhang 20120604
2042 flowtrack->SetSource(AliFlowTrack::kFromAOD); // XZhang 20120604
499fe731 2043 flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2044 }
2045 else if (dynamic_cast<AliMCParticle*>(fTrack))
2046 {
2047 flowtrack->SetSource(AliFlowTrack::kFromMC);
2048 flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2049 }
dd294bb6 2050 flowtrack->SetMass(fTrackMass);
7d27a354 2051 return kTRUE;
2052}
2053
dd294bb6
MK		 2054//-----------------------------------------------------------------------
2055AliFlowTrack* AliFlowTrackCuts::FillFlowTrack(TObjArray* trackCollection, Int_t trackIndex) const
2056{
2057 //fill a flow track constructed from whatever we applied cuts on
2058 //return true on success
2059 switch (fParamType)
2060 {
2061 case kSPDtracklet:
2062 return FillFlowTrackGeneric(trackCollection, trackIndex);
2063 case kPMD:
2064 return FillFlowTrackGeneric(trackCollection, trackIndex);
2065 case kVZERO:
2066 return FillFlowTrackGeneric(trackCollection, trackIndex);
2067 case kKink:
2068 return FillFlowTrackKink(trackCollection, trackIndex);
2069 //case kV0:
2070 // return FillFlowTrackV0(trackCollection, trackIndex);
2071 default:
2072 return FillFlowTrackVParticle(trackCollection, trackIndex);
2073 }
2074}
2075
7d27a354 2076//-----------------------------------------------------------------------
2077Bool_t AliFlowTrackCuts::FillFlowTrack(AliFlowTrack* track) const
2078{
2079 //fill a flow track constructed from whatever we applied cuts on
2080 //return true on success
2081 switch (fParamType)
2082 {
2083 case kSPDtracklet:
2084 return FillFlowTrackGeneric(track);
2085 case kPMD:
2086 return FillFlowTrackGeneric(track);
c27333c2 2087 case kVZERO:
7d27a354 2088 return FillFlowTrackGeneric(track);
2089 default:
2090 return FillFlowTrackVParticle(track);
2091 }
2092}
2093
cac040c5
MK		 2094////-----------------------------------------------------------------------
2095//AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackSPDtracklet() const
2096//{
2097// //make a flow track from tracklet
2098// AliFlowTrack* flowtrack=NULL;
2099// TParticle *tmpTParticle=NULL;
2100// AliMCParticle* tmpAliMCParticle=NULL;
2101// switch (fParamMix)
2102// {
2103// case kPure:
2104// flowtrack = new AliFlowTrack();
2105// flowtrack->SetPhi(fTrackPhi);
2106// flowtrack->SetEta(fTrackEta);
2107// flowtrack->SetWeight(fTrackWeight);
2108// break;
2109// case kTrackWithMCkine:
2110// if (!fMCparticle) return NULL;
2111// flowtrack = new AliFlowTrack();
2112// flowtrack->SetPhi( fMCparticle->Phi() );
2113// flowtrack->SetEta( fMCparticle->Eta() );
2114// flowtrack->SetPt( fMCparticle->Pt() );
2115// break;
2116// case kTrackWithMCpt:
2117// if (!fMCparticle) return NULL;
2118// flowtrack = new AliFlowTrack();
2119// flowtrack->SetPhi(fTrackPhi);
2120// flowtrack->SetEta(fTrackEta);
2121// flowtrack->SetWeight(fTrackWeight);
2122// flowtrack->SetPt(fMCparticle->Pt());
2123// break;
2124// case kTrackWithPtFromFirstMother:
2125// if (!fMCparticle) return NULL;
2126// flowtrack = new AliFlowTrack();
2127// flowtrack->SetPhi(fTrackPhi);
2128// flowtrack->SetEta(fTrackEta);
2129// flowtrack->SetWeight(fTrackWeight);
2130// tmpTParticle = fMCparticle->Particle();
2131// tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2132// flowtrack->SetPt(tmpAliMCParticle->Pt());
2133// break;
2134// default:
2135// flowtrack = new AliFlowTrack();
2136// flowtrack->SetPhi(fTrackPhi);
2137// flowtrack->SetEta(fTrackEta);
2138// flowtrack->SetWeight(fTrackWeight);
2139// break;
2140// }
2141// flowtrack->SetSource(AliFlowTrack::kFromTracklet);
2142// return flowtrack;
2143//}
2144//
2145////-----------------------------------------------------------------------
2146//AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVParticle() const
2147//{
2148// //make flow track from AliVParticle (ESD,AOD,MC)
2149// if (!fTrack) return NULL;
2150// AliFlowTrack* flowtrack=NULL;
2151// TParticle *tmpTParticle=NULL;
2152// AliMCParticle* tmpAliMCParticle=NULL;
2153// AliExternalTrackParam* externalParams=NULL;
2154// AliESDtrack* esdtrack=NULL;
2155// switch(fParamMix)
2156// {
2157// case kPure:
2158// flowtrack = new AliFlowTrack(fTrack);
2159// break;
2160// case kTrackWithMCkine:
2161// flowtrack = new AliFlowTrack(fMCparticle);
2162// break;
2163// case kTrackWithMCPID:
2164// flowtrack = new AliFlowTrack(fTrack);
2165// //flowtrack->setPID(...) from mc, when implemented
2166// break;
2167// case kTrackWithMCpt:
2168// if (!fMCparticle) return NULL;
2169// flowtrack = new AliFlowTrack(fTrack);
2170// flowtrack->SetPt(fMCparticle->Pt());
2171// break;
2172// case kTrackWithPtFromFirstMother:
2173// if (!fMCparticle) return NULL;
2174// flowtrack = new AliFlowTrack(fTrack);
2175// tmpTParticle = fMCparticle->Particle();
2176// tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2177// flowtrack->SetPt(tmpAliMCParticle->Pt());
2178// break;
2179// case kTrackWithTPCInnerParams:
2180// esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2181// if (!esdtrack) return NULL;
2182// externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2183// if (!externalParams) return NULL;
2184// flowtrack = new AliFlowTrack(externalParams);
2185// break;
2186// default:
2187// flowtrack = new AliFlowTrack(fTrack);
2188// break;
2189// }
2190// if (fParamType==kMC)
2191// {
2192// flowtrack->SetSource(AliFlowTrack::kFromMC);
2193// flowtrack->SetID(fTrackLabel);
2194// }
2195// else if (dynamic_cast<AliESDtrack*>(fTrack))
2196// {
2197// flowtrack->SetSource(AliFlowTrack::kFromESD);
2198// flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2199// }
2200// else if (dynamic_cast<AliAODTrack*>(fTrack))
2201// {
2202// flowtrack->SetSource(AliFlowTrack::kFromAOD);
2203// flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2204// }
2205// else if (dynamic_cast<AliMCParticle*>(fTrack))
2206// {
2207// flowtrack->SetSource(AliFlowTrack::kFromMC);
2208// flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2209// }
2210// return flowtrack;
2211//}
2212//
2213////-----------------------------------------------------------------------
2214//AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackPMDtrack() const
2215//{
2216// //make a flow track from PMD track
2217// AliFlowTrack* flowtrack=NULL;
2218// TParticle *tmpTParticle=NULL;
2219// AliMCParticle* tmpAliMCParticle=NULL;
2220// switch (fParamMix)
2221// {
2222// case kPure:
2223// flowtrack = new AliFlowTrack();
2224// flowtrack->SetPhi(fTrackPhi);
2225// flowtrack->SetEta(fTrackEta);
2226// flowtrack->SetWeight(fTrackWeight);
2227// break;
2228// case kTrackWithMCkine:
2229// if (!fMCparticle) return NULL;
2230// flowtrack = new AliFlowTrack();
2231// flowtrack->SetPhi( fMCparticle->Phi() );
2232// flowtrack->SetEta( fMCparticle->Eta() );
2233// flowtrack->SetWeight(fTrackWeight);
2234// flowtrack->SetPt( fMCparticle->Pt() );
2235// break;
2236// case kTrackWithMCpt:
2237// if (!fMCparticle) return NULL;
2238// flowtrack = new AliFlowTrack();
2239// flowtrack->SetPhi(fTrackPhi);
2240// flowtrack->SetEta(fTrackEta);
2241// flowtrack->SetWeight(fTrackWeight);
2242// flowtrack->SetPt(fMCparticle->Pt());
2243// break;
2244// case kTrackWithPtFromFirstMother:
2245// if (!fMCparticle) return NULL;
2246// flowtrack = new AliFlowTrack();
2247// flowtrack->SetPhi(fTrackPhi);
2248// flowtrack->SetEta(fTrackEta);
2249// flowtrack->SetWeight(fTrackWeight);
2250// tmpTParticle = fMCparticle->Particle();
2251// tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2252// flowtrack->SetPt(tmpAliMCParticle->Pt());
2253// break;
2254// default:
2255// flowtrack = new AliFlowTrack();
2256// flowtrack->SetPhi(fTrackPhi);
2257// flowtrack->SetEta(fTrackEta);
2258// flowtrack->SetWeight(fTrackWeight);
2259// break;
2260// }
2261//
2262// flowtrack->SetSource(AliFlowTrack::kFromPMD);
2263// return flowtrack;
2264//}
2265//
2266////-----------------------------------------------------------------------
2267//AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVZERO() const
2268//{
2269// //make a flow track from VZERO
2270// AliFlowTrack* flowtrack=NULL;
2271// TParticle *tmpTParticle=NULL;
2272// AliMCParticle* tmpAliMCParticle=NULL;
2273// switch (fParamMix)
2274// {
2275// case kPure:
2276// flowtrack = new AliFlowTrack();
2277// flowtrack->SetPhi(fTrackPhi);
2278// flowtrack->SetEta(fTrackEta);
2279// flowtrack->SetWeight(fTrackWeight);
2280// break;
2281// case kTrackWithMCkine:
2282// if (!fMCparticle) return NULL;
2283// flowtrack = new AliFlowTrack();
2284// flowtrack->SetPhi( fMCparticle->Phi() );
2285// flowtrack->SetEta( fMCparticle->Eta() );
2286// flowtrack->SetWeight(fTrackWeight);
2287// flowtrack->SetPt( fMCparticle->Pt() );
2288// break;
2289// case kTrackWithMCpt:
2290// if (!fMCparticle) return NULL;
2291// flowtrack = new AliFlowTrack();
2292// flowtrack->SetPhi(fTrackPhi);
2293// flowtrack->SetEta(fTrackEta);
2294// flowtrack->SetWeight(fTrackWeight);
2295// flowtrack->SetPt(fMCparticle->Pt());
2296// break;
2297// case kTrackWithPtFromFirstMother:
2298// if (!fMCparticle) return NULL;
2299// flowtrack = new AliFlowTrack();
2300// flowtrack->SetPhi(fTrackPhi);
2301// flowtrack->SetEta(fTrackEta);
2302// flowtrack->SetWeight(fTrackWeight);
2303// tmpTParticle = fMCparticle->Particle();
2304// tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2305// flowtrack->SetPt(tmpAliMCParticle->Pt());
2306// break;
2307// default:
2308// flowtrack = new AliFlowTrack();
2309// flowtrack->SetPhi(fTrackPhi);
2310// flowtrack->SetEta(fTrackEta);
2311// flowtrack->SetWeight(fTrackWeight);
2312// break;
2313// }
2314//
2315// flowtrack->SetSource(AliFlowTrack::kFromVZERO);
2316// return flowtrack;
2317//}
2318//
2319////-----------------------------------------------------------------------
2320//AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
2321//{
2322// //get a flow track constructed from whatever we applied cuts on
2323// //caller is resposible for deletion
2324// //if construction fails return NULL
2325// //TODO: for tracklets, PMD and VZERO we probably need just one method,
2326// //something like MakeFlowTrackGeneric(), wait with this until
2327// //requirements quirks are known.
2328// switch (fParamType)
2329// {
2330// case kSPDtracklet:
2331// return MakeFlowTrackSPDtracklet();
2332// case kPMD:
2333// return MakeFlowTrackPMDtrack();
2334// case kVZERO:
2335// return MakeFlowTrackVZERO();
2336// default:
2337// return MakeFlowTrackVParticle();
2338// }
2339//}
dd294bb6 2340
127a5825 2341//-----------------------------------------------------------------------
2342Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
2343{
2344 //check if current particle is a physical primary
9a0783cc 2345 if (!fMCevent) return kFALSE;
2346 if (fTrackLabel<0) return kFALSE;
441ea1cf 2347 return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
9a0783cc 2348}
2349
2350//-----------------------------------------------------------------------
441ea1cf 2351Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
9a0783cc 2352{
2353 //check if current particle is a physical primary
2354 Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
9a0783cc 2355 AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
2356 if (!track) return kFALSE;
2357 TParticle* particle = track->Particle();
2358 Bool_t transported = particle->TestBit(kTransportBit);
441ea1cf 2359 //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
2360 //(physprim && (transported || !requiretransported))?"YES":"NO" );
2361 return (physprim && (transported || !requiretransported));
127a5825 2362}
12b2b8bc 2363
924b02b0 2364//-----------------------------------------------------------------------
2365void AliFlowTrackCuts::DefineHistograms()
2366{
2a745a5f 2367 //define qa histograms
1a80f9f6 2368 if (fQA) return;
2369
2281e7c5 2370 const Int_t kNbinsP=200;
499fe731 2371 Double_t binsP[kNbinsP+1];
2372 binsP[0]=0.0;
6a043318 2373 for(int i=1; i<kNbinsP+1; i++)
1a80f9f6 2374 {
2281e7c5 2375 //if(binsP[i-1]+0.05<1.01)
2376 // binsP[i]=binsP[i-1]+0.05;
2377 //else
2378 binsP[i]=binsP[i-1]+0.05;
1a80f9f6 2379 }
2380
cd62a2a8 2381 const Int_t nBinsDCA=1000;
2382 Double_t binsDCA[nBinsDCA+1];
b9cf8f8e 2383 for(int i=0; i<nBinsDCA+1; i++) {binsDCA[i]=0.01*i-5.;}
cd62a2a8 2384 //for(int i=1; i<41; i++) {binsDCA[i+100]=0.1*i+1.0;}
2385
499fe731 2386 Bool_t adddirstatus = TH1::AddDirectoryStatus();
2387 TH1::AddDirectory(kFALSE);
a14b8f3c 2388 fQA=new TList(); fQA->SetOwner();
499fe731 2389 fQA->SetName(Form("%s QA",GetName()));
a14b8f3c 2390 TList* before = new TList(); before->SetOwner();
499fe731 2391 before->SetName("before");
a14b8f3c 2392 TList* after = new TList(); after->SetOwner();
499fe731 2393 after->SetName("after");
2394 fQA->Add(before);
2395 fQA->Add(after);
f21bdf48 2396 before->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2397 after->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2398 before->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2399 after->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2400 before->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
2401 after->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
483a17da 2402 //primary
2403 TH2F* hb = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2404 TH2F* ha = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2405 TAxis* axis = NULL;
2406 axis = hb->GetYaxis();
2407 axis->SetBinLabel(1,"secondary");
2408 axis->SetBinLabel(2,"primary");
2409 axis = ha->GetYaxis();
2410 axis->SetBinLabel(1,"secondary");
2411 axis->SetBinLabel(2,"primary");
2412 before->Add(hb); //3
2413 after->Add(ha); //3
875ee12a 2414 //production process
483a17da 2415 hb = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
875ee12a 2416 -0.5, kMaxMCProcess-0.5);
483a17da 2417 ha = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
875ee12a 2418 -0.5, kMaxMCProcess-0.5);
483a17da 2419 axis = hb->GetYaxis();
875ee12a 2420 for (Int_t i=0; i<kMaxMCProcess; i++)
2421 {
2422 axis->SetBinLabel(i+1,TMCProcessName[i]);
2423 }
644f1368 2424 axis = ha->GetYaxis();
875ee12a 2425 for (Int_t i=0; i<kMaxMCProcess; i++)
2426 {
2427 axis->SetBinLabel(i+1,TMCProcessName[i]);
2428 }
2429 before->Add(hb); //4
2430 after->Add(ha); //4
cd62a2a8 2431 //DCA
2432 before->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2433 after->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2434 before->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
2435 after->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
8eca5d19 2436 //first mother
2281e7c5 2437 hb = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2438 ha = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2439 hb->GetYaxis()->SetBinLabel(1, "#gamma");
2440 ha->GetYaxis()->SetBinLabel(1, "#gamma");
2441 hb->GetYaxis()->SetBinLabel(2, "e^{+}");
2442 ha->GetYaxis()->SetBinLabel(2, "e^{+}");
2443 hb->GetYaxis()->SetBinLabel(3, "e^{-}");
2444 ha->GetYaxis()->SetBinLabel(3, "e^{-}");
2445 hb->GetYaxis()->SetBinLabel(4, "#nu");
2446 ha->GetYaxis()->SetBinLabel(4, "#nu");
2447 hb->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2448 ha->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2449 hb->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2450 ha->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2451 hb->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2452 ha->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2453 hb->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2454 ha->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2455 hb->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2456 ha->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2457 hb->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2458 ha->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2459 hb->GetYaxis()->SetBinLabel(11, "K^{+}");
2460 ha->GetYaxis()->SetBinLabel(11, "K^{+}");
2461 hb->GetYaxis()->SetBinLabel(12, "K^{-}");
2462 ha->GetYaxis()->SetBinLabel(12, "K^{-}");
2463 hb->GetYaxis()->SetBinLabel( 13, "n");
2464 ha->GetYaxis()->SetBinLabel( 13, "n");
2465 hb->GetYaxis()->SetBinLabel( 14, "p");
2466 ha->GetYaxis()->SetBinLabel( 14, "p");
2467 hb->GetYaxis()->SetBinLabel(15, "#bar{p}");
2468 ha->GetYaxis()->SetBinLabel(15, "#bar{p}");
2469 hb->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2470 ha->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2471 hb->GetYaxis()->SetBinLabel(17, "#eta");
2472 ha->GetYaxis()->SetBinLabel(17, "#eta");
2473 hb->GetYaxis()->SetBinLabel(18, "#Lambda");
2474 ha->GetYaxis()->SetBinLabel(18, "#Lambda");
2475 hb->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2476 ha->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2477 hb->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2478 ha->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2479 hb->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2480 ha->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2481 hb->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2482 ha->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2483 hb->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2484 ha->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2485 hb->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2486 ha->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2487 hb->GetYaxis()->SetBinLabel(25, "#bar{n}");
2488 ha->GetYaxis()->SetBinLabel(25, "#bar{n}");
2489 hb->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2490 ha->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2491 hb->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2492 ha->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2493 hb->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2494 ha->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2495 hb->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2496 ha->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2497 hb->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2498 ha->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2499 hb->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2500 ha->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2501 hb->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2502 ha->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2503 hb->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2504 ha->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2505 hb->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2506 ha->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2507 hb->GetYaxis()->SetBinLabel(35, "D^{+}");
2508 ha->GetYaxis()->SetBinLabel(35, "D^{+}");
2509 hb->GetYaxis()->SetBinLabel(36, "D^{-}");
2510 ha->GetYaxis()->SetBinLabel(36, "D^{-}");
2511 hb->GetYaxis()->SetBinLabel(37, "D^{0}");
2512 ha->GetYaxis()->SetBinLabel(37, "D^{0}");
2513 hb->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2514 ha->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2515 hb->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2516 ha->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2517 hb->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2518 ha->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2519 hb->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2520 ha->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2521 hb->GetYaxis()->SetBinLabel(42, "W^{+}");
2522 ha->GetYaxis()->SetBinLabel(42, "W^{+}");
2523 hb->GetYaxis()->SetBinLabel(43, "W^{-}");
2524 ha->GetYaxis()->SetBinLabel(43, "W^{-}");
2525 hb->GetYaxis()->SetBinLabel(44, "Z^{0}");
2526 ha->GetYaxis()->SetBinLabel(44, "Z^{0}");
8eca5d19 2527 before->Add(hb);//7
2528 after->Add(ha);//7
2281e7c5 2529
60875c3c 2530 before->Add(new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2531 after->Add( new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2532
2533 before->Add(new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//9
2534 after->Add( new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//9
2535
2536 before->Add(new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//10
2537 after->Add( new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//10
2538
2539 before->Add(new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//11
2540 after->Add( new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//11
2541
2542 before->Add(new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//12
2543 after->Add( new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//12
2544
dd294bb6
MK		 2545 //kink stuff
2546 before->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2547 after->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2548
2549 before->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2550 after->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2551
2552 before->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2553 after->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2554
2555 before->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2556 after->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2557
2558 before->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
2559 after->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
67006c99 2560 /*
2561 before->Add(new TH2F("Pionbeta",";p_{t}[GeV/c];#beta",kNbinsP,binsP,1000,0,1000 ));//18
2562 after->Add(new TH2F("Pionbeta",";p_{t}[GeV/c];#beta",kNbinsP,binsP,1000,0,1000 ));//18
2563
2564 before->Add(new TH2F("PiondEdX",";p_{t}[GeV/c];dE/dX",kNbinsP,binsP,1000,0,1000 ));//19
2565 after->Add(new TH2F("PiondEdX",";p_{t}[GeV/c];dE/dX",kNbinsP,binsP,1000,0,1000 ));//19
2566*/
499fe731 2567 TH1::AddDirectory(adddirstatus);
924b02b0 2568}
9a0783cc 2569
2570//-----------------------------------------------------------------------
2571Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
2572{
2573 //get the number of tracks in the input event according source
2574 //selection (ESD tracks, tracklets, MC particles etc.)
2575 AliESDEvent* esd=NULL;
2e5052c5 2576 AliAODEvent* aod=NULL; // XZhang 20120615
9a0783cc 2577 switch (fParamType)
2578 {
1a80f9f6 2579 case kSPDtracklet:
2e5052c5 2580 if (!fEvent) return 0; // XZhang 20120615
9a0783cc 2581 esd = dynamic_cast<AliESDEvent*>(fEvent);
2e5052c5 2582 aod = dynamic_cast<AliAODEvent*>(fEvent); // XZhang 20120615
2583// if (!esd) return 0; // XZhang 20120615
2584// return esd->GetMultiplicity()->GetNumberOfTracklets(); // XZhang 20120615
2585 if (esd) return esd->GetMultiplicity()->GetNumberOfTracklets(); // XZhang 20120615
2586 if (aod) return aod->GetTracklets()->GetNumberOfTracklets(); // XZhang 20120615
9a0783cc 2587 case kMC:
2588 if (!fMCevent) return 0;
2589 return fMCevent->GetNumberOfTracks();
d148af7e 2590 case kPMD:
2591 esd = dynamic_cast<AliESDEvent*>(fEvent);
2592 if (!esd) return 0;
2593 return esd->GetNumberOfPmdTracks();
c27333c2
MK		 2594 case kVZERO:
2595 return fgkNumberOfVZEROtracks;
2e5052c5 2596 case kMUON: // XZhang 20120604
2597 if (!fEvent) return 0; // XZhang 20120604
2598 esd = dynamic_cast<AliESDEvent*>(fEvent); // XZhang 20120604
2599 if (esd) return esd->GetNumberOfMuonTracks(); // XZhang 20120604
2600 return fEvent->GetNumberOfTracks(); // if AOD // XZhang 20120604
dd294bb6
MK		 2601 case kKink:
2602 esd = dynamic_cast<AliESDEvent*>(fEvent);
2603 if (!esd) return 0;
2604 return esd->GetNumberOfKinks();
2605 case kV0:
2606 esd = dynamic_cast<AliESDEvent*>(fEvent);
2607 if (!esd) return 0;
2608 return esd->GetNumberOfV0s();
9a0783cc 2609 default:
2610 if (!fEvent) return 0;
2611 return fEvent->GetNumberOfTracks();
2612 }
2613 return 0;
2614}
2615
2616//-----------------------------------------------------------------------
2617TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
2618{
2619 //get the input object according the data source selection:
2620 //(esd tracks, traclets, mc particles,etc...)
2621 AliESDEvent* esd=NULL;
2e5052c5 2622 AliAODEvent* aod=NULL; // XZhang 20120615
9a0783cc 2623 switch (fParamType)
2624 {
1a80f9f6 2625 case kSPDtracklet:
2e5052c5 2626 if (!fEvent) return NULL; // XZhang 20120615
9a0783cc 2627 esd = dynamic_cast<AliESDEvent*>(fEvent);
2e5052c5 2628 aod = dynamic_cast<AliAODEvent*>(fEvent); // XZhang 20120615
2629// if (!esd) return NULL; // XZhang 20120615
2630// return const_cast<AliMultiplicity*>(esd->GetMultiplicity()); // XZhang 20120615
2631 if (esd) return const_cast<AliMultiplicity*>(esd->GetMultiplicity()); // XZhang 20120615
2632 if (aod) return const_cast<AliAODTracklets*>(aod->GetTracklets()); // XZhang 20120615
9a0783cc 2633 case kMC:
2634 if (!fMCevent) return NULL;
2635 return fMCevent->GetTrack(i);
d148af7e 2636 case kPMD:
2637 esd = dynamic_cast<AliESDEvent*>(fEvent);
2638 if (!esd) return NULL;
2639 return esd->GetPmdTrack(i);
c27333c2 2640 case kVZERO:
dd294bb6
MK		 2641 esd = dynamic_cast<AliESDEvent*>(fEvent);
2642 if (!esd) //contributed by G.Ortona
2643 {
2644 aod = dynamic_cast<AliAODEvent*>(fEvent);
2645 if(!aod)return NULL;
2646 return aod->GetVZEROData();
2647 }
2648 return esd->GetVZEROData();
2e5052c5 2649 case kMUON: // XZhang 20120604
2650 if (!fEvent) return NULL; // XZhang 20120604
2651 esd = dynamic_cast<AliESDEvent*>(fEvent); // XZhang 20120604
2652 if (esd) return esd->GetMuonTrack(i); // XZhang 20120604
2653 return fEvent->GetTrack(i); // if AOD // XZhang 20120604
dd294bb6
MK		 2654 case kKink:
2655 esd = dynamic_cast<AliESDEvent*>(fEvent);
2656 if (!esd) return NULL;
2657 return esd->GetKink(i);
2658 case kV0:
2659 esd = dynamic_cast<AliESDEvent*>(fEvent);
2660 if (!esd) return NULL;
2661 return esd->GetV0(i);
9a0783cc 2662 default:
2663 if (!fEvent) return NULL;
2664 return fEvent->GetTrack(i);
2665 }
2666}
a1c43d26 2667
2668//-----------------------------------------------------------------------
2669void AliFlowTrackCuts::Clear(Option_t*)
2670{
2671 //clean up
a1c43d26 2672 fMCevent=NULL;
dd294bb6
MK		 2673 fEvent=NULL;
2674 ClearTrack();
2675}
2676
2677//-----------------------------------------------------------------------
2678void AliFlowTrackCuts::ClearTrack(Option_t*)
2679{
2680 //clean up last track
2681 fKink=NULL;
2682 fV0=NULL;
2683 fTrack=NULL;
a1c43d26 2684 fMCparticle=NULL;
2281e7c5 2685 fTrackLabel=-997;
2686 fTrackWeight=1.0;
a1c43d26 2687 fTrackEta=0.0;
2688 fTrackPhi=0.0;
dd294bb6
MK		 2689 fTrackPt=0.0;
2690 fPOItype=1;
2691 fTrackMass=0.;
a1c43d26 2692}
0d16f553 2693//-----------------------------------------------------------------------
2694Bool_t AliFlowTrackCuts::PassesAODpidCut(const AliAODTrack* track )
2695{
2696 if(!track->GetAODEvent()->GetTOFHeader()){
2697 AliAODPid *pidObj = track->GetDetPid();
2698 if (!pidObj) fESDpid.GetTOFResponse().SetTimeResolution(84.);
2699 else{
2700 Double_t sigmaTOFPidInAOD[10];
2701 pidObj->GetTOFpidResolution(sigmaTOFPidInAOD);
2702 if(sigmaTOFPidInAOD[0] > 84.){
2703 fESDpid.GetTOFResponse().SetTimeResolution(sigmaTOFPidInAOD[0]); // use the electron TOF PID sigma as time resolution (including the T0 used)
2704 }
2705 }
2706 }
2707
2708 //check if passes the selected pid cut for ESDs
2709 Bool_t pass = kTRUE;
2710 switch (fPIDsource)
2711 {
2712 case kTOFbeta:
2713 if (!PassesTOFbetaCut(track)) pass=kFALSE;
2714 break;
2715 case kTOFbayesian:
2716 if (!PassesTOFbayesianCut(track)) pass=kFALSE;
2717 break;
2718 case kTPCbayesian:
2719 if (!PassesTPCbayesianCut(track)) pass=kFALSE;
2720 break;
043e9010 2721 case kTPCTOFNsigma:
2722 if (!PassesTPCTOFNsigmaCut(track)) pass = kFALSE;
2723 break;
67006c99 2724 case kTPCTOFNsigmaPuritybased:
2725 if(!PassesTPCTOFNsigmaCutPuritybased(track)) pass = kFALSE;
2726 break;
0d16f553 2727 default:
2728 return kTRUE;
2729 break;
2730 }
2731 return pass;
32b846cd 2732
0d16f553 2733}
3c67c769 2734//-----------------------------------------------------------------------
2735Bool_t AliFlowTrackCuts::PassesESDpidCut(const AliESDtrack* track )
2736{
2737 //check if passes the selected pid cut for ESDs
2738 Bool_t pass = kTRUE;
2739 switch (fPIDsource)
2740 {
2741 case kTPCpid:
2742 if (!PassesTPCpidCut(track)) pass=kFALSE;
2743 break;
2744 case kTPCdedx:
2745 if (!PassesTPCdedxCut(track)) pass=kFALSE;
2746 break;
2747 case kTOFpid:
2748 if (!PassesTOFpidCut(track)) pass=kFALSE;
2749 break;
2750 case kTOFbeta:
2751 if (!PassesTOFbetaCut(track)) pass=kFALSE;
2752 break;
2753 case kTOFbetaSimple:
2754 if (!PassesTOFbetaSimpleCut(track)) pass=kFALSE;
2755 break;
2756 case kTPCbayesian:
2757 if (!PassesTPCbayesianCut(track)) pass=kFALSE;
2758 break;
2759 // part added by F. Noferini
2760 case kTOFbayesian:
2761 if (!PassesTOFbayesianCut(track)) pass=kFALSE;
2762 break;
2763 // end part added by F. Noferini
2764
2765 //part added by Natasha
2766 case kTPCNuclei:
2767 if (!PassesNucleiSelection(track)) pass=kFALSE;
2768 break;
2769 //end part added by Natasha
d0ce6147 2770
2771 case kTPCTOFNsigma:
2772 if (!PassesTPCTOFNsigmaCut(track)) pass = kFALSE;
2773 break;
3c67c769 2774 default:
2775 printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
2776 pass=kFALSE;
2777 break;
2778 }
2779 return pass;
2780}
2781
32b846cd 2782//-----------------------------------------------------------------------
1a80f9f6 2783Bool_t AliFlowTrackCuts::PassesTOFbetaSimpleCut(const AliESDtrack* track )
32b846cd 2784{
2785 //check if passes PID cut using timing in TOF
8b649d8c 2786 Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
b092f828 2787 (track->GetTOFsignal() > 12000) &&
2788 (track->GetTOFsignal() < 100000) &&
499fe731 2789 (track->GetIntegratedLength() > 365);
2790
875ee12a 2791 if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2792
2793 Bool_t statusMatchingHard = TPCTOFagree(track);
2794 if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2795 return kFALSE;
aab6527a 2796
2797 if (!goodtrack) return kFALSE;
2798
a0241c3a 2799 const Float_t c = 2.99792457999999984e-02;
b092f828 2800 Float_t p = track->GetP();
1a80f9f6 2801 Float_t l = track->GetIntegratedLength();
aab6527a 2802 Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
b092f828 2803 Float_t timeTOF = track->GetTOFsignal()- trackT0;
1a80f9f6 2804 Float_t beta = l/timeTOF/c;
80efa607 2805 Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
1a80f9f6 2806 track->GetIntegratedTimes(integratedTimes);
2807 Float_t betaHypothesis[5] = {0.0,0.0,0.0,0.0,0.0};
2808 Float_t s[5] = {0.0,0.0,0.0,0.0,0.0};
2809 for (Int_t i=0;i<5;i++)
2810 {
2811 betaHypothesis[i] = l/integratedTimes[i]/c;
2812 s[i] = beta-betaHypothesis[i];
2813 }
2814
2815 switch (fParticleID)
2816 {
2817 case AliPID::kPion:
2818 return ( (s[2]<0.015) && (s[2]>-0.015) &&
7288f660 2819 (s[3]>0.025) &&
2820 (s[4]>0.03) );
1a80f9f6 2821 case AliPID::kKaon:
2822 return ( (s[3]<0.015) && (s[3]>-0.015) &&
7288f660 2823 (s[2]<-0.03) &&
2824 (s[4]>0.03) );
1a80f9f6 2825 case AliPID::kProton:
499fe731 2826 return ( (s[4]<0.015) && (s[4]>-0.015) &&
7288f660 2827 (s[3]<-0.025) &&
2828 (s[2]<-0.025) );
1a80f9f6 2829 default:
2830 return kFALSE;
2831 }
2832 return kFALSE;
2833}
2834
499fe731 2835//-----------------------------------------------------------------------
a22c1eb1 2836Float_t AliFlowTrackCuts::GetBeta(const AliVTrack* track, Bool_t QAmode)
499fe731 2837{
2838 //get beta
80efa607 2839 Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
0d16f553 2840 track->GetIntegratedTimes(integratedTimes);
2841
499fe731 2842 const Float_t c = 2.99792457999999984e-02;
0d16f553 2843 Float_t p = track->P();
2844 Float_t l = integratedTimes[0]*c;
499fe731 2845 Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
2846 Float_t timeTOF = track->GetTOFsignal()- trackT0;
a22c1eb1 2847 if(QAmode && timeTOF <= 0) return -999; // avoid division by zero when filling 'before' qa histograms
499fe731 2848 return l/timeTOF/c;
2849}
0d16f553 2850//-----------------------------------------------------------------------
2851Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliAODTrack* track )
2852{
2853 //check if track passes pid selection with an asymmetric TOF beta cut
2854 if (!fTOFpidCuts)
2855 {
2856 //printf("no TOFpidCuts\n");
2857 return kFALSE;
2858 }
2859
2860 //check if passes PID cut using timing in TOF
2861 Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
2862 (track->GetTOFsignal() > 12000) &&
2863 (track->GetTOFsignal() < 100000);
2864
2865 if (!goodtrack) return kFALSE;
2866
2867 const Float_t c = 2.99792457999999984e-02;
80efa607 2868 Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
0d16f553 2869 track->GetIntegratedTimes(integratedTimes);
2870 Float_t l = integratedTimes[0]*c;
499fe731 2871
0d16f553 2872 goodtrack = goodtrack && (l > 365);
2873
2874 if (!goodtrack) return kFALSE;
2875
2876 if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2877
2878 Bool_t statusMatchingHard = TPCTOFagree(track);
2879 if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2880 return kFALSE;
2881
2882
2883 Float_t beta = GetBeta(track);
2884
2885 //construct the pid index because it's not AliPID::EParticleType
2886 Int_t pid = 0;
67006c99 2887 cout<<"TOFbeta: fParticleID = "<<fParticleID<<endl;
0d16f553 2888 switch (fParticleID)
2889 {
2890 case AliPID::kPion:
2891 pid=2;
2892 break;
2893 case AliPID::kKaon:
2894 pid=3;
2895 break;
2896 case AliPID::kProton:
2897 pid=4;
2898 break;
2899 default:
2900 return kFALSE;
2901 }
2902
2903 //signal to cut on
2904 Float_t p = track->P();
2905 Float_t betahypothesis = l/integratedTimes[pid]/c;
2906 Float_t betadiff = beta-betahypothesis;
2907
2908 Float_t* arr = fTOFpidCuts->GetMatrixArray();
2909 Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
2910 if (col<0) return kFALSE;
2911 Float_t min = (*fTOFpidCuts)(1,col);
2912 Float_t max = (*fTOFpidCuts)(2,col);
2913
2914 Bool_t pass = (betadiff>min && betadiff<max);
2915
2916 return pass;
2917}
1a80f9f6 2918//-----------------------------------------------------------------------
2919Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliESDtrack* track )
2920{
3c67c769 2921 //check if track passes pid selection with an asymmetric TOF beta cut
9b692328 2922 if (!fTOFpidCuts)
2923 {
2924 //printf("no TOFpidCuts\n");
2925 return kFALSE;
2926 }
2927
1a80f9f6 2928 //check if passes PID cut using timing in TOF
a14b8f3c 2929 Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
1a80f9f6 2930 (track->GetTOFsignal() > 12000) &&
2931 (track->GetTOFsignal() < 100000) &&
499fe731 2932 (track->GetIntegratedLength() > 365);
2933
875ee12a 2934 if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
2935
0d16f553 2936 if (!goodtrack) return kFALSE;
2937
875ee12a 2938 Bool_t statusMatchingHard = TPCTOFagree(track);
2939 if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
2940 return kFALSE;
1a80f9f6 2941
499fe731 2942 Float_t beta = GetBeta(track);
80efa607 2943 Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
a0241c3a 2944 track->GetIntegratedTimes(integratedTimes);
2945
2946 //construct the pid index because it's not AliPID::EParticleType
32b846cd 2947 Int_t pid = 0;
2b1eaa10 2948 switch (fParticleID)
32b846cd 2949 {
2950 case AliPID::kPion:
2951 pid=2;
2952 break;
2953 case AliPID::kKaon:
2954 pid=3;
2955 break;
2956 case AliPID::kProton:
2957 pid=4;
2958 break;
2959 default:
2960 return kFALSE;
2961 }
a0241c3a 2962
2963 //signal to cut on
499fe731 2964 const Float_t c = 2.99792457999999984e-02;
2965 Float_t l = track->GetIntegratedLength();
2966 Float_t p = track->GetP();
1a80f9f6 2967 Float_t betahypothesis = l/integratedTimes[pid]/c;
2968 Float_t betadiff = beta-betahypothesis;
32b846cd 2969
2970 Float_t* arr = fTOFpidCuts->GetMatrixArray();
a0241c3a 2971 Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
32b846cd 2972 if (col<0) return kFALSE;
2973 Float_t min = (*fTOFpidCuts)(1,col);
2974 Float_t max = (*fTOFpidCuts)(2,col);
2975
1a80f9f6 2976 Bool_t pass = (betadiff>min && betadiff<max);
2977
1a80f9f6 2978 return pass;
32b846cd 2979}
2980
2b1eaa10 2981//-----------------------------------------------------------------------
7d27a354 2982Bool_t AliFlowTrackCuts::PassesTOFpidCut(const AliESDtrack* track) const
2b1eaa10 2983{
2984 //check if passes PID cut using default TOF pid
7d27a354 2985 Double_t pidTOF[AliPID::kSPECIES];
2986 track->GetTOFpid(pidTOF);
2987 if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
2b1eaa10 2988 return kFALSE;
2989}
2990
32b846cd 2991//-----------------------------------------------------------------------
1a80f9f6 2992Bool_t AliFlowTrackCuts::PassesTPCpidCut(const AliESDtrack* track) const
2b1eaa10 2993{
2994 //check if passes PID cut using default TPC pid
2995 Double_t pidTPC[AliPID::kSPECIES];
2996 track->GetTPCpid(pidTPC);
2997 Double_t probablity = 0.;
2998 switch (fParticleID)
2999 {
3000 case AliPID::kPion:
3001 probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
3002 break;
3003 default:
3004 probablity = pidTPC[fParticleID];
3005 }
3006 if (probablity >= fParticleProbability) return kTRUE;
3007 return kFALSE;
3008}
3009
499fe731 3010//-----------------------------------------------------------------------
3c67c769 3011Float_t AliFlowTrackCuts::Getdedx(const AliESDtrack* track) const
499fe731 3012{
3013 //get TPC dedx
3014 return track->GetTPCsignal();
3015}
3016
2b1eaa10 3017//-----------------------------------------------------------------------
1a80f9f6 3018Bool_t AliFlowTrackCuts::PassesTPCdedxCut(const AliESDtrack* track)
32b846cd 3019{
3020 //check if passes PID cut using dedx signal in the TPC
32b846cd 3021 if (!fTPCpidCuts)
3022 {
9b692328 3023 //printf("no TPCpidCuts\n");
32b846cd 3024 return kFALSE;
3025 }
3026
3c75bbd5 3027 const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
3028 if (!tpcparam) return kFALSE;
d79d61fb 3029 Double_t p = tpcparam->GetP();
3030 Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(p, fParticleID);
32b846cd 3031 Float_t sigTPC = track->GetTPCsignal();
3032 Float_t s = (sigTPC-sigExp)/sigExp;
32b846cd 3033
3034 Float_t* arr = fTPCpidCuts->GetMatrixArray();
1a80f9f6 3035 Int_t arrSize = fTPCpidCuts->GetNcols();
3036 Int_t col = TMath::BinarySearch( arrSize, arr, static_cast<Float_t>(p));
32b846cd 3037 if (col<0) return kFALSE;
3038 Float_t min = (*fTPCpidCuts)(1,col);
3039 Float_t max = (*fTPCpidCuts)(2,col);
3040
3041 //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
3042 return (s>min && s<max);
3043}
3044
3045//-----------------------------------------------------------------------
3046void AliFlowTrackCuts::InitPIDcuts()
3047{
3048 //init matrices with PID cuts
3049 TMatrixF* t = NULL;
3050 if (!fTPCpidCuts)
3051 {
2b1eaa10 3052 if (fParticleID==AliPID::kPion)
32b846cd 3053 {
d79d61fb 3054 t = new TMatrixF(3,15);
3055 (*t)(0,0) = 0.20; (*t)(1,0) = -0.4; (*t)(2,0) = 0.0;
3056 (*t)(0,1) = 0.25; (*t)(1,1) = -0.4; (*t)(2,1) = 0.1;
3057 (*t)(0,2) = 0.30; (*t)(1,2) = -0.4; (*t)(2,2) = 0.2;
3058 (*t)(0,3) = 0.35; (*t)(1,3) = -0.4; (*t)(2,3) = 0.2;
32b846cd 3059 (*t)(0,4) = 0.40; (*t)(1,4) = -0.4; (*t)(2,4) = 0.3;
3060 (*t)(0,5) = 0.45; (*t)(1,5) = -0.4; (*t)(2,5) = 0.3;
3061 (*t)(0,6) = 0.50; (*t)(1,6) = -0.4; (*t)(2,6) = 0.25;
3062 (*t)(0,7) = 0.55; (*t)(1,7) = -0.4; (*t)(2,7) = 0.15;
3063 (*t)(0,8) = 0.60; (*t)(1,8) = -0.4; (*t)(2,8) = 0.1;
d79d61fb 3064 (*t)(0,9) = 0.65; (*t)(1,9) = -0.4; (*t)(2,9) = 0.05;
3065 (*t)(0,10) = 0.70; (*t)(1,10) = -0.4; (*t)(2,10) = 0;
3066 (*t)(0,11) = 0.75; (*t)(1,11) = -0.4; (*t)(2,11) = 0;
3067 (*t)(0,12) = 0.80; (*t)(1,12) = -0.4; (*t)(2,12) = -0.05;
3068 (*t)(0,13) = 0.85; (*t)(1,13) = -0.4; (*t)(2,13) = -0.1;
3069 (*t)(0,14) = 0.90; (*t)(1,14) = 0; (*t)(2,14) = 0;
32b846cd 3070 }
3071 else
2b1eaa10 3072 if (fParticleID==AliPID::kKaon)
3ca4688a 3073 {
d79d61fb 3074 t = new TMatrixF(3,12);
3075 (*t)(0,0) = 0.20; (*t)(1,0) = -0.2; (*t)(2,0) = 0.2;
3076 (*t)(0,1) = 0.25; (*t)(1,1) = -0.2; (*t)(2,1) = 0.2;
3077 (*t)(0,2) = 0.30; (*t)(1,2) = -0.2; (*t)(2,2) = 0.2;
3078 (*t)(0,3) = 0.35; (*t)(1,3) = -0.2; (*t)(2,3) = 0.2;
3079 (*t)(0,4) = 0.40; (*t)(1,4) = -0.1; (*t)(2,4) = 0.2;
3080 (*t)(0,5) = 0.45; (*t)(1,5) = -0.1; (*t)(2,5) = 0.2;
3081 (*t)(0,6) = 0.50; (*t)(1,6) =-0.05; (*t)(2,6) = 0.2;
3082 (*t)(0,7) = 0.55; (*t)(1,7) = -0.1; (*t)(2,7) = 0.1;
3083 (*t)(0,8) = 0.60; (*t)(1,8) =-0.05; (*t)(2,8) = 0.1;
3084 (*t)(0,9) = 0.65; (*t)(1,9) = 0; (*t)(2,9) = 0.15;
3085 (*t)(0,10) = 0.70; (*t)(1,10) = 0.05; (*t)(2,10) = 0.2;
1a80f9f6 3086 (*t)(0,11) = 0.75; (*t)(1,11) = 0; (*t)(2,11) = 0;
3ca4688a 3087 }
3088 else
2b1eaa10 3089 if (fParticleID==AliPID::kProton)
32b846cd 3090 {
d79d61fb 3091 t = new TMatrixF(3,9);
3092 (*t)(0,0) = 0.20; (*t)(1,0) = -0.1; (*t)(2,0) = 0.1;
3093 (*t)(0,1) = 0.25; (*t)(1,1) = -0.2; (*t)(2,1) = 0.2;
3094 (*t)(0,2) = 0.80; (*t)(1,2) = -0.1; (*t)(2,2) = 0.2;
3095 (*t)(0,3) = 0.85; (*t)(1,3) =-0.05; (*t)(2,3) = 0.2;
3096 (*t)(0,4) = 0.90; (*t)(1,4) =-0.05; (*t)(2,4) = 0.25;
3097 (*t)(0,5) = 0.95; (*t)(1,5) =-0.05; (*t)(2,5) = 0.25;
3098 (*t)(0,6) = 1.00; (*t)(1,6) = -0.1; (*t)(2,6) = 0.25;
3099 (*t)(0,7) = 1.10; (*t)(1,7) =-0.05; (*t)(2,7) = 0.3;
3100 (*t)(0,8) = 1.20; (*t)(1,8) = 0; (*t)(2,8) = 0;
32b846cd 3101 }
d79d61fb 3102 delete fTPCpidCuts;
32b846cd 3103 fTPCpidCuts=t;
3104 }
7d9ab4fb 3105 t = NULL;
32b846cd 3106 if (!fTOFpidCuts)
3107 {
2b1eaa10 3108 if (fParticleID==AliPID::kPion)
32b846cd 3109 {
a0241c3a 3110 //TOF pions, 0.9 purity
3111 t = new TMatrixF(3,61);
d79d61fb 3112 (*t)(0,0) = 0.000; (*t)(1,0) = 0.000; (*t)(2,0) = 0.000;
3113 (*t)(0,1) = 0.050; (*t)(1,1) = 0.000; (*t)(2,1) = 0.000;
3114 (*t)(0,2) = 0.100; (*t)(1,2) = 0.000; (*t)(2,2) = 0.000;
3115 (*t)(0,3) = 0.150; (*t)(1,3) = 0.000; (*t)(2,3) = 0.000;
1a80f9f6 3116 (*t)(0,4) = 0.200; (*t)(1,4) = -0.030; (*t)(2,4) = 0.030;
d79d61fb 3117 (*t)(0,5) = 0.250; (*t)(1,5) = -0.036; (*t)(2,5) = 0.032;
3118 (*t)(0,6) = 0.300; (*t)(1,6) = -0.038; (*t)(2,6) = 0.032;
3119 (*t)(0,7) = 0.350; (*t)(1,7) = -0.034; (*t)(2,7) = 0.032;
3120 (*t)(0,8) = 0.400; (*t)(1,8) = -0.032; (*t)(2,8) = 0.020;
3121 (*t)(0,9) = 0.450; (*t)(1,9) = -0.030; (*t)(2,9) = 0.020;
3122 (*t)(0,10) = 0.500; (*t)(1,10) = -0.030; (*t)(2,10) = 0.020;
3123 (*t)(0,11) = 0.550; (*t)(1,11) = -0.030; (*t)(2,11) = 0.020;
3124 (*t)(0,12) = 0.600; (*t)(1,12) = -0.030; (*t)(2,12) = 0.020;
3125 (*t)(0,13) = 0.650; (*t)(1,13) = -0.030; (*t)(2,13) = 0.020;
3126 (*t)(0,14) = 0.700; (*t)(1,14) = -0.030; (*t)(2,14) = 0.020;
3127 (*t)(0,15) = 0.750; (*t)(1,15) = -0.030; (*t)(2,15) = 0.020;
3128 (*t)(0,16) = 0.800; (*t)(1,16) = -0.030; (*t)(2,16) = 0.020;
3129 (*t)(0,17) = 0.850; (*t)(1,17) = -0.030; (*t)(2,17) = 0.020;
3130 (*t)(0,18) = 0.900; (*t)(1,18) = -0.030; (*t)(2,18) = 0.020;
3131 (*t)(0,19) = 0.950; (*t)(1,19) = -0.028; (*t)(2,19) = 0.028;
3132 (*t)(0,20) = 1.000; (*t)(1,20) = -0.028; (*t)(2,20) = 0.028;
3133 (*t)(0,21) = 1.100; (*t)(1,21) = -0.028; (*t)(2,21) = 0.028;
3134 (*t)(0,22) = 1.200; (*t)(1,22) = -0.026; (*t)(2,22) = 0.028;
3135 (*t)(0,23) = 1.300; (*t)(1,23) = -0.024; (*t)(2,23) = 0.028;
3136 (*t)(0,24) = 1.400; (*t)(1,24) = -0.020; (*t)(2,24) = 0.028;
3137 (*t)(0,25) = 1.500; (*t)(1,25) = -0.018; (*t)(2,25) = 0.028;
3138 (*t)(0,26) = 1.600; (*t)(1,26) = -0.016; (*t)(2,26) = 0.028;
3139 (*t)(0,27) = 1.700; (*t)(1,27) = -0.014; (*t)(2,27) = 0.028;
3140 (*t)(0,28) = 1.800; (*t)(1,28) = -0.012; (*t)(2,28) = 0.026;
3141 (*t)(0,29) = 1.900; (*t)(1,29) = -0.010; (*t)(2,29) = 0.026;
3142 (*t)(0,30) = 2.000; (*t)(1,30) = -0.008; (*t)(2,30) = 0.026;
3143 (*t)(0,31) = 2.100; (*t)(1,31) = -0.008; (*t)(2,31) = 0.024;
3144 (*t)(0,32) = 2.200; (*t)(1,32) = -0.006; (*t)(2,32) = 0.024;
3145 (*t)(0,33) = 2.300; (*t)(1,33) = -0.004; (*t)(2,33) = 0.024;
3146 (*t)(0,34) = 2.400; (*t)(1,34) = -0.004; (*t)(2,34) = 0.024;
3147 (*t)(0,35) = 2.500; (*t)(1,35) = -0.002; (*t)(2,35) = 0.024;
3148 (*t)(0,36) = 2.600; (*t)(1,36) = -0.002; (*t)(2,36) = 0.024;
3149 (*t)(0,37) = 2.700; (*t)(1,37) = 0.000; (*t)(2,37) = 0.024;
3150 (*t)(0,38) = 2.800; (*t)(1,38) = 0.000; (*t)(2,38) = 0.026;
3151 (*t)(0,39) = 2.900; (*t)(1,39) = 0.000; (*t)(2,39) = 0.024;
3152 (*t)(0,40) = 3.000; (*t)(1,40) = 0.002; (*t)(2,40) = 0.026;
3153 (*t)(0,41) = 3.100; (*t)(1,41) = 0.002; (*t)(2,41) = 0.026;
3154 (*t)(0,42) = 3.200; (*t)(1,42) = 0.002; (*t)(2,42) = 0.026;
3155 (*t)(0,43) = 3.300; (*t)(1,43) = 0.002; (*t)(2,43) = 0.026;
3156 (*t)(0,44) = 3.400; (*t)(1,44) = 0.002; (*t)(2,44) = 0.026;
3157 (*t)(0,45) = 3.500; (*t)(1,45) = 0.002; (*t)(2,45) = 0.026;
3158 (*t)(0,46) = 3.600; (*t)(1,46) = 0.002; (*t)(2,46) = 0.026;
3159 (*t)(0,47) = 3.700; (*t)(1,47) = 0.002; (*t)(2,47) = 0.026;
3160 (*t)(0,48) = 3.800; (*t)(1,48) = 0.002; (*t)(2,48) = 0.026;
3161 (*t)(0,49) = 3.900; (*t)(1,49) = 0.004; (*t)(2,49) = 0.024;
3162 (*t)(0,50) = 4.000; (*t)(1,50) = 0.004; (*t)(2,50) = 0.026;
3163 (*t)(0,51) = 4.100; (*t)(1,51) = 0.004; (*t)(2,51) = 0.026;
3164 (*t)(0,52) = 4.200; (*t)(1,52) = 0.004; (*t)(2,52) = 0.024;
3165 (*t)(0,53) = 4.300; (*t)(1,53) = 0.006; (*t)(2,53) = 0.024;
3166 (*t)(0,54) = 4.400; (*t)(1,54) = 0.000; (*t)(2,54) = 0.000;
3167 (*t)(0,55) = 4.500; (*t)(1,55) = 0.000; (*t)(2,55) = 0.000;
3168 (*t)(0,56) = 4.600; (*t)(1,56) = 0.000; (*t)(2,56) = 0.000;
3169 (*t)(0,57) = 4.700; (*t)(1,57) = 0.000; (*t)(2,57) = 0.000;
3170 (*t)(0,58) = 4.800; (*t)(1,58) = 0.000; (*t)(2,58) = 0.000;
3171 (*t)(0,59) = 4.900; (*t)(1,59) = 0.000; (*t)(2,59) = 0.000;
3172 (*t)(0,60) = 5.900; (*t)(1,60) = 0.000; (*t)(2,60) = 0.000;
32b846cd 3173 }
3174 else
2b1eaa10 3175 if (fParticleID==AliPID::kProton)
32b846cd 3176 {
a0241c3a 3177 //TOF protons, 0.9 purity
3178 t = new TMatrixF(3,61);
d79d61fb 3179 (*t)(0,0) = 0.000; (*t)(1,0) = 0.000; (*t)(2,0) = 0.000;
3180 (*t)(0,1) = 0.050; (*t)(1,1) = 0.000; (*t)(2,1) = 0.000;
3181 (*t)(0,2) = 0.100; (*t)(1,2) = 0.000; (*t)(2,2) = 0.000;
3182 (*t)(0,3) = 0.150; (*t)(1,3) = 0.000; (*t)(2,3) = 0.000;
1a80f9f6 3183 (*t)(0,4) = 0.200; (*t)(1,4) = -0.07; (*t)(2,4) = 0.07;
3184 (*t)(0,5) = 0.200; (*t)(1,5) = -0.07; (*t)(2,5) = 0.07;
3185 (*t)(0,6) = 0.200; (*t)(1,6) = -0.07; (*t)(2,6) = 0.07;
3186 (*t)(0,7) = 0.200; (*t)(1,7) = -0.07; (*t)(2,7) = 0.07;
3187 (*t)(0,8) = 0.200; (*t)(1,8) = -0.07; (*t)(2,8) = 0.07;
3188 (*t)(0,9) = 0.200; (*t)(1,9) = -0.07; (*t)(2,9) = 0.07;
3189 (*t)(0,10) = 0.200; (*t)(1,10) = -0.07; (*t)(2,10) = 0.07;
3190 (*t)(0,11) = 0.200; (*t)(1,11) = -0.07; (*t)(2,11) = 0.07;
3191 (*t)(0,12) = 0.200; (*t)(1,12) = -0.07; (*t)(2,12) = 0.07;
3192 (*t)(0,13) = 0.200; (*t)(1,13) = -0.07; (*t)(2,13) = 0.07;
3193 (*t)(0,14) = 0.200; (*t)(1,14) = -0.07; (*t)(2,14) = 0.07;
3194 (*t)(0,15) = 0.200; (*t)(1,15) = -0.07; (*t)(2,15) = 0.07;
3195 (*t)(0,16) = 0.200; (*t)(1,16) = -0.07; (*t)(2,16) = 0.07;
d79d61fb 3196 (*t)(0,17) = 0.850; (*t)(1,17) = -0.070; (*t)(2,17) = 0.070;
3197 (*t)(0,18) = 0.900; (*t)(1,18) = -0.072; (*t)(2,18) = 0.072;
3198 (*t)(0,19) = 0.950; (*t)(1,19) = -0.072; (*t)(2,19) = 0.072;
3199 (*t)(0,20) = 1.000; (*t)(1,20) = -0.074; (*t)(2,20) = 0.074;
3200 (*t)(0,21) = 1.100; (*t)(1,21) = -0.032; (*t)(2,21) = 0.032;
3201 (*t)(0,22) = 1.200; (*t)(1,22) = -0.026; (*t)(2,22) = 0.026;
3202 (*t)(0,23) = 1.300; (*t)(1,23) = -0.026; (*t)(2,23) = 0.026;
3203 (*t)(0,24) = 1.400; (*t)(1,24) = -0.024; (*t)(2,24) = 0.024;
3204 (*t)(0,25) = 1.500; (*t)(1,25) = -0.024; (*t)(2,25) = 0.024;
3205 (*t)(0,26) = 1.600; (*t)(1,26) = -0.026; (*t)(2,26) = 0.026;
3206 (*t)(0,27) = 1.700; (*t)(1,27) = -0.026; (*t)(2,27) = 0.026;
3207 (*t)(0,28) = 1.800; (*t)(1,28) = -0.026; (*t)(2,28) = 0.026;
3208 (*t)(0,29) = 1.900; (*t)(1,29) = -0.026; (*t)(2,29) = 0.026;
3209 (*t)(0,30) = 2.000; (*t)(1,30) = -0.026; (*t)(2,30) = 0.026;
3210 (*t)(0,31) = 2.100; (*t)(1,31) = -0.026; (*t)(2,31) = 0.026;
3211 (*t)(0,32) = 2.200; (*t)(1,32) = -0.026; (*t)(2,32) = 0.024;
3212 (*t)(0,33) = 2.300; (*t)(1,33) = -0.028; (*t)(2,33) = 0.022;
3213 (*t)(0,34) = 2.400; (*t)(1,34) = -0.028; (*t)(2,34) = 0.020;
3214 (*t)(0,35) = 2.500; (*t)(1,35) = -0.028; (*t)(2,35) = 0.018;
3215 (*t)(0,36) = 2.600; (*t)(1,36) = -0.028; (*t)(2,36) = 0.016;
3216 (*t)(0,37) = 2.700; (*t)(1,37) = -0.028; (*t)(2,37) = 0.016;
3217 (*t)(0,38) = 2.800; (*t)(1,38) = -0.030; (*t)(2,38) = 0.014;
3218 (*t)(0,39) = 2.900; (*t)(1,39) = -0.030; (*t)(2,39) = 0.012;
3219 (*t)(0,40) = 3.000; (*t)(1,40) = -0.030; (*t)(2,40) = 0.012;
3220 (*t)(0,41) = 3.100; (*t)(1,41) = -0.030; (*t)(2,41) = 0.010;
3221 (*t)(0,42) = 3.200; (*t)(1,42) = -0.030; (*t)(2,42) = 0.010;
3222 (*t)(0,43) = 3.300; (*t)(1,43) = -0.030; (*t)(2,43) = 0.010;
3223 (*t)(0,44) = 3.400; (*t)(1,44) = -0.030; (*t)(2,44) = 0.008;
3224 (*t)(0,45) = 3.500; (*t)(1,45) = -0.030; (*t)(2,45) = 0.008;
3225 (*t)(0,46) = 3.600; (*t)(1,46) = -0.030; (*t)(2,46) = 0.008;
3226 (*t)(0,47) = 3.700; (*t)(1,47) = -0.030; (*t)(2,47) = 0.006;
3227 (*t)(0,48) = 3.800; (*t)(1,48) = -0.030; (*t)(2,48) = 0.006;
3228 (*t)(0,49) = 3.900; (*t)(1,49) = -0.030; (*t)(2,49) = 0.006;
3229 (*t)(0,50) = 4.000; (*t)(1,50) = -0.028; (*t)(2,50) = 0.004;
3230 (*t)(0,51) = 4.100; (*t)(1,51) = -0.030; (*t)(2,51) = 0.004;
3231 (*t)(0,52) = 4.200; (*t)(1,52) = -0.030; (*t)(2,52) = 0.004;
3232 (*t)(0,53) = 4.300; (*t)(1,53) = -0.028; (*t)(2,53) = 0.002;
3233 (*t)(0,54) = 4.400; (*t)(1,54) = -0.030; (*t)(2,54) = 0.002;
3234 (*t)(0,55) = 4.500; (*t)(1,55) = -0.028; (*t)(2,55) = 0.002;
3235 (*t)(0,56) = 4.600; (*t)(1,56) = -0.028; (*t)(2,56) = 0.002;
3236 (*t)(0,57) = 4.700; (*t)(1,57) = -0.028; (*t)(2,57) = 0.000;
3237 (*t)(0,58) = 4.800; (*t)(1,58) = -0.028; (*t)(2,58) = 0.002;
3238 (*t)(0,59) = 4.900; (*t)(1,59) = 0.000; (*t)(2,59) = 0.000;
3239 (*t)(0,60) = 5.900; (*t)(1,60) = 0.000; (*t)(2,60) = 0.000;
32b846cd 3240 }
3241 else
2b1eaa10 3242 if (fParticleID==AliPID::kKaon)
32b846cd 3243 {
a0241c3a 3244 //TOF kaons, 0.9 purity
3245 t = new TMatrixF(3,61);
d79d61fb 3246 (*t)(0,0) = 0.000; (*t)(1,0) = 0.000; (*t)(2,0) = 0.000;
3247 (*t)(0,1) = 0.050; (*t)(1,1) = 0.000; (*t)(2,1) = 0.000;
3248 (*t)(0,2) = 0.100; (*t)(1,2) = 0.000; (*t)(2,2) = 0.000;
3249 (*t)(0,3) = 0.150; (*t)(1,3) = 0.000; (*t)(2,3) = 0.000;
1a80f9f6 3250 (*t)(0,4) = 0.200; (*t)(1,4) = -0.05; (*t)(2,4) = 0.05;
3251 (*t)(0,5) = 0.200; (*t)(1,5) = -0.05; (*t)(2,5) = 0.05;
3252 (*t)(0,6) = 0.200; (*t)(1,6) = -0.05; (*t)(2,6) = 0.05;
3253 (*t)(0,7) = 0.200; (*t)(1,7) = -0.05; (*t)(2,7) = 0.05;
3254 (*t)(0,8) = 0.200; (*t)(1,8) = -0.05; (*t)(2,8) = 0.05;
3255 (*t)(0,9) = 0.200; (*t)(1,9) = -0.05; (*t)(2,9) = 0.05;
3256 (*t)(0,10) = 0.200; (*t)(1,10) = -0.05; (*t)(2,10) = 0.05;
d79d61fb 3257 (*t)(0,11) = 0.550; (*t)(1,11) = -0.026; (*t)(2,11) = 0.026;
3258 (*t)(0,12) = 0.600; (*t)(1,12) = -0.026; (*t)(2,12) = 0.026;
3259 (*t)(0,13) = 0.650; (*t)(1,13) = -0.026; (*t)(2,13) = 0.026;
3260 (*t)(0,14) = 0.700; (*t)(1,14) = -0.026; (*t)(2,14) = 0.026;
3261 (*t)(0,15) = 0.750; (*t)(1,15) = -0.026; (*t)(2,15) = 0.026;
3262 (*t)(0,16) = 0.800; (*t)(1,16) = -0.026; (*t)(2,16) = 0.026;
3263 (*t)(0,17) = 0.850; (*t)(1,17) = -0.024; (*t)(2,17) = 0.024;
3264 (*t)(0,18) = 0.900; (*t)(1,18) = -0.024; (*t)(2,18) = 0.024;
3265 (*t)(0,19) = 0.950; (*t)(1,19) = -0.024; (*t)(2,19) = 0.024;
3266 (*t)(0,20) = 1.000; (*t)(1,20) = -0.024; (*t)(2,20) = 0.024;
3267 (*t)(0,21) = 1.100; (*t)(1,21) = -0.024; (*t)(2,21) = 0.024;
3268 (*t)(0,22) = 1.200; (*t)(1,22) = -0.024; (*t)(2,22) = 0.022;
3269 (*t)(0,23) = 1.300; (*t)(1,23) = -0.024; (*t)(2,23) = 0.020;
3270 (*t)(0,24) = 1.400; (*t)(1,24) = -0.026; (*t)(2,24) = 0.016;
3271 (*t)(0,25) = 1.500; (*t)(1,25) = -0.028; (*t)(2,25) = 0.014;
3272 (*t)(0,26) = 1.600; (*t)(1,26) = -0.028; (*t)(2,26) = 0.012;
3273 (*t)(0,27) = 1.700; (*t)(1,27) = -0.028; (*t)(2,27) = 0.010;
3274 (*t)(0,28) = 1.800; (*t)(1,28) = -0.028; (*t)(2,28) = 0.010;
3275 (*t)(0,29) = 1.900; (*t)(1,29) = -0.028; (*t)(2,29) = 0.008;
3276 (*t)(0,30) = 2.000; (*t)(1,30) = -0.028; (*t)(2,30) = 0.006;
3277 (*t)(0,31) = 2.100; (*t)(1,31) = -0.026; (*t)(2,31) = 0.006;
3278 (*t)(0,32) = 2.200; (*t)(1,32) = -0.024; (*t)(2,32) = 0.004;
3279 (*t)(0,33) = 2.300; (*t)(1,33) = -0.020; (*t)(2,33) = 0.002;
3280 (*t)(0,34) = 2.400; (*t)(1,34) = -0.020; (*t)(2,34) = 0.002;
3281 (*t)(0,35) = 2.500; (*t)(1,35) = -0.018; (*t)(2,35) = 0.000;
3282 (*t)(0,36) = 2.600; (*t)(1,36) = -0.016; (*t)(2,36) = 0.000;
3283 (*t)(0,37) = 2.700; (*t)(1,37) = -0.014; (*t)(2,37) = -0.002;
3284 (*t)(0,38) = 2.800; (*t)(1,38) = -0.014; (*t)(2,38) = -0.004;
3285 (*t)(0,39) = 2.900; (*t)(1,39) = -0.012; (*t)(2,39) = -0.004;
3286 (*t)(0,40) = 3.000; (*t)(1,40) = -0.010; (*t)(2,40) = -0.006;
3287 (*t)(0,41) = 3.100; (*t)(1,41) = 0.000; (*t)(2,41) = 0.000;
3288 (*t)(0,42) = 3.200; (*t)(1,42) = 0.000; (*t)(2,42) = 0.000;
3289 (*t)(0,43) = 3.300; (*t)(1,43) = 0.000; (*t)(2,43) = 0.000;
3290 (*t)(0,44) = 3.400; (*t)(1,44) = 0.000; (*t)(2,44) = 0.000;
3291 (*t)(0,45) = 3.500; (*t)(1,45) = 0.000; (*t)(2,45) = 0.000;
3292 (*t)(0,46) = 3.600; (*t)(1,46) = 0.000; (*t)(2,46) = 0.000;
3293 (*t)(0,47) = 3.700; (*t)(1,47) = 0.000; (*t)(2,47) = 0.000;
3294 (*t)(0,48) = 3.800; (*t)(1,48) = 0.000; (*t)(2,48) = 0.000;
3295 (*t)(0,49) = 3.900; (*t)(1,49) = 0.000; (*t)(2,49) = 0.000;
3296 (*t)(0,50) = 4.000; (*t)(1,50) = 0.000; (*t)(2,50) = 0.000;
3297 (*t)(0,51) = 4.100; (*t)(1,51) = 0.000; (*t)(2,51) = 0.000;
3298 (*t)(0,52) = 4.200; (*t)(1,52) = 0.000; (*t)(2,52) = 0.000;
3299 (*t)(0,53) = 4.300; (*t)(1,53) = 0.000; (*t)(2,53) = 0.000;
3300 (*t)(0,54) = 4.400; (*t)(1,54) = 0.000; (*t)(2,54) = 0.000;
3301 (*t)(0,55) = 4.500; (*t)(1,55) = 0.000; (*t)(2,55) = 0.000;
3302 (*t)(0,56) = 4.600; (*t)(1,56) = 0.000; (*t)(2,56) = 0.000;
3303 (*t)(0,57) = 4.700; (*t)(1,57) = 0.000; (*t)(2,57) = 0.000;
3304 (*t)(0,58) = 4.800; (*t)(1,58) = 0.000; (*t)(2,58) = 0.000;
3305 (*t)(0,59) = 4.900; (*t)(1,59) = 0.000; (*t)(2,59) = 0.000;
3306 (*t)(0,60) = 5.900; (*t)(1,60) = 0.000; (*t)(2,60) = 0.000;
32b846cd 3307 }
d79d61fb 3308 delete fTOFpidCuts;
32b846cd 3309 fTOFpidCuts=t;
3310 }
3311}
3abf7ecc 3312
0d16f553 3313//-----------------------------------------------------------------------
3314//-----------------------------------------------------------------------
3315Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliAODTrack* track)
3316{
3317 fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3318 Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3319
3320 Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
3321
3322 if(! kTPC) return kFALSE;
3323
3324 fProbBayes = 0.0;
3325
3326switch (fParticleID)
3327 {
3328 case AliPID::kPion:
3329 fProbBayes = probabilities[2];
3330 break;
3331 case AliPID::kKaon:
3332 fProbBayes = probabilities[3];
3333 break;
3334 case AliPID::kProton:
3335 fProbBayes = probabilities[4];
3336 break;
3337 case AliPID::kElectron:
3338 fProbBayes = probabilities[0];
3339 break;
3340 case AliPID::kMuon:
3341 fProbBayes = probabilities[1];
3342 break;
3343 case AliPID::kDeuteron:
3344 fProbBayes = probabilities[5];
3345 break;
3346 case AliPID::kTriton:
3347 fProbBayes = probabilities[6];
3348 break;
3349 case AliPID::kHe3:
3350 fProbBayes = probabilities[7];
3351 break;
3352 default:
3353 return kFALSE;
3354 }
3355
3356 if(fProbBayes > fParticleProbability){
3357 if(!fCutCharge)
3358 return kTRUE;
3359 else if (fCutCharge && fCharge * track->Charge() > 0)
3360 return kTRUE;
3361 }
3362 return kFALSE;
3363}
cd62a2a8 3364//-----------------------------------------------------------------------
3365Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliESDtrack* track)
3366{
3367 //cut on TPC bayesian pid
cd62a2a8 3368 //if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3369
3370 //Bool_t statusMatchingHard = TPCTOFagree(track);
3371 //if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3372 // return kFALSE;
8b649d8c 3373 fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
9126fd0c 3374 Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
0d16f553 3375 //Int_t kTOF = fBayesianResponse->GetCurrentMask(1); // is TOF on
cd62a2a8 3376
9126fd0c 3377 if(! kTPC) return kFALSE;
3378
0d16f553 3379 // Bool_t statusMatchingHard = 1;
3380 // Float_t mismProb = 0;
3381 // if(kTOF){
3382 // statusMatchingHard = TPCTOFagree(track);
3383 // mismProb = fBayesianResponse->GetTOFMismProb();
3384 // }
3385 // if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3386 // return kFALSE;
9126fd0c 3387
9126fd0c 3388 Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3389
8b649d8c 3390 fProbBayes = 0.0;
cd62a2a8 3391
cd62a2a8 3392 switch (fParticleID)
3393 {
3394 case AliPID::kPion:
8b649d8c 3395 fProbBayes = probabilities[2];
cd62a2a8 3396 break;
3397 case AliPID::kKaon:
8b649d8c 3398 fProbBayes = probabilities[3];
cd62a2a8 3399 break;
3400 case AliPID::kProton:
8b649d8c 3401 fProbBayes = probabilities[4];
cd62a2a8 3402 break;
3403 case AliPID::kElectron:
8b649d8c 3404 fProbBayes = probabilities[0];
3405 break;
3406 case AliPID::kMuon:
3407 fProbBayes = probabilities[1];
3408 break;
3409 case AliPID::kDeuteron:
3410 fProbBayes = probabilities[5];
3411 break;
3412 case AliPID::kTriton:
3413 fProbBayes = probabilities[6];
3414 break;
3415 case AliPID::kHe3:
3416 fProbBayes = probabilities[7];
cd62a2a8 3417 break;
3418 default:
3419 return kFALSE;
3420 }
3421
0d16f553 3422 if(fProbBayes > fParticleProbability)
9126fd0c 3423 {
3424 if(!fCutCharge)
3425 return kTRUE;
3426 else if (fCutCharge && fCharge * track->GetSign() > 0)
3427 return kTRUE;
3428 }
cd62a2a8 3429 return kFALSE;
3430}
0d16f553 3431//-----------------------------------------------------------------------
3432Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliAODTrack* track)
3433{
a22c1eb1 3434 //check is track passes bayesian combined TOF+TPC pid cut
0d16f553 3435 Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3436 (track->GetStatus() & AliESDtrack::kTIME) &&
3437 (track->GetTOFsignal() > 12000) &&
3438 (track->GetTOFsignal() < 100000);
cd62a2a8 3439
0d16f553 3440 if (! goodtrack)
3441 return kFALSE;
3442
3443 Bool_t statusMatchingHard = TPCTOFagree(track);
0d16f553 3444 if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3445 return kFALSE;
3446
3447 fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3448 Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3449
3450 Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3451
3452 fProbBayes = 0.0;
3453
3454 switch (fParticleID)
3455 {
3456 case AliPID::kPion:
3457 fProbBayes = probabilities[2];
3458 break;
3459 case AliPID::kKaon:
3460 fProbBayes = probabilities[3];
3461 break;
3462 case AliPID::kProton:
3463 fProbBayes = probabilities[4];
3464 break;
3465 case AliPID::kElectron:
3466 fProbBayes = probabilities[0];
3467 break;
3468 case AliPID::kMuon:
3469 fProbBayes = probabilities[1];
3470 break;
3471 case AliPID::kDeuteron:
3472 fProbBayes = probabilities[5];
3473 break;
3474 case AliPID::kTriton:
3475 fProbBayes = probabilities[6];
3476 break;
3477 case AliPID::kHe3:
3478 fProbBayes = probabilities[7];
3479 break;
3480 default:
3481 return kFALSE;
3482 }
3483
3484 if(fProbBayes > fParticleProbability && mismProb < 0.5){
3485 if(!fCutCharge)
3486 return kTRUE;
3487 else if (fCutCharge && fCharge * track->Charge() > 0)
3488 return kTRUE;
3489 }
3490 return kFALSE;
3491
3492}
3abf7ecc 3493//-----------------------------------------------------------------------
3494// part added by F. Noferini (some methods)
cd62a2a8 3495Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliESDtrack* track)
875ee12a 3496{
3c67c769 3497 //check is track passes bayesian combined TOF+TPC pid cut
8b649d8c 3498 Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3499 (track->GetStatus() & AliESDtrack::kTIME) &&
875ee12a 3500 (track->GetTOFsignal() > 12000) &&
3501 (track->GetTOFsignal() < 100000) &&
3502 (track->GetIntegratedLength() > 365);
3abf7ecc 3503
3abf7ecc 3504 if (! goodtrack)
3505 return kFALSE;
3506
6a043318 3507 if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3508
875ee12a 3509 Bool_t statusMatchingHard = TPCTOFagree(track);
3510 if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3511 return kFALSE;
3512
9126fd0c 3513 fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
3514 Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3515
9126fd0c 3516 Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3abf7ecc 3517
8b649d8c 3518 fProbBayes = 0.0;
3519
2b1eaa10 3520 switch (fParticleID)
3abf7ecc 3521 {
3522 case AliPID::kPion:
8b649d8c 3523 fProbBayes = probabilities[2];
3abf7ecc 3524 break;
3525 case AliPID::kKaon:
8b649d8c 3526 fProbBayes = probabilities[3];
3abf7ecc 3527 break;
3528 case AliPID::kProton:
8b649d8c 3529 fProbBayes = probabilities[4];
3abf7ecc 3530 break;
3531 case AliPID::kElectron:
8b649d8c 3532 fProbBayes = probabilities[0];
3abf7ecc 3533 break;
8b649d8c 3534 case AliPID::kMuon:
3535 fProbBayes = probabilities[1];
3536 break;
3537 case AliPID::kDeuteron:
3538 fProbBayes = probabilities[5];
3539 break;
3540 case AliPID::kTriton:
3541 fProbBayes = probabilities[6];
3542 break;
3543 case AliPID::kHe3:
3544 fProbBayes = probabilities[7];
3545 break;
3546 default:
3abf7ecc 3547 return kFALSE;
3548 }
3549
3550 // 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);
8b649d8c 3551 if(fProbBayes > fParticleProbability && mismProb < 0.5){
3abf7ecc 3552 if(!fCutCharge)
3553 return kTRUE;
3554 else if (fCutCharge && fCharge * track->GetSign() > 0)
3555 return kTRUE;
3556 }
3557 return kFALSE;
3558}
875ee12a 3559
8eca5d19 3560
3561//-----------------------------------------------------------------------
3562 // part added by Natasha
3563Bool_t AliFlowTrackCuts::PassesNucleiSelection(const AliESDtrack* track)
3564{
3c67c769 3565 //pid selection for heavy nuclei
8eca5d19 3566 Bool_t select=kFALSE;
3567
3568 //if (!track) continue;
6887199d 3569
8eca5d19 3570 if (!track->GetInnerParam())
3571 return kFALSE; //break;
6887199d 3572
8eca5d19 3573 const AliExternalTrackParam* tpcTrack = track->GetInnerParam();
6887199d 3574
8eca5d19 3575 Double_t ptotTPC = tpcTrack->GetP();
3576 Double_t sigTPC = track->GetTPCsignal();
3577 Double_t dEdxBBA = 0.;
3578 Double_t dSigma = 0.;
3579
3580 switch (fParticleID)
6887199d 3581 {
8eca5d19 3582 case AliPID::kDeuteron:
3583 //pid=10;
3584 dEdxBBA = AliExternalTrackParam::BetheBlochAleph(ptotTPC/1.8756,
6887199d 3585 4.60e+00,
3586 8.9684e+00,
3587 1.640e-05,
3588 2.35e+00,
3589 2.35e+00);
8eca5d19 3590 dSigma = (sigTPC - dEdxBBA)/dEdxBBA;
6887199d 3591
8eca5d19 3592 if( ptotTPC<=1.1 && (dSigma < (0.5 - (0.1818*ptotTPC)) ) && (dSigma > ( (0.218*ptotTPC - 0.4) ) ) )
6887199d 3593 {select=kTRUE;}
8eca5d19 3594 break;
6887199d 3595
8eca5d19 3596 case AliPID::kTriton:
3597 //pid=11;
3598 select=kFALSE;
3599 break;
6887199d 3600
8eca5d19 3601 case AliPID::kHe3:
3602 //pid=12;
6887199d 3603 // ----- Pass 2 -------
3604 dEdxBBA = 4.0 * AliExternalTrackParam::BetheBlochAleph( (2.*ptotTPC)/2.8084,
3605 1.74962,
3606 27.4992,
3607 4.00313e-15,
3608 2.42485,
3609 8.31768);
3610 dSigma = (sigTPC - dEdxBBA)/dEdxBBA;
3611 if(ptotTPC<=5.0 && (dSigma >= (-0.03968*ptotTPC - 0.1)) && (dSigma <= (0.31 - 0.0217*ptotTPC)))
3612 {select=kTRUE;}
8eca5d19 3613 break;
6887199d 3614
8eca5d19 3615 case AliPID::kAlpha:
3616 //pid=13;
3617 select=kFALSE;
3618 break;
6887199d 3619
8eca5d19 3620 default:
3621 return kFALSE;
6887199d 3622 }
3623
8eca5d19 3624 return select;
3625}
3626// end part added by Natasha
d0ce6147 3627//-----------------------------------------------------------------------
3628Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCut(const AliAODTrack* track)
3629{
3630 // do a simple combined cut on the n sigma from tpc and tof
3631 // with information of the pid response object (needs pid response task)
3632 // stub, not implemented yet
043e9010 3633 if(!fPIDResponse) return kFALSE;
d0ce6147 3634 if(!track) return kFALSE;
043e9010 3635
3636 // check TOF status
3637 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kFALSE;
3638 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kFALSE;
3639
3640 // check TPC status
3641 if(track->GetTPCsignal() < 10) return kFALSE;
3642
3643 Float_t nsigmaTPC = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC,track,fParticleID);
3644 Float_t nsigmaTOF = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF,track,fParticleID);
3645
3646 Float_t nsigma2 = nsigmaTPC*nsigmaTPC + nsigmaTOF*nsigmaTOF;
3647
3648 return (nsigma2 < fNsigmaCut2);
d0ce6147 3649
3650}
3651//-----------------------------------------------------------------------------
3652Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCut(const AliESDtrack* track)
3653{
3654 // do a simple combined cut on the n sigma from tpc and tof
3655 // with information of the pid response object (needs pid response task)
3656 // stub, not implemented yet
043e9010 3657 if(!fPIDResponse) return kFALSE;
d0ce6147 3658 if(!track) return kFALSE;
043e9010 3659
3660 // check TOF status
3661 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kFALSE;
3662 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kFALSE;
3663
3664 // check TPC status
3665 if(track->GetTPCsignal() < 10) return kFALSE;
3666
3667 Float_t nsigmaTPC = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC,track,fParticleID);
3668 Float_t nsigmaTOF = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF,track,fParticleID);
3669
3670 Float_t nsigma2 = nsigmaTPC*nsigmaTPC + nsigmaTOF*nsigmaTOF;
3671
3672 return (nsigma2 < fNsigmaCut2);
d0ce6147 3673
3674}
67006c99 3675//-----------------------------------------------------------------------------
3676Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCutPuritybased(const AliAODTrack* track)
3677{
3678 // do a combined cut on the n sigma from tpc and tof based on purity of the identified particle. (Particle is selected if Purity>0.8)
3679 // with information of the pid response object (needs pid response task)
d0ce6147 3680
67006c99 3681 if(!fPIDResponse) return kFALSE;
3682 if(!track) return kFALSE;
3683 if(!fUseTPCTOFNsigmaCutContours) return kFALSE;
3684
3685
3686 Bool_t pWithinRange = kFALSE;
3687 Int_t p_bin = -999;
3688 Double_t pBins[50];
3689 for(int b=0;b<50;b++){pBins[b] = 0.1*b;}
c2a07d6a 3690 for(int i=0;i<49;i++){ // fixed from <50 14012015 to avoid out-of-bounds RAB
67006c99 3691 if(track->P()>pBins[i] && track->P()<(pBins[i+1])){
3692 p_bin = i;
3693 }
3694 }
3695 Int_t ispecie = 0;
3696 switch (fParticleID)
3697 {
3698 case AliPID::kPion:
3699 ispecie=0;
3700 break;
3701 case AliPID::kKaon:
3702 ispecie=1;
3703 break;
3704 case AliPID::kProton:
3705 ispecie=2;
3706 break;
3707 default:
3708 return kFALSE;
3709 }
3710
3711 Double_t LowPtPIDTPCnsigLow_Pion[6] = {-3,-3,-3,-3,-3,-3};
3712 Double_t LowPtPIDTPCnsigLow_Kaon[6] = {-3,-2,0,-1.8,-1.2,-0.8}; //for0.4<Pt<0.5 the purity is lower than 0.7
3713 Double_t LowPtPIDTPCnsigHigh_Pion[6] ={2.4,3,3,3,2,1.4};
3714 Double_t LowPtPIDTPCnsigHigh_Kaon[6] ={3,2.2,0,-0.2,1,1.8}; //for 0.4<Pt<0.5 he purity is lower than 0.7
3715
3716 Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
3717 Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
3718
3719
3720 if ( (track->IsOn(AliAODTrack::kITSin)) && (track->IsOn(AliAODTrack::kTOFpid))){
3721 if(p_bin<0) return kFALSE;
3722 if(!fCutContour[p_bin]) cout<<"fCutContour[p_bin] does not exist"<<endl;
3723
3724 if(p_bin==4 || p_bin>7){
3725 if(fCutContour[p_bin]->IsInside(nSigmaTOF,nSigmaTPC)){
3726 return kTRUE;
3727 }
3728 else{
3729 return kFALSE;
3730 }
3731 }
3732
3733 if(p_bin<8 && p_bin!=4){
3734 if(fParticleID==AliPID::kPion && nSigmaTPC>LowPtPIDTPCnsigLow_Pion[p_bin-2] && nSigmaTPC<LowPtPIDTPCnsigHigh_Pion[p_bin-2]){
3735 return kTRUE;
3736 }
3737 if(fParticleID==AliPID::kKaon && nSigmaTPC>LowPtPIDTPCnsigLow_Kaon[p_bin-2] && nSigmaTPC<LowPtPIDTPCnsigHigh_Kaon[p_bin-2]){return kTRUE;}
3738
3739 if(fParticleID==AliPID::kProton && nSigmaTPC>-3 && nSigmaTPC<3){return kTRUE;}
3740 }
3741 }
3742
3743 return kFALSE;
3744}
3abf7ecc 3745//-----------------------------------------------------------------------
875ee12a 3746void AliFlowTrackCuts::SetPriors(Float_t centrCur){
cc0afcfc 3747 //set priors for the bayesian pid selection
9126fd0c 3748 fCurrCentr = centrCur;
3749
875ee12a 3750 fBinLimitPID[0] = 0.300000;
3751 fBinLimitPID[1] = 0.400000;
3752 fBinLimitPID[2] = 0.500000;
3753 fBinLimitPID[3] = 0.600000;
3754 fBinLimitPID[4] = 0.700000;
3755 fBinLimitPID[5] = 0.800000;
3756 fBinLimitPID[6] = 0.900000;
3757 fBinLimitPID[7] = 1.000000;
3758 fBinLimitPID[8] = 1.200000;
3759 fBinLimitPID[9] = 1.400000;
3760 fBinLimitPID[10] = 1.600000;
3761 fBinLimitPID[11] = 1.800000;
3762 fBinLimitPID[12] = 2.000000;
3763 fBinLimitPID[13] = 2.200000;
3764 fBinLimitPID[14] = 2.400000;
3765 fBinLimitPID[15] = 2.600000;
3766 fBinLimitPID[16] = 2.800000;
3767 fBinLimitPID[17] = 3.000000;
3768
3769 // 0-10%
3770 if(centrCur < 10){
3771 fC[0][0] = 0.005;
3772 fC[0][1] = 0.005;
3773 fC[0][2] = 1.0000;
3774 fC[0][3] = 0.0010;
3775 fC[0][4] = 0.0010;
3776
3777 fC[1][0] = 0.005;
3778 fC[1][1] = 0.005;
3779 fC[1][2] = 1.0000;
3780 fC[1][3] = 0.0168;
3781 fC[1][4] = 0.0122;
3782
3783 fC[2][0] = 0.005;
3784 fC[2][1] = 0.005;
3785 fC[2][2] = 1.0000;
3786 fC[2][3] = 0.0272;
3787 fC[2][4] = 0.0070;
3788
3789 fC[3][0] = 0.005;
3790 fC[3][1] = 0.005;
3791 fC[3][2] = 1.0000;
3792 fC[3][3] = 0.0562;
3793 fC[3][4] = 0.0258;
3794
3795 fC[4][0] = 0.005;
3796 fC[4][1] = 0.005;
3797 fC[4][2] = 1.0000;
3798 fC[4][3] = 0.0861;
3799 fC[4][4] = 0.0496;
3800
3801 fC[5][0] = 0.005;
3802 fC[5][1] = 0.005;
3803 fC[5][2] = 1.0000;
3804 fC[5][3] = 0.1168;
3805 fC[5][4] = 0.0740;
3806
3807 fC[6][0] = 0.005;
3808 fC[6][1] = 0.005;
3809 fC[6][2] = 1.0000;
3810 fC[6][3] = 0.1476;
3811 fC[6][4] = 0.0998;
3812
3813 fC[7][0] = 0.005;
3814 fC[7][1] = 0.005;
3815 fC[7][2] = 1.0000;
3816 fC[7][3] = 0.1810;
3817 fC[7][4] = 0.1296;
3818
3819 fC[8][0] = 0.005;
3820 fC[8][1] = 0.005;
3821 fC[8][2] = 1.0000;
3822 fC[8][3] = 0.2240;
3823 fC[8][4] = 0.1827;
3824
3825 fC[9][0] = 0.005;
3826 fC[9][1] = 0.005;
3827 fC[9][2] = 1.0000;
3828 fC[9][3] = 0.2812;
3829 fC[9][4] = 0.2699;
3830
3831 fC[10][0] = 0.005;
3832 fC[10][1] = 0.005;
3833 fC[10][2] = 1.0000;
3834 fC[10][3] = 0.3328;
3835 fC[10][4] = 0.3714;
3836
3837 fC[11][0] = 0.005;
3838 fC[11][1] = 0.005;
3839 fC[11][2] = 1.0000;
3840 fC[11][3] = 0.3780;
3841 fC[11][4] = 0.4810;
3842
3843 fC[12][0] = 0.005;
3844 fC[12][1] = 0.005;
3845 fC[12][2] = 1.0000;
3846 fC[12][3] = 0.4125;
3847 fC[12][4] = 0.5771;
3848
3849 fC[13][0] = 0.005;
3850 fC[13][1] = 0.005;
3851 fC[13][2] = 1.0000;
3852 fC[13][3] = 0.4486;
3853 fC[13][4] = 0.6799;
3854
3855 fC[14][0] = 0.005;
3856 fC[14][1] = 0.005;
3857 fC[14][2] = 1.0000;
3858 fC[14][3] = 0.4840;
3859 fC[14][4] = 0.7668;
3860
3861 fC[15][0] = 0.005;
3862 fC[15][1] = 0.005;
3863 fC[15][2] = 1.0000;
3864 fC[15][3] = 0.4971;
3865 fC[15][4] = 0.8288;
3866
3867 fC[16][0] = 0.005;
3868 fC[16][1] = 0.005;
3869 fC[16][2] = 1.0000;
3870 fC[16][3] = 0.4956;
3871 fC[16][4] = 0.8653;
3872
3873 fC[17][0] = 0.005;
3874 fC[17][1] = 0.005;
3875 fC[17][2] = 1.0000;
3876 fC[17][3] = 0.5173;
3877 fC[17][4] = 0.9059;
3878 }
3879 // 10-20%
3880 else if(centrCur < 20){
3881 fC[0][0] = 0.005;
3882 fC[0][1] = 0.005;
3883 fC[0][2] = 1.0000;
3884 fC[0][3] = 0.0010;
3885 fC[0][4] = 0.0010;
3886
3887 fC[1][0] = 0.005;
3888 fC[1][1] = 0.005;
3889 fC[1][2] = 1.0000;
3890 fC[1][3] = 0.0132;
3891 fC[1][4] = 0.0088;
3892
3893 fC[2][0] = 0.005;
3894 fC[2][1] = 0.005;
3895 fC[2][2] = 1.0000;
3896 fC[2][3] = 0.0283;
3897 fC[2][4] = 0.0068;
3898
3899 fC[3][0] = 0.005;
3900 fC[3][1] = 0.005;
3901 fC[3][2] = 1.0000;
3902 fC[3][3] = 0.0577;
3903 fC[3][4] = 0.0279;
3904
3905 fC[4][0] = 0.005;
3906 fC[4][1] = 0.005;
3907 fC[4][2] = 1.0000;
3908 fC[4][3] = 0.0884;
3909 fC[4][4] = 0.0534;
3910
3911 fC[5][0] = 0.005;
3912 fC[5][1] = 0.005;
3913 fC[5][2] = 1.0000;
3914 fC[5][3] = 0.1179;
3915 fC[5][4] = 0.0794;
3916
3917 fC[6][0] = 0.005;
3918 fC[6][1] = 0.005;
3919 fC[6][2] = 1.0000;
3920 fC[6][3] = 0.1480;
3921 fC[6][4] = 0.1058;
3922
3923 fC[7][0] = 0.005;
3924 fC[7][1] = 0.005;
3925 fC[7][2] = 1.0000;
3926 fC[7][3] = 0.1807;
3927 fC[7][4] = 0.1366;
3928
3929 fC[8][0] = 0.005;
3930 fC[8][1] = 0.005;
3931 fC[8][2] = 1.0000;
3932 fC[8][3] = 0.2219;
3933 fC[8][4] = 0.1891;
3934
3935 fC[9][0] = 0.005;
3936 fC[9][1] = 0.005;
3937 fC[9][2] = 1.0000;
3938 fC[9][3] = 0.2804;
3939 fC[9][4] = 0.2730;
3940
3941 fC[10][0] = 0.005;
3942 fC[10][1] = 0.005;
3943 fC[10][2] = 1.0000;
3944 fC[10][3] = 0.3283;
3945 fC[10][4] = 0.3660;
3946
3947 fC[11][0] = 0.005;
3948 fC[11][1] = 0.005;
3949 fC[11][2] = 1.0000;
3950 fC[11][3] = 0.3710;
3951 fC[11][4] = 0.4647;
3952
3953 fC[12][0] = 0.005;
3954 fC[12][1] = 0.005;
3955 fC[12][2] = 1.0000;
3956 fC[12][3] = 0.4093;
3957 fC[12][4] = 0.5566;
3958
3959 fC[13][0] = 0.005;
3960 fC[13][1] = 0.005;
3961 fC[13][2] = 1.0000;
3962 fC[13][3] = 0.4302;
3963 fC[13][4] = 0.6410;
3964
3965 fC[14][0] = 0.005;
3966 fC[14][1] = 0.005;
3967 fC[14][2] = 1.0000;
3968 fC[14][3] = 0.4649;
3969 fC[14][4] = 0.7055;
3970
3971 fC[15][0] = 0.005;
3972 fC[15][1] = 0.005;
3973 fC[15][2] = 1.0000;
3974 fC[15][3] = 0.4523;
3975 fC[15][4] = 0.7440;
3976
3977 fC[16][0] = 0.005;
3978 fC[16][1] = 0.005;
3979 fC[16][2] = 1.0000;
3980 fC[16][3] = 0.4591;
3981 fC[16][4] = 0.7799;
3982
3983 fC[17][0] = 0.005;
3984 fC[17][1] = 0.005;
3985 fC[17][2] = 1.0000;
3986 fC[17][3] = 0.4804;
3987 fC[17][4] = 0.8218;
3988 }
3989 // 20-30%
3990 else if(centrCur < 30){
3991 fC[0][0] = 0.005;
3992 fC[0][1] = 0.005;
3993 fC[0][2] = 1.0000;
3994 fC[0][3] = 0.0010;
3995 fC[0][4] = 0.0010;
3996
3997 fC[1][0] = 0.005;
3998 fC[1][1] = 0.005;
3999 fC[1][2] = 1.0000;
4000 fC[1][3] = 0.0102;
4001 fC[1][4] = 0.0064;
4002
4003 fC[2][0] = 0.005;
4004 fC[2][1] = 0.005;
4005 fC[2][2] = 1.0000;
4006 fC[2][3] = 0.0292;
4007 fC[2][4] = 0.0066;
4008
4009 fC[3][0] = 0.005;
4010 fC[3][1] = 0.005;
4011 fC[3][2] = 1.0000;
4012 fC[3][3] = 0.0597;
4013 fC[3][4] = 0.0296;
4014
4015 fC[4][0] = 0.005;
4016 fC[4][1] = 0.005;
4017 fC[4][2] = 1.0000;
4018 fC[4][3] = 0.0900;
4019 fC[4][4] = 0.0589;
4020
4021 fC[5][0] = 0.005;
4022 fC[5][1] = 0.005;
4023 fC[5][2] = 1.0000;
4024 fC[5][3] = 0.1199;
4025 fC[5][4] = 0.0859;
4026
4027 fC[6][0] = 0.005;
4028 fC[6][1] = 0.005;
4029 fC[6][2] = 1.0000;
4030 fC[6][3] = 0.1505;
4031 fC[6][4] = 0.1141;
4032
4033 fC[7][0] = 0.005;
4034 fC[7][1] = 0.005;
4035 fC[7][2] = 1.0000;
4036 fC[7][3] = 0.1805;
4037 fC[7][4] = 0.1454;
4038
4039 fC[8][0] = 0.005;
4040 fC[8][1] = 0.005;
4041 fC[8][2] = 1.0000;
4042 fC[8][3] = 0.2221;
4043 fC[8][4] = 0.2004;
4044
4045 fC[9][0] = 0.005;
4046 fC[9][1] = 0.005;
4047 fC[9][2] = 1.0000;
4048 fC[9][3] = 0.2796;
4049 fC[9][4] = 0.2838;
4050
4051 fC[10][0] = 0.005;
4052 fC[10][1] = 0.005;
4053 fC[10][2] = 1.0000;
4054 fC[10][3] = 0.3271;
4055 fC[10][4] = 0.3682;
4056
4057 fC[11][0] = 0.005;
4058 fC[11][1] = 0.005;
4059 fC[11][2] = 1.0000;
4060 fC[11][3] = 0.3648;
4061 fC[11][4] = 0.4509;
4062
4063 fC[12][0] = 0.005;
4064 fC[12][1] = 0.005;
4065 fC[12][2] = 1.0000;
4066 fC[12][3] = 0.3988;
4067 fC[12][4] = 0.5339;
4068
4069 fC[13][0] = 0.005;
4070 fC[13][1] = 0.005;
4071 fC[13][2] = 1.0000;
4072 fC[13][3] = 0.4315;
4073 fC[13][4] = 0.5995;
4074
4075 fC[14][0] = 0.005;
4076 fC[14][1] = 0.005;
4077 fC[14][2] = 1.0000;
4078 fC[14][3] = 0.4548;
4079 fC[14][4] = 0.6612;
4080
4081 fC[15][0] = 0.005;
4082 fC[15][1] = 0.005;
4083 fC[15][2] = 1.0000;
4084 fC[15][3] = 0.4744;
4085 fC[15][4] = 0.7060;
4086
4087 fC[16][0] = 0.005;
4088 fC[16][1] = 0.005;
4089 fC[16][2] = 1.0000;
4090 fC[16][3] = 0.4899;
4091 fC[16][4] = 0.7388;
4092
4093 fC[17][0] = 0.005;
4094 fC[17][1] = 0.005;
4095 fC[17][2] = 1.0000;
4096 fC[17][3] = 0.4411;
4097 fC[17][4] = 0.7293;
4098 }
4099 // 30-40%
4100 else if(centrCur < 40){
4101 fC[0][0] = 0.005;
4102 fC[0][1] = 0.005;
4103 fC[0][2] = 1.0000;
4104 fC[0][3] = 0.0010;
4105 fC[0][4] = 0.0010;
4106
4107 fC[1][0] = 0.005;
4108 fC[1][1] = 0.005;
4109 fC[1][2] = 1.0000;
4110 fC[1][3] = 0.0102;
4111 fC[1][4] = 0.0048;
4112
4113 fC[2][0] = 0.005;
4114 fC[2][1] = 0.005;
4115 fC[2][2] = 1.0000;
4116 fC[2][3] = 0.0306;
4117 fC[2][4] = 0.0079;
4118
4119 fC[3][0] = 0.005;
4120 fC[3][1] = 0.005;
4121 fC[3][2] = 1.0000;
4122 fC[3][3] = 0.0617;
4123 fC[3][4] = 0.0338;
4124
4125 fC[4][0] = 0.005;
4126 fC[4][1] = 0.005;
4127 fC[4][2] = 1.0000;
4128 fC[4][3] = 0.0920;
4129 fC[4][4] = 0.0652;
4130
4131 fC[5][0] = 0.005;
4132 fC[5][1] = 0.005;
4133 fC[5][2] = 1.0000;
4134 fC[5][3] = 0.1211;
4135 fC[5][4] = 0.0955;
4136
4137 fC[6][0] = 0.005;
4138 fC[6][1] = 0.005;
4139 fC[6][2] = 1.0000;
4140 fC[6][3] = 0.1496;
4141 fC[6][4] = 0.1242;
4142
4143 fC[7][0] = 0.005;
4144 fC[7][1] = 0.005;
4145 fC[7][2] = 1.0000;
4146 fC[7][3] = 0.1807;
4147 fC[7][4] = 0.1576;
4148
4149 fC[8][0] = 0.005;
4150 fC[8][1] = 0.005;
4151 fC[8][2] = 1.0000;
4152 fC[8][3] = 0.2195;
4153 fC[8][4] = 0.2097;
4154
4155 fC[9][0] = 0.005;
4156 fC[9][1] = 0.005;
4157 fC[9][2] = 1.0000;
4158 fC[9][3] = 0.2732;
4159 fC[9][4] = 0.2884;
4160
4161 fC[10][0] = 0.005;
4162 fC[10][1] = 0.005;
4163 fC[10][2] = 1.0000;
4164 fC[10][3] = 0.3204;
4165 fC[10][4] = 0.3679;
4166
4167 fC[11][0] = 0.005;
4168 fC[11][1] = 0.005;
4169 fC[11][2] = 1.0000;
4170 fC[11][3] = 0.3564;
4171 fC[11][4] = 0.4449;
4172
4173 fC[12][0] = 0.005;
4174 fC[12][1] = 0.005;
4175 fC[12][2] = 1.0000;
4176 fC[12][3] = 0.3791;
4177 fC[12][4] = 0.5052;
4178
4179 fC[13][0] = 0.005;
4180 fC[13][1] = 0.005;
4181 fC[13][2] = 1.0000;
4182 fC[13][3] = 0.4062;
4183 fC[13][4] = 0.5647;
4184
4185 fC[14][0] = 0.005;
4186 fC[14][1] = 0.005;
4187 fC[14][2] = 1.0000;
4188 fC[14][3] = 0.4234;
4189 fC[14][4] = 0.6203;
4190
4191 fC[15][0] = 0.005;
4192 fC[15][1] = 0.005;
4193 fC[15][2] = 1.0000;
4194 fC[15][3] = 0.4441;
4195 fC[15][4] = 0.6381;
4196
4197 fC[16][0] = 0.005;
4198 fC[16][1] = 0.005;
4199 fC[16][2] = 1.0000;
4200 fC[16][3] = 0.4629;
4201 fC[16][4] = 0.6496;
4202
4203 fC[17][0] = 0.005;
4204 fC[17][1] = 0.005;
4205 fC[17][2] = 1.0000;
4206 fC[17][3] = 0.4293;
4207 fC[17][4] = 0.6491;
4208 }
4209 // 40-50%
4210 else if(centrCur < 50){
4211 fC[0][0] = 0.005;
4212 fC[0][1] = 0.005;
4213 fC[0][2] = 1.0000;
4214 fC[0][3] = 0.0010;
4215 fC[0][4] = 0.0010;
4216
4217 fC[1][0] = 0.005;
4218 fC[1][1] = 0.005;
4219 fC[1][2] = 1.0000;
4220 fC[1][3] = 0.0093;
4221 fC[1][4] = 0.0057;
4222
4223 fC[2][0] = 0.005;
4224 fC[2][1] = 0.005;
4225 fC[2][2] = 1.0000;
4226 fC[2][3] = 0.0319;
4227 fC[2][4] = 0.0075;
4228
4229 fC[3][0] = 0.005;
4230 fC[3][1] = 0.005;
4231 fC[3][2] = 1.0000;
4232 fC[3][3] = 0.0639;
4233 fC[3][4] = 0.0371;
4234
4235 fC[4][0] = 0.005;
4236 fC[4][1] = 0.005;
4237 fC[4][2] = 1.0000;
4238 fC[4][3] = 0.0939;
4239 fC[4][4] = 0.0725;
4240
4241 fC[5][0] = 0.005;
4242 fC[5][1] = 0.005;
4243 fC[5][2] = 1.0000;
4244 fC[5][3] = 0.1224;
4245 fC[5][4] = 0.1045;
4246
4247 fC[6][0] = 0.005;
4248 fC[6][1] = 0.005;
4249 fC[6][2] = 1.0000;
4250 fC[6][3] = 0.1520;
4251 fC[6][4] = 0.1387;
4252
4253 fC[7][0] = 0.005;
4254 fC[7][1] = 0.005;
4255 fC[7][2] = 1.0000;
4256 fC[7][3] = 0.1783;
4257 fC[7][4] = 0.1711;
4258
4259 fC[8][0] = 0.005;
4260 fC[8][1] = 0.005;
4261 fC[8][2] = 1.0000;
4262 fC[8][3] = 0.2202;
4263 fC[8][4] = 0.2269;
4264
4265 fC[9][0] = 0.005;
4266 fC[9][1] = 0.005;
4267 fC[9][2] = 1.0000;
4268 fC[9][3] = 0.2672;
4269 fC[9][4] = 0.2955;
4270
4271 fC[10][0] = 0.005;
4272 fC[10][1] = 0.005;
4273 fC[10][2] = 1.0000;
4274 fC[10][3] = 0.3191;
4275 fC[10][4] = 0.3676;
4276
4277 fC[11][0] = 0.005;
4278 fC[11][1] = 0.005;
4279 fC[11][2] = 1.0000;
4280 fC[11][3] = 0.3434;
4281 fC[11][4] = 0.4321;
4282
4283 fC[12][0] = 0.005;
4284 fC[12][1] = 0.005;
4285 fC[12][2] = 1.0000;
4286 fC[12][3] = 0.3692;
4287 fC[12][4] = 0.4879;
4288
4289 fC[13][0] = 0.005;
4290 fC[13][1] = 0.005;
4291 fC[13][2] = 1.0000;
4292 fC[13][3] = 0.3993;
4293 fC[13][4] = 0.5377;
4294
4295 fC[14][0] = 0.005;
4296 fC[14][1] = 0.005;
4297 fC[14][2] = 1.0000;
4298 fC[14][3] = 0.3818;
4299 fC[14][4] = 0.5547;
4300
4301 fC[15][0] = 0.005;
4302 fC[15][1] = 0.005;
4303 fC[15][2] = 1.0000;
4304 fC[15][3] = 0.4003;
4305 fC[15][4] = 0.5484;
4306
4307 fC[16][0] = 0.005;
4308 fC[16][1] = 0.005;
4309 fC[16][2] = 1.0000;
4310 fC[16][3] = 0.4281;
4311 fC[16][4] = 0.5383;
4312
4313 fC[17][0] = 0.005;
4314 fC[17][1] = 0.005;
4315 fC[17][2] = 1.0000;
4316 fC[17][3] = 0.3960;
4317 fC[17][4] = 0.5374;
4318 }
4319 // 50-60%
4320 else if(centrCur < 60){
4321 fC[0][0] = 0.005;
4322 fC[0][1] = 0.005;
4323 fC[0][2] = 1.0000;
4324 fC[0][3] = 0.0010;
4325 fC[0][4] = 0.0010;
4326
4327 fC[1][0] = 0.005;
4328 fC[1][1] = 0.005;
4329 fC[1][2] = 1.0000;
4330 fC[1][3] = 0.0076;
4331 fC[1][4] = 0.0032;
4332
4333 fC[2][0] = 0.005;
4334 fC[2][1] = 0.005;
4335 fC[2][2] = 1.0000;
4336 fC[2][3] = 0.0329;
4337 fC[2][4] = 0.0085;
4338
4339 fC[3][0] = 0.005;
4340 fC[3][1] = 0.005;
4341 fC[3][2] = 1.0000;
4342 fC[3][3] = 0.0653;
4343 fC[3][4] = 0.0423;
4344
4345 fC[4][0] = 0.005;
4346 fC[4][1] = 0.005;
4347 fC[4][2] = 1.0000;
4348 fC[4][3] = 0.0923;
4349 fC[4][4] = 0.0813;
4350
4351 fC[5][0] = 0.005;
4352 fC[5][1] = 0.005;
4353 fC[5][2] = 1.0000;
4354 fC[5][3] = 0.1219;
4355 fC[5][4] = 0.1161;
4356
4357 fC[6][0] = 0.005;
4358 fC[6][1] = 0.005;
4359 fC[6][2] = 1.0000;
4360 fC[6][3] = 0.1519;
4361 fC[6][4] = 0.1520;
4362
4363 fC[7][0] = 0.005;
4364 fC[7][1] = 0.005;
4365 fC[7][2] = 1.0000;
4366 fC[7][3] = 0.1763;
4367 fC[7][4] = 0.1858;
4368
4369 fC[8][0] = 0.005;
4370 fC[8][1] = 0.005;
4371 fC[8][2] = 1.0000;
4372 fC[8][3] = 0.2178;
4373 fC[8][4] = 0.2385;
4374
4375 fC[9][0] = 0.005;
4376 fC[9][1] = 0.005;
4377 fC[9][2] = 1.0000;
4378 fC[9][3] = 0.2618;
4379 fC[9][4] = 0.3070;
4380
4381 fC[10][0] = 0.005;
4382 fC[10][1] = 0.005;
4383 fC[10][2] = 1.0000;
4384 fC[10][3] = 0.3067;
4385 fC[10][4] = 0.3625;
4386
4387 fC[11][0] = 0.005;
4388 fC[11][1] = 0.005;
4389 fC[11][2] = 1.0000;
4390 fC[11][3] = 0.3336;
4391 fC[11][4] = 0.4188;
4392
4393 fC[12][0] = 0.005;
4394 fC[12][1] = 0.005;
4395 fC[12][2] = 1.0000;
4396 fC[12][3] = 0.3706;
4397 fC[12][4] = 0.4511;
4398
4399 fC[13][0] = 0.005;
4400 fC[13][1] = 0.005;
4401 fC[13][2] = 1.0000;
4402 fC[13][3] = 0.3765;
4403 fC[13][4] = 0.4729;
4404
4405 fC[14][0] = 0.005;
4406 fC[14][1] = 0.005;
4407 fC[14][2] = 1.0000;
4408 fC[14][3] = 0.3942;
4409 fC[14][4] = 0.4855;
4410
4411 fC[15][0] = 0.005;
4412 fC[15][1] = 0.005;
4413 fC[15][2] = 1.0000;
4414 fC[15][3] = 0.4051;
4415 fC[15][4] = 0.4762;
4416
4417 fC[16][0] = 0.005;
4418 fC[16][1] = 0.005;
4419 fC[16][2] = 1.0000;
4420 fC[16][3] = 0.3843;
4421 fC[16][4] = 0.4763;
4422
4423 fC[17][0] = 0.005;
4424 fC[17][1] = 0.005;
4425 fC[17][2] = 1.0000;
4426 fC[17][3] = 0.4237;
4427 fC[17][4] = 0.4773;
4428 }
4429 // 60-70%
4430 else if(centrCur < 70){
4431 fC[0][0] = 0.005;
4432 fC[0][1] = 0.005;
4433 fC[0][2] = 1.0000;
4434 fC[0][3] = 0.0010;
4435 fC[0][4] = 0.0010;
4436
4437 fC[1][0] = 0.005;
4438 fC[1][1] = 0.005;
4439 fC[1][2] = 1.0000;
4440 fC[1][3] = 0.0071;
4441 fC[1][4] = 0.0012;
4442
4443 fC[2][0] = 0.005;
4444 fC[2][1] = 0.005;
4445 fC[2][2] = 1.0000;
4446 fC[2][3] = 0.0336;
4447 fC[2][4] = 0.0097;
4448
4449 fC[3][0] = 0.005;
4450 fC[3][1] = 0.005;
4451 fC[3][2] = 1.0000;
4452 fC[3][3] = 0.0662;
4453 fC[3][4] = 0.0460;
4454
4455 fC[4][0] = 0.005;
4456 fC[4][1] = 0.005;
4457 fC[4][2] = 1.0000;
4458 fC[4][3] = 0.0954;
4459 fC[4][4] = 0.0902;
4460
4461 fC[5][0] = 0.005;
4462 fC[5][1] = 0.005;
4463 fC[5][2] = 1.0000;
4464 fC[5][3] = 0.1181;
4465 fC[5][4] = 0.1306;
4466
4467 fC[6][0] = 0.005;
4468 fC[6][1] = 0.005;
4469 fC[6][2] = 1.0000;
4470 fC[6][3] = 0.1481;
4471 fC[6][4] = 0.1662;
4472
4473 fC[7][0] = 0.005;
4474 fC[7][1] = 0.005;
4475 fC[7][2] = 1.0000;
4476 fC[7][3] = 0.1765;
4477 fC[7][4] = 0.1963;
4478
4479 fC[8][0] = 0.005;
4480 fC[8][1] = 0.005;
4481 fC[8][2] = 1.0000;
4482 fC[8][3] = 0.2155;
4483 fC[8][4] = 0.2433;
4484
4485 fC[9][0] = 0.005;
4486 fC[9][1] = 0.005;
4487 fC[9][2] = 1.0000;
4488 fC[9][3] = 0.2580;
4489 fC[9][4] = 0.3022;
4490
4491 fC[10][0] = 0.005;
4492 fC[10][1] = 0.005;
4493 fC[10][2] = 1.0000;
4494 fC[10][3] = 0.2872;
4495 fC[10][4] = 0.3481;
4496
4497 fC[11][0] = 0.005;
4498 fC[11][1] = 0.005;
4499 fC[11][2] = 1.0000;
4500 fC[11][3] = 0.3170;
4501 fC[11][4] = 0.3847;
4502
4503 fC[12][0] = 0.005;
4504 fC[12][1] = 0.005;
4505 fC[12][2] = 1.0000;
4506 fC[12][3] = 0.3454;
4507 fC[12][4] = 0.4258;
4508
4509 fC[13][0] = 0.005;
4510 fC[13][1] = 0.005;
4511 fC[13][2] = 1.0000;
4512 fC[13][3] = 0.3580;
4513 fC[13][4] = 0.4299;
4514
4515 fC[14][0] = 0.005;
4516 fC[14][1] = 0.005;
4517 fC[14][2] = 1.0000;
4518 fC[14][3] = 0.3903;
4519 fC[14][4] = 0.4326;
4520
4521 fC[15][0] = 0.005;
4522 fC[15][1] = 0.005;
4523 fC[15][2] = 1.0000;
4524 fC[15][3] = 0.3690;
4525 fC[15][4] = 0.4491;
4526
4527 fC[16][0] = 0.005;
4528 fC[16][1] = 0.005;
4529 fC[16][2] = 1.0000;
4530 fC[16][3] = 0.4716;
4531 fC[16][4] = 0.4298;
4532
4533 fC[17][0] = 0.005;
4534 fC[17][1] = 0.005;
4535 fC[17][2] = 1.0000;
4536 fC[17][3] = 0.3875;
4537 fC[17][4] = 0.4083;
4538 }
4539 // 70-80%
4540 else if(centrCur < 80){
4541 fC[0][0] = 0.005;
4542 fC[0][1] = 0.005;
4543 fC[0][2] = 1.0000;
4544 fC[0][3] = 0.0010;
4545 fC[0][4] = 0.0010;
4546
4547 fC[1][0] = 0.005;
4548 fC[1][1] = 0.005;
4549 fC[1][2] = 1.0000;
4550 fC[1][3] = 0.0075;
4551 fC[1][4] = 0.0007;
4552
4553 fC[2][0] = 0.005;
4554 fC[2][1] = 0.005;
4555 fC[2][2] = 1.0000;
4556 fC[2][3] = 0.0313;
4557 fC[2][4] = 0.0124;
4558
4559 fC[3][0] = 0.005;
4560 fC[3][1] = 0.005;
4561 fC[3][2] = 1.0000;
4562 fC[3][3] = 0.0640;
4563 fC[3][4] = 0.0539;
4564
4565 fC[4][0] = 0.005;
4566 fC[4][1] = 0.005;
4567 fC[4][2] = 1.0000;
4568 fC[4][3] = 0.0923;
4569 fC[4][4] = 0.0992;
4570
4571 fC[5][0] = 0.005;
4572 fC[5][1] = 0.005;
4573 fC[5][2] = 1.0000;
4574 fC[5][3] = 0.1202;
4575 fC[5][4] = 0.1417;
4576
4577 fC[6][0] = 0.005;
4578 fC[6][1] = 0.005;
4579 fC[6][2] = 1.0000;
4580 fC[6][3] = 0.1413;
4581 fC[6][4] = 0.1729;
4582
4583 fC[7][0] = 0.005;
4584 fC[7][1] = 0.005;
4585 fC[7][2] = 1.0000;
4586 fC[7][3] = 0.1705;
4587 fC[7][4] = 0.1999;
4588
4589 fC[8][0] = 0.005;
4590 fC[8][1] = 0.005;
4591 fC[8][2] = 1.0000;
4592 fC[8][3] = 0.2103;
4593 fC[8][4] = 0.2472;
4594
4595 fC[9][0] = 0.005;
4596 fC[9][1] = 0.005;
4597 fC[9][2] = 1.0000;
4598 fC[9][3] = 0.2373;
4599 fC[9][4] = 0.2916;
4600
4601 fC[10][0] = 0.005;
4602 fC[10][1] = 0.005;
4603 fC[10][2] = 1.0000;
4604 fC[10][3] = 0.2824;
4605 fC[10][4] = 0.3323;
4606
4607 fC[11][0] = 0.005;
4608 fC[11][1] = 0.005;
4609 fC[11][2] = 1.0000;
4610 fC[11][3] = 0.3046;
4611 fC[11][4] = 0.3576;
4612
4613 fC[12][0] = 0.005;
4614 fC[12][1] = 0.005;
4615 fC[12][2] = 1.0000;
4616 fC[12][3] = 0.3585;
4617 fC[12][4] = 0.4003;
4618
4619 fC[13][0] = 0.005;
4620 fC[13][1] = 0.005;
4621 fC[13][2] = 1.0000;
4622 fC[13][3] = 0.3461;
4623 fC[13][4] = 0.3982;
4624
4625 fC[14][0] = 0.005;
4626 fC[14][1] = 0.005;
4627 fC[14][2] = 1.0000;
4628 fC[14][3] = 0.3362;
4629 fC[14][4] = 0.3776;
4630
4631 fC[15][0] = 0.005;
4632 fC[15][1] = 0.005;
4633 fC[15][2] = 1.0000;
4634 fC[15][3] = 0.3071;
4635 fC[15][4] = 0.3500;
4636
4637 fC[16][0] = 0.005;
4638 fC[16][1] = 0.005;
4639 fC[16][2] = 1.0000;
4640 fC[16][3] = 0.2914;
4641 fC[16][4] = 0.3937;
4642
4643 fC[17][0] = 0.005;
4644 fC[17][1] = 0.005;
4645 fC[17][2] = 1.0000;
4646 fC[17][3] = 0.3727;
4647 fC[17][4] = 0.3877;
4648 }
4649 // 80-100%
4650 else{
4651 fC[0][0] = 0.005;
4652 fC[0][1] = 0.005;
4653 fC[0][2] = 1.0000;
4654 fC[0][3] = 0.0010;
4655 fC[0][4] = 0.0010;
4656
4657 fC[1][0] = 0.005;
4658 fC[1][1] = 0.005;
4659 fC[1][2] = 1.0000;
4660 fC[1][3] = 0.0060;
4661 fC[1][4] = 0.0035;
4662
4663 fC[2][0] = 0.005;
4664 fC[2][1] = 0.005;
4665 fC[2][2] = 1.0000;
4666 fC[2][3] = 0.0323;
4667 fC[2][4] = 0.0113;
4668
4669 fC[3][0] = 0.005;
4670 fC[3][1] = 0.005;
4671 fC[3][2] = 1.0000;
4672 fC[3][3] = 0.0609;
4673 fC[3][4] = 0.0653;
4674
4675 fC[4][0] = 0.005;
4676 fC[4][1] = 0.005;
4677 fC[4][2] = 1.0000;
4678 fC[4][3] = 0.0922;
4679 fC[4][4] = 0.1076;
4680
4681 fC[5][0] = 0.005;
4682 fC[5][1] = 0.005;
4683 fC[5][2] = 1.0000;
4684 fC[5][3] = 0.1096;
4685 fC[5][4] = 0.1328;
4686
4687 fC[6][0] = 0.005;
4688 fC[6][1] = 0.005;
4689 fC[6][2] = 1.0000;
4690 fC[6][3] = 0.1495;
4691 fC[6][4] = 0.1779;
4692
4693 fC[7][0] = 0.005;
4694 fC[7][1] = 0.005;
4695 fC[7][2] = 1.0000;
4696 fC[7][3] = 0.1519;
4697 fC[7][4] = 0.1989;
4698
4699 fC[8][0] = 0.005;
4700 fC[8][1] = 0.005;
4701 fC[8][2] = 1.0000;
4702 fC[8][3] = 0.1817;
4703 fC[8][4] = 0.2472;
4704
4705 fC[9][0] = 0.005;
4706 fC[9][1] = 0.005;
4707 fC[9][2] = 1.0000;
4708 fC[9][3] = 0.2429;
4709 fC[9][4] = 0.2684;
4710
4711 fC[10][0] = 0.005;
4712 fC[10][1] = 0.005;
4713 fC[10][2] = 1.0000;
4714 fC[10][3] = 0.2760;
4715 fC[10][4] = 0.3098;
4716
4717 fC[11][0] = 0.005;
4718 fC[11][1] = 0.005;
4719 fC[11][2] = 1.0000;
4720 fC[11][3] = 0.2673;
4721 fC[11][4] = 0.3198;
4722
4723 fC[12][0] = 0.005;
4724 fC[12][1] = 0.005;
4725 fC[12][2] = 1.0000;
4726 fC[12][3] = 0.3165;
4727 fC[12][4] = 0.3564;
4728
4729 fC[13][0] = 0.005;
4730 fC[13][1] = 0.005;
4731 fC[13][2] = 1.0000;
4732 fC[13][3] = 0.3526;
4733 fC[13][4] = 0.3011;
4734
4735 fC[14][0] = 0.005;
4736 fC[14][1] = 0.005;
4737 fC[14][2] = 1.0000;
4738 fC[14][3] = 0.3788;
4739 fC[14][4] = 0.3011;
4740
4741 fC[15][0] = 0.005;
4742 fC[15][1] = 0.005;
4743 fC[15][2] = 1.0000;
4744 fC[15][3] = 0.3788;
4745 fC[15][4] = 0.3011;
4746
4747 fC[16][0] = 0.005;
4748 fC[16][1] = 0.005;
4749 fC[16][2] = 1.0000;
4750 fC[16][3] = 0.3788;
4751 fC[16][4] = 0.3011;
4752
4753 fC[17][0] = 0.005;
4754 fC[17][1] = 0.005;
4755 fC[17][2] = 1.0000;
4756 fC[17][3] = 0.3788;
4757 fC[17][4] = 0.3011;
4758 }
4759
4760 for(Int_t i=18;i<fgkPIDptBin;i++){
4761 fBinLimitPID[i] = 3.0 + 0.2 * (i-18);
4762 fC[i][0] = fC[17][0];
4763 fC[i][1] = fC[17][1];
4764 fC[i][2] = fC[17][2];
4765 fC[i][3] = fC[17][3];
4766 fC[i][4] = fC[17][4];
4767 }
4768}
4769
4770//---------------------------------------------------------------//
0d16f553 4771Bool_t AliFlowTrackCuts::TPCTOFagree(const AliVTrack *track)
875ee12a 4772{
3c67c769 4773 //check pid agreement between TPC and TOF
875ee12a 4774 Bool_t status = kFALSE;
0d16f553 4775
4776 const Float_t c = 2.99792457999999984e-02;
4777
875ee12a 4778 Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
4779
4780
80efa607 4781 Double_t exptimes[9];
875ee12a 4782 track->GetIntegratedTimes(exptimes);
4783
4784 Float_t dedx = track->GetTPCsignal();
4785
4786 Float_t p = track->P();
4787 Float_t time = track->GetTOFsignal()- fESDpid.GetTOFResponse().GetStartTime(p);
0d16f553 4788 Float_t tl = exptimes[0]*c; // to work both for ESD and AOD
875ee12a 4789
4790 Float_t betagammares = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
4791
4792 Float_t betagamma1 = tl/(time-5 *betagammares) * 33.3564095198152043;
4793
4794// printf("betagamma1 = %f\n",betagamma1);
4795
4796 if(betagamma1 < 0.1) betagamma1 = 0.1;
4797
4798 if(betagamma1 < 0.99999) betagamma1 /= TMath::Sqrt(1-betagamma1*betagamma1);
4799 else betagamma1 = 100;
4800
4801 Float_t betagamma2 = tl/(time+5 *betagammares) * 33.3564095198152043;
4802// printf("betagamma2 = %f\n",betagamma2);
4803
4804 if(betagamma2 < 0.1) betagamma2 = 0.1;
4805
4806 if(betagamma2 < 0.99999) betagamma2 /= TMath::Sqrt(1-betagamma2*betagamma2);
4807 else betagamma2 = 100;
4808
4809
0d16f553 4810 Float_t momtpc=track->GetTPCmomentum();
875ee12a 4811
4812 for(Int_t i=0;i < 5;i++){
4813 Float_t resolutionTOF = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[i], mass[i]);
4814 if(TMath::Abs(exptimes[i] - time) < 5 * resolutionTOF){
4815 Float_t dedxExp = 0;
4816 if(i==0) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
4817 else if(i==1) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
4818 else if(i==2) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
4819 else if(i==3) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
4820 else if(i==4) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
4821
4822 Float_t resolutionTPC = 2;
4823 if(i==0) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kElectron);
4824 else if(i==1) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kMuon);
4825 else if(i==2) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kPion);
4826 else if(i==3) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kKaon);
4827 else if(i==4) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
4828
4829 if(TMath::Abs(dedx - dedxExp) < 3 * resolutionTPC){
4830 status = kTRUE;
4831 }
4832 }
4833 }
4834
4835 Float_t bb1 = fESDpid.GetTPCResponse().Bethe(betagamma1);
4836 Float_t bb2 = fESDpid.GetTPCResponse().Bethe(betagamma2);
4837 Float_t bbM = fESDpid.GetTPCResponse().Bethe((betagamma1+betagamma2)*0.5);
4838
4839
4840 // status = kFALSE;
4841 // for nuclei
4842 Float_t resolutionTOFpr = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
4843 Float_t resolutionTPCpr = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
4844 if(TMath::Abs(dedx-bb1) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4845 status = kTRUE;
4846 }
4847 else if(TMath::Abs(dedx-bb2) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4848 status = kTRUE;
4849 }
4850 else if(TMath::Abs(dedx-bbM) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
4851 status = kTRUE;
4852 }
4853
4854 return status;
3abf7ecc 4855}
4856// end part added by F. Noferini
4857//-----------------------------------------------------------------------
2b1eaa10 4858
2b1eaa10 4859//-----------------------------------------------------------------------
4860const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
4861{
4862 //get the name of the particle id source
4863 switch (s)
4864 {
4865 case kTPCdedx:
4866 return "TPCdedx";
b9cf8f8e 4867 case kTPCbayesian:
4868 return "TPCbayesian";
2b1eaa10 4869 case kTOFbeta:
4870 return "TOFbeta";
4871 case kTPCpid:
4872 return "TPCpid";
4873 case kTOFpid:
4874 return "TOFpid";
4875 case kTOFbayesian:
4876 return "TOFbayesianPID";
1a80f9f6 4877 case kTOFbetaSimple:
4878 return "TOFbetaSimple";
8eca5d19 4879 case kTPCNuclei:
4880 return "TPCnuclei";
d0ce6147 4881 case kTPCTOFNsigma:
4882 return "TPCTOFNsigma";
67006c99 4883 case kTPCTOFNsigmaPuritybased:
4884 return "TPCTOFNsigmaPuritybased";
2b1eaa10 4885 default:
4886 return "NOPID";
4887 }
4888}
4889
4890//-----------------------------------------------------------------------
4891const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
4892{
4893 //return the name of the selected parameter type
4894 switch (type)
4895 {
4896 case kMC:
4897 return "MC";
4898 case kGlobal:
1a80f9f6 4899 return "Global";
4900 case kTPCstandalone:
4901 return "TPCstandalone";
4902 case kSPDtracklet:
4903 return "SPDtracklets";
d148af7e 4904 case kPMD:
4905 return "PMD";
c27333c2
MK		 4906 case kVZERO:
4907 return "VZERO";
2e5052c5 4908 case kMUON: // XZhang 20120604
4909 return "MUON"; // XZhang 20120604
dd294bb6
MK		 4910 case kKink:
4911 return "Kink";
4912 case kV0:
4913 return "V0";
2b1eaa10 4914 default:
4915 return "unknown";
4916 }
4917}
4918
d148af7e 4919//-----------------------------------------------------------------------
3c67c769 4920Bool_t AliFlowTrackCuts::PassesPMDcuts(const AliESDPmdTrack* track )
d148af7e 4921{
4922 //check PMD specific cuts
4923 //clean up from last iteration, and init label
1a80f9f6 4924 Int_t det = track->GetDetector();
4925 //Int_t smn = track->GetSmn();
4926 Float_t clsX = track->GetClusterX();
4927 Float_t clsY = track->GetClusterY();
4928 Float_t clsZ = track->GetClusterZ();
4929 Float_t ncell = track->GetClusterCells();
4930 Float_t adc = track->GetClusterADC();
d148af7e 4931
4932 fTrack = NULL;
4933 fMCparticle=NULL;
2281e7c5 4934 fTrackLabel=-996;
d148af7e 4935
1a80f9f6 4936 fTrackEta = GetPmdEta(clsX,clsY,clsZ);
4937 fTrackPhi = GetPmdPhi(clsX,clsY);
d148af7e 4938 fTrackWeight = 1.0;
4939
4940 Bool_t pass=kTRUE;
4941 if (fCutEta) {if ( fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
4942 if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
1a80f9f6 4943 if (fCutPmdDet) {if(det != fPmdDet) pass = kFALSE;}
4944 if (fCutPmdAdc) {if(adc < fPmdAdc) pass = kFALSE;}
4945 if (fCutPmdNcell) {if(ncell < fPmdNcell) pass = kFALSE;}
d148af7e 4946
1a80f9f6 4947 return pass;
d148af7e 4948}
22289738 4949
4950//-----------------------------------------------------------------------
c27333c2 4951Bool_t AliFlowTrackCuts::PassesVZEROcuts(Int_t id)
22289738 4952{
c27333c2 4953 //check VZERO cuts
b35580e6 4954 if (id<0) return kFALSE;
4955
22289738 4956 //clean up from last iter
dd294bb6 4957 ClearTrack();
22289738 4958
1a80f9f6 4959 fTrackPhi = TMath::PiOver4()*(0.5+id%8);
647a09d3 4960
86a97121 4961 // 10102013 weighting vzero tiles - rbertens@cern.ch
c27333c2 4962 if(!fVZEROgainEqualization) {
86a97121 4963 // if for some reason the equalization is not initialized (e.g. 2011 data)
c27333c2 4964 // the fVZEROxpol[] weights are used to enable or disable vzero rings
54dd223b 4965 if(id<32) { // v0c side
86a97121 4966 fTrackEta = -3.45+0.5*(id/8);
c27333c2
MK		 4967 if(id < 8) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[0];
4968 else if (id < 16 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[1];
4969 else if (id < 24 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[2];
4970 else if (id < 32 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[3];
86a97121 4971 } else { // v0a side
4972 fTrackEta = +4.8-0.6*((id/8)-4);
c27333c2
MK		 4973 if( id < 40) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[0];
4974 else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[1];
4975 else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[2];
4976 else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[3];
86a97121 4977 }
4978 } else { // the equalization is initialized
4979 // note that disabled rings have already been excluded on calibration level in
c27333c2 4980 // AliFlowEvent (so for a disabled ring, fVZEROxpol is zero
86a97121 4981 if(id<32) { // v0c side
4982 fTrackEta = -3.45+0.5*(id/8);
c27333c2
MK		 4983 if(id < 8) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[0]/fVZEROgainEqualization->GetBinContent(1+id);
4984 else if (id < 16 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[1]/fVZEROgainEqualization->GetBinContent(1+id);
4985 else if (id < 24 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[2]/fVZEROgainEqualization->GetBinContent(1+id);
4986 else if (id < 32 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[3]/fVZEROgainEqualization->GetBinContent(1+id);
86a97121 4987 } else { // v0a side
4988 fTrackEta = +4.8-0.6*((id/8)-4);
c27333c2
MK		 4989 if( id < 40) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[0]/fVZEROgainEqualization->GetBinContent(1+id);
4990 else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[1]/fVZEROgainEqualization->GetBinContent(1+id);
4991 else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[2]/fVZEROgainEqualization->GetBinContent(1+id);
4992 else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[3]/fVZEROgainEqualization->GetBinContent(1+id);
86a97121 4993 }
4994 // printf ( " tile %i and weight %.2f \n", id, fTrackWeight);
4995 }
647a09d3 4996
1d78aabb 4997 if (fLinearizeVZEROresponse && id < 64)
92b507bb 4998 {
4999 //this is only needed in pass1 of LHC10h
c27333c2 5000 Float_t multVZERO[fgkNumberOfVZEROtracks];
92b507bb 5001 Float_t dummy=0.0;
cac040c5 5002 AliESDUtils::GetCorrV0((AliESDEvent*)fEvent,dummy,multVZERO);
c27333c2 5003 fTrackWeight = multVZERO[id];
92b507bb 5004 }
22289738 5005
5006 Bool_t pass=kTRUE;
5007 if (fCutEta) {if ( fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
5008 if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
5009
1a80f9f6 5010 return pass;
22289738 5011}
5012
2e5052c5 5013//-----------------------------------------------------------------------------
5014Bool_t AliFlowTrackCuts::PassesMuonCuts(AliVParticle* vparticle)
5015{
5016// XZhang 20120604
dd294bb6 5017 ClearTrack();
2e5052c5 5018 fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
5019 if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
5020 else fMCparticle=NULL;
5021
5022 AliESDMuonTrack *esdTrack = dynamic_cast<AliESDMuonTrack*>(vparticle);
5023 AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(vparticle);
5024 if ((!esdTrack) && (!aodTrack)) return kFALSE;
5025 if (!fMuonTrackCuts->IsSelected(vparticle)) return kFALSE;
5026 HandleVParticle(vparticle); if (!fTrack) return kFALSE;
5027 return kTRUE;
5028}
5029
1a80f9f6 5030//----------------------------------------------------------------------------//
5031Double_t AliFlowTrackCuts::GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
5032{
3c67c769 5033 //get PMD "track" eta coordinate
1a80f9f6 5034 Float_t rpxpy, theta, eta;
5035 rpxpy = TMath::Sqrt(xPos*xPos + yPos*yPos);
5036 theta = TMath::ATan2(rpxpy,zPos);
5037 eta = -TMath::Log(TMath::Tan(0.5*theta));
5038 return eta;
5039}
5040
5041//--------------------------------------------------------------------------//
5042Double_t AliFlowTrackCuts::GetPmdPhi(Float_t xPos, Float_t yPos)
5043{
3c67c769 5044 //Get PMD "track" phi coordinate
1a80f9f6 5045 Float_t pybypx, phi = 0., phi1;
5046 if(xPos==0)
5047 {
5048 if(yPos>0) phi = 90.;
5049 if(yPos<0) phi = 270.;
5050 }
5051 if(xPos != 0)
5052 {
5053 pybypx = yPos/xPos;
5054 if(pybypx < 0) pybypx = - pybypx;
5055 phi1 = TMath::ATan(pybypx)*180./3.14159;
5056
5057 if(xPos > 0 && yPos > 0) phi = phi1; // 1st Quadrant
5058 if(xPos < 0 && yPos > 0) phi = 180 - phi1; // 2nd Quadrant
5059 if(xPos < 0 && yPos < 0) phi = 180 + phi1; // 3rd Quadrant
5060 if(xPos > 0 && yPos < 0) phi = 360 - phi1; // 4th Quadrant
5061
5062 }
5063 phi = phi*3.14159/180.;
5064 return phi;
5065}
499fe731 5066
5067//---------------------------------------------------------------//
5068void AliFlowTrackCuts::Browse(TBrowser* b)
5069{
5070 //some browsing capabilities
5071 if (fQA) b->Add(fQA);
5072}
5073
1a80f9f6 5074//---------------------------------------------------------------//
499fe731 5075Long64_t AliFlowTrackCuts::Merge(TCollection* list)
5076{
5077 //merge
8eca5d19 5078 if (!fQA || !list) return 0;
5079 if (list->IsEmpty()) return 0;
5080 AliFlowTrackCuts* obj=NULL;
5081 TList tmplist;
499fe731 5082 TIter next(list);
5083 while ( (obj = dynamic_cast<AliFlowTrackCuts*>(next())) )
5084 {
5085 if (obj==this) continue;
8eca5d19 5086 tmplist.Add(obj->GetQA());
499fe731 5087 }
8eca5d19 5088 return fQA->Merge(&tmplist);
499fe731 5089}
67006c99 5090//________________________________________________________________//
5091void AliFlowTrackCuts::GetTPCTOFPIDContours()
5092{
5093 fContoursFile = TFile::Open(Form("$ALICE_ROOT/PWGCF/FLOW/database/PIDCutContours_%i-%i.root",fCentralityPercentileMin,fCentralityPercentileMax));
5094 fCutContourList=(TDirectory*)fContoursFile->Get("Filterbit1");
5095 if(!fCutContourList){printf("The contour file is empty"); return;}
5096
5097 Double_t pBinning[50];
5098 for(int b=0;b<50;b++){pBinning[b]=b;}
5099 TString species[3] = {"pion","kaon","proton"};
5100 TList *Species_contours[3];
5101 Int_t ispecie = 0;
5102 for(ispecie = 0; ispecie < 3; ispecie++) {
5103 Species_contours[ispecie] = (TList*)fCutContourList->Get(species[ispecie]);
5104 if(!Species_contours[ispecie]) {
5105 cout<<"Contours for species: "<<species[ispecie]<<" not found!!!"<<endl;
5106 return;
5107 }
5108 }
5109 Int_t iParticle = 0;
5110
5111 switch (fParticleID)
5112 {
5113 case AliPID::kPion:
5114 iParticle=0;
5115 break;
5116 case AliPID::kKaon:
5117 iParticle=1;
5118 break;
5119 case AliPID::kProton:
5120 iParticle=2;
5121 break;
5122 default:
5123 return;
5124 }
5125
5126
5127 for(int i=0;i<50;i++){
5128
5129 TString Graph_Name = "contourlines_";
5130 Graph_Name += species[iParticle];
5131 Graph_Name += Form("%.f%.f-%i%icent",pBinning[i],pBinning[i]+1,fCentralityPercentileMin,fCentralityPercentileMax);
5132 fCutGraph[i] = (TGraph*)Species_contours[iParticle]->FindObject(Graph_Name);
5133
5134 if(!fCutGraph[i]){cout<<"Contour Graph does not exist"<<endl; continue;}
5135
5136 fCutContour[i] = new TCutG(Graph_Name.Data(),fCutGraph[i]->GetN(),fCutGraph[i]->GetX(),fCutGraph[i]->GetY());
5137
5138 }
5139
5140}
5141//__________________
5142Int_t AliFlowTrackCuts::MaxSharedITSClusterCuts(Int_t maxITSclusterShared, AliESDtrack* track){
5143
5144 Int_t counterForSharedCluster = 0;
5145 for(int i =0;i<6;i++){
5146 Bool_t sharedITSCluster = track->HasSharedPointOnITSLayer(i);
5147 Bool_t PointOnITSLayer = track->HasPointOnITSLayer(i);
5148 // cout << "has a point??? " << PointOnITSLayer << " is it shared or not??? " << sharedITSCluster << endl;
5149 if(sharedITSCluster == 1) counterForSharedCluster++;
5150 }
5151 //if(counterForSharedCluster >= maxITSclusterShared) pass=kFALSE;
5152
5153 return counterForSharedCluster;
5154
5155}
5156//___________________
5157Double_t AliFlowTrackCuts::MaxChi2perITSClusterCuts(Double_t maxITSChi2, AliESDtrack* track){
5158
5159 // cout << " Total Shared Cluster == " <<counterForSharedCluster<< endl;
5160 Double_t chi2perClusterITS = track->GetITSchi2()/track->GetNcls(0);
5161 // cout << " chi2perClusterITS == " <<chi2perClusterITS <<endl;
5162
5163
5164 //if(chi2perClusterITS >= maxITSChi2) pass=kFALSE;
5165
5166 return chi2perClusterITS;
5167}
5168
f21bdf48 5169
9126fd0c 5170
2281e7c5 5171