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