protect function PlotNTrackletsTrack against missing ESD info for usage
[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>
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
37ClassImp(AliAnalysisTaskSEDplus)
38
39
40//________________________________________________________________________
41AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
42AliAnalysisTaskSE(),
43fOutput(0),
44fNtupleDplus(0),
45fNtupleDplusbackg(0),
46fHistMass(0),
47fHistSignal(0),
48fHistBackground(0),
49fVHF(0)
50{
51 // Default constructor
d486095a 52}
53
54//________________________________________________________________________
55AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name):
56AliAnalysisTaskSE(name),
57fOutput(0),
58fNtupleDplus(0),
59fNtupleDplusbackg(0),
60fHistMass(0),
61fHistSignal(0),
62fHistBackground(0),
63fVHF(0)
64{
65 // Default constructor
66
67 // Output slot #1 writes into a TList container
68 DefineOutput(1,TList::Class()); //My private output
10bdd1ae 69 // Output slot #2 writes into a TNtuple container
8c34bd86 70 DefineOutput(2,TNtuple::Class()); //My private output
10bdd1ae 71 // Output slot #3 writes into a TNtuple container
8c34bd86 72 DefineOutput(3,TNtuple::Class()); //My private output
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();
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
2a5c77c6 134 OpenFile(2); // 2 is the slot number of the ntuple
8332db59 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");
2a5c77c6 136 OpenFile(3); // 3 is the slot number of the ntuple
8c34bd86 137 fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
138
d486095a 139 return;
140}
141
142//________________________________________________________________________
143void 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();
3b1598da 180 if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
d486095a 181
182
c8112d2f 183 Int_t pdgDgDplustoKpipi[3]={321,211,211};
184
d486095a 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 }
8c34bd86 194 if(d->SelectDplus(fVHF->GetDplusCuts())) {
d486095a 195
196
c8112d2f 197 Int_t labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
8c34bd86 198
199 if(labDp>=0) {
200 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
201 Int_t pdgDp = TMath::Abs(partDp->GetPdgCode());
202 if(pdgDp==411){
413e9336 203
8c34bd86 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);
413e9336 259 fHistMass->Fill(d->InvMassDplus());
fc8d975b 260 }
fc8d975b 261
8c34bd86 262
263 }
264
d486095a 265 if(unsetvtx) d->UnsetOwnPrimaryVtx();
266
267 }
268 // end loop on D+->Kpipi
269
270
271 return;
272}
273
274//________________________________________________________________________
275void 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
d486095a 287 fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
288 fHistSignal = dynamic_cast<TH1F*>(fOutput->FindObject("fHistSignal"));
289 fHistBackground = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBackground"));
8c34bd86 290 fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
291 fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
fc8d975b 292
8c34bd86 293 return;
d486095a 294}
295