AliPhysics  master (3d17d9d)
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 fPIDstrategy(kConservativePID),
51 fnSigma(3.),
52 fDecayChannel(kD0toKpi),
53 fKaonTPCHistoOpt(kKaonTOFtag),
54 fKaonTOFHistoOpt(kSamePionV0tag),
55 fAODProtection(1),
56 fNPtBins(0),
57 fPtLimits(nullptr),
58 fAnalysisCuts(nullptr),
59 fVarForProp(kPt)
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 fPIDstrategy(kConservativePID),
86 fnSigma(3.),
87 fDecayChannel(ch),
88 fKaonTPCHistoOpt(kKaonTOFtag),
89 fKaonTOFHistoOpt(kSamePionV0tag),
90 fAODProtection(1),
91 fNPtBins(0),
92 fPtLimits(nullptr),
93 fAnalysisCuts(cuts),
94 fVarForProp(kPt)
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  }
122  if(fPIDresp) delete fPIDresp;
123  if(fAnalysisCuts) delete fAnalysisCuts;
124  if(fPtLimits) delete[] fPtLimits;
125  delete fOutput;
126 }
127 
128 //________________________________________________________________________
130 {
131  fOutput = new TList();
132  fOutput->SetOwner();
133  fOutput->SetName("OutputHistos");
134 
135  fHistNEvents = new TH1F("hNEvents", "number of events ",15,-0.5,13.5);
136  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsRead");
137  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents Matched dAOD");
138  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents Mismatched dAOD");
139  fHistNEvents->GetXaxis()->SetBinLabel(4,"nEventsAnal");
140  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. passing IsEvSelected");
141  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected due to trigger");
142  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected due to not reco vertex");
143  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for contr vertex");
144  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for vertex out of accept");
145  fHistNEvents->GetXaxis()->SetBinLabel(10,"n. rejected for pileup events");
146  fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of out centrality events");
147  fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of D candidates");
148  fHistNEvents->GetXaxis()->SetBinLabel(13,"no. of D after PID cuts");
149  fHistNEvents->GetXaxis()->SetBinLabel(14,"no. of not on-the-fly rec D");
150 
151  fHistNEvents->GetXaxis()->SetNdivisions(1,false);
152 
153  fHistNEvents->Sumw2();
154  fHistNEvents->SetMinimum(0);
155  fOutput->Add(fHistNEvents);
156 
158  float* ptlims = fAnalysisCuts->GetPtBinLimits();
159  fPtLimits = new double[fNPtBins+1];
160  for(int iPt=0; iPt<=fNPtBins; iPt++) {
161  fPtLimits[iPt] = static_cast<double>(ptlims[iPt]);
162  }
163  if(fPtLimits[fNPtBins]>100) fPtLimits[fNPtBins] = 100.;
164 
165  TString varname = "";
166  if(fVarForProp==kPt)
167  varname = "#it{p}_{T}";
168  else if(fVarForProp==kP)
169  varname = "#it{p}";
170 
171  fHistSystPIDEffD = new TH2F("fHistSystPIDEffD","PID efficiency systematic uncertainty; #it{p}_{T}^{D} (GeV/#it{c}); relative systematic uncertainty",fNPtBins,fPtLimits,500,0.,0.5);
172  fHistPtDauVsD = new TH2F("fHistPtDauVsD",Form("%s Dau vs #it{p}_{T} D; #it{p}_{T}^{D} (GeV/#it{c}); %s^{daugh} (GeV/#it{c})",varname.Data(),varname.Data()),static_cast<int>(fPtLimits[fNPtBins] * 10),0.,fPtLimits[fNPtBins],static_cast<int>(fPtLimits[fNPtBins] * 10),0.,fPtLimits[fNPtBins]);
174  fOutput->Add(fHistPtDauVsD);
175 
176  PostData(1,fOutput);
177 }
178 
179 //________________________________________________________________________
181 {
183  fAnalysisCuts = new AliRDHFCutsDplustoKpipi(*(dynamic_cast<AliRDHFCutsDplustoKpipi*>(fAnalysisCuts)));
184  }
185  else if(fDecayChannel == kDstoKKpi) {
186  fAnalysisCuts = new AliRDHFCutsDstoKKpi(*(dynamic_cast<AliRDHFCutsDstoKKpi*>(fAnalysisCuts)));
187  }
188  else if(fDecayChannel == kD0toKpi) {
189  fAnalysisCuts = new AliRDHFCutsD0toKpi(*(dynamic_cast<AliRDHFCutsD0toKpi*>(fAnalysisCuts)));
190  }
191  else if(fDecayChannel == kDstartoKpipi) {
192  fAnalysisCuts = new AliRDHFCutsDStartoKpipi(*(dynamic_cast<AliRDHFCutsDStartoKpipi*>(fAnalysisCuts)));
193  }
194  else {
195  AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exit.");
196  }
197 }
198 
199 //________________________________________________________________________
201 {
202  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fInputEvent);
203  fHistNEvents->Fill(0); // all events
204 
205  if(fAODProtection>=0){
206  // Protection against different number of events in the AOD and deltaAOD
207  // In case of discrepancy the event is rejected.
208  int matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
209  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
210  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
211  fHistNEvents->Fill(2);
212  PostData(1,fOutput);
213  return;
214  }
215  }
216  TClonesArray *arrayBranch = nullptr, *arrayD0toKpi = nullptr;
217 
218  if(!aod && AODEvent() && IsStandardAOD()) {
219  // In case there is an AOD handler writing a standard AOD, use the AOD
220  // event in memory rather than the input (ESD) event.
221  aod = dynamic_cast<AliAODEvent*> (AODEvent());
222  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
223  // have to taken from the AOD event hold by the AliAODExtension
224  AliAODHandler* aodHandler = (AliAODHandler*)
225  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
226  if(aodHandler->GetExtensions()) {
227  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
228  AliAODEvent *aodFromExt = ext->GetAOD();
230  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
231  else if(fDecayChannel == kD0toKpi)
232  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
233  else if(fDecayChannel == kDstartoKpipi) {
234  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
235  arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
236  }
237  }
238  }
239  else {
241  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
242  else if(fDecayChannel == kD0toKpi)
243  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
244  else if(fDecayChannel == kDstartoKpipi) {
245  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("Dstar");
246  arrayD0toKpi=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
247  }
248  }
249 
250  if (!arrayBranch) {
251  AliError("Could not find array of HF vertices");
252  PostData(1,fOutput);
253  return;
254  }
255 
256  if (!fPIDresp) fPIDresp = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();
257 
258  AliAODVertex *aodVtx = (AliAODVertex*)aod->GetPrimaryVertex();
259  if (!aodVtx || TMath::Abs(aod->GetMagneticField())<0.001) {
260  AliDebug(3, "The event was skipped due to missing vertex or magnetic field issue");
261  PostData(1,fOutput);
262  return;
263  }
264  fHistNEvents->Fill(3); // count event
265 
266  bool isEvSel = fAnalysisCuts->IsEventSelected(aod);
267 
274 
275  int runNumber = aod->GetRunNumber();
276 
277  TClonesArray *arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
278  AliAODMCHeader *mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
279 
280  if(!arrayMC) {
281  AliError("AliAnalysisTaskSEDmesonPIDSysProp::UserExec: MC particles branch not found!\n");
282  PostData(1,fOutput);
283  return;
284  }
285  if(!mcHeader) {
286  AliError("AliAnalysisTaskSEDmesonPIDSysProp::UserExec: MC header branch not found!\n");
287  PostData(1,fOutput);
288  return;
289  }
290 
291  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677)) {
292  // protection for events with empty trigger mask in p-Pb (Run1)
293  PostData(1,fOutput);
294  return;
295  }
297  // events not passing the centrality selection can be removed immediately.
298  PostData(1,fOutput);
299  return;
300  }
301  double zMCVertex = mcHeader->GetVtxZ();
302  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ()) {
303  PostData(1,fOutput);
304  return;
305  }
306  if(!isEvSel){
307  PostData(1,fOutput);
308  return;
309  }
310  fHistNEvents->Fill(4);
311 
312  int nCand = arrayBranch->GetEntriesFast();
313 
314  int nprongs = -1;
315  int pdgcode = -1;
316  int pdgDaughter[3];
317  int pdg2Daughter[2];
319  pdgcode = 411;
320  nprongs = 3;
321  fPartName="Dplus";
322  pdgDaughter[0]=321;
323  pdgDaughter[1]=211;
324  pdgDaughter[2]=211;
325  }else if(fDecayChannel == kDstoKKpi){
326  pdgcode = 431;
327  nprongs = 3;
328  fPartName="Ds";
329  pdgDaughter[0]=321;
330  pdgDaughter[1]=321;
331  pdgDaughter[2]=211;
332  }else if(fDecayChannel == kD0toKpi){
333  pdgcode = 421;
334  nprongs = 2;
335  fPartName="D0";
336  pdgDaughter[0]=321;
337  pdgDaughter[1]=211;
338  }else if(fDecayChannel == kDstartoKpipi){
339  pdgcode = 413;
340  nprongs = 2;
341  fPartName="Dstar";
342  pdgDaughter[0]=421;
343  pdgDaughter[1]=211;
344  pdg2Daughter[0]=321;
345  pdg2Daughter[1]=211;
346  }else{
347  AliError("Wrong decay setting");
348  PostData(1,fOutput);
349  return;
350  }
351 
352  // vHF object is needed to call the method that refills the missing info of the candidates
353  // if they have been deleted in dAOD reconstruction phase
354  // in order to reduce the size of the file
356 
357  for (int iCand = 0; iCand < nCand; iCand++) {
358 
359  AliAODRecoDecayHF* d = nullptr;
360  AliAODRecoDecayHF2Prong* dD0 = nullptr;
361  bool isDStarCand = false;
362  int nDau = 0;
363 
365  d = dynamic_cast<AliAODRecoDecayHF3Prong*>(arrayBranch->UncheckedAt(iCand));
366  nDau = 3;
367  }
368  else if(fDecayChannel == kD0toKpi) {
369  d = dynamic_cast<AliAODRecoDecayHF2Prong*>(arrayBranch->UncheckedAt(iCand));
370  nDau = 2;
371  }
372  else if(fDecayChannel == kDstartoKpipi) {
373  d = dynamic_cast<AliAODRecoCascadeHF*>(arrayBranch->UncheckedAt(iCand));
374  isDStarCand = true;
375  nDau = 3;
376 
377  if(d && d->GetIsFilled()<1)
378  dD0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(arrayD0toKpi->At(d->GetProngID(1)));
379  else
380  dD0 = (dynamic_cast<AliAODRecoCascadeHF*>(d))->Get2Prong();
381  if(!dD0)
382  continue;
383  }
384 
385  if(!d) continue;
386 
387  //Preselection to speed up task
388  TObjArray arrDauTracks(nDau);
389  AliAODTrack *track = nullptr;
391  for(int iDau=0; iDau<nDau; iDau++){
392  AliAODTrack *track = vHF->GetProng(aod,d,iDau);
393  arrDauTracks.AddAt(track,iDau);
394  }
395  }
396  else {
397  for(int iDau=0; iDau<nDau; iDau++){
398  if(iDau == 0)
399  track=vHF->GetProng(aod,d,iDau); //soft pion
400  else
401  track=vHF->GetProng(aod,dD0,iDau-1); //D0 daughters
402  arrDauTracks.AddAt(track,iDau);
403  }
404  }
405  if(!fAnalysisCuts->PreSelect(arrDauTracks)){
406  continue;
407  }
408 
410  if(!vHF->FillRecoCand(aod,dynamic_cast<AliAODRecoDecayHF3Prong*>(d))) {
411  fHistNEvents->Fill(13);
412  continue;
413  }
414  }
415  else if(fDecayChannel == kD0toKpi) {
416  if(!vHF->FillRecoCand(aod,dynamic_cast<AliAODRecoDecayHF2Prong*>(d))) {
417  fHistNEvents->Fill(13);
418  continue;
419  }
420  }
421  else if(fDecayChannel == kDstartoKpipi) {
422  if(!vHF->FillRecoCasc(aod,dynamic_cast<AliAODRecoCascadeHF*>(d),isDStarCand)) {
423  fHistNEvents->Fill(13);
424  continue;
425  }
426  if(!d->GetSecondaryVtx()) continue;
427  }
428 
429  fHistNEvents->Fill(11);
430  bool unsetvtx=false;
431  if(!d->GetOwnPrimaryVtx()){
432  d->SetOwnPrimaryVtx(aodVtx);
433  unsetvtx=true;
434  }
435 
436  bool recVtx=false;
437  AliAODVertex *origownvtx = nullptr;
439  if(d->GetOwnPrimaryVtx())
440  origownvtx = new AliAODVertex(*d->GetOwnPrimaryVtx());
442  recVtx = true;
443  else fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
444  }
445 
446  double ptD = d->Pt();
447  double rapid = d->Y(pdgcode);
448  bool isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(ptD,rapid);
449 
450  if(isFidAcc){
451  int retCodeAnalysisCuts = fAnalysisCuts->IsSelectedPID(d);
452  int retCodeAnalysisTrackCuts = fAnalysisCuts->IsSelected(d,AliRDHFCuts::kTracks,aod); //reject also
453  if(retCodeAnalysisCuts==0 || retCodeAnalysisTrackCuts==0) {
454  if(unsetvtx) d->UnsetOwnPrimaryVtx();
455  continue;
456  }
457  fHistNEvents->Fill(12);
458 
459  int mcLabel=-1;
460  int orig = 0;
461  if(!isDStarCand) mcLabel = d->MatchToMC(pdgcode,arrayMC,nprongs,pdgDaughter);
462  else mcLabel = (dynamic_cast<AliAODRecoCascadeHF*>(d))->MatchToMC(pdgcode,421,pdgDaughter,pdg2Daughter,arrayMC);
463  if(mcLabel<0) {
464  if(unsetvtx) d->UnsetOwnPrimaryVtx();
465  continue;
466  }
467  AliAODMCParticle* partD = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mcLabel));
468  if(partD) orig = AliVertexingHFUtils::CheckOrigin(arrayMC,partD,true);
469  if(orig<4) {
470  if(unsetvtx) d->UnsetOwnPrimaryVtx();
471  continue;
472  }
473 
474  const int nDau = d->GetNDaughters();
475  AliAODTrack* dautrack[nDau];
476  if(fDecayChannel != kDstartoKpipi) {
477  for(int iDau=0; iDau<nDau; iDau++)
478  dautrack[iDau] = dynamic_cast<AliAODTrack*>(d->GetDaughter(iDau));
479  }
480  else {
481  AliAODRecoDecayHF2Prong* D0prong = dynamic_cast<AliAODRecoDecayHF2Prong*>((dynamic_cast<AliAODRecoCascadeHF*>(d))->Get2Prong());
482  if(!D0prong) {
483  if(unsetvtx) d->UnsetOwnPrimaryVtx();
484  continue;
485  }
486  for(int iDau=0; iDau<nDau; iDau++) {
487  dautrack[iDau] = dynamic_cast<AliAODTrack*>(D0prong->GetDaughter(iDau));
488  }
489  }
490 
491  double syst = GetDmesonPIDuncertainty(dautrack,nDau,arrayMC,ptD);
492  if(syst==-999. || syst==0.) {
493  if(unsetvtx) d->UnsetOwnPrimaryVtx();
494  continue;
495  }
496  fHistSystPIDEffD->Fill(ptD,syst);
497  }
498 
499  if(recVtx) fAnalysisCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
500  if(unsetvtx) d->UnsetOwnPrimaryVtx();
501  }
502 
503  delete vHF;
504  vHF = nullptr;
505 
506  PostData(1,fOutput);
507  return;
508 }
509 
510 //________________________________________________________________________
512 {
514  fHistEffPionTPC[0] = (TH1F*)systfile->Get("hEffPionTPCDataV0tag_3sigma");
515  fHistSystPionTPC[0] = (TH1F*)systfile->Get("hRatioEffPionTPCDataV0tag_3sigma");
516  if(!fHistSystPionTPC[0] || !fHistEffPionTPC[0]) return false;
518  fHistEffKaonTPC[0] = (TH1F*)systfile->Get("hEffKaonTPCDataTOFtag_3sigma");
519  fHistSystKaonTPC[0] = (TH1F*)systfile->Get("hRatioEffKaonTPCDataTOFtag_3sigma");
520  }
521  else if(fKaonTPCHistoOpt==kKaonKinkstag) {
522  fHistEffKaonTPC[0] = (TH1F*)systfile->Get("hEffKaonTPCDataKinktag_3sigma");
523  fHistSystKaonTPC[0] = (TH1F*)systfile->Get("hRatioEffKaonTPCDataKinktag_3sigma");
524  }
525  if(!fHistSystKaonTPC[0] || !fHistEffKaonTPC[0]) return false;
526  fHistEffPionTOF = (TH1F*)systfile->Get("hEffPionTOFDataV0tag_3sigma");
527  fHistSystPionTOF = (TH1F*)systfile->Get("hRatioEffPionTOFDataV0tag_3sigma");
528  if(!fHistSystPionTOF || !fHistEffPionTOF) return false;
530  fHistEffKaonTOF = (TH1F*)systfile->Get("hEffKaonTOFDataTPCtag_3sigma");
531  fHistSystKaonTOF = (TH1F*)systfile->Get("hRatioEffKaonTOFDataTPCtag_3sigma");
532  }
533  else if(fKaonTOFHistoOpt==kSamePionV0tag) {
534  fHistEffKaonTOF = (TH1F*)systfile->Get("hEffPionTOFDataV0tag_3sigma");
535  fHistSystKaonTOF = (TH1F*)systfile->Get("hRatioEffPionTOFDataV0tag_3sigma");
536  }
537  if(!fHistSystKaonTOF || !fHistEffKaonTOF) return false;
538  if(fPIDstrategy==kStrongPID) {
539  fHistEffPionTPC[1] = (TH1F*)systfile->Get("hEffPionTPCDataV0tag_2sigma");
540  fHistSystPionTPC[1] = (TH1F*)systfile->Get("hRatioEffPionTPCDataV0tag_2sigma");
541  if(!fHistSystPionTPC[1] || !fHistEffPionTPC[1]) return false;
543  fHistEffKaonTPC[1] = (TH1F*)systfile->Get("hEffKaonTPCDataTOFtag_2sigma");
544  fHistSystKaonTPC[1] = (TH1F*)systfile->Get("hRatioEffKaonTPCDataTOFtag_2sigma");
545  }
546  else if(fKaonTPCHistoOpt==kKaonKinkstag) {
547  fHistEffKaonTPC[1] = (TH1F*)systfile->Get("hEffKaonTPCDataKinktag_2sigma");
548  fHistSystKaonTPC[1] = (TH1F*)systfile->Get("hRatioEffKaonTPCDataKinktag_2sigma");
549  }
550  if(!fHistSystKaonTPC[1] || !fHistEffKaonTPC[1]) return false;
551  }
552  }
553  else if(fPIDstrategy==knSigmaPID) {
554  fHistEffPionTPC[0] = (TH1F*)systfile->Get(Form("hEffPionTPCDataV0tag_%0.fsigma",fnSigma));
555  fHistSystPionTPC[0] = (TH1F*)systfile->Get(Form("hRatioEffPionTPCDataV0tag_%0.fsigma",fnSigma));
556  if(!fHistSystPionTPC[0] || !fHistEffPionTPC[0]) return false;
558  fHistEffKaonTPC[0] = (TH1F*)systfile->Get(Form("hEffKaonTPCDataTOFtag_%0.fsigma",fnSigma));
559  fHistSystKaonTPC[0] = (TH1F*)systfile->Get(Form("hRatioEffKaonTPCDataTOFtag_%0.fsigma",fnSigma));
560  }
561  else if(fKaonTPCHistoOpt==kKaonKinkstag) {
562  fHistEffKaonTPC[0] = (TH1F*)systfile->Get(Form("hEffKaonTPCDataKinktag_%0.fsigma",fnSigma));
563  fHistSystKaonTPC[0] = (TH1F*)systfile->Get(Form("hRatioEffKaonTPCDataKinktag_%0.fsigma",fnSigma));
564  }
565  if(!fHistSystKaonTPC[0] || !fHistEffKaonTPC[0]) return false;
566  fHistEffPionTOF = (TH1F*)systfile->Get(Form("hEffPionTOFDataV0tag_%0.fsigma",fnSigma));
567  fHistSystPionTOF = (TH1F*)systfile->Get(Form("hRatioEffPionTOFDataV0tag_%0.fsigma",fnSigma));
568  if(!fHistSystPionTOF || !fHistEffPionTOF) return false;
570  fHistEffKaonTOF = (TH1F*)systfile->Get(Form("hEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
571  fHistSystKaonTOF = (TH1F*)systfile->Get(Form("hRatioEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
572  }
573  else if(fKaonTOFHistoOpt==kSamePionV0tag) {
574  fHistEffKaonTOF = (TH1F*)systfile->Get(Form("hEffPionTOFDataV0tag_%0.fsigma",fnSigma));
575  fHistSystKaonTOF = (TH1F*)systfile->Get(Form("hRatioEffPionTOFDataV0tag_%0.fsigma",fnSigma));
576  }
577  fHistEffKaonTOF = (TH1F*)systfile->Get(Form("hEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
578  fHistSystKaonTOF = (TH1F*)systfile->Get(Form("hRatioEffKaonTOFDataTPCtag_%0.fsigma",fnSigma));
579  if(!fHistSystKaonTOF || !fHistEffKaonTOF) return false;
580  }
581 
582  return true;
583 }
584 
585 //________________________________________________________________________
586 double AliAnalysisTaskSEDmesonPIDSysProp::GetDmesonPIDuncertainty(AliAODTrack *track[], const int nDau, TClonesArray *arrayMC, double ptD)
587 {
588  double syst = 0.;
589 
590  for(int iDau=0; iDau<nDau; iDau++){
591  if(!track[iDau]){
592  AliWarning("Daughter-particle track not found");
593  return -999.;
594  }
595  int labDau = track[iDau]->GetLabel();
596  if(labDau<0) {
597  AliWarning("Daughter particle not found");
598  return -999.;
599  }
600  AliAODMCParticle* p = dynamic_cast<AliAODMCParticle*>(arrayMC->UncheckedAt(TMath::Abs(labDau)));
601  if(!p) {
602  AliWarning("Daughter particle not found");
603  return -999.;
604  }
605 
606  int daupdgcode = TMath::Abs(p->GetPdgCode());
607  double daupt = track[iDau]->Pt();
608  double daupTPC = track[iDau]->GetTPCmomentum();
609  double dauvar = -1.;
610  if(fVarForProp==kPt)
611  dauvar = daupt;
612  else if(fVarForProp==kP)
613  dauvar = daupTPC;
614 
615  bool isTPCok = false;
616  bool isTOFok = false;
617  if(fPIDresp->CheckPIDStatus(AliPIDResponse::kTPC,track[iDau]) == AliPIDResponse::kDetPidOk) isTPCok = true;
618  if(fPIDresp->CheckPIDStatus(AliPIDResponse::kTOF,track[iDau]) == AliPIDResponse::kDetPidOk) isTOFok = true;
619 
620  int bin = fHistSystPionTPC[0]->GetXaxis()->FindBin(dauvar);
621  double systTPC=0.;
622  double systTOF=0.;
623  double probTPC=0.;
624  double probTOF=0.;
625  double probTPCandTOF=0.;
626  double probTPCorTOF=0.;
627 
629  if(daupdgcode==211) {
630  if(isTPCok) GetSingleTrackSystAndProb(fHistSystPionTPC[0],fHistEffPionTPC[0],bin,systTPC,probTPC);
631  if(isTOFok) GetSingleTrackSystAndProb(fHistSystPionTOF,fHistEffPionTOF,bin,systTOF,probTOF);
632  }
633  else if(daupdgcode==321) {
634  if(isTPCok) GetSingleTrackSystAndProb(fHistSystKaonTPC[0],fHistEffKaonTPC[0],bin,systTPC,probTPC);
635  if(isTOFok) GetSingleTrackSystAndProb(fHistSystKaonTOF,fHistEffKaonTOF,bin,systTOF,probTOF);
636  }
637  probTPCorTOF = probTPC+probTOF-probTPC*probTOF;
638  if(probTPCorTOF>1.e-20) syst += TMath::Sqrt((1-probTPC)*(1-probTPC)*systTOF*systTOF+(1-probTOF)*(1-probTOF)*systTPC*systTPC)/probTPCorTOF;
639  }
640  else if(fPIDstrategy==kStrongPID) {
641  if(daupdgcode==211) {
642  if(isTOFok) {
643  if(isTPCok) GetSingleTrackSystAndProb(fHistSystPionTPC[0],fHistEffPionTPC[0],bin,systTPC,probTPC);
645  }
646  else if(isTPCok && !isTOFok) {//2sisgma cut in case of no TOF
648  }
649  }
650  else if(daupdgcode==321) {
651  if(isTOFok) {
652  if(isTPCok) GetSingleTrackSystAndProb(fHistSystKaonTPC[0],fHistEffKaonTPC[0],bin,systTPC,probTPC);
654  }
655  else if(isTPCok && !isTOFok) {//2sisgma cut in case of no TOF
657  }
658  }
659  probTPCorTOF = probTPC+probTOF-probTPC*probTOF;
660  if(probTPCorTOF>1.e-20) syst += TMath::Sqrt((1-probTPC)*(1-probTPC)*systTOF*systTOF+(1-probTOF)*(1-probTOF)*systTPC*systTPC)/probTPCorTOF;
661  }
662  else if(fPIDstrategy==knSigmaPID) {
663  if(daupdgcode==211) {
664  if(isTPCok && isTOFok) {
667  }
668  }
669  else if(daupdgcode==321) {
670  if(isTPCok && isTOFok) {
673  }
674  }
675  probTPCandTOF = probTPC*probTOF;
676  if(probTPCandTOF>1.e-20) syst += TMath::Sqrt(probTPC*probTPC*systTOF*systTOF+probTOF*probTOF*systTPC*systTPC)/probTPCandTOF;
677  }
678  fHistPtDauVsD->Fill(ptD,dauvar);
679  }
680 
681  return TMath::Abs(syst);
682 }
683 
684 //________________________________________________________________________
685 void AliAnalysisTaskSEDmesonPIDSysProp::GetSingleTrackSystAndProb(TH1F* hSingleTrackSyst, TH1F* hSingleTrackEff, int bin, double &syst, double &prob)
686 {
687  int binmax = hSingleTrackSyst->GetNbinsX();
688  int binmin = 1;
689  if(bin<binmax && bin>binmin) {
690  syst = TMath::Abs(1-hSingleTrackSyst->GetBinContent(bin))*hSingleTrackEff->GetBinContent(bin); //absolute syst
691  prob = hSingleTrackEff->GetBinContent(bin);
692  }
693  else if(bin>=binmax){
694  syst = TMath::Abs(1-hSingleTrackSyst->GetBinContent(binmax))*hSingleTrackEff->GetBinContent(binmax); //absolute syst
695  prob = hSingleTrackEff->GetBinContent(binmax);
696  }
697  else if(bin<=binmin){
698  syst = TMath::Abs(1-hSingleTrackSyst->GetBinContent(binmax))*hSingleTrackEff->GetBinContent(binmin); //absolute syst
699  prob = hSingleTrackEff->GetBinContent(binmin);
700  }
701 }
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:361
AliAODTrack * GetProng(AliVEvent *event, AliAODRecoDecayHF *rd, Int_t iprong)
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:355
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:346
virtual Int_t PreSelect(TObjArray)
Definition: AliRDHFCuts.h:314
Definition: External.C:236
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Double_t ptlims[nPtBins+1]
static Int_t CheckMatchingAODdeltaAODevents()
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:349
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:274
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:321
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_t GetIsFilled() const
AliPIDResponse * fPIDresp
string for particle name
AliAODVertex * GetOwnPrimaryVtx() const
int fKaonTOFHistoOpt
option for syst on kaon TPC PID efficiency
int fDecayChannel
number of sigma PID if nsigma strategy enabled
UShort_t GetProngID(Int_t ip) const
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:358
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 GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:295
Bool_t IsEventSelected(AliVEvent *event)
void CleanOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod, AliAODVertex *origownvtx) const
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:265
TH1F * fHistNEvents
! histo with number of events
TH2F * fHistSystPIDEffD
! histo with PID systematic uncertainty on the D candidate
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:312
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:266
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:343
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:338
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:297
Bool_t RecalcOwnPrimaryVtx(AliAODRecoDecayHF *d, AliAODEvent *aod) const
double GetDmesonPIDuncertainty(AliAODTrack *track[], const int nDau, TClonesArray *arrayMC, double ptD)