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