Included also 3 and 4 prong decays; added variables to ntuple
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSECompareHF.cxx
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 "AliAODRecoCascadeHF.h"
36 #include "AliAnalysisVertexingHF.h"
37 #include "AliAnalysisTaskSE.h"
38 #include "AliAnalysisTaskSECompareHF.h"
39
40 ClassImp(AliAnalysisTaskSECompareHF)
41
42
43 //________________________________________________________________________
44 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
45 AliAnalysisTaskSE(),
46 fOutput(0), 
47 fNtupleCmp(0),
48 fHistMass(0),
49 fVHF(0)
50 {
51   // Default constructor
52
53   // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
54 }
55
56 //________________________________________________________________________
57 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
58 AliAnalysisTaskSE(name),
59 fOutput(0), 
60 fNtupleCmp(0),
61 fHistMass(0),
62 fVHF(0)
63 {
64   // Standard constructor
65
66   // Output slot #1 writes into a TList container
67   DefineOutput(1,TList::Class());  //My private output
68   // Output slot #2 writes into a TNtuple container
69   DefineOutput(2,TNtuple::Class());  //My private output
70 }
71
72 //________________________________________________________________________
73 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
74 {
75   // Destructor
76   if (fOutput) {
77     delete fOutput;
78     fOutput = 0;
79   }
80   if (fVHF) {
81     delete fVHF;
82     fVHF = 0;
83   }
84 }  
85
86 //________________________________________________________________________
87 void AliAnalysisTaskSECompareHF::Init()
88 {
89   // Initialization
90
91   if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
92   
93   gROOT->LoadMacro("ConfigVertexingHF.C");
94
95   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
96   fVHF->PrintStatus();
97   
98   return;
99 }
100
101 //________________________________________________________________________
102 void 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
117   fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec");
118
119   return;
120 }
121
122 //________________________________________________________________________
123 void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
124 {
125   // Execute analysis for current event:
126   // heavy flavor candidates association to MC truth
127
128   
129   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
130   
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
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
148  
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   }
156   
157
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   
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();
222       }
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;
280     }
281
282   } // end loop on vertices
283
284
285   // loop over D*+ candidates
286   /*
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   }
292   */
293
294   // Post the data already here
295   PostData(1,fOutput);
296   PostData(2,fNtupleCmp);
297
298   return;
299 }
300 //________________________________________________________________________
301 void 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
313   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
314
315   fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
316
317   return;
318 }
319