AliPhysics  a6017e1 (a6017e1)
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  if(!d) continue;
414  Double_t ptCand=-999.;
415  Int_t iSpec=-1;
416  Int_t labD=-1;
417  Int_t charge=0;
418  Double_t invMass=0.;
419  if(idCase==0){
420  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF2Prong*)d))continue;
421  ptCand=d->Pt();
422  fHistNCand->Fill(0.,ptCand);
423  labD = d->MatchToMC(421,arrayMC,2,pdg0);
424  if(labD>=0){
425  iSpec=0;
426  AliAODMCParticle *pd0 = (AliAODMCParticle*)arrayMC->At(labD);
427  if(pd0->GetPdgCode()==421) invMass=((AliAODRecoDecayHF2Prong*)d)->InvMassD0();
428  else invMass=((AliAODRecoDecayHF2Prong*)d)->InvMassD0bar();
429  }
430  }else if(idCase==1){
431  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF3Prong*)d))continue;
432  ptCand=d->Pt();
433  if(d->HasSelectionBit(AliRDHFCuts::kDplusCuts)) fHistNCand->Fill(1.,ptCand);
434  if(d->HasSelectionBit(AliRDHFCuts::kDsCuts)) fHistNCand->Fill(3.,ptCand);
435  if(d->HasSelectionBit(AliRDHFCuts::kLcCuts)) fHistNCand->Fill(4.,ptCand);
436  labD = d->MatchToMC(411,arrayMC,3,pdgp);
437  if(labD>=0){
438  AliAODMCParticle *pdp = (AliAODMCParticle*)arrayMC->At(labD);
440  invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDplus();
441  if(pdp->GetPdgCode()==411) charge=1;
442  else if(pdp->GetPdgCode()==-411) charge=-1;
443  }else{
444  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
445  AliAODMCParticle* pdau0=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
446  Int_t pdgCode0=TMath::Abs(pdau0->GetPdgCode());
447  labD = d->MatchToMC(431,arrayMC,3,pdgs);
448  if(labD>=0){
449  if(d->HasSelectionBit(AliRDHFCuts::kDsCuts)) iSpec=3;
450  if(pdgCode0==321) invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDsKKpi();
451  else invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDspiKK();
452  AliAODMCParticle *pds = (AliAODMCParticle*)arrayMC->At(labD);
453  if(pds->GetPdgCode()==431) charge=1;
454  else if(pds->GetPdgCode()==-431) charge=-1;
455  }else{
456  labD = d->MatchToMC(4122,arrayMC,3,pdgl);
457  if(labD>=0){
458  if(d->HasSelectionBit(AliRDHFCuts::kLcCuts)) iSpec=4;
459  if(pdgCode0==2212) invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassLcpKpi();
460  else invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassLcpiKp();
461  AliAODMCParticle *plc = (AliAODMCParticle*)arrayMC->At(labD);
462  if(plc->GetPdgCode()==4122) charge=1;
463  else if(plc->GetPdgCode()==-4122) charge=-1;
464  }
465  }
466  }
467  }else if(idCase==2){
468  if(!vHF->FillRecoCasc(aod,((AliAODRecoCascadeHF*)d),kTRUE))continue;
469  ptCand=d->Pt();
470  fHistNCand->Fill(2.,ptCand);
471  labD = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgst,pdg0,arrayMC);
472  if(labD>=0){
473  iSpec=2;
474  invMass=((AliAODRecoCascadeHF*)d)->InvMassDstarKpipi();
475  AliAODMCParticle *pdst = (AliAODMCParticle*)arrayMC->At(labD);
476  if(pdst->GetPdgCode()==413) charge=1;
477  else if(pdst->GetPdgCode()==-413) charge=-1;
478  }
479  }
480 
481  if(labD>=0 && iSpec>=0){
482  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
483  if(!partD) continue;
484  Double_t ptgen=partD->Pt();
485  Double_t phigen=partD->Phi();
486  Double_t ygen=partD->Y();
487  Double_t ptreco=d->Pt();
488  Double_t dlen=d->DecayLength()*10000.; //um
489  Double_t ndlenxy=d->NormalizedDecayLengthXY();
490  Double_t cosp=d->CosPointingAngle();
491  Int_t iDau=partD->GetFirstDaughter();
492  AliAODMCParticle *dauD=0x0;
493  if(iDau>=0) dauD=(AliAODMCParticle*)arrayMC->At(iDau);
494  if(!dauD) continue;
495  Double_t dx=(d->Xv()-dauD->Xv())*10000.;
496  Double_t dy=(d->Yv()-dauD->Yv())*10000.;
497  Double_t dz=(d->Zv()-dauD->Zv())*10000.;
498  Double_t dlentruex=dauD->Xv()-partD->Xv();
499  Double_t dlentruey=dauD->Yv()-partD->Yv();
500  Double_t dlentruez=dauD->Zv()-partD->Zv();
501  // Double_t dlenxytrue=TMath::Sqrt(dlentruex*dlentruex+dlentruey*dlentruey);
502  Double_t dlentrue=TMath::Sqrt(dlentruex*dlentruex+dlentruey*dlentruey+dlentruez*dlentruez);
503  Double_t phidecvert=TMath::Pi()+TMath::ATan2(-dlentruey,-dlentruex);
504  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partD,kTRUE);//Prompt = 4, FeedDown = 5
505  if(orig<4 || orig>5) continue;
506  Int_t indexh=iSpec*2+(orig-4);
507  if(fEnableExtraHistos && iSpec==1 && charge==-1){
508  //D- kept separate
509  indexh=2*5+(orig-4);
510  }
511  fHistPtYMultRecoFilt[indexh]->Fill(ptgen,ygen,evCentr);
512  if(iSpec==1){
513  if(fRDHFCutsDplus){
514  if(fRDHFCutsDplus->IsEventSelected(aod) &&
516  fHistPtYMultRecoSel[indexh]->Fill(ptgen,ygen,evCentr);
517  }
518  }
519  }
520  fHistXvtxResVsPt[indexh]->Fill(ptreco,dx);
521  fHistYvtxResVsPt[indexh]->Fill(ptreco,dy);
522  fHistZvtxResVsPt[indexh]->Fill(ptreco,dz);
523  fHistInvMassVsPt[indexh]->Fill(ptreco,invMass);
524  fHistDecLenVsPt[indexh]->Fill(ptreco,dlen);
525  fHistNormDLxyVsPt[indexh]->Fill(ptreco,ndlenxy);
526  fHistCosPointVsPt[indexh]->Fill(ptreco,cosp);
527  if(fEnableExtraHistos){
528  fHistXvtxResVsPhi[indexh]->Fill(phigen,dx);
529  fHistYvtxResVsPhi[indexh]->Fill(phigen,dy);
530  fHistZvtxResVsPhi[indexh]->Fill(phigen,dz);
531  fHistXvtxResVsDecLenVsPt[indexh]->Fill(ptreco,dlentrue,dx);
532  fHistYvtxResVsDecLenVsPt[indexh]->Fill(ptreco,dlentrue,dy);
533  fHistZvtxResVsDecLenVsPt[indexh]->Fill(ptreco,dlentrue,dz);
534  // Double_t xorigrot=partD->Xv()*TMath::Cos(phidecvert)+partD->Yv()*TMath::Sin(phidecvert);
535  // Double_t yorigrot=-partD->Xv()*TMath::Sin(phidecvert)+partD->Yv()*TMath::Cos(phidecvert);
536  Double_t xdecgenrot=dauD->Xv()*TMath::Cos(phidecvert)+dauD->Yv()*TMath::Sin(phidecvert);
537  Double_t ydecgenrot=-dauD->Xv()*TMath::Sin(phidecvert)+dauD->Yv()*TMath::Cos(phidecvert);
538  Double_t xdecrecrot=d->Xv()*TMath::Cos(phidecvert)+d->Yv()*TMath::Sin(phidecvert);
539  Double_t ydecrecrot=-d->Xv()*TMath::Sin(phidecvert)+d->Yv()*TMath::Cos(phidecvert);
540  Double_t dxrot=(xdecrecrot-xdecgenrot)*10000.;
541  Double_t dyrot=(ydecrecrot-ydecgenrot)*10000.;
542  fHistXvtxResRotVsPt[indexh]->Fill(ptreco,dxrot);
543  fHistYvtxResRotVsPt[indexh]->Fill(ptreco,dyrot);
544  }
545  }
546  }
547  delete vHF;
548  return;
549 }
550 //_________________________________________________________________
551 Bool_t AliAnalysisTaskDmesonMCPerform::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
553  for (Int_t iProng = 0; iProng<nProng; iProng++){
554  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
555  if(!mcPartDaughter) return kFALSE;
556  Double_t eta = mcPartDaughter->Eta();
557  Double_t pt = mcPartDaughter->Pt();
558  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
559  }
560  return kTRUE;
561 }
562 //________________________________________________________________________
565  //
566  Int_t nTracks=aod->GetNumberOfTracks();
567 
568  for(Int_t it=0; it<nTracks; it++) {
569  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
570  if(!tr) continue;
571  if(tr->GetID()<0) continue;
572  Int_t lab=TMath::Abs(tr->GetLabel());
573  if(lab<kMaxLabel) fMapTrLabel[lab]=it;
574  else printf("Label %d exceeds upper limit\n",lab);
575  }
576  return;
577 }
578 //________________________________________________________________________
580  //
582  //
583  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
584  Int_t nFound=0;
585  for(Int_t jd=0; jd<nProng; jd++){
586  Int_t it=fMapTrLabel[labDau[jd]];
587  if(it>0){
588  AliAODTrack* tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
589  if(tr){
590  if(TMath::Abs(tr->GetLabel())!=labDau[jd]){
591  printf("Mismatched track label\n");
592  continue;
593  }
594  px[jd]=tr->Px();
595  py[jd]=tr->Py();
596  px[jd]=tr->Pz();
597  d0[jd]=0.; // temporary
598  d0err[jd]=0.; // temporary
599  nFound++;
600  }
601  }
602  }
603  if(nFound!=nProng) return 0x0;
604  AliAODRecoDecayHF* dmes=new AliAODRecoDecayHF(0x0,nProng,nProng-2,px,py,pz,d0,d0err);
605  return dmes;
606 }
607 //________________________________________________________________________
609 {
611  //
612  if(fDebug > 1) printf("AnalysisTaskDmesonMCPerform: Terminate() \n");
613 
614  fOutput = dynamic_cast<TList*> (GetOutputData(1));
615  if (!fOutput) {
616  printf("ERROR: fOutput not available\n");
617  return;
618  }
619 
620  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
621  if(fHistNEvents){
622  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(10));
623  }else{
624  printf("ERROR: fHistNEvents not available\n");
625  return;
626  }
627 
628  return;
629 }
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.