72153b8e4d9f71a5cee8e4a20271bb9d20059b6e
[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   // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
54 }
55
56 //________________________________________________________________________
57 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
58 AliAnalysisTaskSE(name),
59 fOutput(0), 
60 fNtupleD0Cmp(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   fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue");
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   // load D0->Kpi candidates                                                   
132   TClonesArray *inputArrayD0toKpi =
133     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
134   if(!inputArrayD0toKpi) {
135     printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
136     return;
137   }
138
139  
140   // load D*+ candidates                                                   
141   TClonesArray *inputArrayDstar =
142     (TClonesArray*)aod->GetList()->FindObject("Dstar");
143   if(!inputArrayDstar) {
144     printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
145     return;
146   }
147   
148
149   // AOD primary vertex
150   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
151   //vtx1->Print();
152
153   // load MC particles
154   TClonesArray *mcArray = 
155     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
156   if(!mcArray) {
157     printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
158     return;
159   }
160
161   // load MC header
162   AliAODMCHeader *mcHeader = 
163     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
164   if(!mcHeader) {
165     printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
166     return;
167   }
168   
169     
170   // loop over D0->Kpi candidates
171   Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
172   printf("Number of D0->Kpi: %d\n",nInD0toKpi);
173
174   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
175     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
176     Bool_t unsetvtx=kFALSE;
177     if(!d->GetOwnPrimaryVtx()) {
178       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
179       unsetvtx=kTRUE;
180     }
181     Int_t okD0=0,okD0bar=0; 
182     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
183       Int_t labD0 = d->MatchToMC(421,mcArray);
184       if(labD0>=0) {
185         AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0);
186         Int_t pdgD0 = partD0->GetPdgCode();
187         Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar());
188         fHistMass->Fill(invmass);
189         // Post the data already here
190         PostData(1,fOutput);
191
192         fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt());
193         PostData(2,fNtupleD0Cmp);
194       }
195     }
196     if(unsetvtx) d->UnsetOwnPrimaryVtx();
197   } // end loop on D0->Kpi
198
199
200   // SIMPLE EXAMPLE OF D*+ MATCHING
201
202     
203   // loop over D*+ candidates
204   /*
205   for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
206     AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
207     Int_t labDstar = c->MatchToMC(413,421,mcArray);
208     if(labDstar>=0) printf("GOOD MATCH FOR D*+\n"); 
209   }
210   */
211   // Post the data already here
212   PostData(1,fOutput);
213   PostData(2,fNtupleD0Cmp);
214
215   return;
216 }
217 //________________________________________________________________________
218 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
219 {
220   // Terminate analysis
221   //
222   if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
223
224   fOutput = dynamic_cast<TList*> (GetOutputData(1));
225   if (!fOutput) {     
226     printf("ERROR: fOutput not available\n");
227     return;
228   }
229
230   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
231
232   fNtupleD0Cmp = dynamic_cast<TNtuple*> (GetOutputData(2));
233
234   return;
235 }
236