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