1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // AliAnalysisTaskSE for the comparison of heavy flavor
19 // decay candidates with the MC truth.
21 // Author: A.Dainese, andrea.dainese@lnl.infn.it
22 /////////////////////////////////////////////////////////////
24 #include <TClonesArray.h>
29 #include "AliAnalysisManager.h"
30 #include "AliAODHandler.h"
31 #include "AliAODEvent.h"
32 #include "AliAODVertex.h"
33 #include "AliESDtrack.h"
34 #include "AliAODTrack.h"
35 #include "AliAODMCHeader.h"
36 #include "AliAODMCParticle.h"
37 #include "AliAODRecoDecayHF2Prong.h"
38 #include "AliAODRecoDecayHF3Prong.h"
39 #include "AliAODRecoDecayHF4Prong.h"
40 #include "AliAODRecoCascadeHF.h"
41 #include "AliAnalysisVertexingHF.h"
42 #include "AliAnalysisTaskSE.h"
43 #include "AliAnalysisTaskSECompareHF.h"
45 ClassImp(AliAnalysisTaskSECompareHF)
48 //________________________________________________________________________
49 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
57 // Default constructor
59 // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
62 //________________________________________________________________________
63 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
64 AliAnalysisTaskSE(name),
71 // Standard constructor
73 // Output slot #1 writes into a TList container
74 DefineOutput(1,TList::Class()); //My private output
75 // Output slot #2 writes into a TNtuple container
76 DefineOutput(2,TNtuple::Class()); //My private output
79 //________________________________________________________________________
80 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
93 //________________________________________________________________________
94 void AliAnalysisTaskSECompareHF::Init()
98 if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
100 gROOT->LoadMacro("ConfigVertexingHF.C");
102 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
108 //________________________________________________________________________
109 void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
111 // Create the output container
113 if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
115 // Several histograms are more conveniently managed in a TList
116 fOutput = new TList();
119 fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
121 fHistMass->SetMinimum(0);
122 fOutput->Add(fHistMass);
124 fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
125 fHistNEvents->Sumw2();
126 fHistNEvents->SetMinimum(0);
127 fOutput->Add(fHistNEvents);
129 OpenFile(2); // 2 is the slot number of the ntuple
130 fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
135 //________________________________________________________________________
136 void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
138 // Execute analysis for current event:
139 // heavy flavor candidates association to MC truth
142 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
144 TClonesArray *inputArrayVertices = 0;
145 TClonesArray *inputArrayD0toKpi = 0;
146 TClonesArray *inputArrayDstar = 0;
148 if(!aod && AODEvent() && IsStandardAOD()) {
149 // In case there is an AOD handler writing a standard AOD, use the AOD
150 // event in memory rather than the input (ESD) event.
151 aod = dynamic_cast<AliAODEvent*> (AODEvent());
152 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
153 // have to taken from the AOD event hold by the AliAODExtension
154 AliAODHandler* aodHandler = (AliAODHandler*)
155 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
156 if(aodHandler->GetExtensions()) {
157 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
158 AliAODEvent *aodFromExt = ext->GetAOD();
160 inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
161 // load D0->Kpi candidates
162 inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
163 // load D*+ candidates
164 inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
168 inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
169 // load D0->Kpi candidates
170 inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
171 // load D*+ candidates
172 inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
176 if(!inputArrayVertices || !aod) {
177 printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
180 if(!inputArrayD0toKpi) {
181 printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
184 if(!inputArrayDstar) {
185 printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
190 fHistNEvents->Fill(0); // count event
191 // Post the data already here
194 // AOD primary vertex
195 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
199 TClonesArray *mcArray =
200 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
202 printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
207 AliAODMCHeader *mcHeader =
208 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
210 printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
214 Int_t nprongs,lab,okD0,okD0bar,pdg;
216 Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
217 AliAODRecoDecayHF2Prong *d2=0;
218 AliAODRecoDecayHF3Prong *d3=0;
219 AliAODRecoDecayHF4Prong *d4=0;
221 Int_t pdgDgD0toKpi[2]={321,211};
222 Int_t pdgDgDplustoKpipi[3]={321,211,211};
223 Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
225 // loop over vertices
226 Int_t nVertices = inputArrayVertices->GetEntriesFast();
227 if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
229 for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
230 AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
233 vtx->GetCovarianceMatrix(covRec);
234 errx=1.; erry=1.; errz=1.;
235 if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
236 if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
237 if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
240 // get parent AliAODRecoDecayHF
241 nprongs= vtx->GetNDaughters();
242 // check that the daughters have kITSrefit and kTPCrefit
243 Bool_t allDgOK=kTRUE;
244 for(Int_t idg=0; idg<nprongs; idg++) {
245 AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
246 if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
247 if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
249 if(!allDgOK) continue;
253 case 2: // look for D0->Kpi
254 d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
255 if(d2->IsLikeSign()) continue;
256 if(d2->Charge() != 0) continue; // these are D*
257 lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
260 if(!d2->GetOwnPrimaryVtx()) {
261 d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
265 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
266 pdg = dMC->GetPdgCode();
267 invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
268 // get a daughter for true pos of decay vertex
269 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
270 dg0MC->XvYvZv(posTrue);
271 fHistMass->Fill(invmass);
272 // Post the data already here
274 Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
275 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
276 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
277 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
278 (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
279 (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
280 fNtupleCmp->Fill(tmp);
281 PostData(2,fNtupleCmp);
283 if(unsetvtx) d2->UnsetOwnPrimaryVtx();
286 case 3: // look for D+
287 d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
288 if(d3->IsLikeSign()) continue;
289 lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
292 if(!d3->GetOwnPrimaryVtx()) {
293 d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
296 //if(d3->SelectDplus(fVHF->GetDplusCuts())) {
297 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
298 pdg = dMC->GetPdgCode();
299 invmass = d3->InvMassDplus();
300 // get a daughter for true pos of decay vertex
301 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
302 dg0MC->XvYvZv(posTrue);
303 Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
304 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
305 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
306 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
307 (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
308 (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
309 fNtupleCmp->Fill(tmp);
310 PostData(2,fNtupleCmp);
312 if(unsetvtx) d3->UnsetOwnPrimaryVtx();
315 case 4: // look for D0->Kpipipi
316 d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
317 if(d4->IsLikeSign()) continue;
318 lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
321 if(!d4->GetOwnPrimaryVtx()) {
322 d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
326 //if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
327 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
328 pdg = dMC->GetPdgCode();
329 //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
330 invmass = 10.; // dummy
331 // get a daughter for true pos of decay vertex
332 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
333 dg0MC->XvYvZv(posTrue);
334 Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
335 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
336 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
337 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
338 (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
339 (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
340 fNtupleCmp->Fill(tmp);
341 PostData(2,fNtupleCmp);
343 if(unsetvtx) d4->UnsetOwnPrimaryVtx();
348 } // end loop on vertices
351 // loop over D*+ candidates
353 for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
354 AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
355 Int_t labDstar = c->MatchToMC(413,421,mcArray);
356 if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
360 // Post the data already here
362 PostData(2,fNtupleCmp);
366 //________________________________________________________________________
367 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
369 // Terminate analysis
371 if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
373 fOutput = dynamic_cast<TList*> (GetOutputData(1));
375 printf("ERROR: fOutput not available\n");
379 fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
380 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
382 //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));