AliPhysics  251aa1e (251aa1e)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSECompareHF.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2008, 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 /* $Id$ */
17 
19 //
20 // AliAnalysisTaskSE for the comparison of heavy flavor
21 // decay candidates with the MC truth.
22 //
23 // Author: A.Dainese, andrea.dainese@lnl.infn.it
25 
26 #include <TClonesArray.h>
27 #include <TNtuple.h>
28 #include <TList.h>
29 #include <TH1F.h>
30 
31 #include "AliAnalysisManager.h"
32 #include "AliAODHandler.h"
33 #include "AliAODEvent.h"
34 #include "AliAODVertex.h"
35 #include "AliESDtrack.h"
36 #include "AliAODTrack.h"
37 #include "AliAODMCHeader.h"
38 #include "AliAODMCParticle.h"
42 #include "AliAODRecoCascadeHF.h"
43 #include "AliAnalysisVertexingHF.h"
44 #include "AliAnalysisTaskSE.h"
46 
50 
51 //________________________________________________________________________
54 fOutput(0),
55 fNtupleCmp(0),
56 fHistMass(0),
57 fHistNEvents(0)
58 {
60 
61  // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
62 }
63 
64 //________________________________________________________________________
66 AliAnalysisTaskSE(name),
67 fOutput(0),
68 fNtupleCmp(0),
69 fHistMass(0),
70 fHistNEvents(0)
71 {
73 
74  // Output slot #1 writes into a TList container
75  DefineOutput(1,TList::Class()); //My private output
76  // Output slot #2 writes into a TNtuple container
77  DefineOutput(2,TNtuple::Class()); //My private output
78 }
79 
80 //________________________________________________________________________
82 {
84  if (fOutput) {
85  delete fOutput;
86  fOutput = 0;
87  }
88 }
89 
90 //________________________________________________________________________
92 {
94 
95  if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
96 
97  return;
98 }
99 
100 //________________________________________________________________________
102 {
104  //
105  if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
106 
107  // Several histograms are more conveniently managed in a TList
108  fOutput = new TList();
109  fOutput->SetOwner();
110 
111  fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
112  fHistMass->Sumw2();
113  fHistMass->SetMinimum(0);
114  fOutput->Add(fHistMass);
115 
116  fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
117  fHistNEvents->Sumw2();
118  fHistNEvents->SetMinimum(0);
119  fOutput->Add(fHistNEvents);
120 
121  OpenFile(2); // 2 is the slot number of the ntuple
122  fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
123 
124  return;
125 }
126 
127 //________________________________________________________________________
129 {
132 
133 
134  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
135 
136  TClonesArray *inputArrayVertices = 0;
137  TClonesArray *inputArrayD0toKpi = 0;
138  TClonesArray *inputArrayDstar = 0;
139 
140  if(!aod && AODEvent() && IsStandardAOD()) {
141  // In case there is an AOD handler writing a standard AOD, use the AOD
142  // event in memory rather than the input (ESD) event.
143  aod = dynamic_cast<AliAODEvent*> (AODEvent());
144  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
145  // have to taken from the AOD event hold by the AliAODExtension
146  AliAODHandler* aodHandler = (AliAODHandler*)
147  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
148  if(aodHandler->GetExtensions()) {
149  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
150  AliAODEvent *aodFromExt = ext->GetAOD();
151  // load HF vertices
152  inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
153  // load D0->Kpi candidates
154  inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
155  // load D*+ candidates
156  inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
157  }
158  } else if(aod) {
159  // load HF vertices
160  inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
161  // load D0->Kpi candidates
162  inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
163  // load D*+ candidates
164  inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
165  }
166 
167 
168  if(!inputArrayVertices || !aod) {
169  printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
170  return;
171  }
172  if(!inputArrayD0toKpi) {
173  printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
174  return;
175  }
176  if(!inputArrayDstar) {
177  printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
178  return;
179  }
180 
181 
182  fHistNEvents->Fill(0); // count event
183  // Post the data already here
184  PostData(1,fOutput);
185 
186  // AOD primary vertex
187  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
188  //vtx1->Print();
189 
190  // load MC particles
191  TClonesArray *mcArray =
192  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
193  if(!mcArray) {
194  printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
195  return;
196  }
197 
198  // load MC header
199  AliAODMCHeader *mcHeader =
200  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
201  if(!mcHeader) {
202  printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
203  return;
204  }
205 
206  Int_t nprongs,lab,okD0,okD0bar,pdg;
207  Bool_t unsetvtx;
208  Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
212 
213  Int_t pdgDgD0toKpi[2]={321,211};
214  Int_t pdgDgDplustoKpipi[3]={321,211,211};
215  Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
216 
217  // loop over vertices
218  Int_t nVertices = inputArrayVertices->GetEntriesFast();
219  if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
220 
221  for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
222  AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
223 
224  vtx->GetXYZ(posRec);
225  vtx->GetCovarianceMatrix(covRec);
226  errx=1.; erry=1.; errz=1.;
227  if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
228  if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
229  if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
230 
231 
232  // get parent AliAODRecoDecayHF
233  nprongs= vtx->GetNDaughters();
234  // check that the daughters have kITSrefit and kTPCrefit
235  Bool_t allDgOK=kTRUE;
236  for(Int_t idg=0; idg<nprongs; idg++) {
237  AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
238  if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
239  if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
240  }
241  if(!allDgOK) continue;
242 
243 
244  switch(nprongs) {
245  case 2: // look for D0->Kpi
246  d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
247  if(d2->IsLikeSign()) continue;
248  if(d2->Charge() != 0) continue; // these are D*
249  lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
250  if(lab>=0) {
251  unsetvtx=kFALSE;
252  if(!d2->GetOwnPrimaryVtx()) {
253  d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
254  unsetvtx=kTRUE;
255  }
256  okD0=1; okD0bar=1;
257  AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
258  pdg = dMC->GetPdgCode();
259  invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
260  // get a daughter for true pos of decay vertex
261  AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
262  dg0MC->XvYvZv(posTrue);
263  fHistMass->Fill(invmass);
264  // Post the data already here
265  PostData(1,fOutput);
266  Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
267  (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
268  (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
269  (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
270  (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
271  (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
272  fNtupleCmp->Fill(tmp);
273  PostData(2,fNtupleCmp);
274 
275  if(unsetvtx) d2->UnsetOwnPrimaryVtx();
276  }
277  break;
278  case 3: // look for D+
279  d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
280  if(d3->IsLikeSign()) continue;
281  lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
282  if(lab>=0) {
283  unsetvtx=kFALSE;
284  if(!d3->GetOwnPrimaryVtx()) {
285  d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
286  unsetvtx=kTRUE;
287  }
288  AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
289  pdg = dMC->GetPdgCode();
290  invmass = d3->InvMassDplus();
291  // get a daughter for true pos of decay vertex
292  AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
293  dg0MC->XvYvZv(posTrue);
294  Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
295  (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
296  (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
297  (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
298  (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
299  (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
300  fNtupleCmp->Fill(tmp);
301  PostData(2,fNtupleCmp);
302  if(unsetvtx) d3->UnsetOwnPrimaryVtx();
303  }
304  break;
305  case 4: // look for D0->Kpipipi
306  d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
307  if(d4->IsLikeSign()) continue;
308  lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
309  if(lab>=0) {
310  unsetvtx=kFALSE;
311  if(!d4->GetOwnPrimaryVtx()) {
312  d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
313  unsetvtx=kTRUE;
314  }
315  okD0=0; okD0bar=0;
316  AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
317  pdg = dMC->GetPdgCode();
318  //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
319  invmass = 10.; // dummy
320  // get a daughter for true pos of decay vertex
321  AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
322  dg0MC->XvYvZv(posTrue);
323  Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
324  (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
325  (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
326  (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
327  (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
328  (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
329  fNtupleCmp->Fill(tmp);
330  PostData(2,fNtupleCmp);
331  if(unsetvtx) d4->UnsetOwnPrimaryVtx();
332  }
333  break;
334  }
335 
336  } // end loop on vertices
337 
338 
339  // loop over D*+ candidates
340  /*
341  for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
342  AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
343  Int_t labDstar = c->MatchToMC(413,421,mcArray);
344  if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
345  }
346  */
347 
348  // Post the data already here
349  PostData(1,fOutput);
350  PostData(2,fNtupleCmp);
351 
352  return;
353 }
354 //________________________________________________________________________
356 {
358  //
359  if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
360 
361  fOutput = dynamic_cast<TList*> (GetOutputData(1));
362  if (!fOutput) {
363  printf("ERROR: fOutput not available\n");
364  return;
365  }
366 
367  fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
368  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
369 
370  //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
371 
372  return;
373 }
374 
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Int_t pdg
double Double_t
Definition: External.C:58
virtual void UserExec(Option_t *option)
virtual void Terminate(Option_t *option)
TH1F * fHistNEvents
! output histogram
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
AliAODVertex * GetOwnPrimaryVtx() const
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
TH1F * fHistMass
! output histogram
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TList * fOutput
! list send on output slot 0
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Bool_t IsLikeSign() const
check if it is like-sign
TNtuple * fNtupleCmp
! output ntuple
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65