AliPhysics  b555aef (b555aef)
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  fAODProtection(1),
50  fRDHFCuts(0x0),
51  fRDHFCutsDplus(0x0)
52 {
53  //
55  //
56 
57  fRDHFCuts=new AliRDHFCutsD0toKpi("EvSelCuts");
61  fRDHFCuts->SetMaxVtxZ(10.);
64 
65  fPartName[0]="Dzero";
66  fPartName[1]="Dplus";
67  fPartName[2]="Dstar";
68  fPartName[3]="Ds";
69  fPartName[4]="Lc2pkpi";
70 
71  // Output slot #1 writes into a TList container
72  DefineOutput(1,TList::Class()); //My private output
73 }
74 
75 //________________________________________________________________________
77 {
78  //
80  //
81  if(fOutput && !fOutput->IsOwner()){
82  // delete histograms stored in fOutput
83  delete fHistNEvents;
84  delete fHistNGenD;
85  delete fHistNCand;
86  for(Int_t j=0; j<2*kDecays; j++){
87  delete fHistPtYMultGen[j];
88  delete fHistPtYMultGenDauInAcc[j];
89  delete fHistPtYMultReco[j];
90  delete fHistPtYMultRecoFilt[j];
91  delete fHistPtYMultRecoSel[j];
92  delete fHistXvtxResVsPt[j];
93  delete fHistYvtxResVsPt[j];
94  delete fHistZvtxResVsPt[j];
95  delete fHistInvMassVsPt[j];
96  delete fHistDecLenVsPt[j];
97  delete fHistNormDLxyVsPt[j];
98  delete fHistCosPointVsPt[j];
99  }
100  }
101  delete fOutput;
102  delete fRDHFCuts;
103 }
104 //________________________________________________________________________
106 {
107  //
109  //
110  if(fDebug > 1) printf("AliAnalysisTaskDmesonMCPerform::UserCreateOutputObjects() \n");
111  fOutput = new TList();
112  fOutput->SetOwner();
113  fOutput->SetName("OutputHistos");
114 
115  fHistNEvents = new TH1F("fHistNEvents", "number of events ",10,-0.5,9.5);
116  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents read");
117  fHistNEvents->GetXaxis()->SetBinLabel(2,"Rejected due to mismatch in trees");
118  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with good AOD");
119  fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
120  fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to vertex reco");
121  fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to pileup");
122  fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to centrality");
123  fHistNEvents->GetXaxis()->SetBinLabel(8,"Rejected due to vtxz");
124  fHistNEvents->GetXaxis()->SetBinLabel(9,"Rejected due to Physics Sel");
125  fHistNEvents->GetXaxis()->SetBinLabel(10,"nEvents accepted");
126  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
127  fHistNEvents->SetMinimum(0);
128  fOutput->Add(fHistNEvents);
129 
130  fHistNGenD=new TH1F("fHistNGenD","number of D mesons",10,-0.5,9.5);
131  TString names[5]={"Dzero","Dplus","Dstar","Ds","Lc"};
132  TString type[2]={"Prompt","Feeddown"};
133  for(Int_t j=0; j<5; j++){
134  for(Int_t i=0; i<2; i++){
135  Int_t index=j*2+i;
136  fHistNGenD->GetXaxis()->SetBinLabel(index+1,Form("%s %s",type[i].Data(),names[j].Data()));
137  }
138  }
139  fHistNGenD->GetXaxis()->SetNdivisions(1,kFALSE);
140  fHistNGenD->SetMinimum(0);
141  fOutput->Add(fHistNGenD);
142 
143  fHistNCand=new TH2F("fHistNCand","number of D meson candidates",5,-0.5,4.5,16,0.,16.);
144  fHistNCand->GetXaxis()->SetBinLabel(1,"D0#rightarrowK#pi");
145  fHistNCand->GetXaxis()->SetBinLabel(2,"D+#rightarrowK#pi#pi");
146  fHistNCand->GetXaxis()->SetBinLabel(3,"D*+#rightarrowD0#pi");
147  fHistNCand->GetXaxis()->SetBinLabel(4,"D_s#rightarrowKK#pi");
148  fHistNCand->GetXaxis()->SetBinLabel(5,"Lc#rightarrowpK#pi");
149  fHistNCand->SetMinimum(0);
150  fOutput->Add(fHistNCand);
151 
152  for(Int_t j=0; j<kDecays; j++){
153  for(Int_t i=0; i<2; i++){
154  Int_t index=j*2+i;
155  fHistPtYMultGen[index]=new TH3F(Form("hPtYMult%sGen%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",16,0.,16.,24,-6.,6.,40,0.,100.);
156  fHistPtYMultGenDauInAcc[index]=new TH3F(Form("hPtYMult%sGenDauInAcc%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",16,0.,16.,24,-6.,6.,40,0.,100.);
157  fHistPtYMultReco[index]=new TH3F(Form("hPtYMult%sReco%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",16,0.,16.,24,-6.,6.,40,0.,100.);
158  fHistPtYMultRecoFilt[index]=new TH3F(Form("hPtYMult%sRecoFilt%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",16,0.,16.,24,-6.,6.,40,0.,100.);
159  fHistPtYMultRecoSel[index]=new TH3F(Form("hPtYMult%sRecoSel%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; y; multPercentile",16,0.,16.,24,-6.,6.,40,0.,100.);
160  fOutput->Add(fHistPtYMultGen[index]);
161  fOutput->Add(fHistPtYMultGenDauInAcc[index]);
162  fOutput->Add(fHistPtYMultReco[index]);
163  fOutput->Add(fHistPtYMultRecoFilt[index]);
164  fOutput->Add(fHistPtYMultRecoSel[index]);
165  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)",16,0.,16.,100,-500.,500.);
166  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)",16,0.,16.,100,-500.,500.);
167  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)",16,0.,16.,100,-500.,500.);
168  fHistInvMassVsPt[index]=new TH2F(Form("hInvMassVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; Inv. Mass (GeV/c^{2})",16,0.,16.,300,1.75,2.35);
169  fHistDecLenVsPt[index]=new TH2F(Form("hDecLenVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; dec. len. (#mum)",16,0.,16.,100,0.,5000.);
170  fHistNormDLxyVsPt[index]=new TH2F(Form("hNormDLxyVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; dec. len. (#mum)",16,0.,16.,100,0.,20.);
171  fHistCosPointVsPt[index]=new TH2F(Form("hCosPointVsPt%s%s",type[i].Data(),fPartName[j].Data())," ; p_{T} (GeV/c) ; cos(#vartheta_{point})",16,0.,16.,100,-1.,1.);
172  fOutput->Add(fHistXvtxResVsPt[index]);
173  fOutput->Add(fHistYvtxResVsPt[index]);
174  fOutput->Add(fHistZvtxResVsPt[index]);
175  fOutput->Add(fHistInvMassVsPt[index]);
176  fOutput->Add(fHistDecLenVsPt[index]);
177  fOutput->Add(fHistNormDLxyVsPt[index]);
178  fOutput->Add(fHistCosPointVsPt[index]);
179  }
180  }
181 
182 
183  PostData(1,fOutput);
184  return;
185 }
186 //________________________________________________________________________
188 {
191 
192  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
193 
194  fHistNEvents->Fill(0); // count event
195 
196  if(fAODProtection>=0){
197  // Protection against different number of events in the AOD and deltaAOD
198  // In case of discrepancy the event is rejected.
199  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
200  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
201  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
202  fHistNEvents->Fill(1);
203  return;
204  }
205  }
206 
207  // Load all the branches of the DeltaAOD - needed for SelectionBit counting
208  TClonesArray *arrayD0toKpi =0;
209  TClonesArray *array3Prong =0;
210  TClonesArray *arrayDstar =0;
211  TClonesArray *arrayCascades =0;
212  if(!aod && AODEvent() && IsStandardAOD()) {
213  // In case there is an AOD handler writing a standard AOD, use the AOD
214  // event in memory rather than the input (ESD) event.
215  aod = dynamic_cast<AliAODEvent*> (AODEvent());
216  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
217  // have to taken from the AOD event hold by the AliAODExtension
218  AliAODHandler* aodHandler = (AliAODHandler*)
219  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
220  if(aodHandler->GetExtensions()) {
221 
222  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
223  AliAODEvent *aodFromExt = ext->GetAOD();
224 
225  array3Prong =(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
226  arrayD0toKpi =(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
227  arrayDstar =(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
228  arrayCascades=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
229  }
230  } else if(aod) {
231  array3Prong =(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
232  arrayD0toKpi =(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
233  arrayDstar =(TClonesArray*)aod->GetList()->FindObject("Dstar");
234  arrayCascades=(TClonesArray*)aod->GetList()->FindObject("CascadesHF");
235  }
236 
237  if(!aod || !array3Prong || !arrayD0toKpi) {
238  printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
239  return;
240  }
241  // fix for temporary bug in ESDfilter
242  // the AODs with null vertex pointer didn't pass the PhysSel
243  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
244  fHistNEvents->Fill(2);
245 
246  //Event selection
247  Bool_t isEvSel=fRDHFCuts->IsEventSelected(aod);
248  if(fRDHFCuts->GetWhyRejection()==5) fHistNEvents->Fill(3);
249  if(!isEvSel && fRDHFCuts->GetWhyRejection()==0) fHistNEvents->Fill(4);
250  if(fRDHFCuts->GetWhyRejection()==1) fHistNEvents->Fill(5);
251  if(fRDHFCuts->GetWhyRejection()==2) fHistNEvents->Fill(6);
252  if(fRDHFCuts->GetWhyRejection()==6)fHistNEvents->Fill(7);
253  if(fRDHFCuts->GetWhyRejection()==7)fHistNEvents->Fill(8);
254 
255  // load MC header
256  AliAODMCHeader *mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
257  if(!mcHeader) {
258  printf("AliAnalysisTaskDmesonMCPerform:UserExec: MC header branch not found!\n");
259  return;
260  }
261 
262  // load MC particles
263  TClonesArray *arrayMC= (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
264  if(!arrayMC) {
265  printf("AliAnalysisTaskDmesonMCPerform::UserExec: MC particles branch not found!\n");
266  return;
267  }
268 
269  Int_t runNumber=aod->GetRunNumber();
270 
271  if(aod->GetTriggerMask()==0 &&
272  (runNumber>=195344 && runNumber<=195677)){
273  // protection for events with empty trigger mask in p-Pb
274  return;
275  }
277 
278  PostData(1,fOutput);
279 
280  for(Int_t i=0; i<kMaxLabel; i++) fMapTrLabel[i]=-999;
281  MapTrackLabels(aod);
282  FillGenLevelHistos(aod,arrayMC,mcHeader);
283  if(!isEvSel)return;
284  fHistNEvents->Fill(9);
285 
286  FillCandLevelHistos(0,aod,arrayD0toKpi,arrayMC);
287  FillCandLevelHistos(1,aod,array3Prong,arrayMC);
288  FillCandLevelHistos(2,aod,arrayDstar,arrayMC);
289 
290  return;
291 }
292 //________________________________________________________________________
293 void AliAnalysisTaskDmesonMCPerform::FillGenLevelHistos(AliAODEvent* aod, TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
294  //
296 
297  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
298  if(TMath::Abs(zMCVertex)>fRDHFCuts->GetMaxVtxZ()) return;
299  Double_t evCentr=1.;
300  if(fRDHFCuts->GetUseCentrality()>0) evCentr=fRDHFCuts->GetCentrality(aod);
301 
302  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
303  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
304  Int_t pdgCode=TMath::Abs(mcPart->GetPdgCode());
305  Int_t iSpec=-1;
306  if(pdgCode == 411) iSpec=1;
307  else if(pdgCode == 413) iSpec=2;
308  else if(pdgCode == 421) iSpec=0;
309  else if(pdgCode == 431) iSpec=3;
310  else if(pdgCode == 4122) iSpec=4;
311  if (iSpec>=0){
312  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,kTRUE);//Prompt = 4, FeedDown = 5
313  if(orig<4 || orig>5) continue;
314  Int_t indexh=iSpec*2+(orig-4);
315  fHistNGenD->Fill(indexh);
316  Int_t deca=0;
317  Bool_t isGoodDecay=kFALSE;
318  Int_t labDau[4]={-1,-1,-1,-1};
319  Int_t nProng=3;
320  if (pdgCode == 411){
321  deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,mcPart,labDau);
322  if(deca>0) isGoodDecay=kTRUE;
323  }else if(pdgCode == 421){
324  deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,mcPart,labDau);
325  if(mcPart->GetNDaughters()!=2) continue;
326  if(deca==1) isGoodDecay=kTRUE;
327  nProng=2;
328  }else if(pdgCode == 431){
329  deca=AliVertexingHFUtils::CheckDsDecay(arrayMC,mcPart,labDau);
330  if(deca==1) isGoodDecay=kTRUE;
331  }else if(pdgCode==413){
332  deca=AliVertexingHFUtils::CheckDstarDecay(arrayMC,mcPart,labDau);
333  if(deca==1) isGoodDecay=kTRUE;
334  }else if(pdgCode==4122){
335  deca=AliVertexingHFUtils::CheckLcpKpiDecay(arrayMC,mcPart,labDau);
336  if(deca>=1) isGoodDecay=kTRUE;
337  }
338  if(labDau[0]==-1){
339  continue; //protection against unfilled array of labels
340  }
341  if(isGoodDecay && indexh>=0){
342  Double_t ptgen=mcPart->Pt();
343  Double_t ygen=mcPart->Y();
344  fHistPtYMultGen[indexh]->Fill(ptgen,ygen,evCentr);
345  if(fRDHFCuts->IsInFiducialAcceptance(ptgen,ygen)){
346  Bool_t dauInAcc=CheckAcceptance(arrayMC,nProng,labDau);
347  if(dauInAcc){
348  fHistPtYMultGenDauInAcc[indexh]->Fill(ptgen,ygen,evCentr);
349  AliAODRecoDecayHF* dmes=GetRecoDecay(aod,nProng,labDau);
350  if(dmes){
351  fHistPtYMultReco[indexh]->Fill(ptgen,ygen,evCentr);
352  }
353  delete dmes;
354  }
355  }
356  }
357  }
358  }
359 }
360 //________________________________________________________________________
361 void AliAnalysisTaskDmesonMCPerform::FillCandLevelHistos(Int_t idCase, AliAODEvent* aod, TClonesArray *arrayDcand, TClonesArray *arrayMC){
362  //
364 
365  Double_t evCentr=1.;
366  if(fRDHFCuts->GetUseCentrality()>0) evCentr=fRDHFCuts->GetCentrality(aod);
367 
368  Int_t pdg0[2]={321,211};
369  Int_t pdgp[3]={321,211,211};
370  Int_t pdgs[3]={321,211,321};
371  Int_t pdgst[2]={421,211};
372  Int_t pdgl[3]={2212,321,211};
373  // vHF object is needed to call the method that refills the missing info of the candidates
374  // if they have been deleted in dAOD reconstruction phase
375  // in order to reduce the size of the file
377 
378  Int_t nCand=arrayDcand->GetEntriesFast();
379  for (Int_t iCand = 0; iCand < nCand; iCand++) {
380  AliAODRecoDecayHF *d=(AliAODRecoDecayHF*)arrayDcand->UncheckedAt(iCand);
381  Double_t ptCand=-999.;
382  Int_t iSpec=-1;
383  Int_t labD=-1;
384  Double_t invMass=0.;
385  if(idCase==0){
386  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF2Prong*)d))continue;
387  ptCand=d->Pt();
388  fHistNCand->Fill(0.,ptCand);
389  labD = d->MatchToMC(421,arrayMC,2,pdg0);
390  if(labD>=0){
391  iSpec=0;
392  AliAODMCParticle *pd0 = (AliAODMCParticle*)arrayMC->At(labD);
393  if(pd0->GetPdgCode()==421) invMass=((AliAODRecoDecayHF2Prong*)d)->InvMassD0();
394  else invMass=((AliAODRecoDecayHF2Prong*)d)->InvMassD0bar();
395  }
396  }else if(idCase==1){
397  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF3Prong*)d))continue;
398  ptCand=d->Pt();
399  if(d->HasSelectionBit(AliRDHFCuts::kDplusCuts)) fHistNCand->Fill(1.,ptCand);
400  if(d->HasSelectionBit(AliRDHFCuts::kDsCuts)) fHistNCand->Fill(3.,ptCand);
401  if(d->HasSelectionBit(AliRDHFCuts::kLcCuts)) fHistNCand->Fill(4.,ptCand);
402  labD = d->MatchToMC(411,arrayMC,3,pdgp);
403  if(labD>=0){
405  invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDplus();
406  }else{
407  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
408  AliAODMCParticle* pdau0=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
409  Int_t pdgCode0=TMath::Abs(pdau0->GetPdgCode());
410  labD = d->MatchToMC(431,arrayMC,3,pdgs);
411  if(labD>=0){
412  if(d->HasSelectionBit(AliRDHFCuts::kDsCuts)) iSpec=3;
413  if(pdgCode0==321) invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDsKKpi();
414  else invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassDspiKK();
415  }else{
416  labD = d->MatchToMC(4122,arrayMC,3,pdgl);
417  if(labD>=0){
418  if(d->HasSelectionBit(AliRDHFCuts::kLcCuts)) iSpec=4;
419  if(pdgCode0==2212) invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassLcpKpi();
420  else invMass=((AliAODRecoDecayHF3Prong*)d)->InvMassLcpiKp();
421  }
422  }
423  }
424  }else if(idCase==2){
425  if(!vHF->FillRecoCasc(aod,((AliAODRecoCascadeHF*)d),kTRUE))continue;
426  ptCand=d->Pt();
427  fHistNCand->Fill(2.,ptCand);
428  labD = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgst,pdg0,arrayMC);
429  if(labD>=0){
430  iSpec=2;
431  invMass=((AliAODRecoCascadeHF*)d)->InvMassDstarKpipi();
432  }
433  }
434 
435  if(labD>=0 && iSpec>=0){
436  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
437  Double_t ptgen=partD->Pt();
438  Double_t ygen=partD->Y();
439  Double_t ptreco=d->Pt();
440  Double_t dlen=d->DecayLength()*10000.; //um
441  Double_t ndlenxy=d->NormalizedDecayLengthXY();
442  Double_t cosp=d->CosPointingAngle();
443  Double_t dx=(d->Xv()-partD->Xv())*10000.;
444  Double_t dy=(d->Yv()-partD->Yv())*10000.;
445  Double_t dz=(d->Zv()-partD->Zv())*10000.;
446  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,partD,kTRUE);//Prompt = 4, FeedDown = 5
447  if(orig<4 || orig>5) continue;
448  Int_t indexh=iSpec*2+(orig-4);
449  fHistPtYMultRecoFilt[indexh]->Fill(ptgen,ygen,evCentr);
450  if(iSpec==1){
451  if(fRDHFCutsDplus){
452  if(fRDHFCutsDplus->IsEventSelected(aod) &&
454  fHistPtYMultRecoSel[indexh]->Fill(ptgen,ygen,evCentr);
455  }
456  }
457  }
458  fHistXvtxResVsPt[indexh]->Fill(ptreco,dx);
459  fHistYvtxResVsPt[indexh]->Fill(ptreco,dy);
460  fHistZvtxResVsPt[indexh]->Fill(ptreco,dz);
461  fHistInvMassVsPt[indexh]->Fill(ptreco,invMass);
462  fHistDecLenVsPt[indexh]->Fill(ptreco,dlen);
463  fHistNormDLxyVsPt[indexh]->Fill(ptreco,ndlenxy);
464  fHistCosPointVsPt[indexh]->Fill(ptreco,cosp);
465  }
466  }
467  delete vHF;
468  return;
469 }
470 //_________________________________________________________________
471 Bool_t AliAnalysisTaskDmesonMCPerform::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
473  for (Int_t iProng = 0; iProng<nProng; iProng++){
474  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
475  if(!mcPartDaughter) return kFALSE;
476  Double_t eta = mcPartDaughter->Eta();
477  Double_t pt = mcPartDaughter->Pt();
478  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
479  }
480  return kTRUE;
481 }
482 //________________________________________________________________________
485  //
486  Int_t nTracks=aod->GetNumberOfTracks();
487 
488  for(Int_t it=0; it<nTracks; it++) {
489  AliAODTrack *tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
490  if(!tr) continue;
491  if(tr->GetID()<0) continue;
492  Int_t lab=TMath::Abs(tr->GetLabel());
493  if(lab<kMaxLabel) fMapTrLabel[lab]=it;
494  else printf("Label %d exceeds upper limit\n",lab);
495  }
496  return;
497 }
498 //________________________________________________________________________
500  //
502  //
503  Double_t px[3],py[3],pz[3],d0[3],d0err[3];
504  Int_t nFound=0;
505  for(Int_t jd=0; jd<nProng; jd++){
506  Int_t it=fMapTrLabel[labDau[jd]];
507  if(it>0){
508  AliAODTrack* tr=dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
509  if(tr){
510  if(TMath::Abs(tr->GetLabel())!=labDau[jd]){
511  printf("Mismatched track label\n");
512  continue;
513  }
514  px[jd]=tr->Px();
515  py[jd]=tr->Py();
516  px[jd]=tr->Pz();
517  d0[jd]=0.; // temporary
518  d0err[jd]=0.; // temporary
519  nFound++;
520  }
521  }
522  }
523  if(nFound!=nProng) return 0x0;
524  AliAODRecoDecayHF* dmes=new AliAODRecoDecayHF(0x0,nProng,nProng-2,px,py,pz,d0,d0err);
525  return dmes;
526 }
527 //________________________________________________________________________
529 {
531  //
532  if(fDebug > 1) printf("AnalysisTaskDmesonMCPerform: Terminate() \n");
533 
534  fOutput = dynamic_cast<TList*> (GetOutputData(1));
535  if (!fOutput) {
536  printf("ERROR: fOutput not available\n");
537  return;
538  }
539 
540  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
541  if(fHistNEvents){
542  printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(10));
543  }else{
544  printf("ERROR: fHistNEvents not available\n");
545  return;
546  }
547 
548  return;
549 }
Double_t NormalizedDecayLengthXY() const
static Int_t CheckD0Decay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
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)
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)
static Int_t CheckLcpKpiDecay(AliMCEvent *mcEvent, Int_t label, Int_t *arrayDauLab)
void SetUseCentrality(Int_t flag=1)
static Int_t CheckMatchingAODdeltaAODevents()
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:304
void SetUseAnyTrigger()
Definition: AliRDHFCuts.h:69
TH2F * fHistDecLenVsPt[2 *kDecays]
! hist. of decay length (meas)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
TH3F * fHistPtYMultGen[2 *kDecays]
! hist. for generated D
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:247
void SetUsePhysicsSelection(Bool_t use=kTRUE)
Definition: AliRDHFCuts.h:342
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
Float_t GetCentrality(AliAODEvent *aodEvent)
Definition: AliRDHFCuts.h:248
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)
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
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
void SetTriggerClass(TString trclass0, TString trclass1="")
Definition: AliRDHFCuts.h:194
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:270
TH2F * fHistZvtxResVsPt[2 *kDecays]
! hist. for sec vert x resol
void SetOptPileup(Int_t opt=0)
Definition: AliRDHFCuts.h:220
Double_t DecayLength() const
Bool_t CheckAcceptance(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
TList * fOutput
! list send on output slot 0
Int_t fMapTrLabel[kMaxLabel]
Names of hadron species.