New task for D0 mass distributions (Chiara B)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSED0Mass.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 D0 candidates invariant mass histogram
19 // and comparison with the MC truth.
20 //
21 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
22 // and Chiara Bianchin, chiara.bianchin@pd.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <Riostream.h>
26 #include <TClonesArray.h>
27 #include <TNtuple.h>
28 #include <TList.h>
29 #include <TH1F.h>
30
31
32 #include "AliAODEvent.h"
33 #include "AliAODVertex.h"
34 #include "AliAODTrack.h"
35 #include "AliAODMCHeader.h"
36 #include "AliAODMCParticle.h"
37 #include "AliAODRecoDecayHF2Prong.h"
38 #include "AliAODRecoCascadeHF.h"
39 #include "AliAnalysisVertexingHF.h"
40 #include "AliAnalysisTaskSE.h"
41 #include "AliAnalysisTaskSED0Mass.h"
42
43 ClassImp(AliAnalysisTaskSED0Mass)
44
45
46 //________________________________________________________________________
47 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
48 AliAnalysisTaskSE(),
49 fOutput(0), 
50 //fNtupleD0Cmp(0),
51 fhistMass(0),
52 fhistSgn(0),
53 fhistBkg(0),
54 fVHF(0)
55 {
56   // Default constructor
57 }
58
59 //________________________________________________________________________
60 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name):
61 AliAnalysisTaskSE(name),
62 fOutput(0), 
63 //fNtupleD0Cmp(0),
64 fhistMass(0),
65 fhistSgn(0),
66 fhistBkg(0),
67 fVHF(0)
68 {
69   // Default constructor
70
71   // Output slot #1 writes into a TList container
72   DefineOutput(1,TList::Class());  //My private output
73 }
74
75 //________________________________________________________________________
76 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
77 {
78   // Destructor
79
80   if (fhistMass){
81     delete fhistMass;
82     fhistMass=0;
83   }
84
85   if (fhistSgn){
86     delete fhistSgn;
87     fhistSgn=0;
88   }
89
90   if (fhistBkg){
91     delete fhistBkg;
92     fhistBkg=0;
93   }
94
95   if (fOutput) {
96     delete fOutput;
97     fOutput = 0;
98   }
99   if (fVHF) {
100     delete fVHF;
101     fVHF = 0;
102   }
103
104 }  
105
106 //________________________________________________________________________
107 void AliAnalysisTaskSED0Mass::Init()
108 {
109   // Initialization
110
111   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
112
113   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
114
115   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
116   // set dedidcated cuts
117   fVHF->SetD0toKpiCuts(0.7,0.03,0.8,0.06,0.06,0.05,0.05,-0.0002,0.6); //2.p-p vertex reconstructed
118   fVHF->PrintStatus();
119
120   return;
121 }
122
123 //________________________________________________________________________
124 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
125 {
126
127   // Create the output container
128   //
129   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
130
131   // Several histograms are more conveniently managed in a TList
132   fOutput = new TList();
133   fOutput->SetOwner();
134
135   const Int_t nhist=3;
136   fhistMass=new TClonesArray("TH1F",nhist);
137   fhistMass->SetName("fhistMass");
138   fhistSgn =new TClonesArray("TH1F",nhist);
139   fhistSgn->SetName("fhistSgn");
140   fhistBkg =new TClonesArray("TH1F",nhist);
141   fhistBkg->SetName("fhistBkg");
142
143   TString nameMass, nameSgn, nameBkg;
144
145   for(Int_t i=0;i<nhist;i++){
146     nameMass="fhistMass_";
147     nameMass+=i+1;
148     cout<<"Creating TH1F with name "<<nameMass<<endl;
149     new ((*fhistMass)[i]) TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
150     nameSgn="fhistSgn_";
151     nameSgn+=i+1;
152     new ((*fhistSgn)[i]) TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
153     nameBkg="fhistBkg_";
154     nameBkg+=i+1;
155     new ((*fhistBkg)[i]) TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
156     printf("Created histograms %s\t%s\t%s\n",fhistMass->UncheckedAt(i)->GetName(),fhistSgn->UncheckedAt(i)->GetName(),fhistBkg->UncheckedAt(i)->GetName());
157   }
158  
159   fOutput->Add(fhistMass);
160   fOutput->Add(fhistSgn);
161   fOutput->Add(fhistBkg);
162
163   return;
164 }
165
166 //________________________________________________________________________
167 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
168 {
169   // Execute analysis for current event:
170   // heavy flavor candidates association to MC truth
171   
172   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
173  
174   // load D0->Kpi candidates                                                   
175   TClonesArray *inputArrayD0toKpi =
176     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
177   if(!inputArrayD0toKpi) {
178     printf("AliAnalysisTaskSECompareHFpt::UserExec: D0toKpi branch not found!\n");
179     return;
180   }
181   
182   // AOD primary vertex
183   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
184   //vtx1->Print();
185   
186   // load MC particles
187   TClonesArray *mcArray = 
188     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
189   if(!mcArray) {
190     printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
191     return;
192   }
193   
194   // load MC header
195   AliAODMCHeader *mcHeader = 
196     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
197   if(!mcHeader) {
198     printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
199     return;
200   }
201   
202   
203   // loop over D0->Kpi candidates
204   Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast();
205   printf("Number of D0->Kpi: %d\n",nInD0toKpi);
206   
207   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
208     
209     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi);
210     Bool_t unsetvtx=kFALSE;
211     if(!d->GetOwnPrimaryVtx()) {
212       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
213       unsetvtx=kTRUE;
214     }
215     
216     //  cout<<"Cuts applied: ";
217     //     Double_t *cutsapplied=(Double_t*)fVHF->GetD0toKpiCuts();
218     //     for (Int_t i=0;i<9;i++){
219     //       printf(cutsapplied[i]\t);
220     //     }
221     //    cout<<endl;
222     
223     
224     Double_t pt = d->Pt();
225     Int_t ptbin=0;
226     
227     //cout<<"P_t = "<<pt<<endl;
228     if (pt>0. && pt<=1.) {
229       ptbin=1; 
230       //printf("I'm in the bin %d\n",ptbin);
231     }
232     else {
233       if(pt>1. && pt<=3.) {
234         ptbin=2;  
235         //printf("I'm in the bin %d\n",ptbin);
236       }
237       else {
238         ptbin=3;  
239         //printf("I'm in the bin %d\n",ptbin);
240       }//if(pt>3)
241     }
242     
243     FillHists(ptbin,d,mcArray);
244    
245     if(unsetvtx) d->UnsetOwnPrimaryVtx();
246  
247   }
248      
249   
250    
251   // Post the data
252   PostData(1,fOutput);
253
254
255   return;
256 }
257 //____________________________________________________________________________*
258 void AliAnalysisTaskSED0Mass::FillHists(Int_t ptbin, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC){
259   //
260   // function used in UserExec:
261   //
262   Int_t okD0=0,okD0bar=0;
263
264   if(part->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {//selected
265     Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
266     //printf("SELECTED\n");
267     Int_t labD0 = part->MatchToMC(421,arrMC); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
268     //printf("labD0 %d",labD0);
269
270     if(labD0>=0) {
271       
272       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
273
274       Int_t pdgD0 = partD0->GetPdgCode();
275       
276       //printf(" pdgD0 %d\n",pdgD0);
277       
278       if (pdgD0==421){ //D0
279         //cout<<"Fill S with D0"<<endl;
280         ((TH1F*)(fhistSgn->UncheckedAt(ptbin-1)))->Fill(invmassD0);
281       }
282       else {//D0bar  
283         //printf("Fill S with D0bar");
284         ((TH1F*)(fhistSgn->UncheckedAt(ptbin-1)))->Fill(invmassD0bar);
285       }
286
287     }
288     else {//background
289           //printf("Fill background");
290       ((TH1F*)(fhistBkg->UncheckedAt(ptbin-1)))->Fill(invmassD0);
291       ((TH1F*)(fhistBkg->UncheckedAt(ptbin-1)))->Fill(invmassD0bar);
292     }
293     
294     //no MC info, just cut selection
295     if (okD0==1) {
296       //printf("Fill mass with D0");
297       ((TH1F*)(fhistMass->UncheckedAt(ptbin-1)))->Fill(invmassD0);
298       
299     }
300     if (okD0bar==1) {
301       ((TH1F*)fhistMass->UncheckedAt(ptbin-1))->Fill(invmassD0bar);
302       //printf("Fill mass with D0bar");
303     }
304       
305   } //else cout<<"NOT SELECTED"<<endl;
306
307 }
308 //________________________________________________________________________
309 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
310 {
311   // Terminate analysis
312   //
313   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
314
315   fOutput = dynamic_cast<TList*> (GetOutputData(1));
316   if (!fOutput) {     
317     printf("ERROR: fOutput not available\n");
318     return;
319   }
320
321   
322   fhistMass = dynamic_cast<TClonesArray*>(fOutput->FindObject("fhistMass"));
323   fhistSgn = dynamic_cast<TClonesArray*>(fOutput->FindObject("fhistSgn"));
324   fhistBkg = dynamic_cast<TClonesArray*>(fOutput->FindObject("fhistBkg"));
325  
326   return;
327 }
328