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