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