double for 4-momentum in Generate()
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliFlowTrackCuts.cxx
CommitLineData
daf66719 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 **************************************************************************/
15
16/* $Id$ */
17
18// AliFlowTrackCuts:
19// ESD track cuts for flow framework
20//
21// origin: Mikolaj Krzewicki (mikolaj.krzewicki@cern.ch)
5559ce24 22//
23// This class gurantees consistency of cut methods, trackparameter
24// selection (global tracks, TPC only, etc..) and parameter mixing
25// in the flow framework. Transparently handles different input types:
26// ESD, MC, AOD.
27// This class works in 2 steps: first the requested track parameters are
28// constructed (to be set by SetParamType() ), then cuts are applied.
29// the constructed track can be requested AFTER checking the cuts by
30// calling GetTrack(), in this case the cut object stays in control,
31// caller does not have to delete the track.
32// Additionally caller can request an AliFlowTrack object to be constructed
33// according the parameter mixing scenario requested by SetParamMix().
34// AliFlowTrack is made using MakeFlowTrack() method, its an 'object factory'
35// so caller needs to take care of the freshly created object.
daf66719 36
37#include <limits.h>
38#include <float.h>
32b846cd 39#include <TMatrix.h>
9a0783cc 40#include "TParticle.h"
2948ac5a 41#include "TObjArray.h"
9a0783cc 42#include "AliStack.h"
daf66719 43#include "AliMCEvent.h"
9a0783cc 44#include "AliESDEvent.h"
daf66719 45#include "AliVParticle.h"
46#include "AliMCParticle.h"
47#include "AliESDtrack.h"
12b2b8bc 48#include "AliMultiplicity.h"
daf66719 49#include "AliAODTrack.h"
50#include "AliFlowTrack.h"
51#include "AliFlowTrackCuts.h"
52#include "AliLog.h"
32b846cd 53#include "AliESDpid.h"
daf66719 54
55ClassImp(AliFlowTrackCuts)
56
57//-----------------------------------------------------------------------
58AliFlowTrackCuts::AliFlowTrackCuts():
59 AliFlowTrackSimpleCuts(),
441ea1cf 60 fAliESDtrackCuts(NULL),
61 fQA(NULL),
62 fCutMC(kFALSE),
63 fCutMCprocessType(kFALSE),
64 fMCprocessType(kPNoProcess),
65 fCutMCPID(kFALSE),
66 fMCPID(0),
67 fIgnoreSignInPID(kFALSE),
68 fCutMCisPrimary(kFALSE),
69 fRequireTransportBitForPrimaries(kTRUE),
70 fMCisPrimary(kFALSE),
71 fRequireCharge(kFALSE),
72 fFakesAreOK(kTRUE),
73 fCutSPDtrackletDeltaPhi(kFALSE),
74 fSPDtrackletDeltaPhiMax(FLT_MAX),
75 fSPDtrackletDeltaPhiMin(-FLT_MAX),
76 fIgnoreTPCzRange(kFALSE),
77 fIgnoreTPCzRangeMax(FLT_MAX),
78 fIgnoreTPCzRangeMin(-FLT_MAX),
79 fCutChi2PerClusterTPC(kFALSE),
80 fMaxChi2PerClusterTPC(FLT_MAX),
81 fMinChi2PerClusterTPC(-FLT_MAX),
82 fCutNClustersTPC(kFALSE),
83 fNClustersTPCMax(INT_MAX),
84 fNClustersTPCMin(INT_MIN),
42655d16 85 fCutNClustersITS(kFALSE),
86 fNClustersITSMax(INT_MAX),
87 fNClustersITSMin(INT_MIN),
5ba02b32 88 fUseAODFilterBit(kFALSE),
89 fAODFilterBit(0),
a63303bd 90 fCutDCAToVertexXY(kFALSE),
91 fCutDCAToVertexZ(kFALSE),
441ea1cf 92 fParamType(kGlobal),
93 fParamMix(kPure),
94 fTrack(NULL),
95 fTrackPhi(0.),
96 fTrackEta(0.),
97 fTrackWeight(0.),
98 fTrackLabel(INT_MIN),
99 fMCevent(NULL),
100 fMCparticle(NULL),
101 fEvent(NULL),
102 fTPCtrack(),
aab6527a 103 fESDpid(),
441ea1cf 104 fPIDsource(kTPCTOFpid),
105 fTPCpidCuts(NULL),
106 fTOFpidCuts(NULL),
107 fTPCTOFpidCrossOverPt(0.4),
108 fAliPID(AliPID::kPion)
109{
110 //io constructor
5104aa35 111 for ( Int_t i=0; i<5; i++ ) { fProbBayes[i]=0.0; }
422c61c7 112 SetPriors(); //init arrays
441ea1cf 113}
114
115//-----------------------------------------------------------------------
116AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
117 AliFlowTrackSimpleCuts(),
daf66719 118 fAliESDtrackCuts(new AliESDtrackCuts()),
a1c43d26 119 fQA(NULL),
1ff4bda1 120 fCutMC(kFALSE),
daf66719 121 fCutMCprocessType(kFALSE),
122 fMCprocessType(kPNoProcess),
123 fCutMCPID(kFALSE),
124 fMCPID(0),
4cbcbead 125 fIgnoreSignInPID(kFALSE),
daf66719 126 fCutMCisPrimary(kFALSE),
441ea1cf 127 fRequireTransportBitForPrimaries(kTRUE),
daf66719 128 fMCisPrimary(kFALSE),
957517fa 129 fRequireCharge(kFALSE),
127a5825 130 fFakesAreOK(kTRUE),
9a0783cc 131 fCutSPDtrackletDeltaPhi(kFALSE),
132 fSPDtrackletDeltaPhiMax(FLT_MAX),
133 fSPDtrackletDeltaPhiMin(-FLT_MAX),
1ff4bda1 134 fIgnoreTPCzRange(kFALSE),
135 fIgnoreTPCzRangeMax(FLT_MAX),
136 fIgnoreTPCzRangeMin(-FLT_MAX),
2948ac5a 137 fCutChi2PerClusterTPC(kFALSE),
138 fMaxChi2PerClusterTPC(FLT_MAX),
139 fMinChi2PerClusterTPC(-FLT_MAX),
32b846cd 140 fCutNClustersTPC(kFALSE),
141 fNClustersTPCMax(INT_MAX),
142 fNClustersTPCMin(INT_MIN),
42655d16 143 fCutNClustersITS(kFALSE),
144 fNClustersITSMax(INT_MAX),
5ba02b32 145 fNClustersITSMin(INT_MIN),
146 fUseAODFilterBit(kFALSE),
147 fAODFilterBit(0),
a63303bd 148 fCutDCAToVertexXY(kFALSE),
149 fCutDCAToVertexZ(kFALSE),
12b2b8bc 150 fParamType(kGlobal),
daf66719 151 fParamMix(kPure),
daf66719 152 fTrack(NULL),
12b2b8bc 153 fTrackPhi(0.),
154 fTrackEta(0.),
155 fTrackWeight(0.),
127a5825 156 fTrackLabel(INT_MIN),
957517fa 157 fMCevent(NULL),
9a0783cc 158 fMCparticle(NULL),
1ff4bda1 159 fEvent(NULL),
32b846cd 160 fTPCtrack(),
aab6527a 161 fESDpid(),
32b846cd 162 fPIDsource(kTPCTOFpid),
163 fTPCpidCuts(NULL),
164 fTOFpidCuts(NULL),
165 fTPCTOFpidCrossOverPt(0.4),
166 fAliPID(AliPID::kPion)
daf66719 167{
168 //constructor
441ea1cf 169 SetName(name);
170 SetTitle("AliFlowTrackCuts");
aab6527a 171 fESDpid.GetTPCResponse().SetBetheBlochParameters( 0.0283086,
172 2.63394e+01,
173 5.04114e-11,
174 2.12543e+00,
175 4.88663e+00 );
5104aa35 176 for ( Int_t i=0; i<5; i++ ) { fProbBayes[i]=0.0; }
422c61c7 177 SetPriors(); //init arrays
daf66719 178}
179
180//-----------------------------------------------------------------------
ee242db3 181AliFlowTrackCuts::AliFlowTrackCuts(const AliFlowTrackCuts& that):
182 AliFlowTrackSimpleCuts(that),
441ea1cf 183 fAliESDtrackCuts(NULL),
a1c43d26 184 fQA(NULL),
1ff4bda1 185 fCutMC(that.fCutMC),
ee242db3 186 fCutMCprocessType(that.fCutMCprocessType),
187 fMCprocessType(that.fMCprocessType),
188 fCutMCPID(that.fCutMCPID),
189 fMCPID(that.fMCPID),
4cbcbead 190 fIgnoreSignInPID(that.fIgnoreSignInPID),
ee242db3 191 fCutMCisPrimary(that.fCutMCisPrimary),
441ea1cf 192 fRequireTransportBitForPrimaries(that.fRequireTransportBitForPrimaries),
ee242db3 193 fMCisPrimary(that.fMCisPrimary),
194 fRequireCharge(that.fRequireCharge),
195 fFakesAreOK(that.fFakesAreOK),
196 fCutSPDtrackletDeltaPhi(that.fCutSPDtrackletDeltaPhi),
197 fSPDtrackletDeltaPhiMax(that.fSPDtrackletDeltaPhiMax),
198 fSPDtrackletDeltaPhiMin(that.fSPDtrackletDeltaPhiMin),
1ff4bda1 199 fIgnoreTPCzRange(that.fIgnoreTPCzRange),
200 fIgnoreTPCzRangeMax(that.fIgnoreTPCzRangeMax),
201 fIgnoreTPCzRangeMin(that.fIgnoreTPCzRangeMin),
2948ac5a 202 fCutChi2PerClusterTPC(that.fCutChi2PerClusterTPC),
203 fMaxChi2PerClusterTPC(that.fMaxChi2PerClusterTPC),
204 fMinChi2PerClusterTPC(that.fMinChi2PerClusterTPC),
32b846cd 205 fCutNClustersTPC(that.fCutNClustersTPC),
206 fNClustersTPCMax(that.fNClustersTPCMax),
207 fNClustersTPCMin(that.fNClustersTPCMin),
42655d16 208 fCutNClustersITS(that.fCutNClustersITS),
209 fNClustersITSMax(that.fNClustersITSMax),
210 fNClustersITSMin(that.fNClustersITSMin),
5ba02b32 211 fUseAODFilterBit(that.fUseAODFilterBit),
212 fAODFilterBit(that.fAODFilterBit),
a63303bd 213 fCutDCAToVertexXY(that.fCutDCAToVertexXY),
214 fCutDCAToVertexZ(that.fCutDCAToVertexZ),
ee242db3 215 fParamType(that.fParamType),
216 fParamMix(that.fParamMix),
daf66719 217 fTrack(NULL),
ee242db3 218 fTrackPhi(0.),
219 fTrackEta(0.),
220 fTrackWeight(0.),
127a5825 221 fTrackLabel(INT_MIN),
957517fa 222 fMCevent(NULL),
9a0783cc 223 fMCparticle(NULL),
1ff4bda1 224 fEvent(NULL),
32b846cd 225 fTPCtrack(),
226 fESDpid(that.fESDpid),
227 fPIDsource(that.fPIDsource),
228 fTPCpidCuts(NULL),
229 fTOFpidCuts(NULL),
230 fTPCTOFpidCrossOverPt(that.fTPCTOFpidCrossOverPt),
231 fAliPID(that.fAliPID)
daf66719 232{
233 //copy constructor
32b846cd 234 if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
235 if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
441ea1cf 236 if (that.fAliESDtrackCuts) fAliESDtrackCuts = new AliESDtrackCuts(*(that.fAliESDtrackCuts));
5104aa35 237 memcpy(fProbBayes,that.fProbBayes,sizeof(fProbBayes));
422c61c7 238 SetPriors(); //init arrays
daf66719 239}
240
241//-----------------------------------------------------------------------
ee242db3 242AliFlowTrackCuts& AliFlowTrackCuts::operator=(const AliFlowTrackCuts& that)
daf66719 243{
244 //assignment
ee242db3 245 AliFlowTrackSimpleCuts::operator=(that);
441ea1cf 246 if (that.fAliESDtrackCuts) *fAliESDtrackCuts=*(that.fAliESDtrackCuts);
a1c43d26 247 fQA=NULL;
1ff4bda1 248 fCutMC=that.fCutMC;
ee242db3 249 fCutMCprocessType=that.fCutMCprocessType;
250 fMCprocessType=that.fMCprocessType;
251 fCutMCPID=that.fCutMCPID;
252 fMCPID=that.fMCPID;
4cbcbead 253 fIgnoreSignInPID=that.fIgnoreSignInPID,
ee242db3 254 fCutMCisPrimary=that.fCutMCisPrimary;
441ea1cf 255 fRequireTransportBitForPrimaries=that.fRequireTransportBitForPrimaries;
ee242db3 256 fMCisPrimary=that.fMCisPrimary;
257 fRequireCharge=that.fRequireCharge;
258 fFakesAreOK=that.fFakesAreOK;
259 fCutSPDtrackletDeltaPhi=that.fCutSPDtrackletDeltaPhi;
260 fSPDtrackletDeltaPhiMax=that.fSPDtrackletDeltaPhiMax;
261 fSPDtrackletDeltaPhiMin=that.fSPDtrackletDeltaPhiMin;
1ff4bda1 262 fIgnoreTPCzRange=that.fIgnoreTPCzRange;
263 fIgnoreTPCzRangeMax=that.fIgnoreTPCzRangeMax;
264 fIgnoreTPCzRangeMin=that.fIgnoreTPCzRangeMin;
2948ac5a 265 fCutChi2PerClusterTPC=that.fCutChi2PerClusterTPC;
266 fMaxChi2PerClusterTPC=that.fMaxChi2PerClusterTPC;
267 fMinChi2PerClusterTPC=that.fMinChi2PerClusterTPC;
32b846cd 268 fCutNClustersTPC=that.fCutNClustersTPC;
269 fNClustersTPCMax=that.fNClustersTPCMax;
270 fNClustersTPCMin=that.fNClustersTPCMin;
42655d16 271 fCutNClustersITS=that.fCutNClustersITS;
272 fNClustersITSMax=that.fNClustersITSMax;
273 fNClustersITSMin=that.fNClustersITSMin;
a63303bd 274 fUseAODFilterBit=that.fUseAODFilterBit;
275 fAODFilterBit=that.fAODFilterBit;
276 fCutDCAToVertexXY=that.fCutDCAToVertexXY;
277 fCutDCAToVertexZ=that.fCutDCAToVertexZ;
ee242db3 278 fParamType=that.fParamType;
279 fParamMix=that.fParamMix;
daf66719 280
daf66719 281 fTrack=NULL;
ee242db3 282 fTrackPhi=0.;
283 fTrackPhi=0.;
284 fTrackWeight=0.;
127a5825 285 fTrackLabel=INT_MIN;
957517fa 286 fMCevent=NULL;
daf66719 287 fMCparticle=NULL;
9a0783cc 288 fEvent=NULL;
daf66719 289
32b846cd 290 fESDpid = that.fESDpid;
291 fPIDsource = that.fPIDsource;
292
293 if (that.fTPCpidCuts) fTPCpidCuts = new TMatrixF(*(that.fTPCpidCuts));
294 if (that.fTOFpidCuts) fTOFpidCuts = new TMatrixF(*(that.fTOFpidCuts));
295 fTPCTOFpidCrossOverPt=that.fTPCTOFpidCrossOverPt;
296
297 fAliPID=that.fAliPID;
5104aa35 298 memcpy(fProbBayes,that.fProbBayes,sizeof(fProbBayes));
32b846cd 299
daf66719 300 return *this;
301}
302
303//-----------------------------------------------------------------------
304AliFlowTrackCuts::~AliFlowTrackCuts()
305{
306 //dtor
daf66719 307 delete fAliESDtrackCuts;
32b846cd 308 delete fTPCpidCuts;
309 delete fTOFpidCuts;
daf66719 310}
311
312//-----------------------------------------------------------------------
aab6527a 313void AliFlowTrackCuts::SetEvent(AliVEvent* event, AliMCEvent* mcEvent)
314{
315 //set the event
316 Clear();
317 fEvent=event;
318 fMCevent=mcEvent;
319
320 //do the magic for ESD
321 AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(event);
322 if (fCutPID && myESD)
323 {
324 //TODO: maybe call it only for the TOF options?
325 // Added by F. Noferini for TOF PID
326 fESDpid.SetTOFResponse(myESD,AliESDpid::kTOF_T0);
327 fESDpid.MakePID(myESD,kFALSE);
328 // End F. Noferini added part
329 }
330
331 //TODO: AOD
332}
333
334//-----------------------------------------------------------------------
335void AliFlowTrackCuts::SetCutMC( Bool_t b )
336{
337 //will we be cutting on MC information?
338 fCutMC=b;
339
340 //if we cut on MC info then also the Bethe Bloch should be the one tuned for MC
341 if (fCutMC)
342 {
343 fESDpid.GetTPCResponse().SetBetheBlochParameters( 2.15898e+00/50.,
344 1.75295e+01,
345 3.40030e-09,
346 1.96178e+00,
347 3.91720e+00);
348 }
349}
350
351//-----------------------------------------------------------------------
12b2b8bc 352Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id)
daf66719 353{
354 //check cuts
355 AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
356 if (vparticle) return PassesCuts(vparticle);
357 AliFlowTrackSimple* flowtrack = dynamic_cast<AliFlowTrackSimple*>(obj);
358 if (flowtrack) return PassesCuts(flowtrack);
12b2b8bc 359 AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
360 if (tracklets) return PassesCuts(tracklets,id);
daf66719 361 return kFALSE; //default when passed wrong type of object
362}
363
364//-----------------------------------------------------------------------
1ff4bda1 365Bool_t AliFlowTrackCuts::IsSelectedMCtruth(TObject* obj, Int_t id)
366{
367 //check cuts
368 AliVParticle* vparticle = dynamic_cast<AliVParticle*>(obj);
369 if (vparticle)
370 {
371 return PassesMCcuts(fMCevent,vparticle->GetLabel());
372 }
373 AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);
374 if (tracklets)
375 {
376 Int_t label0 = tracklets->GetLabel(id,0);
377 Int_t label1 = tracklets->GetLabel(id,1);
378 Int_t label = (label0==label1)?tracklets->GetLabel(id,1):-666;
379 return PassesMCcuts(fMCevent,label);
380 }
381 return kFALSE; //default when passed wrong type of object
382}
383
384//-----------------------------------------------------------------------
daf66719 385Bool_t AliFlowTrackCuts::PassesCuts(AliFlowTrackSimple* track)
386{
387 //check cuts on a flowtracksimple
5559ce24 388
389 //clean up from last iteration
1ff4bda1 390 fTrack = NULL;
daf66719 391 return AliFlowTrackSimpleCuts::PassesCuts(track);
392}
393
394//-----------------------------------------------------------------------
12b2b8bc 395Bool_t AliFlowTrackCuts::PassesCuts(AliMultiplicity* tracklet, Int_t id)
396{
397 //check cuts on a tracklets
398
9a0783cc 399 //clean up from last iteration, and init label
1ff4bda1 400 fTrack = NULL;
12b2b8bc 401 fMCparticle=NULL;
9a0783cc 402 fTrackLabel=-1;
12b2b8bc 403
404 fTrackPhi = tracklet->GetPhi(id);
405 fTrackEta = tracklet->GetEta(id);
406 fTrackWeight = 1.0;
407 if (fCutEta) {if ( fTrackEta < fEtaMin || fTrackEta >= fEtaMax ) return kFALSE;}
408 if (fCutPhi) {if ( fTrackPhi < fPhiMin || fTrackPhi >= fPhiMax ) return kFALSE;}
409
410 //check MC info if available
9a0783cc 411 //if the 2 clusters have different label track cannot be good
412 //and should therefore not pass the mc cuts
413 Int_t label0 = tracklet->GetLabel(id,0);
414 Int_t label1 = tracklet->GetLabel(id,1);
d0471ea0 415 //if possible get label and mcparticle
9a0783cc 416 fTrackLabel = (label0==label1)?tracklet->GetLabel(id,1):-1;
7afa829c 417 if (!fFakesAreOK && fTrackLabel<0) return kFALSE;
d0471ea0 418 if (fTrackLabel>=0 && fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
419 //check MC cuts
1ff4bda1 420 if (fCutMC && !PassesMCcuts()) return kFALSE;
12b2b8bc 421 return kTRUE;
422}
423
424//-----------------------------------------------------------------------
1ff4bda1 425Bool_t AliFlowTrackCuts::PassesMCcuts(AliMCEvent* mcEvent, Int_t label)
12b2b8bc 426{
427 //check the MC info
1ff4bda1 428 if (!mcEvent) return kFALSE;
429 if (label<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
430 AliMCParticle* mcparticle = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
431 if (!mcparticle) {AliError("no MC track"); return kFALSE;}
12b2b8bc 432
433 if (fCutMCisPrimary)
434 {
441ea1cf 435 if (IsPhysicalPrimary(mcEvent,label,fRequireTransportBitForPrimaries) != fMCisPrimary) return kFALSE;
12b2b8bc 436 }
437 if (fCutMCPID)
438 {
1ff4bda1 439 Int_t pdgCode = mcparticle->PdgCode();
4cbcbead 440 if (fIgnoreSignInPID)
441 {
442 if (TMath::Abs(fMCPID) != TMath::Abs(pdgCode)) return kFALSE;
443 }
444 else
445 {
446 if (fMCPID != pdgCode) return kFALSE;
447 }
12b2b8bc 448 }
449 if ( fCutMCprocessType )
450 {
1ff4bda1 451 TParticle* particle = mcparticle->Particle();
12b2b8bc 452 Int_t processID = particle->GetUniqueID();
453 if (processID != fMCprocessType ) return kFALSE;
454 }
455 return kTRUE;
456}
1ff4bda1 457//-----------------------------------------------------------------------
458Bool_t AliFlowTrackCuts::PassesMCcuts()
459{
460 if (!fMCevent) return kFALSE;
461 if (fTrackLabel<0) return kFALSE;//otherwise AliCMevent prints a warning before returning NULL
462 fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
463 return PassesMCcuts(fMCevent,fTrackLabel);
464}
12b2b8bc 465
466//-----------------------------------------------------------------------
daf66719 467Bool_t AliFlowTrackCuts::PassesCuts(AliVParticle* vparticle)
468{
469 //check cuts for an ESD vparticle
470
127a5825 471 ////////////////////////////////////////////////////////////////
472 // start by preparing the track parameters to cut on //////////
473 ////////////////////////////////////////////////////////////////
5559ce24 474 //clean up from last iteration
1ff4bda1 475 fTrack=NULL;
5559ce24 476
957517fa 477 //get the label and the mc particle
127a5825 478 fTrackLabel = (fFakesAreOK)?TMath::Abs(vparticle->GetLabel()):vparticle->GetLabel();
479 if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
daf66719 480 else fMCparticle=NULL;
481
957517fa 482 Bool_t isMCparticle = kFALSE; //some things are different for MC particles, check!
daf66719 483 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*>(vparticle);
42655d16 484 AliAODTrack* aodTrack = NULL;
daf66719 485 if (esdTrack)
42655d16 486 //for an ESD track we do some magic sometimes like constructing TPC only parameters
487 //or doing some hybrid, handle that here
daf66719 488 HandleESDtrack(esdTrack);
489 else
957517fa 490 {
daf66719 491 HandleVParticle(vparticle);
957517fa 492 //now check if produced particle is MC
493 isMCparticle = (dynamic_cast<AliMCParticle*>(fTrack))!=NULL;
42655d16 494 aodTrack = dynamic_cast<AliAODTrack*>(vparticle); //keep the additional dynamic cast out of the way for ESDs
957517fa 495 }
127a5825 496 ////////////////////////////////////////////////////////////////
497 ////////////////////////////////////////////////////////////////
498
1ff4bda1 499 if (!fTrack) return kFALSE;
42655d16 500 //because it may be different from global, not needed for aodTrack because we dont do anything funky there
501 if (esdTrack) esdTrack = static_cast<AliESDtrack*>(fTrack);
2948ac5a 502
924b02b0 503 Bool_t pass=kTRUE;
9a0783cc 504 //check the common cuts for the current particle fTrack (MC,AOD,ESD)
32b846cd 505 Double_t pt = fTrack->Pt();
7afa829c 506 if (!fFakesAreOK) {if (fTrackLabel<0) pass=kFALSE;}
32b846cd 507 if (fCutPt) {if (pt < fPtMin || pt >= fPtMax ) pass=kFALSE;}
924b02b0 508 if (fCutEta) {if (fTrack->Eta() < fEtaMin || fTrack->Eta() >= fEtaMax ) pass=kFALSE;}
509 if (fCutPhi) {if (fTrack->Phi() < fPhiMin || fTrack->Phi() >= fPhiMax ) pass=kFALSE;}
510 if (fRequireCharge) {if (fTrack->Charge() == 0) pass=kFALSE;}
511 if (fCutCharge && !isMCparticle) {if (fTrack->Charge() != fCharge) pass=kFALSE;}
957517fa 512 if (fCutCharge && isMCparticle)
513 {
514 //in case of an MC particle the charge is stored in units of 1/3|e|
515 Int_t charge = TMath::Nint(fTrack->Charge()/3.0); //mc particles have charge in units of 1/3e
32b846cd 516 if (charge!=fCharge) pass=kFALSE;
957517fa 517 }
924b02b0 518 //if(fCutPID) {if (fTrack->PID() != fPID) pass=kFALSE;}
daf66719 519
957517fa 520 //when additionally MC info is required
1ff4bda1 521 if (fCutMC && !PassesMCcuts()) pass=kFALSE;
daf66719 522
42655d16 523 //the case of ESD or AOD
bb6ca457 524 if (esdTrack) { if (!PassesESDcuts(esdTrack)) { pass=kFALSE; } }
525 if (aodTrack) { if (!PassesAODcuts(aodTrack)) { pass=kFALSE; } }
42655d16 526
527 //true by default, if we didn't set any cuts
528 return pass;
529}
530
531//_______________________________________________________________________
532Bool_t AliFlowTrackCuts::PassesAODcuts(AliAODTrack* track)
533{
534 Bool_t pass = kTRUE;
535
536 if (fCutNClustersTPC)
1ff4bda1 537 {
42655d16 538 Int_t ntpccls = track->GetTPCNcls();
539 if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
540 }
541
542 if (fCutNClustersITS)
543 {
544 Int_t nitscls = track->GetITSNcls();
545 if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
546 }
5ba02b32 547
548 if (fCutChi2PerClusterTPC)
549 {
550 Double_t chi2tpc = track->Chi2perNDF();
551 if (chi2tpc < fMinChi2PerClusterTPC || chi2tpc > fMaxChi2PerClusterTPC) pass=kFALSE;
552 }
553
554 if (GetRequireTPCRefit() && !(track->GetStatus() & AliESDtrack::kTPCrefit) ) pass=kFALSE;
555 if (GetRequireITSRefit() && !(track->GetStatus() & AliESDtrack::kITSrefit) ) pass=kFALSE;
556
557 if (fUseAODFilterBit && !track->TestFilterBit(fAODFilterBit)) pass=kFALSE;
558
a63303bd 559 if (fCutDCAToVertexXY && track->DCA()>GetMaxDCAToVertexXY()) pass=kFALSE;
5ba02b32 560
42655d16 561
562 return pass;
563}
564
565//_______________________________________________________________________
566Bool_t AliFlowTrackCuts::PassesESDcuts(AliESDtrack* track)
567{
568 Bool_t pass=kTRUE;
569 if (fIgnoreTPCzRange)
570 {
571 const AliExternalTrackParam* pin = track->GetOuterParam();
572 const AliExternalTrackParam* pout = track->GetInnerParam();
573 if (pin&&pout)
441ea1cf 574 {
42655d16 575 Double_t zin = pin->GetZ();
576 Double_t zout = pout->GetZ();
577 if (zin*zout<0) pass=kFALSE; //reject if cross the membrane
578 if (zin < fIgnoreTPCzRangeMin || zin > fIgnoreTPCzRangeMax) pass=kFALSE;
579 if (zout < fIgnoreTPCzRangeMin || zout > fIgnoreTPCzRangeMax) pass=kFALSE;
441ea1cf 580 }
42655d16 581 }
32b846cd 582
42655d16 583 Int_t ntpccls = ( fParamType==kESD_TPConly )?
584 track->GetTPCNclsIter1():track->GetTPCNcls();
585 if (fCutChi2PerClusterTPC)
586 {
587 Float_t tpcchi2 = (fParamType==kESD_TPConly)?
588 track->GetTPCchi2Iter1():track->GetTPCchi2();
589 tpcchi2 = (ntpccls>0)?tpcchi2/ntpccls:-FLT_MAX;
590 if (tpcchi2<fMinChi2PerClusterTPC || tpcchi2 >=fMaxChi2PerClusterTPC)
591 pass=kFALSE;
592 }
593
594 if (fCutNClustersTPC)
595 {
596 if (ntpccls < fNClustersTPCMin || ntpccls > fNClustersTPCMax) pass=kFALSE;
597 }
598
599 Int_t nitscls = track->GetNcls(0);
600 if (fCutNClustersITS)
601 {
602 if (nitscls < fNClustersITSMin || nitscls > fNClustersITSMax) pass=kFALSE;
603 }
32b846cd 604
42655d16 605 Double_t pt = fTrack->Pt();
606 if (fCutPID)
607 {
608 switch (fPIDsource)
32b846cd 609 {
42655d16 610 case kTPCpid:
611 if (!PassesTPCpidCut(track)) pass=kFALSE;
612 break;
613 case kTOFpid:
614 if (!PassesTOFpidCut(track)) pass=kFALSE;
615 break;
616 case kTPCTOFpid:
617 if (pt< fTPCTOFpidCrossOverPt)
618 {
619 if (!PassesTPCpidCut(track)) pass=kFALSE;
620 }
621 else //if (pt>=fTPCTOFpidCrossOverPt)
622 {
623 if (!PassesTOFpidCut(track)) pass=kFALSE;
624 }
625 break;
626 // part added by F. Noferini
627 case kTOFbayesian:
628 if (!PassesTOFbayesianCut(track)) pass=kFALSE;
629 break;
630 // end part added by F. Noferini
631 default:
632 printf("AliFlowTrackCuts::PassesCuts() this should never be called!\n");
633 pass=kFALSE;
634 break;
32b846cd 635 }
42655d16 636 }
32b846cd 637
42655d16 638 //some stuff is still handled by AliESDtrackCuts class - delegate
639 if (fAliESDtrackCuts)
640 {
641 if (!fAliESDtrackCuts->IsSelected(track)) pass=kFALSE;
1ff4bda1 642 }
42655d16 643
644 return pass;
daf66719 645}
646
647//-----------------------------------------------------------------------
648void AliFlowTrackCuts::HandleVParticle(AliVParticle* track)
649{
650 //handle the general case
daf66719 651 switch (fParamType)
652 {
daf66719 653 default:
daf66719 654 fTrack = track;
32b846cd 655 break;
daf66719 656 }
657}
658
659//-----------------------------------------------------------------------
660void AliFlowTrackCuts::HandleESDtrack(AliESDtrack* track)
661{
662 //handle esd track
f98ea3f1 663 AliExternalTrackParam* ip=NULL;
daf66719 664 switch (fParamType)
665 {
12b2b8bc 666 case kGlobal:
daf66719 667 fTrack = track;
daf66719 668 break;
669 case kESD_TPConly:
1ff4bda1 670 if (!track->FillTPCOnlyTrack(fTPCtrack))
671 {
672 fTrack=NULL;
673 fMCparticle=NULL;
674 fTrackLabel=-1;
675 return;
676 }
f98ea3f1 677 ip = const_cast<AliExternalTrackParam*>(fTPCtrack.GetInnerParam());
678 if (!ip) { ip = new AliExternalTrackParam(*(track->GetInnerParam())); }
679 else { *ip = *(track->GetInnerParam()); }
1ff4bda1 680 fTrack = &fTPCtrack;
957517fa 681 //recalculate the label and mc particle, they may differ as TPClabel != global label
127a5825 682 fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
683 if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
957517fa 684 else fMCparticle=NULL;
daf66719 685 break;
daf66719 686 default:
687 fTrack = track;
32b846cd 688 break;
daf66719 689 }
690}
691
692//-----------------------------------------------------------------------
693AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts()
694{
695 //get standard cuts
441ea1cf 696 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly cuts");
daf66719 697 delete cuts->fAliESDtrackCuts;
698 cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
5559ce24 699 cuts->SetParamType(kESD_TPConly);
daf66719 700 return cuts;
701}
702
703//-----------------------------------------------------------------------
704AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
705{
706 //get standard cuts
441ea1cf 707 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
daf66719 708 delete cuts->fAliESDtrackCuts;
709 cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
12b2b8bc 710 cuts->SetParamType(kGlobal);
daf66719 711 return cuts;
712}
713
714//-----------------------------------------------------------------------
41dc4195 715AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack() const
daf66719 716{
717 //get a flow track constructed from whatever we applied cuts on
718 //caller is resposible for deletion
9a0783cc 719 //if construction fails return NULL
daf66719 720 AliFlowTrack* flowtrack=NULL;
d0471ea0 721 TParticle *tmpTParticle=NULL;
722 AliMCParticle* tmpAliMCParticle=NULL;
12b2b8bc 723 if (fParamType==kESD_SPDtracklet)
daf66719 724 {
9a0783cc 725 switch (fParamMix)
726 {
727 case kPure:
d0471ea0 728 flowtrack = new AliFlowTrack();
9a0783cc 729 flowtrack->SetPhi(fTrackPhi);
730 flowtrack->SetEta(fTrackEta);
731 break;
732 case kTrackWithMCkine:
733 if (!fMCparticle) return NULL;
d0471ea0 734 flowtrack = new AliFlowTrack();
9a0783cc 735 flowtrack->SetPhi( fMCparticle->Phi() );
736 flowtrack->SetEta( fMCparticle->Eta() );
737 flowtrack->SetPt( fMCparticle->Pt() );
738 break;
739 case kTrackWithMCpt:
740 if (!fMCparticle) return NULL;
d0471ea0 741 flowtrack = new AliFlowTrack();
9a0783cc 742 flowtrack->SetPhi(fTrackPhi);
743 flowtrack->SetEta(fTrackEta);
744 flowtrack->SetPt(fMCparticle->Pt());
745 break;
d0471ea0 746 case kTrackWithPtFromFirstMother:
747 if (!fMCparticle) return NULL;
748 flowtrack = new AliFlowTrack();
749 flowtrack->SetPhi(fTrackPhi);
750 flowtrack->SetEta(fTrackEta);
751 tmpTParticle = fMCparticle->Particle();
752 tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
753 flowtrack->SetPt(tmpAliMCParticle->Pt());
32b846cd 754 break;
9a0783cc 755 default:
d0471ea0 756 flowtrack = new AliFlowTrack();
9a0783cc 757 flowtrack->SetPhi(fTrackPhi);
758 flowtrack->SetEta(fTrackEta);
32b846cd 759 break;
9a0783cc 760 }
12b2b8bc 761 flowtrack->SetSource(AliFlowTrack::kFromTracklet);
762 }
763 else
764 {
1ff4bda1 765 if (!fTrack) return NULL;
12b2b8bc 766 switch(fParamMix)
767 {
768 case kPure:
769 flowtrack = new AliFlowTrack(fTrack);
770 break;
771 case kTrackWithMCkine:
772 flowtrack = new AliFlowTrack(fMCparticle);
773 break;
774 case kTrackWithMCPID:
775 flowtrack = new AliFlowTrack(fTrack);
9a0783cc 776 //flowtrack->setPID(...) from mc, when implemented
12b2b8bc 777 break;
9a0783cc 778 case kTrackWithMCpt:
779 if (!fMCparticle) return NULL;
780 flowtrack = new AliFlowTrack(fTrack);
781 flowtrack->SetPt(fMCparticle->Pt());
32b846cd 782 break;
d0471ea0 783 case kTrackWithPtFromFirstMother:
784 if (!fMCparticle) return NULL;
785 flowtrack = new AliFlowTrack(fTrack);
786 tmpTParticle = fMCparticle->Particle();
787 tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
788 flowtrack->SetPt(tmpAliMCParticle->Pt());
32b846cd 789 break;
12b2b8bc 790 default:
791 flowtrack = new AliFlowTrack(fTrack);
32b846cd 792 break;
12b2b8bc 793 }
794 if (fParamType==kMC) flowtrack->SetSource(AliFlowTrack::kFromMC);
795 else if (dynamic_cast<AliESDtrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromESD);
796 else if (dynamic_cast<AliAODTrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromAOD);
797 else if (dynamic_cast<AliMCParticle*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromMC);
daf66719 798 }
6e214c87 799
daf66719 800 return flowtrack;
801}
127a5825 802
803//-----------------------------------------------------------------------
804Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
805{
806 //check if current particle is a physical primary
9a0783cc 807 if (!fMCevent) return kFALSE;
808 if (fTrackLabel<0) return kFALSE;
441ea1cf 809 return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
9a0783cc 810}
811
812//-----------------------------------------------------------------------
441ea1cf 813Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
9a0783cc 814{
815 //check if current particle is a physical primary
816 Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
9a0783cc 817 AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
818 if (!track) return kFALSE;
819 TParticle* particle = track->Particle();
820 Bool_t transported = particle->TestBit(kTransportBit);
441ea1cf 821 //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
822 //(physprim && (transported || !requiretransported))?"YES":"NO" );
823 return (physprim && (transported || !requiretransported));
127a5825 824}
12b2b8bc 825
826//-----------------------------------------------------------------------
827const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
828{
829 //return the name of the selected parameter type
830 switch (type)
831 {
832 case kMC:
833 return "MC";
834 case kGlobal:
835 return "ESD global";
836 case kESD_TPConly:
837 return "TPC only";
838 case kESD_SPDtracklet:
2a745a5f 839 return "SPD tracklet";
12b2b8bc 840 default:
2a745a5f 841 return "unknown";
12b2b8bc 842 }
12b2b8bc 843}
924b02b0 844
845//-----------------------------------------------------------------------
846void AliFlowTrackCuts::DefineHistograms()
847{
2a745a5f 848 //define qa histograms
924b02b0 849}
9a0783cc 850
851//-----------------------------------------------------------------------
852Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
853{
854 //get the number of tracks in the input event according source
855 //selection (ESD tracks, tracklets, MC particles etc.)
856 AliESDEvent* esd=NULL;
857 switch (fParamType)
858 {
859 case kESD_SPDtracklet:
860 esd = dynamic_cast<AliESDEvent*>(fEvent);
861 if (!esd) return 0;
862 return esd->GetMultiplicity()->GetNumberOfTracklets();
863 case kMC:
864 if (!fMCevent) return 0;
865 return fMCevent->GetNumberOfTracks();
866 default:
867 if (!fEvent) return 0;
868 return fEvent->GetNumberOfTracks();
869 }
870 return 0;
871}
872
873//-----------------------------------------------------------------------
874TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
875{
876 //get the input object according the data source selection:
877 //(esd tracks, traclets, mc particles,etc...)
878 AliESDEvent* esd=NULL;
879 switch (fParamType)
880 {
881 case kESD_SPDtracklet:
882 esd = dynamic_cast<AliESDEvent*>(fEvent);
883 if (!esd) return NULL;
884 return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
885 case kMC:
886 if (!fMCevent) return NULL;
887 return fMCevent->GetTrack(i);
888 default:
889 if (!fEvent) return NULL;
890 return fEvent->GetTrack(i);
891 }
892}
a1c43d26 893
894//-----------------------------------------------------------------------
895void AliFlowTrackCuts::Clear(Option_t*)
896{
897 //clean up
898 fTrack=NULL;
899 fMCevent=NULL;
900 fMCparticle=NULL;
901 fTrackLabel=0;
902 fTrackWeight=0.0;
903 fTrackEta=0.0;
904 fTrackPhi=0.0;
905}
32b846cd 906
907//-----------------------------------------------------------------------
b092f828 908Bool_t AliFlowTrackCuts::PassesTOFpidCut(AliESDtrack* track )
32b846cd 909{
910 //check if passes PID cut using timing in TOF
b092f828 911 Bool_t goodtrack = (track) &&
912 (track->GetStatus() & AliESDtrack::kTOFpid) &&
913 (track->GetTOFsignal() > 12000) &&
914 (track->GetTOFsignal() < 100000) &&
915 (track->GetIntegratedLength() > 365) &&
916 !(track->GetStatus() & AliESDtrack::kTOFmismatch);
aab6527a 917
918 if (!goodtrack) return kFALSE;
919
b092f828 920 Float_t pt = track->Pt();
921 Float_t p = track->GetP();
aab6527a 922 Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
b092f828 923 Float_t timeTOF = track->GetTOFsignal()- trackT0;
32b846cd 924 //2=pion 3=kaon 4=protons
925 Double_t inttimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
b092f828 926 track->GetIntegratedTimes(inttimes);
32b846cd 927 //construct the pid index because it's screwed up in TOF
928 Int_t pid = 0;
929 switch (fAliPID)
930 {
931 case AliPID::kPion:
932 pid=2;
933 break;
934 case AliPID::kKaon:
935 pid=3;
936 break;
937 case AliPID::kProton:
938 pid=4;
939 break;
940 default:
941 return kFALSE;
942 }
943 Float_t s = timeTOF-inttimes[pid];
944
945 Float_t* arr = fTOFpidCuts->GetMatrixArray();
946 Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
947 if (col<0) return kFALSE;
948 Float_t min = (*fTOFpidCuts)(1,col);
949 Float_t max = (*fTOFpidCuts)(2,col);
950
2a745a5f 951 //printf("--------------TOF pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
32b846cd 952 return (s>min && s<max);
953}
954
955//-----------------------------------------------------------------------
956Bool_t AliFlowTrackCuts::PassesTPCpidCut(AliESDtrack* track)
957{
958 //check if passes PID cut using dedx signal in the TPC
32b846cd 959 if (!fTPCpidCuts)
960 {
961 printf("no TPCpidCuts\n");
962 return kFALSE;
963 }
964
3c75bbd5 965 const AliExternalTrackParam* tpcparam = track->GetInnerParam(); //tpc only params at the inner wall
966 if (!tpcparam) return kFALSE;
967 Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(tpcparam->GetP(), fAliPID);
32b846cd 968 Float_t sigTPC = track->GetTPCsignal();
969 Float_t s = (sigTPC-sigExp)/sigExp;
970 Double_t pt = track->Pt();
971
972 Float_t* arr = fTPCpidCuts->GetMatrixArray();
973 Int_t col = TMath::BinarySearch(fTPCpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
974 if (col<0) return kFALSE;
975 Float_t min = (*fTPCpidCuts)(1,col);
976 Float_t max = (*fTPCpidCuts)(2,col);
977
978 //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
979 return (s>min && s<max);
980}
981
982//-----------------------------------------------------------------------
983void AliFlowTrackCuts::InitPIDcuts()
984{
985 //init matrices with PID cuts
986 TMatrixF* t = NULL;
987 if (!fTPCpidCuts)
988 {
989 if (fAliPID==AliPID::kPion)
990 {
991 t = new TMatrixF(3,10);
992 (*t)(0,0) = 0.20; (*t)(1,0) = -0.4; (*t)(2,0) = 0.2;
993 (*t)(0,1) = 0.25; (*t)(1,1) = -0.4; (*t)(2,1) = 0.2;
994 (*t)(0,2) = 0.30; (*t)(1,2) = -0.4; (*t)(2,2) = 0.25;
995 (*t)(0,3) = 0.35; (*t)(1,3) = -0.4; (*t)(2,3) = 0.25;
996 (*t)(0,4) = 0.40; (*t)(1,4) = -0.4; (*t)(2,4) = 0.3;
997 (*t)(0,5) = 0.45; (*t)(1,5) = -0.4; (*t)(2,5) = 0.3;
998 (*t)(0,6) = 0.50; (*t)(1,6) = -0.4; (*t)(2,6) = 0.25;
999 (*t)(0,7) = 0.55; (*t)(1,7) = -0.4; (*t)(2,7) = 0.15;
1000 (*t)(0,8) = 0.60; (*t)(1,8) = -0.4; (*t)(2,8) = 0.1;
1001 (*t)(0,9) = 0.65; (*t)(1,9) = 0; (*t)(2,9) = 0;
1002 }
1003 else
3ca4688a 1004 if (fAliPID==AliPID::kKaon)
1005 {
1006 t = new TMatrixF(3,7);
1007 (*t)(0,0) = 0.20; (*t)(1,0) = -0.2; (*t)(2,0) = 0.4;
1008 (*t)(0,1) = 0.25; (*t)(1,1) =-0.15; (*t)(2,1) = 0.4;
1009 (*t)(0,2) = 0.30; (*t)(1,2) = -0.1; (*t)(2,2) = 0.4;
1010 (*t)(0,3) = 0.35; (*t)(1,3) = -0.1; (*t)(2,3) = 0.4;
1011 (*t)(0,4) = 0.40; (*t)(1,4) = -0.1; (*t)(2,4) = 0.6;
1012 (*t)(0,5) = 0.45; (*t)(1,5) = -0.1; (*t)(2,5) = 0.6;
1013 (*t)(0,6) = 0.50; (*t)(1,6) = 0; (*t)(2,6) = 0;
1014 }
1015 else
32b846cd 1016 if (fAliPID==AliPID::kProton)
1017 {
1018 t = new TMatrixF(3,16);
3ca4688a 1019 (*t)(0,0) = 0.20; (*t)(1,0) = 0; (*t)(2,0) = 0;
1020 (*t)(0,1) = 0.25; (*t)(1,1) = -0.2; (*t)(2,1) = 0.3;
32b846cd 1021 (*t)(0,2) = 0.30; (*t)(1,2) = -0.2; (*t)(2,2) = 0.6;
1022 (*t)(0,3) = 0.35; (*t)(1,3) = -0.2; (*t)(2,3) = 0.6;
1023 (*t)(0,4) = 0.40; (*t)(1,4) = -0.2; (*t)(2,4) = 0.6;
1024 (*t)(0,5) = 0.45; (*t)(1,5) = -0.15; (*t)(2,5) = 0.6;
1025 (*t)(0,6) = 0.50; (*t)(1,6) = -0.1; (*t)(2,6) = 0.6;
3ca4688a 1026 (*t)(0,7) = 0.55; (*t)(1,7) = -0.05; (*t)(2,7) = 0.6;
32b846cd 1027 (*t)(0,8) = 0.60; (*t)(1,8) = -0.05; (*t)(2,8) = 0.45;
1028 (*t)(0,9) = 0.65; (*t)(1,9) = -0.05; (*t)(2,9) = 0.45;
1029 (*t)(0,10) = 0.70; (*t)(1,10) = -0.05; (*t)(2,10) = 0.45;
1030 (*t)(0,11) = 0.75; (*t)(1,11) = -0.05; (*t)(2,11) = 0.45;
1031 (*t)(0,12) = 0.80; (*t)(1,12) = 0; (*t)(2,12) = 0.45;
1032 (*t)(0,13) = 0.85; (*t)(1,13) = 0; (*t)(2,13) = 0.45;
1033 (*t)(0,14) = 0.90; (*t)(1,14) = 0; (*t)(2,14) = 0.45;
1034 (*t)(0,15) = 0.95; (*t)(1,15) = 0; (*t)(2,15) = 0;
1035 }
32b846cd 1036 fTPCpidCuts=t;
1037 }
7d9ab4fb 1038 t = NULL;
32b846cd 1039 if (!fTOFpidCuts)
1040 {
1041 if (fAliPID==AliPID::kPion)
1042 {
3abf7ecc 1043 t = new TMatrixF(3,31);
32b846cd 1044 (*t)(0,0) = 0.3; (*t)(1,0) = -700; (*t)(2,0) = 700;
1045 (*t)(0,1) = 0.35; (*t)(1,1) = -800; (*t)(2,1) = 800;
1046 (*t)(0,2) = 0.40; (*t)(1,2) = -600; (*t)(2,2) = 800;
1047 (*t)(0,3) = 0.45; (*t)(1,3) = -500; (*t)(2,3) = 700;
1048 (*t)(0,4) = 0.50; (*t)(1,4) = -400; (*t)(2,4) = 700;
1049 (*t)(0,5) = 0.55; (*t)(1,5) = -400; (*t)(2,5) = 700;
1050 (*t)(0,6) = 0.60; (*t)(1,6) = -400; (*t)(2,6) = 700;
1051 (*t)(0,7) = 0.65; (*t)(1,7) = -400; (*t)(2,7) = 700;
1052 (*t)(0,8) = 0.70; (*t)(1,8) = -400; (*t)(2,8) = 700;
1053 (*t)(0,9) = 0.75; (*t)(1,9) = -400; (*t)(2,9) = 700;
1054 (*t)(0,10) = 0.80; (*t)(1,10) = -400; (*t)(2,10) = 600;
1055 (*t)(0,11) = 0.85; (*t)(1,11) = -400; (*t)(2,11) = 600;
1056 (*t)(0,12) = 0.90; (*t)(1,12) = -400; (*t)(2,12) = 600;
1057 (*t)(0,13) = 0.95; (*t)(1,13) = -400; (*t)(2,13) = 600;
1058 (*t)(0,14) = 1.00; (*t)(1,14) = -400; (*t)(2,14) = 550;
1059 (*t)(0,15) = 1.10; (*t)(1,15) = -400; (*t)(2,15) = 450;
1060 (*t)(0,16) = 1.20; (*t)(1,16) = -400; (*t)(2,16) = 400;
1061 (*t)(0,17) = 1.30; (*t)(1,17) = -400; (*t)(2,17) = 300;
1062 (*t)(0,18) = 1.40; (*t)(1,18) = -400; (*t)(2,18) = 300;
1063 (*t)(0,19) = 1.50; (*t)(1,19) = -400; (*t)(2,19) = 250;
1064 (*t)(0,20) = 1.60; (*t)(1,20) = -400; (*t)(2,20) = 200;
1065 (*t)(0,21) = 1.70; (*t)(1,21) = -400; (*t)(2,21) = 150;
1066 (*t)(0,22) = 1.80; (*t)(1,22) = -400; (*t)(2,22) = 100;
1067 (*t)(0,23) = 1.90; (*t)(1,23) = -400; (*t)(2,23) = 70;
1068 (*t)(0,24) = 2.00; (*t)(1,24) = -400; (*t)(2,24) = 50;
441ea1cf 1069 (*t)(0,25) = 2.10; (*t)(1,25) = -400; (*t)(2,25) = 0;
3abf7ecc 1070 (*t)(0,26) = 2.20; (*t)(1,26) = -400; (*t)(2,26) = 0;
1071 (*t)(0,27) = 2.30; (*t)(1,26) = -400; (*t)(2,26) = 0;
1072 (*t)(0,28) = 2.40; (*t)(1,26) = -400; (*t)(2,26) = -50;
1073 (*t)(0,29) = 2.50; (*t)(1,26) = -400; (*t)(2,26) = -50;
1074 (*t)(0,30) = 2.60; (*t)(1,26) = 0; (*t)(2,26) = 0;
32b846cd 1075 }
1076 else
1077 if (fAliPID==AliPID::kProton)
1078 {
1079 t = new TMatrixF(3,39);
3ca4688a 1080 (*t)(0,0) = 0.3; (*t)(1,0) = 0; (*t)(2,0) = 0;
1081 (*t)(0,1) = 0.35; (*t)(1,1) = 0; (*t)(2,1) = 0;
1082 (*t)(0,2) = 0.40; (*t)(1,2) = 0; (*t)(2,2) = 0;
1083 (*t)(0,3) = 0.45; (*t)(1,3) = 0; (*t)(2,3) = 0;
1084 (*t)(0,4) = 0.50; (*t)(1,4) = 0; (*t)(2,4) = 0;
32b846cd 1085 (*t)(0,5) = 0.55; (*t)(1,5) = -900; (*t)(2,5) = 600;
1086 (*t)(0,6) = 0.60; (*t)(1,6) = -800; (*t)(2,6) = 600;
1087 (*t)(0,7) = 0.65; (*t)(1,7) = -800; (*t)(2,7) = 600;
1088 (*t)(0,8) = 0.70; (*t)(1,8) = -800; (*t)(2,8) = 600;
1089 (*t)(0,9) = 0.75; (*t)(1,9) = -700; (*t)(2,9) = 500;
1090 (*t)(0,10) = 0.80; (*t)(1,10) = -700; (*t)(2,10) = 500;
1091 (*t)(0,11) = 0.85; (*t)(1,11) = -700; (*t)(2,11) = 500;
1092 (*t)(0,12) = 0.90; (*t)(1,12) = -600; (*t)(2,12) = 500;
1093 (*t)(0,13) = 0.95; (*t)(1,13) = -600; (*t)(2,13) = 500;
1094 (*t)(0,14) = 1.00; (*t)(1,14) = -600; (*t)(2,14) = 500;
1095 (*t)(0,15) = 1.10; (*t)(1,15) = -600; (*t)(2,15) = 500;
1096 (*t)(0,16) = 1.20; (*t)(1,16) = -500; (*t)(2,16) = 500;
1097 (*t)(0,17) = 1.30; (*t)(1,17) = -500; (*t)(2,17) = 500;
1098 (*t)(0,18) = 1.40; (*t)(1,18) = -500; (*t)(2,18) = 500;
1099 (*t)(0,19) = 1.50; (*t)(1,19) = -500; (*t)(2,19) = 500;
1100 (*t)(0,20) = 1.60; (*t)(1,20) = -400; (*t)(2,20) = 500;
1101 (*t)(0,21) = 1.70; (*t)(1,21) = -400; (*t)(2,21) = 500;
1102 (*t)(0,22) = 1.80; (*t)(1,22) = -400; (*t)(2,22) = 500;
1103 (*t)(0,23) = 1.90; (*t)(1,23) = -400; (*t)(2,23) = 500;
1104 (*t)(0,24) = 2.00; (*t)(1,24) = -400; (*t)(2,24) = 500;
1105 (*t)(0,25) = 2.10; (*t)(1,25) = -350; (*t)(2,25) = 500;
1106 (*t)(0,26) = 2.20; (*t)(1,26) = -350; (*t)(2,26) = 500;
1107 (*t)(0,27) = 2.30; (*t)(1,27) = -300; (*t)(2,27) = 500;
1108 (*t)(0,28) = 2.40; (*t)(1,28) = -300; (*t)(2,28) = 500;
2a745a5f 1109 (*t)(0,29) = 2.50; (*t)(1,29) = -300; (*t)(2,29) = 500;
1110 (*t)(0,30) = 2.60; (*t)(1,30) = -250; (*t)(2,30) = 500;
1111 (*t)(0,31) = 2.70; (*t)(1,31) = -200; (*t)(2,31) = 500;
1112 (*t)(0,32) = 2.80; (*t)(1,32) = -150; (*t)(2,32) = 500;
1113 (*t)(0,33) = 2.90; (*t)(1,33) = -150; (*t)(2,33) = 500;
1114 (*t)(0,34) = 3.00; (*t)(1,34) = -100; (*t)(2,34) = 400;
1115 (*t)(0,35) = 3.10; (*t)(1,35) = -100; (*t)(2,35) = 400;
3abf7ecc 1116 (*t)(0,36) = 3.50; (*t)(1,36) = -100; (*t)(2,36) = 400;
1117 (*t)(0,37) = 3.60; (*t)(1,37) = 0; (*t)(2,37) = 0;
1118 (*t)(0,38) = 3.70; (*t)(1,38) = 0; (*t)(2,38) = 0;
32b846cd 1119 }
1120 else
1121 if (fAliPID==AliPID::kKaon)
1122 {
3abf7ecc 1123 t = new TMatrixF(3,26);
32b846cd 1124 (*t)(0,0) = 0.3; (*t)(1,0) = 0; (*t)(2,0) = 0;
1125 (*t)(0,1) = 0.35; (*t)(1,1) = 0; (*t)(2,1) = 0;
1126 (*t)(0,2) = 0.40; (*t)(1,2) = -800; (*t)(2,2) = 600;
1127 (*t)(0,3) = 0.45; (*t)(1,3) = -800; (*t)(2,3) = 600;
1128 (*t)(0,4) = 0.50; (*t)(1,4) = -800; (*t)(2,4) = 600;
1129 (*t)(0,5) = 0.55; (*t)(1,5) = -800; (*t)(2,5) = 600;
1130 (*t)(0,6) = 0.60; (*t)(1,6) = -800; (*t)(2,6) = 600;
1131 (*t)(0,7) = 0.65; (*t)(1,7) = -700; (*t)(2,7) = 600;
1132 (*t)(0,8) = 0.70; (*t)(1,8) = -600; (*t)(2,8) = 600;
1133 (*t)(0,9) = 0.75; (*t)(1,9) = -600; (*t)(2,9) = 500;
1134 (*t)(0,10) = 0.80; (*t)(1,10) = -500; (*t)(2,10) = 500;
1135 (*t)(0,11) = 0.85; (*t)(1,11) = -500; (*t)(2,11) = 500;
1136 (*t)(0,12) = 0.90; (*t)(1,12) = -400; (*t)(2,12) = 500;
1137 (*t)(0,13) = 0.95; (*t)(1,13) = -400; (*t)(2,13) = 500;
1138 (*t)(0,14) = 1.00; (*t)(1,14) = -400; (*t)(2,14) = 500;
1139 (*t)(0,15) = 1.10; (*t)(1,15) = -350; (*t)(2,15) = 450;
1140 (*t)(0,16) = 1.20; (*t)(1,16) = -300; (*t)(2,16) = 400;
1141 (*t)(0,17) = 1.30; (*t)(1,17) = -300; (*t)(2,17) = 400;
1142 (*t)(0,18) = 1.40; (*t)(1,18) = -250; (*t)(2,18) = 400;
1143 (*t)(0,19) = 1.50; (*t)(1,19) = -200; (*t)(2,19) = 400;
1144 (*t)(0,20) = 1.60; (*t)(1,20) = -150; (*t)(2,20) = 400;
1145 (*t)(0,21) = 1.70; (*t)(1,21) = -100; (*t)(2,21) = 400;
3abf7ecc 1146 (*t)(0,22) = 1.80; (*t)(1,22) = -50; (*t)(2,22) = 400;
1147 (*t)(0,23) = 1.90; (*t)(1,22) = 0; (*t)(2,22) = 400;
1148 (*t)(0,24) = 2.00; (*t)(1,22) = 50; (*t)(2,22) = 400;
1149 (*t)(0,25) = 2.10; (*t)(1,22) = 0; (*t)(2,22) = 0;
32b846cd 1150 }
1151 fTOFpidCuts=t;
1152 }
1153}
3abf7ecc 1154
1155//-----------------------------------------------------------------------
1156// part added by F. Noferini (some methods)
1157Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(AliESDtrack* track){
3abf7ecc 1158
1159 Bool_t goodtrack = track && (track->GetStatus() & AliESDtrack::kTOFpid) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000) && (track->GetIntegratedLength() > 365) && !(track->GetStatus() & AliESDtrack::kTOFmismatch);
1160
1161 if (! goodtrack)
1162 return kFALSE;
1163
1164 Int_t pdg = GetESDPdg(track,"bayesianALL");
1165 // printf("pdg set to %i\n",pdg);
1166
1167 Int_t pid = 0;
1168 Float_t prob = 0;
1169 switch (fAliPID)
1170 {
1171 case AliPID::kPion:
1172 pid=211;
1173 prob = fProbBayes[2];
1174 break;
1175 case AliPID::kKaon:
1176 pid=321;
1177 prob = fProbBayes[3];
1178 break;
1179 case AliPID::kProton:
1180 pid=2212;
1181 prob = fProbBayes[4];
1182 break;
1183 case AliPID::kElectron:
1184 pid=-11;
1185 prob = fProbBayes[0];
1186 break;
1187 default:
1188 return kFALSE;
1189 }
1190
1191 // 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);
1192 if(TMath::Abs(pdg) == TMath::Abs(pid) && prob > 0.8){
1193 if(!fCutCharge)
1194 return kTRUE;
1195 else if (fCutCharge && fCharge * track->GetSign() > 0)
1196 return kTRUE;
1197 }
1198 return kFALSE;
1199}
1200//-----------------------------------------------------------------------
1201Int_t AliFlowTrackCuts::GetESDPdg(AliESDtrack *track,Option_t *option,Int_t ipart,Float_t cPi,Float_t cKa,Float_t cPr){
1202 Int_t pdg = 0;
1203 Int_t pdgvalues[5] = {-11,-13,211,321,2212};
1204 Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
1205
1206 if(strstr(option,"bayesianTOF")){ // Bayesian TOF PID
1207 Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1208 Double_t rcc=0.;
1209
1210 Float_t pt = track->Pt();
1211
1212 Int_t iptesd = 0;
1213 while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1214
1215 if(cPi < 0){
1216 c[0] = fC[iptesd][0];
1217 c[1] = fC[iptesd][1];
1218 c[2] = fC[iptesd][2];
1219 c[3] = fC[iptesd][3];
1220 c[4] = fC[iptesd][4];
1221 }
1222 else{
1223 c[0] = 0.0;
1224 c[1] = 0.0;
1225 c[2] = cPi;
1226 c[3] = cKa;
1227 c[4] = cPr;
1228 }
1229
1230 Double_t r1[10]; track->GetTOFpid(r1);
1231
1232 Int_t i;
1233 for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1234
1235 Double_t w[10];
1236 for (i=0; i<5; i++){
1237 w[i]=c[i]*r1[i]/rcc;
1238 fProbBayes[i] = w[i];
1239 }
1240 if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1241 pdg = 211*Int_t(track->GetSign());
1242 }
1243 else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1244 pdg = 2212*Int_t(track->GetSign());
1245 }
1246 else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1247 pdg = 321*Int_t(track->GetSign());
1248 }
1249 else if (w[0]>=w[1]) { //electrons
1250 pdg = -11*Int_t(track->GetSign());
1251 }
1252 else{ // muon
1253 pdg = -13*Int_t(track->GetSign());
1254 }
1255 }
1256
1257 else if(strstr(option,"bayesianTPC")){ // Bayesian TPC PID
1258 Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1259 Double_t rcc=0.;
1260
1261 Float_t pt = track->Pt();
1262
1263 Int_t iptesd = 0;
1264 while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1265
1266 if(cPi < 0){
1267 c[0] = fC[iptesd][0];
1268 c[1] = fC[iptesd][1];
1269 c[2] = fC[iptesd][2];
1270 c[3] = fC[iptesd][3];
1271 c[4] = fC[iptesd][4];
1272 }
1273 else{
1274 c[0] = 0.0;
1275 c[1] = 0.0;
1276 c[2] = cPi;
1277 c[3] = cKa;
1278 c[4] = cPr;
1279 }
1280
1281 Double_t r1[10]; track->GetTPCpid(r1);
1282
1283 Int_t i;
1284 for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1285
1286 Double_t w[10];
1287 for (i=0; i<5; i++){
1288 w[i]=c[i]*r1[i]/rcc;
1289 fProbBayes[i] = w[i];
1290 }
1291 if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1292 pdg = 211*Int_t(track->GetSign());
1293 }
1294 else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1295 pdg = 2212*Int_t(track->GetSign());
1296 }
1297 else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1298 pdg = 321*Int_t(track->GetSign());
1299 }
1300 else if (w[0]>=w[1]) { //electrons
1301 pdg = -11*Int_t(track->GetSign());
1302 }
1303 else{ // muon
1304 pdg = -13*Int_t(track->GetSign());
1305 }
1306 }
1307
1308 else if(strstr(option,"bayesianALL")){
1309 Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1310 Double_t rcc=0.;
1311
1312 Float_t pt = track->Pt();
1313
1314 Int_t iptesd = 0;
1315 while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1316
1317 if(cPi < 0){
1318 c[0] = fC[iptesd][0];
1319 c[1] = fC[iptesd][1];
1320 c[2] = fC[iptesd][2];
1321 c[3] = fC[iptesd][3];
1322 c[4] = fC[iptesd][4];
1323 }
1324 else{
1325 c[0] = 0.0;
1326 c[1] = 0.0;
1327 c[2] = cPi;
1328 c[3] = cKa;
1329 c[4] = cPr;
1330 }
1331
1332 Double_t r1[10]; track->GetTOFpid(r1);
1333 Double_t r2[10]; track->GetTPCpid(r2);
1334
1335 Int_t i;
1336 for (i=0; i<5; i++) rcc+=(c[i]*r1[i]*r2[i]);
1337
1338
1339 Double_t w[10];
1340 for (i=0; i<5; i++){
1341 w[i]=c[i]*r1[i]*r2[i]/rcc;
1342 fProbBayes[i] = w[i];
1343 }
1344
1345 if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1346 pdg = 211*Int_t(track->GetSign());
1347 }
1348 else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1349 pdg = 2212*Int_t(track->GetSign());
1350 }
1351 else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1352 pdg = 321*Int_t(track->GetSign());
1353 }
1354 else if (w[0]>=w[1]) { //electrons
1355 pdg = -11*Int_t(track->GetSign());
1356 }
1357 else{ // muon
1358 pdg = -13*Int_t(track->GetSign());
1359 }
1360 }
1361
1362 else if(strstr(option,"sigmacutTOF")){
1363 printf("PID not implemented yet: %s\nNO PID!!!!\n",option);
1364 Float_t p = track->P();
1365
1366 // Take expected times
1367 Double_t exptimes[5];
1368 track->GetIntegratedTimes(exptimes);
1369
1370 // Take resolution for TOF response
aab6527a 1371 // like fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1372 Float_t resolution = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
3abf7ecc 1373
1374 if(TMath::Abs(exptimes[ipart] - track->GetTOFsignal()) < 3 * resolution){
1375 pdg = pdgvalues[ipart] * Int_t(track->GetSign());
1376 }
1377 }
1378
1379 else{
1380 printf("Invalid PID option: %s\nNO PID!!!!\n",option);
1381 }
1382
1383 return pdg;
1384}
1385//-----------------------------------------------------------------------
1386void AliFlowTrackCuts::SetPriors(){
1387 // set abbundancies
1388 fBinLimitPID[0] = 0.30;
1389 fC[0][0] = 0.015;
1390 fC[0][1] = 0.015;
1391 fC[0][2] = 1;
1392 fC[0][3] = 0.0025;
1393 fC[0][4] = 0.000015;
1394 fBinLimitPID[1] = 0.35;
1395 fC[1][0] = 0.015;
1396 fC[1][1] = 0.015;
1397 fC[1][2] = 1;
1398 fC[1][3] = 0.01;
1399 fC[1][4] = 0.001;
1400 fBinLimitPID[2] = 0.40;
1401 fC[2][0] = 0.015;
1402 fC[2][1] = 0.015;
1403 fC[2][2] = 1;
1404 fC[2][3] = 0.026;
1405 fC[2][4] = 0.004;
1406 fBinLimitPID[3] = 0.45;
1407 fC[3][0] = 0.015;
1408 fC[3][1] = 0.015;
1409 fC[3][2] = 1;
1410 fC[3][3] = 0.026;
1411 fC[3][4] = 0.004;
1412 fBinLimitPID[4] = 0.50;
1413 fC[4][0] = 0.015;
1414 fC[4][1] = 0.015;
1415 fC[4][2] = 1.000000;
1416 fC[4][3] = 0.05;
1417 fC[4][4] = 0.01;
1418 fBinLimitPID[5] = 0.60;
1419 fC[5][0] = 0.012;
1420 fC[5][1] = 0.012;
1421 fC[5][2] = 1;
1422 fC[5][3] = 0.085;
1423 fC[5][4] = 0.022;
1424 fBinLimitPID[6] = 0.70;
1425 fC[6][0] = 0.01;
1426 fC[6][1] = 0.01;
1427 fC[6][2] = 1;
1428 fC[6][3] = 0.12;
1429 fC[6][4] = 0.036;
1430 fBinLimitPID[7] = 0.80;
1431 fC[7][0] = 0.0095;
1432 fC[7][1] = 0.0095;
1433 fC[7][2] = 1;
1434 fC[7][3] = 0.15;
1435 fC[7][4] = 0.05;
1436 fBinLimitPID[8] = 0.90;
1437 fC[8][0] = 0.0085;
1438 fC[8][1] = 0.0085;
1439 fC[8][2] = 1;
1440 fC[8][3] = 0.18;
1441 fC[8][4] = 0.074;
1442 fBinLimitPID[9] = 1;
1443 fC[9][0] = 0.008;
1444 fC[9][1] = 0.008;
1445 fC[9][2] = 1;
1446 fC[9][3] = 0.22;
1447 fC[9][4] = 0.1;
1448 fBinLimitPID[10] = 1.20;
1449 fC[10][0] = 0.007;
1450 fC[10][1] = 0.007;
1451 fC[10][2] = 1;
1452 fC[10][3] = 0.28;
1453 fC[10][4] = 0.16;
1454 fBinLimitPID[11] = 1.40;
1455 fC[11][0] = 0.0066;
1456 fC[11][1] = 0.0066;
1457 fC[11][2] = 1;
1458 fC[11][3] = 0.35;
1459 fC[11][4] = 0.23;
1460 fBinLimitPID[12] = 1.60;
1461 fC[12][0] = 0.0075;
1462 fC[12][1] = 0.0075;
1463 fC[12][2] = 1;
1464 fC[12][3] = 0.40;
1465 fC[12][4] = 0.31;
1466 fBinLimitPID[13] = 1.80;
1467 fC[13][0] = 0.0062;
1468 fC[13][1] = 0.0062;
1469 fC[13][2] = 1;
1470 fC[13][3] = 0.45;
1471 fC[13][4] = 0.39;
1472 fBinLimitPID[14] = 2.00;
1473 fC[14][0] = 0.005;
1474 fC[14][1] = 0.005;
1475 fC[14][2] = 1;
1476 fC[14][3] = 0.46;
1477 fC[14][4] = 0.47;
1478 fBinLimitPID[15] = 2.20;
1479 fC[15][0] = 0.0042;
1480 fC[15][1] = 0.0042;
1481 fC[15][2] = 1;
1482 fC[15][3] = 0.5;
1483 fC[15][4] = 0.55;
1484 fBinLimitPID[16] = 2.40;
1485 fC[16][0] = 0.007;
1486 fC[16][1] = 0.007;
1487 fC[16][2] = 1;
1488 fC[16][3] = 0.5;
1489 fC[16][4] = 0.6;
1490
1491 for(Int_t i=17;i<fnPIDptBin;i++){
1492 fBinLimitPID[i] = 2.0 + 0.2 * (i-14);
1493 fC[i][0] = fC[13][0];
1494 fC[i][1] = fC[13][1];
1495 fC[i][2] = fC[13][2];
1496 fC[i][3] = fC[13][3];
1497 fC[i][4] = fC[13][4];
1498 }
1499}
1500// end part added by F. Noferini
1501//-----------------------------------------------------------------------