some coverty fixes
[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
422c61c7 111 SetPriors(); //init arrays
441ea1cf 112}
113
114//-----------------------------------------------------------------------
115AliFlowTrackCuts::AliFlowTrackCuts(const char* name):
116 AliFlowTrackSimpleCuts(),
daf66719 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
308//-----------------------------------------------------------------------
aab6527a 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
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
360//-----------------------------------------------------------------------
1ff4bda1 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
380//-----------------------------------------------------------------------
daf66719 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
390//-----------------------------------------------------------------------
12b2b8bc 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
462//-----------------------------------------------------------------------
daf66719 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
daf66719 659 switch (fParamType)
660 {
12b2b8bc 661 case kGlobal:
daf66719 662 fTrack = track;
daf66719 663 break;
664 case kESD_TPConly:
1ff4bda1 665 if (!track->FillTPCOnlyTrack(fTPCtrack))
666 {
667 fTrack=NULL;
668 fMCparticle=NULL;
669 fTrackLabel=-1;
670 return;
671 }
672 fTrack = &fTPCtrack;
957517fa 673 //recalculate the label and mc particle, they may differ as TPClabel != global label
127a5825 674 fTrackLabel = (fFakesAreOK)?TMath::Abs(fTrack->GetLabel()):fTrack->GetLabel();
675 if (fMCevent) fMCparticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(fTrackLabel));
957517fa 676 else fMCparticle=NULL;
daf66719 677 break;
daf66719 678 default:
679 fTrack = track;
32b846cd 680 break;
daf66719 681 }
682}
683
684//-----------------------------------------------------------------------
685AliFlowTrackCuts* AliFlowTrackCuts::GetStandardTPCOnlyTrackCuts()
686{
687 //get standard cuts
441ea1cf 688 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard TPConly cuts");
daf66719 689 delete cuts->fAliESDtrackCuts;
690 cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
5559ce24 691 cuts->SetParamType(kESD_TPConly);
daf66719 692 return cuts;
693}
694
695//-----------------------------------------------------------------------
696AliFlowTrackCuts* AliFlowTrackCuts::GetStandardITSTPCTrackCuts2009(Bool_t selPrimaries)
697{
698 //get standard cuts
441ea1cf 699 AliFlowTrackCuts* cuts = new AliFlowTrackCuts("standard global track cuts 2009");
daf66719 700 delete cuts->fAliESDtrackCuts;
701 cuts->fAliESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selPrimaries);
12b2b8bc 702 cuts->SetParamType(kGlobal);
daf66719 703 return cuts;
704}
705
706//-----------------------------------------------------------------------
6e214c87 707AliFlowTrack* AliFlowTrackCuts::MakeFlowTrack(int index) const
daf66719 708{
709 //get a flow track constructed from whatever we applied cuts on
710 //caller is resposible for deletion
9a0783cc 711 //if construction fails return NULL
daf66719 712 AliFlowTrack* flowtrack=NULL;
d0471ea0 713 TParticle *tmpTParticle=NULL;
714 AliMCParticle* tmpAliMCParticle=NULL;
12b2b8bc 715 if (fParamType==kESD_SPDtracklet)
daf66719 716 {
9a0783cc 717 switch (fParamMix)
718 {
719 case kPure:
d0471ea0 720 flowtrack = new AliFlowTrack();
9a0783cc 721 flowtrack->SetPhi(fTrackPhi);
722 flowtrack->SetEta(fTrackEta);
723 break;
724 case kTrackWithMCkine:
725 if (!fMCparticle) return NULL;
d0471ea0 726 flowtrack = new AliFlowTrack();
9a0783cc 727 flowtrack->SetPhi( fMCparticle->Phi() );
728 flowtrack->SetEta( fMCparticle->Eta() );
729 flowtrack->SetPt( fMCparticle->Pt() );
730 break;
731 case kTrackWithMCpt:
732 if (!fMCparticle) return NULL;
d0471ea0 733 flowtrack = new AliFlowTrack();
9a0783cc 734 flowtrack->SetPhi(fTrackPhi);
735 flowtrack->SetEta(fTrackEta);
736 flowtrack->SetPt(fMCparticle->Pt());
737 break;
d0471ea0 738 case kTrackWithPtFromFirstMother:
739 if (!fMCparticle) return NULL;
740 flowtrack = new AliFlowTrack();
741 flowtrack->SetPhi(fTrackPhi);
742 flowtrack->SetEta(fTrackEta);
743 tmpTParticle = fMCparticle->Particle();
744 tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
745 flowtrack->SetPt(tmpAliMCParticle->Pt());
32b846cd 746 break;
9a0783cc 747 default:
d0471ea0 748 flowtrack = new AliFlowTrack();
9a0783cc 749 flowtrack->SetPhi(fTrackPhi);
750 flowtrack->SetEta(fTrackEta);
32b846cd 751 break;
9a0783cc 752 }
12b2b8bc 753 flowtrack->SetSource(AliFlowTrack::kFromTracklet);
754 }
755 else
756 {
1ff4bda1 757 if (!fTrack) return NULL;
12b2b8bc 758 switch(fParamMix)
759 {
760 case kPure:
761 flowtrack = new AliFlowTrack(fTrack);
762 break;
763 case kTrackWithMCkine:
764 flowtrack = new AliFlowTrack(fMCparticle);
765 break;
766 case kTrackWithMCPID:
767 flowtrack = new AliFlowTrack(fTrack);
9a0783cc 768 //flowtrack->setPID(...) from mc, when implemented
12b2b8bc 769 break;
9a0783cc 770 case kTrackWithMCpt:
771 if (!fMCparticle) return NULL;
772 flowtrack = new AliFlowTrack(fTrack);
773 flowtrack->SetPt(fMCparticle->Pt());
32b846cd 774 break;
d0471ea0 775 case kTrackWithPtFromFirstMother:
776 if (!fMCparticle) return NULL;
777 flowtrack = new AliFlowTrack(fTrack);
778 tmpTParticle = fMCparticle->Particle();
779 tmpAliMCParticle = static_cast<AliMCParticle*>(fMCevent->GetTrack(tmpTParticle->GetFirstMother()));
780 flowtrack->SetPt(tmpAliMCParticle->Pt());
32b846cd 781 break;
12b2b8bc 782 default:
783 flowtrack = new AliFlowTrack(fTrack);
32b846cd 784 break;
12b2b8bc 785 }
786 if (fParamType==kMC) flowtrack->SetSource(AliFlowTrack::kFromMC);
787 else if (dynamic_cast<AliESDtrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromESD);
788 else if (dynamic_cast<AliAODTrack*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromAOD);
789 else if (dynamic_cast<AliMCParticle*>(fTrack)) flowtrack->SetSource(AliFlowTrack::kFromMC);
daf66719 790 }
6e214c87 791
792 if(flowtrack)
793 flowtrack->SetIndexOnCollection(index);
daf66719 794 return flowtrack;
795}
127a5825 796
797//-----------------------------------------------------------------------
798Bool_t AliFlowTrackCuts::IsPhysicalPrimary() const
799{
800 //check if current particle is a physical primary
9a0783cc 801 if (!fMCevent) return kFALSE;
802 if (fTrackLabel<0) return kFALSE;
441ea1cf 803 return IsPhysicalPrimary(fMCevent, fTrackLabel, fRequireTransportBitForPrimaries);
9a0783cc 804}
805
806//-----------------------------------------------------------------------
441ea1cf 807Bool_t AliFlowTrackCuts::IsPhysicalPrimary(AliMCEvent* mcEvent, Int_t label, Bool_t requiretransported)
9a0783cc 808{
809 //check if current particle is a physical primary
810 Bool_t physprim=mcEvent->IsPhysicalPrimary(label);
9a0783cc 811 AliMCParticle* track = static_cast<AliMCParticle*>(mcEvent->GetTrack(label));
812 if (!track) return kFALSE;
813 TParticle* particle = track->Particle();
814 Bool_t transported = particle->TestBit(kTransportBit);
441ea1cf 815 //printf("label: %i prim: %s, transp: %s, pass: %s\n",label, (physprim)?"YES":"NO ",(transported)?"YES":"NO ",
816 //(physprim && (transported || !requiretransported))?"YES":"NO" );
817 return (physprim && (transported || !requiretransported));
127a5825 818}
12b2b8bc 819
820//-----------------------------------------------------------------------
821const char* AliFlowTrackCuts::GetParamTypeName(trackParameterType type)
822{
823 //return the name of the selected parameter type
824 switch (type)
825 {
826 case kMC:
827 return "MC";
828 case kGlobal:
829 return "ESD global";
830 case kESD_TPConly:
831 return "TPC only";
832 case kESD_SPDtracklet:
2a745a5f 833 return "SPD tracklet";
12b2b8bc 834 default:
2a745a5f 835 return "unknown";
12b2b8bc 836 }
12b2b8bc 837}
924b02b0 838
839//-----------------------------------------------------------------------
840void AliFlowTrackCuts::DefineHistograms()
841{
2a745a5f 842 //define qa histograms
924b02b0 843}
9a0783cc 844
845//-----------------------------------------------------------------------
846Int_t AliFlowTrackCuts::GetNumberOfInputObjects() const
847{
848 //get the number of tracks in the input event according source
849 //selection (ESD tracks, tracklets, MC particles etc.)
850 AliESDEvent* esd=NULL;
851 switch (fParamType)
852 {
853 case kESD_SPDtracklet:
854 esd = dynamic_cast<AliESDEvent*>(fEvent);
855 if (!esd) return 0;
856 return esd->GetMultiplicity()->GetNumberOfTracklets();
857 case kMC:
858 if (!fMCevent) return 0;
859 return fMCevent->GetNumberOfTracks();
860 default:
861 if (!fEvent) return 0;
862 return fEvent->GetNumberOfTracks();
863 }
864 return 0;
865}
866
867//-----------------------------------------------------------------------
868TObject* AliFlowTrackCuts::GetInputObject(Int_t i)
869{
870 //get the input object according the data source selection:
871 //(esd tracks, traclets, mc particles,etc...)
872 AliESDEvent* esd=NULL;
873 switch (fParamType)
874 {
875 case kESD_SPDtracklet:
876 esd = dynamic_cast<AliESDEvent*>(fEvent);
877 if (!esd) return NULL;
878 return const_cast<AliMultiplicity*>(esd->GetMultiplicity());
879 case kMC:
880 if (!fMCevent) return NULL;
881 return fMCevent->GetTrack(i);
882 default:
883 if (!fEvent) return NULL;
884 return fEvent->GetTrack(i);
885 }
886}
a1c43d26 887
888//-----------------------------------------------------------------------
889void AliFlowTrackCuts::Clear(Option_t*)
890{
891 //clean up
892 fTrack=NULL;
893 fMCevent=NULL;
894 fMCparticle=NULL;
895 fTrackLabel=0;
896 fTrackWeight=0.0;
897 fTrackEta=0.0;
898 fTrackPhi=0.0;
899}
32b846cd 900
901//-----------------------------------------------------------------------
902Bool_t AliFlowTrackCuts::PassesTOFpidCut(AliESDtrack* t )
903{
904 //check if passes PID cut using timing in TOF
aab6527a 905 Bool_t goodtrack = (t) &&
906 (t->GetStatus() & AliESDtrack::kTOFpid) &&
907 (t->GetTOFsignal() > 12000) &&
908 (t->GetTOFsignal() < 100000) &&
909 (t->GetIntegratedLength() > 365) &&
910 !(t->GetStatus() & AliESDtrack::kTOFmismatch);
911
912 if (!goodtrack) return kFALSE;
913
32b846cd 914 Float_t pt = t->Pt();
915 Float_t p = t->GetP();
aab6527a 916 Float_t trackT0 = fESDpid.GetTOFResponse().GetStartTime(p);
32b846cd 917 Float_t timeTOF = t->GetTOFsignal()- trackT0;
918 //2=pion 3=kaon 4=protons
919 Double_t inttimes[5] = {-1.0,-1.0,-1.0,-1.0,-1.0};
920 t->GetIntegratedTimes(inttimes);
921 //construct the pid index because it's screwed up in TOF
922 Int_t pid = 0;
923 switch (fAliPID)
924 {
925 case AliPID::kPion:
926 pid=2;
927 break;
928 case AliPID::kKaon:
929 pid=3;
930 break;
931 case AliPID::kProton:
932 pid=4;
933 break;
934 default:
935 return kFALSE;
936 }
937 Float_t s = timeTOF-inttimes[pid];
938
939 Float_t* arr = fTOFpidCuts->GetMatrixArray();
940 Int_t col = TMath::BinarySearch(fTOFpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
941 if (col<0) return kFALSE;
942 Float_t min = (*fTOFpidCuts)(1,col);
943 Float_t max = (*fTOFpidCuts)(2,col);
944
2a745a5f 945 //printf("--------------TOF pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
32b846cd 946 return (s>min && s<max);
947}
948
949//-----------------------------------------------------------------------
950Bool_t AliFlowTrackCuts::PassesTPCpidCut(AliESDtrack* track)
951{
952 //check if passes PID cut using dedx signal in the TPC
32b846cd 953 if (!fTPCpidCuts)
954 {
955 printf("no TPCpidCuts\n");
956 return kFALSE;
957 }
958
959 const AliExternalTrackParam* tpcparam = track->GetInnerParam();
960 if (!tpcparam) return kFALSE;
aab6527a 961 Float_t sigExp = fESDpid.GetTPCResponse().GetExpectedSignal(tpcparam->GetP(), fAliPID);
32b846cd 962 Float_t sigTPC = track->GetTPCsignal();
963 Float_t s = (sigTPC-sigExp)/sigExp;
964 Double_t pt = track->Pt();
965
966 Float_t* arr = fTPCpidCuts->GetMatrixArray();
967 Int_t col = TMath::BinarySearch(fTPCpidCuts->GetNcols(),arr,static_cast<Float_t>(pt));
968 if (col<0) return kFALSE;
969 Float_t min = (*fTPCpidCuts)(1,col);
970 Float_t max = (*fTPCpidCuts)(2,col);
971
972 //printf("------------TPC pid cut %s\n",(s>min && s<max)?"PASS":"FAIL");
973 return (s>min && s<max);
974}
975
976//-----------------------------------------------------------------------
977void AliFlowTrackCuts::InitPIDcuts()
978{
979 //init matrices with PID cuts
980 TMatrixF* t = NULL;
981 if (!fTPCpidCuts)
982 {
983 if (fAliPID==AliPID::kPion)
984 {
985 t = new TMatrixF(3,10);
986 (*t)(0,0) = 0.20; (*t)(1,0) = -0.4; (*t)(2,0) = 0.2;
987 (*t)(0,1) = 0.25; (*t)(1,1) = -0.4; (*t)(2,1) = 0.2;
988 (*t)(0,2) = 0.30; (*t)(1,2) = -0.4; (*t)(2,2) = 0.25;
989 (*t)(0,3) = 0.35; (*t)(1,3) = -0.4; (*t)(2,3) = 0.25;
990 (*t)(0,4) = 0.40; (*t)(1,4) = -0.4; (*t)(2,4) = 0.3;
991 (*t)(0,5) = 0.45; (*t)(1,5) = -0.4; (*t)(2,5) = 0.3;
992 (*t)(0,6) = 0.50; (*t)(1,6) = -0.4; (*t)(2,6) = 0.25;
993 (*t)(0,7) = 0.55; (*t)(1,7) = -0.4; (*t)(2,7) = 0.15;
994 (*t)(0,8) = 0.60; (*t)(1,8) = -0.4; (*t)(2,8) = 0.1;
995 (*t)(0,9) = 0.65; (*t)(1,9) = 0; (*t)(2,9) = 0;
996 }
997 else
3ca4688a 998 if (fAliPID==AliPID::kKaon)
999 {
1000 t = new TMatrixF(3,7);
1001 (*t)(0,0) = 0.20; (*t)(1,0) = -0.2; (*t)(2,0) = 0.4;
1002 (*t)(0,1) = 0.25; (*t)(1,1) =-0.15; (*t)(2,1) = 0.4;
1003 (*t)(0,2) = 0.30; (*t)(1,2) = -0.1; (*t)(2,2) = 0.4;
1004 (*t)(0,3) = 0.35; (*t)(1,3) = -0.1; (*t)(2,3) = 0.4;
1005 (*t)(0,4) = 0.40; (*t)(1,4) = -0.1; (*t)(2,4) = 0.6;
1006 (*t)(0,5) = 0.45; (*t)(1,5) = -0.1; (*t)(2,5) = 0.6;
1007 (*t)(0,6) = 0.50; (*t)(1,6) = 0; (*t)(2,6) = 0;
1008 }
1009 else
32b846cd 1010 if (fAliPID==AliPID::kProton)
1011 {
1012 t = new TMatrixF(3,16);
3ca4688a 1013 (*t)(0,0) = 0.20; (*t)(1,0) = 0; (*t)(2,0) = 0;
1014 (*t)(0,1) = 0.25; (*t)(1,1) = -0.2; (*t)(2,1) = 0.3;
32b846cd 1015 (*t)(0,2) = 0.30; (*t)(1,2) = -0.2; (*t)(2,2) = 0.6;
1016 (*t)(0,3) = 0.35; (*t)(1,3) = -0.2; (*t)(2,3) = 0.6;
1017 (*t)(0,4) = 0.40; (*t)(1,4) = -0.2; (*t)(2,4) = 0.6;
1018 (*t)(0,5) = 0.45; (*t)(1,5) = -0.15; (*t)(2,5) = 0.6;
1019 (*t)(0,6) = 0.50; (*t)(1,6) = -0.1; (*t)(2,6) = 0.6;
3ca4688a 1020 (*t)(0,7) = 0.55; (*t)(1,7) = -0.05; (*t)(2,7) = 0.6;
32b846cd 1021 (*t)(0,8) = 0.60; (*t)(1,8) = -0.05; (*t)(2,8) = 0.45;
1022 (*t)(0,9) = 0.65; (*t)(1,9) = -0.05; (*t)(2,9) = 0.45;
1023 (*t)(0,10) = 0.70; (*t)(1,10) = -0.05; (*t)(2,10) = 0.45;
1024 (*t)(0,11) = 0.75; (*t)(1,11) = -0.05; (*t)(2,11) = 0.45;
1025 (*t)(0,12) = 0.80; (*t)(1,12) = 0; (*t)(2,12) = 0.45;
1026 (*t)(0,13) = 0.85; (*t)(1,13) = 0; (*t)(2,13) = 0.45;
1027 (*t)(0,14) = 0.90; (*t)(1,14) = 0; (*t)(2,14) = 0.45;
1028 (*t)(0,15) = 0.95; (*t)(1,15) = 0; (*t)(2,15) = 0;
1029 }
32b846cd 1030 fTPCpidCuts=t;
1031 }
7d9ab4fb 1032 t = NULL;
32b846cd 1033 if (!fTOFpidCuts)
1034 {
1035 if (fAliPID==AliPID::kPion)
1036 {
3abf7ecc 1037 t = new TMatrixF(3,31);
32b846cd 1038 (*t)(0,0) = 0.3; (*t)(1,0) = -700; (*t)(2,0) = 700;
1039 (*t)(0,1) = 0.35; (*t)(1,1) = -800; (*t)(2,1) = 800;
1040 (*t)(0,2) = 0.40; (*t)(1,2) = -600; (*t)(2,2) = 800;
1041 (*t)(0,3) = 0.45; (*t)(1,3) = -500; (*t)(2,3) = 700;
1042 (*t)(0,4) = 0.50; (*t)(1,4) = -400; (*t)(2,4) = 700;
1043 (*t)(0,5) = 0.55; (*t)(1,5) = -400; (*t)(2,5) = 700;
1044 (*t)(0,6) = 0.60; (*t)(1,6) = -400; (*t)(2,6) = 700;
1045 (*t)(0,7) = 0.65; (*t)(1,7) = -400; (*t)(2,7) = 700;
1046 (*t)(0,8) = 0.70; (*t)(1,8) = -400; (*t)(2,8) = 700;
1047 (*t)(0,9) = 0.75; (*t)(1,9) = -400; (*t)(2,9) = 700;
1048 (*t)(0,10) = 0.80; (*t)(1,10) = -400; (*t)(2,10) = 600;
1049 (*t)(0,11) = 0.85; (*t)(1,11) = -400; (*t)(2,11) = 600;
1050 (*t)(0,12) = 0.90; (*t)(1,12) = -400; (*t)(2,12) = 600;
1051 (*t)(0,13) = 0.95; (*t)(1,13) = -400; (*t)(2,13) = 600;
1052 (*t)(0,14) = 1.00; (*t)(1,14) = -400; (*t)(2,14) = 550;
1053 (*t)(0,15) = 1.10; (*t)(1,15) = -400; (*t)(2,15) = 450;
1054 (*t)(0,16) = 1.20; (*t)(1,16) = -400; (*t)(2,16) = 400;
1055 (*t)(0,17) = 1.30; (*t)(1,17) = -400; (*t)(2,17) = 300;
1056 (*t)(0,18) = 1.40; (*t)(1,18) = -400; (*t)(2,18) = 300;
1057 (*t)(0,19) = 1.50; (*t)(1,19) = -400; (*t)(2,19) = 250;
1058 (*t)(0,20) = 1.60; (*t)(1,20) = -400; (*t)(2,20) = 200;
1059 (*t)(0,21) = 1.70; (*t)(1,21) = -400; (*t)(2,21) = 150;
1060 (*t)(0,22) = 1.80; (*t)(1,22) = -400; (*t)(2,22) = 100;
1061 (*t)(0,23) = 1.90; (*t)(1,23) = -400; (*t)(2,23) = 70;
1062 (*t)(0,24) = 2.00; (*t)(1,24) = -400; (*t)(2,24) = 50;
441ea1cf 1063 (*t)(0,25) = 2.10; (*t)(1,25) = -400; (*t)(2,25) = 0;
3abf7ecc 1064 (*t)(0,26) = 2.20; (*t)(1,26) = -400; (*t)(2,26) = 0;
1065 (*t)(0,27) = 2.30; (*t)(1,26) = -400; (*t)(2,26) = 0;
1066 (*t)(0,28) = 2.40; (*t)(1,26) = -400; (*t)(2,26) = -50;
1067 (*t)(0,29) = 2.50; (*t)(1,26) = -400; (*t)(2,26) = -50;
1068 (*t)(0,30) = 2.60; (*t)(1,26) = 0; (*t)(2,26) = 0;
32b846cd 1069 }
1070 else
1071 if (fAliPID==AliPID::kProton)
1072 {
1073 t = new TMatrixF(3,39);
3ca4688a 1074 (*t)(0,0) = 0.3; (*t)(1,0) = 0; (*t)(2,0) = 0;
1075 (*t)(0,1) = 0.35; (*t)(1,1) = 0; (*t)(2,1) = 0;
1076 (*t)(0,2) = 0.40; (*t)(1,2) = 0; (*t)(2,2) = 0;
1077 (*t)(0,3) = 0.45; (*t)(1,3) = 0; (*t)(2,3) = 0;
1078 (*t)(0,4) = 0.50; (*t)(1,4) = 0; (*t)(2,4) = 0;
32b846cd 1079 (*t)(0,5) = 0.55; (*t)(1,5) = -900; (*t)(2,5) = 600;
1080 (*t)(0,6) = 0.60; (*t)(1,6) = -800; (*t)(2,6) = 600;
1081 (*t)(0,7) = 0.65; (*t)(1,7) = -800; (*t)(2,7) = 600;
1082 (*t)(0,8) = 0.70; (*t)(1,8) = -800; (*t)(2,8) = 600;
1083 (*t)(0,9) = 0.75; (*t)(1,9) = -700; (*t)(2,9) = 500;
1084 (*t)(0,10) = 0.80; (*t)(1,10) = -700; (*t)(2,10) = 500;
1085 (*t)(0,11) = 0.85; (*t)(1,11) = -700; (*t)(2,11) = 500;
1086 (*t)(0,12) = 0.90; (*t)(1,12) = -600; (*t)(2,12) = 500;
1087 (*t)(0,13) = 0.95; (*t)(1,13) = -600; (*t)(2,13) = 500;
1088 (*t)(0,14) = 1.00; (*t)(1,14) = -600; (*t)(2,14) = 500;
1089 (*t)(0,15) = 1.10; (*t)(1,15) = -600; (*t)(2,15) = 500;
1090 (*t)(0,16) = 1.20; (*t)(1,16) = -500; (*t)(2,16) = 500;
1091 (*t)(0,17) = 1.30; (*t)(1,17) = -500; (*t)(2,17) = 500;
1092 (*t)(0,18) = 1.40; (*t)(1,18) = -500; (*t)(2,18) = 500;
1093 (*t)(0,19) = 1.50; (*t)(1,19) = -500; (*t)(2,19) = 500;
1094 (*t)(0,20) = 1.60; (*t)(1,20) = -400; (*t)(2,20) = 500;
1095 (*t)(0,21) = 1.70; (*t)(1,21) = -400; (*t)(2,21) = 500;
1096 (*t)(0,22) = 1.80; (*t)(1,22) = -400; (*t)(2,22) = 500;
1097 (*t)(0,23) = 1.90; (*t)(1,23) = -400; (*t)(2,23) = 500;
1098 (*t)(0,24) = 2.00; (*t)(1,24) = -400; (*t)(2,24) = 500;
1099 (*t)(0,25) = 2.10; (*t)(1,25) = -350; (*t)(2,25) = 500;
1100 (*t)(0,26) = 2.20; (*t)(1,26) = -350; (*t)(2,26) = 500;
1101 (*t)(0,27) = 2.30; (*t)(1,27) = -300; (*t)(2,27) = 500;
1102 (*t)(0,28) = 2.40; (*t)(1,28) = -300; (*t)(2,28) = 500;
2a745a5f 1103 (*t)(0,29) = 2.50; (*t)(1,29) = -300; (*t)(2,29) = 500;
1104 (*t)(0,30) = 2.60; (*t)(1,30) = -250; (*t)(2,30) = 500;
1105 (*t)(0,31) = 2.70; (*t)(1,31) = -200; (*t)(2,31) = 500;
1106 (*t)(0,32) = 2.80; (*t)(1,32) = -150; (*t)(2,32) = 500;
1107 (*t)(0,33) = 2.90; (*t)(1,33) = -150; (*t)(2,33) = 500;
1108 (*t)(0,34) = 3.00; (*t)(1,34) = -100; (*t)(2,34) = 400;
1109 (*t)(0,35) = 3.10; (*t)(1,35) = -100; (*t)(2,35) = 400;
3abf7ecc 1110 (*t)(0,36) = 3.50; (*t)(1,36) = -100; (*t)(2,36) = 400;
1111 (*t)(0,37) = 3.60; (*t)(1,37) = 0; (*t)(2,37) = 0;
1112 (*t)(0,38) = 3.70; (*t)(1,38) = 0; (*t)(2,38) = 0;
32b846cd 1113 }
1114 else
1115 if (fAliPID==AliPID::kKaon)
1116 {
3abf7ecc 1117 t = new TMatrixF(3,26);
32b846cd 1118 (*t)(0,0) = 0.3; (*t)(1,0) = 0; (*t)(2,0) = 0;
1119 (*t)(0,1) = 0.35; (*t)(1,1) = 0; (*t)(2,1) = 0;
1120 (*t)(0,2) = 0.40; (*t)(1,2) = -800; (*t)(2,2) = 600;
1121 (*t)(0,3) = 0.45; (*t)(1,3) = -800; (*t)(2,3) = 600;
1122 (*t)(0,4) = 0.50; (*t)(1,4) = -800; (*t)(2,4) = 600;
1123 (*t)(0,5) = 0.55; (*t)(1,5) = -800; (*t)(2,5) = 600;
1124 (*t)(0,6) = 0.60; (*t)(1,6) = -800; (*t)(2,6) = 600;
1125 (*t)(0,7) = 0.65; (*t)(1,7) = -700; (*t)(2,7) = 600;
1126 (*t)(0,8) = 0.70; (*t)(1,8) = -600; (*t)(2,8) = 600;
1127 (*t)(0,9) = 0.75; (*t)(1,9) = -600; (*t)(2,9) = 500;
1128 (*t)(0,10) = 0.80; (*t)(1,10) = -500; (*t)(2,10) = 500;
1129 (*t)(0,11) = 0.85; (*t)(1,11) = -500; (*t)(2,11) = 500;
1130 (*t)(0,12) = 0.90; (*t)(1,12) = -400; (*t)(2,12) = 500;
1131 (*t)(0,13) = 0.95; (*t)(1,13) = -400; (*t)(2,13) = 500;
1132 (*t)(0,14) = 1.00; (*t)(1,14) = -400; (*t)(2,14) = 500;
1133 (*t)(0,15) = 1.10; (*t)(1,15) = -350; (*t)(2,15) = 450;
1134 (*t)(0,16) = 1.20; (*t)(1,16) = -300; (*t)(2,16) = 400;
1135 (*t)(0,17) = 1.30; (*t)(1,17) = -300; (*t)(2,17) = 400;
1136 (*t)(0,18) = 1.40; (*t)(1,18) = -250; (*t)(2,18) = 400;
1137 (*t)(0,19) = 1.50; (*t)(1,19) = -200; (*t)(2,19) = 400;
1138 (*t)(0,20) = 1.60; (*t)(1,20) = -150; (*t)(2,20) = 400;
1139 (*t)(0,21) = 1.70; (*t)(1,21) = -100; (*t)(2,21) = 400;
3abf7ecc 1140 (*t)(0,22) = 1.80; (*t)(1,22) = -50; (*t)(2,22) = 400;
1141 (*t)(0,23) = 1.90; (*t)(1,22) = 0; (*t)(2,22) = 400;
1142 (*t)(0,24) = 2.00; (*t)(1,22) = 50; (*t)(2,22) = 400;
1143 (*t)(0,25) = 2.10; (*t)(1,22) = 0; (*t)(2,22) = 0;
32b846cd 1144 }
1145 fTOFpidCuts=t;
1146 }
1147}
3abf7ecc 1148
1149//-----------------------------------------------------------------------
1150// part added by F. Noferini (some methods)
1151Bool_t AliFlowTrackCuts::PassesTOFbayesianCut(AliESDtrack* track){
3abf7ecc 1152
1153 Bool_t goodtrack = track && (track->GetStatus() & AliESDtrack::kTOFpid) && (track->GetTOFsignal() > 12000) && (track->GetTOFsignal() < 100000) && (track->GetIntegratedLength() > 365) && !(track->GetStatus() & AliESDtrack::kTOFmismatch);
1154
1155 if (! goodtrack)
1156 return kFALSE;
1157
1158 Int_t pdg = GetESDPdg(track,"bayesianALL");
1159 // printf("pdg set to %i\n",pdg);
1160
1161 Int_t pid = 0;
1162 Float_t prob = 0;
1163 switch (fAliPID)
1164 {
1165 case AliPID::kPion:
1166 pid=211;
1167 prob = fProbBayes[2];
1168 break;
1169 case AliPID::kKaon:
1170 pid=321;
1171 prob = fProbBayes[3];
1172 break;
1173 case AliPID::kProton:
1174 pid=2212;
1175 prob = fProbBayes[4];
1176 break;
1177 case AliPID::kElectron:
1178 pid=-11;
1179 prob = fProbBayes[0];
1180 break;
1181 default:
1182 return kFALSE;
1183 }
1184
1185 // 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);
1186 if(TMath::Abs(pdg) == TMath::Abs(pid) && prob > 0.8){
1187 if(!fCutCharge)
1188 return kTRUE;
1189 else if (fCutCharge && fCharge * track->GetSign() > 0)
1190 return kTRUE;
1191 }
1192 return kFALSE;
1193}
1194//-----------------------------------------------------------------------
1195Int_t AliFlowTrackCuts::GetESDPdg(AliESDtrack *track,Option_t *option,Int_t ipart,Float_t cPi,Float_t cKa,Float_t cPr){
1196 Int_t pdg = 0;
1197 Int_t pdgvalues[5] = {-11,-13,211,321,2212};
1198 Float_t mass[5] = {5.10998909999999971e-04,1.05658000000000002e-01,1.39570000000000000e-01,4.93676999999999977e-01,9.38271999999999995e-01};
1199
1200 if(strstr(option,"bayesianTOF")){ // Bayesian TOF PID
1201 Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1202 Double_t rcc=0.;
1203
1204 Float_t pt = track->Pt();
1205
1206 Int_t iptesd = 0;
1207 while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1208
1209 if(cPi < 0){
1210 c[0] = fC[iptesd][0];
1211 c[1] = fC[iptesd][1];
1212 c[2] = fC[iptesd][2];
1213 c[3] = fC[iptesd][3];
1214 c[4] = fC[iptesd][4];
1215 }
1216 else{
1217 c[0] = 0.0;
1218 c[1] = 0.0;
1219 c[2] = cPi;
1220 c[3] = cKa;
1221 c[4] = cPr;
1222 }
1223
1224 Double_t r1[10]; track->GetTOFpid(r1);
1225
1226 Int_t i;
1227 for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1228
1229 Double_t w[10];
1230 for (i=0; i<5; i++){
1231 w[i]=c[i]*r1[i]/rcc;
1232 fProbBayes[i] = w[i];
1233 }
1234 if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1235 pdg = 211*Int_t(track->GetSign());
1236 }
1237 else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1238 pdg = 2212*Int_t(track->GetSign());
1239 }
1240 else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1241 pdg = 321*Int_t(track->GetSign());
1242 }
1243 else if (w[0]>=w[1]) { //electrons
1244 pdg = -11*Int_t(track->GetSign());
1245 }
1246 else{ // muon
1247 pdg = -13*Int_t(track->GetSign());
1248 }
1249 }
1250
1251 else if(strstr(option,"bayesianTPC")){ // Bayesian TPC PID
1252 Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1253 Double_t rcc=0.;
1254
1255 Float_t pt = track->Pt();
1256
1257 Int_t iptesd = 0;
1258 while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1259
1260 if(cPi < 0){
1261 c[0] = fC[iptesd][0];
1262 c[1] = fC[iptesd][1];
1263 c[2] = fC[iptesd][2];
1264 c[3] = fC[iptesd][3];
1265 c[4] = fC[iptesd][4];
1266 }
1267 else{
1268 c[0] = 0.0;
1269 c[1] = 0.0;
1270 c[2] = cPi;
1271 c[3] = cKa;
1272 c[4] = cPr;
1273 }
1274
1275 Double_t r1[10]; track->GetTPCpid(r1);
1276
1277 Int_t i;
1278 for (i=0; i<5; i++) rcc+=(c[i]*r1[i]);
1279
1280 Double_t w[10];
1281 for (i=0; i<5; i++){
1282 w[i]=c[i]*r1[i]/rcc;
1283 fProbBayes[i] = w[i];
1284 }
1285 if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1286 pdg = 211*Int_t(track->GetSign());
1287 }
1288 else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1289 pdg = 2212*Int_t(track->GetSign());
1290 }
1291 else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1292 pdg = 321*Int_t(track->GetSign());
1293 }
1294 else if (w[0]>=w[1]) { //electrons
1295 pdg = -11*Int_t(track->GetSign());
1296 }
1297 else{ // muon
1298 pdg = -13*Int_t(track->GetSign());
1299 }
1300 }
1301
1302 else if(strstr(option,"bayesianALL")){
1303 Double_t c[5]={0.01, 0.01, 0.85, 0.1, 0.05};
1304 Double_t rcc=0.;
1305
1306 Float_t pt = track->Pt();
1307
1308 Int_t iptesd = 0;
1309 while(pt > fBinLimitPID[iptesd] && iptesd < fnPIDptBin-1) iptesd++;
1310
1311 if(cPi < 0){
1312 c[0] = fC[iptesd][0];
1313 c[1] = fC[iptesd][1];
1314 c[2] = fC[iptesd][2];
1315 c[3] = fC[iptesd][3];
1316 c[4] = fC[iptesd][4];
1317 }
1318 else{
1319 c[0] = 0.0;
1320 c[1] = 0.0;
1321 c[2] = cPi;
1322 c[3] = cKa;
1323 c[4] = cPr;
1324 }
1325
1326 Double_t r1[10]; track->GetTOFpid(r1);
1327 Double_t r2[10]; track->GetTPCpid(r2);
1328
1329 Int_t i;
1330 for (i=0; i<5; i++) rcc+=(c[i]*r1[i]*r2[i]);
1331
1332
1333 Double_t w[10];
1334 for (i=0; i<5; i++){
1335 w[i]=c[i]*r1[i]*r2[i]/rcc;
1336 fProbBayes[i] = w[i];
1337 }
1338
1339 if (w[2]>=w[3] && w[2]>=w[4] && w[2]>=w[1] && w[2]>=w[0]) {//pion
1340 pdg = 211*Int_t(track->GetSign());
1341 }
1342 else if (w[4]>=w[3] && w[4]>=w[1] && w[4]>=w[0]) {//proton
1343 pdg = 2212*Int_t(track->GetSign());
1344 }
1345 else if (w[3]>=w[1] && w[3]>=w[0]){//kaon
1346 pdg = 321*Int_t(track->GetSign());
1347 }
1348 else if (w[0]>=w[1]) { //electrons
1349 pdg = -11*Int_t(track->GetSign());
1350 }
1351 else{ // muon
1352 pdg = -13*Int_t(track->GetSign());
1353 }
1354 }
1355
1356 else if(strstr(option,"sigmacutTOF")){
1357 printf("PID not implemented yet: %s\nNO PID!!!!\n",option);
1358 Float_t p = track->P();
1359
1360 // Take expected times
1361 Double_t exptimes[5];
1362 track->GetIntegratedTimes(exptimes);
1363
1364 // Take resolution for TOF response
aab6527a 1365 // like fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
1366 Float_t resolution = fESDpid.GetTOFResponse().GetExpectedSigma(p, exptimes[ipart], mass[ipart]);
3abf7ecc 1367
1368 if(TMath::Abs(exptimes[ipart] - track->GetTOFsignal()) < 3 * resolution){
1369 pdg = pdgvalues[ipart] * Int_t(track->GetSign());
1370 }
1371 }
1372
1373 else{
1374 printf("Invalid PID option: %s\nNO PID!!!!\n",option);
1375 }
1376
1377 return pdg;
1378}
1379//-----------------------------------------------------------------------
1380void AliFlowTrackCuts::SetPriors(){
1381 // set abbundancies
1382 fBinLimitPID[0] = 0.30;
1383 fC[0][0] = 0.015;
1384 fC[0][1] = 0.015;
1385 fC[0][2] = 1;
1386 fC[0][3] = 0.0025;
1387 fC[0][4] = 0.000015;
1388 fBinLimitPID[1] = 0.35;
1389 fC[1][0] = 0.015;
1390 fC[1][1] = 0.015;
1391 fC[1][2] = 1;
1392 fC[1][3] = 0.01;
1393 fC[1][4] = 0.001;
1394 fBinLimitPID[2] = 0.40;
1395 fC[2][0] = 0.015;
1396 fC[2][1] = 0.015;
1397 fC[2][2] = 1;
1398 fC[2][3] = 0.026;
1399 fC[2][4] = 0.004;
1400 fBinLimitPID[3] = 0.45;
1401 fC[3][0] = 0.015;
1402 fC[3][1] = 0.015;
1403 fC[3][2] = 1;
1404 fC[3][3] = 0.026;
1405 fC[3][4] = 0.004;
1406 fBinLimitPID[4] = 0.50;
1407 fC[4][0] = 0.015;
1408 fC[4][1] = 0.015;
1409 fC[4][2] = 1.000000;
1410 fC[4][3] = 0.05;
1411 fC[4][4] = 0.01;
1412 fBinLimitPID[5] = 0.60;
1413 fC[5][0] = 0.012;
1414 fC[5][1] = 0.012;
1415 fC[5][2] = 1;
1416 fC[5][3] = 0.085;
1417 fC[5][4] = 0.022;
1418 fBinLimitPID[6] = 0.70;
1419 fC[6][0] = 0.01;
1420 fC[6][1] = 0.01;
1421 fC[6][2] = 1;
1422 fC[6][3] = 0.12;
1423 fC[6][4] = 0.036;
1424 fBinLimitPID[7] = 0.80;
1425 fC[7][0] = 0.0095;
1426 fC[7][1] = 0.0095;
1427 fC[7][2] = 1;
1428 fC[7][3] = 0.15;
1429 fC[7][4] = 0.05;
1430 fBinLimitPID[8] = 0.90;
1431 fC[8][0] = 0.0085;
1432 fC[8][1] = 0.0085;
1433 fC[8][2] = 1;
1434 fC[8][3] = 0.18;
1435 fC[8][4] = 0.074;
1436 fBinLimitPID[9] = 1;
1437 fC[9][0] = 0.008;
1438 fC[9][1] = 0.008;
1439 fC[9][2] = 1;
1440 fC[9][3] = 0.22;
1441 fC[9][4] = 0.1;
1442 fBinLimitPID[10] = 1.20;
1443 fC[10][0] = 0.007;
1444 fC[10][1] = 0.007;
1445 fC[10][2] = 1;
1446 fC[10][3] = 0.28;
1447 fC[10][4] = 0.16;
1448 fBinLimitPID[11] = 1.40;
1449 fC[11][0] = 0.0066;
1450 fC[11][1] = 0.0066;
1451 fC[11][2] = 1;
1452 fC[11][3] = 0.35;
1453 fC[11][4] = 0.23;
1454 fBinLimitPID[12] = 1.60;
1455 fC[12][0] = 0.0075;
1456 fC[12][1] = 0.0075;
1457 fC[12][2] = 1;
1458 fC[12][3] = 0.40;
1459 fC[12][4] = 0.31;
1460 fBinLimitPID[13] = 1.80;
1461 fC[13][0] = 0.0062;
1462 fC[13][1] = 0.0062;
1463 fC[13][2] = 1;
1464 fC[13][3] = 0.45;
1465 fC[13][4] = 0.39;
1466 fBinLimitPID[14] = 2.00;
1467 fC[14][0] = 0.005;
1468 fC[14][1] = 0.005;
1469 fC[14][2] = 1;
1470 fC[14][3] = 0.46;
1471 fC[14][4] = 0.47;
1472 fBinLimitPID[15] = 2.20;
1473 fC[15][0] = 0.0042;
1474 fC[15][1] = 0.0042;
1475 fC[15][2] = 1;
1476 fC[15][3] = 0.5;
1477 fC[15][4] = 0.55;
1478 fBinLimitPID[16] = 2.40;
1479 fC[16][0] = 0.007;
1480 fC[16][1] = 0.007;
1481 fC[16][2] = 1;
1482 fC[16][3] = 0.5;
1483 fC[16][4] = 0.6;
1484
1485 for(Int_t i=17;i<fnPIDptBin;i++){
1486 fBinLimitPID[i] = 2.0 + 0.2 * (i-14);
1487 fC[i][0] = fC[13][0];
1488 fC[i][1] = fC[13][1];
1489 fC[i][2] = fC[13][2];
1490 fC[i][3] = fC[13][3];
1491 fC[i][4] = fC[13][4];
1492 }
1493}
1494// end part added by F. Noferini
1495//-----------------------------------------------------------------------