Fill ntuple only on request; more mass histos (Francesco, Renu)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.cxx
CommitLineData
d486095a 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.
8c34bd86 20// Author: Renu Bala, bala@to.infn.it
d486095a 21/////////////////////////////////////////////////////////////
22
23#include <TClonesArray.h>
24#include <TNtuple.h>
25#include <TList.h>
1f4e9722 26#include <TString.h>
d486095a 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
38ClassImp(AliAnalysisTaskSEDplus)
39
40
41//________________________________________________________________________
42AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
43AliAnalysisTaskSE(),
44fOutput(0),
45fNtupleDplus(0),
46fNtupleDplusbackg(0),
1f4e9722 47fFillNtuple(kFALSE),
d486095a 48fVHF(0)
49{
50 // Default constructor
d486095a 51}
52
53//________________________________________________________________________
1f4e9722 54AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
d486095a 55AliAnalysisTaskSE(name),
56fOutput(0),
57fNtupleDplus(0),
58fNtupleDplusbackg(0),
1f4e9722 59fFillNtuple(fillNtuple),
d486095a 60fVHF(0)
61{
62 // Default constructor
63
64 // Output slot #1 writes into a TList container
65 DefineOutput(1,TList::Class()); //My private output
1f4e9722 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 }
d486095a 73}
74
75//________________________________________________________________________
76AliAnalysisTaskSEDplus::~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//________________________________________________________________________
90void 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//________________________________________________________________________
105void 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();
1f4e9722 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 }
d486095a 147
8c34bd86 148
1f4e9722 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
d486095a 156 return;
157}
158
159//________________________________________________________________________
160void 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
1f4e9722 178 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
179 // vtx1->Print();
180
d486095a 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 }
1f4e9722 188
d486095a 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 }
1f4e9722 196 Int_t n3Prong = array3Prong->GetEntriesFast();
197 if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
d486095a 198
199
1f4e9722 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 }
d486095a 251
252
1f4e9722 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){
8c34bd86 259 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
1f4e9722 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);
8c34bd86 306 fNtupleDplusbackg->Fill(tmpbkg);
307 PostData(3,fNtupleDplusbackg);
fc8d975b 308 }
8c34bd86 309 }
1f4e9722 310 PostData(1,fOutput);
d486095a 311 }
1f4e9722 312 if(unsetvtx) d->UnsetOwnPrimaryVtx();
313 }
314
d486095a 315 return;
316}
317
318//________________________________________________________________________
319void 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
1f4e9722 331 if(fFillNtuple){
332 fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
333 fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
334 }
fc8d975b 335
8c34bd86 336 return;
d486095a 337}
338