Updated cuts classes (Rossella, Renu, Francesco, Giacomo, Chiara)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDs.cxx
CommitLineData
25c94fc8 1/**************************************************************************
2 * Copyright(c) 2007-2009, 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/* $Id: AliITSCorrMapSDD.cxx 32906 2009-06-12 16:56:53Z prino $ */
17
18///////////////////////////////////////////////////////////////////
19// //
20// Analysis task to produce Ds candidates mass spectra //
21// Origin: F.Prino, Torino, prino@to.infn.it //
22// //
23///////////////////////////////////////////////////////////////////
24
25#include <TClonesArray.h>
26#include <TNtuple.h>
27#include <TList.h>
28#include <TString.h>
29#include <TH1F.h>
30#include <TMath.h>
31#include <TDatabasePDG.h>
32
33#include "AliAnalysisManager.h"
34#include "AliAODHandler.h"
35#include "AliAODEvent.h"
36#include "AliAODVertex.h"
37#include "AliAODTrack.h"
38#include "AliAODMCHeader.h"
39#include "AliAODMCParticle.h"
40#include "AliAODRecoDecayHF3Prong.h"
41#include "AliAnalysisVertexingHF.h"
42#include "AliAnalysisTaskSE.h"
43#include "AliAnalysisTaskSEDs.h"
44
45ClassImp(AliAnalysisTaskSEDs)
46
47
48//________________________________________________________________________
49AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
50 AliAnalysisTaskSE(),
51 fOutput(0),
52 fHistNEvents(0),
53 fReadMC(kFALSE),
54 fNPtBins(0),
55 fMassRange(0.2),
56 fVHF(0)
57{
58 // Default constructor
59}
60
61//________________________________________________________________________
62AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name):
63 AliAnalysisTaskSE(name),
64 fOutput(0),
65 fHistNEvents(0),
66 fReadMC(kFALSE),
67 fNPtBins(0),
68 fMassRange(0.2),
69 fVHF(0)
70{
71 // Default constructor
72 // Output slot #1 writes into a TList container
73 Double_t ptlim[6]={0.,1.,3.,5.,10.,9999999.};
74 SetPtBins(5,ptlim);
75
76 DefineOutput(1,TList::Class()); //My private output
77
78}
79
80//________________________________________________________________________
81void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Double_t* lim){
82 // define pt bins for analysis
83 if(n>kMaxPtBins){
84 printf("Max. number of Pt bins = %d\n",kMaxPtBins);
85 fNPtBins=kMaxPtBins;
86 fPtLimits[0]=0.;
87 fPtLimits[1]=1.;
88 fPtLimits[2]=3.;
89 fPtLimits[3]=5.;
90 fPtLimits[4]=10.;
91 for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
92 }else{
93 fNPtBins=n;
94 for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
95 for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
96 }
97 if(fDebug > 1){
98 printf("Number of Pt bins = %d\n",fNPtBins);
99 for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
100 }
101}
102//________________________________________________________________________
103AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
104{
105 // Destructor
106 if (fOutput) {
107 delete fOutput;
108 fOutput = 0;
109 }
110 if (fVHF) {
111 delete fVHF;
112 fVHF = 0;
113 }
114}
115
116//________________________________________________________________________
117void AliAnalysisTaskSEDs::Init()
118{
119 // Initialization
120
121 if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
122
123 gROOT->LoadMacro("ConfigVertexingHF.C");
124
125 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
126 fVHF->PrintStatus();
127
128 return;
129}
130
131//________________________________________________________________________
132void AliAnalysisTaskSEDs::UserCreateOutputObjects()
133{
134 // Create the output container
135 //
136 if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
137
138 // Several histograms are more conveniently managed in a TList
139 fOutput = new TList();
140 fOutput->SetOwner();
141 fOutput->SetName("OutputHistos");
142
143 Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
144 Double_t minMass=massDs-0.5*fMassRange;
145 Double_t maxMass=massDs+0.5*fMassRange;
146 TString hisname;
147 Int_t index;
148 for(Int_t i=0;i<fNPtBins;i++){
149 index=GetHistoIndex(i);
150 hisname.Form("hMassAllPt%d",i);
151 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
152 fMassHist[index]->Sumw2();
153 hisname.Form("hMassAllPt%dCuts",i);
154 fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
155 fMassHistCuts[index]->Sumw2();
156 hisname.Form("hCosPAllPt%d",i);
157 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
158 fCosPHist[index]->Sumw2();
159 hisname.Form("hDLenAllPt%d",i);
160 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
161 fDLenHist[index]->Sumw2();
162 hisname.Form("hDalitzAllPt%d",i);
163 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
164 fDalitz[index]->Sumw2();
165
166 index=GetSignalHistoIndex(i);
167 hisname.Form("hMassSigPt%d",i);
168 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
169 fMassHist[index]->Sumw2();
170 hisname.Form("hMassSigPt%dCuts",i);
171 fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
172 fMassHistCuts[index]->Sumw2();
173 hisname.Form("hCosPSigPt%d",i);
174 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
175 fCosPHist[index]->Sumw2();
176 hisname.Form("hDLenSigPt%d",i);
177 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
178 fDLenHist[index]->Sumw2();
179 hisname.Form("hDalitzSigPt%d",i);
180 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
181 fDalitz[index]->Sumw2();
182
183 hisname.Form("hMassBkgPt%d",i);
184 index=GetBackgroundHistoIndex(i);
185 fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
186 fMassHist[index]->Sumw2();
187 hisname.Form("hMassBkgPt%dCuts",i);
188 fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
189 fMassHistCuts[index]->Sumw2();
190 hisname.Form("hCosPBkgPt%d",i);
191 fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
192 fCosPHist[index]->Sumw2();
193 hisname.Form("hDLenBkgPt%d",i);
194 fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
195 fDLenHist[index]->Sumw2();
196 hisname.Form("hDalitzBkgPt%d",i);
197 fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
198 fDalitz[index]->Sumw2();
199 }
200
201 for(Int_t i=0; i<3*fNPtBins; i++){
202 fOutput->Add(fMassHist[i]);
203 fOutput->Add(fMassHistCuts[i]);
204 fOutput->Add(fCosPHist[i]);
205 fOutput->Add(fDLenHist[i]);
206 fOutput->Add(fDalitz[i]);
207 }
208
209 fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",4,-0.5,3.5);
210 fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",4,-0.5,3.5);
211 fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",4,-0.5,3.5);
212 fChanHistCuts[0] = new TH1F("hChanAllCuts", "KKpi and piKK candidates",4,-0.5,3.5);
213 fChanHistCuts[1] = new TH1F("hChanSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
214 fChanHistCuts[2] = new TH1F("hChanBkgCuts", "KKpi and piKK candidates",4,-0.5,3.5);
215 for(Int_t i=0;i<3;i++){
216 fChanHist[i]->Sumw2();
217 fChanHist[i]->SetMinimum(0);
218 fChanHistCuts[i]->Sumw2();
219 fChanHistCuts[i]->SetMinimum(0);
220 fOutput->Add(fChanHist[i]);
221 fOutput->Add(fChanHistCuts[i]);
222 }
223
224 fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
225 fHistNEvents->Sumw2();
226 fHistNEvents->SetMinimum(0);
227 fOutput->Add(fHistNEvents);
228
229
230 return;
231}
232
233//________________________________________________________________________
234void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
235{
236 // Ds selection for current event, fill mass histos and selecetion variable histo
237 // separate signal and backgound if fReadMC is activated
238
239 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
240 fHistNEvents->Fill(0); // count event
241 // Post the data already here
242 PostData(1,fOutput);
243
244
245 TClonesArray *array3Prong = 0;
246 if(!aod && AODEvent() && IsStandardAOD()) {
247 // In case there is an AOD handler writing a standard AOD, use the AOD
248 // event in memory rather than the input (ESD) event.
249 aod = dynamic_cast<AliAODEvent*> (AODEvent());
250 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
251 // have to taken from the AOD event hold by the AliAODExtension
252 AliAODHandler* aodHandler = (AliAODHandler*)
253 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
254 if(aodHandler->GetExtensions()) {
255 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
256 AliAODEvent *aodFromExt = ext->GetAOD();
257 array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
258 }
259 } else {
260 array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
261 }
262
263 if(!array3Prong) {
264 printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
265 return;
266 }
267
268
269 TClonesArray *arrayMC=0;
270 AliAODMCHeader *mcHeader=0;
271
272 // AOD primary vertex
273 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
274 // vtx1->Print();
275
276 // load MC particles
277 if(fReadMC){
278
279 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
280 if(!arrayMC) {
281 printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
282 return;
283 }
284
285 // load MC header
286 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
287 if(!mcHeader) {
288 printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
289 return;
290 }
291 }
292
293 Int_t n3Prong = array3Prong->GetEntriesFast();
294 if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
295
296
297 Int_t pdgDstoKKpi[3]={321,321,211};
298 Double_t cutsDs[14]={0.2,0.4,0.4,0.,0.,0.005,0.038,0.,0.,0.95,0.,0.1,0.004,0.035};
299 for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
300 AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
301
302
303 Bool_t unsetvtx=kFALSE;
304 if(!d->GetOwnPrimaryVtx()){
305 d->SetOwnPrimaryVtx(vtx1);
306 unsetvtx=kTRUE;
307 }
308 Int_t isKKpi,ispiKK;
309 Int_t isPhi,isK0star;
310 Int_t isKKpiTC,ispiKKTC;
311 Int_t isPhiTC,isK0starTC;
312 if(d->SelectDs(fVHF->GetDsCuts(),isKKpi,ispiKK,isPhi,isK0star)) {
313 Double_t ptCand = d->Pt();
314 Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,ptCand);
315 Bool_t passTightCuts=d->SelectDs(cutsDs,isKKpiTC,ispiKKTC,isPhiTC,isK0starTC);
316 Int_t labDs=-1;
317 if(fReadMC){
318 labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
319 }
320
321 Double_t dlen=d->DecayLength();
322 Double_t cosp=d->CosPointingAngle();
323 Int_t index=GetHistoIndex(iPtBin);
324 Int_t type=0;
325 if(isKKpi) type+=1;
326 if(ispiKK) type+=2;
327 Int_t typeTC=0;
328 if(isKKpiTC) typeTC+=1;
329 if(ispiKKTC) typeTC+=2;
330 fCosPHist[index]->Fill(cosp);
331 fDLenHist[index]->Fill(dlen);
332 fChanHist[0]->Fill(type);
333 if(passTightCuts) fChanHistCuts[0]->Fill(typeTC);
334 if(fReadMC){
335 if(labDs>=0) {
336 index=GetSignalHistoIndex(iPtBin);
337 fCosPHist[index]->Fill(cosp);
338 fDLenHist[index]->Fill(dlen);
339 fChanHist[1]->Fill(type);
340 if(passTightCuts) fChanHistCuts[1]->Fill(typeTC);
341 }else{
342 index=GetBackgroundHistoIndex(iPtBin);
343 fCosPHist[index]->Fill(cosp);
344 fDLenHist[index]->Fill(dlen);
345 fChanHist[2]->Fill(type);
346 if(passTightCuts) fChanHistCuts[2]->Fill(typeTC);
347 }
348 }
349 if(isKKpi){
350 index=GetHistoIndex(iPtBin);
351 Double_t invMass=d->InvMassDsKKpi();
352 fMassHist[index]->Fill(invMass);
353 Double_t mass01=d->InvMass2Prongs(0,1,321,321);
354 Double_t mass12=d->InvMass2Prongs(1,2,321,211);
355 fDalitz[index]->Fill(mass01,mass12);
356 if(passTightCuts && isKKpiTC) fMassHistCuts[index]->Fill(invMass);
357 if(fReadMC){
358 if(labDs>=0) {
359 index=GetSignalHistoIndex(iPtBin);
360 fMassHist[index]->Fill(invMass);
361 fDalitz[index]->Fill(mass01,mass12);
362 if(passTightCuts&& isKKpiTC) fMassHistCuts[index]->Fill(invMass);
363 }else{
364 index=GetBackgroundHistoIndex(iPtBin);
365 fMassHist[index]->Fill(invMass);
366 fDalitz[index]->Fill(mass01,mass12);
367 if(passTightCuts&& isKKpiTC) fMassHistCuts[index]->Fill(invMass);
368 }
369 }
370 }
371 if(ispiKK){
372 index=GetHistoIndex(iPtBin);
373 Double_t invMass=d->InvMassDspiKK();
374 fMassHist[index]->Fill(invMass);
375 if(passTightCuts && ispiKKTC) fMassHistCuts[index]->Fill(invMass);
376 if(fReadMC){
377 if(labDs>=0) {
378 index=GetSignalHistoIndex(iPtBin);
379 fMassHist[index]->Fill(invMass);
380 if(passTightCuts && ispiKKTC) fMassHistCuts[index]->Fill(invMass);
381 }else{
382 index=GetBackgroundHistoIndex(iPtBin);
383 fMassHist[index]->Fill(invMass);
384 if(passTightCuts && ispiKKTC) fMassHistCuts[index]->Fill(invMass);
385 }
386 }
387 }
388 }
389 if(unsetvtx) d->UnsetOwnPrimaryVtx();
390 }
391
392
393 PostData(1,fOutput);
394 return;
395}
396
397//_________________________________________________________________
398
399void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
400{
401 // Terminate analysis
402 //
403 if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
404 fOutput = dynamic_cast<TList*> (GetOutputData(1));
405 if (!fOutput) {
406 printf("ERROR: fOutput not available\n");
407 return;
408 }
409 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
410 fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
411 fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
412 fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
413 fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
414 fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
415 fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));
416
417
418 TString hisname;
419 Int_t index;
420 for(Int_t i=0;i<fNPtBins;i++){
421
422 index=GetHistoIndex(i);
423 hisname.Form("hMassAllPt%d",i);
424 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
425 hisname.Form("hMassAllPt%dCuts",i);
426 fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
427 hisname.Form("hCosPAllPt%d",i);
428 fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
429 hisname.Form("hDLenAllPt%d",i);
430 fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
431 hisname.Form("hDalitzAllPt%d",i);
432 fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
433
434 index=GetSignalHistoIndex(i);
435 hisname.Form("hMassSigPt%d",i);
436 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
437 hisname.Form("hMassSigPt%dCuts",i);
438 fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
439 hisname.Form("hCosPSigPt%d",i);
440 fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
441 hisname.Form("hDLenSigPt%d",i);
442 fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
443 hisname.Form("hDalitzSigPt%d",i);
444 fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
445
446 index=GetBackgroundHistoIndex(i);
447 hisname.Form("hMassBkgPt%d",i);
448 fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
449 hisname.Form("hMassBkgPt%dCuts",i);
450 fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
451 hisname.Form("hCosPBkgPt%d",i);
452 fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
453 hisname.Form("hDLenBkgPt%d",i);
454 fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
455 hisname.Form("hDalitzBkgPt%d",i);
456 fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
457
458 }
459 return;
460}