New file to configure TOF preprocessor for NOISE runs:
[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>
b557eb43 28
29#include "AliAnalysisManager.h"
30#include "AliAODHandler.h"
d486095a 31#include "AliAODEvent.h"
32#include "AliAODVertex.h"
33#include "AliAODTrack.h"
34#include "AliAODMCHeader.h"
35#include "AliAODMCParticle.h"
36#include "AliAODRecoDecayHF3Prong.h"
37#include "AliAnalysisVertexingHF.h"
38#include "AliAnalysisTaskSE.h"
39#include "AliAnalysisTaskSEDplus.h"
40
41ClassImp(AliAnalysisTaskSEDplus)
42
43
44//________________________________________________________________________
45AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
46AliAnalysisTaskSE(),
47fOutput(0),
48fNtupleDplus(0),
49fNtupleDplusbackg(0),
1f4e9722 50fFillNtuple(kFALSE),
d486095a 51fVHF(0)
52{
53 // Default constructor
d486095a 54}
55
56//________________________________________________________________________
1f4e9722 57AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
d486095a 58AliAnalysisTaskSE(name),
59fOutput(0),
60fNtupleDplus(0),
61fNtupleDplusbackg(0),
1f4e9722 62fFillNtuple(fillNtuple),
d486095a 63fVHF(0)
64{
65 // Default constructor
66
67 // Output slot #1 writes into a TList container
68 DefineOutput(1,TList::Class()); //My private output
1f4e9722 69
70 if(fFillNtuple){
71 // Output slot #2 writes into a TNtuple container
72 DefineOutput(2,TNtuple::Class()); //My private output
73 // Output slot #3 writes into a TNtuple container
74 DefineOutput(3,TNtuple::Class()); //My private output
75 }
d486095a 76}
77
78//________________________________________________________________________
79AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
80{
81 // Destructor
82 if (fOutput) {
83 delete fOutput;
84 fOutput = 0;
85 }
86 if (fVHF) {
87 delete fVHF;
88 fVHF = 0;
89 }
90}
91
92//________________________________________________________________________
93void AliAnalysisTaskSEDplus::Init()
94{
95 // Initialization
96
97 if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
98
99 gROOT->LoadMacro("ConfigVertexingHF.C");
100
101 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
102 fVHF->PrintStatus();
103
104 return;
105}
106
107//________________________________________________________________________
108void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
109{
110 // Create the output container
111 //
112 if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
113
114 // Several histograms are more conveniently managed in a TList
115 fOutput = new TList();
116 fOutput->SetOwner();
1f4e9722 117 fOutput->SetName("OutputHistos");
118
119 Int_t nPtBins=4;
120 TString hisname;
121 for(Int_t i=0;i<nPtBins;i++){
122 hisname.Form("hMassPt%d",i);
123 TH1F* hm=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
124 hm->Sumw2();
125 hisname.Form("hSigPt%d",i);
126 TH1F* hs=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
127 hs->Sumw2();
128 hisname.Form("hBkgPt%d",i);
129 TH1F* hb=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
130 hb->Sumw2();
131
132 hisname.Form("hMassPt%dTC",i);
133 TH1F* hmtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
134 hmtc->Sumw2();
135 hisname.Form("hSigPt%dTC",i);
136 TH1F* hstc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
137 hstc->Sumw2();
138 hisname.Form("hBkgPt%dTC",i);
139 TH1F* hbtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
140 hbtc->Sumw2();
141
142
143 fOutput->Add(hm);
144 fOutput->Add(hs);
145 fOutput->Add(hb);
146 fOutput->Add(hmtc);
147 fOutput->Add(hstc);
148 fOutput->Add(hbtc);
149 }
d486095a 150
8c34bd86 151
1f4e9722 152 if(fFillNtuple){
153 OpenFile(2); // 2 is the slot number of the ntuple
154 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");
155 OpenFile(3); // 3 is the slot number of the ntuple
156 fNtupleDplusbackg = new TNtuple("fNtupleDplusbackg","D + backg","Ptpibkg:Ptkbkg:Ptpi2bkg:PtRecbkg:PointingAnglebkg:DLbkg:VxRecbkg:VyRecbkg:VzRecbkg:InvMassbkg:sigvertbkg:d0Pibkg:d0Kbkg:d0Pi2bkg");
157 }
158
d486095a 159 return;
160}
161
162//________________________________________________________________________
163void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
164{
165 // Execute analysis for current event:
166 // heavy flavor candidates association to MC truth
167
168
169
170 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
171
b557eb43 172 TClonesArray *array3Prong = 0;
173
174 if(!aod && AODEvent() && IsStandardAOD()) {
175 // In case there is an AOD handler writing a standard AOD, use the AOD
176 // event in memory rather than the input (ESD) event.
177 aod = dynamic_cast<AliAODEvent*> (AODEvent());
178 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
179 // have to taken from the AOD event hold by the AliAODExtension
180 AliAODHandler* aodHandler = (AliAODHandler*)
181 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
182 if(aodHandler->GetExtensions()) {
183 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
184 AliAODEvent *aodFromExt = ext->GetAOD();
185 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
186 }
187 } else {
188 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
189 }
8931c313 190
d486095a 191 if(!array3Prong) {
192 printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
193 return;
194 }
195
196 // AOD primary vertex
1f4e9722 197 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
198 // vtx1->Print();
199
d486095a 200 // load MC particles
201 TClonesArray *arrayMC =
202 (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
203 if(!arrayMC) {
204 printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
205 return;
206 }
1f4e9722 207
d486095a 208 // load MC header
209 AliAODMCHeader *mcHeader =
210 (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
211 if(!mcHeader) {
212 printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
213 return;
214 }
1f4e9722 215 Int_t n3Prong = array3Prong->GetEntriesFast();
216 if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
d486095a 217
218
1f4e9722 219 Int_t pdgDgDplustoKpipi[3]={321,211,211};
220 Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
221 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
222 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
223
224
225 Bool_t unsetvtx=kFALSE;
226 if(!d->GetOwnPrimaryVtx()){
227 d->SetOwnPrimaryVtx(vtx1);
228 unsetvtx=kTRUE;
229 }
230 if(d->SelectDplus(fVHF->GetDplusCuts())) {
231 Int_t iPtBin=0;
232 Double_t ptCand = d->Pt();
233 if(ptCand<2.){
234 iPtBin=0;
235 cutsDplus[7]=0.08;
236 cutsDplus[8]=0.5;
237 cutsDplus[10]=0.979;
238 }
239 else if(ptCand>2. && ptCand<3){
240 iPtBin=1;
241 cutsDplus[7]=0.08;
242 cutsDplus[8]=0.5;
243 cutsDplus[9]=0.991;
244 }else if(ptCand>3. && ptCand<5){
245 iPtBin=2;
246 cutsDplus[7]=0.1;
247 cutsDplus[8]=0.5;
248 cutsDplus[9]=0.9955;
249 }else{
250 iPtBin=3;
251 cutsDplus[7]=0.1;
252 cutsDplus[8]=0.5;
253 cutsDplus[9]=0.997;
254 }
255 Bool_t passTightCuts=d->SelectDplus(cutsDplus);
256 Int_t labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
257 Double_t invMass=d->InvMassDplus();
258
259 TString hisNameA(Form("hMassPt%d",iPtBin));
260 TString hisNameS(Form("hSigPt%d",iPtBin));
261 TString hisNameB(Form("hBkgPt%d",iPtBin));
262 TString hisNameATC(Form("hMassPt%dTC",iPtBin));
263 TString hisNameSTC(Form("hSigPt%dTC",iPtBin));
264 TString hisNameBTC(Form("hBkgPt%dTC",iPtBin));
265
266 ((TH1F*)(fOutput->FindObject(hisNameA)))->Fill(invMass);
267 if(passTightCuts){
268 ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
269 }
d486095a 270
271
1f4e9722 272 if(labDp>=0) {
273 ((TH1F*)(fOutput->FindObject(hisNameS)))->Fill(invMass);
274 if(passTightCuts){
275 ((TH1F*)(fOutput->FindObject(hisNameSTC)))->Fill(invMass);
276 }
277 if(fFillNtuple){
8c34bd86 278 AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
1f4e9722 279 AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
280 Float_t tmp[20];
281 tmp[0]=TMath::Abs(partDp->GetPdgCode());
282 tmp[1]=partDp->Px()-d->Px();
283 tmp[2]=partDp->Py()-d->Py();
284 tmp[3]=partDp->Pz()-d->Pz();
285 tmp[4]=d->PtProng(0);
286 tmp[5]=d->PtProng(1);
287 tmp[6]=d->PtProng(2);
288 tmp[7]=d->Pt();
289 tmp[8]=partDp->Pt();
290 tmp[9]=d->CosPointingAngle();
291 tmp[10]=d->DecayLength();
292 tmp[11]=dg0->Xv();
293 tmp[12]=d->Xv();
294 tmp[13]=d->Yv();
295 tmp[14]=d->Zv();
296 tmp[15]=d->InvMassDplus();
297 tmp[16]=d->GetSigmaVert();
298 tmp[17]=d->Getd0Prong(0);
299 tmp[18]=d->Getd0Prong(1);
300 tmp[19]=d->Getd0Prong(2);
301 fNtupleDplus->Fill(tmp);
302 PostData(2,fNtupleDplus);
303 }
304 }else{
305 ((TH1F*)(fOutput->FindObject(hisNameB)))->Fill(invMass);
306 if(passTightCuts){
307 ((TH1F*)(fOutput->FindObject(hisNameBTC)))->Fill(invMass);
308 }
309 if(fFillNtuple){
310 Float_t tmpbkg[14];
311 tmpbkg[0]=d->PtProng(0);
312 tmpbkg[1]=d->PtProng(1);
313 tmpbkg[2]=d->PtProng(2);
314 tmpbkg[3]=d->Pt();
315 tmpbkg[4]=d->CosPointingAngle();
316 tmpbkg[5]=d->DecayLength();
317 tmpbkg[6]=d->Xv();
318 tmpbkg[7]=d->Yv();
319 tmpbkg[8]=d->Zv();
320 tmpbkg[9]=d->InvMassDplus();
321 tmpbkg[10]=d->GetSigmaVert();
322 tmpbkg[11]=d->Getd0Prong(0);
323 tmpbkg[12]=d->Getd0Prong(1);
324 tmpbkg[13]=d->Getd0Prong(2);
8c34bd86 325 fNtupleDplusbackg->Fill(tmpbkg);
326 PostData(3,fNtupleDplusbackg);
fc8d975b 327 }
8c34bd86 328 }
1f4e9722 329 PostData(1,fOutput);
d486095a 330 }
1f4e9722 331 if(unsetvtx) d->UnsetOwnPrimaryVtx();
332 }
333
d486095a 334 return;
335}
336
337//________________________________________________________________________
338void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
339{
340 // Terminate analysis
341 //
342 if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
343
344 fOutput = dynamic_cast<TList*> (GetOutputData(1));
345 if (!fOutput) {
346 printf("ERROR: fOutput not available\n");
347 return;
348 }
349
1f4e9722 350 if(fFillNtuple){
351 fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
352 fNtupleDplusbackg = dynamic_cast<TNtuple*>(GetOutputData(3));
353 }
fc8d975b 354
8c34bd86 355 return;
d486095a 356}
357