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