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