AliPhysics  cc5d4fb (cc5d4fb)
AliAnalysisTaskDmesonMCPerform.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2018, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 //*************************************************************************
17 // Class AliAnalysisTaskDmesonMCPerform
18 // AliAnalysisTaskSE for performance studies of
19 // D meson hadronic decays in MC simulations
20 // F. Prino, prino@to.infn.it
21 //*************************************************************************
22 
23 #include <TList.h>
24 
25 #include "AliAnalysisManager.h"
26 #include "AliInputEventHandler.h"
27 #include "AliAODHandler.h"
28 #include "AliAODEvent.h"
29 #include "AliAODMCParticle.h"
30 #include "AliAODMCHeader.h"
31 #include "AliAODVertex.h"
32 #include "AliAnalysisVertexingHF.h"
33 #include "AliVertexingHFUtils.h"
34 #include "AliAODRecoDecayHF.h"
37 
41 
42 //________________________________________________________________________
44  AliAnalysisTaskSE("taskDperf"),
45  fOutput(0x0),
46  fHistNEvents(0x0),
47  fHistNGenD(0x0),
48  fHistNCand(0x0),
49  fNPtBins(16),
50  fMinPt(0.),
51  fMaxPt(16.),
52  fEnableExtraHistos(kFALSE),
53  fAODProtection(1),
54  fRDHFCuts(0x0),
55  fRDHFCutsDplus(0x0)
56 {
57  //
59  //
60 
61  fRDHFCuts=new AliRDHFCutsD0toKpi("EvSelCuts");
65  fRDHFCuts->SetMaxVtxZ(10.);
68 
69  fPartName[0]="Dzero";
70  fPartName[1]="Dplus";
71  fPartName[2]="Dstar";
72  fPartName[3]="Ds";
73  fPartName[4]="Lc2pkpi";
74  fPartName[5]="Dminus";
75 
76  // Output slot #1 writes into a TList container
77  DefineOutput(1,TList::Class()); //My private output
78 }
79 
80 //________________________________________________________________________
82 {
83  //
85  //
86  if(fOutput && !fOutput->IsOwner()){
87  // delete histograms stored in fOutput
88  delete fHistNEvents;
89  delete fHistNGenD;
90  delete fHistNCand;
91  for(Int_t j=0; j<2*kDecays; j++){
92  delete fHistPtYMultGen[j];
93  delete fHistPtYMultGenDauInAcc[j];
94  delete fHistPtYMultReco[j];
95  delete fHistPtYMultRecoFilt[j];
96  delete fHistPtYMultRecoSel[j];
97  delete fHistXvtxResVsPt[j];
98  delete fHistYvtxResVsPt[j];
99  delete fHistZvtxResVsPt[j];
100  delete fHistXvtxResRotVsPt[j];
101  delete fHistYvtxResRotVsPt[j];
102  delete fHistXvtxResVsPhi[j];
103  delete fHistYvtxResVsPhi[j];
104  delete fHistZvtxResVsPhi[j];
105  delete fHistXvtxResVsDecLenVsPt[j];
106  delete fHistYvtxResVsDecLenVsPt[j];
107  delete fHistZvtxResVsDecLenVsPt[j];
108  delete fHistInvMassVsPt[j];
109  delete fHistDecLenVsPt[j];
110  delete fHistNormDLxyVsPt[j];
111  delete fHistCosPointVsPt[j];
112  }
113  }
114  delete fOutput;
115  delete fRDHFCuts;
116 }
117 //________________________________________________________________________
119 {
120  //
122  //
123  if(fDebug > 1) printf("AliAnalysisTaskDmesonMCPerform::UserCreateOutputObjects() \n");
124  fOutput = new TList();
125  fOutput->SetOwner();
126  fOutput->SetName("OutputHistos");
127 
128  fHistNEvents = new TH1F("fHistNEvents", "number of events ",10,-0.5,9.5);
129  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents read");
130  fHistNEvents->GetXaxis()->SetBinLabel(2,"Rejected due to mismatch in trees");
131  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with good AOD");
132  fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
133  fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to vertex reco");
134  fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to pileup");
135  fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to centrality");
136  fHistNEvents->GetXaxis()->SetBinLabel(8,"Rejected due to vtxz");
137  fHistNEvents->GetXaxis()->SetBinLabel(9,"Rejected due to Physics Sel");
138  fHistNEvents->GetXaxis()->SetBinLabel(10,"nEvents accepted");
139  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
140  fHistNEvents->SetMinimum(0);
141  fOutput->Add(fHistNEvents);
142 
143  fHistNGenD=new TH1F("fHistNGenD","number of D mesons",10,-0.5,9.5);
144  TString names[5]={"Dzero","Dplus","Dstar","Ds","Lc"};
145  TString type[2]={"Prompt","Feeddown"};
146  for(Int_t j=0; j<5; j++){
147  for(Int_t i=0; i<2; i++){
148  Int_t index=j*2+i;
149  fHistNGenD->GetXaxis()->SetBinLabel(index+1,Form("%s %s",type[i].Data(),names[j].Data()));
150  }
151  }
152  fHistNGenD->GetXaxis()->SetNdivisions(1,kFALSE);
153  fHistNGenD->SetMinimum(0);
154  fOutput->Add(fHistNGenD);
155 
156  fHistNCand=new TH2F("fHistNCand","number of D meson candidates",5,-0.5,4.5,fNPtBins,fMinPt,fMaxPt);
157  fHistNCand->GetXaxis()->SetBinLabel(1,"D0#rightarrowK#pi");
158  fHistNCand->GetXaxis()->SetBinLabel(2,"D+#rightarrowK#pi#pi");
159  fHistNCand->GetXaxis()->SetBinLabel(3,"D*+#rightarrowD0#pi");
160  fHistNCand->GetXaxis()->SetBinLabel(4,"D_s#rightarrowKK#pi");
161  fHistNCand->GetXaxis()->SetBinLabel(5,"Lc#rightarrowpK#pi");
162  fHistNCand->SetMinimum(0);
163  fOutput->Add(fHistNCand);
164 
165  for(Int_t j=0; j<kDecays; j++){
166  for(Int_t i=0; i<2; i++){
167  Int_t index=j*2+i;
168  fHistPtYMultGen[index]=new TH3F(Form("hPtYMult%sGen%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",fNPtBins,fMinPt,fMaxPt,24,-6.,6.,40,0.,100.);
169  fHistPtYMultGenDauInAcc[index]=new TH3F(Form("hPtYMult%sGenDauInAcc%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",fNPtBins,fMinPt,fMaxPt,24,-6.,6.,40,0.,100.);
170  fHistPtYMultReco[index]=new TH3F(Form("hPtYMult%sReco%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",fNPtBins,fMinPt,fMaxPt,24,-6.,6.,40,0.,100.);
171  fHistPtYMultRecoFilt[index]=new TH3F(Form("hPtYMult%sRecoFilt%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",fNPtBins,fMinPt,fMaxPt,24,-6.,6.,40,0.,100.);
172  fHistPtYMultRecoSel[index]=new TH3F(Form("hPtYMult%sRecoSel%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",fNPtBins,fMinPt,fMaxPt,24,-6.,6.,40,0.,100.);
173  fOutput->Add(fHistPtYMultGen[index]);
174  fOutput->Add(fHistPtYMultGenDauInAcc[index]);
175  fOutput->Add(fHistPtYMultReco[index]);
176  fOutput->Add(fHistPtYMultRecoFilt[index]);
177  fOutput->Add(fHistPtYMultRecoSel[index]);
178  fHistXvtxResVsPt[index]=new TH2F(Form("hXvtxResVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; x_{v}(rec)-x_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,-500.,500.);
179  fHistYvtxResVsPt[index]=new TH2F(Form("hYvtxResVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y_{v}(rec)-y_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,-500.,500.);
180  fHistZvtxResVsPt[index]=new TH2F(Form("hZvtxResVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; z_{v}(rec)-z_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,-500.,500.);
181  fHistInvMassVsPt[index]=new TH2F(Form("hInvMassVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; Inv. Mass (GeV/c^{2})",fNPtBins,fMinPt,fMaxPt,300,1.75,2.35);
182  fHistDecLenVsPt[index]=new TH2F(Form("hDecLenVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; dec. len. (#mum)",fNPtBins,fMinPt,fMaxPt,100,0.,5000.);
183  fHistNormDLxyVsPt[index]=new TH2F(Form("hNormDLxyVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; dec. len. (#mum)",fNPtBins,fMinPt,fMaxPt,100,0.,20.);
184  fHistCosPointVsPt[index]=new TH2F(Form("hCosPointVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; cos(#vartheta_{point})",fNPtBins,fMinPt,fMaxPt,100,-1.,1.);
185  fOutput->Add(fHistXvtxResVsPt[index]);
186  fOutput->Add(fHistYvtxResVsPt[index]);
187  fOutput->Add(fHistZvtxResVsPt[index]);
188  fOutput->Add(fHistInvMassVsPt[index]);
189  fOutput->Add(fHistDecLenVsPt[index]);
190  fOutput->Add(fHistNormDLxyVsPt[index]);
191  fOutput->Add(fHistCosPointVsPt[index]);
192 
193  if(fEnableExtraHistos){
194  fHistXvtxResRotVsPt[index]=new TH2F(Form("hXvtxResRotVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; x'_{v}(rec)-x'_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,-500.,500.);
195  fHistYvtxResRotVsPt[index]=new TH2F(Form("hYvtxResRotVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y'_{v}(rec)-y'_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,-500.,500.);
196  fHistXvtxResVsPhi[index]=new TH2F(Form("hXvtxResVsPhi%s%s",type[i].Data(),fPartName[j].Data())," ; #varphi (rad) ; x_{v}(rec)-x_{v}(gen) (#mum)",36,0.,2*TMath::Pi(),100,-500.,500.);
197  fHistYvtxResVsPhi[index]=new TH2F(Form("hYvtxResVsPhi%s%s",type[i].Data(),fPartName[j].Data())," ; #varphi (rad) ; y_{v}(rec)-y_{v}(gen) (#mum)",36,0.,2*TMath::Pi(),100,-500.,500.);
198  fHistZvtxResVsPhi[index]=new TH2F(Form("hZvtxResVsPhi%s%s",type[i].Data(),fPartName[j].Data())," ; #varphi (rad) ; z_{v}(rec)-z_{v}(gen) (#mum)",36,0.,2*TMath::Pi(),100,-500.,500.);
199  fOutput->Add(fHistXvtxResRotVsPt[index]);
200  fOutput->Add(fHistYvtxResRotVsPt[index]);
201  fOutput->Add(fHistXvtxResVsPhi[index]);
202  fOutput->Add(fHistYvtxResVsPhi[index]);
203  fOutput->Add(fHistZvtxResVsPhi[index]);
204  fHistXvtxResVsDecLenVsPt[index]=new TH3F(Form("hXvtxResVsDecLenVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; dec. len. (#mum) ; x_{v}(rec)-x_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,0.,5000.,100,-500.,500.);
205  fHistYvtxResVsDecLenVsPt[index]=new TH3F(Form("hYvtxResVsDecLenVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; dec. len. (#mum) ; y_{v}(rec)-y_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,0.,5000.,100,-500.,500.);
206  fHistZvtxResVsDecLenVsPt[index]=new TH3F(Form("hZvtxResVsDecLenVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; dec. len. (#mum) ; z_{v}(rec)-z_{v}(gen) (#mum)",fNPtBins,fMinPt,fMaxPt,100,0.,5000.,100,-500.,500.);
207  fOutput->Add(fHistXvtxResVsDecLenVsPt[index]);
208  fOutput->Add(fHistYvtxResVsDecLenVsPt[index]);
209  fOutput->Add(fHistZvtxResVsDecLenVsPt[index]);
210  }
211  }
212  }
213 
214 
215  PostData(1,fOutput);
216  return;
217 }
218 //________________________________________________________________________
220 {
223 
224  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
225 
226  fHistNEvents->Fill(0); // count event
227 
228  if(fAODProtection>=0){
229  // Protection against different number of events in the AOD and deltaAOD
230  // In case of discrepancy the event is rejected.
231  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
232  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
233  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
234  fHistNEvents->Fill(1);
235  return;
236  }
237  }
238 
239  // Load all the branches of the DeltaAOD - needed for SelectionBit counting
240  TClonesArray *arrayD0toKpi =0;
241  TClonesArray *array3Prong =0;
242  TClonesArray *arrayDstar =0;
243  TClonesArray *arrayCascades =0;
244  if(!aod && AODEvent() && IsStandardAOD()) {
245  // In case there is an AOD handler writing a standard AOD, use the AOD
246  // event in memory rather than the input (ESD) event.
247  aod = dynamic_cast<AliAODEvent*> (AODEvent());
248  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
249  // have to taken from the AOD event hold by the AliAODExtension
250  AliAODHandler* aodHandler = (AliAODHandler*)
251  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
252  if(aodHandler->GetExtensions()) {
253 
254  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
255  AliAODEvent *aodFromExt = ext->GetAOD();
256 
257  array3Prong =(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
258  arrayD0toKpi =(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
259  arrayDstar =(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
260  arrayCascades=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
261  }
262  } else if(aod) {
263  array3Prong =(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
264  arrayD0toKpi =(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
265  arrayDstar =(TClonesArray*)aod->GetList()->FindObject("Dstar");
266  arrayCascades=(TClonesArray*)aod->GetList()->FindObject("CascadesHF");
267  }
268 
269  if(!aod || !array3Prong || !arrayD0toKpi) {
270  printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
271  return;
272  }
273  // fix for temporary bug in ESDfilter
274  // the AODs with null vertex pointer didn't pass the PhysSel
275  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
276  fHistNEvents->Fill(2);
277 
278  //Event selection
279  Bool_t isEvSel=fRDHFCuts->IsEventSelected(aod);
280  if(fRDHFCuts->GetWhyRejection()==5) fHistNEvents->Fill(3);
281  if(!isEvSel && fRDHFCuts->GetWhyRejection()==0) fHistNEvents->Fill(4);
282  if(fRDHFCuts->GetWhyRejection()==1) fHistNEvents->Fill(5);
283  if(fRDHFCuts->GetWhyRejection()==2) fHistNEvents->Fill(6);
284  if(fRDHFCuts->GetWhyRejection()==6)fHistNEvents->Fill(7);
285  if(fRDHFCuts->GetWhyRejection()==7)fHistNEvents->Fill(8);
286 
287  // load MC header
288  AliAODMCHeader *mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
289  if(!mcHeader) {
290  printf("AliAnalysisTaskDmesonMCPerform:UserExec: MC header branch not found!\n");
291  return;
292  }
293 
294  // load MC particles
295  TClonesArray *arrayMC= (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
296  if(!arrayMC) {
297  printf("AliAnalysisTaskDmesonMCPerform::UserExec: MC particles branch not found!\n");
298  return;
299  }
300 
301  Int_t runNumber=aod->GetRunNumber();
302 
303  if(aod->GetTriggerMask()==0 &&
304  (runNumber>=195344 && runNumber<=195677)){
305  // protection for events with empty trigger mask in p-Pb
306  return;
307  }
309 
310  PostData(1,fOutput);
311 
312  for(Int_t i=0; i<kMaxLabel; i++) fMapTrLabel[i]=-999;
313  MapTrackLabels(aod);
314  FillGenLevelHistos(aod,arrayMC,mcHeader);
315  if(!isEvSel)return;
316  fHistNEvents->Fill(9);
317 
318  FillCandLevelHistos(0,aod,arrayD0toKpi,arrayMC);
319  FillCandLevelHistos(1,aod,array3Prong,arrayMC);
320  FillCandLevelHistos(2,aod,arrayDstar,arrayMC);
321 
322  return;
323 }
324 //________________________________________________________________________
325 void AliAnalysisTaskDmesonMCPerform::FillGenLevelHistos(AliAODEvent* aod, TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
326  //
328 
329  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
330  if(TMath::Abs(zMCVertex)>fRDHFCuts->GetMaxVtxZ()) return;
331  Double_t evCentr=1.;
332  if(fRDHFCuts->GetUseCentrality()>0) evCentr=fRDHFCuts->GetCentrality(aod);
333 
334  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
335  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
336  Int_t pdgCode=TMath::Abs(mcPart->GetPdgCode());
337  Int_t iSpec=-1;
338  if(pdgCode == 411) iSpec=1;
339  else if(pdgCode == 413) iSpec=2;
340  else if(pdgCode == 421) iSpec=0;
341  else if(pdgCode == 431) iSpec=3;
342  else if(pdgCode == 4122) iSpec=4;
343  if (iSpec>=0){
344  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
345  if(orig<4 || orig>5) continue;
346  Int_t indexh=iSpec*2+(orig-4);
347  fHistNGenD->Fill(indexh);
348  Int_t deca=0;
349  Bool_t isGoodDecay=kFALSE;
350  Int_t labDau[4]={-1,-1,-1,-1};
351  Int_t nProng=3;
352  if (pdgCode == 411){
353  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,mcPart,labDau);
354  if(deca>0) isGoodDecay=kTRUE;
355  }else if(pdgCode == 421){
356  deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,mcPart,labDau);
357  if(mcPart->GetNDaughters()!=2) continue;
358  if(deca==1) isGoodDecay=kTRUE;
359  nProng=2;
360  }else if(pdgCode == 431){
361  deca=AliVertexingHFUtils::CheckDsDecay(arrayMC,mcPart,labDau);
362  if(deca==1) isGoodDecay=kTRUE;
363  }else if(pdgCode==413){
364  deca=AliVertexingHFUtils::CheckDstarDecay(arrayMC,mcPart,labDau);
365  if(deca==1) isGoodDecay=kTRUE;
366  }else if(pdgCode==4122){
367  deca=AliVertexingHFUtils::CheckLcpKpiDecay(arrayMC,mcPart,labDau);
368  if(deca>=1) isGoodDecay=kTRUE;
369  }
370  if(labDau[0]==-1){
371  continue; //protection against unfilled array of labels
372  }
373  if(isGoodDecay && indexh>=0){
374  Double_t ptgen=mcPart->Pt();
375  Double_t ygen=mcPart->Y();
376  fHistPtYMultGen[indexh]->Fill(ptgen,ygen,evCentr);
377  if(fRDHFCuts->IsInFiducialAcceptance(ptgen,ygen)){
378  Bool_t dauInAcc=CheckAcceptance(arrayMC,nProng,labDau);
379  if(dauInAcc){
380  fHistPtYMultGenDauInAcc[indexh]->Fill(ptgen,ygen,evCentr);
381  AliAODRecoDecayHF* dmes=GetRecoDecay(aod,nProng,labDau);
382  if(dmes){
383  fHistPtYMultReco[indexh]->Fill(ptgen,ygen,evCentr);
384  }
385  delete dmes;
386  }
387  }
388  }
389  }
390  }
391 }
392 //________________________________________________________________________
393 void AliAnalysisTaskDmesonMCPerform::FillCandLevelHistos(Int_t idCase, AliAODEvent* aod, TClonesArray *arrayDcand, TClonesArray *arrayMC){
394  //
396 
397  Double_t evCentr=1.;
398  if(fRDHFCuts->GetUseCentrality()>0) evCentr=fRDHFCuts->GetCentrality(aod);
399 
400  Int_t pdg0[2]={321,211};
401  Int_t pdgp[3]={321,211,211};
402  Int_t pdgs[3]={321,211,321};
403  Int_t pdgst[2]={421,211};
404  Int_t pdgl[3]={2212,321,211};
405  // vHF object is needed to call the method that refills the missing info of the candidates
406  // if they have been deleted in dAOD reconstruction phase
407  // in order to reduce the size of the file
409 
410  Int_t nCand=arrayDcand->GetEntriesFast();
411  for (Int_t iCand = 0; iCand < nCand; iCand++) {
412  AliAODRecoDecayHF *d=(AliAODRecoDecayHF*)arrayDcand->UncheckedAt(iCand);
413  Double_t ptCand=-999.;
414  Int_t iSpec=-1;
415  Int_t labD=-1;
416  Int_t charge=0;
417  Double_t invMass=0.;
418  if(idCase==0){
419  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF2Prong*)d))continue;
420  ptCand=d->Pt();
421  fHistNCand->Fill(0.,ptCand);
422  labD = d->MatchToMC(421,arrayMC,2,pdg0);
423  if(labD>=0){
424  iSpec=0;
425  AliAODMCParticle *pd0 = (AliAODMCParticle*)arrayMC->At(labD);
426  if(pd0->GetPdgCode()==421) invMass=((AliAODRecoDecayHF2Prong*)d)->InvMassD0();
427  else invMass=((AliAODRecoDecayHF2Prong*)d)->InvMassD0bar();
428  }
429  }else if(idCase==1){
430  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF3Prong*)d))continue;
431  ptCand=d->Pt();
432  if(d->HasSelectionBit(AliRDHFCuts::kDplusCuts)) fHistNCand->Fill(1.,ptCand);
433  if(d->HasSelectionBit(AliRDHFCuts::kDsCuts)) fHistNCand->Fill(3.,ptCand);
434  if(d->HasSelectionBit(AliRDHFCuts::kLcCuts)) fHistNCand->Fill(4.,ptCand);
435  labD = d->MatchToMC(411,arrayMC,3,pdgp);
436  if(labD>=0){
437  AliAODMCParticle *pdp = (AliAODMCParticle*)arrayMC->At(labD);
439  invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDplus();
440  if(pdp->GetPdgCode()==411) charge=1;
441  else if(pdp->GetPdgCode()==-411) charge=-1;
442  }else{
443  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
444  AliAODMCParticle* pdau0=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
445  Int_t pdgCode0=TMath::Abs(pdau0->GetPdgCode());
446  labD = d->MatchToMC(431,arrayMC,3,pdgs);
447  if(labD>=0){
448  if(d->HasSelectionBit(AliRDHFCuts::kDsCuts)) iSpec=3;
449  if(pdgCode0==321) invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDsKKpi();
450  else invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDspiKK();
451  AliAODMCParticle *pds = (AliAODMCParticle*)arrayMC->At(labD);
452  if(pds->GetPdgCode()==431) charge=1;
453  else if(pds->GetPdgCode()==-431) charge=-1;
454  }else{
455  labD = d->MatchToMC(4122,arrayMC,3,pdgl);
456  if(labD>=0){
457  if(d->HasSelectionBit(AliRDHFCuts::kLcCuts)) iSpec=4;
458  if(pdgCode0==2212) invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassLcpKpi();
459  else invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassLcpiKp();
460  AliAODMCParticle *plc = (AliAODMCParticle*)arrayMC->At(labD);
461  if(plc->GetPdgCode()==4122) charge=1;
462  else if(plc->GetPdgCode()==-4122) charge=-1;
463  }
464  }
465  }
466  }else if(idCase==2){
467  if(!vHF->FillRecoCasc(aod,((AliAODRecoCascadeHF*)d),kTRUE))continue;
468  ptCand=d->Pt();
469  fHistNCand->Fill(2.,ptCand);
470  labD = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgst,pdg0,arrayMC);
471  if(labD>=0){
472  iSpec=2;
473  invMass=((AliAODRecoCascadeHF*)d)->InvMassDstarKpipi();
474  AliAODMCParticle *pdst = (AliAODMCParticle*)arrayMC->At(labD);
475  if(pdst->GetPdgCode()==413) charge=1;
476  else if(pdst->GetPdgCode()==-413) charge=-1;
477  }
478  }
479 
480  if(labD>=0 && iSpec>=0){
481  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
482  Double_t ptgen=partD->Pt();
483  Double_t phigen=partD->Phi();
484  Double_t ygen=partD->Y();
485  Double_t ptreco=d->Pt();
486  Double_t dlen=d->DecayLength()*10000.; //um
487  Double_t ndlenxy=d->NormalizedDecayLengthXY();
488  Double_t cosp=d->CosPointingAngle();
489  Double_t dx=(d->Xv()-partD->Xv())*10000.;
490  Double_t dy=(d->Yv()-partD->Yv())*10000.;
491  Double_t dz=(d->Zv()-partD->Zv())*10000.;
492  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partD,kTRUE);//Prompt = 4, FeedDown = 5
493  if(orig<4 || orig>5) continue;
494  Int_t indexh=iSpec*2+(orig-4);
495  if(iSpec==1 && charge==-1){
496  //D- kept separate
497  indexh=2*5+(orig-4);
498  }
499  fHistPtYMultRecoFilt[indexh]->Fill(ptgen,ygen,evCentr);
500  if(iSpec==1){
501  if(fRDHFCutsDplus){
502  if(fRDHFCutsDplus->IsEventSelected(aod) &&
504  fHistPtYMultRecoSel[indexh]->Fill(ptgen,ygen,evCentr);
505  }
506  }
507  }
508  fHistXvtxResVsPt[indexh]->Fill(ptreco,dx);
509  fHistYvtxResVsPt[indexh]->Fill(ptreco,dy);
510  fHistZvtxResVsPt[indexh]->Fill(ptreco,dz);
511  fHistInvMassVsPt[indexh]->Fill(ptreco,invMass);
512  fHistDecLenVsPt[indexh]->Fill(ptreco,dlen);
513  fHistNormDLxyVsPt[indexh]->Fill(ptreco,ndlenxy);
514  fHistCosPointVsPt[indexh]->Fill(ptreco,cosp);
515  if(fEnableExtraHistos){
516  fHistXvtxResVsPhi[indexh]->Fill(phigen,dx);
517  fHistYvtxResVsPhi[indexh]->Fill(phigen,dy);
518  fHistZvtxResVsPhi[indexh]->Fill(phigen,dz);
519  fHistXvtxResVsDecLenVsPt[indexh]->Fill(ptreco,dlen,dx);
520  fHistYvtxResVsDecLenVsPt[indexh]->Fill(ptreco,dlen,dy);
521  fHistZvtxResVsDecLenVsPt[indexh]->Fill(ptreco,dlen,dz);
522  Double_t xgenrot=partD->Xv()*TMath::Cos(phigen)-partD->Yv()*TMath::Sin(phigen);
523  Double_t ygenrot=partD->Xv()*TMath::Sin(phigen)+partD->Yv()*TMath::Cos(phigen);
524  Double_t xrecrot=d->Xv()*TMath::Cos(phigen)-d->Yv()*TMath::Sin(phigen);
525  Double_t yrecrot=d->Xv()*TMath::Sin(phigen)+d->Yv()*TMath::Cos(phigen);
526  Double_t dxrot=(xrecrot-xgenrot)*10000.;
527  Double_t dyrot=(yrecrot-ygenrot)*10000.;
528  fHistXvtxResRotVsPt[indexh]->Fill(ptreco,dxrot);
529  fHistYvtxResRotVsPt[indexh]->Fill(ptreco,dyrot);
530  }
531  }
532  }
533  delete vHF;
534  return;
535 }
536 //_________________________________________________________________
537 Bool_t AliAnalysisTaskDmesonMCPerform::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
539  for (Int_t iProng = 0; iProng<nProng; iProng++){
540  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
541  if(!mcPartDaughter) return kFALSE;
542  Double_t eta = mcPartDaughter->Eta();
543  Double_t pt = mcPartDaughter->Pt();
544  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
545  }
546  return kTRUE;
547 }
548 //________________________________________________________________________
551  //
552  Int_t nTracks=aod->GetNumberOfTracks();
553 
554  for(Int_t it=0; it<nTracks; it++) {
555  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
556  if(!tr) continue;
557  if(tr->GetID()<0) continue;
558  Int_t lab=TMath::Abs(tr->GetLabel());
559  if(lab<kMaxLabel) fMapTrLabel[lab]=it;
560  else printf("Label %d exceeds upper limit\n",lab);
561  }
562  return;
563 }
564 //________________________________________________________________________
566  //
568  //
569  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
570  Int_t nFound=0;
571  for(Int_t jd=0; jd<nProng; jd++){
572  Int_t it=fMapTrLabel[labDau[jd]];
573  if(it>0){
574  AliAODTrack* tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
575  if(tr){
576  if(TMath::Abs(tr->GetLabel())!=labDau[jd]){
577  printf("Mismatched track label\n");
578  continue;
579  }
580  px[jd]=tr->Px();
581  py[jd]=tr->Py();
582  px[jd]=tr->Pz();
583  d0[jd]=0.; // temporary
584  d0err[jd]=0.; // temporary
585  nFound++;
586  }
587  }
588  }
589  if(nFound!=nProng) return 0x0;
590  AliAODRecoDecayHF* dmes=new AliAODRecoDecayHF(0x0,nProng,nProng-2,px,py,pz,d0,d0err);
591  return dmes;
592 }
593 //________________________________________________________________________
595 {
597  //
598  if(fDebug > 1) printf("AnalysisTaskDmesonMCPerform: Terminate() \n");
599 
600  fOutput = dynamic_cast<TList*> (GetOutputData(1));
601  if (!fOutput) {
602  printf("ERROR: fOutput not available\n");
603  return;
604  }
605 
606  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
607  if(fHistNEvents){
608  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(10));
609  }else{
610  printf("ERROR: fHistNEvents not available\n");
611  return;
612  }
613 
614  return;
615 }
Int_t charge
Double_t NormalizedDecayLengthXY() const
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH2F * fHistYvtxResRotVsPt[2 *kDecays]
! hist. for sec vert x resol
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
static Int_t CheckDplusDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH3F * fHistZvtxResVsDecLenVsPt[2 *kDecays]
! hist. for sec vert x resol
void SetMaxVtxZ(Float_t z=1e6)
Definition: AliRDHFCuts.h:61
Int_t IsEventSelectedInCentrality(AliVEvent *event)
Bool_t HasSelectionBit(Int_t i) const
AliAODRecoDecayHF * GetRecoDecay(AliAODEvent *aod, Int_t nProng, Int_t *labDau)
static Int_t CheckDsDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
TH3F * fHistYvtxResVsDecLenVsPt[2 *kDecays]
! hist. for sec vert x resol
static Int_t CheckLcpKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
void SetUseCentrality(Int_t flag=1)
static Int_t CheckMatchingAODdeltaAODevents()
TH2F * fHistXvtxResRotVsPt[2 *kDecays]
! hist. for sec vert x resol
TH2F * fHistXvtxResVsPt[2 *kDecays]
! hist. for sec vert x resol
TH3F * fHistPtYMultRecoFilt[2 *kDecays]
! hist. for D candidates
TH3F * fHistPtYMultRecoSel[2 *kDecays]
! hist. for D candidates (sel)
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:315
void SetUseAnyTrigger()
Definition: AliRDHFCuts.h:72
TH2F * fHistDecLenVsPt[2 *kDecays]
! hist. of decay length (meas)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
TH3F * fHistPtYMultGen[2 *kDecays]
! hist. for generated D
Int_t fAODProtection
flag to control the additional info
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
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:257
void SetUsePhysicsSelection(Bool_t use=kTRUE)
Definition: AliRDHFCuts.h:362
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
AliRDHFCutsD0toKpi * fRDHFCuts
flag to activate protection against AOD-dAOD mismatch.
TString fPartName[kDecays]
Cuts for Dplus.
AliRDHFCutsDplustoKpipi * fRDHFCutsDplus
Cuts for event selection.
TH2F * fHistNormDLxyVsPt[2 *kDecays]
! hist. of decay length (meas)
int Int_t
Definition: External.C:63
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
TH2F * fHistNCand
! hist. for N. of D candidates
TH2F * fHistZvtxResVsPhi[2 *kDecays]
! hist. for sec vert x resol
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:258
static Int_t CheckDstarDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
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)
Double_t fMinPt
number of pt bins in histos
TH2F * fHistInvMassVsPt[2 *kDecays]
! hist. of inv mass (meas)
Bool_t IsEventSelected(AliVEvent *event)
void FillCandLevelHistos(Int_t idCase, AliAODEvent *aod, TClonesArray *arrayDcand, TClonesArray *arrayMC)
TH3F * fHistPtYMultGenDauInAcc[2 *kDecays]
! hist. for generated D in acc
void FillGenLevelHistos(AliAODEvent *aod, TClonesArray *arrayMC, AliAODMCHeader *mcHeader)
TH1F * fHistNEvents
! hist. for N. of events
TH3F * fHistPtYMultReco[2 *kDecays]
! hist. for D with reco daught.
TH2F * fHistYvtxResVsPt[2 *kDecays]
! hist. for sec vert x resol
const char Option_t
Definition: External.C:48
TH2F * fHistCosPointVsPt[2 *kDecays]
! hist. of cos(theta_p) (meas)
TH1F * fHistNGenD
! hist. for N. of D&#39;s
TH2F * fHistYvtxResVsPhi[2 *kDecays]
! hist. for sec vert x resol
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
void SetTriggerClass(TString trclass0, TString trclass1="")
Definition: AliRDHFCuts.h:197
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:280
TH2F * fHistZvtxResVsPt[2 *kDecays]
! hist. for sec vert x resol
void SetOptPileup(Int_t opt=0)
Definition: AliRDHFCuts.h:223
Double_t DecayLength() const
TH2F * fHistXvtxResVsPhi[2 *kDecays]
! hist. for sec vert x resol
Bool_t CheckAcceptance(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
TH3F * fHistXvtxResVsDecLenVsPt[2 *kDecays]
! hist. for sec vert x resol
TList * fOutput
! list send on output slot 0
Int_t fMapTrLabel[kMaxLabel]
Names of hadron species.