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