AliPhysics  9b6b435 (9b6b435)
AliAnalysisTaskSEDmesonPIDSysProp.cxx
Go to the documentation of this file.
1 /* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. */
2 
4 // \class AliAnalysisTaskSEDmesonPIDSysProp //
5 // \brief analysis task for PID Systematic uncertainty propagation from the single track to the D mesons //
6 // \author: A. M. Barbano, anastasia.maria.barbano@cern.ch //
7 // \author: F. Grosa, fabrizio.grosa@cern.ch //
9 
11 #include "AliAnalysisTaskSE.h"
12 #include "AliAnalysisManager.h"
13 #include "AliAODHandler.h"
14 #include "AliLog.h"
15 #include "AliMCEventHandler.h"
16 #include "AliMCEvent.h"
18 #include "AliRDHFCutsDstoKKpi.h"
19 #include "AliRDHFCutsD0toKpi.h"
21 #include "AliAODMCParticle.h"
22 #include "AliAODMCHeader.h"
23 #include "AliAODRecoDecay.h"
26 #include "AliAODRecoCascadeHF.h"
27 #include "AliAnalysisVertexingHF.h"
28 #include "AliVertexingHFUtils.h"
29 #include "AliPIDResponse.h"
30 #include <TFile.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TChain.h>
34 
36 
37 //________________________________________________________________________
39 AliAnalysisTaskSE("taskPIDSysProp"),
40 fOutput(nullptr),
41 fHistNEvents(nullptr),
42 fHistPtDauVsD(nullptr),
43 fHistSystPIDEffD(nullptr),
44 fHistEffPionTOF(nullptr),
45 fHistEffKaonTOF(nullptr),
46 fHistSystPionTOF(nullptr),
47 fHistSystKaonTOF(nullptr),
48 fPartName(""),
49 fPIDresp(nullptr),
50 fSystFileName(""),
51 fPIDstrategy(kConservativePID),
52 fnSigma(3.),
53 fDecayChannel(kD0toKpi),
54 fKaonTPCHistoOpt(kKaonTOFtag),
55 fKaonTOFHistoOpt(kSamePionV0tag),
56 fAODProtection(1),
57 fNPtBins(0),
58 fPtLimits(nullptr),
59 fAnalysisCuts(nullptr)
60 {
61  for(int iHist=0; iHist<2; iHist++) {
62  fHistEffPionTPC[iHist]=nullptr;
63  fHistEffKaonTPC[iHist]=nullptr;
64  fHistSystPionTPC[iHist]=nullptr;
65  fHistSystKaonTPC[iHist]=nullptr;
66  }
67 
68  DefineInput(0, TChain::Class());
69  DefineOutput(1, TList::Class());
70 }
71 
72 //________________________________________________________________________
74 AliAnalysisTaskSE("taskPIDSysProp"),
75 fOutput(nullptr),
76 fHistNEvents(nullptr),
77 fHistPtDauVsD(nullptr),
78 fHistSystPIDEffD(nullptr),
79 fHistEffPionTOF(nullptr),
80 fHistEffKaonTOF(nullptr),
81 fHistSystPionTOF(nullptr),
82 fHistSystKaonTOF(nullptr),
83 fPartName(""),
84 fPIDresp(nullptr),
85 fSystFileName(systfilename),
86 fPIDstrategy(kConservativePID),
87 fnSigma(3.),
88 fDecayChannel(ch),
89 fKaonTPCHistoOpt(kKaonTOFtag),
90 fKaonTOFHistoOpt(kSamePionV0tag),
91 fAODProtection(1),
92 fNPtBins(0),
93 fPtLimits(nullptr),
94 fAnalysisCuts(cuts)
95 {
96  for(int iHist=0; iHist<2; iHist++) {
97  fHistEffPionTPC[iHist]=nullptr;
98  fHistEffKaonTPC[iHist]=nullptr;
99  fHistSystPionTPC[iHist]=nullptr;
100  fHistSystKaonTPC[iHist]=nullptr;
101  }
102 
103  DefineInput(0, TChain::Class());
104  DefineOutput(1, TList::Class());
105 }
106 
107 //___________________________________________________________________________
109 {
110  if(fOutput && !fOutput->IsOwner()){
111  delete fHistNEvents;
112  for(int iHist=0; iHist<2; iHist++) {
113  if(fHistSystPionTPC[iHist]) delete fHistSystPionTPC[iHist];
114  if(fHistSystKaonTPC[iHist]) delete fHistSystKaonTPC[iHist];
115  }
116  delete fHistSystPionTOF;
117  delete fHistSystKaonTOF;
118 
119  delete fHistPtDauVsD;
120  delete fHistSystPIDEffD;
121  delete fOutput;
122  }
123  if(fPIDresp) delete fPIDresp;
124  if(fAnalysisCuts) delete fAnalysisCuts;
125  if(fPtLimits) delete[] fPtLimits;
126  fOutput = nullptr;
127 }
128 
129 //________________________________________________________________________
131 {
132  int load = LoadEffSystFile();
133  if(load>0) AliFatal("Impossible to load single track systematic file, check if it is correct! Exit.");
134 
135  fOutput = new TList();
136  fOutput->SetOwner();
137  fOutput->SetName("OutputHistos");
138 
139  fHistNEvents = new TH1F("hNEvents", "number of events ",15,-0.5,13.5);
140  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsRead");
141  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents Matched dAOD");
142  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents Mismatched dAOD");
143  fHistNEvents->GetXaxis()->SetBinLabel(4,"nEventsAnal");
144  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. passing IsEvSelected");
145  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected due to trigger");
146  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected due to not reco vertex");
147  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for contr vertex");
148  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for vertex out of accept");
149  fHistNEvents->GetXaxis()->SetBinLabel(10,"n. rejected for pileup events");
150  fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of out centrality events");
151  fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of D candidates");
152  fHistNEvents->GetXaxis()->SetBinLabel(13,"no. of D after PID cuts");
153  fHistNEvents->GetXaxis()->SetBinLabel(14,"no. of not on-the-fly rec D");
154 
155  fHistNEvents->GetXaxis()->SetNdivisions(1,false);
156 
157  fHistNEvents->Sumw2();
158  fHistNEvents->SetMinimum(0);
159  fOutput->Add(fHistNEvents);
160 
162  float* ptlims = fAnalysisCuts->GetPtBinLimits();
163  fPtLimits = new double[fNPtBins+1];
164  for(int iPt=0; iPt<=fNPtBins; iPt++) {
165  fPtLimits[iPt] = static_cast<double>(ptlims[iPt]);
166  }
167  if(fPtLimits[fNPtBins]>100) fPtLimits[fNPtBins] = 100.;
168 
169  fHistSystPIDEffD = new TH2F("fHistSystPIDEffD","PID efficiency systematic uncertainty; #it{p}_{T}^{D} (GeV/#it{c}); relative systematic uncertainty",fNPtBins,fPtLimits,150,0.,0.15);
170  fHistPtDauVsD = new TH2F("fHistPtDauVsD","#it{p}_{T} Dau vs #it{p}_{T} D; #it{p}_{T}^{D} (GeV/#it{c}); #it{p}_{T}^{daugh} (GeV/#it{c})",static_cast<int>(fPtLimits[fNPtBins] * 10),0.,fPtLimits[fNPtBins],static_cast<int>(fPtLimits[fNPtBins] * 10),0.,fPtLimits[fNPtBins]);
172  fOutput->Add(fHistPtDauVsD);
173 
174  PostData(1,fOutput);
175 }
176 
177 //________________________________________________________________________
179 {
181  fAnalysisCuts = new AliRDHFCutsDplustoKpipi(*(dynamic_cast<AliRDHFCutsDplustoKpipi*>(fAnalysisCuts)));
182  }
183  else if(fDecayChannel == kDstoKKpi) {
184  fAnalysisCuts = new AliRDHFCutsDstoKKpi(*(dynamic_cast<AliRDHFCutsDstoKKpi*>(fAnalysisCuts)));
185  }
186  else if(fDecayChannel == kD0toKpi) {
187  fAnalysisCuts = new AliRDHFCutsD0toKpi(*(dynamic_cast<AliRDHFCutsD0toKpi*>(fAnalysisCuts)));
188  }
189  else if(fDecayChannel == kDstartoKpipi) {
190  fAnalysisCuts = new AliRDHFCutsDStartoKpipi(*(dynamic_cast<AliRDHFCutsDStartoKpipi*>(fAnalysisCuts)));
191  }
192  else {
193  AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exit.");
194  }
195 }
196 
197 //________________________________________________________________________
199 {
200  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fInputEvent);
201  fHistNEvents->Fill(0); // all events
202 
203  if(fAODProtection>=0){
204  // Protection against different number of events in the AOD and deltaAOD
205  // In case of discrepancy the event is rejected.
206  int matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
207  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
208  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
209  fHistNEvents->Fill(2);
210  PostData(1,fOutput);
211  return;
212  }
213  }
214  TClonesArray *arrayBranch=0;
215 
216  if(!aod && AODEvent() && IsStandardAOD()) {
217  // In case there is an AOD handler writing a standard AOD, use the AOD
218  // event in memory rather than the input (ESD) event.
219  aod = dynamic_cast<AliAODEvent*> (AODEvent());
220  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
221  // have to taken from the AOD event hold by the AliAODExtension
222  AliAODHandler* aodHandler = (AliAODHandler*)
223  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
224  if(aodHandler->GetExtensions()) {
225  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
226  AliAODEvent *aodFromExt = ext->GetAOD();
228  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
229  else if(fDecayChannel == kD0toKpi)
230  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
231  else if(fDecayChannel == kDstartoKpipi)
232  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
233  }
234  }
235  else {
237  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
238  else if(fDecayChannel == kD0toKpi)
239  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
240  else if(fDecayChannel == kDstartoKpipi)
241  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("Dstar");
242  }
243 
244  if (!arrayBranch) {
245  AliError("Could not find array of HF vertices");
246  PostData(1,fOutput);
247  return;
248  }
249 
250  if (!fPIDresp) fPIDresp = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
251 
252  AliAODVertex *aodVtx = (AliAODVertex*)aod->GetPrimaryVertex();
253  if (!aodVtx || TMath::Abs(aod->GetMagneticField())<0.001) {
254  AliDebug(3, "The event was skipped due to missing vertex or magnetic field issue");
255  PostData(1,fOutput);
256  return;
257  }
258  fHistNEvents->Fill(3); // count event
259 
260  bool isEvSel = fAnalysisCuts->IsEventSelected(aod);
261  float ntracks = aod->GetNumberOfTracks();
262 
269 
270  int runNumber = aod->GetRunNumber();
271 
272  TClonesArray *arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
273  AliAODMCHeader *mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
274 
275  if(!arrayMC) {
276  AliError("AliAnalysisTaskSEDmesonPIDSysProp::UserExec: MC particles branch not found!\n");
277  PostData(1,fOutput);
278  return;
279  }
280  if(!mcHeader) {
281  AliError("AliAnalysisTaskSEDmesonPIDSysProp::UserExec: MC header branch not found!\n");
282  PostData(1,fOutput);
283  return;
284  }
285 
286  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677)) {
287  // protection for events with empty trigger mask in p-Pb (Run1)
288  PostData(1,fOutput);
289  return;
290  }
292  // events not passing the centrality selection can be removed immediately.
293  PostData(1,fOutput);
294  return;
295  }
296  double zMCVertex = mcHeader->GetVtxZ();
297  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ()) {
298  PostData(1,fOutput);
299  return;
300  }
301  if(!isEvSel){
302  PostData(1,fOutput);
303  return;
304  }
305  fHistNEvents->Fill(4);
306 
307  int nCand = arrayBranch->GetEntriesFast();
308 
309  int nprongs = -1;
310  int pdgcode = -1;
311  int pdgDaughter[3];
312  int pdg2Daughter[2];
314  pdgcode = 411;
315  nprongs = 3;
316  fPartName="Dplus";
317  pdgDaughter[0]=321;
318  pdgDaughter[1]=211;
319  pdgDaughter[2]=211;
320  }else if(fDecayChannel == kDstoKKpi){
321  pdgcode = 431;
322  nprongs = 3;
323  fPartName="Ds";
324  pdgDaughter[0]=321;
325  pdgDaughter[1]=321;
326  pdgDaughter[2]=211;
327  }else if(fDecayChannel == kD0toKpi){
328  pdgcode = 421;
329  nprongs = 2;
330  fPartName="D0";
331  pdgDaughter[0]=321;
332  pdgDaughter[1]=211;
333  }else if(fDecayChannel == kDstartoKpipi){
334  pdgcode = 413;
335  nprongs = 2;
336  fPartName="Dstar";
337  pdgDaughter[0]=421;
338  pdgDaughter[1]=211;
339  pdg2Daughter[0]=321;
340  pdg2Daughter[1]=211;
341  }else{
342  AliError("Wrong decay setting");
343  PostData(1,fOutput);
344  return;
345  }
346 
347  // vHF object is needed to call the method that refills the missing info of the candidates
348  // if they have been deleted in dAOD reconstruction phase
349  // in order to reduce the size of the file
351 
352  for (int iCand = 0; iCand < nCand; iCand++) {
353 
354  AliAODRecoDecayHF* d = nullptr;
355  bool isDStarCand = false;
356 
358  d = dynamic_cast<AliAODRecoDecayHF3Prong*>(arrayBranch->UncheckedAt(iCand));
359  if(!vHF->FillRecoCand(aod,dynamic_cast<AliAODRecoDecayHF3Prong*>(d))) {
360  fHistNEvents->Fill(13);
361  continue;
362  }
363  }
364  else if(fDecayChannel == kD0toKpi) {
365  d = dynamic_cast<AliAODRecoDecayHF2Prong*>(arrayBranch->UncheckedAt(iCand));
366  if(!vHF->FillRecoCand(aod,dynamic_cast<AliAODRecoDecayHF2Prong*>(d))) {
367  fHistNEvents->Fill(13);
368  continue;
369  }
370  }
371  else if(fDecayChannel == kDstartoKpipi) {
372  d = dynamic_cast<AliAODRecoCascadeHF*>(arrayBranch->UncheckedAt(iCand));
373  if(!d) continue;
374  isDStarCand = true;
375  if(!vHF->FillRecoCasc(aod,dynamic_cast<AliAODRecoCascadeHF*>(d),isDStarCand)) {
376  fHistNEvents->Fill(13);
377  continue;
378  }
379  if(!d->GetSecondaryVtx()) continue;
380  }
381 
382  fHistNEvents->Fill(11);
383  bool unsetvtx=false;
384  if(!d->GetOwnPrimaryVtx()){
385  d->SetOwnPrimaryVtx(aodVtx);
386  unsetvtx=true;
387  }
388 
389  bool recVtx=false;
390  AliAODVertex *origownvtx = nullptr;
391 
392  double ptD = d->Pt();
393  double rapid = d->Y(pdgcode);
394  bool isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(ptD,rapid);
395 
396  if(isFidAcc){
397  int retCodeAnalysisCuts = fAnalysisCuts->IsSelectedPID(d);
398  if(retCodeAnalysisCuts==0) continue;
399  fHistNEvents->Fill(12);
400 
401  int mcLabel=-1;
402  int orig = 0;
403  if(!isDStarCand) mcLabel = d->MatchToMC(pdgcode,arrayMC,nprongs,pdgDaughter);
404  else mcLabel = (dynamic_cast<AliAODRecoCascadeHF*>(d))->MatchToMC(pdgcode,421,pdgDaughter,pdg2Daughter,arrayMC);
405  if(mcLabel<0) continue;
406  AliAODMCParticle* partD = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mcLabel));
407  if(partD) orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partD,kTRUE);
408  if(orig<4) continue;
409 
410  const int nDau = d->GetNDaughters();
411  AliAODTrack* dautrack[nDau];
412  if(fDecayChannel != kDstartoKpipi) {
413  for(int iDau=0; iDau<nDau; iDau++)
414  dautrack[iDau] = dynamic_cast<AliAODTrack*>(d->GetDaughter(iDau));
415  }
416  else {
417  AliAODRecoDecayHF2Prong* D0prong = dynamic_cast<AliAODRecoDecayHF2Prong*>((dynamic_cast<AliAODRecoCascadeHF*>(d))->Get2Prong());
418  if(!D0prong) {
419  continue;
420  }
421  for(int iDau=0; iDau<nDau; iDau++) {
422  dautrack[iDau] = dynamic_cast<AliAODTrack*>(D0prong->GetDaughter(iDau));
423  }
424  }
425 
426  double syst = GetDmesonPIDuncertainty(dautrack,nDau,arrayMC,ptD);
427  if(syst==-999. || syst==0.) continue;
428  fHistSystPIDEffD->Fill(ptD,syst);
429  }
430  if(unsetvtx) d->UnsetOwnPrimaryVtx();
431  }
432 
433  delete vHF;
434  vHF = nullptr;
435 
436  PostData(1,fOutput);
437  return;
438 }
439 
440 //________________________________________________________________________
442 {
443  TFile* infile = TFile::Open(fSystFileName.Data());
444  if(!infile) return 1;
445 
447  fHistEffPionTPC[0] = (TH1F*)infile->Get("hEffPionTPCDataV0tag_3sigma");
448  fHistSystPionTPC[0] = (TH1F*)infile->Get("hRatioEffPionTPCDataV0tag_3sigma");
449  if(!fHistSystPionTPC[0] || !fHistEffPionTPC[0]) return 2;
451  fHistEffKaonTPC[0] = (TH1F*)infile->Get("hEffKaonTPCDataTOFtag_3sigma");
452  fHistSystKaonTPC[0] = (TH1F*)infile->Get("hRatioEffKaonTPCDataTOFtag_3sigma");
453  }
454  else if(fKaonTPCHistoOpt==kKaonKinkstag) {
455  fHistEffKaonTPC[0] = (TH1F*)infile->Get("hEffKaonTPCDataKinktag_3sigma");
456  fHistSystKaonTPC[0] = (TH1F*)infile->Get("hRatioEffKaonTPCDataKinktag_3sigma");
457  }
458  if(!fHistSystKaonTPC[0] || !fHistEffKaonTPC[0]) return 3;
459  fHistEffPionTOF = (TH1F*)infile->Get("hEffPionTOFDataV0tag_3sigma");
460  fHistSystPionTOF = (TH1F*)infile->Get("hRatioEffPionTOFDataV0tag_3sigma");
461  if(!fHistSystPionTOF || !fHistEffPionTOF) return 4;
463  fHistEffKaonTOF = (TH1F*)infile->Get("hEffKaonTOFDataTPCtag_3sigma");
464  fHistSystKaonTOF = (TH1F*)infile->Get("hRatioEffKaonTOFDataTPCtag_3sigma");
465  }
466  else if(fKaonTOFHistoOpt==kSamePionV0tag) {
467  fHistEffKaonTOF = (TH1F*)infile->Get("hEffPionTOFDataV0tag_3sigma");
468  fHistSystKaonTOF = (TH1F*)infile->Get("hRatioEffPionTOFDataV0tag_3sigma");
469  }
470  if(!fHistSystKaonTOF || !fHistEffKaonTOF) return 5;
471  if(fPIDstrategy==kStrongPID) {
472  fHistEffPionTPC[1] = (TH1F*)infile->Get("hEffPionTPCDataV0tag_2sigma");
473  fHistSystPionTPC[1] = (TH1F*)infile->Get("hRatioEffPionTPCDataV0tag_2sigma");
474  if(!fHistSystPionTPC[1] || !fHistEffPionTPC[1]) return 6;
476  fHistEffKaonTPC[1] = (TH1F*)infile->Get("hEffKaonTPCDataTOFtag_2sigma");
477  fHistSystKaonTPC[1] = (TH1F*)infile->Get("hRatioEffKaonTPCDataTOFtag_2sigma");
478  }
479  else if(fKaonTPCHistoOpt==kKaonKinkstag) {
480  fHistEffKaonTPC[1] = (TH1F*)infile->Get("hEffKaonTPCDataKinktag_2sigma");
481  fHistSystKaonTPC[1] = (TH1F*)infile->Get("hRatioEffKaonTPCDataKinktag_2sigma");
482  }
483  if(!fHistSystKaonTPC[1] || !fHistEffKaonTPC[1]) return 7;
484  }
485  }
486  else if(fPIDstrategy==knSigmaPID) {
487  fHistEffPionTPC[0] = (TH1F*)infile->Get(Form("hEffPionTPCDataV0tag_%0.fsigma",fnSigma));
488  fHistSystPionTPC[0] = (TH1F*)infile->Get(Form("hRatioEffPionTPCDataV0tag_%0.fsigma",fnSigma));
489  if(!fHistSystPionTPC[0] || !fHistEffPionTPC[0]) return 8;
491  fHistEffKaonTPC[0] = (TH1F*)infile->Get(Form("hEffKaonTPCDataTOFtag_%0.fsigma",fnSigma));
492  fHistSystKaonTPC[0] = (TH1F*)infile->Get(Form("hRatioEffKaonTPCDataTOFtag_%0.fsigma",fnSigma));
493  }
494  else if(fKaonTPCHistoOpt==kKaonKinkstag) {
495  fHistEffKaonTPC[0] = (TH1F*)infile->Get(Form("hEffKaonTPCDataKinktag_%0.fsigma",fnSigma));
496  fHistSystKaonTPC[0] = (TH1F*)infile->Get(Form("hRatioEffKaonTPCDataKinktag_%0.fsigma",fnSigma));
497  }
498  if(!fHistSystKaonTPC[0] || !fHistEffKaonTPC[0]) return 9;
499  fHistEffPionTOF = (TH1F*)infile->Get(Form("hEffPionTOFDataV0tag_%0.fsigma",fnSigma));
500  fHistSystPionTOF = (TH1F*)infile->Get(Form("hRatioEffPionTOFDataV0tag_%0.fsigma",fnSigma));
501  if(!fHistSystPionTOF || !fHistEffPionTOF) return 10;
503  fHistEffKaonTOF = (TH1F*)infile->Get(Form("hEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
504  fHistSystKaonTOF = (TH1F*)infile->Get(Form("hRatioEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
505  }
506  else if(fKaonTOFHistoOpt==kSamePionV0tag) {
507  fHistEffKaonTOF = (TH1F*)infile->Get(Form("hEffPionTOFDataV0tag_%0.fsigma",fnSigma));
508  fHistSystKaonTOF = (TH1F*)infile->Get(Form("hRatioEffPionTOFDataV0tag_%0.fsigma",fnSigma));
509  }
510  fHistEffKaonTOF = (TH1F*)infile->Get(Form("hEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
511  fHistSystKaonTOF = (TH1F*)infile->Get(Form("hRatioEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
512  if(!fHistSystKaonTOF || !fHistEffKaonTOF) return 11;
513  }
514 
515  return 0;
516 }
517 
518 //________________________________________________________________________
519 double AliAnalysisTaskSEDmesonPIDSysProp::GetDmesonPIDuncertainty(AliAODTrack *track[], const int nDau, TClonesArray *arrayMC, double ptD)
520 {
521  double syst = 0.;
522 
523  for(int iDau=0; iDau<nDau; iDau++){
524  if(!track[iDau]){
525  AliWarning("Daughter-particle track not found");
526  return -999.;
527  }
528  int labDau = track[iDau]->GetLabel();
529  if(labDau<0) {
530  AliWarning("Daughter particle not found");
531  return -999.;
532  }
533  AliAODMCParticle* p = dynamic_cast<AliAODMCParticle*>(arrayMC->UncheckedAt(TMath::Abs(labDau)));
534  if(!p) {
535  AliWarning("Daughter particle not found");
536  return -999.;
537  }
538 
539  int daupdgcode = TMath::Abs(p->GetPdgCode());
540  double daupt = track[iDau]->Pt();
541 
542  bool isTPCok = false;
543  bool isTOFok = false;
544  if(fPIDresp->CheckPIDStatus(AliPIDResponse::kTPC,track[iDau]) == AliPIDResponse::kDetPidOk) isTPCok = true;
545  if(fPIDresp->CheckPIDStatus(AliPIDResponse::kTOF,track[iDau]) == AliPIDResponse::kDetPidOk) isTOFok = true;
546 
547  int bin = fHistSystPionTPC[0]->GetXaxis()->FindBin(daupt);
548  double systTPC=0.;
549  double systTOF=0.;
550  double probTPC=0.;
551  double probTOF=0.;
552  double probTPCandTOF=0.;
553  double probTPCorTOF=0.;
554 
556  if(daupdgcode==211) {
557  if(isTPCok) GetSingleTrackSystAndProb(fHistSystPionTPC[0],fHistEffPionTPC[0],bin,systTPC,probTPC);
558  if(isTOFok) GetSingleTrackSystAndProb(fHistSystPionTOF,fHistEffPionTOF,bin,systTOF,probTOF);
559  }
560  else if(daupdgcode==321) {
561  if(isTPCok) GetSingleTrackSystAndProb(fHistSystKaonTPC[0],fHistEffKaonTPC[0],bin,systTPC,probTPC);
562  if(isTOFok) GetSingleTrackSystAndProb(fHistSystKaonTOF,fHistEffKaonTOF,bin,systTOF,probTOF);
563  }
564  probTPCorTOF = probTPC+probTOF-probTPC*probTOF;
565  if(probTPCorTOF>1.e-20) syst += TMath::Sqrt((1-probTPC)*(1-probTPC)*systTOF*systTOF+(1-probTOF)*(1-probTOF)*systTPC*systTPC)/probTPCorTOF;
566  }
567  else if(fPIDstrategy==kStrongPID) {
568  if(daupdgcode==211) {
569  if(isTOFok) {
570  if(isTPCok) GetSingleTrackSystAndProb(fHistSystPionTPC[0],fHistEffPionTPC[0],bin,systTPC,probTPC);
572  }
573  else if(isTPCok && !isTOFok) {//2sisgma cut in case of no TOF
575  }
576  }
577  else if(daupdgcode==321) {
578  if(isTOFok) {
579  if(isTPCok) GetSingleTrackSystAndProb(fHistSystKaonTPC[0],fHistEffKaonTPC[0],bin,systTPC,probTPC);
581  }
582  else if(isTPCok && !isTOFok) {//2sisgma cut in case of no TOF
584  }
585  }
586  probTPCorTOF = probTPC+probTOF-probTPC*probTOF;
587  if(probTPCorTOF>1.e-20) syst += TMath::Sqrt((1-probTPC)*(1-probTPC)*systTOF*systTOF+(1-probTOF)*(1-probTOF)*systTPC*systTPC)/probTPCorTOF;
588  }
589  else if(fPIDstrategy==knSigmaPID) {
590  if(daupdgcode==211) {
591  if(isTPCok && isTOFok) {
594  }
595  }
596  else if(daupdgcode==321) {
597  if(isTPCok && isTOFok) {
600  }
601  }
602  probTPCandTOF = probTPC*probTOF;
603  if(probTPCandTOF>1.e-20) syst += TMath::Sqrt(probTPC*probTPC*systTOF*systTOF+probTOF*probTOF*systTPC*systTPC)/probTPCandTOF;
604  }
605  fHistPtDauVsD->Fill(ptD,daupt);
606  }
607 
608  return TMath::Abs(syst);
609 }
610 
611 //________________________________________________________________________
612 void AliAnalysisTaskSEDmesonPIDSysProp::GetSingleTrackSystAndProb(TH1F* hSingleTrackSyst, TH1F* hSingleTrackEff, int bin, double &syst, double &prob)
613 {
614  int binmax = hSingleTrackSyst->GetNbinsX();
615  int binmin = 1;
616  if(bin<binmax && bin>binmin) {
617  syst = TMath::Abs(1-hSingleTrackSyst->GetBinContent(bin))*hSingleTrackEff->GetBinContent(bin); //absolute syst
618  prob = hSingleTrackEff->GetBinContent(bin);
619  }
620  else if(bin>=binmax){
621  syst = TMath::Abs(1-hSingleTrackSyst->GetBinContent(binmax))*hSingleTrackEff->GetBinContent(binmax); //absolute syst
622  prob = hSingleTrackEff->GetBinContent(binmax);
623  }
624  else if(bin<=binmin){
625  syst = TMath::Abs(1-hSingleTrackSyst->GetBinContent(binmax))*hSingleTrackEff->GetBinContent(binmin); //absolute syst
626  prob = hSingleTrackEff->GetBinContent(binmin);
627  }
628 }
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:337
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:331
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:322
Definition: External.C:236
Int_t IsEventSelectedInCentrality(AliVEvent *event)
static Int_t CheckMatchingAODdeltaAODevents()
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:325
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:258
double fnSigma
PID strategy (conservative, strong, nsigma..)
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *)
Definition: AliRDHFCuts.h:301
int fAODProtection
option for syst on kaon TOF PID efficiency
Class for cuts on AOD reconstructed D+->Kpipi.
int fNPtBins
flag to activate protection against AOD-dAOD mismatch.
int fPIDstrategy
Name of file with single track syst. unc.
AliPIDResponse * fPIDresp
string for particle name
AliAODVertex * GetOwnPrimaryVtx() const
const double ptlims[]
int fKaonTOFHistoOpt
option for syst on kaon TPC PID efficiency
int fDecayChannel
number of sigma PID if nsigma strategy enabled
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:334
void GetSingleTrackSystAndProb(TH1F *hSingleTrackSyst, TH1F *hSingleTrackEff, int bin, double &syst, double &prob)
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
TH2F * fHistPtDauVsD
! histo with pT daughters vs pT candidate
Bool_t IsEventSelected(AliVEvent *event)
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:249
TH1F * fHistNEvents
! histo with number of events
TH2F * fHistSystPIDEffD
! histo with PID systematic uncertainty on the D candidate
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:250
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:319
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:314
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:281
double GetDmesonPIDuncertainty(AliAODTrack *track[], const int nDau, TClonesArray *arrayMC, double ptD)