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 "AliAODEvent.h"
30 #include "AliAODVertex.h"
31 #include "AliAODTrack.h"
32 #include "AliAODMCHeader.h"
33 #include "AliAODMCParticle.h"
34 #include "AliAODRecoDecayHF2Prong.h"
35 #include "AliAODRecoCascadeHF.h"
36 #include "AliAnalysisVertexingHF.h"
37 #include "AliAnalysisTaskSE.h"
38 #include "AliAnalysisTaskSECompareHF.h"
40 ClassImp(AliAnalysisTaskSECompareHF)
43 //________________________________________________________________________
44 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
52 // Default constructor
54 // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
57 //________________________________________________________________________
58 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
59 AliAnalysisTaskSE(name),
66 // Standard constructor
68 // Output slot #1 writes into a TList container
69 DefineOutput(1,TList::Class()); //My private output
70 // Output slot #2 writes into a TNtuple container
71 DefineOutput(2,TNtuple::Class()); //My private output
74 //________________________________________________________________________
75 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
88 //________________________________________________________________________
89 void AliAnalysisTaskSECompareHF::Init()
93 if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
95 gROOT->LoadMacro("ConfigVertexingHF.C");
97 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
103 //________________________________________________________________________
104 void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
106 // Create the output container
108 if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
110 // Several histograms are more conveniently managed in a TList
111 fOutput = new TList();
114 fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
116 fHistMass->SetMinimum(0);
117 fOutput->Add(fHistMass);
119 fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
120 fHistNEvents->Sumw2();
121 fHistNEvents->SetMinimum(0);
122 fOutput->Add(fHistNEvents);
124 fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec");
129 //________________________________________________________________________
130 void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
132 // Execute analysis for current event:
133 // heavy flavor candidates association to MC truth
136 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
138 fHistNEvents->Fill(0); // count event
139 // Post the data already here
143 TClonesArray *inputArrayVertices =
144 (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
145 if(!inputArrayVertices) {
146 printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
150 // load D0->Kpi candidates
151 TClonesArray *inputArrayD0toKpi =
152 (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
153 if(!inputArrayD0toKpi) {
154 printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
159 // load D*+ candidates
160 TClonesArray *inputArrayDstar =
161 (TClonesArray*)aod->GetList()->FindObject("Dstar");
162 if(!inputArrayDstar) {
163 printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
168 // AOD primary vertex
169 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
173 TClonesArray *mcArray =
174 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
176 printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
181 AliAODMCHeader *mcHeader =
182 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
184 printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
188 Int_t nprongs,lab,okD0,okD0bar,pdg;
190 Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
191 AliAODRecoDecayHF2Prong *d2=0;
192 AliAODRecoDecayHF3Prong *d3=0;
193 AliAODRecoDecayHF4Prong *d4=0;
195 // loop over vertices
196 Int_t nVertices = inputArrayVertices->GetEntriesFast();
197 printf("Number of vertices: %d\n",nVertices);
199 for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
200 AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
203 vtx->GetCovarianceMatrix(covRec);
204 errx=1.; erry=1.; errz=1.;
205 if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
206 if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
207 if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
210 // get parent AliAODRecoDecayHF
211 nprongs= vtx->GetNDaughters();
214 case 2: // look for D0->Kpi
215 d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
216 if(d2->Charge() != 0) continue; // these are D*
217 lab = d2->MatchToMC(421,mcArray);
220 if(!d2->GetOwnPrimaryVtx()) {
221 d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
225 //if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) { // no cuts, otherwise we bias the resolution
226 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
227 pdg = dMC->GetPdgCode();
228 invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
229 fHistMass->Fill(invmass);
230 // Post the data already here
232 // get a daughter for true pos of decay vertex
233 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
234 dg0MC->XvYvZv(posTrue);
235 fNtupleCmp->Fill(pdg,nprongs,
236 posRec[0],posTrue[0],errx,
237 posRec[1],posTrue[1],erry,
238 posRec[2],posTrue[2],errz,
239 d2->GetReducedChi2(),d2->Pt(),invmass);
240 PostData(2,fNtupleCmp);
242 if(unsetvtx) d2->UnsetOwnPrimaryVtx();
245 case 3: // look for D+
246 d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
247 lab = d3->MatchToMC(411,mcArray);
250 if(!d3->GetOwnPrimaryVtx()) {
251 d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
254 //if(d3->SelectDplus(fVHF->GetDplusCuts())) {
255 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
256 pdg = dMC->GetPdgCode();
257 invmass = d3->InvMassDplus();
258 // get a daughter for true pos of decay vertex
259 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
260 dg0MC->XvYvZv(posTrue);
261 fNtupleCmp->Fill(pdg,nprongs,
262 posRec[0],posTrue[0],errx,
263 posRec[1],posTrue[1],erry,
264 posRec[2],posTrue[2],errz,
265 d3->GetReducedChi2(),d3->Pt(),invmass);
266 PostData(2,fNtupleCmp);
268 if(unsetvtx) d3->UnsetOwnPrimaryVtx();
271 case 4: // look for D0->Kpipipi
272 d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
273 lab = d4->MatchToMC(421,mcArray);
276 if(!d4->GetOwnPrimaryVtx()) {
277 d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
281 //if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
282 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
283 pdg = dMC->GetPdgCode();
284 //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
285 invmass = 10.; // dummy
286 // get a daughter for true pos of decay vertex
287 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
288 dg0MC->XvYvZv(posTrue);
289 fNtupleCmp->Fill(pdg,nprongs,
290 posRec[0],posTrue[0],errx,
291 posRec[1],posTrue[1],erry,
292 posRec[2],posTrue[2],errz,
293 d4->GetReducedChi2(),d4->Pt(),invmass);
294 PostData(2,fNtupleCmp);
296 if(unsetvtx) d4->UnsetOwnPrimaryVtx();
301 } // end loop on vertices
304 // loop over D*+ candidates
306 for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
307 AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
308 Int_t labDstar = c->MatchToMC(413,421,mcArray);
309 if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
313 // Post the data already here
315 PostData(2,fNtupleCmp);
319 //________________________________________________________________________
320 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
322 // Terminate analysis
324 if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
326 fOutput = dynamic_cast<TList*> (GetOutputData(1));
328 printf("ERROR: fOutput not available\n");
332 fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
333 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
335 fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));