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