AliPhysics  master (3d17d9d)
AliAnalysisTaskSEHFSystPID.cxx
Go to the documentation of this file.
1 /* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. */
2 
4 // \class AliAnalysisTaskSEHFSystPID //
5 // \brief analysis task for the study of PID systematic uncertainties of HF particles //
6 // \author: A. M. Barbano, anastasia.maria.barbano@cern.ch //
7 // \author: F. Grosa, fabrizio.grosa@cern.ch //
9 
10 #include <cstddef>
11 #include <TVector3.h>
12 #include <TLorentzVector.h>
13 #include <TDatabasePDG.h>
14 #include <TChain.h>
15 
17 #include "AliAODHandler.h"
18 #include "AliInputEventHandler.h"
19 #include "AliAODMCParticle.h"
20 #include "AliAnalysisManager.h"
21 #include "AliMultSelection.h"
22 #include "AliESDtrack.h"
23 #include "AliRDHFCuts.h"
24 
26 
27 //________________________________________________________________________
29 AliAnalysisTaskSE("TaskNsigmaPID"),
30 fOutputList(nullptr),
31 fHistNEvents(nullptr),
32 fHistQtVsMassKinks(nullptr),
33 fHistPDaughterVsMotherKink(nullptr),
34 fHistdEdxVsPMotherKink(nullptr),
35 fHistOpeningAngleVsPMotherKink(nullptr),
36 fHistNTPCclsVsRadius(nullptr),
37 fPIDtree(nullptr),
38 fPTPC(0),
39 fPTOF(0),
40 fPHMPID(0),
41 fdEdxTPC(0),
42 fdEdxITS(0),
43 fToF(0),
44 fPt(0),
45 fTPCNcls(0),
46 fTPCNclsPID(0),
47 fTrackLength(0),
48 fStartTimeRes(0),
49 fTPCNcrossed(0),
50 fTPCFindable(0),
51 fITSclsMap(0),
52 fHMPIDsignal(0),
53 fHMPIDoccupancy(0),
54 fTrackInfoMap(0),
55 fEta(-9999),
56 fPhi(9999),
57 fPDGcode(-1),
58 fTag(0),
59 fNsigmaMaxForKaonTag(0.2),
60 fNsigmaMaxForNucleiTag(3.),
61 fQtMinKinks(0.15),
62 fRMinKinks(120),
63 fRMaxKinks(210),
64 fDeadZoneWidth(3.),
65 fCutGeoNcrNclLength(130.),
66 fCutGeoNcrNclGeom1Pt(1.5),
67 fCutGeoNcrNclFractionNcr(0.85),
68 fCutGeoNcrNclFractionNcl(0.7),
69 fCentMin(0.),
70 fCentMax(100.),
71 fCentEstimator(kCentOff),
72 fTriggerClass(""),
73 fTriggerMask(AliVEvent::kINT7),
74 fIsMC(false),
75 fSystem(0),
76 fConversionFactordEdx(100.),
77 fESDtrackCuts(nullptr),
78 fAOD(nullptr),
79 fPIDresp(nullptr),
80 fV0cuts(nullptr),
81 fFillTreeWithPIDInfo(true),
82 fFillTreeWithNsigmaPIDOnly(false),
83 fFillTreeWithRawPIDOnly(false),
84 fFillTreeWithTrackQualityInfo(false),
85 fEnabledDownSampling(false),
86 fFracToKeepDownSampling(0.1),
87 fPtMaxDownSampling(1.5),
88 fDownSamplingOpt(0),
89 fAODProtection(1),
90 fRunNumberPrevEvent(-1),
91 fEnableNsigmaTPCDataCorr(false),
92 fSystNsigmaTPCDataCorr(AliAODPidHF::kNone),
93 fMeanNsigmaTPCPionData{},
94 fMeanNsigmaTPCKaonData{},
95 fMeanNsigmaTPCProtonData{},
96 fSigmaNsigmaTPCPionData{},
97 fSigmaNsigmaTPCKaonData{},
98 fSigmaNsigmaTPCProtonData{},
99 fPlimitsNsigmaTPCDataCorr{},
100 fNPbinsNsigmaTPCDataCorr(0),
101 fEtalimitsNsigmaTPCDataCorr{},
102 fNEtabinsNsigmaTPCDataCorr(0),
103 fUseAliEventCuts(false),
104 fAliEventCuts(),
105 fApplyPbPbOutOfBunchPileupCuts(),
106 fUseTimeRangeCutForPbPb2018(true),
107 fTimeRangeCut()
108 {
109  //
110  // default constructur
111  //
112 
113  for(int iHisto=0; iHisto<5; iHisto++) fHistArmenteroPlot[iHisto] = nullptr;
114 
115  for(int iDet=0; iDet<kNMaxDet; iDet++) {
116  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++) {
117  fPIDNsigma[iDet][iHypo] = numeric_limits<short>::min();
118  fHistNsigmaVsPt[iDet][iHypo] = nullptr;
119  }
120  }
121 
122  fEnabledSpecies[kPion] = true;
123  fEnabledSpecies[kKaon] = true;
124  fEnabledSpecies[kProton] = true;
125  fEnabledSpecies[kElectron] = false;
126  fEnabledSpecies[kDeuteron] = false;
127  fEnabledSpecies[kTriton] = false;
128  fEnabledSpecies[kHe3] = false;
129 
130  fEnabledDet[kITS] = false;
131  fEnabledDet[kTPC] = true;
132  fEnabledDet[kTOF] = true;
133  fEnabledDet[kHMPID] = false;
134 
135  for(int iP=0; iP<=AliAODPidHF::kMaxPBins; iP++) {
136  fPlimitsNsigmaTPCDataCorr[iP] = 0.;
137  }
138  for(int iEta=0; iEta<=AliAODPidHF::kMaxEtaBins; iEta++) {
139  fEtalimitsNsigmaTPCDataCorr[iEta] = 0.;
140  }
141 
142  fAliEventCuts.SetManualMode();
143 }
144 
145 //________________________________________________________________________
147 AliAnalysisTaskSE(name),
148 fOutputList(nullptr),
149 fHistNEvents(nullptr),
150 fHistQtVsMassKinks(nullptr),
151 fHistPDaughterVsMotherKink(nullptr),
152 fHistdEdxVsPMotherKink(nullptr),
153 fHistOpeningAngleVsPMotherKink(nullptr),
154 fHistNTPCclsVsRadius(nullptr),
155 fPIDtree(nullptr),
156 fPTPC(0),
157 fPTOF(0),
158 fPHMPID(0),
159 fdEdxTPC(0),
160 fdEdxITS(0),
161 fToF(0),
162 fPt(0),
163 fTPCNcls(0),
164 fTPCNclsPID(0),
165 fTrackLength(0),
166 fStartTimeRes(0),
167 fTPCNcrossed(0),
168 fTPCFindable(0),
169 fITSclsMap(0),
170 fHMPIDsignal(0),
171 fHMPIDoccupancy(0),
172 fTrackInfoMap(0),
173 fEta(-9999),
174 fPhi(9999),
175 fPDGcode(-1),
176 fTag(0),
177 fNsigmaMaxForKaonTag(0.2),
178 fNsigmaMaxForNucleiTag(3.),
179 fQtMinKinks(0.15),
180 fRMinKinks(120),
181 fRMaxKinks(210),
182 fDeadZoneWidth(3.),
183 fCutGeoNcrNclLength(130.),
184 fCutGeoNcrNclGeom1Pt(1.5),
185 fCutGeoNcrNclFractionNcr(0.85),
186 fCutGeoNcrNclFractionNcl(0.7),
187 fCentMin(0.),
188 fCentMax(100.),
189 fCentEstimator(kCentOff),
190 fTriggerClass(""),
191 fTriggerMask(AliVEvent::kINT7),
192 fIsMC(false),
193 fSystem(system),
194 fConversionFactordEdx(100.),
195 fESDtrackCuts(nullptr),
196 fAOD(nullptr),
197 fPIDresp(nullptr),
198 fV0cuts(nullptr),
199 fFillTreeWithPIDInfo(true),
200 fFillTreeWithNsigmaPIDOnly(false),
201 fFillTreeWithRawPIDOnly(false),
202 fFillTreeWithTrackQualityInfo(false),
203 fEnabledDownSampling(false),
204 fFracToKeepDownSampling(0.1),
205 fPtMaxDownSampling(1.5),
206 fDownSamplingOpt(0),
207 fAODProtection(1),
208 fRunNumberPrevEvent(-1),
209 fEnableNsigmaTPCDataCorr(false),
210 fSystNsigmaTPCDataCorr(AliAODPidHF::kNone),
211 fMeanNsigmaTPCPionData{},
221 fUseAliEventCuts(false),
222 fAliEventCuts(),
226 {
227  //
228  // standard constructur
229  //
230 
231  for(int iHisto=0; iHisto<5; iHisto++) fHistArmenteroPlot[iHisto] = nullptr;
232 
233  for(int iDet=0; iDet<kNMaxDet; iDet++) {
234  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++) {
235  fPIDNsigma[iDet][iHypo] = numeric_limits<short>::min();
236  fHistNsigmaVsPt[iDet][iHypo] = nullptr;
237  }
238  }
239 
240  fEnabledSpecies[kPion] = true;
241  fEnabledSpecies[kKaon] = true;
242  fEnabledSpecies[kProton] = true;
243  fEnabledSpecies[kElectron] = false;
244  fEnabledSpecies[kDeuteron] = false;
245  fEnabledSpecies[kTriton] = false;
246  fEnabledSpecies[kHe3] = false;
247 
248  fEnabledDet[kITS] = false;
249  fEnabledDet[kTPC] = true;
250  fEnabledDet[kTOF] = true;
251  fEnabledDet[kHMPID] = false;
252 
253  for(int iP=0; iP<=AliAODPidHF::kMaxPBins; iP++) {
254  fPlimitsNsigmaTPCDataCorr[iP] = 0.;
255  }
256  for(int iEta=0; iEta<=AliAODPidHF::kMaxEtaBins; iEta++) {
257  fEtalimitsNsigmaTPCDataCorr[iEta] = 0.;
258  }
259 
260  fAliEventCuts.SetManualMode();
261 
262  DefineInput(0, TChain::Class());
263  DefineOutput(1, TList::Class());
264  DefineOutput(2, TTree::Class());
265 }
266 
267 //________________________________________________________________________
269 {
270  // Destructor
271  if (fOutputList) {
272  delete fHistNEvents;
273  for(int iHisto=0; iHisto<5; iHisto++) delete fHistArmenteroPlot[iHisto];
274  delete fHistQtVsMassKinks;
277  delete fHistdEdxVsPMotherKink;
278  delete fHistNTPCclsVsRadius;
279  delete fOutputList;
280  }
281 
282  if(fPIDtree) delete fPIDtree;
283  if(fESDtrackCuts) delete fESDtrackCuts;
284  if(fV0cuts) delete fV0cuts;
285 }
286 
287 //________________________________________________________________________
289 {
290  // create track cuts
291  if(!fESDtrackCuts) {
292  fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
293  fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE,0);
294  fESDtrackCuts->SetEtaRange(-0.8, 0.8);
295  fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
296  }
297 
298  // V0 Kine cuts
299  fV0cuts = new AliAODv0KineCuts;
300  fV0cuts->SetNoKinks(false);
301 
302  fOutputList = new TList();
303  fOutputList->SetOwner(true);
304 
305  fHistNEvents = new TH1F("fHistNEvents","Number of processed events;;Number of events",13,-1.5,11.5);
306  fHistNEvents->Sumw2();
307  fHistNEvents->SetMinimum(0);
308  fHistNEvents->GetXaxis()->SetBinLabel(1,"Read from AOD");
309  fHistNEvents->GetXaxis()->SetBinLabel(2,"Mismatch AOD");
310  fHistNEvents->GetXaxis()->SetBinLabel(3,"Pass Phys. Sel. + Trig");
311  fHistNEvents->GetXaxis()->SetBinLabel(4,"No vertex");
312  fHistNEvents->GetXaxis()->SetBinLabel(5,"Vertex contributors < 1");
313  fHistNEvents->GetXaxis()->SetBinLabel(6,"Without SPD vertex");
314  fHistNEvents->GetXaxis()->SetBinLabel(7,"Error on zVertex>0.5");
315  fHistNEvents->GetXaxis()->SetBinLabel(8,"|zVertex|>10");
316  fHistNEvents->GetXaxis()->SetBinLabel(9,"Good Z vertex");
317  fHistNEvents->GetXaxis()->SetBinLabel(10,"Cent corr cuts");
318  fHistNEvents->GetXaxis()->SetBinLabel(11,"V0mult vs. nTPC cls");
319  fHistNEvents->GetXaxis()->SetBinLabel(12,"excluded time range");
320  fHistNEvents->GetXaxis()->SetBinLabel(13,"Selected events");
322 
323  TString armenteronames[5] = {"All","K0s","Lambda","AntiLambda","Gamma"};
324  TString armenterolabels[5] = {"all","K_{s}^{0}","#Lambda","#bar{#Lambda}","#gamma"};
325  for(int iHisto=0; iHisto<5; iHisto++) {
326  fHistArmenteroPlot[iHisto] = new TH2F(Form("fHistArmenteroPlot%s",armenteronames[iHisto].Data()),Form("Armenteros-Podolanski plot - %s;#it{#alpha};#it{q}_{T} (GeV/#it{c})",armenterolabels[iHisto].Data()),200,-1.,1.,200,0.,0.4);
327  fOutputList->Add(fHistArmenteroPlot[iHisto]);
328  }
329 
330  fHistQtVsMassKinks = new TH2F("fHistQtVsMassKinks","#it{q}_{T} vs. #it{M}(#mu#nu) kink mother;#it{M}(#mu#nu) (GeV/#it{c}^{2});#it{q}_{T} (GeV/#it{c})",250,0.1,0.6,150,0.,0.3);
332 
333  fHistPDaughterVsMotherKink = new TH2F("fHistPDaughterVsMotherKink","#it{p} (daughter) vs. #it{p} (mother) kinks;#it{p}^{mother} (GeV/#it{c});#it{p}^{dau} (GeV/#it{c})",100,0,20,100,0,20);
335 
336  fHistOpeningAngleVsPMotherKink = new TH2F("fHistOpeningAngleVsPMotherKink","#it{#alpha} vs. #it{p} (mother) kinks;#it{p}^{mother} (GeV/#it{c});#it{#alpha} (deg)",100,0,20,90,0,90);
338 
339  fHistdEdxVsPMotherKink = new TH2F("fHistdEdxVsPMotherKink","dE/dx vs. #it{p} (mother) kinks;#it{p}^{mother} (GeV/#it{c});d#it{E}/d#it{x} (a.u.)",100,0,20,125,0,250);
341 
342  fHistNTPCclsVsRadius = new TH2F("fHistNTPCclsVsRadius","N TPC clusters (mother) vs. #it{R} kinks;#it{R} (cm);N TPC clusters (mother)",50,0,250,160,-0.5,159.5);
344 
345  TString detnames[kNMaxDet] = {"ITS","TPC","TOF","HMPID"};
346  TString partnameshort[kNMaxHypo] = {"pi","K","p","e","d","t","He3"};
347  TString hyponames[kNMaxHypo] = {"Pion","Kaon","Proton","Electron","Deuteron","Triton","He3"};
348 
349  if(fIsMC) {
350  for(int iDet=0; iDet<kNMaxDet; iDet++) {
351  if(!fEnabledDet[iDet])
352  continue;
353  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++) {
354  if(!fEnabledSpecies[iHypo])
355  continue;
356  fHistNsigmaVsPt[iDet][iHypo] = new TH2F(Form("fHistNsigma%svsPt_%s",detnames[iDet].Data(),hyponames[iHypo].Data()),Form(";#it{p}_{T} (GeV/#it{c});N_{#sigma}^{%s}(%s)",detnames[iDet].Data(),hyponames[iHypo].Data()),500,0,50,1000,-50,50);
357  fOutputList->Add(fHistNsigmaVsPt[iDet][iHypo]);
358  }
359  }
360  }
361 
362  fPIDtree = new TTree("fPIDtree","fPIDtree");
363  fPIDtree->Branch("pT",&fPt,"pT/s");
364  fPIDtree->Branch("eta",&fEta,"eta/S");
365  fPIDtree->Branch("phi",&fPhi,"phi/s");
367  if(fEnabledDet[kITS])
368  fPIDtree->Branch("p",&fP,"p/s");
369  if(fEnabledDet[kTPC])
370  fPIDtree->Branch("pTPC",&fPTPC,"pTPC/s");
371  if(fEnabledDet[kTOF])
372  fPIDtree->Branch("pTOF",&fPTOF,"pTOF/s");
373  if(fEnabledDet[kHMPID])
374  fPIDtree->Branch("pHMPID",&fPHMPID,"pHMPID/s");
376  for(int iDet=0; iDet<kNMaxDet; iDet++) {
377  if(!fEnabledDet[iDet])
378  continue;
379  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++) {
380  if(!fEnabledSpecies[iHypo])
381  continue;
382  fPIDtree->Branch(Form("n_sigma_%s_%s",detnames[iDet].Data(),partnameshort[iHypo].Data()),&fPIDNsigma[iDet][iHypo],Form("n_sigma_%s_%s/S",detnames[iDet].Data(),partnameshort[iHypo].Data()));
383  }
384  }
385  }
387  if(fEnabledDet[kITS]) {
388  fPIDtree->Branch("dEdxITS",&fdEdxITS,"dEdxITS/s");
389  fPIDtree->Branch("NclusterPIDTPC",&fTPCNclsPID,"NclusterPIDTPC/b");
390  }
391  if(fEnabledDet[kTPC]) {
392  fPIDtree->Branch("dEdxTPC",&fdEdxTPC,"dEdxTPC/s");
393  fPIDtree->Branch("ITSclsMap",&fITSclsMap,"ITSclsMap/b");
394  }
395  if(fEnabledDet[kTOF]) {
396  fPIDtree->Branch("ToF",&fToF,"ToF/s");
397  fPIDtree->Branch("TrackLength",&fTrackLength,"TrackLength/s");
398  fPIDtree->Branch("StartTimeRes",&fStartTimeRes,"StartTimeRes/s");
399  }
400  if(fEnabledDet[kHMPID]) {
401  fPIDtree->Branch("HMPIDsig",&fHMPIDsignal,"HMPIDsig/s");
402  fPIDtree->Branch("HMPIDocc",&fHMPIDoccupancy,"HMPIDocc/s");
403  }
404  }
405  }
407  fPIDtree->Branch("NclusterTPC",&fTPCNcls,"NclusterTPC/b");
408  fPIDtree->Branch("NcrossedRowsTPC",&fTPCNcrossed,"NcrossedRowsTPC/b");
409  fPIDtree->Branch("NFindableTPC",&fTPCFindable,"NFindableClustersTPC/b");
410  }
411  fPIDtree->Branch("trackbits",&fTrackInfoMap,"trackbits/b"); // basic track info always filled
412  fPIDtree->Branch("tag",&fTag,"tag/s");
413  if(fIsMC) fPIDtree->Branch("PDGcode",&fPDGcode,"PDGcode/I");
414 
415  if(fUseAliEventCuts) { //add QA plots if event cuts used
416  fAliEventCuts.AddQAplotsToList(fOutputList,true);
417  }
418 
419  // post data
420  PostData(1, fOutputList);
421  PostData(2, fPIDtree);
422 }
423 
424 //________________________________________________________________________
426 {
427  // main event loop
428  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
429  if (!fAOD && AODEvent() && IsStandardAOD()) {
430  // In case there is an AOD handler writing a standard AOD, use the AOD
431  // event in memory rather than the input (ESD) event.
432  fAOD = dynamic_cast<AliAODEvent*> (AODEvent());
433  }
434  if (!fAOD) {
435  AliWarning("AliAnalysisTaskSEHFSystPID::Exec(): bad AOD");
436  PostData(1, fOutputList);
437  return;
438  }
439  fHistNEvents->Fill(-1);
440 
441  if(fAODProtection>=0){
442  // Protection against different number of events in the AOD and deltaAOD
443  // In case of discrepancy the event is rejected.
444  int matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
445  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
446  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
447  fHistNEvents->Fill(0);
448  PostData(1, fOutputList);
449  return;
450  }
451  }
452 
453  if(fUseAliEventCuts)
454  fAliEventCuts.AcceptEvent(fAOD); //for QA plots
455 
456  if(TMath::Abs(fAOD->GetMagneticField())<0.001) return;
457 
458  AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
459  if(!aodHandler) {
460  AliWarning("AliAnalysisTaskSEHFSystPID::Exec(): No AliInputEventHandler!");
461  return;
462  }
463 
464  unsigned int maskPhysSel = aodHandler->IsEventSelected();
465  TString firedTriggerClasses = fAOD->GetFiredTriggerClasses();
466  if(!fIsMC && (fAOD->GetRunNumber()<136851 || fAOD->GetRunNumber()>139517)) {
467  if(!(firedTriggerClasses.Contains(fTriggerClass.Data()))) return;
468  }
469  if((maskPhysSel & fTriggerMask)==0) return;
470 
471  if(!IsCentralitySelected()) {
472  PostData(1, fOutputList);
473  return;
474  }
475 
476  fHistNEvents->Fill(1);
477 
478  if (!IsVertexAccepted()) {
479  fHistNEvents->Fill(2);
480  PostData(1, fOutputList);
481  return;
482  }
483 
484  if(fUseAliEventCuts) {
486  if(sel>0) {
487  fHistNEvents->Fill(sel);
488  PostData(1, fOutputList);
489  return;
490  }
491  }
492  else {
494  if(fAOD->GetRunNumber() != fRunNumberPrevEvent){
495  fTimeRangeCut.InitFromRunNumber(fAOD->GetRunNumber());
496  }
497  if(fTimeRangeCut.CutEvent(fAOD)){
498  fHistNEvents->Fill(10);
499  PostData(1, fOutputList);
500  return;
501  }
502  }
503  }
504 
505  fHistNEvents->Fill(11);
506  // load MC particles
507  TClonesArray *arrayMC=0;
508  if(fIsMC){
509  arrayMC = (TClonesArray*)fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
510  if(!arrayMC) {
511  AliWarning("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
512  return;
513  }
514  }
515 
516  //get pid response
517  if(!fPIDresp) fPIDresp = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
518 
519  //load data correction for NsigmaTPC, if enabled
521  if(fAOD->GetRunNumber()!=fRunNumberPrevEvent) {
523  }
524  }
525 
526  // V0 selection
527  if(fSystem==0) fV0cuts->SetMode(AliAODv0KineCuts::kPurity,AliAODv0KineCuts::kPP);
528  else if(fSystem==1) fV0cuts->SetMode(AliAODv0KineCuts::kPurity,AliAODv0KineCuts::kPbPb);
529  fV0cuts->SetEvent(fAOD);
530 
531  vector<short> idPionFromK0s;
532  vector<short> idPionFromL;
533  vector<short> idProtonFromL;
534  vector<short> idElectronFromGamma;
535  vector<short> idKaonFromKinks;
536  GetTaggedV0s(idPionFromK0s,idPionFromL,idProtonFromL,idElectronFromGamma);
537  GetTaggedKaonsFromKinks(idKaonFromKinks);
538 
539  vector<short>::iterator it;
540 
541  const int nTracks = fAOD->GetNumberOfTracks();
542  //loop on tracks
543 
544  for(int iTrack=0; iTrack<nTracks; iTrack++) {
545  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
546  if(!track) continue;
547  //applying ESDtrackCut
548  if(!fESDtrackCuts->IsSelected(track)) continue;
549 
550  if(fEnabledDownSampling && fFracToKeepDownSampling<1. && track->Pt()<fPtMaxDownSampling) {
551  double pseudoRand = track->Pt()*1000.-(long)(track->Pt()*1000);
552  if(fDownSamplingOpt==0 && pseudoRand>fFracToKeepDownSampling) continue; //keep tracks with pseudorand < fFracToKeepDownSampling
553  else if(fDownSamplingOpt==1 && pseudoRand<(1-fFracToKeepDownSampling)) continue; //keep tracks with pseudorand > 1-fFracToKeepDownSampling
554  }
555 
556  fPt = ConvertFloatToUnsignedShort(track->Pt()*1000);
557  fP = ConvertFloatToUnsignedShort(track->P()*1000);
558  fPTPC = ConvertFloatToUnsignedShort(track->GetTPCmomentum()*1000);
560  double momHMPID[3];
561  bool momHMPIDok = track->GetOuterHmpPxPyPz(momHMPID);
562  if(momHMPIDok)
563  fPHMPID = ConvertFloatToUnsignedShort(TMath::Sqrt(momHMPID[0]*momHMPID[0]+momHMPID[1]*momHMPID[1]+momHMPID[2]*momHMPID[2])*1000);
564  else
565  fPHMPID = numeric_limits<short>::min();
566 
567  fEta = ConvertFloatToShort(track->Eta()*1000);
568  fPhi = ConvertFloatToUnsignedShort(track->Phi()*1000);
569 
570  //charge
571  if(track->Charge()>0) {
572  fTag |= kPositiveTrack;
573  fTag &= ~kNegativeTrack;
574  }
575  else if(track->Charge()<0) {
576  fTag |= kNegativeTrack;
577  fTag &= ~kPositiveTrack;
578  }
579 
580  if(track->HasPointOnITSLayer(0) || track->HasPointOnITSLayer(1))
582  else
584  if(track->HasPointOnITSLayer(1))
586  else
588  if(track->GetStatus() & AliESDtrack::kITSrefit)
590  else
592  if(track->GetStatus() & AliESDtrack::kTPCrefit)
594  else
596  if(IsSelectedByGeometricalCut(track))
598  else
600 
602  fTPCNcls = static_cast<unsigned char>(track->GetTPCNcls());
603  fTPCNcrossed = static_cast<unsigned char>(track->GetTPCNCrossedRows());
604  fTPCFindable = static_cast<unsigned char>(track->GetTPCNclsF());
605  }
606 
607  //PID
608  bool isDetOk[kNMaxDet] = {false,false,false};
609  if (fPIDresp->CheckPIDStatus(AliPIDResponse::kITS,track) == AliPIDResponse::kDetPidOk) {
611  isDetOk[kITS] = true;
612  }
613  else {
615  }
616  if (fPIDresp->CheckPIDStatus(AliPIDResponse::kTPC,track) == AliPIDResponse::kDetPidOk) {
618  isDetOk[kTPC] = true;
619  }
620  else {
622  }
623  if (fPIDresp->CheckPIDStatus(AliPIDResponse::kTOF,track) == AliPIDResponse::kDetPidOk) {
625  isDetOk[kTOF] = true;
626  }
627  else {
629  }
630  if (fPIDresp->CheckPIDStatus(AliPIDResponse::kHMPID,track) == AliPIDResponse::kDetPidOk) {
631  isDetOk[kHMPID] = true;
632  }
633 
635  if(!fFillTreeWithNsigmaPIDOnly) { //raw variables
636  //ITS variables
638  fITSclsMap = track->GetITSClusterMap();
639 
640  //TPC variables
642  fTPCNclsPID = static_cast<unsigned char>(track->GetTPCsignalN());
643 
644  //TOF variables
645  fTrackLength = ConvertFloatToUnsignedShort(track->GetIntegratedLength()*10);
646  fStartTimeRes = ConvertFloatToUnsignedShort(fPIDresp->GetTOFResponse().GetStartTimeRes(track->P())*100);
647 
648  if (!(track->GetStatus() & AliESDtrack::kTOFout) || !(track->GetStatus() & AliESDtrack::kTIME)) {
649  fToF = 0;
650  }
651  else {
652  if (fTrackLength < 3500) fToF = 0;
653  else {
654  float tof = track->GetTOFsignal();
655  float time0 = fPIDresp->GetTOFResponse().GetStartTime(track->P());
656  fToF = ConvertFloatToUnsignedShort((tof-time0)/10);
657  }
658  }
659 
660  //HMPID variables
661  fHMPIDsignal = ConvertFloatToUnsignedShort(track->GetHMPIDsignal()*100);
662  fHMPIDoccupancy = ConvertFloatToUnsignedShort(track->GetHMPIDoccupancy()*100);
663  }
664  if(!fFillTreeWithRawPIDOnly) { // nsigma
665  for(int iDet=0; iDet<kNMaxDet; iDet++) {
666  if(isDetOk[iDet]) {
667  FillNsigma(iDet, track);
668  }
669  else {
670  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++) {
671  fPIDNsigma[iDet][iHypo] = numeric_limits<short>::min();
672  }
673  }
674  }
675  }
676  }
677 
678  short trackid = track->GetID();
679 
680  bool filltree = false;
681  it = find(idPionFromK0s.begin(),idPionFromK0s.end(),trackid);
682  if(it!=idPionFromK0s.end()) {
683  fTag |= kIsPionFromK0s;
684  filltree = true;
685  }
686  else
687  fTag &= ~kIsPionFromK0s;
688 
689  it = find(idPionFromL.begin(),idPionFromL.end(),trackid);
690  if(it!=idPionFromL.end()) {
691  fTag |= kIsPionFromL;
692  filltree = true;
693  }
694  else
695  fTag &= ~kIsPionFromL;
696 
697  it = find(idProtonFromL.begin(),idProtonFromL.end(),trackid);
698  if(it!=idProtonFromL.end()) {
699  filltree = true;
700  fTag |= kIsProtonFromL;
701  }
702  else
703  fTag &= ~kIsProtonFromL;
704 
705  it = find(idElectronFromGamma.begin(),idElectronFromGamma.end(),trackid);
706  if(it!=idElectronFromGamma.end()) {
707  filltree = true;
709  }
710  else
712 
713  it = find(idKaonFromKinks.begin(),idKaonFromKinks.end(),trackid);
714  if(it!=idKaonFromKinks.end()) {
715  filltree = true;
717  }
718  else
720 
721  if(TMath::Abs(fPIDresp->NumberOfSigmasTOF(track,AliPID::kKaon))<fNsigmaMaxForKaonTag && isDetOk[kTOF]) {
722  filltree = true;
723  fTag |= kIsKaonFromTOF;
724  }
725  else
726  fTag &= ~kIsKaonFromTOF;
727 
728  if(TMath::Abs(fPIDresp->NumberOfSigmasTPC(track,AliPID::kKaon))<fNsigmaMaxForKaonTag && isDetOk[kTPC]) {
729  filltree = true;
730  fTag |= kIsKaonFromTPC;
731  }
732  else
733  fTag &= ~kIsKaonFromTPC;
734 
735  if(TMath::Abs(fPIDresp->NumberOfSigmasHMPID(track,AliPID::kKaon))<fNsigmaMaxForKaonTag && isDetOk[kHMPID]) {
736  filltree = true;
738  }
739  else
741 
742  if(isDetOk[kTPC] && isDetOk[kTOF] && TMath::Abs(fPIDresp->NumberOfSigmasTPC(track,AliPID::kDeuteron))<fNsigmaMaxForNucleiTag && TMath::Abs(fPIDresp->NumberOfSigmasTOF(track,AliPID::kDeuteron))<fNsigmaMaxForNucleiTag) {
743  filltree = true;
745  }
746  else
748 
749  if(isDetOk[kTPC] && isDetOk[kTOF] && TMath::Abs(fPIDresp->NumberOfSigmasTPC(track,AliPID::kTriton))<fNsigmaMaxForNucleiTag && TMath::Abs(fPIDresp->NumberOfSigmasTOF(track,AliPID::kTriton))<fNsigmaMaxForNucleiTag) {
750  filltree = true;
752  }
753  else
755 
756  if(isDetOk[kTPC] && isDetOk[kTOF] && TMath::Abs(fPIDresp->NumberOfSigmasTPC(track,AliPID::kHe3))<fNsigmaMaxForNucleiTag && TMath::Abs(fPIDresp->NumberOfSigmasTOF(track,AliPID::kHe3))<fNsigmaMaxForNucleiTag) {
757  filltree = true;
759  }
760  else
762 
763  if(fIsMC) {
764  fPDGcode = GetPDGcodeFromMC(track,arrayMC);
765  AliPIDResponse::EDetector det[4] = {AliPIDResponse::kITS,AliPIDResponse::kTPC,AliPIDResponse::kTOF,AliPIDResponse::kHMPID};
766  for(int iDet=0; iDet<kNMaxDet; iDet++) {
767  if(!fEnabledDet[iDet])
768  continue;
769 
770  switch(fPDGcode) {
771  case 211:
773  fHistNsigmaVsPt[iDet][kPion]->Fill(track->Pt(),fPIDresp->NumberOfSigmas(det[iDet],track,AliPID::kPion));
774  break;
775  case 321:
777  fHistNsigmaVsPt[iDet][kKaon]->Fill(track->Pt(),fPIDresp->NumberOfSigmas(det[iDet],track,AliPID::kKaon));
778  break;
779  case 2212:
781  fHistNsigmaVsPt[iDet][kProton]->Fill(track->Pt(),fPIDresp->NumberOfSigmas(det[iDet],track,AliPID::kProton));
782  break;
783  case 11:
785  fHistNsigmaVsPt[iDet][kElectron]->Fill(track->Pt(),fPIDresp->NumberOfSigmas(det[iDet],track,AliPID::kElectron));
786  break;
787  case 1000010020:
789  fHistNsigmaVsPt[iDet][kDeuteron]->Fill(track->Pt(),fPIDresp->NumberOfSigmas(det[iDet],track,AliPID::kDeuteron));
790  break;
791  case 1000010030:
793  fHistNsigmaVsPt[iDet][kTriton]->Fill(track->Pt(),fPIDresp->NumberOfSigmas(det[iDet],track,AliPID::kTriton));
794  break;
795  case 1000020030:
796  if(fEnabledSpecies[kHe3])
797  fHistNsigmaVsPt[iDet][kHe3]->Fill(track->Pt(),fPIDresp->NumberOfSigmas(det[iDet],track,AliPID::kHe3));
798  break;
799  }
800  }
801  }
802  else fPDGcode = 0;
803 
804  if(filltree && ((fIsMC && fPDGcode>=0) || !fIsMC)) fPIDtree->Fill();
805 
806  fTag = 0;
807  fTrackInfoMap = 0;
808 
809  // Post output data
810  PostData(2, fPIDtree);
811  }
812 
813  //clear vectors of ids
814  idPionFromK0s.clear();
815  idPionFromL.clear();
816  idProtonFromL.clear();
817  idElectronFromGamma.clear();
818  idKaonFromKinks.clear();
819 
820  fRunNumberPrevEvent = fAOD->GetRunNumber();
821 
822  // Post output data
823  PostData(1, fOutputList);
824 }
825 
826 //________________________________________________________________________
827 void AliAnalysisTaskSEHFSystPID::EnableParticleSpecies(bool pi, bool kao, bool pr, bool el, bool deu, bool tr, bool He3)
828 {
829  // function to enable/disable different particle species
830  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++)
831  fEnabledSpecies[iHypo] = false;
832 
833  if(pi) fEnabledSpecies[kPion] = true;
834  if(kao) fEnabledSpecies[kKaon] = true;
835  if(pr) fEnabledSpecies[kProton] = true;
836  if(el) fEnabledSpecies[kElectron] = true;
837  if(deu) fEnabledSpecies[kDeuteron] = true;
838  if(tr) fEnabledSpecies[kTriton] = true;
839  if(He3) fEnabledSpecies[kHe3] = true;
840 }
841 
842 //________________________________________________________________________
843 void AliAnalysisTaskSEHFSystPID::EnableDetectors(bool ITS, bool TPC, bool TOF, bool HMPID)
844 {
845  // function to enable/disable different PID detectors
846  for(int iDet=0; iDet<kNMaxDet; iDet++)
847  fEnabledDet[iDet] = false;
848 
849  if(ITS) fEnabledDet[kITS] = true;
850  if(TPC) fEnabledDet[kTPC] = true;
851  if(TOF) fEnabledDet[kTOF] = true;
852  if(HMPID) fEnabledDet[kHMPID] = true;
853 }
854 
855 //________________________________________________________________________
857 {
858  // function to check if a proper vertex is reconstructed and write z-position in vertexZ
859  const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
860  if(!vertex){
861  fHistNEvents->Fill(2);
862  return false;
863  }
864  else{
865  TString title=vertex->GetTitle();
866  if(title.Contains("Z") || title.Contains("3D")) return false;
867  if(vertex->GetNContributors()<1) {
868  fHistNEvents->Fill(3);
869  return false;
870  }
871  }
872 
873  const AliVVertex *vSPD = fAOD->GetPrimaryVertexSPD();
874  if(!vSPD || (vSPD && vSPD->GetNContributors()<1)){
875  fHistNEvents->Fill(4);
876  return false;
877  }
878  else{
879  double dz = vSPD->GetZ()-vertex->GetZ();
880  if(TMath::Abs(dz)>0.5) {
881  fHistNEvents->Fill(5);
882  return false;
883  }
884  double covTrc[6],covSPD[6];
885  vertex->GetCovarianceMatrix(covTrc);
886  vSPD->GetCovarianceMatrix(covSPD);
887  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
888  double errTrc = TMath::Sqrt(covTrc[5]);
889  double nsigTot = TMath::Abs(dz)/errTot, nsigTrc = TMath::Abs(dz)/errTrc;
890  if (TMath::Abs(dz)>0.2 || nsigTot>10 || nsigTrc>20) {
891  fHistNEvents->Fill(5);
892  return false;
893  }
894  }
895 
896  if(TMath::Abs(vertex->GetZ())>10.) {
897  fHistNEvents->Fill(6);
898  return false;
899  }
900  fHistNEvents->Fill(7);
901 
902  return true;
903 }
904 
905 //________________________________________________________________________
907 
908  if(fCentEstimator==kCentOff) return true;
909 
910  AliMultSelection *multSelection = (AliMultSelection*)fAOD->FindListObject("MultSelection");
911  if(!multSelection){
912  AliWarning("AliMultSelection could not be found in the aod event list of objects");
913  return false;
914  }
915 
916  float cent=-999;
917  if(fCentEstimator==kCentV0M) cent=multSelection->GetMultiplicityPercentile("V0M");
918  else if(fCentEstimator==kCentV0A) cent=multSelection->GetMultiplicityPercentile("V0A");
919  else if(fCentEstimator==kCentZNA) cent=multSelection->GetMultiplicityPercentile("ZNA");
920  else if(fCentEstimator==kCentCL1) cent=multSelection->GetMultiplicityPercentile("CL1");
921  else if(fCentEstimator==kCentCL0) cent=multSelection->GetMultiplicityPercentile("CL0");
922  else {
923  AliWarning(Form("CENTRALITY ESTIMATE WITH ESTIMATOR %d NOT YET IMPLEMENTED FOR NEW FRAMEWORK",fCentEstimator));
924  return false;
925  }
926 
927  if(cent>=fCentMin && cent<=fCentMax) return true;
928  return false;
929 }
930 
931 //________________________________________________________________________
932 void AliAnalysisTaskSEHFSystPID::GetTaggedV0s(vector<short> &idPionFromK0s, vector<short> &idPionFromL, vector<short> &idProtonFromL, vector<short> &idElectronFromGamma) {
933  // tag tracks from V0 decays
934  const int nV0s = fAOD->GetNumberOfV0s();
935  AliAODv0 *V0=nullptr;
936 
937  for (int iV0=0; iV0<nV0s; iV0++){
938  V0 = (AliAODv0*)fAOD->GetV0(iV0);
939  if(!V0) continue;
940  //if(V0->GetOnFlyStatus()) continue;
941 
942  AliAODTrack* pTrack=dynamic_cast<AliAODTrack*>(V0->GetDaughter(0));
943  AliAODTrack* nTrack=dynamic_cast<AliAODTrack*>(V0->GetDaughter(1));
944  if(!fESDtrackCuts->IsSelected(pTrack) || !fESDtrackCuts->IsSelected(nTrack)) continue;
945 
946  // Get the particle selection
947  bool foundV0 = false;
948  int pdgV0=-1, pdgP=-1, pdgN=-1;
949  foundV0 = fV0cuts->ProcessV0(V0, pdgV0, pdgP, pdgN);
950  if(!foundV0) continue;
951 
952  // v0 Armenteros plot (QA)
953  fHistArmenteroPlot[0]->Fill(V0->AlphaV0(),V0->PtArmV0());
954 
955  short iTrackP = V0->GetPosID(); // positive track
956  short iTrackN = V0->GetNegID(); // negative track
957 
958  if(pdgP==211 && pdgN==-211) {
959  idPionFromK0s.push_back(iTrackP);
960  idPionFromK0s.push_back(iTrackN);
961  fHistArmenteroPlot[1]->Fill(V0->AlphaV0(),V0->PtArmV0());
962  }
963  else if(pdgP==2212 && pdgN==-211) {
964  idProtonFromL.push_back(iTrackP);
965  idPionFromL.push_back(iTrackN);
966  fHistArmenteroPlot[2]->Fill(V0->AlphaV0(),V0->PtArmV0());
967  }
968  else if(pdgP==211 && pdgN==-2212) {
969  idPionFromL.push_back(iTrackP);
970  idProtonFromL.push_back(iTrackN);
971  fHistArmenteroPlot[3]->Fill(V0->AlphaV0(),V0->PtArmV0());
972  }
973  else if(pdgP==-11 && pdgN==11) {
974  idElectronFromGamma.push_back(iTrackP);
975  idElectronFromGamma.push_back(iTrackN);
976  fHistArmenteroPlot[4]->Fill(V0->AlphaV0(),V0->PtArmV0());
977  }
978  }
979 }
980 
981 //________________________________________________________________________
982 int AliAnalysisTaskSEHFSystPID::GetPDGcodeFromMC(AliAODTrack* track, TClonesArray* arrayMC)
983 {
984  // Get pdg code
985  int pdg = -1;
986  if(!track) return pdg;
987  int label = track->GetLabel();
988  if(label<0) return pdg;
989  AliAODMCParticle* partMC = dynamic_cast<AliAODMCParticle*>(arrayMC->At(label));
990  if(!partMC) return pdg;
991  pdg = TMath::Abs(partMC->GetPdgCode());
992 
993  return pdg;
994 }
995 
996 //___________________________________________________________________
997 AliAODTrack* AliAnalysisTaskSEHFSystPID::IsKinkDaughter(AliAODTrack* track)
998 {
999  // Check if track is a kink daughter --> in this case returns the mother track
1000  AliAODVertex *maybeKink=track->GetProdVertex();
1001  if(!maybeKink) return nullptr;
1002  if(maybeKink->GetType()==AliAODVertex::kKink) return ((AliAODTrack *)maybeKink->GetParent());
1003 
1004  return nullptr;
1005 }
1006 
1007 //___________________________________________________________________
1008 void AliAnalysisTaskSEHFSystPID::GetTaggedKaonsFromKinks(vector<short> &idKaonFromKinks)
1009 {
1010  // Tag kink mother tracks
1011 
1012  const int nTracks = fAOD->GetNumberOfTracks();
1013  for(int iTrack=0; iTrack<nTracks; iTrack++) {
1014 
1015  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
1016  if(track->GetTPCNcls()<30 || track->IsPrimaryCandidate()) continue;
1017  AliAODTrack* kinkmothertrack = IsKinkDaughter(track);
1018  if(!kinkmothertrack || !fESDtrackCuts->IsSelected(kinkmothertrack) || !kinkmothertrack->IsPrimaryCandidate()) continue;
1019  AliAODVertex *prodvtx=track->GetProdVertex();
1020  float R = TMath::Sqrt(prodvtx->GetX()*prodvtx->GetX()+prodvtx->GetY()*prodvtx->GetY());
1021  if(R<fRMinKinks || R>fRMaxKinks) continue;
1022  if(fPIDresp->NumberOfSigmasTPC(track,AliPID::kMuon)>3.5 && fPIDresp->NumberOfSigmasTOF(track,AliPID::kMuon)>3.5) continue;
1023 
1024  //kinks at vertex
1025  float bfield = fAOD->GetMagneticField();
1026  AliESDtrack kinkMotherESD(kinkmothertrack);
1027  AliESDtrack kinkDaughterESD(track);
1028  double xKinkMother=-999;
1029  double xKinkDaughter=-999;
1030  double minDist = kinkMotherESD.GetDCA(&kinkDaughterESD, bfield, xKinkMother, xKinkDaughter);
1031  if (minDist > 2.0) continue;
1032  double mother3MomentumKinkVtxArray[3] = {0, 0, 0};
1033  double daughter3MomentumKinkVtxArray[3] = {0, 0, 0};
1034  if (!kinkMotherESD.GetPxPyPzAt(xKinkMother, bfield, mother3MomentumKinkVtxArray)) continue;
1035  if (!kinkDaughterESD.GetPxPyPzAt(xKinkDaughter, bfield, daughter3MomentumKinkVtxArray)) continue;
1036  TVector3 mother3MomentumDCA(mother3MomentumKinkVtxArray);
1037  TVector3 daughter3MomentumDCA(daughter3MomentumKinkVtxArray);
1038  float qt = daughter3MomentumDCA.Perp(mother3MomentumDCA);
1039  float openingAngle = daughter3MomentumDCA.Angle(mother3MomentumDCA);
1040  openingAngle *= 180/TMath::Pi();
1041  if(qt<fQtMinKinks || qt>0.300) continue;
1042  if(openingAngle<2 || openingAngle>MaxOpeningAngleKnu(kinkmothertrack->P())) continue;
1043 
1044  //inv-mass
1045  float massmu = TDatabasePDG::Instance()->GetParticle(13)->Mass(); //muon
1046  TVector3 transferedMom = mother3MomentumDCA-daughter3MomentumDCA;
1047  float energyDaughterMu = TMath::Sqrt(daughter3MomentumDCA.Mag()*daughter3MomentumDCA.Mag()+massmu*massmu);
1048  float invmassMuNu = (energyDaughterMu+transferedMom.Mag())*(energyDaughterMu+transferedMom.Mag())-mother3MomentumDCA.Mag()*mother3MomentumDCA.Mag();
1049  if(invmassMuNu>0) invmassMuNu = TMath::Sqrt(invmassMuNu);
1050  if(invmassMuNu<0 || invmassMuNu>0.8) continue;
1051 
1052  //N-clusters vs. R
1053  int TPCnclsMother = kinkmothertrack->GetTPCNcls();
1054  int TPCnclsMax = -31.67+(11./12.)*R;
1055  int TPCnclsMin = -85.5 +(65./95.)*R;
1056  if(TPCnclsMother<TPCnclsMin || TPCnclsMother>TPCnclsMax) continue;
1057 
1058  idKaonFromKinks.push_back(kinkmothertrack->GetID());
1059  fHistQtVsMassKinks->Fill(invmassMuNu,qt);
1060  fHistPDaughterVsMotherKink->Fill(kinkmothertrack->P(),track->P());
1061  fHistOpeningAngleVsPMotherKink->Fill(kinkmothertrack->P(),openingAngle);
1062  fHistdEdxVsPMotherKink->Fill(kinkmothertrack->P(),kinkmothertrack->GetTPCsignal());
1063  fHistNTPCclsVsRadius->Fill(R,TPCnclsMother);
1064  }
1065 }
1066 
1067 //___________________________________________________________________
1069 
1070  float par0 = 0.493677;
1071  float par1 = 0.9127037;
1072  float par2 = TMath::Pi();
1073 
1074  return TMath::ATan(par0*par1*1./TMath::Sqrt(p*p*(1-par1*par1)-(par0*par0*par1*par1)))*180/par2;
1075 }
1076 
1077 //________________________________________________________________
1079 {
1080  float t_d = fPIDresp->GetTOFResponse().GetExpectedSignal(track, AliPID::kTriton); //largest mass possible with Z=1
1081  float len = track->GetIntegratedLength();
1082  float beta_d = len / (t_d * kCSPEED);
1083  float mass = AliPID::ParticleMassZ(AliPID::kTriton); //largest mass possible with Z=1
1084 
1085  if(TMath::Abs(beta_d-1.) < 1.e-12) return track->GetTPCmomentum();
1086  else return mass*beta_d/sqrt(1.-(beta_d*beta_d));
1087 }
1088 
1089 //___________________________________________________________________
1091  if(num>=static_cast<float>(numeric_limits<short>::max())) return numeric_limits<short>::max();
1092  else if(num<=static_cast<float>(numeric_limits<short>::min())) return numeric_limits<short>::min();
1093 
1094  num = round(num);
1095  return static_cast<short>(num);
1096 }
1097 
1098 //___________________________________________________________________
1100  if(num>=static_cast<float>(numeric_limits<unsigned short>::max())) return numeric_limits<unsigned short>::max();
1101  else if(num<=static_cast<float>(numeric_limits<unsigned short>::min())) return numeric_limits<unsigned short>::min();
1102 
1103  num = round(num);
1104  return static_cast<unsigned short>(num);
1105 }
1106 
1107 //________________________________________________________________
1108 void AliAnalysisTaskSEHFSystPID::GetNsigmaTPCMeanSigmaData(float &mean, float &sigma, AliPID::EParticleType species, float pTPC, float eta) {
1109 
1110  int bin = TMath::BinarySearch(fNPbinsNsigmaTPCDataCorr,fPlimitsNsigmaTPCDataCorr,pTPC);
1111  if(bin<0) bin=0; //underflow --> equal to min value
1112  else if(bin>fNPbinsNsigmaTPCDataCorr-1) bin=fNPbinsNsigmaTPCDataCorr-1; //overflow --> equal to max value
1113 
1114  int etabin = TMath::BinarySearch(fNEtabinsNsigmaTPCDataCorr,fEtalimitsNsigmaTPCDataCorr,TMath::Abs(eta));
1115  if(etabin<0) etabin=0; //underflow --> equal to min value
1116  else if(etabin>fNEtabinsNsigmaTPCDataCorr-1) etabin=fNEtabinsNsigmaTPCDataCorr-1; //overflow --> equal to max value
1117 
1118  switch(species) {
1119  case AliPID::kPion:
1120  {
1121  mean = fMeanNsigmaTPCPionData[etabin][bin];
1122  sigma = fSigmaNsigmaTPCPionData[etabin][bin];
1123  break;
1124  }
1125  case AliPID::kKaon:
1126  {
1127  mean = fMeanNsigmaTPCKaonData[etabin][bin];
1128  sigma = fSigmaNsigmaTPCKaonData[etabin][bin];
1129  break;
1130  }
1131  case AliPID::kProton:
1132  {
1133  mean = fMeanNsigmaTPCProtonData[etabin][bin];
1134  sigma = fSigmaNsigmaTPCProtonData[etabin][bin];
1135  break;
1136  }
1137  default:
1138  {
1139  mean = 0.;
1140  sigma = 1.;
1141  break;
1142  }
1143  }
1144 }
1145 
1146 //________________________________________________________________
1148 
1149  int run = fAOD->GetRunNumber();
1150  if(fAOD->GetRunNumber()!=fRunNumberPrevEvent) {
1151  if(run >= 244917 && run <= 246994)
1152  fAliEventCuts.SetupRun2PbPb();
1153  else if(run >= 295369 && run <= 297624)
1154  fAliEventCuts.SetupPbPb2018();
1155  }
1156 
1157  // cut on correlations for out of bunch pileup in PbPb run2
1159  if(!fAliEventCuts.PassedCut(AliEventCuts::kCorrelations))
1160  return 8;
1161  }
1162  else if(fApplyPbPbOutOfBunchPileupCuts==2 && run >= 295369 && run <= 297624){
1163  // Ionut cut on V0multiplicity vs. n TPC clusters (Pb-Pb 2018)
1164  AliAODVZERO* v0data=(AliAODVZERO*)((AliAODEvent*)fAOD)->GetVZEROData();
1165  float mTotV0=v0data->GetMTotV0A()+v0data->GetMTotV0C();
1166  int nTPCcls=((AliAODEvent*)fAOD)->GetNumberOfTPCClusters();
1167  float mV0TPCclsCut=-2000.+(0.013*nTPCcls)+(1.25e-9*nTPCcls*nTPCcls);
1168  if(mTotV0<mV0TPCclsCut){
1169  return 9;
1170  }
1171  }
1172 
1174  fAliEventCuts.UseTimeRangeCut();
1175  if(!fAliEventCuts.PassedCut(AliEventCuts::kTriggerClasses))
1176  return 10;
1177  }
1178 
1179  return 0;
1180 }
1181 
1182 //________________________________________________________________
1184 
1185  // convert to ESD track here
1186  AliESDtrack esdTrack(track);
1187  // set the TPC cluster info
1188  esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
1189  esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
1190  esdTrack.SetTPCPointsF(track->GetTPCNclsF());
1191 
1192  float nCrossedRowsTPC = esdTrack.GetTPCCrossedRows();
1193  float lengthInActiveZoneTPC=esdTrack.GetLengthInActiveZone(0,fDeadZoneWidth,220.,fAOD->GetMagneticField());
1194  double cutGeoNcrNclLength=fCutGeoNcrNclLength-TMath::Power(TMath::Abs(esdTrack.GetSigned1Pt()),fCutGeoNcrNclGeom1Pt);
1195  if (lengthInActiveZoneTPC<cutGeoNcrNclLength)
1196  return false;
1197  if (nCrossedRowsTPC<fCutGeoNcrNclFractionNcr*cutGeoNcrNclLength)
1198  return false;
1199  if (esdTrack.GetTPCncls()<fCutGeoNcrNclFractionNcl*cutGeoNcrNclLength)
1200  return false;
1201 
1202  return true;
1203 }
1204 
1205 //________________________________________________________________
1206 bool AliAnalysisTaskSEHFSystPID::FillNsigma(int iDet, AliAODTrack* track) {
1207 
1208  if(iDet>=kNMaxDet)
1209  return false;
1210 
1211  AliPID::EParticleType hypopidresp[7] = {AliPID::kPion,AliPID::kKaon,AliPID::kProton,AliPID::kElectron,AliPID::kDeuteron,AliPID::kTriton,AliPID::kHe3};
1212 
1213  switch(iDet) {
1214  case kITS:
1215  {
1216  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++)
1217  fPIDNsigma[iDet][iHypo] = ConvertFloatToShort(fPIDresp->NumberOfSigmasITS(track,hypopidresp[iHypo])*100);
1218  break;
1219  }
1220  case kTPC:
1221  {
1222  float nSigmaTPC[kNMaxHypo];
1223 
1224  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++)
1225  nSigmaTPC[iHypo] = fPIDresp->NumberOfSigmasTPC(track,hypopidresp[iHypo]);
1226 
1227  if(fEnableNsigmaTPCDataCorr) { //only pion, kaon and protons
1228  float mean[kProton+1], sigma[kProton+1];
1229  for(int iHypo=0; iHypo<=kProton; iHypo++) {
1230  GetNsigmaTPCMeanSigmaData(mean[iHypo], sigma[iHypo], hypopidresp[iHypo], track->GetTPCmomentum(), track->Eta());
1231 
1232  if(nSigmaTPC[iHypo]>-990.)
1233  nSigmaTPC[iHypo] = (nSigmaTPC[iHypo]-mean[iHypo]) / sigma[iHypo];
1234  }
1235  }
1236 
1237  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++)
1238  fPIDNsigma[iDet][iHypo] = ConvertFloatToShort(nSigmaTPC[iHypo]*100);
1239 
1240  break;
1241  }
1242  case kTOF:
1243  {
1244  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++)
1245  fPIDNsigma[iDet][iHypo] = ConvertFloatToShort(fPIDresp->NumberOfSigmasTOF(track,hypopidresp[iHypo])*100);
1246 
1247  break;
1248  }
1249  case kHMPID:
1250  {
1251  for(int iHypo=0; iHypo<kNMaxHypo; iHypo++)
1252  fPIDNsigma[iDet][iHypo] = ConvertFloatToShort(fPIDresp->NumberOfSigmasHMPID(track,hypopidresp[iHypo])*100);
1253 
1254  break;
1255  }
1256  }
1257 
1258  return true;
1259 }
float fNsigmaMaxForNucleiTag
max nSigma value to tag kaons
Int_t pdg
int fApplyPbPbOutOfBunchPileupCuts
event-cut object for centrality correlation event cuts
void GetTaggedV0s(vector< short > &idPionFromK0s, vector< short > &idPionFromL, vector< short > &idProtonFromL, vector< short > &idElectronFromGamma)
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
float fCutGeoNcrNclLength
1st parameter geometrical cut
AliPIDResponse * fPIDresp
AOD object.
TH1F * fHistNEvents
! histo with number of events
unsigned short fPt
ToF signal to fill the tree.
unsigned short fPTOF
TPC momentum to fill the tree.
void GetTaggedKaonsFromKinks(vector< short > &idKaonFromKinks)
vector< vector< float > > fSigmaNsigmaTPCProtonData
array of NsigmaTPC kaon mean in data
bool fEnabledDet[kNMaxDet]
array of flags to enable particle species
Double_t mass
float fNsigmaMaxForKaonTag
bit map for tag (see enum above)
unsigned short fTrackLength
number of PID clusters in TPC to fill the tree
double fFracToKeepDownSampling
flag to enable/disable downsampling
static Int_t CheckMatchingAODdeltaAODevents()
AliEventCuts fAliEventCuts
flag to enable usage of AliEventCuts foe event-selection
static void SetNsigmaTPCDataDrivenCorrection(Int_t run, Int_t system, Int_t &nPbins, Float_t Plims[kMaxPBins+1], Int_t &nEtabins, Float_t absEtalims[kMaxEtaBins+1], vector< vector< Float_t > > &meanNsigmaTPCpion, vector< vector< Float_t > > &meanNsigmaTPCkaon, vector< vector< Float_t > > &meanNsigmaTPCproton, vector< vector< Float_t > > &sigmaNsigmaTPCpion, vector< vector< Float_t > > &sigmaNsigmaTPCkaon, vector< vector< Float_t > > &sigmaNsigmaTPCproton)
int fNEtabinsNsigmaTPCDataCorr
vector of eta limits for data-driven NsigmaTPC correction
int fSystNsigmaTPCDataCorr
flag to enable data-driven NsigmaTPC correction
TH2F * fHistNsigmaVsPt[kNMaxDet][kNMaxHypo]
! array of histos for nsigma vs pt (MC truth)
vector< vector< float > > fMeanNsigmaTPCPionData
system for data-driven NsigmaTPC correction
unsigned char fTPCNcls
transverse momentum to fill the tree
int GetPDGcodeFromMC(AliAODTrack *track, TClonesArray *arrayMC)
void GetNsigmaTPCMeanSigmaData(float &mean, float &sigma, AliPID::EParticleType species, float pTPC, float eta)
bool FillNsigma(int iDet, AliAODTrack *track)
TH2F * fHistOpeningAngleVsPMotherKink
! histo for opening angle vs. pT mother kink
unsigned char fTPCNcrossed
start time resolution for TOF PID
int fAODProtection
option for downsampling
TH2F * fHistNTPCclsVsRadius
! histo for nTPC clusters vs. R mother kink
TH2F * fHistQtVsMassKinks
! histo for mother-kink qt vs. mass distribution
unsigned char fTPCNclsPID
number of clusters in TPC to fill the tree
unsigned short fStartTimeRes
track length for TOF PID
Double_t * sigma
int fNPbinsNsigmaTPCDataCorr
array of p limits for data-driven NsigmaTPC correction
float fEtalimitsNsigmaTPCDataCorr[AliAODPidHF::kMaxEtaBins+1]
number of p bins for data-driven NsigmaTPC correction
virtual void UserExec(Option_t *option)
unsigned short fdEdxTPC
HMPID momentum to fill the tree.
unsigned short fToF
ITS dEdX to fill the tree.
double fPtMaxDownSampling
fraction to keep when downsampling activated
float fPlimitsNsigmaTPCDataCorr[AliAODPidHF::kMaxPBins+1]
array of NsigmaTPC proton mean in data
AliTimeRangeCut fTimeRangeCut
flag to enable time-range cut in PbPb 2018
float fCentMin
5th parameter geometrical cut
unsigned short fHMPIDsignal
ITS cluster map.
float fRMaxKinks
min radius in XY for kinks
int fDownSamplingOpt
pT max of tracks to downsample
AliAODEvent * fAOD
single-track cut set
unsigned short fPTPC
Momentum at primary vertex to fill the tree.
static const int kMaxPBins
Definition: AliAODPidHF.h:45
bool IsSelectedByGeometricalCut(AliAODTrack *track)
unsigned char fTPCFindable
number of TPC crossed rows
vector< vector< float > > fSigmaNsigmaTPCKaonData
array of NsigmaTPC pion mean in data
TH2F * fHistdEdxVsPMotherKink
! histo for mother kink TPC dEdx vs. p
float GetTOFmomentum(AliAODTrack *track)
TH2F * fHistPDaughterVsMotherKink
! histo for pT daughter vs. pT mother kink
AliESDtrackCuts * fESDtrackCuts
conversion factor from float to short for dE/dx
unsigned short fPhi
pseudorapidity of the track
bool fUseAliEventCuts
number of eta bins for data-driven NsigmaTPC correction
unsigned short fHMPIDoccupancy
HMPID signal.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
bool fEnableNsigmaTPCDataCorr
run number of previous event
float fCutGeoNcrNclFractionNcl
4th parameter geometrical cut
unsigned short fTag
PDG code in case of MC to fill the tree.
vector< vector< float > > fMeanNsigmaTPCKaonData
array of NsigmaTPC pion mean in data
void EnableParticleSpecies(bool pi=true, bool kao=true, bool pr=true, bool el=false, bool deu=false, bool tr=false, bool He3=false)
TH2F * fHistArmenteroPlot[5]
! histo for armenteros-podolanski plot
bool fUseTimeRangeCutForPbPb2018
option for Pb-Pb out-of bunch pileup cuts with AliEventCuts
TList * fOutputList
! output list for histograms
TTree * fPIDtree
! tree with PID info
AliAODv0KineCuts * fV0cuts
basic pid object
unsigned long long fTriggerMask
trigger class
int fSystem
flag to switch on the MC analysis for the efficiency estimation
unsigned short fdEdxITS
TPC dEdX to fill the tree.
static const int kMaxEtaBins
Definition: AliAODPidHF.h:44
unsigned short ConvertFloatToUnsignedShort(float num)
float fConversionFactordEdx
system: 0->pp,pPb 1->PbPb
void EnableDetectors(bool ITS=false, bool TPC=true, bool TOF=true, bool HMPID=false)
bool fFillTreeWithNsigmaPIDOnly
flag to enable filling of the tree with PID variables
const char Option_t
Definition: External.C:48
unsigned short fPHMPID
TOF momentum to fill the tree.
vector< vector< float > > fSigmaNsigmaTPCPionData
array of NsigmaTPC proton mean in data
unsigned char fITSclsMap
number of TPC findable clusters
const Double_t pi
float fDeadZoneWidth
max radius in XY for kink
TString fTriggerClass
centrality estimator
bool fFillTreeWithTrackQualityInfo
flag to enable filling of the tree with only raw variables for the PID
int fRunNumberPrevEvent
flag to activate protection against AOD-dAOD mismatch
unsigned char fTrackInfoMap
HMPID occupancy.
bool fEnabledDownSampling
flag to enable filling of the tree with track selections
short fEta
bit map with some track info (see enum above)
short fPIDNsigma[kNMaxDet][kNMaxHypo]
array of flags to enable detectors
float fCutGeoNcrNclGeom1Pt
2nd parameter geometrical cut
int fPDGcode
azimuthal angle of the track
unsigned short fP
Nsigma PID to fill the tree.
bool fFillTreeWithRawPIDOnly
flag to enable filling of the tree with only Nsigma variables for the PID
AliAODTrack * IsKinkDaughter(AliAODTrack *track)
float fCutGeoNcrNclFractionNcr
3rd parameter geometrical cut
vector< vector< float > > fMeanNsigmaTPCProtonData
array of NsigmaTPC kaon mean in data