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