Added LHC10h run list for flow analysis (Giacomo)
[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
27de2dfb 16/* $Id$ */
17
3a219f60 18/////////////////////////////////////////////////////////////
19//
20// AliAnalysisTaskSE for the comparison of heavy flavor
21// decay candidates with the MC truth.
22//
23// Author: A.Dainese, andrea.dainese@lnl.infn.it
24/////////////////////////////////////////////////////////////
25
26#include <TClonesArray.h>
27#include <TNtuple.h>
28#include <TList.h>
29#include <TH1F.h>
30
b557eb43 31#include "AliAnalysisManager.h"
32#include "AliAODHandler.h"
3a219f60 33#include "AliAODEvent.h"
34#include "AliAODVertex.h"
607719af 35#include "AliESDtrack.h"
3a219f60 36#include "AliAODTrack.h"
37#include "AliAODMCHeader.h"
38#include "AliAODMCParticle.h"
39#include "AliAODRecoDecayHF2Prong.h"
c4cdf105 40#include "AliAODRecoDecayHF3Prong.h"
41#include "AliAODRecoDecayHF4Prong.h"
b3999999 42#include "AliAODRecoCascadeHF.h"
2ae9c7cb 43#include "AliAnalysisVertexingHF.h"
3a219f60 44#include "AliAnalysisTaskSE.h"
45#include "AliAnalysisTaskSECompareHF.h"
46
47ClassImp(AliAnalysisTaskSECompareHF)
48
49
50//________________________________________________________________________
51AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
52AliAnalysisTaskSE(),
53fOutput(0),
2a041947 54fNtupleCmp(0),
12bfc069 55fHistMass(0),
8887af0b 56fHistNEvents(0)
3a219f60 57{
58 // Default constructor
3a219f60 59
9d6d35b0 60 // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
3a219f60 61}
62
63//________________________________________________________________________
64AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
65AliAnalysisTaskSE(name),
66fOutput(0),
2a041947 67fNtupleCmp(0),
12bfc069 68fHistMass(0),
8887af0b 69fHistNEvents(0)
3a219f60 70{
9d6d35b0 71 // Standard constructor
3a219f60 72
73 // Output slot #1 writes into a TList container
74 DefineOutput(1,TList::Class()); //My private output
10bdd1ae 75 // Output slot #2 writes into a TNtuple container
76 DefineOutput(2,TNtuple::Class()); //My private output
3a219f60 77}
78
79//________________________________________________________________________
80AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
81{
82 // Destructor
83 if (fOutput) {
84 delete fOutput;
85 fOutput = 0;
86 }
87}
88
89//________________________________________________________________________
90void AliAnalysisTaskSECompareHF::Init()
91{
92 // Initialization
93
94 if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
9d6d35b0 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
a273622f 115 fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
116 fHistNEvents->Sumw2();
117 fHistNEvents->SetMinimum(0);
118 fOutput->Add(fHistNEvents);
119
2a5c77c6 120 OpenFile(2); // 2 is the slot number of the ntuple
607719af 121 fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
3a219f60 122
123 return;
124}
125
126//________________________________________________________________________
127void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
128{
129 // Execute analysis for current event:
130 // heavy flavor candidates association to MC truth
9d6d35b0 131
3a219f60 132
133 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
8931c313 134
b557eb43 135 TClonesArray *inputArrayVertices = 0;
136 TClonesArray *inputArrayD0toKpi = 0;
137 TClonesArray *inputArrayDstar = 0;
138
139 if(!aod && AODEvent() && IsStandardAOD()) {
140 // In case there is an AOD handler writing a standard AOD, use the AOD
141 // event in memory rather than the input (ESD) event.
142 aod = dynamic_cast<AliAODEvent*> (AODEvent());
143 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
144 // have to taken from the AOD event hold by the AliAODExtension
145 AliAODHandler* aodHandler = (AliAODHandler*)
146 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
147 if(aodHandler->GetExtensions()) {
148 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
149 AliAODEvent *aodFromExt = ext->GetAOD();
150 // load HF vertices
151 inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
152 // load D0->Kpi candidates
153 inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
154 // load D*+ candidates
155 inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
156 }
dc222f77 157 } else if(aod) {
b557eb43 158 // load HF vertices
159 inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
160 // load D0->Kpi candidates
161 inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
162 // load D*+ candidates
163 inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
164 }
165
2a041947 166
20c2e031 167 if(!inputArrayVertices || !aod) {
2a041947 168 printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
169 return;
170 }
3a219f60 171 if(!inputArrayD0toKpi) {
172 printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
173 return;
174 }
b3999999 175 if(!inputArrayDstar) {
176 printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
177 return;
178 }
9d6d35b0 179
b3999999 180
b557eb43 181 fHistNEvents->Fill(0); // count event
182 // Post the data already here
183 PostData(1,fOutput);
184
3a219f60 185 // AOD primary vertex
186 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
187 //vtx1->Print();
188
189 // load MC particles
190 TClonesArray *mcArray =
191 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
192 if(!mcArray) {
193 printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
194 return;
195 }
196
197 // load MC header
198 AliAODMCHeader *mcHeader =
199 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
200 if(!mcHeader) {
201 printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
202 return;
203 }
204
2a041947 205 Int_t nprongs,lab,okD0,okD0bar,pdg;
206 Bool_t unsetvtx;
90426506 207 Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
952f1e33 208 AliAODRecoDecayHF2Prong *d2=0;
209 AliAODRecoDecayHF3Prong *d3=0;
210 AliAODRecoDecayHF4Prong *d4=0;
2a041947 211
c8112d2f 212 Int_t pdgDgD0toKpi[2]={321,211};
213 Int_t pdgDgDplustoKpipi[3]={321,211,211};
214 Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
215
2a041947 216 // loop over vertices
217 Int_t nVertices = inputArrayVertices->GetEntriesFast();
3b1598da 218 if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
2a041947 219
220 for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
221 AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
222
90426506 223 vtx->GetXYZ(posRec);
224 vtx->GetCovarianceMatrix(covRec);
225 errx=1.; erry=1.; errz=1.;
226 if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
227 if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
228 if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
229
230
2a041947 231 // get parent AliAODRecoDecayHF
232 nprongs= vtx->GetNDaughters();
607719af 233 // check that the daughters have kITSrefit and kTPCrefit
234 Bool_t allDgOK=kTRUE;
235 for(Int_t idg=0; idg<nprongs; idg++) {
236 AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
237 if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
238 if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
239 }
240 if(!allDgOK) continue;
241
2a041947 242
243 switch(nprongs) {
244 case 2: // look for D0->Kpi
952f1e33 245 d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
27917274 246 if(d2->IsLikeSign()) continue;
a273622f 247 if(d2->Charge() != 0) continue; // these are D*
c8112d2f 248 lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
2a041947 249 if(lab>=0) {
250 unsetvtx=kFALSE;
251 if(!d2->GetOwnPrimaryVtx()) {
252 d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
253 unsetvtx=kTRUE;
254 }
cb8088a2 255 okD0=1; okD0bar=1;
256 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
257 pdg = dMC->GetPdgCode();
258 invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
259 // get a daughter for true pos of decay vertex
260 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
261 dg0MC->XvYvZv(posTrue);
262 fHistMass->Fill(invmass);
263 // Post the data already here
264 PostData(1,fOutput);
265 Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
266 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
267 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
268 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
269 (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
270 (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
271 fNtupleCmp->Fill(tmp);
272 PostData(2,fNtupleCmp);
273
2a041947 274 if(unsetvtx) d2->UnsetOwnPrimaryVtx();
3a219f60 275 }
2a041947 276 break;
277 case 3: // look for D+
952f1e33 278 d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
27917274 279 if(d3->IsLikeSign()) continue;
c8112d2f 280 lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
2a041947 281 if(lab>=0) {
282 unsetvtx=kFALSE;
283 if(!d3->GetOwnPrimaryVtx()) {
284 d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
285 unsetvtx=kTRUE;
286 }
90426506 287 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
288 pdg = dMC->GetPdgCode();
289 invmass = d3->InvMassDplus();
290 // get a daughter for true pos of decay vertex
291 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
292 dg0MC->XvYvZv(posTrue);
607719af 293 Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
294 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
295 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
296 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
297 (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
298 (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
299 fNtupleCmp->Fill(tmp);
a7c1934a 300 PostData(2,fNtupleCmp);
2a041947 301 if(unsetvtx) d3->UnsetOwnPrimaryVtx();
302 }
303 break;
304 case 4: // look for D0->Kpipipi
952f1e33 305 d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
27917274 306 if(d4->IsLikeSign()) continue;
c8112d2f 307 lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
2a041947 308 if(lab>=0) {
309 unsetvtx=kFALSE;
310 if(!d4->GetOwnPrimaryVtx()) {
311 d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
312 unsetvtx=kTRUE;
313 }
314 okD0=0; okD0bar=0;
90426506 315 AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
316 pdg = dMC->GetPdgCode();
317 //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
318 invmass = 10.; // dummy
319 // get a daughter for true pos of decay vertex
320 AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
321 dg0MC->XvYvZv(posTrue);
607719af 322 Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
323 (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
324 (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
325 (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
326 (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
327 (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
328 fNtupleCmp->Fill(tmp);
a7c1934a 329 PostData(2,fNtupleCmp);
2a041947 330 if(unsetvtx) d4->UnsetOwnPrimaryVtx();
331 }
332 break;
3a219f60 333 }
3a219f60 334
2a041947 335 } // end loop on vertices
3a219f60 336
9d6d35b0 337
b3999999 338 // loop over D*+ candidates
4dd17b11 339 /*
b3999999 340 for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
341 AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
342 Int_t labDstar = c->MatchToMC(413,421,mcArray);
343 if(labDstar>=0) printf("GOOD MATCH FOR D*+\n");
344 }
4dd17b11 345 */
2a041947 346
4dd17b11 347 // Post the data already here
348 PostData(1,fOutput);
2a041947 349 PostData(2,fNtupleCmp);
b3999999 350
3a219f60 351 return;
352}
3a219f60 353//________________________________________________________________________
354void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
355{
356 // Terminate analysis
357 //
358 if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
359
360 fOutput = dynamic_cast<TList*> (GetOutputData(1));
361 if (!fOutput) {
362 printf("ERROR: fOutput not available\n");
363 return;
364 }
365
3a219f60 366 fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
a273622f 367 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
3a219f60 368
27917274 369 //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
10bdd1ae 370
3a219f60 371 return;
372}
373