AliPhysics  c2ade29 (c2ade29)
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 
25 
26 //________________________________________________________________________
28 AliAnalysisTaskSE("TaskNsigmaPID"),
29 fOutputList(nullptr),
30 fHistNEvents(nullptr),
31 fHistQtVsMassKinks(nullptr),
32 fHistPDaughterVsMotherKink(nullptr),
33 fHistdEdxVsPMotherKink(nullptr),
34 fHistOpeningAngleVsPMotherKink(nullptr),
35 fHistNTPCclsVsRadius(nullptr),
36 fPIDtree(nullptr),
37 fPTPC(0),
38 fPTOF(0),
39 fdEdxTPC(0),
40 fToF(0),
41 fPt(0),
42 fTPCNcls(0),
43 fTPCNclsPID(0),
44 fTrackLength(0),
45 fStartTimeRes(0),
46 fPDGcode(-1),
47 fTag(0),
48 fNsigmaMaxForTag(0.02),
49 fQtMinKinks(0.15),
50 fRMinKinks(120),
51 fRMaxKinks(210),
52 fCentMin(0.),
53 fCentMax(100.),
54 fCentEstimator(kCentOff),
55 fTriggerClass(""),
56 fTriggerMask(AliVEvent::kINT7),
57 fIsMC(false),
58 fSystem(0),
59 fESDtrackCuts(nullptr),
60 fAOD(nullptr),
61 fPIDresp(nullptr),
62 fV0cuts(nullptr),
63 fFillTreeWithNsigmaPIDOnly(false),
64 fEnabledDownSampling(false),
65 fFracToKeepDownSampling(0.1),
66 fPtMaxDownSampling(1.5)
67 {
68  //
69  // default constructur
70  //
71 
72  for(int iVar=0; iVar<6; iVar++) fPIDNsigma[iVar] = -999.;
73  for(int iHisto=0; iHisto<5; iHisto++) fHistArmenteroPlot[iHisto] = nullptr;
74  for(int iHisto=0; iHisto<kNHypo; iHisto++) {
75  fHistNsigmaTPCvsPt[iHisto] = nullptr;
76  fHistNsigmaTOFvsPt[iHisto] = nullptr;
77  }
78 }
79 
80 //________________________________________________________________________
82 AliAnalysisTaskSE(name),
83 fOutputList(nullptr),
84 fHistNEvents(nullptr),
85 fHistQtVsMassKinks(nullptr),
86 fHistPDaughterVsMotherKink(nullptr),
87 fHistdEdxVsPMotherKink(nullptr),
88 fHistOpeningAngleVsPMotherKink(nullptr),
89 fHistNTPCclsVsRadius(nullptr),
90 fPIDtree(nullptr),
91 fPTPC(0),
92 fPTOF(0),
93 fdEdxTPC(0),
94 fToF(0),
95 fPt(0),
96 fTPCNcls(0),
97 fTPCNclsPID(0),
98 fTrackLength(0),
99 fStartTimeRes(0),
100 fPDGcode(-1),
101 fTag(0),
102 fNsigmaMaxForTag(0.02),
103 fQtMinKinks(0.15),
104 fRMinKinks(120),
105 fRMaxKinks(210),
106 fCentMin(0.),
107 fCentMax(100.),
108 fCentEstimator(kCentOff),
109 fTriggerClass(""),
110 fTriggerMask(AliVEvent::kINT7),
111 fIsMC(false),
112 fSystem(system),
113 fESDtrackCuts(nullptr),
114 fAOD(nullptr),
115 fPIDresp(nullptr),
116 fV0cuts(nullptr),
117 fFillTreeWithNsigmaPIDOnly(false),
118 fEnabledDownSampling(false),
119 fFracToKeepDownSampling(0.1),
120 fPtMaxDownSampling(1.5)
121 {
122  //
123  // standard constructur
124  //
125 
126  for(int iVar=0; iVar<6; iVar++) fPIDNsigma[iVar] = -999.;
127  for(int iHisto=0; iHisto<5; iHisto++) fHistArmenteroPlot[iHisto] = nullptr;
128  for(int iHisto=0; iHisto<kNHypo; iHisto++) {
129  fHistNsigmaTPCvsPt[iHisto] = nullptr;
130  fHistNsigmaTOFvsPt[iHisto] = nullptr;
131  }
132 
133  DefineInput(0, TChain::Class());
134  DefineOutput(1, TList::Class());
135  DefineOutput(2, TTree::Class());
136 }
137 
138 //________________________________________________________________________
140 {
141  // Destructor
142  if (fOutputList) {
143  delete fHistNEvents;
144  for(int iHisto=0; iHisto<5; iHisto++) delete fHistArmenteroPlot[iHisto];
145  delete fHistQtVsMassKinks;
148  delete fHistdEdxVsPMotherKink;
149  delete fHistNTPCclsVsRadius;
150  delete fOutputList;
151  }
152 
153  if(fPIDtree) delete fPIDtree;
154  if(fESDtrackCuts) delete fESDtrackCuts;
155  if(fV0cuts) delete fV0cuts;
156 }
157 
158 //________________________________________________________________________
160 {
161  // create track cuts
162  if(!fESDtrackCuts) {
163  fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
164  fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE,0);
165  fESDtrackCuts->SetEtaRange(-0.8, 0.8);
166  fESDtrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
167  }
168 
169  // V0 Kine cuts
170  fV0cuts = new AliAODv0KineCuts;
171  fV0cuts->SetNoKinks(false);
172 
173  fOutputList = new TList();
174  fOutputList->SetOwner(true);
175 
176  fHistNEvents = new TH1F("fHistNEvents","Number of processed events;;Number of events",8,-1.5,6.5);
177  fHistNEvents->Sumw2();
178  fHistNEvents->SetMinimum(0);
179  fHistNEvents->GetXaxis()->SetBinLabel(1,"Read from AOD");
180  fHistNEvents->GetXaxis()->SetBinLabel(2,"Pass Phys. Sel. + Trig");
181  fHistNEvents->GetXaxis()->SetBinLabel(3,"No vertex");
182  fHistNEvents->GetXaxis()->SetBinLabel(4,"Vertex contributors < 1");
183  fHistNEvents->GetXaxis()->SetBinLabel(5,"Without SPD vertex");
184  fHistNEvents->GetXaxis()->SetBinLabel(6,"Error on zVertex>0.5");
185  fHistNEvents->GetXaxis()->SetBinLabel(7,"|zVertex|>10");
186  fHistNEvents->GetXaxis()->SetBinLabel(8,"Good Z vertex");
188 
189  TString armenteronames[5] = {"All","K0s","Lambda","AntiLambda","Gamma"};
190  TString armenterolabels[5] = {"all","K_{s}^{0}","#Lambda","#bar{#Lambda}","#gamma"};
191  for(int iHisto=0; iHisto<5; iHisto++) {
192  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);
193  fOutputList->Add(fHistArmenteroPlot[iHisto]);
194  }
195 
196  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);
198 
199  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);
201 
202  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);
204 
205  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);
207 
208  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);
210 
211  if(fIsMC) {
212  for(int iHisto=0; iHisto<kNHypo; iHisto++) {
213  fHistNsigmaTPCvsPt[iHisto] = new TH2F(Form("fHistNsigmaTPCvsPt_%s",hyponames[iHisto].Data()),Form(";#it{p}_{T} (GeV/#it{c});N_{#sigma}^{TPC}(%s)",hyponames[iHisto].Data()),500,0,50,1000,-50,50);
214  fHistNsigmaTOFvsPt[iHisto] = new TH2F(Form("fHistNsigmaTOFvsPt_%s",hyponames[iHisto].Data()),Form(";#it{p}_{T} (GeV/#it{c});N_{#sigma}^{TOF}(%s)",hyponames[iHisto].Data()),500,0,50,1000,-50,50);
215  fOutputList->Add(fHistNsigmaTPCvsPt[iHisto]);
216  fOutputList->Add(fHistNsigmaTOFvsPt[iHisto]);
217  }
218  }
219 
220  fPIDtree = new TTree("fPIDtree","fPIDtree");
221  TString PIDbranchnames[6] = {"n_sigma_TPC_pi","n_sigma_TPC_K","n_sigma_TPC_p","n_sigma_TOF_pi","n_sigma_TOF_K","n_sigma_TOF_p"};
222  for(int iVar=0; iVar<6; iVar++) {
223  fPIDtree->Branch(PIDbranchnames[iVar].Data(),&fPIDNsigma[iVar],Form("%s/S",PIDbranchnames[iVar].Data()));
224  }
225  fPIDtree->Branch("pT",&fPt,"pT/s");
227  fPIDtree->Branch("pTPC",&fPTPC,"pTPC/s");
228  fPIDtree->Branch("pTOF",&fPTOF,"pTOF/s");
229  fPIDtree->Branch("dEdx",&fdEdxTPC,"dEdx/s");
230  fPIDtree->Branch("ToF",&fToF,"ToF/s");
231  fPIDtree->Branch("NclusterTPC",&fTPCNcls,"NclusterTPC/b");
232  fPIDtree->Branch("NclusterPIDTPC",&fTPCNclsPID,"NclusterPIDTPC/b");
233  fPIDtree->Branch("TrackLength",&fTrackLength,"TrackLength/s");
234  fPIDtree->Branch("StartTimeRes",&fStartTimeRes,"StartTimeRes/s");
235  }
236  fPIDtree->Branch("tag",&fTag,"tag/b");
237  if(fIsMC) fPIDtree->Branch("PDGcode",&fPDGcode,"PDGcode/S");
238 
239  // post data
240  PostData(1, fOutputList);
241  PostData(2, fPIDtree);
242 }
243 
244 //________________________________________________________________________
246 {
247  // main event loop
248  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
249  if (!fAOD && AODEvent() && IsStandardAOD()) {
250  // In case there is an AOD handler writing a standard AOD, use the AOD
251  // event in memory rather than the input (ESD) event.
252  fAOD = dynamic_cast<AliAODEvent*> (AODEvent());
253  }
254  if (!fAOD) {
255  AliWarning("AliAnalysisTaskSEHFSystPID::Exec(): bad AOD");
256  PostData(1, fOutputList);
257  return;
258  }
259  fHistNEvents->Fill(-1);
260 
261  if(TMath::Abs(fAOD->GetMagneticField())<0.001) return;
262 
263  AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
264  if(!aodHandler) {
265  AliWarning("AliAnalysisTaskSEHFSystPID::Exec(): No AliInputEventHandler!");
266  return;
267  }
268 
269  unsigned int maskPhysSel = aodHandler->IsEventSelected();
270  TString firedTriggerClasses = fAOD->GetFiredTriggerClasses();
271  if(!fIsMC && (fAOD->GetRunNumber()<136851 || fAOD->GetRunNumber()>139517)) {
272  if(!(firedTriggerClasses.Contains(fTriggerClass.Data()))) return;
273  }
274  if((maskPhysSel & fTriggerMask)==0) return;
275 
276  if(!IsCentralitySelected()) {
277  PostData(1, fOutputList);
278  return;
279  }
280 
281  fHistNEvents->Fill(0);
282 
283  if (!IsVertexAccepted()) {
284  fHistNEvents->Fill(1);
285  PostData(1, fOutputList);
286  return;
287  }
288 
289  // load MC particles
290  TClonesArray *arrayMC=0;
291  if(fIsMC){
292  arrayMC = (TClonesArray*)fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
293  if(!arrayMC) {
294  AliWarning("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
295  return;
296  }
297  }
298 
299  //get pid response
300  if(!fPIDresp) fPIDresp = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
301 
302  // V0 selection
303  if(fSystem==0) fV0cuts->SetMode(AliAODv0KineCuts::kPurity,AliAODv0KineCuts::kPP);
304  else if(fSystem==1) fV0cuts->SetMode(AliAODv0KineCuts::kPurity,AliAODv0KineCuts::kPbPb);
305  fV0cuts->SetEvent(fAOD);
306 
307  vector<short> idPionFromK0s;
308  vector<short> idPionFromL;
309  vector<short> idProtonFromL;
310  vector<short> idElectronFromGamma;
311  vector<short> idKaonFromKinks;
312  GetTaggedV0s(idPionFromK0s,idPionFromL,idProtonFromL,idElectronFromGamma);
313  GetTaggedKaonsFromKinks(idKaonFromKinks);
314 
315  vector<short>::iterator it;
316 
317  const int nTracks = fAOD->GetNumberOfTracks();
318  //loop on tracks
319 
320  for(int iTrack=0; iTrack<nTracks; iTrack++) {
321  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
322  if(!track) continue;
323  //applying ESDtrackCut
324  if(!fESDtrackCuts->IsSelected(track)) continue;
325 
326  if(fEnabledDownSampling && track->Pt()<fPtMaxDownSampling) {
327  double pseudoRand = track->Pt()*1000.-(long)(track->Pt()*1000.);
328  if(pseudoRand>fFracToKeepDownSampling) continue;
329  }
330 
331  fPt = ConvertFloatToUnsignedShort(track->Pt()*1000);
332 
334  //TPC variables
335  fTPCNcls = static_cast<unsigned char>(track->GetTPCNcls());
336  fPTPC = ConvertFloatToUnsignedShort(track->GetTPCmomentum()*1000);
337  fdEdxTPC = ConvertFloatToUnsignedShort(track->GetTPCsignal()*100);
338  fTPCNclsPID = static_cast<unsigned char>(track->GetTPCsignalN());
339 
340  //TOF variables
342  fTrackLength = ConvertFloatToUnsignedShort(track->GetIntegratedLength()*10);
343  fStartTimeRes = ConvertFloatToUnsignedShort(fPIDresp->GetTOFResponse().GetStartTimeRes(track->P())*100);
344 
345  if (!(track->GetStatus() & AliESDtrack::kTOFout) || !(track->GetStatus() & AliESDtrack::kTIME)) {
346  fToF = 0;
347  }
348  else {
349  if (fTrackLength < 3500) fToF = 0;
350  else {
351  float tof = track->GetTOFsignal();
352  float time0 = fPIDresp->GetTOFResponse().GetStartTime(track->P());
353  fToF = ConvertFloatToUnsignedShort((tof-time0)/10);
354  }
355  }
356  }
357 
358  bool isTPCok = false;
359  bool isTOFok = false;
360  if (fPIDresp->CheckPIDStatus(AliPIDResponse::kTPC,track) == AliPIDResponse::kDetPidOk) isTPCok = true;
361  if (fPIDresp->CheckPIDStatus(AliPIDResponse::kTOF,track) == AliPIDResponse::kDetPidOk) isTOFok = true;
362 
363  if(isTPCok) {
364  fPIDNsigma[0] = ConvertFloatToShort(fPIDresp->NumberOfSigmasTPC(track,AliPID::kPion)*100);
365  fPIDNsigma[1] = ConvertFloatToShort(fPIDresp->NumberOfSigmasTPC(track,AliPID::kKaon)*100);
366  fPIDNsigma[2] = ConvertFloatToShort(fPIDresp->NumberOfSigmasTPC(track,AliPID::kProton)*100);
367  }
368  else for(int iVar=0; iVar<3; iVar++) fPIDNsigma[iVar] = numeric_limits<short>::min();
369  if(isTOFok) {
370  fPIDNsigma[3] = ConvertFloatToShort(fPIDresp->NumberOfSigmasTOF(track,AliPID::kPion)*100);
371  fPIDNsigma[4] = ConvertFloatToShort(fPIDresp->NumberOfSigmasTOF(track,AliPID::kKaon)*100);
372  fPIDNsigma[5] = ConvertFloatToShort(fPIDresp->NumberOfSigmasTOF(track,AliPID::kProton)*100);
373  }
374  else for(int iVar=3; iVar<6; iVar++) fPIDNsigma[iVar] = numeric_limits<short>::min();
375 
376  short trackid = track->GetID();
377 
378  it = find(idPionFromK0s.begin(),idPionFromK0s.end(),trackid);
379  if(it!=idPionFromK0s.end()) fTag |= kIsPionFromK0s;
380  else fTag &= ~kIsPionFromK0s;
381  it = find(idPionFromL.begin(),idPionFromL.end(),trackid);
382  if(it!=idPionFromL.end()) fTag |= kIsPionFromL;
383  else fTag &= ~kIsPionFromL;
384  it = find(idProtonFromL.begin(),idProtonFromL.end(),trackid);
385  if(it!=idProtonFromL.end()) fTag |= kIsProtonFromL;
386  else fTag &= ~kIsProtonFromL;
387  it = find(idElectronFromGamma.begin(),idElectronFromGamma.end(),trackid);
388  if(it!=idElectronFromGamma.end()) fTag |= kIsElectronFromGamma;
389  else fTag &= ~kIsElectronFromGamma;
390  it = find(idKaonFromKinks.begin(),idKaonFromKinks.end(),trackid);
391  if(it!=idKaonFromKinks.end()) fTag |= kIsKaonFromKinks;
392  else fTag &= ~kIsKaonFromKinks;
393  if(TMath::Abs(fPIDresp->NumberOfSigmasTOF(track,AliPID::kKaon))<fNsigmaMaxForTag && isTOFok && isTPCok) fTag |= kIsKaonFromTOF;
394  else fTag &= ~kIsKaonFromTOF;
395  if(TMath::Abs(fPIDresp->NumberOfSigmasTPC(track,AliPID::kKaon))<fNsigmaMaxForTag && isTOFok && isTPCok) fTag |= kIsKaonFromTPC;
396  else fTag &= ~kIsKaonFromTPC;
397 
398  if(fIsMC) {
399  fPDGcode = GetPDGcodeFromMC(track,arrayMC);
400  if(fPDGcode==211) {
401  fHistNsigmaTPCvsPt[kPion]->Fill(track->Pt(),fPIDresp->NumberOfSigmasTPC(track,AliPID::kPion));
402  fHistNsigmaTOFvsPt[kPion]->Fill(track->Pt(),fPIDresp->NumberOfSigmasTOF(track,AliPID::kPion));
403  }
404  else if(fPDGcode==321) {
405  fHistNsigmaTPCvsPt[kKaon]->Fill(track->Pt(),fPIDresp->NumberOfSigmasTPC(track,AliPID::kKaon));
406  fHistNsigmaTOFvsPt[kKaon]->Fill(track->Pt(),fPIDresp->NumberOfSigmasTOF(track,AliPID::kKaon));
407  }
408  else if(fPDGcode==2212) {
409  fHistNsigmaTPCvsPt[kProton]->Fill(track->Pt(),fPIDresp->NumberOfSigmasTPC(track,AliPID::kProton));
410  fHistNsigmaTOFvsPt[kProton]->Fill(track->Pt(),fPIDresp->NumberOfSigmasTOF(track,AliPID::kProton));
411  }
412  }
413  else fPDGcode = 0;
414 
415  if(fTag!=0 && ((fIsMC && fPDGcode>=0) || !fIsMC)) fPIDtree->Fill();
416 
417  fTag = 0;
418 
419  // Post output data
420  PostData(2, fPIDtree);
421  }
422 
423  //clear vectors of ids
424  idPionFromK0s.clear();
425  idPionFromL.clear();
426  idProtonFromL.clear();
427  idElectronFromGamma.clear();
428  idKaonFromKinks.clear();
429 
430  // Post output data
431  PostData(1, fOutputList);
432 }
433 
434 //________________________________________________________________________
436 {
437  // function to check if a proper vertex is reconstructed and write z-position in vertexZ
438  const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
439  if(!vertex){
440  fHistNEvents->Fill(1);
441  return false;
442  }
443  else{
444  TString title=vertex->GetTitle();
445  if(title.Contains("Z") || title.Contains("3D")) return false;
446  if(vertex->GetNContributors()<1) {
447  fHistNEvents->Fill(2);
448  return false;
449  }
450  }
451 
452  const AliVVertex *vSPD = fAOD->GetPrimaryVertexSPD();
453  if(!vSPD || (vSPD && vSPD->GetNContributors()<1)){
454  fHistNEvents->Fill(3);
455  return false;
456  }
457  else{
458  double dz = vSPD->GetZ()-vertex->GetZ();
459  if(TMath::Abs(dz)>0.5) {
460  fHistNEvents->Fill(4);
461  return false;
462  }
463  double covTrc[6],covSPD[6];
464  vertex->GetCovarianceMatrix(covTrc);
465  vSPD->GetCovarianceMatrix(covSPD);
466  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
467  double errTrc = TMath::Sqrt(covTrc[5]);
468  double nsigTot = TMath::Abs(dz)/errTot, nsigTrc = TMath::Abs(dz)/errTrc;
469  if (TMath::Abs(dz)>0.2 || nsigTot>10 || nsigTrc>20) {
470  fHistNEvents->Fill(4);
471  return false;
472  }
473  }
474 
475  if(TMath::Abs(vertex->GetZ())>10.) {
476  fHistNEvents->Fill(5);
477  return false;
478  }
479  fHistNEvents->Fill(6);
480 
481  return true;
482 }
483 
484 //________________________________________________________________________
486 
487  if(fCentEstimator==kCentOff) return true;
488 
489  AliMultSelection *multSelection = (AliMultSelection*)fAOD->FindListObject("MultSelection");
490  if(!multSelection){
491  AliWarning("AliMultSelection could not be found in the aod event list of objects");
492  return false;
493  }
494 
495  float cent=-999;
496  if(fCentEstimator==kCentV0M) cent=multSelection->GetMultiplicityPercentile("V0M");
497  else if(fCentEstimator==kCentV0A) cent=multSelection->GetMultiplicityPercentile("V0A");
498  else if(fCentEstimator==kCentZNA) cent=multSelection->GetMultiplicityPercentile("ZNA");
499  else if(fCentEstimator==kCentCL1) cent=multSelection->GetMultiplicityPercentile("CL1");
500  else if(fCentEstimator==kCentCL0) cent=multSelection->GetMultiplicityPercentile("CL0");
501  else {
502  AliWarning(Form("CENTRALITY ESTIMATE WITH ESTIMATOR %d NOT YET IMPLEMENTED FOR NEW FRAMEWORK",fCentEstimator));
503  return false;
504  }
505 
506  if(cent>=fCentMin && cent<=fCentMax) return true;
507  return false;
508 }
509 
510 //________________________________________________________________________
511 void AliAnalysisTaskSEHFSystPID::GetTaggedV0s(vector<short> &idPionFromK0s, vector<short> &idPionFromL, vector<short> &idProtonFromL, vector<short> &idElectronFromGamma) {
512  // tag tracks from V0 decays
513  const int nV0s = fAOD->GetNumberOfV0s();
514  AliAODv0 *V0=nullptr;
515 
516  for (int iV0=0; iV0<nV0s; iV0++){
517  V0 = (AliAODv0*)fAOD->GetV0(iV0);
518  if(!V0) continue;
519  //if(V0->GetOnFlyStatus()) continue;
520 
521  AliAODTrack* pTrack=dynamic_cast<AliAODTrack*>(V0->GetDaughter(0));
522  AliAODTrack* nTrack=dynamic_cast<AliAODTrack*>(V0->GetDaughter(1));
523  if(!fESDtrackCuts->IsSelected(pTrack) || !fESDtrackCuts->IsSelected(nTrack)) continue;
524 
525  // Get the particle selection
526  bool foundV0 = false;
527  int pdgV0=-1, pdgP=-1, pdgN=-1;
528  foundV0 = fV0cuts->ProcessV0(V0, pdgV0, pdgP, pdgN);
529  if(!foundV0) continue;
530 
531  // v0 Armenteros plot (QA)
532  fHistArmenteroPlot[0]->Fill(V0->AlphaV0(),V0->PtArmV0());
533 
534  short iTrackP = V0->GetPosID(); // positive track
535  short iTrackN = V0->GetNegID(); // negative track
536 
537  if(pdgP==211 && pdgN==-211) {
538  idPionFromK0s.push_back(iTrackP);
539  idPionFromK0s.push_back(iTrackN);
540  fHistArmenteroPlot[1]->Fill(V0->AlphaV0(),V0->PtArmV0());
541  }
542  else if(pdgP==2212 && pdgN==-211) {
543  idProtonFromL.push_back(iTrackP);
544  idPionFromL.push_back(iTrackN);
545  fHistArmenteroPlot[2]->Fill(V0->AlphaV0(),V0->PtArmV0());
546  }
547  else if(pdgP==211 && pdgN==-2212) {
548  idPionFromL.push_back(iTrackP);
549  idProtonFromL.push_back(iTrackN);
550  fHistArmenteroPlot[3]->Fill(V0->AlphaV0(),V0->PtArmV0());
551  }
552  else if(pdgP==-11 && pdgN==11) {
553  idElectronFromGamma.push_back(iTrackP);
554  idElectronFromGamma.push_back(iTrackN);
555  fHistArmenteroPlot[4]->Fill(V0->AlphaV0(),V0->PtArmV0());
556  }
557  }
558 }
559 
560 //________________________________________________________________________
561 short AliAnalysisTaskSEHFSystPID::GetPDGcodeFromMC(AliAODTrack* track, TClonesArray* arrayMC)
562 {
563  // Get pdg code
564  short pdg = -1;
565  int label = track->GetLabel();
566  if(label<0) return pdg;
567  AliAODMCParticle* partMC = dynamic_cast<AliAODMCParticle*>(arrayMC->At(label));
568  if(!partMC) return pdg;
569  pdg = TMath::Abs(partMC->GetPdgCode());
570  if(partMC->GetPdgCode()>numeric_limits<short>::max()) pdg = -1;
571 
572  return pdg;
573 }
574 
575 //___________________________________________________________________
576 AliAODTrack* AliAnalysisTaskSEHFSystPID::IsKinkDaughter(AliAODTrack* track)
577 {
578  // Check if track is a kink daughter --> in this case returns the mother track
579  AliAODVertex *maybeKink=track->GetProdVertex();
580  if(!maybeKink) return nullptr;
581  if(maybeKink->GetType()==AliAODVertex::kKink) return ((AliAODTrack *)maybeKink->GetParent());
582 
583  return nullptr;
584 }
585 
586 //___________________________________________________________________
587 void AliAnalysisTaskSEHFSystPID::GetTaggedKaonsFromKinks(vector<short> &idKaonFromKinks)
588 {
589  // Tag kink mother tracks
590 
591  const int nTracks = fAOD->GetNumberOfTracks();
592  for(int iTrack=0; iTrack<nTracks; iTrack++) {
593 
594  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
595  if(track->GetTPCNcls()<30 || track->IsPrimaryCandidate()) continue;
596  AliAODTrack* kinkmothertrack = IsKinkDaughter(track);
597  if(!kinkmothertrack || !fESDtrackCuts->IsSelected(kinkmothertrack) || !kinkmothertrack->IsPrimaryCandidate()) continue;
598  AliAODVertex *prodvtx=track->GetProdVertex();
599  float R = TMath::Sqrt(prodvtx->GetX()*prodvtx->GetX()+prodvtx->GetY()*prodvtx->GetY());
600  if(R<fRMinKinks || R>fRMaxKinks) continue;
601  if(fPIDresp->NumberOfSigmasTPC(track,AliPID::kMuon)>3.5 && fPIDresp->NumberOfSigmasTOF(track,AliPID::kMuon)>3.5) continue;
602 
603  //kinks at vertex
604  float bfield = fAOD->GetMagneticField();
605  AliESDtrack kinkMotherESD(kinkmothertrack);
606  AliESDtrack kinkDaughterESD(track);
607  double xKinkMother=-999;
608  double xKinkDaughter=-999;
609  double minDist = kinkMotherESD.GetDCA(&kinkDaughterESD, bfield, xKinkMother, xKinkDaughter);
610  if (minDist > 2.0) continue;
611  double mother3MomentumKinkVtxArray[3] = {0, 0, 0};
612  double daughter3MomentumKinkVtxArray[3] = {0, 0, 0};
613  if (!kinkMotherESD.GetPxPyPzAt(xKinkMother, bfield, mother3MomentumKinkVtxArray)) continue;
614  if (!kinkDaughterESD.GetPxPyPzAt(xKinkDaughter, bfield, daughter3MomentumKinkVtxArray)) continue;
615  TVector3 mother3MomentumDCA(mother3MomentumKinkVtxArray);
616  TVector3 daughter3MomentumDCA(daughter3MomentumKinkVtxArray);
617  float qt = daughter3MomentumDCA.Perp(mother3MomentumDCA);
618  float openingAngle = daughter3MomentumDCA.Angle(mother3MomentumDCA);
619  openingAngle *= 180/TMath::Pi();
620  if(qt<fQtMinKinks || qt>0.300) continue;
621  if(openingAngle<2 || openingAngle>MaxOpeningAngleKnu(kinkmothertrack->P())) continue;
622 
623  //inv-mass
624  float massmu = TDatabasePDG::Instance()->GetParticle(13)->Mass(); //muon
625  TVector3 transferedMom = mother3MomentumDCA-daughter3MomentumDCA;
626  float energyDaughterMu = TMath::Sqrt(daughter3MomentumDCA.Mag()*daughter3MomentumDCA.Mag()+massmu*massmu);
627  float invmassMuNu = (energyDaughterMu+transferedMom.Mag())*(energyDaughterMu+transferedMom.Mag())-mother3MomentumDCA.Mag()*mother3MomentumDCA.Mag();
628  if(invmassMuNu>0) invmassMuNu = TMath::Sqrt(invmassMuNu);
629  if(invmassMuNu<0 || invmassMuNu>0.8) continue;
630 
631  //N-clusters vs. R
632  int TPCnclsMother = kinkmothertrack->GetTPCNcls();
633  int TPCnclsMax = -31.67+(11./12.)*R;
634  int TPCnclsMin = -85.5 +(65./95.)*R;
635  if(TPCnclsMother<TPCnclsMin || TPCnclsMother>TPCnclsMax) continue;
636 
637  idKaonFromKinks.push_back(kinkmothertrack->GetID());
638  fHistQtVsMassKinks->Fill(invmassMuNu,qt);
639  fHistPDaughterVsMotherKink->Fill(kinkmothertrack->P(),track->P());
640  fHistOpeningAngleVsPMotherKink->Fill(kinkmothertrack->P(),openingAngle);
641  fHistdEdxVsPMotherKink->Fill(kinkmothertrack->P(),kinkmothertrack->GetTPCsignal());
642  fHistNTPCclsVsRadius->Fill(R,TPCnclsMother);
643  }
644 }
645 
646 //___________________________________________________________________
648 
649  float par0 = 0.493677;
650  float par1 = 0.9127037;
651  float par2 = TMath::Pi();
652 
653  return TMath::ATan(par0*par1*1./TMath::Sqrt(p*p*(1-par1*par1)-(par0*par0*par1*par1)))*180/par2;
654 }
655 
656 //________________________________________________________________
658 {
659  float t_d = fPIDresp->GetTOFResponse().GetExpectedSignal(track, AliPID::kTriton); //largest mass possible with Z=1
660  float len = track->GetIntegratedLength();
661  float beta_d = len / (t_d * kCSPEED);
662  float mass = AliPID::ParticleMassZ(AliPID::kTriton); //largest mass possible with Z=1
663 
664  if(TMath::Abs(beta_d-1.) < 1.e-12) return track->GetTPCmomentum();
665  else return mass*beta_d/sqrt(1.-(beta_d*beta_d));
666 }
667 
668 //___________________________________________________________________
670  if(num>=static_cast<float>(numeric_limits<short>::max())) return numeric_limits<short>::max();
671  else if(num<=static_cast<float>(numeric_limits<short>::min())) return numeric_limits<short>::min();
672 
673  return static_cast<short>(num);
674 }
675 
676 //___________________________________________________________________
678  if(num>=static_cast<float>(numeric_limits<unsigned short>::max())) return numeric_limits<unsigned short>::max();
679  else if(num<=static_cast<float>(numeric_limits<unsigned short>::min())) return numeric_limits<unsigned short>::min();
680 
681  return static_cast<unsigned short>(num);
682 }
Int_t pdg
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
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)
Double_t mass
unsigned short fTrackLength
number of PID clusters in TPC to fill the tree
double fFracToKeepDownSampling
flag to enable/disable downsampling
unsigned char fTPCNcls
transverse momentum to fill the tree
TH2F * fHistOpeningAngleVsPMotherKink
! histo for opening angle vs. pT mother kink
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
TH2F * fHistNsigmaTOFvsPt[kNHypo]
! array of histos for nsigmaTPC vs pt (MC truth)
virtual void UserExec(Option_t *option)
unsigned short fdEdxTPC
TOF momentum to fill the tree.
unsigned short fToF
TPC dEdX to fill the tree.
double fPtMaxDownSampling
fraction to keep when downsampling activated
float fCentMin
max radius in XY for kinks
float fRMaxKinks
min radius in XY for kinks
AliAODEvent * fAOD
single-track cut set
TH2F * fHistNsigmaTPCvsPt[kNHypo]
! array of histos for nsigmaTPC vs pt (MC truth)
unsigned short fPTPC
Nsigma PID to fill the tree.
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
system: 0->pp,pPb 1->PbPb
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)
short fPDGcode
start time resolution for TOF PID
TH2F * fHistArmenteroPlot[5]
! histo for armenteros-podolanski plot
unsigned char fTag
PDG code in case of MC to fill the tree.
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 ConvertFloatToUnsignedShort(float num)
const char Option_t
Definition: External.C:48
float fNsigmaMaxForTag
bit map for tag (see enum above)
TString fTriggerClass
centrality estimator
short GetPDGcodeFromMC(AliAODTrack *track, TClonesArray *arrayMC)
bool fEnabledDownSampling
flag to enable filling of the tree with only Nsigma variables for the PID
AliAODTrack * IsKinkDaughter(AliAODTrack *track)