A. Dainese added a method to match the reconstructed particles to
[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"
2ae9c7cb 35#include "AliAnalysisVertexingHF.h"
3a219f60 36#include "AliAnalysisTaskSE.h"
37#include "AliAnalysisTaskSECompareHF.h"
38
39ClassImp(AliAnalysisTaskSECompareHF)
40
41
42//________________________________________________________________________
43AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
44AliAnalysisTaskSE(),
45fOutput(0),
46fNtupleD0Cmp(0),
12bfc069 47fHistMass(0),
48fVHF(0)
3a219f60 49{
50 // Default constructor
3a219f60 51
52 // Output slot #1 writes into a TList container
53 DefineOutput(1,TList::Class()); //My private output
54}
55
56//________________________________________________________________________
57AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
58AliAnalysisTaskSE(name),
59fOutput(0),
60fNtupleD0Cmp(0),
12bfc069 61fHistMass(0),
62fVHF(0)
3a219f60 63{
64 // Default constructor
3a219f60 65
66 // Output slot #1 writes into a TList container
67 DefineOutput(1,TList::Class()); //My private output
68}
69
70//________________________________________________________________________
71AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
72{
73 // Destructor
74 if (fOutput) {
75 delete fOutput;
76 fOutput = 0;
77 }
12bfc069 78 if (fVHF) {
79 delete fVHF;
80 fVHF = 0;
81 }
3a219f60 82}
83
84//________________________________________________________________________
85void AliAnalysisTaskSECompareHF::Init()
86{
87 // Initialization
88
89 if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
90
12bfc069 91 gROOT->LoadMacro("ConfigVertexingHF.C");
92
93 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
94 fVHF->PrintStatus();
95
3a219f60 96 return;
97}
98
99//________________________________________________________________________
100void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
101{
102 // Create the output container
103 //
104 if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
105
106 // Several histograms are more conveniently managed in a TList
107 fOutput = new TList();
108 fOutput->SetOwner();
109
110 fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
111 fHistMass->Sumw2();
112 fHistMass->SetMinimum(0);
113 fOutput->Add(fHistMass);
114
115 fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
116 fOutput->Add(fNtupleD0Cmp);
117
118 return;
119}
120
121//________________________________________________________________________
122void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
123{
124 // Execute analysis for current event:
125 // heavy flavor candidates association to MC truth
126
127 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
128
129 // load D0->Kpi candidates
130 TClonesArray *inputArrayD0toKpi =
131 (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
132 if(!inputArrayD0toKpi) {
133 printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
134 return;
135 }
136
137 // AOD primary vertex
138 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
139 //vtx1->Print();
140
141 // load MC particles
142 TClonesArray *mcArray =
143 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
144 if(!mcArray) {
145 printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
146 return;
147 }
148
149 // load MC header
150 AliAODMCHeader *mcHeader =
151 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
152 if(!mcHeader) {
153 printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
154 return;
155 }
156
157
158
159 // loop over D0->Kpi candidates
160 Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
161 printf("Number of D0->Kpi: %d\n",nInD0toKpi);
162
163 Int_t lab0,lab1,labMother,labD0daugh0,labD0daugh1,pdgMother,pdgD0;
164
165 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
166 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
167 labD0daugh0=-1;
168 labD0daugh1=-1;
169 Bool_t unsetvtx=kFALSE;
170 if(!d->GetOwnPrimaryVtx()) {
171 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
172 unsetvtx=kTRUE;
173 }
174 Int_t okD0=0,okD0bar=0;
12bfc069 175 if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
3a219f60 176 // get daughter AOD tracks
177 AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
178 AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
179
180 lab0 = trk0->GetLabel();
181 lab1 = trk1->GetLabel();
182
183
184 AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(lab0);
185 if(!part0) {
186 printf("no MC particle\n");
187 continue;
188 }
189 while(part0->GetMother()>=0) {
190 labMother=part0->GetMother();
191 part0 = (AliAODMCParticle*)mcArray->At(labMother);
192 if(!part0) {
193 printf("no MC mother particle\n");
194 break;
195 }
196 pdgMother = TMath::Abs(part0->GetPdgCode());
197 if(pdgMother==421) {
198 labD0daugh0=labMother;
199 break;
200 }
201 }
202
203 AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(lab1);
204 if(!part1) {
205 printf("no MC particle\n");
206 continue;
207 }
208 while(part1->GetMother()>=0) {
209 labMother=part1->GetMother();
210 part1 = (AliAODMCParticle*)mcArray->At(labMother);
211 if(!part1) {
212 printf("no MC mother particle\n");
213 break;
214 }
215 pdgMother = TMath::Abs(part1->GetPdgCode());
216 if(pdgMother==421) {
217 labD0daugh1=labMother;
218 break;
219 }
220 }
221
222 // check if the candidate is signal
223 if(labD0daugh0>=0 && labD0daugh1>=0 && labD0daugh0==labD0daugh1) {
224
225 AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0daugh0);
12bfc069 226 // check that the D0 decays in 2 prongs
227 if (TMath::Abs(partD0->GetDaughter(1)-partD0->GetDaughter(0))==1) {
228
229 pdgD0 = partD0->GetPdgCode();
230 Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
231 fHistMass->Fill(invmass);
3a219f60 232
12bfc069 233 // Post the data already here
234 PostData(1,fOutput);
3a219f60 235
12bfc069 236 fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
237 }
3a219f60 238 }
239
240
241 }
242 if(unsetvtx) d->UnsetOwnPrimaryVtx();
243 } // end loop on D0->Kpi
244
245
246 return;
247}
248
249//________________________________________________________________________
250void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
251{
252 // Terminate analysis
253 //
254 if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
255
256 fOutput = dynamic_cast<TList*> (GetOutputData(1));
257 if (!fOutput) {
258 printf("ERROR: fOutput not available\n");
259 return;
260 }
261
262 fNtupleD0Cmp = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleD0Cmp"));
263 fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
264
265 return;
266}
267