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