AliPhysics  095eea3 (095eea3)
 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 Int_t AliFlowTrackCuts::GetITStype(const AliAODTrack* track) const
1384 {
1385  // set ITStype (test purpose only)
1386  Int_t ITStype=0;
1387  Bool_t bHITs[6] = {kFALSE};
1388  Int_t hitcnt=0;
1389  Bool_t isSPD=kFALSE, isSDD=kFALSE, isSSD=kFALSE;
1390  for (Int_t i=0; i<6; i++) {
1391  if(track->HasPointOnITSLayer(i)) {
1392  if(i<2) isSPD=kTRUE;
1393  if(i>2 && i<4) isSDD=kTRUE;
1394  if(i>4) isSSD=kTRUE;
1395  bHITs[i]=kTRUE;
1396  hitcnt++;
1397  }
1398  }
1399  if(hitcnt==1) {
1400  for (Int_t i=0; i<6; i++) {
1401  if(bHITs[i]) ITStype = i+1;
1402  }
1403  }
1404  if(hitcnt>=2) {
1405  if(isSPD & !isSDD & !isSSD) ITStype = 7;
1406  if(isSPD & isSDD & !isSSD) ITStype = 8;
1407  if(isSPD & !isSDD & isSSD) ITStype = 9;
1408  if(!isSPD & isSDD & !isSSD) ITStype = 10;
1409  if(!isSPD & isSDD & isSSD) ITStype = 11;
1410  if(!isSPD & !isSDD & isSSD) ITStype = 12;
1411  if(isSPD & isSDD & isSSD) {
1412  if(bHITs[0] && !bHITs[1]) ITStype = 13;
1413  if(!bHITs[0] && bHITs[1]) ITStype = 14;
1414  if(bHITs[0] && bHITs[1]) ITStype = 15;
1415  }
1416  }
1417  return ITStype;
1418 }
1419 
1420 //_______________________________________________________________________
1421 Bool_t AliFlowTrackCuts::PassesAODcuts(const AliAODTrack* track, Bool_t passedFid)
1422 {
1423  //check cuts for AOD
1424  Bool_t pass = passedFid;
1425 // AliAODPid *pidObj = track->GetDetPid();
1426 
1427  if (fCutNClustersTPC)
1428  {
1429  Int_t ntpccls = track->GetTPCNcls();
1430  if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
1431  }
1432 
1433  if (fCutNClustersITS)
1434  {
1435  Int_t nitscls = track->GetITSNcls();
1436  if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
1437  }
1438 
1439  if (fCutChi2PerClusterTPC)
1440  {
1441  Double_t chi2tpc = track->Chi2perNDF();
1442  if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
1443  }
1444 
1445  if (fCutChi2PerClusterITS)
1446  {
1447  Int_t nitscls = track->GetITSNcls();
1448  if(nitscls>0) {
1449  Double_t chi2TIS = track->GetITSchi2()/track->GetITSNcls();
1450  if (chi2TIS >= fMaxChi2PerClusterITS) pass=kFALSE;
1451  }
1452  }
1453 
1454  if (fCutITSClusterGlobal)
1455  {
1456  Bool_t bSPDCl = kFALSE;
1457  for (Int_t i=0; i<2; i++) {
1458  if(track->HasPointOnITSLayer(i)) bSPDCl = kTRUE;
1459  }
1460  Bool_t bSDDCl = track->HasPointOnITSLayer(2);
1461  Bool_t temppass = kTRUE;
1462  if (bSPDCl || (!bSPDCl && bSDDCl)) temppass = kTRUE;
1463  else temppass = kFALSE;
1464  if(!temppass) pass=kFALSE;
1465  }
1466 
1467  if (fCutFracSharedTPCCluster)
1468  {
1469  Int_t ntpccls = track->GetTPCncls();
1470  if(ntpccls>0) {
1471  Int_t ntpcclsS = track->GetTPCnclsS();
1472  Double_t fshtpccls = 1.*ntpcclsS/ntpccls;
1473  if (fshtpccls > fMaxFracSharedTPCCluster) pass=kFALSE;
1474  }
1475  }
1476 
1477  if (fCutCrossedTPCRows)
1478  {
1479  Int_t nCrossedRows = track->GetTPCNCrossedRows();
1480  if (nCrossedRows <= fMinNCrossedRows) pass=kFALSE;
1481  Float_t CrossedRowsOverFindableClusters = track->GetTPCFoundFraction();
1482  if (CrossedRowsOverFindableClusters < fMinCrossedRowsOverFindableClusters) pass=kFALSE;
1483  }
1484 
1485  if (fCutGoldenChi2)
1486  {
1487  Double_t GoldenChi2 = track->GetChi2TPCConstrainedVsGlobal();
1488  if (GoldenChi2 >= fMaxGoldenChi2) pass=kFALSE;
1489  }
1490 
1491  if (fRequireTOFSignal)
1492  {
1493  if(TMath::Abs(track->GetTOFsignalDz())>10.) pass=kFALSE;
1494  if(track->GetTOFsignal() < 12000.) pass=kFALSE;
1495  if(track->GetTOFsignal() > 25000.) pass=kFALSE;
1496  }
1497 
1498  if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
1499  if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
1500 
1501  if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
1502  Double_t DCAxy = track->DCA();
1503  Double_t DCAz = track->ZAtDCA();
1504  if(fCutDCAToVertexXYAOD || fCutDCAToVertexZAOD || fCutDCAToVertexXYPtDepAOD) {
1505  if (std::abs((Int_t)DCAxy)==999 || std::abs((Int_t)DCAz)==999) {
1506  // re-evaluate the dca as it seems to not be natively present
1507  // allowed only for tracks inside the beam pipe
1508  Double_t pos[3] = {-99., -99., -99.};
1509  track->GetPosition(pos);
1510  if(pos[0]*pos[0]+pos[1]*pos[1] <= 3.*3.) {
1511  AliAODTrack copy(*track); // stack copy
1512  Double_t b[2] = {-99., -99.};
1513  Double_t bCov[3] = {-99., -99., -99.};
1514  if(copy.PropagateToDCA(fEvent->GetPrimaryVertex(), fEvent->GetMagneticField(), 100., b, bCov)) {
1515  DCAxy = b[0];
1516  DCAz = b[1];
1517  }
1518  }
1519  }
1520  if (fCutDCAToVertexXYAOD) {
1521  if (TMath::Abs(DCAxy)>fMaxDCAxyAOD) pass=kFALSE;
1522  }
1523  if (fCutDCAToVertexXYPtDepAOD) {
1524  Double_t MaxDCAPtDep = 2.4;
1525  if (fAODFilterBit!=128) MaxDCAPtDep = 0.0182+0.0350/pow(track->Pt(),1.01);
1526  else MaxDCAPtDep = 0.4+0.2/pow(track->Pt(),0.3);
1527  if (TMath::Abs(DCAxy)>MaxDCAPtDep) pass=kFALSE;
1528  }
1529  if (fCutDCAToVertexZAOD) {
1530  if (TMath::Abs(DCAz)>fMaxDCAzAOD) pass=kFALSE;
1531  }
1532  }
1533  Double_t dedx = track->GetTPCsignal();
1534  if(fCutMinimalTPCdedx) {
1535  if (dedx < fMinimalTPCdedx) pass=kFALSE;
1536  }
1537  Double_t time[9];
1538  track->GetIntegratedTimes(time);
1539 
1540  if(fCutTPCSecbound && track->Pt()>fCutTPCSecboundMinpt) {
1541  Double_t xyz[3]={-9999.,-9999.,-9999.};
1542  const double r = 84.; // TPC IROC radius
1543  if (!track->GetXYZatR(r, fEvent->GetMagneticField(), xyz, NULL)) pass=kFALSE;
1544  Double_t cra = TMath::ATan2(xyz[1],xyz[0]); // crossing angle
1545  Double_t dpe = 3.*TMath::TwoPi()/360.;// excluded region (\pm 3 degree)
1546  for(Int_t nb=-9; nb<=9; nb++) {
1547  Double_t bnp = nb*TMath::Pi()/9.; // TPC boundaries azimuthal angle
1548  if(cra<bnp+dpe && cra>bnp-dpe) pass=kFALSE;
1549  }
1550  }
1551  if(fCutTPCSecboundVar) {
1552  Double_t phimod = track->Phi();
1553  if(fEvent->GetMagneticField() < 0) // for negative polarity field
1554  phimod = TMath::TwoPi() - phimod;
1555  if(track->Charge() < 0) // for negative charge
1556  phimod = TMath::TwoPi() - phimod;
1557  if(phimod < 0)
1558  cout << "Warning!!!!! phi < 0: " << phimod << endl;
1559 
1560  phimod += TMath::Pi()/18.0; // to center gap in the middle
1561  phimod = fmod(phimod, TMath::Pi()/9.0);
1562  if(phimod < fPhiCutHigh->Eval(track->Pt()) && phimod > fPhiCutLow->Eval(track->Pt()))
1563  pass=kFALSE; // reject track
1564  }
1565 
1566  if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
1567  {
1568  if (!PassesAODpidCut(track)) pass=kFALSE;
1569  }
1570 
1571  if (fQA) {
1572  // changed 04062014 used to be filled before possible PID cut
1573  Double_t momTPC = track->GetTPCmomentum();
1574  QAbefore( 0)->Fill(momTPC,GetBeta(track, kTRUE));
1575  if(pass) QAafter( 0)->Fill(momTPC, GetBeta(track, kTRUE));
1576  QAbefore( 1)->Fill(momTPC,dedx);
1577  QAbefore( 5)->Fill(track->Pt(),DCAxy);
1578  QAbefore( 6)->Fill(track->Pt(),DCAz);
1579  if (pass) QAafter( 1)->Fill(momTPC,dedx);
1580  if (pass) QAafter( 5)->Fill(track->Pt(),DCAxy);
1581  if (pass) QAafter( 6)->Fill(track->Pt(),DCAz);
1582  QAbefore( 8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1583  if (pass) QAafter( 8)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kElectron]));
1584  QAbefore( 9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1585  if (pass) QAafter( 9)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kMuon]));
1586  QAbefore(10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1587  if (pass) QAafter( 10)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kPion]));
1588  QAbefore(11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1589  if (pass) QAafter( 11)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kKaon]));
1590  QAbefore(12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
1591  if (pass) QAafter( 12)->Fill(track->P(),(track->GetTOFsignal()-time[AliPID::kProton]));
1592  QAbefore(18)->Fill(track->P(),track->Chi2perNDF());
1593  if (pass) QAafter( 18)->Fill(track->P(),track->Chi2perNDF());
1594  Int_t ntpccls = track->GetTPCncls();
1595  if(ntpccls>0) {
1596  Int_t ntpcclsS = track->GetTPCnclsS();
1597  Double_t fshtpccls = 1.*ntpcclsS/ntpccls;
1598  QAbefore(19)->Fill(track->P(),fshtpccls);
1599  if (pass) QAafter(19)->Fill(track->P(),fshtpccls);
1600  }
1601  QAbefore(20)->Fill(track->P(),ntpccls);
1602  if (pass) QAafter(20)->Fill(track->P(),ntpccls);
1603  Int_t nitscls = track->GetITSNcls();
1604  if(nitscls>0) {
1605  Double_t chi2TIS = track->GetITSchi2()/track->GetITSNcls();
1606  QAbefore(21)->Fill(track->P(),chi2TIS);
1607  if (pass) QAafter(21)->Fill(track->P(),chi2TIS);
1608  }
1609  }
1610 
1611  return pass;
1612 }
1613 
1614 //_______________________________________________________________________
1615 Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
1616 {
1617  //check cuts on ESD tracks
1618  Bool_t pass=kTRUE;
1619  Float_t dcaxy=0.0;
1620  Float_t dcaz=0.0;
1621  track->GetImpactParameters(dcaxy,dcaz);
1622  const AliExternalTrackParam* pout = track->GetOuterParam();
1623  const AliExternalTrackParam* pin = track->GetInnerParam();
1624 
1625  if (fIgnoreTPCzRange)
1626  {
1627  if (pin&&pout)
1628  {
1629  Double_t zin = pin->GetZ();
1630  Double_t zout = pout->GetZ();
1631  if (zin*zout<0) pass=kFALSE; //reject if cross the membrane
1632  if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
1633  if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
1634  }
1635  }
1636 
1637  Int_t ntpccls = ( fParamType==kTPCstandalone )?
1638  track->GetTPCNclsIter1():track->GetTPCNcls();
1639  if (fCutChi2PerClusterTPC)
1640  {
1641  Float_t tpcchi2 = (fParamType==kTPCstandalone)?
1642  track->GetTPCchi2Iter1():track->GetTPCchi2();
1643  tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
1644  if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
1645  pass=kFALSE;
1646  }
1647 
1648  if (fCutMinimalTPCdedx)
1649  {
1650  if (track->GetTPCsignal() < fMinimalTPCdedx) pass=kFALSE;
1651  }
1652 
1653  if (fCutNClustersTPC)
1654  {
1655  if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
1656  }
1657 
1658  Int_t nitscls = track->GetNcls(0);
1659  if (fCutNClustersITS)
1660  {
1661  if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
1662  }
1663 
1664  //some stuff is still handled by AliESDtrackCuts class - delegate
1665  if (fAliESDtrackCuts)
1666  {
1667  if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
1668  }
1669 
1670  //PID part with pid QA
1671  Double_t beta = GetBeta(track);
1672  Double_t dedx = Getdedx(track);
1673  if (fQA)
1674  {
1675  if (pass) QAbefore(0)->Fill(track->GetP(),beta);
1676  if (pass) QAbefore(1)->Fill(pin->GetP(),dedx);
1677  }
1678  if (fCutPID && (fParticleID!=AliPID::kUnknown)) //if kUnknown don't cut on PID
1679  {
1680  if (!PassesESDpidCut(track)) pass=kFALSE;
1681  }
1682  if (fCutRejectElectronsWithTPCpid)
1683  {
1684  //reject electrons using standard TPC pid
1685  //TODO this should be rethought....
1686  Double_t pidTPC[AliPID::kSPECIES];
1687  track->GetTPCpid(pidTPC);
1688  if (pidTPC[AliPID::kElectron]<fParticleProbability) pass=kFALSE;
1689  }
1690  if(fCutITSclusterShared) {
1691  Int_t counterForSharedCluster=MaxSharedITSClusterCuts(track);
1692  if(counterForSharedCluster >= fMaxITSclusterShared) pass=kFALSE;
1693  }
1694 
1695  if(fCutITSChi2) {
1696  Double_t chi2perClusterITS = MaxChi2perITSClusterCuts(track);
1697  if(chi2perClusterITS >= fMaxITSChi2) pass=kFALSE;
1698  }
1699 
1700  if (fQA)
1701  {
1702  if (pass) QAafter(0)->Fill(track->GetP(),beta);
1703  if (pass) QAafter(1)->Fill(pin->GetP(),dedx);
1704  }
1705  //end pid part with qa
1706 
1707  if (fQA)
1708  {
1709  Double_t pt = track->Pt();
1710  QAbefore(5)->Fill(pt,dcaxy);
1711  QAbefore(6)->Fill(pt,dcaz);
1712  if (pass) QAafter(5)->Fill(pt,dcaxy);
1713  if (pass) QAafter(6)->Fill(pt,dcaz);
1714  QAbefore(17)->Fill(Float_t(track->GetKinkIndex(0)));
1715  if (pass) QAafter(17)->Fill(Float_t(track->GetKinkIndex(0)));
1716  }
1717 
1718  return pass;
1719 }
1720 
1721 //-----------------------------------------------------------------------
1722 void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
1723 {
1724  //handle the general case
1725  switch (fParamType)
1726  {
1727  default:
1728  fTrack = track;
1729  break;
1730  }
1731 }
1732 
1733 //-----------------------------------------------------------------------
1734 void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
1735 {
1736  //handle esd track
1737  switch (fParamType)
1738  {
1739  case kGlobal:
1740  fTrack = track;
1741  break;
1742  case kTPCstandalone:
1743  if (!track->FillTPCOnlyTrack(fTPCtrack))
1744  {
1745  fTrack=NULL;
1746  fMCparticle=NULL;
1747  fTrackLabel=-998;
1748  return;
1749  }
1750  fTrack = &fTPCtrack;
1751  //recalculate the label and mc particle, they may differ as TPClabel != global label
1752  fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1753  if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1754  else fMCparticle=NULL;
1755  break;
1756  default:
1757  if (fForceTPCstandalone)
1758  {
1759  if (!track->FillTPCOnlyTrack(fTPCtrack))
1760  {
1761  fTrack=NULL;
1762  fMCparticle=NULL;
1763  fTrackLabel=-998;
1764  return;
1765  }
1766  fTrack = &fTPCtrack;
1767  //recalculate the label and mc particle, they may differ as TPClabel != global label
1768  fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
1769  if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
1770  else fMCparticle=NULL;
1771  }
1772  else
1773  fTrack = track;
1774  break;
1775  }
1776 }
1777 
1778 //-----------------------------------------------------------------------
1779 Int_t AliFlowTrackCuts::Count(AliVEvent* event)
1780 {
1781  //calculate the number of track in given event.
1782  //if argument is NULL(default) take the event attached
1783  //by SetEvent()
1784 
1785  Int_t multiplicity = 0;
1786  if (!event)
1787  {
1788  for (Int_t i=0; i<GetNumberOfInputObjects(); i++)
1789  {
1790  if (IsSelected(GetInputObject(i))) multiplicity++;
1791  }
1792  }
1793  else
1794  {
1795  for (Int_t i=0; i<event->GetNumberOfTracks(); i++)
1796  {
1797  if (IsSelected(event->GetTrack(i))) multiplicity++;
1798  }
1799  }
1800  return multiplicity;
1801 }
1802 
1803 //-----------------------------------------------------------------------
1804 AliFlowTrackCuts* AliFlowTrackCuts::GetAODTrackCutsForFilterBit(UInt_t bit, TString suffix)
1805 {
1806  // object which in its default form only cuts on filterbit (for AOD analysis)
1807  // NOTE : the cut object is defined in such a way that ONLY filterbit is tested,
1808  // all other paramters are NOT checked
1809  // the exception is TPCdecx which is always cut for reasons of backward compatibility
1810  // but by setting the threshold to larger than -99999999 effectively there is no check
1811  AliFlowTrackCuts* cuts = new AliFlowTrackCuts(Form("AOD fitlerbit %i, %s", (int)bit, suffix.Data()));
1812  cuts->SetMinimalTPCdedx(-999999999);
1813  cuts->SetAODfilterBit(bit);
1814  cuts->SetParamType(AliFlowTrackCuts::kAODFilterBit);
1815  return cuts;
1816 }
1817 //-----------------------------------------------------------------------
1818 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts()
1819 {
1820  //returns the lhc10h vzero track cuts, this function
1821  //is left here for backward compatibility
1822  //if a run is recognized as 11h, the calibration method will
1823  //switch to 11h calbiration, which means that the cut
1824  //object is updated but not replaced.
1825  //calibratin is only available for PbPb runs
1826  return GetStandardVZEROOnlyTrackCuts2010();
1827 }
1828 //-----------------------------------------------------------------------
1829 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2010()
1830 {
1831  //get standard VZERO cuts
1832  //DISCLAIMER: LHC10h VZERO calibration consists (by default) of two steps
1833  //1) re-weigting of signal
1834  //2) re-centering of q-vectors
1835  //step 2 is available only for n==2 and n==3, for the higher harmonics the user
1836  //is repsonsible for making sure the q-sub distributions are (sufficiently) flat
1837  //or a sensible NUA procedure is applied !
1838  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts");
1839  cuts->SetParamType(AliFlowTrackCuts::kVZERO);
1840  cuts->SetEtaRange( -10, +10 );
1841  cuts->SetEtaGap(-1., 1.);
1842  cuts->SetPhiMin( 0 );
1843  cuts->SetPhiMax( TMath::TwoPi() );
1844  // options for the reweighting
1845  cuts->SetVZEROgainEqualizationPerRing(kFALSE);
1846  cuts->SetApplyRecentering(kTRUE);
1847  // to exclude a ring , do e.g.
1848  // cuts->SetUseVZERORing(7, kFALSE);
1849  // excluding a ring will break the re-centering as re-centering relies on a
1850  // database file which tuned to receiving info from all rings
1851  return cuts;
1852 }
1853 //-----------------------------------------------------------------------
1854 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts2011()
1855 {
1856  //get standard VZERO cuts for 2011 data
1857  //in this case, the vzero segments will be weighted by
1858  //VZEROEqMultiplicity,
1859  //if recentering is enableded, the sub-q vectors
1860  //will be taken from the event header, so make sure to run
1861  //the VZERO event plane selection task before this task !
1862  //DISCLAIMER: recentering is only available for n==2
1863  //for the higher harmonics the user
1864  //is repsonsible for making sure the q-sub distributions are (sufficiently) flat
1865  //or a sensible NUA procedure is applied !
1866  //recentering replaces the already evaluated q-vectors, so
1867  //when chosen, additional settings (e.g. excluding rings)
1868  //have no effect. recentering is true by default
1869  //
1870  //NOTE user is responsible for running the vzero event plane
1871  //selection task in advance, e.g. add to your launcher macro
1872  //
1873  // gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskVZEROEPSelection.C");
1874  // AddTaskVZEROEPSelection();
1875  //
1876  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard vzero flow cuts 2011");
1877  cuts->SetParamType(kVZERO);
1878  cuts->SetEtaRange( -10, +10 );
1879  cuts->SetEtaGap(-1., 1.);
1880  cuts->SetPhiMin( 0 );
1881  cuts->SetPhiMax( TMath::TwoPi() );
1882  cuts->SetApplyRecentering(kTRUE);
1883  cuts->SetVZEROgainEqualizationPerRing(kFALSE);
1884  return cuts;
1885 }
1886 //-----------------------------------------------------------------------
1887 AliFlowTrackCuts* AliFlowTrackCuts::GetBetaVZEROOnlyTrackCuts()
1888 {
1889  // test patch for VZERO track cuts
1890  // april 20 2015.
1891  // DISCLAIMER: BETA TEST OF CODE, DO NOT USE FOR PHYSICS RESULTS !!!
1892  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("beta version of vzero ");
1893  cuts->SetParamType(AliFlowTrackCuts::kBetaVZERO);
1894  cuts->SetEtaRange( -10, +10 );
1895  cuts->SetEtaGap(-1., 1.);
1896  cuts->SetPhiMin( 0 );
1897  cuts->SetPhiMax( TMath::TwoPi() );
1898  cuts->SetApplyRecentering(kTRUE);
1899  // idea of the cuts is that calibration is done per ring
1900  // and that it is transparent for different data taking periods
1901  return cuts;
1902 }
1903 
1904 
1905 //-----------------------------------------------------------------------
1906 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardGlobalTrackCuts2010()
1907 {
1908  //get standard cuts
1909  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard Global tracks");
1910  cuts->SetParamType(kGlobal);
1911  cuts->SetPtRange(0.2,5.);
1912  cuts->SetEtaRange(-0.8,0.8);
1913  cuts->SetMinNClustersTPC(70);
1914  cuts->SetMinChi2PerClusterTPC(0.1);
1915  cuts->SetMaxChi2PerClusterTPC(4.0);
1916  cuts->SetMinNClustersITS(2);
1917  cuts->SetRequireITSRefit(kTRUE);
1918  cuts->SetRequireTPCRefit(kTRUE);
1919  cuts->SetMaxDCAToVertexXY(0.3);
1920  cuts->SetMaxDCAToVertexZ(0.3);
1921  cuts->SetAcceptKinkDaughters(kFALSE);
1922  cuts->SetMinimalTPCdedx(10.);
1923 
1924  return cuts;
1925 }
1926 
1927 //-----------------------------------------------------------------------
1928 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts2010()
1929 {
1930  //get standard cuts
1931  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone 2010");
1932  cuts->SetParamType(kTPCstandalone);
1933  cuts->SetPtRange(0.2,5.);
1934  cuts->SetEtaRange(-0.8,0.8);
1935  cuts->SetMinNClustersTPC(70);
1936  cuts->SetMinChi2PerClusterTPC(0.2);
1937  cuts->SetMaxChi2PerClusterTPC(4.0);
1938  cuts->SetMaxDCAToVertexXY(3.0);
1939  cuts->SetMaxDCAToVertexZ(3.0);
1940  cuts->SetDCAToVertex2D(kTRUE);
1941  cuts->SetAcceptKinkDaughters(kFALSE);
1942  cuts->SetMinimalTPCdedx(10.);
1943 
1944  return cuts;
1945 }
1946 
1947 //-----------------------------------------------------------------------
1948 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCStandaloneTrackCuts()
1949 {
1950  //get standard cuts
1951  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPC standalone");
1952  cuts->SetParamType(kTPCstandalone);
1953  cuts->SetPtRange(0.2,5.);
1954  cuts->SetEtaRange(-0.8,0.8);
1955  cuts->SetMinNClustersTPC(70);
1956  cuts->SetMinChi2PerClusterTPC(0.2);
1957  cuts->SetMaxChi2PerClusterTPC(4.0);
1958  cuts->SetMaxDCAToVertexXY(3.0);
1959  cuts->SetMaxDCAToVertexZ(3.0);
1960  cuts->SetDCAToVertex2D(kTRUE);
1961  cuts->SetAcceptKinkDaughters(kFALSE);
1962  cuts->SetMinimalTPCdedx(10.);
1963 
1964  return cuts;
1965 }
1966 
1967 //-----------------------------------------------------------------------
1968 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
1969 {
1970  //get standard cuts
1971  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
1972  delete cuts->fAliESDtrackCuts;
1973  cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
1974  cuts->SetParamType(kGlobal);
1975  return cuts;
1976 }
1977 
1978 //-----------------------------------------------------------------------------
1979 AliFlowTrackCuts* AliFlowTrackCuts::GetStandardMuonTrackCuts(Bool_t isMC, Int_t passN)
1980 {
1981 // XZhang 20120604
1982  AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard muon track cuts");
1983  cuts->SetParamType(kMUON);
1984  cuts->SetStandardMuonTrackCuts();
1985  cuts->SetIsMuonMC(isMC);
1986  cuts->SetMuonPassNumber(passN);
1987  return cuts;
1988 }
1989 
1990 //-----------------------------------------------------------------------
1991 //AliFlowTrack* AliFlowTrackCuts::FillFlowTrackV0(TObjArray* trackCollection, Int_t trackIndex) const
1992 //{
1993 // //fill flow track from a reconstructed V0 (topological)
1994 // AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
1995 // if (flowtrack)
1996 // {
1997 // flowtrack->Clear();
1998 // }
1999 // else
2000 // {
2001 // flowtrack = new AliFlowCandidateTrack();
2002 // trackCollection->AddAtAndExpand(flowtrack,trackIndex);
2003 // }
2004 //
2005 // TParticle *tmpTParticle=NULL;
2006 // AliMCParticle* tmpAliMCParticle=NULL;
2007 // AliExternalTrackParam* externalParams=NULL;
2008 // AliESDtrack* esdtrack=NULL;
2009 // switch(fParamMix)
2010 // {
2011 // case kPure:
2012 // flowtrack->Set(fTrack);
2013 // break;
2014 // case kTrackWithMCkine:
2015 // flowtrack->Set(fMCparticle);
2016 // break;
2017 // case kTrackWithMCPID:
2018 // flowtrack->Set(fTrack);
2019 // //flowtrack->setPID(...) from mc, when implemented
2020 // break;
2021 // case kTrackWithMCpt:
2022 // if (!fMCparticle) return NULL;
2023 // flowtrack->Set(fTrack);
2024 // flowtrack->SetPt(fMCparticle->Pt());
2025 // break;
2026 // case kTrackWithPtFromFirstMother:
2027 // if (!fMCparticle) return NULL;
2028 // flowtrack->Set(fTrack);
2029 // tmpTParticle = fMCparticle->Particle();
2030 // tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2031 // flowtrack->SetPt(tmpAliMCParticle->Pt());
2032 // break;
2033 // case kTrackWithTPCInnerParams:
2034 // esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2035 // if (!esdtrack) return NULL;
2036 // externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2037 // if (!externalParams) return NULL;
2038 // flowtrack->Set(externalParams);
2039 // break;
2040 // default:
2041 // flowtrack->Set(fTrack);
2042 // break;
2043 // }
2044 // if (fParamType==kMC)
2045 // {
2046 // flowtrack->SetSource(AliFlowTrack::kFromMC);
2047 // flowtrack->SetID(fTrackLabel);
2048 // }
2049 // else if (dynamic_cast<AliAODTrack*>(fTrack))
2050 // {
2051 // if (fParamType==kMUON) // XZhang 20120604
2052 // flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
2053 // else // XZhang 20120604
2054 // flowtrack->SetSource(AliFlowTrack::kFromAOD); // XZhang 20120604
2055 // flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2056 // }
2057 // else if (dynamic_cast<AliMCParticle*>(fTrack))
2058 // {
2059 // flowtrack->SetSource(AliFlowTrack::kFromMC);
2060 // flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2061 // }
2062 //
2063 // if (fV0)
2064 // {
2065 // //add daughter indices
2066 // }
2067 //
2068 // return flowtrack;
2069 //}
2070 
2071 //-----------------------------------------------------------------------
2072 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackKink(TObjArray* trackCollection, Int_t trackIndex) const
2073 {
2074  //fill flow track from AliVParticle (ESD,AOD,MC)
2075  AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
2076  if (flowtrack)
2077  {
2078  flowtrack->Clear();
2079  }
2080  else
2081  {
2082  flowtrack = new AliFlowCandidateTrack();
2083  trackCollection->AddAtAndExpand(flowtrack,trackIndex);
2084  }
2085 
2086  TParticle *tmpTParticle=NULL;
2087  AliMCParticle* tmpAliMCParticle=NULL;
2088  AliExternalTrackParam* externalParams=NULL;
2089  AliESDtrack* esdtrack=NULL;
2090  switch(fParamMix)
2091  {
2092  case kPure:
2093  flowtrack->Set(fTrack);
2094  break;
2095  case kTrackWithMCkine:
2096  flowtrack->Set(fMCparticle);
2097  break;
2098  case kTrackWithMCPID:
2099  flowtrack->Set(fTrack);
2100  //flowtrack->setPID(...) from mc, when implemented
2101  break;
2102  case kTrackWithMCpt:
2103  if (!fMCparticle) return NULL;
2104  flowtrack->Set(fTrack);
2105  flowtrack->SetPt(fMCparticle->Pt());
2106  break;
2107  case kTrackWithPtFromFirstMother:
2108  if (!fMCparticle) return NULL;
2109  flowtrack->Set(fTrack);
2110  tmpTParticle = fMCparticle->Particle();
2111  tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2112  flowtrack->SetPt(tmpAliMCParticle->Pt());
2113  break;
2114  case kTrackWithTPCInnerParams:
2115  esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2116  if (!esdtrack) return NULL;
2117  externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2118  if (!externalParams) return NULL;
2119  flowtrack->Set(externalParams);
2120  break;
2121  default:
2122  flowtrack->Set(fTrack);
2123  break;
2124  }
2125  if (fParamType==kMC)
2126  {
2127  flowtrack->SetSource(AliFlowTrack::kFromMC);
2128  flowtrack->SetID(fTrackLabel);
2129  }
2130  else if (dynamic_cast<AliESDtrack*>(fTrack))
2131  {
2132  flowtrack->SetSource(AliFlowTrack::kFromESD);
2133  flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2134  }
2135  else if (dynamic_cast<AliESDMuonTrack*>(fTrack)) // XZhang 20120604
2136  { // XZhang 20120604
2137  flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
2138  flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID()); // XZhang 20120604
2139  } // XZhang 20120604
2140  else if (dynamic_cast<AliAODTrack*>(fTrack))
2141  {
2142  if (fParamType==kMUON) // XZhang 20120604
2143  flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
2144  else // XZhang 20120604
2145  flowtrack->SetSource(AliFlowTrack::kFromAOD); // XZhang 20120604
2146  flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2147  }
2148  else if (dynamic_cast<AliMCParticle*>(fTrack))
2149  {
2150  flowtrack->SetSource(AliFlowTrack::kFromMC);
2151  flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2152  }
2153 
2154  if (fKink)
2155  {
2156  Int_t indexMother = fKink->GetIndex(0);
2157  Int_t indexDaughter = fKink->GetIndex(1);
2158  flowtrack->SetID(indexMother);
2159  flowtrack->AddDaughter(indexDaughter);
2160  flowtrack->SetMass(1.);
2161  flowtrack->SetSource(AliFlowTrack::kFromKink);
2162  }
2163 
2164  flowtrack->SetMass(fTrackMass);
2165 
2166  return flowtrack;
2167 }
2168 
2169 //-----------------------------------------------------------------------
2170 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackGeneric(TObjArray* trackCollection, Int_t trackIndex) const
2171 {
2172  //fill a flow track from tracklet,vzero,pmd,...
2173  AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
2174  if (flowtrack)
2175  {
2176  flowtrack->Clear();
2177  }
2178  else
2179  {
2180  flowtrack = new AliFlowTrack();
2181  trackCollection->AddAtAndExpand(flowtrack,trackIndex);
2182  }
2183 
2184  if (FillFlowTrackGeneric(flowtrack)) return flowtrack;
2185  else
2186  {
2187  trackCollection->RemoveAt(trackIndex);
2188  delete flowtrack;
2189  return NULL;
2190  }
2191 }
2192 
2193 //-----------------------------------------------------------------------
2194 Bool_t AliFlowTrackCuts::FillFlowTrackGeneric(AliFlowTrack* flowtrack) const
2195 {
2196  //fill a flow track from tracklet,vzero,pmd,...
2197  TParticle *tmpTParticle=NULL;
2198  AliMCParticle* tmpAliMCParticle=NULL;
2199  switch (fParamMix)
2200  {
2201  case kPure:
2202  flowtrack->SetPhi(fTrackPhi);
2203  flowtrack->SetEta(fTrackEta);
2204  flowtrack->SetWeight(fTrackWeight);
2205  flowtrack->SetPt(fTrackPt);
2206  flowtrack->SetMass(fTrackMass);
2207  break;
2208  case kTrackWithMCkine:
2209  if (!fMCparticle) return kFALSE;
2210  flowtrack->SetPhi( fMCparticle->Phi() );
2211  flowtrack->SetEta( fMCparticle->Eta() );
2212  flowtrack->SetPt( fMCparticle->Pt() );
2213  flowtrack->SetMass(fTrackMass);
2214  break;
2215  case kTrackWithMCpt:
2216  if (!fMCparticle) return kFALSE;
2217  flowtrack->SetPhi(fTrackPhi);
2218  flowtrack->SetEta(fTrackEta);
2219  flowtrack->SetWeight(fTrackWeight);
2220  flowtrack->SetPt(fMCparticle->Pt());
2221  flowtrack->SetMass(fTrackMass);
2222  break;
2223  case kTrackWithPtFromFirstMother:
2224  if (!fMCparticle) return kFALSE;
2225  flowtrack->SetPhi(fTrackPhi);
2226  flowtrack->SetEta(fTrackEta);
2227  flowtrack->SetWeight(fTrackWeight);
2228  tmpTParticle = fMCparticle->Particle();
2229  tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2230  flowtrack->SetPt(tmpAliMCParticle->Pt());
2231  flowtrack->SetMass(fTrackMass);
2232  break;
2233  default:
2234  flowtrack->SetPhi(fTrackPhi);
2235  flowtrack->SetEta(fTrackEta);
2236  flowtrack->SetWeight(fTrackWeight);
2237  flowtrack->SetPt(fTrackPt);
2238  flowtrack->SetMass(fTrackMass);
2239  break;
2240  }
2241  flowtrack->SetSource(AliFlowTrack::kFromTracklet);
2242  return kTRUE;
2243 }
2244 
2245 //-----------------------------------------------------------------------
2246 AliFlowTrack* AliFlowTrackCuts::FillFlowTrackVParticle(TObjArray* trackCollection, Int_t trackIndex) const
2247 {
2248  //fill flow track from AliVParticle (ESD,AOD,MC)
2249 
2250  AliFlowTrack* flowtrack = static_cast<AliFlowTrack*>(trackCollection->At(trackIndex));
2251  if (flowtrack)
2252  {
2253  flowtrack->Clear();
2254  }
2255  else
2256  {
2257  flowtrack = new AliFlowTrack();
2258  trackCollection->AddAtAndExpand(flowtrack,trackIndex);
2259  }
2260 
2261  if (FillFlowTrackVParticle(flowtrack)) return flowtrack;
2262  else
2263  {
2264  trackCollection->RemoveAt(trackIndex);
2265  delete flowtrack;
2266  return NULL;
2267  }
2268 }
2269 
2270 //-----------------------------------------------------------------------
2271 Bool_t AliFlowTrackCuts::FillFlowTrackVParticle(AliFlowTrack* flowtrack) const
2272 {
2273  //fill flow track from AliVParticle (ESD,AOD,MC)
2274  if (!fTrack) return kFALSE;
2275  TParticle *tmpTParticle=NULL;
2276  AliMCParticle* tmpAliMCParticle=NULL;
2277  AliExternalTrackParam* externalParams=NULL;
2278  AliESDtrack* esdtrack=NULL;
2279  switch(fParamMix)
2280  {
2281  case kPure:
2282  flowtrack->Set(fTrack);
2283  break;
2284  case kTrackWithMCkine:
2285  flowtrack->Set(fMCparticle);
2286  break;
2287  case kTrackWithMCPID:
2288  flowtrack->Set(fTrack);
2289  //flowtrack->setPID(...) from mc, when implemented
2290  break;
2291  case kTrackWithMCpt:
2292  if (!fMCparticle) return kFALSE;
2293  flowtrack->Set(fTrack);
2294  flowtrack->SetPt(fMCparticle->Pt());
2295  break;
2296  case kTrackWithPtFromFirstMother:
2297  if (!fMCparticle) return kFALSE;
2298  flowtrack->Set(fTrack);
2299  tmpTParticle = fMCparticle->Particle();
2300  tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2301  flowtrack->SetPt(tmpAliMCParticle->Pt());
2302  break;
2303  case kTrackWithTPCInnerParams:
2304  esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2305  if (!esdtrack) return kFALSE;
2306  externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2307  if (!externalParams) return kFALSE;
2308  flowtrack->Set(externalParams);
2309  break;
2310  default:
2311  flowtrack->Set(fTrack);
2312  break;
2313  }
2314  if (fParamType==kMC)
2315  {
2316  flowtrack->SetSource(AliFlowTrack::kFromMC);
2317  flowtrack->SetID(fTrackLabel);
2318  }
2319  else if (dynamic_cast<AliESDtrack*>(fTrack))
2320  {
2321  flowtrack->SetSource(AliFlowTrack::kFromESD);
2322  flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2323  }
2324  else if (dynamic_cast<AliESDMuonTrack*>(fTrack)) // XZhang 20120604
2325  { // XZhang 20120604
2326  flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
2327  flowtrack->SetID((Int_t)static_cast<AliESDMuonTrack*>(fTrack)->GetUniqueID()); // XZhang 20120604
2328  } // XZhang 20120604
2329  else if (dynamic_cast<AliAODTrack*>(fTrack))
2330  {
2331  if (fParamType==kMUON) // XZhang 20120604
2332  flowtrack->SetSource(AliFlowTrack::kFromMUON); // XZhang 20120604
2333  else // XZhang 20120604
2334  flowtrack->SetSource(AliFlowTrack::kFromAOD); // XZhang 20120604
2335  flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2336  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fTrack);
2337  flowtrack->SetITStype(GetITStype(aodTrack));
2338  }
2339  else if (dynamic_cast<AliMCParticle*>(fTrack))
2340  {
2341  flowtrack->SetSource(AliFlowTrack::kFromMC);
2342  flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2343  }
2344  flowtrack->SetMass(fTrackMass);
2345  return kTRUE;
2346 }
2347 
2348 //-----------------------------------------------------------------------
2349 AliFlowTrack* AliFlowTrackCuts::FillFlowTrack(TObjArray* trackCollection, Int_t trackIndex) const
2350 {
2351  //fill a flow track constructed from whatever we applied cuts on
2352  //return true on success
2353  switch (fParamType)
2354  {
2355  case kSPDtracklet:
2356  return FillFlowTrackGeneric(trackCollection, trackIndex);
2357  case kPMD:
2358  return FillFlowTrackGeneric(trackCollection, trackIndex);
2359  case kVZERO:
2360  return FillFlowTrackGeneric(trackCollection, trackIndex);
2361  case kBetaVZERO:
2362  return FillFlowTrackGeneric(trackCollection, trackIndex);
2363  case kDeltaVZERO:
2364  return FillFlowTrackGeneric(trackCollection, trackIndex);
2365  case kKappaVZERO:
2366  return FillFlowTrackGeneric(trackCollection, trackIndex);
2367  case kHotfixHI:
2368  return FillFlowTrackGeneric(trackCollection, trackIndex);
2369  case kKink:
2370  return FillFlowTrackKink(trackCollection, trackIndex);
2371  //case kV0:
2372  // return FillFlowTrackV0(trackCollection, trackIndex);
2373  default:
2374  return FillFlowTrackVParticle(trackCollection, trackIndex);
2375  }
2376 }
2377 
2378 //-----------------------------------------------------------------------
2379 Bool_t AliFlowTrackCuts::FillFlowTrack(AliFlowTrack* track) const
2380 {
2381  //fill a flow track constructed from whatever we applied cuts on
2382  //return true on success
2383  switch (fParamType)
2384  {
2385  case kSPDtracklet:
2386  return FillFlowTrackGeneric(track);
2387  case kPMD:
2388  return FillFlowTrackGeneric(track);
2389  case kVZERO:
2390  return FillFlowTrackGeneric(track);
2391  case kBetaVZERO:
2392  return FillFlowTrackGeneric(track);
2393  case kDeltaVZERO:
2394  return FillFlowTrackGeneric(track);
2395  case kKappaVZERO:
2396  return FillFlowTrackGeneric(track);
2397  case kHotfixHI:
2398  return FillFlowTrackGeneric(track);
2399  default:
2400  return FillFlowTrackVParticle(track);
2401  }
2402 }
2403 
2405 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackSPDtracklet() const
2406 //{
2407 // //make a flow track from tracklet
2408 // AliFlowTrack* flowtrack=NULL;
2409 // TParticle *tmpTParticle=NULL;
2410 // AliMCParticle* tmpAliMCParticle=NULL;
2411 // switch (fParamMix)
2412 // {
2413 // case kPure:
2414 // flowtrack = new AliFlowTrack();
2415 // flowtrack->SetPhi(fTrackPhi);
2416 // flowtrack->SetEta(fTrackEta);
2417 // flowtrack->SetWeight(fTrackWeight);
2418 // break;
2419 // case kTrackWithMCkine:
2420 // if (!fMCparticle) return NULL;
2421 // flowtrack = new AliFlowTrack();
2422 // flowtrack->SetPhi( fMCparticle->Phi() );
2423 // flowtrack->SetEta( fMCparticle->Eta() );
2424 // flowtrack->SetPt( fMCparticle->Pt() );
2425 // break;
2426 // case kTrackWithMCpt:
2427 // if (!fMCparticle) return NULL;
2428 // flowtrack = new AliFlowTrack();
2429 // flowtrack->SetPhi(fTrackPhi);
2430 // flowtrack->SetEta(fTrackEta);
2431 // flowtrack->SetWeight(fTrackWeight);
2432 // flowtrack->SetPt(fMCparticle->Pt());
2433 // break;
2434 // case kTrackWithPtFromFirstMother:
2435 // if (!fMCparticle) return NULL;
2436 // flowtrack = new AliFlowTrack();
2437 // flowtrack->SetPhi(fTrackPhi);
2438 // flowtrack->SetEta(fTrackEta);
2439 // flowtrack->SetWeight(fTrackWeight);
2440 // tmpTParticle = fMCparticle->Particle();
2441 // tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2442 // flowtrack->SetPt(tmpAliMCParticle->Pt());
2443 // break;
2444 // default:
2445 // flowtrack = new AliFlowTrack();
2446 // flowtrack->SetPhi(fTrackPhi);
2447 // flowtrack->SetEta(fTrackEta);
2448 // flowtrack->SetWeight(fTrackWeight);
2449 // break;
2450 // }
2451 // flowtrack->SetSource(AliFlowTrack::kFromTracklet);
2452 // return flowtrack;
2453 //}
2454 //
2456 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVParticle() const
2457 //{
2458 // //make flow track from AliVParticle (ESD,AOD,MC)
2459 // if (!fTrack) return NULL;
2460 // AliFlowTrack* flowtrack=NULL;
2461 // TParticle *tmpTParticle=NULL;
2462 // AliMCParticle* tmpAliMCParticle=NULL;
2463 // AliExternalTrackParam* externalParams=NULL;
2464 // AliESDtrack* esdtrack=NULL;
2465 // switch(fParamMix)
2466 // {
2467 // case kPure:
2468 // flowtrack = new AliFlowTrack(fTrack);
2469 // break;
2470 // case kTrackWithMCkine:
2471 // flowtrack = new AliFlowTrack(fMCparticle);
2472 // break;
2473 // case kTrackWithMCPID:
2474 // flowtrack = new AliFlowTrack(fTrack);
2475 // //flowtrack->setPID(...) from mc, when implemented
2476 // break;
2477 // case kTrackWithMCpt:
2478 // if (!fMCparticle) return NULL;
2479 // flowtrack = new AliFlowTrack(fTrack);
2480 // flowtrack->SetPt(fMCparticle->Pt());
2481 // break;
2482 // case kTrackWithPtFromFirstMother:
2483 // if (!fMCparticle) return NULL;
2484 // flowtrack = new AliFlowTrack(fTrack);
2485 // tmpTParticle = fMCparticle->Particle();
2486 // tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2487 // flowtrack->SetPt(tmpAliMCParticle->Pt());
2488 // break;
2489 // case kTrackWithTPCInnerParams:
2490 // esdtrack = dynamic_cast<AliESDtrack*>(fTrack);
2491 // if (!esdtrack) return NULL;
2492 // externalParams = const_cast<AliExternalTrackParam*>(esdtrack->GetTPCInnerParam());
2493 // if (!externalParams) return NULL;
2494 // flowtrack = new AliFlowTrack(externalParams);
2495 // break;
2496 // default:
2497 // flowtrack = new AliFlowTrack(fTrack);
2498 // break;
2499 // }
2500 // if (fParamType==kMC)
2501 // {
2502 // flowtrack->SetSource(AliFlowTrack::kFromMC);
2503 // flowtrack->SetID(fTrackLabel);
2504 // }
2505 // else if (dynamic_cast<AliESDtrack*>(fTrack))
2506 // {
2507 // flowtrack->SetSource(AliFlowTrack::kFromESD);
2508 // flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2509 // }
2510 // else if (dynamic_cast<AliAODTrack*>(fTrack))
2511 // {
2512 // flowtrack->SetSource(AliFlowTrack::kFromAOD);
2513 // flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2514 // }
2515 // else if (dynamic_cast<AliMCParticle*>(fTrack))
2516 // {
2517 // flowtrack->SetSource(AliFlowTrack::kFromMC);
2518 // flowtrack->SetID(static_cast<AliVTrack*>(fTrack)->GetID());
2519 // }
2520 // return flowtrack;
2521 //}
2522 //
2524 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackPMDtrack() const
2525 //{
2526 // //make a flow track from PMD track
2527 // AliFlowTrack* flowtrack=NULL;
2528 // TParticle *tmpTParticle=NULL;
2529 // AliMCParticle* tmpAliMCParticle=NULL;
2530 // switch (fParamMix)
2531 // {
2532 // case kPure:
2533 // flowtrack = new AliFlowTrack();
2534 // flowtrack->SetPhi(fTrackPhi);
2535 // flowtrack->SetEta(fTrackEta);
2536 // flowtrack->SetWeight(fTrackWeight);
2537 // break;
2538 // case kTrackWithMCkine:
2539 // if (!fMCparticle) return NULL;
2540 // flowtrack = new AliFlowTrack();
2541 // flowtrack->SetPhi( fMCparticle->Phi() );
2542 // flowtrack->SetEta( fMCparticle->Eta() );
2543 // flowtrack->SetWeight(fTrackWeight);
2544 // flowtrack->SetPt( fMCparticle->Pt() );
2545 // break;
2546 // case kTrackWithMCpt:
2547 // if (!fMCparticle) return NULL;
2548 // flowtrack = new AliFlowTrack();
2549 // flowtrack->SetPhi(fTrackPhi);
2550 // flowtrack->SetEta(fTrackEta);
2551 // flowtrack->SetWeight(fTrackWeight);
2552 // flowtrack->SetPt(fMCparticle->Pt());
2553 // break;
2554 // case kTrackWithPtFromFirstMother:
2555 // if (!fMCparticle) return NULL;
2556 // flowtrack = new AliFlowTrack();
2557 // flowtrack->SetPhi(fTrackPhi);
2558 // flowtrack->SetEta(fTrackEta);
2559 // flowtrack->SetWeight(fTrackWeight);
2560 // tmpTParticle = fMCparticle->Particle();
2561 // tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2562 // flowtrack->SetPt(tmpAliMCParticle->Pt());
2563 // break;
2564 // default:
2565 // flowtrack = new AliFlowTrack();
2566 // flowtrack->SetPhi(fTrackPhi);
2567 // flowtrack->SetEta(fTrackEta);
2568 // flowtrack->SetWeight(fTrackWeight);
2569 // break;
2570 // }
2571 //
2572 // flowtrack->SetSource(AliFlowTrack::kFromPMD);
2573 // return flowtrack;
2574 //}
2575 //
2577 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrackVZERO() const
2578 //{
2579 // //make a flow track from VZERO
2580 // AliFlowTrack* flowtrack=NULL;
2581 // TParticle *tmpTParticle=NULL;
2582 // AliMCParticle* tmpAliMCParticle=NULL;
2583 // switch (fParamMix)
2584 // {
2585 // case kPure:
2586 // flowtrack = new AliFlowTrack();
2587 // flowtrack->SetPhi(fTrackPhi);
2588 // flowtrack->SetEta(fTrackEta);
2589 // flowtrack->SetWeight(fTrackWeight);
2590 // break;
2591 // case kTrackWithMCkine:
2592 // if (!fMCparticle) return NULL;
2593 // flowtrack = new AliFlowTrack();
2594 // flowtrack->SetPhi( fMCparticle->Phi() );
2595 // flowtrack->SetEta( fMCparticle->Eta() );
2596 // flowtrack->SetWeight(fTrackWeight);
2597 // flowtrack->SetPt( fMCparticle->Pt() );
2598 // break;
2599 // case kTrackWithMCpt:
2600 // if (!fMCparticle) return NULL;
2601 // flowtrack = new AliFlowTrack();
2602 // flowtrack->SetPhi(fTrackPhi);
2603 // flowtrack->SetEta(fTrackEta);
2604 // flowtrack->SetWeight(fTrackWeight);
2605 // flowtrack->SetPt(fMCparticle->Pt());
2606 // break;
2607 // case kTrackWithPtFromFirstMother:
2608 // if (!fMCparticle) return NULL;
2609 // flowtrack = new AliFlowTrack();
2610 // flowtrack->SetPhi(fTrackPhi);
2611 // flowtrack->SetEta(fTrackEta);
2612 // flowtrack->SetWeight(fTrackWeight);
2613 // tmpTParticle = fMCparticle->Particle();
2614 // tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
2615 // flowtrack->SetPt(tmpAliMCParticle->Pt());
2616 // break;
2617 // default:
2618 // flowtrack = new AliFlowTrack();
2619 // flowtrack->SetPhi(fTrackPhi);
2620 // flowtrack->SetEta(fTrackEta);
2621 // flowtrack->SetWeight(fTrackWeight);
2622 // break;
2623 // }
2624 //
2625 // flowtrack->SetSource(AliFlowTrack::kFromVZERO);
2626 // return flowtrack;
2627 //}
2628 //
2630 //AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
2631 //{
2632 // //get a flow track constructed from whatever we applied cuts on
2633 // //caller is resposible for deletion
2634 // //if construction fails return NULL
2635 // //TODO: for tracklets, PMD and VZERO we probably need just one method,
2636 // //something like MakeFlowTrackGeneric(), wait with this until
2637 // //requirements quirks are known.
2638 // switch (fParamType)
2639 // {
2640 // case kSPDtracklet:
2641 // return MakeFlowTrackSPDtracklet();
2642 // case kPMD:
2643 // return MakeFlowTrackPMDtrack();
2644 // case kVZERO:
2645 // return MakeFlowTrackVZERO();
2646 // default:
2647 // return MakeFlowTrackVParticle();
2648 // }
2649 //}
2650 
2651 //-----------------------------------------------------------------------
2652 Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
2653 {
2654  //check if current particle is a physical primary
2655  if (!fMCevent) return kFALSE;
2656  if (fTrackLabel<0) return kFALSE;
2657  return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
2658 }
2659 
2660 //-----------------------------------------------------------------------
2661 Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
2662 {
2663  //check if current particle is a physical primary
2664  Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
2665  AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
2666  if (!track) return kFALSE;
2667  TParticle* particle = track->Particle();
2668  Bool_t transported = particle->TestBit(kTransportBit);
2669  //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
2670  //(physprim && (transported || !requiretransported))?"YES":"NO" );
2671  return (physprim && (transported || !requiretransported));
2672 }
2673 
2674 //-----------------------------------------------------------------------
2675 void AliFlowTrackCuts::DefineHistograms()
2676 {
2677  //define qa histograms
2678  if (fQA) return;
2679 
2680  const Int_t kNbinsP=200;
2681  Double_t binsP[kNbinsP+1];
2682  binsP[0]=0.0;
2683  for(int i=1; i<kNbinsP+1; i++)
2684  {
2685  //if(binsP[i-1]+0.05<1.01)
2686  // binsP[i]=binsP[i-1]+0.05;
2687  //else
2688  binsP[i]=binsP[i-1]+0.05;
2689  }
2690 
2691  const Int_t nBinsDCA=1000;
2692  Double_t binsDCA[nBinsDCA+1];
2693  for(int i=0; i<nBinsDCA+1; i++) {binsDCA[i]=0.01*i-5.;}
2694  //for(int i=1; i<41; i++) {binsDCA[i+100]=0.1*i+1.0;}
2695 
2696  Bool_t adddirstatus = TH1::AddDirectoryStatus();
2697  TH1::AddDirectory(kFALSE);
2698  fQA=new TList(); fQA->SetOwner();
2699  fQA->SetName(Form("%s QA",GetName()));
2700  TList* before = new TList(); before->SetOwner();
2701  before->SetName("before");
2702  TList* after = new TList(); after->SetOwner();
2703  after->SetName("after");
2704  fQA->Add(before);
2705  fQA->Add(after);
2706  before->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2707  after->Add(new TH2F("TOFbeta",";p [GeV/c];#beta",kNbinsP,binsP,1000,0.4,1.1)); //0
2708  before->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2709  after->Add(new TH2F("TPCdedx",";p [GeV/c];dEdx",kNbinsP,binsP,500,0,500)); //1
2710  before->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
2711  after->Add(new TH2F("MC pid",";p[GeV/c];species",kNbinsP,binsP,10,-5, 5)); //2
2712  //primary
2713  TH2F* hb = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2714  TH2F* ha = new TH2F("MC primary",";p[GeV/c];primary",kNbinsP,binsP,2,-1,1);
2715  TAxis* axis = NULL;
2716  axis = hb->GetYaxis();
2717  axis->SetBinLabel(1,"secondary");
2718  axis->SetBinLabel(2,"primary");
2719  axis = ha->GetYaxis();
2720  axis->SetBinLabel(1,"secondary");
2721  axis->SetBinLabel(2,"primary");
2722  before->Add(hb); //3
2723  after->Add(ha); //3
2724  //production process
2725  hb = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
2726  -0.5, kMaxMCProcess-0.5);
2727  ha = new TH2F("MC production process",";p[GeV/c];",kNbinsP,binsP,kMaxMCProcess,
2728  -0.5, kMaxMCProcess-0.5);
2729  axis = hb->GetYaxis();
2730  for (Int_t i=0; i<kMaxMCProcess; i++)
2731  {
2732  axis->SetBinLabel(i+1,TMCProcessName[i]);
2733  }
2734  axis = ha->GetYaxis();
2735  for (Int_t i=0; i<kMaxMCProcess; i++)
2736  {
2737  axis->SetBinLabel(i+1,TMCProcessName[i]);
2738  }
2739  before->Add(hb); //4
2740  after->Add(ha); //4
2741  //DCA
2742  before->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2743  after->Add(new TH2F("DCAxy",";p_{t}[GeV/c];DCAxy[cm]", 100, 0., 10., nBinsDCA, binsDCA));//5
2744  before->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
2745  after->Add(new TH2F("DCAz",";p_{t}[GeV/c];DCAz[cm]", 100, 0., 10., nBinsDCA, binsDCA));//6
2746  //first mother
2747  hb = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2748  ha = new TH2F("MC first mother",";p[GeV/c];",kNbinsP,binsP,44,1.,45.);
2749  hb->GetYaxis()->SetBinLabel(1, "#gamma");
2750  ha->GetYaxis()->SetBinLabel(1, "#gamma");
2751  hb->GetYaxis()->SetBinLabel(2, "e^{+}");
2752  ha->GetYaxis()->SetBinLabel(2, "e^{+}");
2753  hb->GetYaxis()->SetBinLabel(3, "e^{-}");
2754  ha->GetYaxis()->SetBinLabel(3, "e^{-}");
2755  hb->GetYaxis()->SetBinLabel(4, "#nu");
2756  ha->GetYaxis()->SetBinLabel(4, "#nu");
2757  hb->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2758  ha->GetYaxis()->SetBinLabel(5, "#mu^{+}");
2759  hb->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2760  ha->GetYaxis()->SetBinLabel(6, "#mu^{-}");
2761  hb->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2762  ha->GetYaxis()->SetBinLabel(7, "#pi^{0}");
2763  hb->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2764  ha->GetYaxis()->SetBinLabel(8, "#pi^{+}");
2765  hb->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2766  ha->GetYaxis()->SetBinLabel(9, "#pi^{-}");
2767  hb->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2768  ha->GetYaxis()->SetBinLabel(10, "K^{0}_{L}");
2769  hb->GetYaxis()->SetBinLabel(11, "K^{+}");
2770  ha->GetYaxis()->SetBinLabel(11, "K^{+}");
2771  hb->GetYaxis()->SetBinLabel(12, "K^{-}");
2772  ha->GetYaxis()->SetBinLabel(12, "K^{-}");
2773  hb->GetYaxis()->SetBinLabel( 13, "n");
2774  ha->GetYaxis()->SetBinLabel( 13, "n");
2775  hb->GetYaxis()->SetBinLabel( 14, "p");
2776  ha->GetYaxis()->SetBinLabel( 14, "p");
2777  hb->GetYaxis()->SetBinLabel(15, "#bar{p}");
2778  ha->GetYaxis()->SetBinLabel(15, "#bar{p}");
2779  hb->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2780  ha->GetYaxis()->SetBinLabel(16, "K^{0}_{S}");
2781  hb->GetYaxis()->SetBinLabel(17, "#eta");
2782  ha->GetYaxis()->SetBinLabel(17, "#eta");
2783  hb->GetYaxis()->SetBinLabel(18, "#Lambda");
2784  ha->GetYaxis()->SetBinLabel(18, "#Lambda");
2785  hb->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2786  ha->GetYaxis()->SetBinLabel(19, "#Sigma^{+}");
2787  hb->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2788  ha->GetYaxis()->SetBinLabel(20, "#Sigma^{0}");
2789  hb->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2790  ha->GetYaxis()->SetBinLabel(21, "#Sigma^{-}");
2791  hb->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2792  ha->GetYaxis()->SetBinLabel(22, "#Theta^{0}");
2793  hb->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2794  ha->GetYaxis()->SetBinLabel(23, "#Theta^{-}");
2795  hb->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2796  ha->GetYaxis()->SetBinLabel(24, "#Omega^{-}");
2797  hb->GetYaxis()->SetBinLabel(25, "#bar{n}");
2798  ha->GetYaxis()->SetBinLabel(25, "#bar{n}");
2799  hb->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2800  ha->GetYaxis()->SetBinLabel(26, "#bar{#Lambda}");
2801  hb->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2802  ha->GetYaxis()->SetBinLabel(27, "#bar{#Sigma}^{-}");
2803  hb->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2804  ha->GetYaxis()->SetBinLabel(28, "#bar{#Sigma}^{0}");
2805  hb->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2806  ha->GetYaxis()->SetBinLabel(29, "#bar{#Sigma}^{+}");
2807  hb->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2808  ha->GetYaxis()->SetBinLabel(30, "#bar{#Theta}^{0}");
2809  hb->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2810  ha->GetYaxis()->SetBinLabel(31, "#bar{#Theta}^{+}");
2811  hb->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2812  ha->GetYaxis()->SetBinLabel(32, "#bar{#Omega}^{+}");
2813  hb->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2814  ha->GetYaxis()->SetBinLabel(33, "#tau^{+}");
2815  hb->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2816  ha->GetYaxis()->SetBinLabel(34, "#tau^{-}");
2817  hb->GetYaxis()->SetBinLabel(35, "D^{+}");
2818  ha->GetYaxis()->SetBinLabel(35, "D^{+}");
2819  hb->GetYaxis()->SetBinLabel(36, "D^{-}");
2820  ha->GetYaxis()->SetBinLabel(36, "D^{-}");
2821  hb->GetYaxis()->SetBinLabel(37, "D^{0}");
2822  ha->GetYaxis()->SetBinLabel(37, "D^{0}");
2823  hb->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2824  ha->GetYaxis()->SetBinLabel(38, "#bar{D}^{0}");
2825  hb->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2826  ha->GetYaxis()->SetBinLabel(39, "D_{s}^{+}");
2827  hb->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2828  ha->GetYaxis()->SetBinLabel(40, "#bar{D_{s}}^{-}");
2829  hb->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2830  ha->GetYaxis()->SetBinLabel(41, "#Lambda_{c}^{+}");
2831  hb->GetYaxis()->SetBinLabel(42, "W^{+}");
2832  ha->GetYaxis()->SetBinLabel(42, "W^{+}");
2833  hb->GetYaxis()->SetBinLabel(43, "W^{-}");
2834  ha->GetYaxis()->SetBinLabel(43, "W^{-}");
2835  hb->GetYaxis()->SetBinLabel(44, "Z^{0}");
2836  ha->GetYaxis()->SetBinLabel(44, "Z^{0}");
2837  before->Add(hb);//7
2838  after->Add(ha);//7
2839 
2840  before->Add(new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2841  after->Add( new TH2F("TOFkElectron",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//8
2842 
2843  before->Add(new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//9
2844  after->Add( new TH2F("TOFkMuon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//9
2845 
2846  before->Add(new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//10
2847  after->Add( new TH2F("TOFkPion",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//10
2848 
2849  before->Add(new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//11
2850  after->Add( new TH2F("TOFkKaon",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//11
2851 
2852  before->Add(new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//12
2853  after->Add( new TH2F("TOFkProton",";p_{t}[GeV/c];TOF signal - IT", kNbinsP,binsP,1000,-2e4, 2e4));//12
2854 
2855  //kink stuff
2856  before->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2857  after->Add(new TH1F("KinkQt",";q_{t}[GeV/c];counts", 200, 0., 0.3));//13
2858 
2859  before->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2860  after->Add(new TH1F("KinkMinv",";m_{inv}(#mu#nu)[GeV/c^{2}];counts;Kink M_{inv}", 200, 0., 0.7));//14
2861 
2862  before->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2863  after->Add(new TH2F("KinkVertex",";x[cm];y[cm];Kink vertex position",250,-250.,250., 250,-250.,250.));//15
2864 
2865  before->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2866  after->Add(new TH2F("KinkAngleMp",";p_{mother}[GeV/c];Kink decay angle [deg];", 100,0.,6., 100,0.,80.));//16
2867 
2868  before->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
2869  after->Add(new TH1F("KinkIndex",";Kink index;counts", 2, 0., 2.));//17
2870 
2871  before->Add(new TH2F("Chi2PerClusterTPC",";p[GeV/c];species",kNbinsP,binsP,50,0.,5.)); //18
2872  after->Add(new TH2F("Chi2PerClusterTPC",";p[GeV/c];species",kNbinsP,binsP,50,0.,5.)); //18
2873 
2874  before->Add(new TH2F("FractionSharedClusterTPC",";p[GeV/c];species",kNbinsP,binsP,50,0.,1.)); //19
2875  after->Add(new TH2F("FractionSharedClusterTPC",";p[GeV/c];species",kNbinsP,binsP,50,0.,1.)); //19
2876 
2877  before->Add(new TH2F("NClustersTPC",";p[GeV/c];species",kNbinsP,binsP,150,50.,200.)); //20
2878  after->Add(new TH2F("NClustersTPC",";p[GeV/c];species",kNbinsP,binsP,150,50.,200.)); //20
2879 
2880  before->Add(new TH2F("Chi2PerClusterITS",";p[GeV/c];species",kNbinsP,binsP,100,0.,50.)); //21
2881  after->Add(new TH2F("Chi2PerClusterITS",";p[GeV/c];species",kNbinsP,binsP,100,0.,50.)); //21
2882 
2883  TH1::AddDirectory(adddirstatus);
2884 }
2885 
2886 //-----------------------------------------------------------------------
2887 Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
2888 {
2889  //get the number of tracks in the input event according source
2890  //selection (ESD tracks, tracklets, MC particles etc.)
2891  AliESDEvent* esd=NULL;
2892  AliAODEvent* aod=NULL; // XZhang 20120615
2893  switch (fParamType)
2894  {
2895  case kSPDtracklet:
2896  if (!fEvent) return 0; // XZhang 20120615
2897  esd = dynamic_cast<AliESDEvent*>(fEvent);
2898  aod = dynamic_cast<AliAODEvent*>(fEvent); // XZhang 20120615
2899 // if (!esd) return 0; // XZhang 20120615
2900 // return esd->GetMultiplicity()->GetNumberOfTracklets(); // XZhang 20120615
2901  if (esd) return esd->GetMultiplicity()->GetNumberOfTracklets(); // XZhang 20120615
2902  if (aod) return aod->GetTracklets()->GetNumberOfTracklets(); // XZhang 20120615
2903  case kMC:
2904  if (!fMCevent) return 0;
2905  return fMCevent->GetNumberOfTracks();
2906  case kPMD:
2907  esd = dynamic_cast<AliESDEvent*>(fEvent);
2908  if (!esd) return 0;
2909  return esd->GetNumberOfPmdTracks();
2910  case kVZERO:
2911  return fgkNumberOfVZEROtracks;
2912  case kBetaVZERO:
2913  return fgkNumberOfVZEROtracks;
2914  case kDeltaVZERO:
2915  return fgkNumberOfVZEROtracks;
2916  case kKappaVZERO:
2917  return fgkNumberOfVZEROtracks;
2918  case kHotfixHI:
2919  return fgkNumberOfVZEROtracks;
2920  case kMUON: // XZhang 20120604
2921  if (!fEvent) return 0; // XZhang 20120604
2922  esd = dynamic_cast<AliESDEvent*>(fEvent); // XZhang 20120604
2923  if (esd) return esd->GetNumberOfMuonTracks(); // XZhang 20120604
2924  return fEvent->GetNumberOfTracks(); // if AOD // XZhang 20120604
2925  case kKink:
2926  esd = dynamic_cast<AliESDEvent*>(fEvent);
2927  if (!esd) return 0;
2928  return esd->GetNumberOfKinks();
2929  case kV0:
2930  esd = dynamic_cast<AliESDEvent*>(fEvent);
2931  if (!esd) return 0;
2932  return esd->GetNumberOfV0s();
2933  default:
2934  if (!fEvent) return 0;
2935  return fEvent->GetNumberOfTracks();
2936  }
2937  return 0;
2938 }
2939 
2940 //-----------------------------------------------------------------------
2941 TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
2942 {
2943  //get the input object according the data source selection:
2944  //(esd tracks, traclets, mc particles,etc...)
2945  AliESDEvent* esd=NULL;
2946  AliAODEvent* aod=NULL; // XZhang 20120615
2947  switch (fParamType)
2948  {
2949  case kSPDtracklet:
2950  if (!fEvent) return NULL; // XZhang 20120615
2951  esd = dynamic_cast<AliESDEvent*>(fEvent);
2952  aod = dynamic_cast<AliAODEvent*>(fEvent); // XZhang 20120615
2953 // if (!esd) return NULL; // XZhang 20120615
2954 // return const_cast<AliMultiplicity*>(esd->GetMultiplicity()); // XZhang 20120615
2955  if (esd) return const_cast<AliMultiplicity*>(esd->GetMultiplicity()); // XZhang 20120615
2956  if (aod) return const_cast<AliAODTracklets*>(aod->GetTracklets()); // XZhang 20120615
2957  case kMC:
2958  if (!fMCevent) return NULL;
2959  return fMCevent->GetTrack(i);
2960  case kPMD:
2961  esd = dynamic_cast<AliESDEvent*>(fEvent);
2962  if (!esd) return NULL;
2963  return esd->GetPmdTrack(i);
2964  case kVZERO:
2965  esd = dynamic_cast<AliESDEvent*>(fEvent);
2966  if (!esd) //contributed by G.Ortona
2967  {
2968  aod = dynamic_cast<AliAODEvent*>(fEvent);
2969  if(!aod)return NULL;
2970  return aod->GetVZEROData();
2971  }
2972  return esd->GetVZEROData();
2973  case kBetaVZERO:
2974  esd = dynamic_cast<AliESDEvent*>(fEvent);
2975  if (!esd) //contributed by G.Ortona
2976  {
2977  aod = dynamic_cast<AliAODEvent*>(fEvent);
2978  if(!aod)return NULL;
2979  return aod->GetVZEROData();
2980  }
2981  return esd->GetVZEROData();
2982  case kDeltaVZERO:
2983  esd = dynamic_cast<AliESDEvent*>(fEvent);
2984  if (!esd) //contributed by G.Ortona
2985  {
2986  aod = dynamic_cast<AliAODEvent*>(fEvent);
2987  if(!aod)return NULL;
2988  return aod->GetVZEROData();
2989  }
2990  return esd->GetVZEROData();
2991  case kKappaVZERO:
2992  esd = dynamic_cast<AliESDEvent*>(fEvent);
2993  if (!esd) //contributed by G.Ortona
2994  {
2995  aod = dynamic_cast<AliAODEvent*>(fEvent);
2996  if(!aod)return NULL;
2997  return aod->GetVZEROData();
2998  }
2999  return esd->GetVZEROData();
3000  case kHotfixHI:
3001  esd = dynamic_cast<AliESDEvent*>(fEvent);
3002  if (!esd) //contributed by G.Ortona
3003  {
3004  aod = dynamic_cast<AliAODEvent*>(fEvent);
3005  if(!aod)return NULL;
3006  return aod->GetVZEROData();
3007  }
3008  return esd->GetVZEROData();
3009  case kMUON: // XZhang 20120604
3010  if (!fEvent) return NULL; // XZhang 20120604
3011  esd = dynamic_cast<AliESDEvent*>(fEvent); // XZhang 20120604
3012  if (esd) return esd->GetMuonTrack(i); // XZhang 20120604
3013  return fEvent->GetTrack(i); // if AOD // XZhang 20120604
3014  case kKink:
3015  esd = dynamic_cast<AliESDEvent*>(fEvent);
3016  if (!esd) return NULL;
3017  return esd->GetKink(i);
3018  case kV0:
3019  esd = dynamic_cast<AliESDEvent*>(fEvent);
3020  if (!esd) return NULL;
3021  return esd->GetV0(i);
3022  default:
3023  if (!fEvent) return NULL;
3024  return fEvent->GetTrack(i);
3025  }
3026 }
3027 
3028 //-----------------------------------------------------------------------
3029 void AliFlowTrackCuts::Clear(Option_t*)
3030 {
3031  //clean up
3032  fMCevent=NULL;
3033  fEvent=NULL;
3034  ClearTrack();
3035 }
3036 
3037 //-----------------------------------------------------------------------
3038 void AliFlowTrackCuts::ClearTrack(Option_t*)
3039 {
3040  //clean up last track
3041  fKink=NULL;
3042  fV0=NULL;
3043  fTrack=NULL;
3044  fMCparticle=NULL;
3045  fTrackLabel=-997;
3046  fTrackWeight=1.0;
3047  fTrackEta=0.0;
3048  fTrackPhi=0.0;
3049  fTrackPt=0.0;
3050  fPOItype=1;
3051  fTrackMass=0.;
3052 }
3053 //-----------------------------------------------------------------------
3054 Bool_t AliFlowTrackCuts::PassesAODpidCut(const AliAODTrack* track )
3055 {
3056  if(!track->GetAODEvent()->GetTOFHeader()){
3057  AliAODPid *pidObj = track->GetDetPid();
3058  if (!pidObj) fESDpid.GetTOFResponse().SetTimeResolution(84.);
3059  else{
3060  Double_t sigmaTOFPidInAOD[10];
3061  pidObj->GetTOFpidResolution(sigmaTOFPidInAOD);
3062  if(sigmaTOFPidInAOD[0] > 84.){
3063  fESDpid.GetTOFResponse().SetTimeResolution(sigmaTOFPidInAOD[0]); // use the electron TOF PID sigma as time resolution (including the T0 used)
3064  }
3065  }
3066  }
3067 
3068  //check if passes the selected pid cut for ESDs
3069  Bool_t pass = kTRUE;
3070  switch (fPIDsource)
3071  {
3072  case kTOFbeta:
3073  if (!PassesTOFbetaCut(track)) pass=kFALSE;
3074  break;
3075  case kTOFbayesian:
3076  if (!PassesTOFbayesianCut(track)) pass=kFALSE;
3077  break;
3078  case kTPCbayesian:
3079  if (!PassesTPCbayesianCut(track)) pass=kFALSE;
3080  break;
3081  case kTPCTOFNsigma:
3082  if (!PassesTPCTOFNsigmaCut(track)) pass = kFALSE;
3083  break;
3084  case kTPCTOFNsigmaPurity:
3085  if(!PassesTPCTOFNsigmaPurityCut(track)) pass = kFALSE;
3086  break;
3087  case kTPCTPCTOFNsigma:
3088  if(!PassesTPCTPCTOFNsigmaCut(track)) pass = kFALSE;
3089  break;
3090  default:
3091  return kTRUE;
3092  break;
3093  }
3094  return pass;
3095 
3096 }
3097 //-----------------------------------------------------------------------
3098 Bool_t AliFlowTrackCuts::PassesESDpidCut(const AliESDtrack* track )
3099 {
3100  //check if passes the selected pid cut for ESDs
3101  Bool_t pass = kTRUE;
3102  switch (fPIDsource)
3103  {
3104  case kTPCpid:
3105  if (!PassesTPCpidCut(track)) pass=kFALSE;
3106  break;
3107  case kTPCdedx:
3108  if (!PassesTPCdedxCut(track)) pass=kFALSE;
3109  break;
3110  case kTOFpid:
3111  if (!PassesTOFpidCut(track)) pass=kFALSE;
3112  break;
3113  case kTOFbeta:
3114  if (!PassesTOFbetaCut(track)) pass=kFALSE;
3115  break;
3116  case kTOFbetaSimple:
3117  if (!PassesTOFbetaSimpleCut(track)) pass=kFALSE;
3118  break;
3119  case kTPCbayesian:
3120  if (!PassesTPCbayesianCut(track)) pass=kFALSE;
3121  break;
3122  // part added by F. Noferini
3123  case kTOFbayesian:
3124  if (!PassesTOFbayesianCut(track)) pass=kFALSE;
3125  break;
3126  // end part added by F. Noferini
3127 
3128  //part added by Natasha
3129  case kTPCNuclei:
3130  if (!PassesNucleiSelection(track)) pass=kFALSE;
3131  break;
3132  //end part added by Natasha
3133 
3134  case kTPCTOFNsigma:
3135  if (!PassesTPCTOFNsigmaCut(track)) pass = kFALSE;
3136  break;
3137  default:
3138  printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
3139  pass=kFALSE;
3140  break;
3141  }
3142  return pass;
3143 }
3144 
3145 //-----------------------------------------------------------------------
3146 Bool_t AliFlowTrackCuts::PassesTOFbetaSimpleCut(const AliESDtrack* track )
3147 {
3148  //check if passes PID cut using timing in TOF
3149  Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
3150  (track->GetTOFsignal() > 12000) &&
3151  (track->GetTOFsignal() < 100000) &&
3152  (track->GetIntegratedLength() > 365);
3153 
3154  if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
3155 
3156  Bool_t statusMatchingHard = TPCTOFagree(track);
3157  if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3158  return kFALSE;
3159 
3160  if (!goodtrack) return kFALSE;
3161 
3162  const Float_t c = 2.99792457999999984e-02;
3163  Float_t p = track->GetP();
3164  Float_t l = track->GetIntegratedLength();
3165  Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
3166  Float_t timeTOF = track->GetTOFsignal()- trackT0;
3167  Float_t beta = l/timeTOF/c;
3168  Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
3169  track->GetIntegratedTimes(integratedTimes);
3170  Float_t betaHypothesis[5] = {0.0,0.0,0.0,0.0,0.0};
3171  Float_t s[5] = {0.0,0.0,0.0,0.0,0.0};
3172  for (Int_t i=0;i<5;i++)
3173  {
3174  betaHypothesis[i] = l/integratedTimes[i]/c;
3175  s[i] = beta-betaHypothesis[i];
3176  }
3177 
3178  switch (fParticleID)
3179  {
3180  case AliPID::kPion:
3181  return ( (s[2]<0.015) && (s[2]>-0.015) &&
3182  (s[3]>0.025) &&
3183  (s[4]>0.03) );
3184  case AliPID::kKaon:
3185  return ( (s[3]<0.015) && (s[3]>-0.015) &&
3186  (s[2]<-0.03) &&
3187  (s[4]>0.03) );
3188  case AliPID::kProton:
3189  return ( (s[4]<0.015) && (s[4]>-0.015) &&
3190  (s[3]<-0.025) &&
3191  (s[2]<-0.025) );
3192  default:
3193  return kFALSE;
3194  }
3195  return kFALSE;
3196 }
3197 
3198 //-----------------------------------------------------------------------
3199 Float_t AliFlowTrackCuts::GetBeta(const AliVTrack* track, Bool_t QAmode)
3200 {
3201  //get beta
3202  Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
3203  track->GetIntegratedTimes(integratedTimes);
3204 
3205  const Float_t c = 2.99792457999999984e-02;
3206  Float_t p = track->P();
3207  Float_t l = integratedTimes[0]*c;
3208  Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
3209  Float_t timeTOF = track->GetTOFsignal()- trackT0;
3210  if(QAmode && timeTOF <= 0) return -999; // avoid division by zero when filling 'before' qa histograms
3211  return l/timeTOF/c;
3212 }
3213 //-----------------------------------------------------------------------
3214 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliAODTrack* track )
3215 {
3216  //check if track passes pid selection with an asymmetric TOF beta cut
3217  if (!fTOFpidCuts)
3218  {
3219  //printf("no TOFpidCuts\n");
3220  return kFALSE;
3221  }
3222 
3223  //check if passes PID cut using timing in TOF
3224  Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
3225  (track->GetTOFsignal() > 12000) &&
3226  (track->GetTOFsignal() < 100000);
3227 
3228  if (!goodtrack) return kFALSE;
3229 
3230  const Float_t c = 2.99792457999999984e-02;
3231  Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
3232  track->GetIntegratedTimes(integratedTimes);
3233  Float_t l = integratedTimes[0]*c;
3234 
3235  goodtrack = goodtrack && (l > 365);
3236 
3237  if (!goodtrack) return kFALSE;
3238 
3239  if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
3240 
3241  Bool_t statusMatchingHard = TPCTOFagree(track);
3242  if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3243  return kFALSE;
3244 
3245 
3246  Float_t beta = GetBeta(track);
3247 
3248  //construct the pid index because it's not AliPID::EParticleType
3249  Int_t pid = 0;
3250  cout<<"TOFbeta: fParticleID = "<<fParticleID<<endl;
3251  switch (fParticleID)
3252  {
3253  case AliPID::kPion:
3254  pid=2;
3255  break;
3256  case AliPID::kKaon:
3257  pid=3;
3258  break;
3259  case AliPID::kProton:
3260  pid=4;
3261  break;
3262  default:
3263  return kFALSE;
3264  }
3265 
3266  //signal to cut on
3267  Float_t p = track->P();
3268  Float_t betahypothesis = l/integratedTimes[pid]/c;
3269  Float_t betadiff = beta-betahypothesis;
3270 
3271  Float_t* arr = fTOFpidCuts->GetMatrixArray();
3272  Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
3273  if (col<0) return kFALSE;
3274  Float_t min = (*fTOFpidCuts)(1,col);
3275  Float_t max = (*fTOFpidCuts)(2,col);
3276 
3277  Bool_t pass = (betadiff>min && betadiff<max);
3278 
3279  return pass;
3280 }
3281 //-----------------------------------------------------------------------
3282 Bool_t AliFlowTrackCuts::PassesTOFbetaCut(const AliESDtrack* track )
3283 {
3284  //check if track passes pid selection with an asymmetric TOF beta cut
3285  if (!fTOFpidCuts)
3286  {
3287  //printf("no TOFpidCuts\n");
3288  return kFALSE;
3289  }
3290 
3291  //check if passes PID cut using timing in TOF
3292  Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFpid) &&
3293  (track->GetTOFsignal() > 12000) &&
3294  (track->GetTOFsignal() < 100000) &&
3295  (track->GetIntegratedLength() > 365);
3296 
3297  if (!fAllowTOFmismatchFlag) {if ((track->GetStatus() & AliESDtrack::kTOFmismatch)) return kFALSE;}
3298 
3299  if (!goodtrack) return kFALSE;
3300 
3301  Bool_t statusMatchingHard = TPCTOFagree(track);
3302  if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3303  return kFALSE;
3304 
3305  Float_t beta = GetBeta(track);
3306  Double_t integratedTimes[9] = {-1.0,-1.0,-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0};
3307  track->GetIntegratedTimes(integratedTimes);
3308 
3309  //construct the pid index because it's not AliPID::EParticleType
3310  Int_t pid = 0;
3311  switch (fParticleID)
3312  {
3313  case AliPID::kPion:
3314  pid=2;
3315  break;
3316  case AliPID::kKaon:
3317  pid=3;
3318  break;
3319  case AliPID::kProton:
3320  pid=4;
3321  break;
3322  default:
3323  return kFALSE;
3324  }
3325 
3326  //signal to cut on
3327  const Float_t c = 2.99792457999999984e-02;
3328  Float_t l = track->GetIntegratedLength();
3329  Float_t p = track->GetP();
3330  Float_t betahypothesis = l/integratedTimes[pid]/c;
3331  Float_t betadiff = beta-betahypothesis;
3332 
3333  Float_t* arr = fTOFpidCuts->GetMatrixArray();
3334  Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(p));
3335  if (col<0) return kFALSE;
3336  Float_t min = (*fTOFpidCuts)(1,col);
3337  Float_t max = (*fTOFpidCuts)(2,col);
3338 
3339  Bool_t pass = (betadiff>min && betadiff<max);
3340 
3341  return pass;
3342 }
3343 
3344 //-----------------------------------------------------------------------
3345 Bool_t AliFlowTrackCuts::PassesTOFpidCut(const AliESDtrack* track) const
3346 {
3347  //check if passes PID cut using default TOF pid
3348  Double_t pidTOF[AliPID::kSPECIES];
3349  track->GetTOFpid(pidTOF);
3350  if (pidTOF[fParticleID]>=fParticleProbability) return kTRUE;
3351  return kFALSE;
3352 }
3353 
3354 //-----------------------------------------------------------------------
3355 Bool_t AliFlowTrackCuts::PassesTPCpidCut(const AliESDtrack* track) const
3356 {
3357  //check if passes PID cut using default TPC pid
3358  Double_t pidTPC[AliPID::kSPECIES];
3359  track->GetTPCpid(pidTPC);
3360  Double_t probablity = 0.;
3361  switch (fParticleID)
3362  {
3363  case AliPID::kPion:
3364  probablity = pidTPC[AliPID::kPion] + pidTPC[AliPID::kMuon];
3365  break;
3366  default:
3367  probablity = pidTPC[fParticleID];
3368  }
3369  if (probablity >= fParticleProbability) return kTRUE;
3370  return kFALSE;
3371 }
3372 
3373 //-----------------------------------------------------------------------
3374 Float_t AliFlowTrackCuts::Getdedx(const AliESDtrack* track) const
3375 {
3376  //get TPC dedx
3377  return track->GetTPCsignal();
3378 }
3379 
3380 //-----------------------------------------------------------------------
3381 Bool_t AliFlowTrackCuts::PassesTPCdedxCut(const AliESDtrack* track)
3382 {
3383  //check if passes PID cut using dedx signal in the TPC
3384  if (!fTPCpidCuts)
3385  {
3386  //printf("no TPCpidCuts\n");
3387  return kFALSE;
3388  }
3389 
3390  const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
3391  if (!tpcparam) return kFALSE;
3392  Double_t p = tpcparam->GetP();
3393  Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(p, fParticleID);
3394  Float_t sigTPC = track->GetTPCsignal();
3395  Float_t s = (sigTPC-sigExp)/sigExp;
3396 
3397  Float_t* arr = fTPCpidCuts->GetMatrixArray();
3398  Int_t arrSize = fTPCpidCuts->GetNcols();
3399  Int_t col = TMath::BinarySearch( arrSize, arr, static_cast<Float_t>(p));
3400  if (col<0) return kFALSE;
3401  Float_t min = (*fTPCpidCuts)(1,col);
3402  Float_t max = (*fTPCpidCuts)(2,col);
3403 
3404  //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
3405  return (s>min && s<max);
3406 }
3407 
3408 //-----------------------------------------------------------------------
3409 void AliFlowTrackCuts::InitPIDcuts()
3410 {
3411  //init matrices with PID cuts
3412  TMatrixF* t = NULL;
3413  if (!fTPCpidCuts)
3414  {
3415  if (fParticleID==AliPID::kPion)
3416  {
3417  t = new TMatrixF(3,15);
3418  (*t)(0,0) = 0.20; (*t)(1,0) = -0.4; (*t)(2,0) = 0.0;
3419  (*t)(0,1) = 0.25; (*t)(1,1) = -0.4; (*t)(2,1) = 0.1;
3420  (*t)(0,2) = 0.30; (*t)(1,2) = -0.4; (*t)(2,2) = 0.2;
3421  (*t)(0,3) = 0.35; (*t)(1,3) = -0.4; (*t)(2,3) = 0.2;
3422  (*t)(0,4) = 0.40; (*t)(1,4) = -0.4; (*t)(2,4) = 0.3;
3423  (*t)(0,5) = 0.45; (*t)(1,5) = -0.4; (*t)(2,5) = 0.3;
3424  (*t)(0,6) = 0.50; (*t)(1,6) = -0.4; (*t)(2,6) = 0.25;
3425  (*t)(0,7) = 0.55; (*t)(1,7) = -0.4; (*t)(2,7) = 0.15;
3426  (*t)(0,8) = 0.60; (*t)(1,8) = -0.4; (*t)(2,8) = 0.1;
3427  (*t)(0,9) = 0.65; (*t)(1,9) = -0.4; (*t)(2,9) = 0.05;
3428  (*t)(0,10) = 0.70; (*t)(1,10) = -0.4; (*t)(2,10) = 0;
3429  (*t)(0,11) = 0.75; (*t)(1,11) = -0.4; (*t)(2,11) = 0;
3430  (*t)(0,12) = 0.80; (*t)(1,12) = -0.4; (*t)(2,12) = -0.05;
3431  (*t)(0,13) = 0.85; (*t)(1,13) = -0.4; (*t)(2,13) = -0.1;
3432  (*t)(0,14) = 0.90; (*t)(1,14) = 0; (*t)(2,14) = 0;
3433  }
3434  else
3435  if (fParticleID==AliPID::kKaon)
3436  {
3437  t = new TMatrixF(3,12);
3438  (*t)(0,0) = 0.20; (*t)(1,0) = -0.2; (*t)(2,0) = 0.2;
3439  (*t)(0,1) = 0.25; (*t)(1,1) = -0.2; (*t)(2,1) = 0.2;
3440  (*t)(0,2) = 0.30; (*t)(1,2) = -0.2; (*t)(2,2) = 0.2;
3441  (*t)(0,3) = 0.35; (*t)(1,3) = -0.2; (*t)(2,3) = 0.2;
3442  (*t)(0,4) = 0.40; (*t)(1,4) = -0.1; (*t)(2,4) = 0.2;
3443  (*t)(0,5) = 0.45; (*t)(1,5) = -0.1; (*t)(2,5) = 0.2;
3444  (*t)(0,6) = 0.50; (*t)(1,6) =-0.05; (*t)(2,6) = 0.2;
3445  (*t)(0,7) = 0.55; (*t)(1,7) = -0.1; (*t)(2,7) = 0.1;
3446  (*t)(0,8) = 0.60; (*t)(1,8) =-0.05; (*t)(2,8) = 0.1;
3447  (*t)(0,9) = 0.65; (*t)(1,9) = 0; (*t)(2,9) = 0.15;
3448  (*t)(0,10) = 0.70; (*t)(1,10) = 0.05; (*t)(2,10) = 0.2;
3449  (*t)(0,11) = 0.75; (*t)(1,11) = 0; (*t)(2,11) = 0;
3450  }
3451  else
3452  if (fParticleID==AliPID::kProton)
3453  {
3454  t = new TMatrixF(3,9);
3455  (*t)(0,0) = 0.20; (*t)(1,0) = -0.1; (*t)(2,0) = 0.1;
3456  (*t)(0,1) = 0.25; (*t)(1,1) = -0.2; (*t)(2,1) = 0.2;
3457  (*t)(0,2) = 0.80; (*t)(1,2) = -0.1; (*t)(2,2) = 0.2;
3458  (*t)(0,3) = 0.85; (*t)(1,3) =-0.05; (*t)(2,3) = 0.2;
3459  (*t)(0,4) = 0.90; (*t)(1,4) =-0.05; (*t)(2,4) = 0.25;
3460  (*t)(0,5) = 0.95; (*t)(1,5) =-0.05; (*t)(2,5) = 0.25;
3461  (*t)(0,6) = 1.00; (*t)(1,6) = -0.1; (*t)(2,6) = 0.25;
3462  (*t)(0,7) = 1.10; (*t)(1,7) =-0.05; (*t)(2,7) = 0.3;
3463  (*t)(0,8) = 1.20; (*t)(1,8) = 0; (*t)(2,8) = 0;
3464  }
3465  delete fTPCpidCuts;
3466  fTPCpidCuts=t;
3467  }
3468  t = NULL;
3469  if (!fTOFpidCuts)
3470  {
3471  if (fParticleID==AliPID::kPion)
3472  {
3473  //TOF pions, 0.9 purity
3474  t = new TMatrixF(3,61);
3475  (*t)(0,0) = 0.000; (*t)(1,0) = 0.000; (*t)(2,0) = 0.000;
3476  (*t)(0,1) = 0.050; (*t)(1,1) = 0.000; (*t)(2,1) = 0.000;
3477  (*t)(0,2) = 0.100; (*t)(1,2) = 0.000; (*t)(2,2) = 0.000;
3478  (*t)(0,3) = 0.150; (*t)(1,3) = 0.000; (*t)(2,3) = 0.000;
3479  (*t)(0,4) = 0.200; (*t)(1,4) = -0.030; (*t)(2,4) = 0.030;
3480  (*t)(0,5) = 0.250; (*t)(1,5) = -0.036; (*t)(2,5) = 0.032;
3481  (*t)(0,6) = 0.300; (*t)(1,6) = -0.038; (*t)(2,6) = 0.032;
3482  (*t)(0,7) = 0.350; (*t)(1,7) = -0.034; (*t)(2,7) = 0.032;
3483  (*t)(0,8) = 0.400; (*t)(1,8) = -0.032; (*t)(2,8) = 0.020;
3484  (*t)(0,9) = 0.450; (*t)(1,9) = -0.030; (*t)(2,9) = 0.020;
3485  (*t)(0,10) = 0.500; (*t)(1,10) = -0.030; (*t)(2,10) = 0.020;
3486  (*t)(0,11) = 0.550; (*t)(1,11) = -0.030; (*t)(2,11) = 0.020;
3487  (*t)(0,12) = 0.600; (*t)(1,12) = -0.030; (*t)(2,12) = 0.020;
3488  (*t)(0,13) = 0.650; (*t)(1,13) = -0.030; (*t)(2,13) = 0.020;
3489  (*t)(0,14) = 0.700; (*t)(1,14) = -0.030; (*t)(2,14) = 0.020;
3490  (*t)(0,15) = 0.750; (*t)(1,15) = -0.030; (*t)(2,15) = 0.020;
3491  (*t)(0,16) = 0.800; (*t)(1,16) = -0.030; (*t)(2,16) = 0.020;
3492  (*t)(0,17) = 0.850; (*t)(1,17) = -0.030; (*t)(2,17) = 0.020;
3493  (*t)(0,18) = 0.900; (*t)(1,18) = -0.030; (*t)(2,18) = 0.020;
3494  (*t)(0,19) = 0.950; (*t)(1,19) = -0.028; (*t)(2,19) = 0.028;
3495  (*t)(0,20) = 1.000; (*t)(1,20) = -0.028; (*t)(2,20) = 0.028;
3496  (*t)(0,21) = 1.100; (*t)(1,21) = -0.028; (*t)(2,21) = 0.028;
3497  (*t)(0,22) = 1.200; (*t)(1,22) = -0.026; (*t)(2,22) = 0.028;
3498  (*t)(0,23) = 1.300; (*t)(1,23) = -0.024; (*t)(2,23) = 0.028;
3499  (*t)(0,24) = 1.400; (*t)(1,24) = -0.020; (*t)(2,24) = 0.028;
3500  (*t)(0,25) = 1.500; (*t)(1,25) = -0.018; (*t)(2,25) = 0.028;
3501  (*t)(0,26) = 1.600; (*t)(1,26) = -0.016; (*t)(2,26) = 0.028;
3502  (*t)(0,27) = 1.700; (*t)(1,27) = -0.014; (*t)(2,27) = 0.028;
3503  (*t)(0,28) = 1.800; (*t)(1,28) = -0.012; (*t)(2,28) = 0.026;
3504  (*t)(0,29) = 1.900; (*t)(1,29) = -0.010; (*t)(2,29) = 0.026;
3505  (*t)(0,30) = 2.000; (*t)(1,30) = -0.008; (*t)(2,30) = 0.026;
3506  (*t)(0,31) = 2.100; (*t)(1,31) = -0.008; (*t)(2,31) = 0.024;
3507  (*t)(0,32) = 2.200; (*t)(1,32) = -0.006; (*t)(2,32) = 0.024;
3508  (*t)(0,33) = 2.300; (*t)(1,33) = -0.004; (*t)(2,33) = 0.024;
3509  (*t)(0,34) = 2.400; (*t)(1,34) = -0.004; (*t)(2,34) = 0.024;
3510  (*t)(0,35) = 2.500; (*t)(1,35) = -0.002; (*t)(2,35) = 0.024;
3511  (*t)(0,36) = 2.600; (*t)(1,36) = -0.002; (*t)(2,36) = 0.024;
3512  (*t)(0,37) = 2.700; (*t)(1,37) = 0.000; (*t)(2,37) = 0.024;
3513  (*t)(0,38) = 2.800; (*t)(1,38) = 0.000; (*t)(2,38) = 0.026;
3514  (*t)(0,39) = 2.900; (*t)(1,39) = 0.000; (*t)(2,39) = 0.024;
3515  (*t)(0,40) = 3.000; (*t)(1,40) = 0.002; (*t)(2,40) = 0.026;
3516  (*t)(0,41) = 3.100; (*t)(1,41) = 0.002; (*t)(2,41) = 0.026;
3517  (*t)(0,42) = 3.200; (*t)(1,42) = 0.002; (*t)(2,42) = 0.026;
3518  (*t)(0,43) = 3.300; (*t)(1,43) = 0.002; (*t)(2,43) = 0.026;
3519  (*t)(0,44) = 3.400; (*t)(1,44) = 0.002; (*t)(2,44) = 0.026;
3520  (*t)(0,45) = 3.500; (*t)(1,45) = 0.002; (*t)(2,45) = 0.026;
3521  (*t)(0,46) = 3.600; (*t)(1,46) = 0.002; (*t)(2,46) = 0.026;
3522  (*t)(0,47) = 3.700; (*t)(1,47) = 0.002; (*t)(2,47) = 0.026;
3523  (*t)(0,48) = 3.800; (*t)(1,48) = 0.002; (*t)(2,48) = 0.026;
3524  (*t)(0,49) = 3.900; (*t)(1,49) = 0.004; (*t)(2,49) = 0.024;
3525  (*t)(0,50) = 4.000; (*t)(1,50) = 0.004; (*t)(2,50) = 0.026;
3526  (*t)(0,51) = 4.100; (*t)(1,51) = 0.004; (*t)(2,51) = 0.026;
3527  (*t)(0,52) = 4.200; (*t)(1,52) = 0.004; (*t)(2,52) = 0.024;
3528  (*t)(0,53) = 4.300; (*t)(1,53) = 0.006; (*t)(2,53) = 0.024;
3529  (*t)(0,54) = 4.400; (*t)(1,54) = 0.000; (*t)(2,54) = 0.000;
3530  (*t)(0,55) = 4.500; (*t)(1,55) = 0.000; (*t)(2,55) = 0.000;
3531  (*t)(0,56) = 4.600; (*t)(1,56) = 0.000; (*t)(2,56) = 0.000;
3532  (*t)(0,57) = 4.700; (*t)(1,57) = 0.000; (*t)(2,57) = 0.000;
3533  (*t)(0,58) = 4.800; (*t)(1,58) = 0.000; (*t)(2,58) = 0.000;
3534  (*t)(0,59) = 4.900; (*t)(1,59) = 0.000; (*t)(2,59) = 0.000;
3535  (*t)(0,60) = 5.900; (*t)(1,60) = 0.000; (*t)(2,60) = 0.000;
3536  }
3537  else
3538  if (fParticleID==AliPID::kProton)
3539  {
3540  //TOF protons, 0.9 purity
3541  t = new TMatrixF(3,61);
3542  (*t)(0,0) = 0.000; (*t)(1,0) = 0.000; (*t)(2,0) = 0.000;
3543  (*t)(0,1) = 0.050; (*t)(1,1) = 0.000; (*t)(2,1) = 0.000;
3544  (*t)(0,2) = 0.100; (*t)(1,2) = 0.000; (*t)(2,2) = 0.000;
3545  (*t)(0,3) = 0.150; (*t)(1,3) = 0.000; (*t)(2,3) = 0.000;
3546  (*t)(0,4) = 0.200; (*t)(1,4) = -0.07; (*t)(2,4) = 0.07;
3547  (*t)(0,5) = 0.200; (*t)(1,5) = -0.07; (*t)(2,5) = 0.07;
3548  (*t)(0,6) = 0.200; (*t)(1,6) = -0.07; (*t)(2,6) = 0.07;
3549  (*t)(0,7) = 0.200; (*t)(1,7) = -0.07; (*t)(2,7) = 0.07;
3550  (*t)(0,8) = 0.200; (*t)(1,8) = -0.07; (*t)(2,8) = 0.07;
3551  (*t)(0,9) = 0.200; (*t)(1,9) = -0.07; (*t)(2,9) = 0.07;
3552  (*t)(0,10) = 0.200; (*t)(1,10) = -0.07; (*t)(2,10) = 0.07;
3553  (*t)(0,11) = 0.200; (*t)(1,11) = -0.07; (*t)(2,11) = 0.07;
3554  (*t)(0,12) = 0.200; (*t)(1,12) = -0.07; (*t)(2,12) = 0.07;
3555  (*t)(0,13) = 0.200; (*t)(1,13) = -0.07; (*t)(2,13) = 0.07;
3556  (*t)(0,14) = 0.200; (*t)(1,14) = -0.07; (*t)(2,14) = 0.07;
3557  (*t)(0,15) = 0.200; (*t)(1,15) = -0.07; (*t)(2,15) = 0.07;
3558  (*t)(0,16) = 0.200; (*t)(1,16) = -0.07; (*t)(2,16) = 0.07;
3559  (*t)(0,17) = 0.850; (*t)(1,17) = -0.070; (*t)(2,17) = 0.070;
3560  (*t)(0,18) = 0.900; (*t)(1,18) = -0.072; (*t)(2,18) = 0.072;
3561  (*t)(0,19) = 0.950; (*t)(1,19) = -0.072; (*t)(2,19) = 0.072;
3562  (*t)(0,20) = 1.000; (*t)(1,20) = -0.074; (*t)(2,20) = 0.074;
3563  (*t)(0,21) = 1.100; (*t)(1,21) = -0.032; (*t)(2,21) = 0.032;
3564  (*t)(0,22) = 1.200; (*t)(1,22) = -0.026; (*t)(2,22) = 0.026;
3565  (*t)(0,23) = 1.300; (*t)(1,23) = -0.026; (*t)(2,23) = 0.026;
3566  (*t)(0,24) = 1.400; (*t)(1,24) = -0.024; (*t)(2,24) = 0.024;
3567  (*t)(0,25) = 1.500; (*t)(1,25) = -0.024; (*t)(2,25) = 0.024;
3568  (*t)(0,26) = 1.600; (*t)(1,26) = -0.026; (*t)(2,26) = 0.026;
3569  (*t)(0,27) = 1.700; (*t)(1,27) = -0.026; (*t)(2,27) = 0.026;
3570  (*t)(0,28) = 1.800; (*t)(1,28) = -0.026; (*t)(2,28) = 0.026;
3571  (*t)(0,29) = 1.900; (*t)(1,29) = -0.026; (*t)(2,29) = 0.026;
3572  (*t)(0,30) = 2.000; (*t)(1,30) = -0.026; (*t)(2,30) = 0.026;
3573  (*t)(0,31) = 2.100; (*t)(1,31) = -0.026; (*t)(2,31) = 0.026;
3574  (*t)(0,32) = 2.200; (*t)(1,32) = -0.026; (*t)(2,32) = 0.024;
3575  (*t)(0,33) = 2.300; (*t)(1,33) = -0.028; (*t)(2,33) = 0.022;
3576  (*t)(0,34) = 2.400; (*t)(1,34) = -0.028; (*t)(2,34) = 0.020;
3577  (*t)(0,35) = 2.500; (*t)(1,35) = -0.028; (*t)(2,35) = 0.018;
3578  (*t)(0,36) = 2.600; (*t)(1,36) = -0.028; (*t)(2,36) = 0.016;
3579  (*t)(0,37) = 2.700; (*t)(1,37) = -0.028; (*t)(2,37) = 0.016;
3580  (*t)(0,38) = 2.800; (*t)(1,38) = -0.030; (*t)(2,38) = 0.014;
3581  (*t)(0,39) = 2.900; (*t)(1,39) = -0.030; (*t)(2,39) = 0.012;
3582  (*t)(0,40) = 3.000; (*t)(1,40) = -0.030; (*t)(2,40) = 0.012;
3583  (*t)(0,41) = 3.100; (*t)(1,41) = -0.030; (*t)(2,41) = 0.010;
3584  (*t)(0,42) = 3.200; (*t)(1,42) = -0.030; (*t)(2,42) = 0.010;
3585  (*t)(0,43) = 3.300; (*t)(1,43) = -0.030; (*t)(2,43) = 0.010;
3586  (*t)(0,44) = 3.400; (*t)(1,44) = -0.030; (*t)(2,44) = 0.008;
3587  (*t)(0,45) = 3.500; (*t)(1,45) = -0.030; (*t)(2,45) = 0.008;
3588  (*t)(0,46) = 3.600; (*t)(1,46) = -0.030; (*t)(2,46) = 0.008;
3589  (*t)(0,47) = 3.700; (*t)(1,47) = -0.030; (*t)(2,47) = 0.006;
3590  (*t)(0,48) = 3.800; (*t)(1,48) = -0.030; (*t)(2,48) = 0.006;
3591  (*t)(0,49) = 3.900; (*t)(1,49) = -0.030; (*t)(2,49) = 0.006;
3592  (*t)(0,50) = 4.000; (*t)(1,50) = -0.028; (*t)(2,50) = 0.004;
3593  (*t)(0,51) = 4.100; (*t)(1,51) = -0.030; (*t)(2,51) = 0.004;
3594  (*t)(0,52) = 4.200; (*t)(1,52) = -0.030; (*t)(2,52) = 0.004;
3595  (*t)(0,53) = 4.300; (*t)(1,53) = -0.028; (*t)(2,53) = 0.002;
3596  (*t)(0,54) = 4.400; (*t)(1,54) = -0.030; (*t)(2,54) = 0.002;
3597  (*t)(0,55) = 4.500; (*t)(1,55) = -0.028; (*t)(2,55) = 0.002;
3598  (*t)(0,56) = 4.600; (*t)(1,56) = -0.028; (*t)(2,56) = 0.002;
3599  (*t)(0,57) = 4.700; (*t)(1,57) = -0.028; (*t)(2,57) = 0.000;
3600  (*t)(0,58) = 4.800; (*t)(1,58) = -0.028; (*t)(2,58) = 0.002;
3601  (*t)(0,59) = 4.900; (*t)(1,59) = 0.000; (*t)(2,59) = 0.000;
3602  (*t)(0,60) = 5.900; (*t)(1,60) = 0.000; (*t)(2,60) = 0.000;
3603  }
3604  else
3605  if (fParticleID==AliPID::kKaon)
3606  {
3607  //TOF kaons, 0.9 purity
3608  t = new TMatrixF(3,61);
3609  (*t)(0,0) = 0.000; (*t)(1,0) = 0.000; (*t)(2,0) = 0.000;
3610  (*t)(0,1) = 0.050; (*t)(1,1) = 0.000; (*t)(2,1) = 0.000;
3611  (*t)(0,2) = 0.100; (*t)(1,2) = 0.000; (*t)(2,2) = 0.000;
3612  (*t)(0,3) = 0.150; (*t)(1,3) = 0.000; (*t)(2,3) = 0.000;
3613  (*t)(0,4) = 0.200; (*t)(1,4) = -0.05; (*t)(2,4) = 0.05;
3614  (*t)(0,5) = 0.200; (*t)(1,5) = -0.05; (*t)(2,5) = 0.05;
3615  (*t)(0,6) = 0.200; (*t)(1,6) = -0.05; (*t)(2,6) = 0.05;
3616  (*t)(0,7) = 0.200; (*t)(1,7) = -0.05; (*t)(2,7) = 0.05;
3617  (*t)(0,8) = 0.200; (*t)(1,8) = -0.05; (*t)(2,8) = 0.05;
3618  (*t)(0,9) = 0.200; (*t)(1,9) = -0.05; (*t)(2,9) = 0.05;
3619  (*t)(0,10) = 0.200; (*t)(1,10) = -0.05; (*t)(2,10) = 0.05;
3620  (*t)(0,11) = 0.550; (*t)(1,11) = -0.026; (*t)(2,11) = 0.026;
3621  (*t)(0,12) = 0.600; (*t)(1,12) = -0.026; (*t)(2,12) = 0.026;
3622  (*t)(0,13) = 0.650; (*t)(1,13) = -0.026; (*t)(2,13) = 0.026;
3623  (*t)(0,14) = 0.700; (*t)(1,14) = -0.026; (*t)(2,14) = 0.026;
3624  (*t)(0,15) = 0.750; (*t)(1,15) = -0.026; (*t)(2,15) = 0.026;
3625  (*t)(0,16) = 0.800; (*t)(1,16) = -0.026; (*t)(2,16) = 0.026;
3626  (*t)(0,17) = 0.850; (*t)(1,17) = -0.024; (*t)(2,17) = 0.024;
3627  (*t)(0,18) = 0.900; (*t)(1,18) = -0.024; (*t)(2,18) = 0.024;
3628  (*t)(0,19) = 0.950; (*t)(1,19) = -0.024; (*t)(2,19) = 0.024;
3629  (*t)(0,20) = 1.000; (*t)(1,20) = -0.024; (*t)(2,20) = 0.024;
3630  (*t)(0,21) = 1.100; (*t)(1,21) = -0.024; (*t)(2,21) = 0.024;
3631  (*t)(0,22) = 1.200; (*t)(1,22) = -0.024; (*t)(2,22) = 0.022;
3632  (*t)(0,23) = 1.300; (*t)(1,23) = -0.024; (*t)(2,23) = 0.020;
3633  (*t)(0,24) = 1.400; (*t)(1,24) = -0.026; (*t)(2,24) = 0.016;
3634  (*t)(0,25) = 1.500; (*t)(1,25) = -0.028; (*t)(2,25) = 0.014;
3635  (*t)(0,26) = 1.600; (*t)(1,26) = -0.028; (*t)(2,26) = 0.012;
3636  (*t)(0,27) = 1.700; (*t)(1,27) = -0.028; (*t)(2,27) = 0.010;
3637  (*t)(0,28) = 1.800; (*t)(1,28) = -0.028; (*t)(2,28) = 0.010;
3638  (*t)(0,29) = 1.900; (*t)(1,29) = -0.028; (*t)(2,29) = 0.008;
3639  (*t)(0,30) = 2.000; (*t)(1,30) = -0.028; (*t)(2,30) = 0.006;
3640  (*t)(0,31) = 2.100; (*t)(1,31) = -0.026; (*t)(2,31) = 0.006;
3641  (*t)(0,32) = 2.200; (*t)(1,32) = -0.024; (*t)(2,32) = 0.004;
3642  (*t)(0,33) = 2.300; (*t)(1,33) = -0.020; (*t)(2,33) = 0.002;
3643  (*t)(0,34) = 2.400; (*t)(1,34) = -0.020; (*t)(2,34) = 0.002;
3644  (*t)(0,35) = 2.500; (*t)(1,35) = -0.018; (*t)(2,35) = 0.000;
3645  (*t)(0,36) = 2.600; (*t)(1,36) = -0.016; (*t)(2,36) = 0.000;
3646  (*t)(0,37) = 2.700; (*t)(1,37) = -0.014; (*t)(2,37) = -0.002;
3647  (*t)(0,38) = 2.800; (*t)(1,38) = -0.014; (*t)(2,38) = -0.004;
3648  (*t)(0,39) = 2.900; (*t)(1,39) = -0.012; (*t)(2,39) = -0.004;
3649  (*t)(0,40) = 3.000; (*t)(1,40) = -0.010; (*t)(2,40) = -0.006;
3650  (*t)(0,41) = 3.100; (*t)(1,41) = 0.000; (*t)(2,41) = 0.000;
3651  (*t)(0,42) = 3.200; (*t)(1,42) = 0.000; (*t)(2,42) = 0.000;
3652  (*t)(0,43) = 3.300; (*t)(1,43) = 0.000; (*t)(2,43) = 0.000;
3653  (*t)(0,44) = 3.400; (*t)(1,44) = 0.000; (*t)(2,44) = 0.000;
3654  (*t)(0,45) = 3.500; (*t)(1,45) = 0.000; (*t)(2,45) = 0.000;
3655  (*t)(0,46) = 3.600; (*t)(1,46) = 0.000; (*t)(2,46) = 0.000;
3656  (*t)(0,47) = 3.700; (*t)(1,47) = 0.000; (*t)(2,47) = 0.000;
3657  (*t)(0,48) = 3.800; (*t)(1,48) = 0.000; (*t)(2,48) = 0.000;
3658  (*t)(0,49) = 3.900; (*t)(1,49) = 0.000; (*t)(2,49) = 0.000;
3659  (*t)(0,50) = 4.000; (*t)(1,50) = 0.000; (*t)(2,50) = 0.000;
3660  (*t)(0,51) = 4.100; (*t)(1,51) = 0.000; (*t)(2,51) = 0.000;
3661  (*t)(0,52) = 4.200; (*t)(1,52) = 0.000; (*t)(2,52) = 0.000;
3662  (*t)(0,53) = 4.300; (*t)(1,53) = 0.000; (*t)(2,53) = 0.000;
3663  (*t)(0,54) = 4.400; (*t)(1,54) = 0.000; (*t)(2,54) = 0.000;
3664  (*t)(0,55) = 4.500; (*t)(1,55) = 0.000; (*t)(2,55) = 0.000;
3665  (*t)(0,56) = 4.600; (*t)(1,56) = 0.000; (*t)(2,56) = 0.000;
3666  (*t)(0,57) = 4.700; (*t)(1,57) = 0.000; (*t)(2,57) = 0.000;
3667  (*t)(0,58) = 4.800; (*t)(1,58) = 0.000; (*t)(2,58) = 0.000;
3668  (*t)(0,59) = 4.900; (*t)(1,59) = 0.000; (*t)(2,59) = 0.000;
3669  (*t)(0,60) = 5.900; (*t)(1,60) = 0.000; (*t)(2,60) = 0.000;
3670  }
3671  delete fTOFpidCuts;
3672  fTOFpidCuts=t;
3673  }
3674 }
3675 
3676 //-----------------------------------------------------------------------
3677 //-----------------------------------------------------------------------
3678 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliAODTrack* track)
3679 {
3680  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3681  Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3682 
3683  Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
3684 
3685  if(! kTPC) return kFALSE;
3686 
3687  fProbBayes = 0.0;
3688 
3689 switch (fParticleID)
3690  {
3691  case AliPID::kPion:
3692  fProbBayes = probabilities[2];
3693  break;
3694  case AliPID::kKaon:
3695  fProbBayes = probabilities[3];
3696  break;
3697  case AliPID::kProton:
3698  fProbBayes = probabilities[4];
3699  break;
3700  case AliPID::kElectron:
3701  fProbBayes = probabilities[0];
3702  break;
3703  case AliPID::kMuon:
3704  fProbBayes = probabilities[1];
3705  break;
3706  case AliPID::kDeuteron:
3707  fProbBayes = probabilities[5];
3708  break;
3709  case AliPID::kTriton:
3710  fProbBayes = probabilities[6];
3711  break;
3712  case AliPID::kHe3:
3713  fProbBayes = probabilities[7];
3714  break;
3715  default:
3716  return kFALSE;
3717  }
3718 
3719  if(fProbBayes > fParticleProbability){
3720  if(!fCutCharge)
3721  return kTRUE;
3722  else if (fCutCharge && fCharge * track->Charge() > 0)
3723  return kTRUE;
3724  }
3725  return kFALSE;
3726 }
3727 //-----------------------------------------------------------------------
3728 Bool_t AliFlowTrackCuts::PassesTPCbayesianCut(const AliESDtrack* track)
3729 {
3730  //cut on TPC bayesian pid
3731  //if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3732 
3733  //Bool_t statusMatchingHard = TPCTOFagree(track);
3734  //if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3735  // return kFALSE;
3736  fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
3737  Int_t kTPC = fBayesianResponse->GetCurrentMask(0); // is TPC on
3738  //Int_t kTOF = fBayesianResponse->GetCurrentMask(1); // is TOF on
3739 
3740  if(! kTPC) return kFALSE;
3741 
3742  // Bool_t statusMatchingHard = 1;
3743  // Float_t mismProb = 0;
3744  // if(kTOF){
3745  // statusMatchingHard = TPCTOFagree(track);
3746  // mismProb = fBayesianResponse->GetTOFMismProb();
3747  // }
3748  // if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3749  // return kFALSE;
3750 
3751  Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3752 
3753  fProbBayes = 0.0;
3754 
3755  switch (fParticleID)
3756  {
3757  case AliPID::kPion:
3758  fProbBayes = probabilities[2];
3759  break;
3760  case AliPID::kKaon:
3761  fProbBayes = probabilities[3];
3762  break;
3763  case AliPID::kProton:
3764  fProbBayes = probabilities[4];
3765  break;
3766  case AliPID::kElectron:
3767  fProbBayes = probabilities[0];
3768  break;
3769  case AliPID::kMuon:
3770  fProbBayes = probabilities[1];
3771  break;
3772  case AliPID::kDeuteron:
3773  fProbBayes = probabilities[5];
3774  break;
3775  case AliPID::kTriton:
3776  fProbBayes = probabilities[6];
3777  break;
3778  case AliPID::kHe3:
3779  fProbBayes = probabilities[7];
3780  break;
3781  default:
3782  return kFALSE;
3783  }
3784 
3785  if(fProbBayes > fParticleProbability)
3786  {
3787  if(!fCutCharge)
3788  return kTRUE;
3789  else if (fCutCharge && fCharge * track->GetSign() > 0)
3790  return kTRUE;
3791  }
3792  return kFALSE;
3793 }
3794 //-----------------------------------------------------------------------
3795 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliAODTrack* track)
3796 {
3797  //check is track passes bayesian combined TOF+TPC pid cut
3798  Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3799  (track->GetStatus() & AliESDtrack::kTIME) &&
3800  (track->GetTOFsignal() > 12000) &&
3801  (track->GetTOFsignal() < 100000);
3802 
3803  if (! goodtrack)
3804  return kFALSE;
3805 
3806  Bool_t statusMatchingHard = TPCTOFagree(track);
3807  if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3808  return kFALSE;
3809 
3810  fBayesianResponse->ComputeProb(track,track->GetAODEvent()); // fCurrCentr is needed for mismatch fraction
3811  Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mism$
3812 
3813  Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3814 
3815  fProbBayes = 0.0;
3816 
3817  switch (fParticleID)
3818  {
3819  case AliPID::kPion:
3820  fProbBayes = probabilities[2];
3821  break;
3822  case AliPID::kKaon:
3823  fProbBayes = probabilities[3];
3824  break;
3825  case AliPID::kProton:
3826  fProbBayes = probabilities[4];
3827  break;
3828  case AliPID::kElectron:
3829  fProbBayes = probabilities[0];
3830  break;
3831  case AliPID::kMuon:
3832  fProbBayes = probabilities[1];
3833  break;
3834  case AliPID::kDeuteron:
3835  fProbBayes = probabilities[5];
3836  break;
3837  case AliPID::kTriton:
3838  fProbBayes = probabilities[6];
3839  break;
3840  case AliPID::kHe3:
3841  fProbBayes = probabilities[7];
3842  break;
3843  default:
3844  return kFALSE;
3845  }
3846 
3847  if(fProbBayes > fParticleProbability && mismProb < 0.5){
3848  if(!fCutCharge)
3849  return kTRUE;
3850  else if (fCutCharge && fCharge * track->Charge() > 0)
3851  return kTRUE;
3852  }
3853  return kFALSE;
3854 
3855 }
3856 //-----------------------------------------------------------------------
3857 // part added by F. Noferini (some methods)
3858 Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(const AliESDtrack* track)
3859 {
3860  //check is track passes bayesian combined TOF+TPC pid cut
3861  Bool_t goodtrack = (track->GetStatus() & AliESDtrack::kTOFout) &&
3862  (track->GetStatus() & AliESDtrack::kTIME) &&
3863  (track->GetTOFsignal() > 12000) &&
3864  (track->GetTOFsignal() < 100000) &&
3865  (track->GetIntegratedLength() > 365);
3866 
3867  if (! goodtrack)
3868  return kFALSE;
3869 
3870  if (! fAllowTOFmismatchFlag) {if (track->GetStatus() & AliESDtrack::kTOFmismatch) return kFALSE;}
3871 
3872  Bool_t statusMatchingHard = TPCTOFagree(track);
3873  if (fRequireStrictTOFTPCagreement && (!statusMatchingHard))
3874  return kFALSE;
3875 
3876  fBayesianResponse->ComputeProb(track,fCurrCentr); // fCurrCentr is needed for mismatch fraction
3877  Float_t *probabilities = fBayesianResponse->GetProb(); // Bayesian Probability (from 0 to 4) (Combined TPC || TOF) including a tuning of priors and TOF mismatch parameterization
3878 
3879  Float_t mismProb = fBayesianResponse->GetTOFMismProb(); // mismatch Bayesian probabilities
3880 
3881  fProbBayes = 0.0;
3882 
3883  switch (fParticleID)
3884  {
3885  case AliPID::kPion:
3886  fProbBayes = probabilities[2];
3887  break;
3888  case AliPID::kKaon:
3889  fProbBayes = probabilities[3];
3890  break;
3891  case AliPID::kProton:
3892  fProbBayes = probabilities[4];
3893  break;
3894  case AliPID::kElectron:
3895  fProbBayes = probabilities[0];
3896  break;
3897  case AliPID::kMuon:
3898  fProbBayes = probabilities[1];
3899  break;
3900  case AliPID::kDeuteron:
3901  fProbBayes = probabilities[5];
3902  break;
3903  case AliPID::kTriton:
3904  fProbBayes = probabilities[6];
3905  break;
3906  case AliPID::kHe3:
3907  fProbBayes = probabilities[7];
3908  break;
3909  default:
3910  return kFALSE;
3911  }
3912 
3913  // 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);
3914  if(fProbBayes > fParticleProbability && mismProb < 0.5){
3915  if(!fCutCharge)
3916  return kTRUE;
3917  else if (fCutCharge && fCharge * track->GetSign() > 0)
3918  return kTRUE;
3919  }
3920  return kFALSE;
3921 }
3922 
3923 
3924 //-----------------------------------------------------------------------
3925  // part added by Natasha
3926 Bool_t AliFlowTrackCuts::PassesNucleiSelection(const AliESDtrack* track)
3927 {
3928  //pid selection for heavy nuclei
3929  Bool_t select=kFALSE;
3930 
3931  //if (!track) continue;
3932 
3933  if (!track->GetInnerParam())
3934  return kFALSE; //break;
3935 
3936  const AliExternalTrackParam* tpcTrack = track->GetInnerParam();
3937 
3938  Double_t ptotTPC = tpcTrack->GetP();
3939  Double_t sigTPC = track->GetTPCsignal();
3940  Double_t dEdxBBA = 0.;
3941  Double_t dSigma = 0.;
3942 
3943  switch (fParticleID)
3944  {
3945  case AliPID::kDeuteron:
3946  //pid=10;
3947  dEdxBBA = AliExternalTrackParam::BetheBlochAleph(ptotTPC/1.8756,
3948  4.60e+00,
3949  8.9684e+00,
3950  1.640e-05,
3951  2.35e+00,
3952  2.35e+00);
3953  dSigma = (sigTPC - dEdxBBA)/dEdxBBA;
3954 
3955  if( ptotTPC<=1.1 && (dSigma < (0.5 - (0.1818*ptotTPC)) ) && (dSigma > ( (0.218*ptotTPC - 0.4) ) ) )
3956  {select=kTRUE;}
3957  break;
3958 
3959  case AliPID::kTriton:
3960  //pid=11;
3961  select=kFALSE;
3962  break;
3963 
3964  case AliPID::kHe3:
3965  //pid=12;
3966  // ----- Pass 2 -------
3967  dEdxBBA = 4.0 * AliExternalTrackParam::BetheBlochAleph( (2.*ptotTPC)/2.8084,
3968  1.74962,
3969  27.4992,
3970  4.00313e-15,
3971  2.42485,
3972  8.31768);
3973  dSigma = (sigTPC - dEdxBBA)/dEdxBBA;
3974  if(ptotTPC<=5.0 && (dSigma >= (-0.03968*ptotTPC - 0.1)) && (dSigma <= (0.31 - 0.0217*ptotTPC)))
3975  {select=kTRUE;}
3976  break;
3977 
3978  case AliPID::kAlpha:
3979  //pid=13;
3980  select=kFALSE;
3981  break;
3982 
3983  default:
3984  return kFALSE;
3985  }
3986 
3987  return select;
3988 }
3989 // end part added by Natasha
3990 //-----------------------------------------------------------------------
3991 Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCut(const AliAODTrack* track)
3992 {
3993  // do a simple combined cut on the n sigma from tpc and tof
3994  // with information of the pid response object (needs pid response task)
3995  // stub, not implemented yet
3996  if(!fPIDResponse) return kFALSE;
3997  if(!track) return kFALSE;
3998 
3999  // check TOF status
4000  if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kFALSE;
4001  if ((track->GetStatus()&AliVTrack::kTIME)==0) return kFALSE;
4002 
4003  // check TPC status
4004  if(track->GetTPCsignal() < 10) return kFALSE;
4005 
4006  Float_t nsigmaTPC = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC,track,fParticleID);
4007  Float_t nsigmaTOF = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF,track,fParticleID);
4008 
4009  Float_t nsigma2 = nsigmaTPC*nsigmaTPC + nsigmaTOF*nsigmaTOF;
4010 
4011  return (nsigma2 < fNsigmaCut2);
4012 
4013 }
4014 //-----------------------------------------------------------------------------
4015 Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaCut(const AliESDtrack* track)
4016 {
4017  // do a simple combined cut on the n sigma from tpc and tof
4018  // with information of the pid response object (needs pid response task)
4019  // stub, not implemented yet
4020  if(!fPIDResponse) return kFALSE;
4021  if(!track) return kFALSE;
4022 
4023  // check TOF status
4024  if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kFALSE;
4025  if ((track->GetStatus()&AliVTrack::kTIME)==0) return kFALSE;
4026 
4027  // check TPC status
4028  if(track->GetTPCsignal() < 10) return kFALSE;
4029 
4030  Float_t nsigmaTPC = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC,track,fParticleID);
4031  Float_t nsigmaTOF = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF,track,fParticleID);
4032 
4033  Float_t nsigma2 = nsigmaTPC*nsigmaTPC + nsigmaTOF*nsigmaTOF;
4034 
4035  return (nsigma2 < fNsigmaCut2);
4036 
4037 }
4038 //-----------------------------------------------------------------------------
4039 Bool_t AliFlowTrackCuts::PassesTPCTOFNsigmaPurityCut(const AliAODTrack* track)
4040 {
4041  // 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)
4042  // with information of the pid response object (needs pid response task)
4043 
4044  if(!fPIDResponse) return kFALSE;
4045  if(!track) return kFALSE;
4046 
4047  Int_t p_int = -999;
4048  Double_t pInterval=0;
4049  for(int i=0;i<60;i++){
4050  pInterval=0.1*i;
4051  if(track->P()>pInterval && track->P()<pInterval+0.1){p_int = i;}
4052  }
4053  /*
4054  Double_t LowPtPIDTPCnsigLow_Pion[2] = {-3,-3};
4055  Double_t LowPtPIDTPCnsigLow_Kaon[2] = {-3,-2};
4056  Double_t LowPtPIDTPCnsigHigh_Pion[2] ={2.4,3};
4057  Double_t LowPtPIDTPCnsigHigh_Kaon[2] ={3,2.2};
4058  */
4059 
4060  Float_t nsigmaTPC = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC,track,fParticleID);
4061  Float_t nsigmaTOF = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF,track,fParticleID);
4062 
4063  int index = (fParticleID-2)*60 + p_int;
4064  if ( (track->IsOn(AliAODTrack::kITSin))){
4065  if(p_int<2) return kFALSE;
4066 
4067  if(!fPurityFunction[index]){ cout<<"fPurityFunction[index] does not exist"<<endl; return kFALSE;}
4068  if(p_int>1){
4069  if((track->IsOn(AliAODTrack::kTOFpid))){
4070  if(fPurityFunction[index]->Eval(nsigmaTOF,nsigmaTPC)>fPurityLevel){
4071  if(TMath::Sqrt(TMath::Power(nsigmaTPC,2)+TMath::Power(nsigmaTOF,2))<3){
4072  return kTRUE;
4073  }
4074  }
4075  }
4076  }
4077  /*if(p_int<4){
4078  if(fParticleID==AliPID::kKaon && nsigmaTPC>LowPtPIDTPCnsigLow_Kaon[p_int-2] && nsigmaTPC<LowPtPIDTPCnsigHigh_Kaon[p_int-2]){return kTRUE;}
4079  if(fParticleID==AliPID::kPion && nsigmaTPC>LowPtPIDTPCnsigLow_Pion[p_int-2] && nsigmaTPC<LowPtPIDTPCnsigHigh_Pion[p_int-2]){return kTRUE;}
4080  if(fParticleID==AliPID::kProton && nsigmaTPC>-3 && nsigmaTPC<3){return kTRUE;}
4081  }*/
4082  }
4083  return kFALSE;
4084 }
4085 //-----------------------------------------------------------------------
4086 void AliFlowTrackCuts::SetPriors(Float_t centrCur){
4087  //set priors for the bayesian pid selection
4088  fCurrCentr = centrCur;
4089 
4090  fBinLimitPID[0] = 0.300000;
4091  fBinLimitPID[1] = 0.400000;
4092  fBinLimitPID[2] = 0.500000;
4093  fBinLimitPID[3] = 0.600000;
4094  fBinLimitPID[4] = 0.700000;
4095  fBinLimitPID[5] = 0.800000;
4096  fBinLimitPID[6] = 0.900000;
4097  fBinLimitPID[7] = 1.000000;
4098  fBinLimitPID[8] = 1.200000;
4099  fBinLimitPID[9] = 1.400000;
4100  fBinLimitPID[10] = 1.600000;
4101  fBinLimitPID[11] = 1.800000;
4102  fBinLimitPID[12] = 2.000000;
4103  fBinLimitPID[13] = 2.200000;
4104  fBinLimitPID[14] = 2.400000;
4105  fBinLimitPID[15] = 2.600000;
4106  fBinLimitPID[16] = 2.800000;
4107  fBinLimitPID[17] = 3.000000;
4108 
4109  // 0-10%
4110  if(centrCur < 10){
4111  fC[0][0] = 0.005;
4112  fC[0][1] = 0.005;
4113  fC[0][2] = 1.0000;
4114  fC[0][3] = 0.0010;
4115  fC[0][4] = 0.0010;
4116 
4117  fC[1][0] = 0.005;
4118  fC[1][1] = 0.005;
4119  fC[1][2] = 1.0000;
4120  fC[1][3] = 0.0168;
4121  fC[1][4] = 0.0122;
4122 
4123  fC[2][0] = 0.005;
4124  fC[2][1] = 0.005;
4125  fC[2][2] = 1.0000;
4126  fC[2][3] = 0.0272;
4127  fC[2][4] = 0.0070;
4128 
4129  fC[3][0] = 0.005;
4130  fC[3][1] = 0.005;
4131  fC[3][2] = 1.0000;
4132  fC[3][3] = 0.0562;
4133  fC[3][4] = 0.0258;
4134 
4135  fC[4][0] = 0.005;
4136  fC[4][1] = 0.005;
4137  fC[4][2] = 1.0000;
4138  fC[4][3] = 0.0861;
4139  fC[4][4] = 0.0496;
4140 
4141  fC[5][0] = 0.005;
4142  fC[5][1] = 0.005;
4143  fC[5][2] = 1.0000;
4144  fC[5][3] = 0.1168;
4145  fC[5][4] = 0.0740;
4146 
4147  fC[6][0] = 0.005;
4148  fC[6][1] = 0.005;
4149  fC[6][2] = 1.0000;
4150  fC[6][3] = 0.1476;
4151  fC[6][4] = 0.0998;
4152 
4153  fC[7][0] = 0.005;
4154  fC[7][1] = 0.005;
4155  fC[7][2] = 1.0000;
4156  fC[7][3] = 0.1810;
4157  fC[7][4] = 0.1296;
4158 
4159  fC[8][0] = 0.005;
4160  fC[8][1] = 0.005;
4161  fC[8][2] = 1.0000;
4162  fC[8][3] = 0.2240;
4163  fC[8][4] = 0.1827;
4164 
4165  fC[9][0] = 0.005;
4166  fC[9][1] = 0.005;
4167  fC[9][2] = 1.0000;
4168  fC[9][3] = 0.2812;
4169  fC[9][4] = 0.2699;
4170 
4171  fC[10][0] = 0.005;
4172  fC[10][1] = 0.005;
4173  fC[10][2] = 1.0000;
4174  fC[10][3] = 0.3328;
4175  fC[10][4] = 0.3714;
4176 
4177  fC[11][0] = 0.005;
4178  fC[11][1] = 0.005;
4179  fC[11][2] = 1.0000;
4180  fC[11][3] = 0.3780;
4181  fC[11][4] = 0.4810;
4182 
4183  fC[12][0] = 0.005;
4184  fC[12][1] = 0.005;
4185  fC[12][2] = 1.0000;
4186  fC[12][3] = 0.4125;
4187  fC[12][4] = 0.5771;
4188 
4189  fC[13][0] = 0.005;
4190  fC[13][1] = 0.005;
4191  fC[13][2] = 1.0000;
4192  fC[13][3] = 0.4486;
4193  fC[13][4] = 0.6799;
4194 
4195  fC[14][0] = 0.005;
4196  fC[14][1] = 0.005;
4197  fC[14][2] = 1.0000;
4198  fC[14][3] = 0.4840;
4199  fC[14][4] = 0.7668;
4200 
4201  fC[15][0] = 0.005;
4202  fC[15][1] = 0.005;
4203  fC[15][2] = 1.0000;
4204  fC[15][3] = 0.4971;
4205  fC[15][4] = 0.8288;
4206 
4207  fC[16][0] = 0.005;
4208  fC[16][1] = 0.005;
4209  fC[16][2] = 1.0000;
4210  fC[16][3] = 0.4956;
4211  fC[16][4] = 0.8653;
4212 
4213  fC[17][0] = 0.005;
4214  fC[17][1] = 0.005;
4215  fC[17][2] = 1.0000;
4216  fC[17][3] = 0.5173;
4217  fC[17][4] = 0.9059;
4218  }
4219  // 10-20%
4220  else if(centrCur < 20){
4221  fC[0][0] = 0.005;
4222  fC[0][1] = 0.005;
4223  fC[0][2] = 1.0000;
4224  fC[0][3] = 0.0010;
4225  fC[0][4] = 0.0010;
4226 
4227  fC[1][0] = 0.005;
4228  fC[1][1] = 0.005;
4229  fC[1][2] = 1.0000;
4230  fC[1][3] = 0.0132;
4231  fC[1][4] = 0.0088;
4232 
4233  fC[2][0] = 0.005;
4234  fC[2][1] = 0.005;
4235  fC[2][2] = 1.0000;
4236  fC[2][3] = 0.0283;
4237  fC[2][4] = 0.0068;
4238 
4239  fC[3][0] = 0.005;
4240  fC[3][1] = 0.005;
4241  fC[3][2] = 1.0000;
4242  fC[3][3] = 0.0577;
4243  fC[3][4] = 0.0279;
4244 
4245  fC[4][0] = 0.005;
4246  fC[4][1] = 0.005;
4247  fC[4][2] = 1.0000;
4248  fC[4][3] = 0.0884;
4249  fC[4][4] = 0.0534;
4250 
4251  fC[5][0] = 0.005;
4252  fC[5][1] = 0.005;
4253  fC[5][2] = 1.0000;
4254  fC[5][3] = 0.1179;
4255  fC[5][4] = 0.0794;
4256 
4257  fC[6][0] = 0.005;
4258  fC[6][1] = 0.005;
4259  fC[6][2] = 1.0000;
4260  fC[6][3] = 0.1480;
4261  fC[6][4] = 0.1058;
4262 
4263  fC[7][0] = 0.005;
4264  fC[7][1] = 0.005;
4265  fC[7][2] = 1.0000;
4266  fC[7][3] = 0.1807;
4267  fC[7][4] = 0.1366;
4268 
4269  fC[8][0] = 0.005;
4270  fC[8][1] = 0.005;
4271  fC[8][2] = 1.0000;
4272  fC[8][3] = 0.2219;
4273  fC[8][4] = 0.1891;
4274 
4275  fC[9][0] = 0.005;
4276  fC[9][1] = 0.005;
4277  fC[9][2] = 1.0000;
4278  fC[9][3] = 0.2804;
4279  fC[9][4] = 0.2730;
4280 
4281  fC[10][0] = 0.005;
4282  fC[10][1] = 0.005;
4283  fC[10][2] = 1.0000;
4284  fC[10][3] = 0.3283;
4285  fC[10][4] = 0.3660;
4286 
4287  fC[11][0] = 0.005;
4288  fC[11][1] = 0.005;
4289  fC[11][2] = 1.0000;
4290  fC[11][3] = 0.3710;
4291  fC[11][4] = 0.4647;
4292 
4293  fC[12][0] = 0.005;
4294  fC[12][1] = 0.005;
4295  fC[12][2] = 1.0000;
4296  fC[12][3] = 0.4093;
4297  fC[12][4] = 0.5566;
4298 
4299  fC[13][0] = 0.005;
4300  fC[13][1] = 0.005;
4301  fC[13][2] = 1.0000;
4302  fC[13][3] = 0.4302;
4303  fC[13][4] = 0.6410;
4304 
4305  fC[14][0] = 0.005;
4306  fC[14][1] = 0.005;
4307  fC[14][2] = 1.0000;
4308  fC[14][3] = 0.4649;
4309  fC[14][4] = 0.7055;
4310 
4311  fC[15][0] = 0.005;
4312  fC[15][1] = 0.005;
4313  fC[15][2] = 1.0000;
4314  fC[15][3] = 0.4523;
4315  fC[15][4] = 0.7440;
4316 
4317  fC[16][0] = 0.005;
4318  fC[16][1] = 0.005;
4319  fC[16][2] = 1.0000;
4320  fC[16][3] = 0.4591;
4321  fC[16][4] = 0.7799;
4322 
4323  fC[17][0] = 0.005;
4324  fC[17][1] = 0.005;
4325  fC[17][2] = 1.0000;
4326  fC[17][3] = 0.4804;
4327  fC[17][4] = 0.8218;
4328  }
4329  // 20-30%
4330  else if(centrCur < 30){
4331  fC[0][0] = 0.005;
4332  fC[0][1] = 0.005;
4333  fC[0][2] = 1.0000;
4334  fC[0][3] = 0.0010;
4335  fC[0][4] = 0.0010;
4336 
4337  fC[1][0] = 0.005;
4338  fC[1][1] = 0.005;
4339  fC[1][2] = 1.0000;
4340  fC[1][3] = 0.0102;
4341  fC[1][4] = 0.0064;
4342 
4343  fC[2][0] = 0.005;
4344  fC[2][1] = 0.005;
4345  fC[2][2] = 1.0000;
4346  fC[2][3] = 0.0292;
4347  fC[2][4] = 0.0066;
4348 
4349  fC[3][0] = 0.005;
4350  fC[3][1] = 0.005;
4351  fC[3][2] = 1.0000;
4352  fC[3][3] = 0.0597;
4353  fC[3][4] = 0.0296;
4354 
4355  fC[4][0] = 0.005;
4356  fC[4][1] = 0.005;
4357  fC[4][2] = 1.0000;
4358  fC[4][3] = 0.0900;
4359  fC[4][4] = 0.0589;
4360 
4361  fC[5][0] = 0.005;
4362  fC[5][1] = 0.005;
4363  fC[5][2] = 1.0000;
4364  fC[5][3] = 0.1199;
4365  fC[5][4] = 0.0859;
4366 
4367  fC[6][0] = 0.005;
4368  fC[6][1] = 0.005;
4369  fC[6][2] = 1.0000;
4370  fC[6][3] = 0.1505;
4371  fC[6][4] = 0.1141;
4372 
4373  fC[7][0] = 0.005;
4374  fC[7][1] = 0.005;
4375  fC[7][2] = 1.0000;
4376  fC[7][3] = 0.1805;
4377  fC[7][4] = 0.1454;
4378 
4379  fC[8][0] = 0.005;
4380  fC[8][1] = 0.005;
4381  fC[8][2] = 1.0000;
4382  fC[8][3] = 0.2221;
4383  fC[8][4] = 0.2004;
4384 
4385  fC[9][0] = 0.005;
4386  fC[9][1] = 0.005;
4387  fC[9][2] = 1.0000;
4388  fC[9][3] = 0.2796;
4389  fC[9][4] = 0.2838;
4390 
4391  fC[10][0] = 0.005;
4392  fC[10][1] = 0.005;
4393  fC[10][2] = 1.0000;
4394  fC[10][3] = 0.3271;
4395  fC[10][4] = 0.3682;
4396 
4397  fC[11][0] = 0.005;
4398  fC[11][1] = 0.005;
4399  fC[11][2] = 1.0000;
4400  fC[11][3] = 0.3648;
4401  fC[11][4] = 0.4509;
4402 
4403  fC[12][0] = 0.005;
4404  fC[12][1] = 0.005;
4405  fC[12][2] = 1.0000;
4406  fC[12][3] = 0.3988;
4407  fC[12][4] = 0.5339;
4408 
4409  fC[13][0] = 0.005;
4410  fC[13][1] = 0.005;
4411  fC[13][2] = 1.0000;
4412  fC[13][3] = 0.4315;
4413  fC[13][4] = 0.5995;
4414 
4415  fC[14][0] = 0.005;
4416  fC[14][1] = 0.005;
4417  fC[14][2] = 1.0000;
4418  fC[14][3] = 0.4548;
4419  fC[14][4] = 0.6612;
4420 
4421  fC[15][0] = 0.005;
4422  fC[15][1] = 0.005;
4423  fC[15][2] = 1.0000;
4424  fC[15][3] = 0.4744;
4425  fC[15][4] = 0.7060;
4426 
4427  fC[16][0] = 0.005;
4428  fC[16][1] = 0.005;
4429  fC[16][2] = 1.0000;
4430  fC[16][3] = 0.4899;
4431  fC[16][4] = 0.7388;
4432 
4433  fC[17][0] = 0.005;
4434  fC[17][1] = 0.005;
4435  fC[17][2] = 1.0000;
4436  fC[17][3] = 0.4411;
4437  fC[17][4] = 0.7293;
4438  }
4439  // 30-40%
4440  else if(centrCur < 40){
4441  fC[0][0] = 0.005;
4442  fC[0][1] = 0.005;
4443  fC[0][2] = 1.0000;
4444  fC[0][3] = 0.0010;
4445  fC[0][4] = 0.0010;
4446 
4447  fC[1][0] = 0.005;
4448  fC[1][1] = 0.005;
4449  fC[1][2] = 1.0000;
4450  fC[1][3] = 0.0102;
4451  fC[1][4] = 0.0048;
4452 
4453  fC[2][0] = 0.005;
4454  fC[2][1] = 0.005;
4455  fC[2][2] = 1.0000;
4456  fC[2][3] = 0.0306;
4457  fC[2][4] = 0.0079;
4458 
4459  fC[3][0] = 0.005;
4460  fC[3][1] = 0.005;
4461  fC[3][2] = 1.0000;
4462  fC[3][3] = 0.0617;
4463  fC[3][4] = 0.0338;
4464 
4465  fC[4][0] = 0.005;
4466  fC[4][1] = 0.005;
4467  fC[4][2] = 1.0000;
4468  fC[4][3] = 0.0920;
4469  fC[4][4] = 0.0652;
4470 
4471  fC[5][0] = 0.005;
4472  fC[5][1] = 0.005;
4473  fC[5][2] = 1.0000;
4474  fC[5][3] = 0.1211;
4475  fC[5][4] = 0.0955;
4476 
4477  fC[6][0] = 0.005;
4478  fC[6][1] = 0.005;
4479  fC[6][2] = 1.0000;
4480  fC[6][3] = 0.1496;
4481  fC[6][4] = 0.1242;
4482 
4483  fC[7][0] = 0.005;
4484  fC[7][1] = 0.005;
4485  fC[7][2] = 1.0000;
4486  fC[7][3] = 0.1807;
4487  fC[7][4] = 0.1576;
4488 
4489  fC[8][0] = 0.005;
4490  fC[8][1] = 0.005;
4491  fC[8][2] = 1.0000;
4492  fC[8][3] = 0.2195;
4493  fC[8][4] = 0.2097;
4494 
4495  fC[9][0] = 0.005;
4496  fC[9][1] = 0.005;
4497  fC[9][2] = 1.0000;
4498  fC[9][3] = 0.2732;
4499  fC[9][4] = 0.2884;
4500 
4501  fC[10][0] = 0.005;
4502  fC[10][1] = 0.005;
4503  fC[10][2] = 1.0000;
4504  fC[10][3] = 0.3204;
4505  fC[10][4] = 0.3679;
4506 
4507  fC[11][0] = 0.005;
4508  fC[11][1] = 0.005;
4509  fC[11][2] = 1.0000;
4510  fC[11][3] = 0.3564;
4511  fC[11][4] = 0.4449;
4512 
4513  fC[12][0] = 0.005;
4514  fC[12][1] = 0.005;
4515  fC[12][2] = 1.0000;
4516  fC[12][3] = 0.3791;
4517  fC[12][4] = 0.5052;
4518 
4519  fC[13][0] = 0.005;
4520  fC[13][1] = 0.005;
4521  fC[13][2] = 1.0000;
4522  fC[13][3] = 0.4062;
4523  fC[13][4] = 0.5647;
4524 
4525  fC[14][0] = 0.005;
4526  fC[14][1] = 0.005;
4527  fC[14][2] = 1.0000;
4528  fC[14][3] = 0.4234;
4529  fC[14][4] = 0.6203;
4530 
4531  fC[15][0] = 0.005;
4532  fC[15][1] = 0.005;
4533  fC[15][2] = 1.0000;
4534  fC[15][3] = 0.4441;
4535  fC[15][4] = 0.6381;
4536 
4537  fC[16][0] = 0.005;
4538  fC[16][1] = 0.005;
4539  fC[16][2] = 1.0000;
4540  fC[16][3] = 0.4629;
4541  fC[16][4] = 0.6496;
4542 
4543  fC[17][0] = 0.005;
4544  fC[17][1] = 0.005;
4545  fC[17][2] = 1.0000;
4546  fC[17][3] = 0.4293;
4547  fC[17][4] = 0.6491;
4548  }
4549  // 40-50%
4550  else if(centrCur < 50){
4551  fC[0][0] = 0.005;
4552  fC[0][1] = 0.005;
4553  fC[0][2] = 1.0000;
4554  fC[0][3] = 0.0010;
4555  fC[0][4] = 0.0010;
4556 
4557  fC[1][0] = 0.005;
4558  fC[1][1] = 0.005;
4559  fC[1][2] = 1.0000;
4560  fC[1][3] = 0.0093;
4561  fC[1][4] = 0.0057;
4562 
4563  fC[2][0] = 0.005;
4564  fC[2][1] = 0.005;
4565  fC[2][2] = 1.0000;
4566  fC[2][3] = 0.0319;
4567  fC[2][4] = 0.0075;
4568 
4569  fC[3][0] = 0.005;
4570  fC[3][1] = 0.005;
4571  fC[3][2] = 1.0000;
4572  fC[3][3] = 0.0639;
4573  fC[3][4] = 0.0371;
4574 
4575  fC[4][0] = 0.005;
4576  fC[4][1] = 0.005;
4577  fC[4][2] = 1.0000;
4578  fC[4][3] = 0.0939;
4579  fC[4][4] = 0.0725;
4580 
4581  fC[5][0] = 0.005;
4582  fC[5][1] = 0.005;
4583  fC[5][2] = 1.0000;
4584  fC[5][3] = 0.1224;
4585  fC[5][4] = 0.1045;
4586 
4587  fC[6][0] = 0.005;
4588  fC[6][1] = 0.005;
4589  fC[6][2] = 1.0000;
4590  fC[6][3] = 0.1520;
4591  fC[6][4] = 0.1387;
4592 
4593  fC[7][0] = 0.005;
4594  fC[7][1] = 0.005;
4595  fC[7][2] = 1.0000;
4596  fC[7][3] = 0.1783;
4597  fC[7][4] = 0.1711;
4598 
4599  fC[8][0] = 0.005;
4600  fC[8][1] = 0.005;
4601  fC[8][2] = 1.0000;
4602  fC[8][3] = 0.2202;
4603  fC[8][4] = 0.2269;
4604 
4605  fC[9][0] = 0.005;
4606  fC[9][1] = 0.005;
4607  fC[9][2] = 1.0000;
4608  fC[9][3] = 0.2672;
4609  fC[9][4] = 0.2955;
4610 
4611  fC[10][0] = 0.005;
4612  fC[10][1] = 0.005;
4613  fC[10][2] = 1.0000;
4614  fC[10][3] = 0.3191;
4615  fC[10][4] = 0.3676;
4616 
4617  fC[11][0] = 0.005;
4618  fC[11][1] = 0.005;
4619  fC[11][2] = 1.0000;
4620  fC[11][3] = 0.3434;
4621  fC[11][4] = 0.4321;
4622 
4623  fC[12][0] = 0.005;
4624  fC[12][1] = 0.005;
4625  fC[12][2] = 1.0000;
4626  fC[12][3] = 0.3692;
4627  fC[12][4] = 0.4879;
4628 
4629  fC[13][0] = 0.005;
4630  fC[13][1] = 0.005;
4631  fC[13][2] = 1.0000;
4632  fC[13][3] = 0.3993;
4633  fC[13][4] = 0.5377;
4634 
4635  fC[14][0] = 0.005;
4636  fC[14][1] = 0.005;
4637  fC[14][2] = 1.0000;
4638  fC[14][3] = 0.3818;
4639  fC[14][4] = 0.5547;
4640 
4641  fC[15][0] = 0.005;
4642  fC[15][1] = 0.005;
4643  fC[15][2] = 1.0000;
4644  fC[15][3] = 0.4003;
4645  fC[15][4] = 0.5484;
4646 
4647  fC[16][0] = 0.005;
4648  fC[16][1] = 0.005;
4649  fC[16][2] = 1.0000;
4650  fC[16][3] = 0.4281;
4651  fC[16][4] = 0.5383;
4652 
4653  fC[17][0] = 0.005;
4654  fC[17][1] = 0.005;
4655  fC[17][2] = 1.0000;
4656  fC[17][3] = 0.3960;
4657  fC[17][4] = 0.5374;
4658  }
4659  // 50-60%
4660  else if(centrCur < 60){
4661  fC[0][0] = 0.005;
4662  fC[0][1] = 0.005;
4663  fC[0][2] = 1.0000;
4664  fC[0][3] = 0.0010;
4665  fC[0][4] = 0.0010;
4666 
4667  fC[1][0] = 0.005;
4668  fC[1][1] = 0.005;
4669  fC[1][2] = 1.0000;
4670  fC[1][3] = 0.0076;
4671  fC[1][4] = 0.0032;
4672 
4673  fC[2][0] = 0.005;
4674  fC[2][1] = 0.005;
4675  fC[2][2] = 1.0000;
4676  fC[2][3] = 0.0329;
4677  fC[2][4] = 0.0085;
4678 
4679  fC[3][0] = 0.005;
4680  fC[3][1] = 0.005;
4681  fC[3][2] = 1.0000;
4682  fC[3][3] = 0.0653;
4683  fC[3][4] = 0.0423;
4684 
4685  fC[4][0] = 0.005;
4686  fC[4][1] = 0.005;
4687  fC[4][2] = 1.0000;
4688  fC[4][3] = 0.0923;
4689  fC[4][4] = 0.0813;
4690 
4691  fC[5][0] = 0.005;
4692  fC[5][1] = 0.005;
4693  fC[5][2] = 1.0000;
4694  fC[5][3] = 0.1219;
4695  fC[5][4] = 0.1161;
4696 
4697  fC[6][0] = 0.005;
4698  fC[6][1] = 0.005;
4699  fC[6][2] = 1.0000;
4700  fC[6][3] = 0.1519;
4701  fC[6][4] = 0.1520;
4702 
4703  fC[7][0] = 0.005;
4704  fC[7][1] = 0.005;
4705  fC[7][2] = 1.0000;
4706  fC[7][3] = 0.1763;
4707  fC[7][4] = 0.1858;
4708 
4709  fC[8][0] = 0.005;
4710  fC[8][1] = 0.005;
4711  fC[8][2] = 1.0000;
4712  fC[8][3] = 0.2178;
4713  fC[8][4] = 0.2385;
4714 
4715  fC[9][0] = 0.005;
4716  fC[9][1] = 0.005;
4717  fC[9][2] = 1.0000;
4718  fC[9][3] = 0.2618;
4719  fC[9][4] = 0.3070;
4720 
4721  fC[10][0] = 0.005;
4722  fC[10][1] = 0.005;
4723  fC[10][2] = 1.0000;
4724  fC[10][3] = 0.3067;
4725  fC[10][4] = 0.3625;
4726 
4727  fC[11][0] = 0.005;
4728  fC[11][1] = 0.005;
4729  fC[11][2] = 1.0000;
4730  fC[11][3] = 0.3336;
4731  fC[11][4] = 0.4188;
4732 
4733  fC[12][0] = 0.005;
4734  fC[12][1] = 0.005;
4735  fC[12][2] = 1.0000;
4736  fC[12][3] = 0.3706;
4737  fC[12][4] = 0.4511;
4738 
4739  fC[13][0] = 0.005;
4740  fC[13][1] = 0.005;
4741  fC[13][2] = 1.0000;
4742  fC[13][3] = 0.3765;
4743  fC[13][4] = 0.4729;
4744 
4745  fC[14][0] = 0.005;
4746  fC[14][1] = 0.005;
4747  fC[14][2] = 1.0000;
4748  fC[14][3] = 0.3942;
4749  fC[14][4] = 0.4855;
4750 
4751  fC[15][0] = 0.005;
4752  fC[15][1] = 0.005;
4753  fC[15][2] = 1.0000;
4754  fC[15][3] = 0.4051;
4755  fC[15][4] = 0.4762;
4756 
4757  fC[16][0] = 0.005;
4758  fC[16][1] = 0.005;
4759  fC[16][2] = 1.0000;
4760  fC[16][3] = 0.3843;
4761  fC[16][4] = 0.4763;
4762 
4763  fC[17][0] = 0.005;
4764  fC[17][1] = 0.005;
4765  fC[17][2] = 1.0000;
4766  fC[17][3] = 0.4237;
4767  fC[17][4] = 0.4773;
4768  }
4769  // 60-70%
4770  else if(centrCur < 70){
4771  fC[0][0] = 0.005;
4772  fC[0][1] = 0.005;
4773  fC[0][2] = 1.0000;
4774  fC[0][3] = 0.0010;
4775  fC[0][4] = 0.0010;
4776 
4777  fC[1][0] = 0.005;
4778  fC[1][1] = 0.005;
4779  fC[1][2] = 1.0000;
4780  fC[1][3] = 0.0071;
4781  fC[1][4] = 0.0012;
4782 
4783  fC[2][0] = 0.005;
4784  fC[2][1] = 0.005;
4785  fC[2][2] = 1.0000;
4786  fC[2][3] = 0.0336;
4787  fC[2][4] = 0.0097;
4788 
4789  fC[3][0] = 0.005;
4790  fC[3][1] = 0.005;
4791  fC[3][2] = 1.0000;
4792  fC[3][3] = 0.0662;
4793  fC[3][4] = 0.0460;
4794 
4795  fC[4][0] = 0.005;
4796  fC[4][1] = 0.005;
4797  fC[4][2] = 1.0000;
4798  fC[4][3] = 0.0954;
4799  fC[4][4] = 0.0902;
4800 
4801  fC[5][0] = 0.005;
4802  fC[5][1] = 0.005;
4803  fC[5][2] = 1.0000;
4804  fC[5][3] = 0.1181;
4805  fC[5][4] = 0.1306;
4806 
4807  fC[6][0] = 0.005;
4808  fC[6][1] = 0.005;
4809  fC[6][2] = 1.0000;
4810  fC[6][3] = 0.1481;
4811  fC[6][4] = 0.1662;
4812 
4813  fC[7][0] = 0.005;
4814  fC[7][1] = 0.005;
4815  fC[7][2] = 1.0000;
4816  fC[7][3] = 0.1765;
4817  fC[7][4] = 0.1963;
4818 
4819  fC[8][0] = 0.005;
4820  fC[8][1] = 0.005;
4821  fC[8][2] = 1.0000;
4822  fC[8][3] = 0.2155;
4823  fC[8][4] = 0.2433;
4824 
4825  fC[9][0] = 0.005;
4826  fC[9][1] = 0.005;
4827  fC[9][2] = 1.0000;
4828  fC[9][3] = 0.2580;
4829  fC[9][4] = 0.3022;
4830 
4831  fC[10][0] = 0.005;
4832  fC[10][1] = 0.005;
4833  fC[10][2] = 1.0000;
4834  fC[10][3] = 0.2872;
4835  fC[10][4] = 0.3481;
4836 
4837  fC[11][0] = 0.005;
4838  fC[11][1] = 0.005;
4839  fC[11][2] = 1.0000;
4840  fC[11][3] = 0.3170;
4841  fC[11][4] = 0.3847;
4842 
4843  fC[12][0] = 0.005;
4844  fC[12][1] = 0.005;
4845  fC[12][2] = 1.0000;
4846  fC[12][3] = 0.3454;
4847  fC[12][4] = 0.4258;
4848 
4849  fC[13][0] = 0.005;
4850  fC[13][1] = 0.005;
4851  fC[13][2] = 1.0000;
4852  fC[13][3] = 0.3580;
4853  fC[13][4] = 0.4299;
4854 
4855  fC[14][0] = 0.005;
4856  fC[14][1] = 0.005;
4857  fC[14][2] = 1.0000;
4858  fC[14][3] = 0.3903;
4859  fC[14][4] = 0.4326;
4860 
4861  fC[15][0] = 0.005;
4862  fC[15][1] = 0.005;
4863  fC[15][2] = 1.0000;
4864  fC[15][3] = 0.3690;
4865  fC[15][4] = 0.4491;
4866 
4867  fC[16][0] = 0.005;
4868  fC[16][1] = 0.005;
4869  fC[16][2] = 1.0000;
4870  fC[16][3] = 0.4716;
4871  fC[16][4] = 0.4298;
4872 
4873  fC[17][0] = 0.005;
4874  fC[17][1] = 0.005;
4875  fC[17][2] = 1.0000;
4876  fC[17][3] = 0.3875;
4877  fC[17][4] = 0.4083;
4878  }
4879  // 70-80%
4880  else if(centrCur < 80){
4881  fC[0][0] = 0.005;
4882  fC[0][1] = 0.005;
4883  fC[0][2] = 1.0000;
4884  fC[0][3] = 0.0010;
4885  fC[0][4] = 0.0010;
4886 
4887  fC[1][0] = 0.005;
4888  fC[1][1] = 0.005;
4889  fC[1][2] = 1.0000;
4890  fC[1][3] = 0.0075;
4891  fC[1][4] = 0.0007;
4892 
4893  fC[2][0] = 0.005;
4894  fC[2][1] = 0.005;
4895  fC[2][2] = 1.0000;
4896  fC[2][3] = 0.0313;
4897  fC[2][4] = 0.0124;
4898 
4899  fC[3][0] = 0.005;
4900  fC[3][1] = 0.005;
4901  fC[3][2] = 1.0000;
4902  fC[3][3] = 0.0640;
4903  fC[3][4] = 0.0539;
4904 
4905  fC[4][0] = 0.005;
4906  fC[4][1] = 0.005;
4907  fC[4][2] = 1.0000;
4908  fC[4][3] = 0.0923;
4909  fC[4][4] = 0.0992;
4910 
4911  fC[5][0] = 0.005;
4912  fC[5][1] = 0.005;
4913  fC[5][2] = 1.0000;
4914  fC[5][3] = 0.1202;
4915  fC[5][4] = 0.1417;
4916 
4917  fC[6][0] = 0.005;
4918  fC[6][1] = 0.005;
4919  fC[6][2] = 1.0000;
4920  fC[6][3] = 0.1413;
4921  fC[6][4] = 0.1729;
4922 
4923  fC[7][0] = 0.005;
4924  fC[7][1] = 0.005;
4925  fC[7][2] = 1.0000;
4926  fC[7][3] = 0.1705;
4927  fC[7][4] = 0.1999;
4928 
4929  fC[8][0] = 0.005;
4930  fC[8][1] = 0.005;
4931  fC[8][2] = 1.0000;
4932  fC[8][3] = 0.2103;
4933  fC[8][4] = 0.2472;
4934 
4935  fC[9][0] = 0.005;
4936  fC[9][1] = 0.005;
4937  fC[9][2] = 1.0000;
4938  fC[9][3] = 0.2373;
4939  fC[9][4] = 0.2916;
4940 
4941  fC[10][0] = 0.005;
4942  fC[10][1] = 0.005;
4943  fC[10][2] = 1.0000;
4944  fC[10][3] = 0.2824;
4945  fC[10][4] = 0.3323;
4946 
4947  fC[11][0] = 0.005;
4948  fC[11][1] = 0.005;
4949  fC[11][2] = 1.0000;
4950  fC[11][3] = 0.3046;
4951  fC[11][4] = 0.3576;
4952 
4953  fC[12][0] = 0.005;
4954  fC[12][1] = 0.005;
4955  fC[12][2] = 1.0000;
4956  fC[12][3] = 0.3585;
4957  fC[12][4] = 0.4003;
4958 
4959  fC[13][0] = 0.005;
4960  fC[13][1] = 0.005;
4961  fC[13][2] = 1.0000;
4962  fC[13][3] = 0.3461;
4963  fC[13][4] = 0.3982;
4964 
4965  fC[14][0] = 0.005;
4966  fC[14][1] = 0.005;
4967  fC[14][2] = 1.0000;
4968  fC[14][3] = 0.3362;
4969  fC[14][4] = 0.3776;
4970 
4971  fC[15][0] = 0.005;
4972  fC[15][1] = 0.005;
4973  fC[15][2] = 1.0000;
4974  fC[15][3] = 0.3071;
4975  fC[15][4] = 0.3500;
4976 
4977  fC[16][0] = 0.005;
4978  fC[16][1] = 0.005;
4979  fC[16][2] = 1.0000;
4980  fC[16][3] = 0.2914;
4981  fC[16][4] = 0.3937;
4982 
4983  fC[17][0] = 0.005;
4984  fC[17][1] = 0.005;
4985  fC[17][2] = 1.0000;
4986  fC[17][3] = 0.3727;
4987  fC[17][4] = 0.3877;
4988  }
4989  // 80-100%
4990  else{
4991  fC[0][0] = 0.005;
4992  fC[0][1] = 0.005;
4993  fC[0][2] = 1.0000;
4994  fC[0][3] = 0.0010;
4995  fC[0][4] = 0.0010;
4996 
4997  fC[1][0] = 0.005;
4998  fC[1][1] = 0.005;
4999  fC[1][2] = 1.0000;
5000  fC[1][3] = 0.0060;
5001  fC[1][4] = 0.0035;
5002 
5003  fC[2][0] = 0.005;
5004  fC[2][1] = 0.005;
5005  fC[2][2] = 1.0000;
5006  fC[2][3] = 0.0323;
5007  fC[2][4] = 0.0113;
5008 
5009  fC[3][0] = 0.005;
5010  fC[3][1] = 0.005;
5011  fC[3][2] = 1.0000;
5012  fC[3][3] = 0.0609;
5013  fC[3][4] = 0.0653;
5014 
5015  fC[4][0] = 0.005;
5016  fC[4][1] = 0.005;
5017  fC[4][2] = 1.0000;
5018  fC[4][3] = 0.0922;
5019  fC[4][4] = 0.1076;
5020 
5021  fC[5][0] = 0.005;
5022  fC[5][1] = 0.005;
5023  fC[5][2] = 1.0000;
5024  fC[5][3] = 0.1096;
5025  fC[5][4] = 0.1328;
5026 
5027  fC[6][0] = 0.005;
5028  fC[6][1] = 0.005;
5029  fC[6][2] = 1.0000;
5030  fC[6][3] = 0.1495;
5031  fC[6][4] = 0.1779;
5032 
5033  fC[7][0] = 0.005;
5034  fC[7][1] = 0.005;
5035  fC[7][2] = 1.0000;
5036  fC[7][3] = 0.1519;
5037  fC[7][4] = 0.1989;
5038 
5039  fC[8][0] = 0.005;
5040  fC[8][1] = 0.005;
5041  fC[8][2] = 1.0000;
5042  fC[8][3] = 0.1817;
5043  fC[8][4] = 0.2472;
5044 
5045  fC[9][0] = 0.005;
5046  fC[9][1] = 0.005;
5047  fC[9][2] = 1.0000;
5048  fC[9][3] = 0.2429;
5049  fC[9][4] = 0.2684;
5050 
5051  fC[10][0] = 0.005;
5052  fC[10][1] = 0.005;
5053  fC[10][2] = 1.0000;
5054  fC[10][3] = 0.2760;
5055  fC[10][4] = 0.3098;
5056 
5057  fC[11][0] = 0.005;
5058  fC[11][1] = 0.005;
5059  fC[11][2] = 1.0000;
5060  fC[11][3] = 0.2673;
5061  fC[11][4] = 0.3198;
5062 
5063  fC[12][0] = 0.005;
5064  fC[12][1] = 0.005;
5065  fC[12][2] = 1.0000;
5066  fC[12][3] = 0.3165;
5067  fC[12][4] = 0.3564;
5068 
5069  fC[13][0] = 0.005;
5070  fC[13][1] = 0.005;
5071  fC[13][2] = 1.0000;
5072  fC[13][3] = 0.3526;
5073  fC[13][4] = 0.3011;
5074 
5075  fC[14][0] = 0.005;
5076  fC[14][1] = 0.005;
5077  fC[14][2] = 1.0000;
5078  fC[14][3] = 0.3788;
5079  fC[14][4] = 0.3011;
5080 
5081  fC[15][0] = 0.005;
5082  fC[15][1] = 0.005;
5083  fC[15][2] = 1.0000;
5084  fC[15][3] = 0.3788;
5085  fC[15][4] = 0.3011;
5086 
5087  fC[16][0] = 0.005;
5088  fC[16][1] = 0.005;
5089  fC[16][2] = 1.0000;
5090  fC[16][3] = 0.3788;
5091  fC[16][4] = 0.3011;
5092 
5093  fC[17][0] = 0.005;
5094  fC[17][1] = 0.005;
5095  fC[17][2] = 1.0000;
5096  fC[17][3] = 0.3788;
5097  fC[17][4] = 0.3011;
5098  }
5099 
5100  for(Int_t i=18;i<fgkPIDptBin;i++){
5101  fBinLimitPID[i] = 3.0 + 0.2 * (i-18);
5102  fC[i][0] = fC[17][0];
5103  fC[i][1] = fC[17][1];
5104  fC[i][2] = fC[17][2];
5105  fC[i][3] = fC[17][3];
5106  fC[i][4] = fC[17][4];
5107  }
5108 }
5109 
5110 //---------------------------------------------------------------//
5111 Bool_t AliFlowTrackCuts::TPCTOFagree(const AliVTrack *track)
5112 {
5113  //check pid agreement between TPC and TOF
5114  Bool_t status = kFALSE;
5115 
5116  const Float_t c = 2.99792457999999984e-02;
5117 
5118  Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
5119 
5120 
5121  Double_t exptimes[9];
5122  track->GetIntegratedTimes(exptimes);
5123 
5124  Float_t dedx = track->GetTPCsignal();
5125 
5126  Float_t p = track->P();
5127  Float_t time = track->GetTOFsignal()- fESDpid.GetTOFResponse().GetStartTime(p);
5128  Float_t tl = exptimes[0]*c; // to work both for ESD and AOD
5129 
5130  Float_t betagammares = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
5131 
5132  Float_t betagamma1 = tl/(time-5 *betagammares) * 33.3564095198152043;
5133 
5134 // printf("betagamma1 = %f\n",betagamma1);
5135 
5136  if(betagamma1 < 0.1) betagamma1 = 0.1;
5137 
5138  if(betagamma1 < 0.99999) betagamma1 /= TMath::Sqrt(1-betagamma1*betagamma1);
5139  else betagamma1 = 100;
5140 
5141  Float_t betagamma2 = tl/(time+5 *betagammares) * 33.3564095198152043;
5142 // printf("betagamma2 = %f\n",betagamma2);
5143 
5144  if(betagamma2 < 0.1) betagamma2 = 0.1;
5145 
5146  if(betagamma2 < 0.99999) betagamma2 /= TMath::Sqrt(1-betagamma2*betagamma2);
5147  else betagamma2 = 100;
5148 
5149 
5150  Float_t momtpc=track->GetTPCmomentum();
5151 
5152  for(Int_t i=0;i < 5;i++){
5153  Float_t resolutionTOF = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[i], mass[i]);
5154  if(TMath::Abs(exptimes[i] - time) < 5 * resolutionTOF){
5155  Float_t dedxExp = 0;
5156  if(i==0) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kElectron);
5157  else if(i==1) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kMuon);
5158  else if(i==2) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kPion);
5159  else if(i==3) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kKaon);
5160  else if(i==4) dedxExp = fESDpid.GetTPCResponse().GetExpectedSignal(momtpc,AliPID::kProton);
5161 
5162  Float_t resolutionTPC = 2;
5163  if(i==0) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kElectron);
5164  else if(i==1) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kMuon);
5165  else if(i==2) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kPion);
5166  else if(i==3) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kKaon);
5167  else if(i==4) resolutionTPC = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
5168 
5169  if(TMath::Abs(dedx - dedxExp) < 3 * resolutionTPC){
5170  status = kTRUE;
5171  }
5172  }
5173  }
5174 
5175  Float_t bb1 = fESDpid.GetTPCResponse().Bethe(betagamma1);
5176  Float_t bb2 = fESDpid.GetTPCResponse().Bethe(betagamma2);
5177  Float_t bbM = fESDpid.GetTPCResponse().Bethe((betagamma1+betagamma2)*0.5);
5178 
5179 
5180  // status = kFALSE;
5181  // for nuclei
5182  Float_t resolutionTOFpr = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[4], mass[4]);
5183  Float_t resolutionTPCpr = fESDpid.GetTPCResponse().GetExpectedSigma(momtpc,track->GetTPCsignalN(),AliPID::kProton);
5184  if(TMath::Abs(dedx-bb1) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
5185  status = kTRUE;
5186  }
5187  else if(TMath::Abs(dedx-bb2) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
5188  status = kTRUE;
5189  }
5190  else if(TMath::Abs(dedx-bbM) < resolutionTPCpr*3 && exptimes[4] < time-7*resolutionTOFpr){
5191  status = kTRUE;
5192  }
5193 
5194  return status;
5195 }
5196 // end part added by F. Noferini
5197 //-----------------------------------------------------------------------
5198 Bool_t AliFlowTrackCuts::PassesTPCTPCTOFNsigmaCut(const AliAODTrack* track){
5199 
5200  //Added by Bernhard Hohlweger (bernhard.hohlweger@cern.ch)
5201  //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.
5202  //Default values tuned on 2010h.
5203  //Pt threshold can be set by calling AliFlowTrackCuts::SetPtTOFPIDoff(Double_t pt).
5204 
5205  Bool_t pass = kTRUE;
5206 
5207  if(!fPIDResponse) pass = kFALSE;
5208  if(!track) pass = kFALSE;
5209 
5210  // check TPC status
5211  if(track->GetTPCsignal() < 10) pass = kFALSE;
5212 
5213  if(fPtTOFPIDoff == 0.){
5214  if(fParticleID == AliPID::kProton){
5215  fPtTOFPIDoff = 0.75;
5216  }else if(fParticleID == AliPID::kPion){
5217  fPtTOFPIDoff = 0.55;
5218  }else if(fParticleID == AliPID::kKaon){
5219  fPtTOFPIDoff = 0.45;
5220  }else{
5221  //by default just use the standard TPCTOFNsigmaCut
5222  fPtTOFPIDoff = 0.;
5223  }
5224  }
5225  if(pass){
5226  Double_t Pt = track->Pt();
5227  Float_t nsigmaTPC = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC,track,fParticleID);
5228  Float_t nsigma2 = 999.;
5229  if(Pt < fPtTOFPIDoff){
5230  nsigma2 = nsigmaTPC*nsigmaTPC;
5231  }else{
5232  // check TOF status
5233  if (((track->GetStatus()&AliVTrack::kTOFout)==0)&&((track->GetStatus()&AliVTrack::kTIME)==0)){
5234  pass = kFALSE;
5235  }else{
5236  Float_t nsigmaTOF = fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF,track,fParticleID);
5237  nsigma2 = nsigmaTPC*nsigmaTPC + nsigmaTOF*nsigmaTOF;
5238  }
5239  }
5240  pass = pass && (nsigma2 < fNsigmaCut2);
5241  }
5242 
5243  return pass;
5244 
5245 }
5246 //-----------------------------------------------------------------------
5247 const char* AliFlowTrackCuts::PIDsourceName(PIDsource s)
5248 {
5249  //get the name of the particle id source
5250  switch (s)
5251  {
5252  case kTPCdedx:
5253  return "TPCdedx";
5254  case kTPCbayesian:
5255  return "TPCbayesian";
5256  case kTOFbeta:
5257  return "TOFbeta";
5258  case kTPCpid:
5259  return "TPCpid";
5260  case kTOFpid:
5261  return "TOFpid";
5262  case kTOFbayesian:
5263  return "TOFbayesianPID";
5264  case kTOFbetaSimple:
5265  return "TOFbetaSimple";
5266  case kTPCNuclei:
5267  return "TPCnuclei";
5268  case kTPCTOFNsigma:
5269  return "TPCTOFNsigma";
5270  case kTPCTOFNsigmaPurity:
5271  return "TPCTOFNsigmaPurity";
5272  default:
5273  return "NOPID";
5274  }
5275 }
5276 
5277 //-----------------------------------------------------------------------
5278 const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
5279 {
5280  //return the name of the selected parameter type
5281  switch (type)
5282  {
5283  case kMC:
5284  return "MC";
5285  case kGlobal:
5286  return "Global";
5287  case kTPCstandalone:
5288  return "TPCstandalone";
5289  case kSPDtracklet:
5290  return "SPDtracklets";
5291  case kPMD:
5292  return "PMD";
5293  case kVZERO:
5294  return "VZERO";
5295  case kBetaVZERO:
5296  return "BetaVZERO";
5297  case kDeltaVZERO:
5298  return "DeltaVZERO";
5299  case kKappaVZERO:
5300  return "KappaVZERO";
5301  case kHotfixHI:
5302  return "HotfixHI";
5303  case kMUON: // XZhang 20120604
5304  return "MUON"; // XZhang 20120604
5305  case kKink:
5306  return "Kink";
5307  case kV0:
5308  return "V0";
5309  default:
5310  return "unknown";
5311  }
5312 }
5313 
5314 //-----------------------------------------------------------------------
5315 Bool_t AliFlowTrackCuts::PassesPMDcuts(const AliESDPmdTrack* track )
5316 {
5317  //check PMD specific cuts
5318  //clean up from last iteration, and init label
5319  Int_t det = track->GetDetector();
5320  //Int_t smn = track->GetSmn();
5321  Float_t clsX = track->GetClusterX();
5322  Float_t clsY = track->GetClusterY();
5323  Float_t clsZ = track->GetClusterZ();
5324  Float_t ncell = track->GetClusterCells();
5325  Float_t adc = track->GetClusterADC();
5326 
5327  fTrack = NULL;
5328  fMCparticle=NULL;
5329  fTrackLabel=-996;
5330 
5331  fTrackEta = GetPmdEta(clsX,clsY,clsZ);
5332  fTrackPhi = GetPmdPhi(clsX,clsY);
5333  fTrackWeight = 1.0;
5334 
5335  Bool_t pass=kTRUE;
5336  if (fCutEta) {if ( fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
5337  if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
5338  if (fCutPmdDet) {if(det != fPmdDet) pass = kFALSE;}
5339  if (fCutPmdAdc) {if(adc < fPmdAdc) pass = kFALSE;}
5340  if (fCutPmdNcell) {if(ncell < fPmdNcell) pass = kFALSE;}
5341 
5342  return pass;
5343 }
5344 
5345 //-----------------------------------------------------------------------
5346 Bool_t AliFlowTrackCuts::PassesVZEROcuts(Int_t id)
5347 {
5348  //check VZERO cuts
5349  if (id<0) return kFALSE;
5350 
5351  //clean up from last iter
5352  ClearTrack();
5353 
5354  fTrackPhi = TMath::PiOver4()*(0.5+id%8);
5355 
5356  // 10102013 weighting vzero tiles - rbertens@cern.ch
5357  if(!fVZEROgainEqualization) {
5358  // if for some reason the equalization is not initialized (e.g. 2011 data)
5359  // the fVZEROxpol[] weights are used to enable or disable vzero rings
5360  if(id<32) { // v0c side
5361  fTrackEta = -3.45+0.5*(id/8);
5362  if(id < 8) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[0];
5363  else if (id < 16 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[1];
5364  else if (id < 24 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[2];
5365  else if (id < 32 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[3];
5366  } else { // v0a side
5367  fTrackEta = +4.8-0.6*((id/8)-4);
5368  if( id < 40) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[0];
5369  else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[1];
5370  else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[2];
5371  else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[3];
5372  }
5373  } else { // the equalization is initialized
5374  // note that disabled rings have already been excluded on calibration level in
5375  // AliFlowEvent (so for a disabled ring, fVZEROxpol is zero
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]/fVZEROgainEqualization->GetBinContent(1+id);
5379  else if (id < 16 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[1]/fVZEROgainEqualization->GetBinContent(1+id);
5380  else if (id < 24 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[2]/fVZEROgainEqualization->GetBinContent(1+id);
5381  else if (id < 32 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[3]/fVZEROgainEqualization->GetBinContent(1+id);
5382  } else { // v0a side
5383  fTrackEta = +4.8-0.6*((id/8)-4);
5384  if( id < 40) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[0]/fVZEROgainEqualization->GetBinContent(1+id);
5385  else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[1]/fVZEROgainEqualization->GetBinContent(1+id);
5386  else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[2]/fVZEROgainEqualization->GetBinContent(1+id);
5387  else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[3]/fVZEROgainEqualization->GetBinContent(1+id);
5388  }
5389  // printf ( " tile %i and weight %.2f \n", id, fTrackWeight);
5390  }
5391 
5392  // 29052015 weighting vzero tiles - jacopo.margutti@cern.ch
5393  if(fVZEROgainEqualizationCen) {
5394  //Added by Bernhard Hohlweger - bhohlweg@cern.ch
5395  Double_t EventCentrality = -999;
5396  if(fEvent->GetRunNumber() < 209122){
5397  //For Run1 Data the Old Centrality Percentile Method is available whereas for Run2 a new method was implemented
5398  //Cut was done for the first run of the LHC15a period
5399  EventCentrality = fEvent->GetCentrality()->GetCentralityPercentile("V0M");
5400  //09062016 Bernhard Hohlweger
5401  Double_t CorrectionFactor = fVZEROgainEqualizationCen->GetBinContent(fVZEROgainEqualizationCen->FindBin(id,EventCentrality));
5402  // the fVZEROxpol[] weights are used to enable or disable vzero rings
5403  if(id<32) { // v0c side
5404  fTrackEta = -3.45+0.5*(id/8);
5405  if(id < 8) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[0];
5406  else if (id < 16 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[1];
5407  else if (id < 24 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[2];
5408  else if (id < 32 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROCpol[3];
5409  } else { // v0a side
5410  fTrackEta = +4.8-0.6*((id/8)-4);
5411  if( id < 40) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[0];
5412  else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[1];
5413  else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[2];
5414  else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROEqMultiplicity(id)*fVZEROApol[3];
5415  }
5416  if(CorrectionFactor) {
5417  fTrackWeight /= CorrectionFactor;
5418  }
5419  }else{//for Run2 Data
5420  Double_t CorrectionFactor = -1.;
5421  AliMultSelection *MultSelection = 0x0;
5422  MultSelection = (AliMultSelection * ) fEvent->FindListObject("MultSelection");
5423  if( !MultSelection) {
5424  //If you get this warning (and EventCentrality -999) please check that the AliMultSelectionTask actually ran (before your task)
5425  AliFatal("AliMultSelection not found, did you Run AliMultSelectionTask? \n");
5426  }else{
5427  EventCentrality = MultSelection->GetMultiplicityPercentile("V0M");
5428  if(EventCentrality < 0 && EventCentrality > 100){
5429  AliWarning("No Correction Available for this Centrality \n");
5430  }else{
5431  CorrectionFactor = fVZEROgainEqualizationCen->GetBinContent(fVZEROgainEqualizationCen->FindBin(id,EventCentrality));
5432  }
5433  }
5434  if(id<32) { // v0c side
5435  fTrackEta = -3.45+0.5*(id/8);
5436  if(id < 8) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[0];
5437  else if (id < 16 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[1];
5438  else if (id < 24 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[2];
5439  else if (id < 32 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROCpol[3];
5440  } else { // v0a side
5441  fTrackEta = +4.8-0.6*((id/8)-4);
5442  if( id < 40) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[0];
5443  else if ( id < 48 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[1];
5444  else if ( id < 56 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[2];
5445  else if ( id < 64 ) fTrackWeight = fEvent->GetVZEROData()->GetMultiplicity(id)*fVZEROApol[3];
5446  }
5447  if(CorrectionFactor>0) {
5448  fTrackWeight /= CorrectionFactor;
5449  }
5450  }
5451  }//end of if(fVZEROgainEqualizationCen)
5452 
5453  if (fLinearizeVZEROresponse && id < 64)
5454  {
5455  //this is only needed in pass1 of LHC10h
5456  Float_t multVZERO[fgkNumberOfVZEROtracks];
5457  Float_t dummy=0.0;
5458  AliESDUtils::GetCorrV0((AliESDEvent*)fEvent,dummy,multVZERO);
5459  fTrackWeight = multVZERO[id];
5460  }
5461 
5462  Bool_t pass=kTRUE;
5463  if (fCutEta) {if ( fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) pass = kFALSE;}
5464  if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) pass = kFALSE;}
5465 
5466  return pass;
5467 }
5468 
5469 //-----------------------------------------------------------------------------
5470 Bool_t AliFlowTrackCuts::PassesMuonCuts(AliVParticle* vparticle)
5471 {
5472 // XZhang 20120604
5473  ClearTrack();
5474  fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
5475  if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
5476  else fMCparticle=NULL;
5477 
5478  AliESDMuonTrack *esdTrack = dynamic_cast<AliESDMuonTrack*>(vparticle);
5479  AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(vparticle);
5480  if ((!esdTrack) && (!aodTrack)) return kFALSE;
5481  if (!fMuonTrackCuts->IsSelected(vparticle)) return kFALSE;
5482  HandleVParticle(vparticle); if (!fTrack) return kFALSE;
5483  return kTRUE;
5484 }
5485 
5486 //----------------------------------------------------------------------------//
5487 Double_t AliFlowTrackCuts::GetPmdEta(Float_t xPos, Float_t yPos, Float_t zPos)
5488 {
5489  //get PMD "track" eta coordinate
5490  Float_t rpxpy, theta, eta;
5491  rpxpy = TMath::Sqrt(xPos*xPos + yPos*yPos);
5492  theta = TMath::ATan2(rpxpy,zPos);
5493  eta = -TMath::Log(TMath::Tan(0.5*theta));
5494  return eta;
5495 }
5496 
5497 //--------------------------------------------------------------------------//
5498 Double_t AliFlowTrackCuts::GetPmdPhi(Float_t xPos, Float_t yPos)
5499 {
5500  //Get PMD "track" phi coordinate
5501  Float_t pybypx, phi = 0., phi1;
5502  if(xPos==0)
5503  {
5504  if(yPos>0) phi = 90.;
5505  if(yPos<0) phi = 270.;
5506  }
5507  if(xPos != 0)
5508  {
5509  pybypx = yPos/xPos;
5510  if(pybypx < 0) pybypx = - pybypx;
5511  phi1 = TMath::ATan(pybypx)*180./3.14159;
5512 
5513  if(xPos > 0 && yPos > 0) phi = phi1; // 1st Quadrant
5514  if(xPos < 0 && yPos > 0) phi = 180 - phi1; // 2nd Quadrant
5515  if(xPos < 0 && yPos < 0) phi = 180 + phi1; // 3rd Quadrant
5516  if(xPos > 0 && yPos < 0) phi = 360 - phi1; // 4th Quadrant
5517 
5518  }
5519  phi = phi*3.14159/180.;
5520  return phi;
5521 }
5522 
5523 //---------------------------------------------------------------//
5524 void AliFlowTrackCuts::Browse(TBrowser* b)
5525 {
5526  //some browsing capabilities
5527  if (fQA) b->Add(fQA);
5528 }
5529 
5530 //---------------------------------------------------------------//
5531 Long64_t AliFlowTrackCuts::Merge(TCollection* list)
5532 {
5533  //merge
5534  if (!fQA || !list) return 0;
5535  if (list->IsEmpty()) return 0;
5536  AliFlowTrackCuts* obj=NULL;
5537  TList tmplist;
5538  TIter next(list);
5539  while ( (obj = dynamic_cast<AliFlowTrackCuts*>(next())) )
5540  {
5541  if (obj==this) continue;
5542  tmplist.Add(obj->GetQA());
5543  }
5544  return fQA->Merge(&tmplist);
5545 }
5546 //________________________________________________________________//
5547 void AliFlowTrackCuts::SetTPCTOFNsigmaPIDPurityFunctions(Float_t PurityLevel)
5548 {
5549  fPurityLevel = PurityLevel;
5550  cout<<"fPurityLevel = "<<fPurityLevel<<endl;
5551 
5552  fPurityFunctionsFile = TFile::Open(Form("alien:///alice/cern.ch/user/n/nmohamma/PurityFunctions_%i-%icent.root",fCentralityPercentileMin,fCentralityPercentileMax));
5553  if((!fPurityFunctionsFile) || (!fPurityFunctionsFile->IsOpen())) {
5554  printf("The purity functions file does not exist");
5555  return;
5556  }
5557  fPurityFunctionsList=(TDirectory*)fPurityFunctionsFile->Get("Filterbit1");
5558  if(!fPurityFunctionsList){printf("The purity functions list is empty"); return;}
5559 
5560  TString species[3] = {"pion","kaon","proton"};