added option to calculate dca relative to TPC vertex
[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.
20
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
fc8d975b 52
d486095a 53}
54
55//________________________________________________________________________
56AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name):
57AliAnalysisTaskSE(name),
58fOutput(0),
59fNtupleDplus(0),
60fNtupleDplusbackg(0),
61fHistMass(0),
62fHistSignal(0),
63fHistBackground(0),
64fVHF(0)
65{
66 // Default constructor
67
68 // Output slot #1 writes into a TList container
69 DefineOutput(1,TList::Class()); //My private output
10bdd1ae 70 // Output slot #2 writes into a TNtuple container
fc8d975b 71 DefineOutput(2,TNtuple::Class()); //My private output //AD
10bdd1ae 72 // Output slot #3 writes into a TNtuple container
fc8d975b 73 DefineOutput(3,TNtuple::Class()); //My private output //AD
d486095a 74}
75
76//________________________________________________________________________
77AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
78{
79 // Destructor
80 if (fOutput) {
81 delete fOutput;
82 fOutput = 0;
83 }
84 if (fVHF) {
85 delete fVHF;
86 fVHF = 0;
87 }
88}
89
90//________________________________________________________________________
91void AliAnalysisTaskSEDplus::Init()
92{
93 // Initialization
94
95 if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
96
97 gROOT->LoadMacro("ConfigVertexingHF.C");
98
99 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
100 fVHF->PrintStatus();
101
102 return;
103}
104
105//________________________________________________________________________
106void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
107{
108 // Create the output container
109 //
110 if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
111
112 // Several histograms are more conveniently managed in a TList
113 fOutput = new TList();
114 fOutput->SetOwner();
115
116 fHistMass = new TH1F("fHistMass", "D^{+} invariant mass; M [GeV]; Entries",100,1.765,1.965);
117 fHistSignal = new TH1F("fHistSignal", "D^{+} invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
118 fHistBackground =new TH1F("fHistBackground", "Background invariant mass - MC; M [GeV]; Entries",100,1.765,1.965);
119
120
121
122
123 fHistMass->Sumw2();
124 fHistSignal->Sumw2();
125 fHistBackground->Sumw2();
126
127
128 //fHistMass->SetMinimum(0);
129 //fHistSignal->SetMinimum(0);
130 //fHistBackground->SetMinimum(0);
131 fOutput->Add(fHistMass);
132 fOutput->Add(fHistSignal);
133 fOutput->Add(fHistBackground);
134
135
136 fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Ptpi:Ptpi2:PtK:PtRec:PtTrue:PointingAngle:DecLeng:VxTrue:VxRec:InvMass:sigvert");
137 fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptpi2bkg:PtKbkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:InvMassbkg:sigvertbkg");
fc8d975b 138
139 //fOutput->Add(fNtupleDplus); // AD
140 //fOutput->Add(fNtupleDplusbackg); //AD
141
d486095a 142 return;
143}
144
145//________________________________________________________________________
146void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
147{
148 // Execute analysis for current event:
149 // heavy flavor candidates association to MC truth
150
151
152
153 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
154
155 // load Dplus->Kpipi candidates
156 TClonesArray *array3Prong =
157 (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
158 if(!array3Prong) {
159 printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
160 return;
161 }
162
163 // AOD primary vertex
164 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
165 // vtx1->Print();
166
167 // load MC particles
168 TClonesArray *arrayMC =
169 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
170 if(!arrayMC) {
171 printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
172 return;
173 }
174
175 // load MC header
176 AliAODMCHeader *mcHeader =
177 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
178 if(!mcHeader) {
179 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
180 return;
181 }
182 Int_t n3Prong = array3Prong->GetEntriesFast();
183 printf("Number of D+->Kpipi: %d\n",n3Prong);
184
185
186 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
187 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
188
189
190 Bool_t unsetvtx=kFALSE;
191 if(!d->GetOwnPrimaryVtx()){
192 d->SetOwnPrimaryVtx(vtx1);
193 unsetvtx=kTRUE;
194 }
fc8d975b 195 if(d->SelectDplus(fVHF->GetDplusCuts()))
196 {
d486095a 197
198
fc8d975b 199 Int_t labDp = d->MatchToMC(411,arrayMC);
200 // if(labDp>=0) {
201 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
202 if(labDp>=0) {
203 Int_t pdgDp = TMath::Abs(partDp->GetPdgCode());
204 if(pdgDp==411){
205
206 fHistSignal->Fill(d->InvMassDplus());
207 fHistMass->Fill(d->InvMassDplus());
413e9336 208
fc8d975b 209 // Post the data already here
210 PostData(1,fOutput);
211
212 fNtupleDplus->Fill(pdgDp,partDp->Px()-d->Px(),partDp->Py()-d->Py(),partDp->Pz()-d->Pz(),d->PtProng(0),d->PtProng(2),d->PtProng(1),d->Pt(),partDp->Pt(),d->CosPointingAngle(),d->DecayLength(),partDp->Xv(),d->Xv(),d->InvMassDplus(),d->GetSigmaVert());
213 PostData(2,fNtupleDplus); //AD
214 }
215 }
216else {
217 fHistBackground->Fill(d->InvMassDplus());
218 fNtupleDplusbackg->Fill(d->PtProng(0),d->PtProng(2),d->PtProng(1),d->Pt(),d->CosPointingAngle(),d->DecayLength(),d->Xv(),d->InvMassDplus(),d->GetSigmaVert());
219 PostData(3,fNtupleDplusbackg); //AD
220
413e9336 221 fHistMass->Fill(d->InvMassDplus());
fc8d975b 222 }
223
413e9336 224
fc8d975b 225 }
d486095a 226
fc8d975b 227
228
d486095a 229 if(unsetvtx) d->UnsetOwnPrimaryVtx();
230
231 }
232 // end loop on D+->Kpipi
233
234
235 return;
236}
237
238//________________________________________________________________________
239void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
240{
241 // Terminate analysis
242 //
243 if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
244
245 fOutput = dynamic_cast<TList*> (GetOutputData(1));
246 if (!fOutput) {
247 printf("ERROR: fOutput not available\n");
248 return;
249 }
250
fc8d975b 251 //fNtupleDplus = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleDplus"));//AD
d486095a 252 fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
253 fHistSignal = dynamic_cast<TH1F*>(fOutput->FindObject("fHistSignal"));
254 fHistBackground = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBackground"));
fc8d975b 255 //fNtupleDplusbackg = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleDplusbackg")); //AD
10bdd1ae 256
fc8d975b 257 fNtupleDplus = dynamic_cast<TNtuple*> (GetOutputData(2)); //AD
258 fNtupleDplusbackg = dynamic_cast<TNtuple*> (GetOutputData(3)); //AD
259
260
261 return;
d486095a 262}
263