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