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