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