c1ddc39c953301d54cc3f8d9309286808413e262
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.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 extraction of signal(e.g D+) of heavy flavor
19 // decay candidates with the MC truth.
20 // Author: Renu Bala, bala@to.infn.it
21 /////////////////////////////////////////////////////////////
22
23 #include <TClonesArray.h>
24 #include <TNtuple.h>
25 #include <TList.h>
26 #include <TH1F.h>
27 #include "AliAODEvent.h"
28 #include "AliAODVertex.h"
29 #include "AliAODTrack.h"
30 #include "AliAODMCHeader.h"
31 #include "AliAODMCParticle.h"
32 #include "AliAODRecoDecayHF3Prong.h"
33 #include "AliAnalysisVertexingHF.h"
34 #include "AliAnalysisTaskSE.h"
35 #include "AliAnalysisTaskSEDplus.h"
36
37 ClassImp(AliAnalysisTaskSEDplus)
38
39
40 //________________________________________________________________________
41 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
42 AliAnalysisTaskSE(),
43 fOutput(0), 
44 fNtupleDplus(0),
45 fNtupleDplusbackg(0),
46 fHistMass(0),
47 fHistSignal(0),
48 fHistBackground(0),
49 fVHF(0)
50 {
51   // Default constructor
52 }
53
54 //________________________________________________________________________
55 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name):
56 AliAnalysisTaskSE(name),
57 fOutput(0), 
58 fNtupleDplus(0),
59 fNtupleDplusbackg(0),
60 fHistMass(0),
61 fHistSignal(0),
62 fHistBackground(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   // Output slot #2 writes into a TNtuple container
70   DefineOutput(2,TNtuple::Class());  //My private output
71   // Output slot #3 writes into a TNtuple container
72   DefineOutput(3,TNtuple::Class());  //My private output
73 }
74
75 //________________________________________________________________________
76 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
77 {
78   // Destructor
79   if (fOutput) {
80     delete fOutput;
81     fOutput = 0;
82   }
83   if (fVHF) {
84     delete fVHF;
85     fVHF = 0;
86   }
87 }  
88
89 //________________________________________________________________________
90 void AliAnalysisTaskSEDplus::Init()
91 {
92   // Initialization
93
94   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
95
96   gROOT->LoadMacro("ConfigVertexingHF.C");
97
98   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
99   fVHF->PrintStatus();
100
101   return;
102 }
103
104 //________________________________________________________________________
105 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
106 {
107   // Create the output container
108   //
109   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
110
111   // Several histograms are more conveniently managed in a TList
112   fOutput = new TList();
113   fOutput->SetOwner();
114
115   fHistMass = new TH1F("fHistMass", "D^{+} invariant mass; M [GeV]; Entries",100,1.765,1.965);
116   fHistSignal = new TH1F("fHistSignal", "D^{+} invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
117   fHistBackground =new TH1F("fHistBackground", "Background invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
118
119
120
121
122   fHistMass->Sumw2();
123   fHistSignal->Sumw2();
124   fHistBackground->Sumw2();
125
126
127   //fHistMass->SetMinimum(0);
128   //fHistSignal->SetMinimum(0);
129   //fHistBackground->SetMinimum(0);
130   fOutput->Add(fHistMass);
131   fOutput->Add(fHistSignal);
132   fOutput->Add(fHistBackground);
133   
134   OpenFile(2); // 2 is the slot number of the ntuple
135   fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Ptpi:PtK:Ptpi2:PtRec:PtTrue:PointingAngle:DecLeng:VxTrue:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2");
136   OpenFile(3); // 3 is the slot number of the ntuple
137   fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
138   
139   return;
140 }
141
142 //________________________________________________________________________
143 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
144 {
145   // Execute analysis for current event:
146   // heavy flavor candidates association to MC truth
147
148   
149   
150   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
151
152   // load Dplus->Kpipi candidates                                                   
153   TClonesArray *array3Prong =
154     (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
155   if(!array3Prong) {
156     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
157     return;
158   }
159
160   // AOD primary vertex
161     AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
162     //    vtx1->Print();
163
164   // load MC particles
165   TClonesArray *arrayMC = 
166     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
167   if(!arrayMC) {
168     printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
169     return;
170   }
171
172   // load MC header
173   AliAODMCHeader *mcHeader = 
174     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
175   if(!mcHeader) {
176     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
177     return;
178   }
179     Int_t n3Prong = array3Prong->GetEntriesFast();
180     if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
181
182
183     Int_t pdgDgDplustoKpipi[3]={321,211,211};
184
185     for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
186       AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
187
188
189       Bool_t unsetvtx=kFALSE;
190       if(!d->GetOwnPrimaryVtx()){
191         d->SetOwnPrimaryVtx(vtx1);
192         unsetvtx=kTRUE;
193       }
194       if(d->SelectDplus(fVHF->GetDplusCuts())) {
195           
196           
197         Int_t labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
198     
199         if(labDp>=0) {
200           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
201           Int_t pdgDp = TMath::Abs(partDp->GetPdgCode());
202           if(pdgDp==411){
203             
204             fHistSignal->Fill(d->InvMassDplus());
205             fHistMass->Fill(d->InvMassDplus());
206             
207             // Post the data already here
208             PostData(1,fOutput);
209             
210             AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
211
212             Float_t tmp[20];
213             tmp[0]=pdgDp;
214             tmp[1]=partDp->Px()-d->Px();
215             tmp[2]=partDp->Py()-d->Py();
216             tmp[3]=partDp->Pz()-d->Pz();
217             tmp[4]=d->PtProng(0);
218             tmp[5]=d->PtProng(1);
219             tmp[6]=d->PtProng(2);
220             tmp[7]=d->Pt();
221             tmp[8]=partDp->Pt();
222             tmp[9]=d->CosPointingAngle();
223             tmp[10]=d->DecayLength();
224             tmp[11]=dg0->Xv();
225             tmp[12]=d->Xv();
226             tmp[13]=d->Yv();
227             tmp[14]=d->Zv();
228             tmp[15]=d->InvMassDplus();
229             tmp[16]=d->GetSigmaVert();
230             tmp[17]=d->Getd0Prong(0);
231             tmp[18]=d->Getd0Prong(1);
232             tmp[19]=d->Getd0Prong(2);
233
234
235             fNtupleDplus->Fill(tmp);
236             PostData(2,fNtupleDplus);
237           }
238         } else {     
239
240             Float_t tmpbkg[14];
241             tmpbkg[0]=d->PtProng(0);
242             tmpbkg[1]=d->PtProng(1);
243             tmpbkg[2]=d->PtProng(2);
244             tmpbkg[3]=d->Pt();
245             tmpbkg[4]=d->CosPointingAngle();
246             tmpbkg[5]=d->DecayLength();
247             tmpbkg[6]=d->Xv();
248             tmpbkg[7]=d->Yv();
249             tmpbkg[8]=d->Zv();
250             tmpbkg[9]=d->InvMassDplus();
251             tmpbkg[10]=d->GetSigmaVert();
252             tmpbkg[11]=d->Getd0Prong(0);
253             tmpbkg[12]=d->Getd0Prong(1);
254             tmpbkg[13]=d->Getd0Prong(2);
255
256           fHistBackground->Fill(d->InvMassDplus());
257           fNtupleDplusbackg->Fill(tmpbkg);          
258           PostData(3,fNtupleDplusbackg);
259           fHistMass->Fill(d->InvMassDplus());
260         }
261         
262         
263       }
264
265       if(unsetvtx) d->UnsetOwnPrimaryVtx();
266
267     }
268      // end loop on D+->Kpipi
269
270
271   return;
272 }
273
274 //________________________________________________________________________
275 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
276 {
277   // Terminate analysis
278   //
279   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
280
281   fOutput = dynamic_cast<TList*> (GetOutputData(1));
282   if (!fOutput) {     
283     printf("ERROR: fOutput not available\n");
284     return;
285   }
286
287   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
288   fHistSignal = dynamic_cast<TH1F*>(fOutput->FindObject("fHistSignal"));
289   fHistBackground = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBackground"));
290   fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
291   fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
292
293   return;
294 }
295