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