Included also 3 and 4 prong decays; added variables to ntuple
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSECompareHF.cxx
CommitLineData
3a219f60 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/////////////////////////////////////////////////////////////
17//
18// AliAnalysisTaskSE for the comparison of heavy flavor
19// decay candidates with the MC truth.
20//
21// Author: A.Dainese, andrea.dainese@lnl.infn.it
22/////////////////////////////////////////////////////////////
23
24#include <TClonesArray.h>
25#include <TNtuple.h>
26#include <TList.h>
27#include <TH1F.h>
28
29#include "AliAODEvent.h"
30#include "AliAODVertex.h"
31#include "AliAODTrack.h"
32#include "AliAODMCHeader.h"
33#include "AliAODMCParticle.h"
34#include "AliAODRecoDecayHF2Prong.h"
b3999999 35#include "AliAODRecoCascadeHF.h"
2ae9c7cb 36#include "AliAnalysisVertexingHF.h"
3a219f60 37#include "AliAnalysisTaskSE.h"
38#include "AliAnalysisTaskSECompareHF.h"
39
40ClassImp(AliAnalysisTaskSECompareHF)
41
42
43//________________________________________________________________________
44AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
45AliAnalysisTaskSE(),
46fOutput(0),
2a041947 47fNtupleCmp(0),
12bfc069 48fHistMass(0),
49fVHF(0)
3a219f60 50{
51 // Default constructor
3a219f60 52
9d6d35b0 53 // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
3a219f60 54}
55
56//________________________________________________________________________
57AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
58AliAnalysisTaskSE(name),
59fOutput(0),
2a041947 60fNtupleCmp(0),
12bfc069 61fHistMass(0),
62fVHF(0)
3a219f60 63{
9d6d35b0 64 // Standard constructor
3a219f60 65
66 // Output slot #1 writes into a TList container
67 DefineOutput(1,TList::Class()); //My private output
10bdd1ae 68 // Output slot #2 writes into a TNtuple container
69 DefineOutput(2,TNtuple::Class()); //My private output
3a219f60 70}
71
72//________________________________________________________________________
73AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
74{
75 // Destructor
76 if (fOutput) {
77 delete fOutput;
78 fOutput = 0;
79 }
12bfc069 80 if (fVHF) {
81 delete fVHF;
82 fVHF = 0;
83 }
3a219f60 84}
85
86//________________________________________________________________________
87void AliAnalysisTaskSECompareHF::Init()
88{
89 // Initialization
90
91 if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
9d6d35b0 92
12bfc069 93 gROOT->LoadMacro("ConfigVertexingHF.C");
94
95 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
96 fVHF->PrintStatus();
9d6d35b0 97
3a219f60 98 return;
99}
100
101//________________________________________________________________________
102void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
103{
104 // Create the output container
105 //
106 if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
107
108 // Several histograms are more conveniently managed in a TList
109 fOutput = new TList();
110 fOutput->SetOwner();
111
112 fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
113 fHistMass->Sumw2();
114 fHistMass->SetMinimum(0);
115 fOutput->Add(fHistMass);
116
2a041947 117 fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec");
3a219f60 118
119 return;
120}
121
122//________________________________________________________________________
123void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
124{
125 // Execute analysis for current event:
126 // heavy flavor candidates association to MC truth
9d6d35b0 127
3a219f60 128
129 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
9d6d35b0 130
2a041947 131
132 // load HF vertices
133 TClonesArray *inputArrayVertices =
134 (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
135 if(!inputArrayVertices) {
136 printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
137 return;
138 }
139
3a219f60 140 // load D0->Kpi candidates
141 TClonesArray *inputArrayD0toKpi =
142 (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
143 if(!inputArrayD0toKpi) {
144 printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
145 return;
146 }
147
9d6d35b0 148
b3999999 149 // load D*+ candidates
150 TClonesArray *inputArrayDstar =
151 (TClonesArray*)aod->GetList()->FindObject("Dstar");
152 if(!inputArrayDstar) {
153 printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
154 return;
155 }
9d6d35b0 156
b3999999 157
3a219f60 158 // AOD primary vertex
159 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
160 //vtx1->Print();
161
162 // load MC particles
163 TClonesArray *mcArray =
164 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
165 if(!mcArray) {
166 printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
167 return;
168 }
169
170 // load MC header
171 AliAODMCHeader *mcHeader =
172 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
173 if(!mcHeader) {
174 printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
175 return;
176 }
177
2a041947 178 Int_t nprongs,lab,okD0,okD0bar,pdg;
179 Bool_t unsetvtx;
180 Double_t invmass,cov[6],errx,erry,errz;
181
182 // loop over vertices
183 Int_t nVertices = inputArrayVertices->GetEntriesFast();
184 printf("Number of vertices: %d\n",nVertices);
185
186 for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
187 AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
188
189 // get parent AliAODRecoDecayHF
190 nprongs= vtx->GetNDaughters();
191
192 switch(nprongs) {
193 case 2: // look for D0->Kpi
194 AliAODRecoDecayHF2Prong *d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
195 lab = d2->MatchToMC(421,mcArray);
196 if(lab>=0) {
197 unsetvtx=kFALSE;
198 if(!d2->GetOwnPrimaryVtx()) {
199 d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
200 unsetvtx=kTRUE;
201 }
202 okD0=0; okD0bar=0;
203 if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
204 AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
205 pdg = part->GetPdgCode();
206 invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
207 fHistMass->Fill(invmass);
208 // Post the data already here
209 PostData(1,fOutput);
210
211 AliAODVertex *secVtx = d2->GetSecondaryVtx();
212 secVtx->GetCovarianceMatrix(cov);
213 errx=1.; erry=1.; errz=1.;
214 if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
215 if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
216 if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
217
218 fNtupleCmp->Fill(pdg,nprongs,d2->Xv(),part->Xv(),errx,d2->Yv(),part->Yv(),erry,d2->Zv(),part->Zv(),errz,d2->GetReducedChi2(),d2->Pt(),invmass);
219 PostData(2,fNtupleCmp);
220 }
221 if(unsetvtx) d2->UnsetOwnPrimaryVtx();
3a219f60 222 }
2a041947 223 break;
224 case 3: // look for D+
225 AliAODRecoDecayHF3Prong *d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
226 lab = d3->MatchToMC(411,mcArray);
227 if(lab>=0) {
228 unsetvtx=kFALSE;
229 if(!d3->GetOwnPrimaryVtx()) {
230 d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
231 unsetvtx=kTRUE;
232 }
233 if(d3->SelectDplus(fVHF->GetDplusCuts())) {
234 AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
235 pdg = part->GetPdgCode();
236 invmass = d3->InvMassDplus();
237 AliAODVertex *secVtx = d3->GetSecondaryVtx();
238 secVtx->GetCovarianceMatrix(cov);
239 errx=1.; erry=1.; errz=1.;
240 if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
241 if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
242 if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
243
244 fNtupleCmp->Fill(pdg,nprongs,d3->Xv(),part->Xv(),errx,d3->Yv(),part->Yv(),erry,d3->Zv(),part->Zv(),errz,d3->GetReducedChi2(),d3->Pt(),invmass);
245 PostData(2,fNtupleCmp);
246 }
247 if(unsetvtx) d3->UnsetOwnPrimaryVtx();
248 }
249 break;
250 case 4: // look for D0->Kpipipi
251 AliAODRecoDecayHF4Prong *d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
252 lab = d4->MatchToMC(421,mcArray);
253 if(lab>=0) {
254 unsetvtx=kFALSE;
255 if(!d4->GetOwnPrimaryVtx()) {
256 d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
257 unsetvtx=kTRUE;
258 }
259 okD0=0; okD0bar=0;
260 if(d4->SelectD0(fVHF->GetD0to4ProngsCuts(),okD0,okD0bar)) {
261 AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(lab);
262 pdg = part->GetPdgCode();
263 //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
264 invmass = 10.;
265
266 AliAODVertex *secVtx = d4->GetSecondaryVtx();
267 secVtx->GetCovarianceMatrix(cov);
268 errx=1.; erry=1.; errz=1.;
269 if(cov[0]>0) errx = TMath::Sqrt(cov[0]);
270 if(cov[2]>0) erry = TMath::Sqrt(cov[2]);
271 if(cov[5]>0) errz = TMath::Sqrt(cov[5]);
272
273 fNtupleCmp->Fill(pdg,nprongs,d4->Xv(),part->Xv(),errx,d4->Yv(),part->Yv(),erry,d4->Zv(),part->Zv(),errz,d4->GetReducedChi2(),d4->Pt(),invmass);
274 PostData(2,fNtupleCmp);
275
276 }
277 if(unsetvtx) d4->UnsetOwnPrimaryVtx();
278 }
279 break;
3a219f60 280 }
3a219f60 281
2a041947 282 } // end loop on vertices
3a219f60 283
9d6d35b0 284
b3999999 285 // loop over D*+ candidates
4dd17b11 286 /*
b3999999 287 for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
288 AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
289 Int_t labDstar = c->MatchToMC(413,421,mcArray);
290 if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
291 }
4dd17b11 292 */
2a041947 293
4dd17b11 294 // Post the data already here
295 PostData(1,fOutput);
2a041947 296 PostData(2,fNtupleCmp);
b3999999 297
3a219f60 298 return;
299}
3a219f60 300//________________________________________________________________________
301void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
302{
303 // Terminate analysis
304 //
305 if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
306
307 fOutput = dynamic_cast<TList*> (GetOutputData(1));
308 if (!fOutput) {
309 printf("ERROR: fOutput not available\n");
310 return;
311 }
312
3a219f60 313 fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
314
2a041947 315 fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
10bdd1ae 316
3a219f60 317 return;
318}
319