Use new method AliAODRecoDecay::MatchToMC that simplifies matching to MC
[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 fNtupleD0Cmp(0),
48 fHistMass(0),
49 fVHF(0)
50 {
51   // Default constructor
52
53   // Output slot #1 writes into a TList container
54   DefineOutput(1,TList::Class());  //My private output
55 }
56
57 //________________________________________________________________________
58 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
59 AliAnalysisTaskSE(name),
60 fOutput(0), 
61 fNtupleD0Cmp(0),
62 fHistMass(0),
63 fVHF(0)
64 {
65   // Default constructor
66
67   // Output slot #1 writes into a TList container
68   DefineOutput(1,TList::Class());  //My private output
69 }
70
71 //________________________________________________________________________
72 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
73 {
74   // Destructor
75   if (fOutput) {
76     delete fOutput;
77     fOutput = 0;
78   }
79   if (fVHF) {
80     delete fVHF;
81     fVHF = 0;
82   }
83 }  
84
85 //________________________________________________________________________
86 void AliAnalysisTaskSECompareHF::Init()
87 {
88   // Initialization
89
90   if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
91
92   gROOT->LoadMacro("ConfigVertexingHF.C");
93
94   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
95   fVHF->PrintStatus();
96
97   return;
98 }
99
100 //________________________________________________________________________
101 void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
102 {
103   // Create the output container
104   //
105   if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
106
107   // Several histograms are more conveniently managed in a TList
108   fOutput = new TList();
109   fOutput->SetOwner();
110
111   fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
112   fHistMass->Sumw2();
113   fHistMass->SetMinimum(0);
114   fOutput->Add(fHistMass);
115
116   fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
117   fOutput->Add(fNtupleD0Cmp);
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   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
129  
130   // load D0->Kpi candidates                                                   
131   TClonesArray *inputArrayD0toKpi =
132     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
133   if(!inputArrayD0toKpi) {
134     printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
135     return;
136   }
137
138   /*
139   // load D*+ candidates                                                   
140   TClonesArray *inputArrayDstar =
141     (TClonesArray*)aod->GetList()->FindObject("Dstar");
142   if(!inputArrayDstar) {
143     printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
144     return;
145   }
146   */
147
148   // AOD primary vertex
149   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
150   //vtx1->Print();
151
152   // load MC particles
153   TClonesArray *mcArray = 
154     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
155   if(!mcArray) {
156     printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
157     return;
158   }
159
160   // load MC header
161   AliAODMCHeader *mcHeader = 
162     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
163   if(!mcHeader) {
164     printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
165     return;
166   }
167   
168     
169   // loop over D0->Kpi candidates
170   Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
171   printf("Number of D0->Kpi: %d\n",nInD0toKpi);
172
173   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
174     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
175     Bool_t unsetvtx=kFALSE;
176     if(!d->GetOwnPrimaryVtx()) {
177       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
178       unsetvtx=kTRUE;
179     }
180     Int_t okD0=0,okD0bar=0; 
181     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
182       Int_t labD0 = d->MatchToMC(421,mcArray);
183       if(labD0>=0) {
184         AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0);
185         Int_t pdgD0 = partD0->GetPdgCode();
186         Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
187         fHistMass->Fill(invmass);
188         // Post the data already here
189         PostData(1,fOutput);
190
191         fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
192       }
193     }
194     if(unsetvtx) d->UnsetOwnPrimaryVtx();
195   } // end loop on D0->Kpi
196
197
198   // SIMPLE EXAMPLE OF D*+ MATCHING
199   /*
200   // loop over D*+ candidates
201   for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
202     AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
203     Int_t labDstar = c->MatchToMC(413,421,mcArray);
204     if(labDstar>=0) printf("GOOD MATCH FOR D*+\n"); 
205   }
206   */
207
208   return;
209 }
210
211 //________________________________________________________________________
212 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
213 {
214   // Terminate analysis
215   //
216   if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
217
218   fOutput = dynamic_cast<TList*> (GetOutputData(1));
219   if (!fOutput) {     
220     printf("ERROR: fOutput not available\n");
221     return;
222   }
223
224   fNtupleD0Cmp = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleD0Cmp"));
225   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
226
227   return;
228 }
229