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