AliPhysics  1811c8f (1811c8f)
AliAnalysisTaskTrackingSysPropagation.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 
17 // AliAnalysisTask for Tracking Systematics (matching efficiency //
18 // + tracking efficiency) porpagation at the D-meson level //
19 // (including daughter's kinematics) //
21 
22 #include <TFile.h>
23 #include <TH1F.h>
24 #include <TF1.h>
25 #include <TH2F.h>
26 #include <TH3F.h>
27 #include <TChain.h>
28 #include <Rtypes.h>
29 #include "AliAnalysisTaskSE.h"
30 #include "AliAnalysisManager.h"
31 #include "AliAnalysisDataContainer.h"
32 #include "AliAODHandler.h"
33 #include "AliLog.h"
34 #include "AliMCEventHandler.h"
35 #include "AliMCEvent.h"
37 #include "AliMultiplicity.h"
38 #include "AliRDHFCuts.h"
40 #include "AliRDHFCutsDstoKKpi.h"
41 #include "AliRDHFCutsD0toKpi.h"
43 #include "AliAODMCParticle.h"
44 #include "AliAODMCHeader.h"
45 #include "AliAODRecoDecay.h"
48 #include "AliAODRecoCascadeHF.h"
49 #include "AliAnalysisVertexingHF.h"
50 
51 
53 /* $Id$ */
54 
55 //________________________________________________________________________
57 AliAnalysisTaskSE("taskTrackingSysProp"),
58  fPartName(""),
59  fOutput(0x0),
60  fAnalysisCuts(0),
61  fHistNEvents(0x0),
62  fHistMESyst(0x0),
63  fHistTrEffSyst(0x0),
64  fhPtDauVsD(0x0),
65  fhSystMatchEffD(0x0),
66  fDecayChannel(AliAnalysisTaskTrackingSysPropagation::kDplustoKpipi),
67  fPDGcode(411),
68  fAODProtection(1),
69  fMaxPt(60.)
70 {
71  DefineInput(0, TChain::Class());
72  DefineOutput(1, TList::Class());
73 }
74 
75 //________________________________________________________________________
77  AliAnalysisTaskSE("taskTrackingSysProp"),
78  fPartName(""),
79  fOutput(0x0),
80  fAnalysisCuts(cuts),
81  fHistNEvents(0x0),
82  fHistMESyst(0x0),
83  fHistTrEffSyst(0x0),
84  fhPtDauVsD(0x0),
85  fhSystMatchEffD(0x0),
86  fDecayChannel(ch),
87  fPDGcode(411),
88  fAODProtection(1),
89  fMaxPt(60.)
90 {
91  fHistMESyst = new TH1F(*(static_cast<TH1F*>(HistMESys)));
92  fHistTrEffSyst = new TH1F(*(static_cast<TH1F*>(HistTrEffSys)));
93 
94  DefineInput(0, TChain::Class());
95  DefineOutput(1, TList::Class());
96 }
97 
98 //___________________________________________________________________________
100 
101  if(fOutput && !fOutput->IsOwner()){
102  delete fHistNEvents;
103  delete fHistMESyst;
104  delete fHistTrEffSyst;
105  delete fhPtDauVsD;
106  delete fhSystMatchEffD;
107  }
108  delete fOutput;
109  delete fAnalysisCuts;
110  fOutput = 0;
111 }
112 
113 //________________________________________________________________________
115 
116  fOutput = new TList();
117  fOutput->SetOwner();
118  fOutput->SetName("OutputHistos");
119 
120  fHistNEvents = new TH1F("hNEvents", "number of events ",15,-0.5,14.5);
121  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsRead");
122  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents Matched dAOD");
123  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents Mismatched dAOD");
124  fHistNEvents->GetXaxis()->SetBinLabel(4,"nEventsAnal");
125  fHistNEvents->GetXaxis()->SetBinLabel(5,"n. passing IsEvSelected");
126  fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected due to trigger");
127  fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected due to not reco vertex");
128  fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for contr vertex");
129  fHistNEvents->GetXaxis()->SetBinLabel(9,"n. rejected for vertex out of accept");
130  fHistNEvents->GetXaxis()->SetBinLabel(10,"n. rejected for pileup events");
131  fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of out centrality events");
132  fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of 3 prong candidates");
133  fHistNEvents->GetXaxis()->SetBinLabel(13,"no. of D after filtering cuts");
134  fHistNEvents->GetXaxis()->SetBinLabel(14,"no. of D after selection cuts");
135  fHistNEvents->GetXaxis()->SetBinLabel(15,"no. of not on-the-fly rec D");
136 
137  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
138 
139  fHistNEvents->Sumw2();
140  fHistNEvents->SetMinimum(0);
141  fOutput->Add(fHistNEvents);
142 
143  int nbins = (int)(10*fMaxPt);
144  fhSystMatchEffD = new TH2F("fhSystMatchEffD","Matching Efficiency; p_{T} D; syst. (\%)",nbins,0.,fMaxPt,50,0.,25.);
145  fhPtDauVsD = new TH2F("fhPtDauVsD","Pt Dau vs D; p_{T} D; p_{T} daugh",nbins,0.,fMaxPt,nbins,0.,fMaxPt);
146 
147  fOutput->Add(fHistMESyst);
148  fOutput->Add(fHistTrEffSyst);
149  fOutput->Add(fhSystMatchEffD);
150  fOutput->Add(fhPtDauVsD);
151 
152  PostData(1,fOutput);
153 }
154 //________________________________________________________________________
156 
158  fAnalysisCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fAnalysisCuts)));
159  }
160  else if(fDecayChannel == kDstoKKpi) {
161  fAnalysisCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fAnalysisCuts)));
162  }
163  else if(fDecayChannel == kD0toKpi) {
164  fAnalysisCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fAnalysisCuts)));
165  }
166  else if(fDecayChannel == kDstartoKpipi) {
167  fAnalysisCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fAnalysisCuts)));
168  }
169  else {
170  AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
171  }
172 }
173 //________________________________________________________________________
175 
176  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fInputEvent);
177  fHistNEvents->Fill(0); // all events
178 
179  if(fAODProtection>=0){
180  // Protection against different number of events in the AOD and deltaAOD
181  // In case of discrepancy the event is rejected.
182  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
183  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
184  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
185  fHistNEvents->Fill(2);
186  PostData(1,fOutput);
187  return;
188  }
189  }
190  TClonesArray *arrayBranch=0;
191 
192  if(!aod && AODEvent() && IsStandardAOD()) {
193  // In case there is an AOD handler writing a standard AOD, use the AOD
194  // event in memory rather than the input (ESD) event.
195  aod = dynamic_cast<AliAODEvent*> (AODEvent());
196  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
197  // have to taken from the AOD event hold by the AliAODExtension
198  AliAODHandler* aodHandler = (AliAODHandler*)
199  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
200  if(aodHandler->GetExtensions()) {
201  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
202  AliAODEvent *aodFromExt = ext->GetAOD();
204  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
205  else if(fDecayChannel == kD0toKpi)
206  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
207  else if(fDecayChannel == kDstartoKpipi)
208  arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
209  }
210  }
211  else {
213  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
214  else if(fDecayChannel == kD0toKpi)
215  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
216  else if(fDecayChannel == kDstartoKpipi)
217  arrayBranch=(TClonesArray*)aod->GetList()->FindObject("Dstar");
218  }
219 
220  if (!arrayBranch) {
221  AliError("Could not find array of HF vertices");
222  PostData(1,fOutput);
223  return;
224  }
225 
226  AliAODVertex *aodVtx = (AliAODVertex*)aod->GetPrimaryVertex();
227  if (!aodVtx || TMath::Abs(aod->GetMagneticField())<0.001) {
228  AliDebug(3, "The event was skipped due to missing vertex or magnetic field issue");
229  PostData(1,fOutput);
230  return;
231  }
232  fHistNEvents->Fill(3); // count event
234  fPDGcode = 411;
235  }else if(fDecayChannel == kDstoKKpi){
236  fPDGcode = 431;
237  }else if(fDecayChannel == kD0toKpi){
238  fPDGcode = 421;
239  }else if(fDecayChannel == kDstartoKpipi){
240  fPDGcode = 413;
241  }
242  Bool_t isEvSel = fAnalysisCuts->IsEventSelected(aod);
243  Float_t ntracks = aod->GetNumberOfTracks();
244 
251 
252  Int_t runNumber = aod->GetRunNumber();
253 
254  TClonesArray *arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
255  AliAODMCHeader *mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
256 
257  if(!arrayMC) {
258  AliError("AliAnalysisTaskTrackingSysPropagation::UserExec: MC particles branch not found!\n");
259  PostData(1,fOutput);
260  return;
261  }
262  if(!mcHeader) {
263  AliError("AliAnalysisTaskTrackingSysPropagation::UserExec: MC header branch not found!\n");
264  PostData(1,fOutput);
265  return;
266  }
267 
268  // arrayMC->Print();
269  if(aod->GetTriggerMask()==0 && (runNumber>=195344 && runNumber<=195677)) {
270  // protection for events with empty trigger mask in p-Pb
271  PostData(1,fOutput);
272  return;
273  }
275  // events not passing the centrality selection can be removed immediately.
276  PostData(1,fOutput);
277  return;
278  }
279  Double_t zMCVertex = mcHeader->GetVtxZ();
280  if (TMath::Abs(zMCVertex) > fAnalysisCuts->GetMaxVtxZ()) {
281  PostData(1,fOutput);
282  return;
283  }
284  if(!isEvSel){
285  PostData(1,fOutput);
286  return;
287  }
288  fHistNEvents->Fill(4);
289 
290  Int_t nCand = arrayBranch->GetEntriesFast();
291 
292  Int_t nprongs = -1;
293  Int_t pdgDaughter[3];
294  Int_t pdg2Daughter[2];
296  fPDGcode = 411;
297  nprongs = 3;
298  fPartName="Dplus";
299  pdgDaughter[0]=321;
300  pdgDaughter[1]=211;
301  pdgDaughter[2]=211;
302  }else if(fDecayChannel == kDstoKKpi){
303  fPDGcode = 431;
304  nprongs = 3;
305  fPartName="Ds";
306  pdgDaughter[0]=321;
307  pdgDaughter[1]=321;
308  pdgDaughter[2]=211;
309  }else if(fDecayChannel == kD0toKpi){
310  fPDGcode = 421;
311  nprongs = 2;
312  fPartName="D0";
313  pdgDaughter[0]=321;
314  pdgDaughter[1]=211;
315  }else if(fDecayChannel == kDstartoKpipi){
316  fPDGcode = 413;
317  nprongs = 2;
318  fPartName="Dstar";
319  pdgDaughter[0]=421;
320  pdgDaughter[1]=211;
321  pdg2Daughter[0]=321;
322  pdg2Daughter[1]=211;
323  }else{
324  AliError("WRONG DECAY SETTING");
325  PostData(1,fOutput);
326  return;
327  }
328 
329  // vHF object is needed to call the method that refills the missing info of the candidates
330  // if they have been deleted in dAOD reconstruction phase
331  // in order to reduce the size of the file
333 
334  for (Int_t iCand = 0; iCand < nCand; iCand++) {
335 
336  AliAODRecoDecayHF* d = 0x0;
337 
339  d = (AliAODRecoDecayHF3Prong*)arrayBranch->UncheckedAt(iCand);
340  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF3Prong*)d)) {
341  fHistNEvents->Fill(14);
342  continue;
343  }
344  }
345  else if(fDecayChannel == kD0toKpi) {
346  d = (AliAODRecoDecayHF2Prong*)arrayBranch->UncheckedAt(iCand);
347  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF2Prong*)d)) {
348  fHistNEvents->Fill(14);
349  continue;
350  }
351  }
352  else if(fDecayChannel == kDstartoKpipi) {
353  d = (AliAODRecoCascadeHF*)arrayBranch->At(iCand);
354  if(!d) continue;
355  Bool_t isDStarCand =kTRUE;
356  if(!vHF->FillRecoCasc(aod,(AliAODRecoCascadeHF*)d,isDStarCand)) {
357  fHistNEvents->Fill(14);
358  continue;
359  }
360  if(!d->GetSecondaryVtx()) continue;
361  }
362 
363  fHistNEvents->Fill(11);
364  Bool_t unsetvtx=kFALSE;
365  if(!d->GetOwnPrimaryVtx()){
366  d->SetOwnPrimaryVtx(aodVtx);
367  unsetvtx=kTRUE;
368  }
369 
370  Bool_t recVtx=kFALSE;
371  AliAODVertex *origownvtx=0x0;
372 
373  Double_t ptD = d->Pt();
374  Double_t rapid = d->Y(fPDGcode);
375  Bool_t isFidAcc = fAnalysisCuts->IsInFiducialAcceptance(ptD,rapid);
376 
377  if(isFidAcc){
378  Int_t retCodeAnalysisCuts = fAnalysisCuts->IsSelected(d,AliRDHFCuts::kAll,aod);
379  if(retCodeAnalysisCuts > 0) {
380  fHistNEvents->Fill(13);
381  if(fDecayChannel == kDstoKKpi){
382  Int_t isPhiKKpi = retCodeAnalysisCuts&4; //for Ds
383  Int_t isPhipiKK = retCodeAnalysisCuts&8; //for Ds
384  if(!(isPhiKKpi || isPhipiKK)) continue;
385  }
386 
387  double syst = 0.;
388  Int_t nDau = d->GetNDaughters();
389  AliAODTrack *track;
390 
391  Int_t mcLabel=-1;
392  int nbinsME = fHistMESyst->GetNbinsX();
393  int bin = 0;
394  if(!(fDecayChannel == kDstartoKpipi)) {
395  mcLabel = d->MatchToMC(fPDGcode,arrayMC,nprongs,pdgDaughter);
396  if (mcLabel < 0) continue;
397  for(Int_t iDau=0; iDau < nDau; iDau++){
398  track = (AliAODTrack*)d->GetDaughter(iDau);
399  if(!track){
400  AliError("Daughter particle track not found");
401  PostData(1,fOutput);
402  return;
403  }
404  Int_t labDau = track->GetLabel();
405  AliAODMCParticle* p = (AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau));
406  double ptdau = track->Pt();
407 
408  bin = fHistMESyst->FindBin(ptdau);
409  if(bin > 0 && bin <= nbinsME) syst += fHistMESyst->GetBinContent(bin);
410  else if(bin > nbinsME) syst += fHistMESyst->GetBinContent(nbinsME);
411  else if(bin == 0) {
412  AliError("Check input histo at low pt!");
413  PostData(1,fOutput);
414  return;
415  }
416  fhPtDauVsD->Fill(ptD,ptdau);
417  }
418  }
419  else { //D*
420  mcLabel = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDaughter,pdg2Daughter,arrayMC,kFALSE);
421  if (mcLabel < 0) continue;
422 
423  //soft pion
424  track = (AliAODTrack*)(((AliAODRecoCascadeHF*)d)->GetBachelor());
425  if(!track) {
426  PostData(1,fOutput);
427  return;
428  }
429  Int_t labDau = track->GetLabel();
430  AliAODMCParticle* p = (AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau));
431  double ptdau = track->Pt();
432 
433  bin = fHistMESyst->FindBin(ptdau);
434  if(bin > 0 && bin <= nbinsME) syst += fHistMESyst->GetBinContent(bin);
435  else if(bin > nbinsME) syst += fHistMESyst->GetBinContent(nbinsME);
436  else if(bin == 0) {
437  AliError("Check input histo at low pt!");
438  PostData(1,fOutput);
439  return;
440  }
441  fhPtDauVsD->Fill(ptD,ptdau);
442  //D0
443  AliAODRecoDecayHF2Prong *D0 = ((AliAODRecoCascadeHF*)d)->Get2Prong();
444  for(Int_t iDau=0; iDau < 2; iDau++){
445  track = 0x0;
446  track = (AliAODTrack*)D0->GetDaughter(iDau);
447  if(!track) {
448  PostData(1,fOutput);
449  return;
450  }
451  labDau = track->GetLabel();
452  p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau));
453  ptdau = track->Pt();
454 
455  bin = fHistMESyst->FindBin(ptdau);
456  if(bin > 0 && bin <= nbinsME) syst += fHistMESyst->GetBinContent(bin);
457  else if(bin > nbinsME) syst += fHistMESyst->GetBinContent(nbinsME);
458  else if(bin == 0) {
459  AliError("Check input histo at low pt!");
460  PostData(1,fOutput);
461  return;
462  }
463  fhPtDauVsD->Fill(ptD,ptdau);
464  }
465  }
466  int nbinsTE = fHistTrEffSyst->GetNbinsX();
467  bin = fHistTrEffSyst->FindBin(ptD);
468  double trackEffSys = 0.;
469  if(bin > 0 && bin <= nbinsTE) trackEffSys = fHistTrEffSyst->GetBinContent(bin);
470  else if(bin > nbinsTE) trackEffSys = fHistTrEffSyst->GetBinContent(nbinsTE);
471  else if(bin == 0) {
472  AliError("Check input histo at low pt!");
473  PostData(1,fOutput);
474  return;
475  }
476  syst = TMath::Sqrt(syst*syst+trackEffSys*trackEffSys);
477  fhSystMatchEffD->Fill(ptD,syst);
478  }
479  }
480  if(unsetvtx) d->UnsetOwnPrimaryVtx();
481  }
482  PostData(1,fOutput);
483  return;
484 }
485 
486 //________________________________________________________________________
488 
489  fOutput = dynamic_cast<TList*>(GetOutputData(1));
490  if (!fOutput) {
491  AliError("ERROR: fOutput not available\n");
492  return;
493  }
494  if(fHistNEvents){
495  Printf("Number of Analyzed Events = %f",fHistNEvents->GetBinContent(1));
496  }else{
497  AliError("ERROR: fHistNEvents not available\n");
498  }
499 
500  Printf("end of Terminate");
501  return;
502 }
Bool_t IsEventRejectedDueToCentrality() const
Definition: AliRDHFCuts.h:330
Bool_t IsEventRejectedDueToZVertexOutsideFiducialRegion() const
Definition: AliRDHFCuts.h:324
Bool_t IsEventRejectedDueToNotRecoVertex() const
Definition: AliRDHFCuts.h:318
double Double_t
Definition: External.C:58
Definition: External.C:236
Int_t IsEventSelectedInCentrality(AliVEvent *event)
static Int_t CheckMatchingAODdeltaAODevents()
Bool_t IsEventRejectedDueToVertexContributors() const
Definition: AliRDHFCuts.h:321
TH2F * fhSystMatchEffD
histo with Pt daughters vs pt candidate
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
TH1F * fHistTrEffSyst
histo with match. eff. systematics vs pt (need to be passed as input)
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:256
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
DecChannel fDecayChannel
histo with systematic uncertainty on the candidate
Class for cuts on AOD reconstructed D+->Kpipi.
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Double_t fMaxPt
flag to activate protection against AOD-dAOD mismatch.
AliAODVertex * GetOwnPrimaryVtx() const
TH2F * fhPtDauVsD
histo with track. eff. systematics vs pt (need to be passed as input)
Bool_t IsEventRejectedDueToPileup() const
Definition: AliRDHFCuts.h:327
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Bool_t IsEventSelected(AliVEvent *event)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:290
const char Option_t
Definition: External.C:48
Bool_t IsEventRejectedDueToTrigger() const
Definition: AliRDHFCuts.h:315
const Int_t nbins
bool Bool_t
Definition: External.C:53
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:310
Int_t GetUseCentrality() const
Definition: AliRDHFCuts.h:279
AliRDHFCuts * fAnalysisCuts
max pt in the outputs histos