3f2558626b51f56a166bc30f2df5435c7661ae27
[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 <TString.h>
27 #include <TH1F.h>
28 #include "AliAODEvent.h"
29 #include "AliAODVertex.h"
30 #include "AliAODTrack.h"
31 #include "AliAODMCHeader.h"
32 #include "AliAODMCParticle.h"
33 #include "AliAODRecoDecayHF3Prong.h"
34 #include "AliAnalysisVertexingHF.h"
35 #include "AliAnalysisTaskSE.h"
36 #include "AliAnalysisTaskSEDplus.h"
37
38 ClassImp(AliAnalysisTaskSEDplus)
39
40
41 //________________________________________________________________________
42 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
43 AliAnalysisTaskSE(),
44 fOutput(0), 
45 fNtupleDplus(0),
46 fNtupleDplusbackg(0),
47 fFillNtuple(kFALSE),
48 fVHF(0)
49 {
50   // Default constructor
51 }
52
53 //________________________________________________________________________
54 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
55 AliAnalysisTaskSE(name),
56 fOutput(0), 
57 fNtupleDplus(0),
58 fNtupleDplusbackg(0),
59 fFillNtuple(fillNtuple),
60 fVHF(0)
61 {
62   // Default constructor
63
64   // Output slot #1 writes into a TList container
65   DefineOutput(1,TList::Class());  //My private output
66
67   if(fFillNtuple){
68     // Output slot #2 writes into a TNtuple container
69     DefineOutput(2,TNtuple::Class());  //My private output
70     // Output slot #3 writes into a TNtuple container
71     DefineOutput(3,TNtuple::Class());  //My private output
72   }
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   fOutput->SetName("OutputHistos");
115
116   Int_t nPtBins=4;
117   TString hisname;
118   for(Int_t i=0;i<nPtBins;i++){
119     hisname.Form("hMassPt%d",i);
120     TH1F* hm=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
121     hm->Sumw2();
122     hisname.Form("hSigPt%d",i);
123     TH1F* hs=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
124     hs->Sumw2();
125     hisname.Form("hBkgPt%d",i);
126     TH1F* hb=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
127     hb->Sumw2();
128
129     hisname.Form("hMassPt%dTC",i);
130     TH1F* hmtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
131     hmtc->Sumw2();
132     hisname.Form("hSigPt%dTC",i);
133     TH1F* hstc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
134     hstc->Sumw2();
135     hisname.Form("hBkgPt%dTC",i);
136     TH1F* hbtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
137     hbtc->Sumw2();
138
139
140     fOutput->Add(hm);
141     fOutput->Add(hs);
142     fOutput->Add(hb);
143     fOutput->Add(hmtc);
144     fOutput->Add(hstc);
145     fOutput->Add(hbtc);
146   }
147
148   
149   if(fFillNtuple){
150     OpenFile(2); // 2 is the slot number of the ntuple
151     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");
152     OpenFile(3); // 3 is the slot number of the ntuple
153     fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
154   }
155
156   return;
157 }
158
159 //________________________________________________________________________
160 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
161 {
162   // Execute analysis for current event:
163   // heavy flavor candidates association to MC truth
164
165   
166   
167   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
168
169   // load Dplus->Kpipi candidates                                                   
170   TClonesArray *array3Prong =
171     (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
172   if(!array3Prong) {
173     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
174     return;
175   }
176
177   // AOD primary vertex
178   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
179   //    vtx1->Print();
180   
181   // load MC particles
182   TClonesArray *arrayMC = 
183     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
184   if(!arrayMC) {
185     printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
186     return;
187   }
188     
189   // load MC header
190   AliAODMCHeader *mcHeader = 
191     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
192   if(!mcHeader) {
193     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
194     return;
195   }
196   Int_t n3Prong = array3Prong->GetEntriesFast();
197   if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
198
199
200   Int_t pdgDgDplustoKpipi[3]={321,211,211};
201   Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
202   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
203     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
204     
205     
206     Bool_t unsetvtx=kFALSE;
207     if(!d->GetOwnPrimaryVtx()){
208       d->SetOwnPrimaryVtx(vtx1);
209       unsetvtx=kTRUE;
210     }
211     if(d->SelectDplus(fVHF->GetDplusCuts())) {
212       Int_t iPtBin=0;
213       Double_t ptCand = d->Pt();
214       if(ptCand<2.){
215         iPtBin=0;
216         cutsDplus[7]=0.08;
217         cutsDplus[8]=0.5;
218         cutsDplus[10]=0.979;
219       }
220       else if(ptCand>2. && ptCand<3){ 
221         iPtBin=1;
222         cutsDplus[7]=0.08;
223         cutsDplus[8]=0.5;
224         cutsDplus[9]=0.991;
225       }else if(ptCand>3. && ptCand<5){ 
226         iPtBin=2;
227         cutsDplus[7]=0.1;
228         cutsDplus[8]=0.5;
229         cutsDplus[9]=0.9955;
230       }else{
231         iPtBin=3;
232         cutsDplus[7]=0.1;
233         cutsDplus[8]=0.5;
234         cutsDplus[9]=0.997;
235       }
236       Bool_t passTightCuts=d->SelectDplus(cutsDplus);
237       Int_t labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
238       Double_t invMass=d->InvMassDplus();
239
240       TString hisNameA(Form("hMassPt%d",iPtBin));
241       TString hisNameS(Form("hSigPt%d",iPtBin));
242       TString hisNameB(Form("hBkgPt%d",iPtBin));
243       TString hisNameATC(Form("hMassPt%dTC",iPtBin));
244       TString hisNameSTC(Form("hSigPt%dTC",iPtBin));
245       TString hisNameBTC(Form("hBkgPt%dTC",iPtBin));
246       
247       ((TH1F*)(fOutput->FindObject(hisNameA)))->Fill(invMass);
248       if(passTightCuts){
249         ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
250       }
251
252
253       if(labDp>=0) {
254         ((TH1F*)(fOutput->FindObject(hisNameS)))->Fill(invMass);
255         if(passTightCuts){
256           ((TH1F*)(fOutput->FindObject(hisNameSTC)))->Fill(invMass);
257         }
258         if(fFillNtuple){
259           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
260           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));         
261           Float_t tmp[20];
262           tmp[0]=TMath::Abs(partDp->GetPdgCode());
263           tmp[1]=partDp->Px()-d->Px();
264           tmp[2]=partDp->Py()-d->Py();
265           tmp[3]=partDp->Pz()-d->Pz();
266           tmp[4]=d->PtProng(0);
267           tmp[5]=d->PtProng(1);
268           tmp[6]=d->PtProng(2);
269           tmp[7]=d->Pt();
270           tmp[8]=partDp->Pt();
271           tmp[9]=d->CosPointingAngle();
272           tmp[10]=d->DecayLength();
273           tmp[11]=dg0->Xv();
274           tmp[12]=d->Xv();
275           tmp[13]=d->Yv();
276           tmp[14]=d->Zv();
277           tmp[15]=d->InvMassDplus();
278           tmp[16]=d->GetSigmaVert();
279           tmp[17]=d->Getd0Prong(0);
280           tmp[18]=d->Getd0Prong(1);
281           tmp[19]=d->Getd0Prong(2);       
282           fNtupleDplus->Fill(tmp);
283           PostData(2,fNtupleDplus);
284         }
285       }else{
286         ((TH1F*)(fOutput->FindObject(hisNameB)))->Fill(invMass);
287         if(passTightCuts){
288           ((TH1F*)(fOutput->FindObject(hisNameBTC)))->Fill(invMass);
289         }
290         if(fFillNtuple){
291           Float_t tmpbkg[14];
292           tmpbkg[0]=d->PtProng(0);
293           tmpbkg[1]=d->PtProng(1);
294           tmpbkg[2]=d->PtProng(2);
295           tmpbkg[3]=d->Pt();
296           tmpbkg[4]=d->CosPointingAngle();
297           tmpbkg[5]=d->DecayLength();
298           tmpbkg[6]=d->Xv();
299           tmpbkg[7]=d->Yv();
300           tmpbkg[8]=d->Zv();
301           tmpbkg[9]=d->InvMassDplus();
302           tmpbkg[10]=d->GetSigmaVert();
303           tmpbkg[11]=d->Getd0Prong(0);
304           tmpbkg[12]=d->Getd0Prong(1);
305           tmpbkg[13]=d->Getd0Prong(2);
306           fNtupleDplusbackg->Fill(tmpbkg);          
307           PostData(3,fNtupleDplusbackg);
308         }
309       }
310       PostData(1,fOutput);    
311     }
312     if(unsetvtx) d->UnsetOwnPrimaryVtx();
313   }
314   
315   return;
316 }
317
318 //________________________________________________________________________
319 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
320 {
321   // Terminate analysis
322   //
323   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
324
325   fOutput = dynamic_cast<TList*> (GetOutputData(1));
326   if (!fOutput) {     
327     printf("ERROR: fOutput not available\n");
328     return;
329   }
330
331   if(fFillNtuple){
332     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
333     fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
334   }
335
336   return;
337 }
338