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