]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSESignificance.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New task for cut optimisation (Giacomo, Renu, Francesco, Chiara)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSESignificance.cxx
CommitLineData
cbedddce 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// AliAnalysisTaskSESignificane to calculate effects on
19// significance of D mesons cut
20// Authors: G. Ortona, ortona@to.infn.it
21// F. Prino, prino@to.infn.it
22// Renu Bala, bala@to.infn.it
23// Chiara Bianchin, cbianchi@pd.infn.it
24/////////////////////////////////////////////////////////////
25
26#include <Riostream.h>
27#include <TClonesArray.h>
28#include <TCanvas.h>
29#include <TList.h>
30#include <TFile.h>
31#include <TString.h>
32#include <TH1F.h>
33#include <TDatabasePDG.h>
34
35#include "AliAnalysisManager.h"
36#include "AliAODHandler.h"
37#include "AliAODEvent.h"
38#include "AliAODVertex.h"
39#include "AliAODTrack.h"
40#include "AliAODMCHeader.h"
41#include "AliAODMCParticle.h"
42#include "AliAODRecoDecayHF3Prong.h"
43#include "AliAODRecoDecayHF.h"
44#include "AliAODRecoDecayHF2Prong.h"
45#include "AliAODRecoDecayHF4Prong.h"
46#include "AliAODRecoCascadeHF.h"
47
48#include "AliAnalysisTaskSE.h"
49#include "AliRDHFCutsDplustoKpipi.h"
50#include "AliRDHFCutsD0toKpi.h"
51#include "AliMultiDimVector.h"
52
53#include "AliAnalysisTaskSESignificance.h"
54
55ClassImp(AliAnalysisTaskSESignificance)
56
57
58//________________________________________________________________________
59AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance():
60 AliAnalysisTaskSE(),
61 fOutput(0),
62 fCutList(0),
63 fHistNEvents(0),
64 fUpmasslimit(1.965),
65 fLowmasslimit(1.765),
66 fRDCuts(0),
67 fNPtBins(0),
68 fReadMC(kFALSE),
69 fDecChannel(0),
70 fSelectionlevel(0)
71{
72 // Default constructor
73}
74
75//________________________________________________________________________
76AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
77 AliAnalysisTaskSE(name),
78 fOutput(0),
79 fCutList(listMDV),
80 fHistNEvents(0),
81 fUpmasslimit(0),
82 fLowmasslimit(0),
83 fRDCuts(rdCuts),
84 fNPtBins(0),
85 fReadMC(kFALSE),
86 fDecChannel(decaychannel),
87 fSelectionlevel(selectionlevel)
88{
89
90 Int_t pdg=421;
91 switch(fDecChannel){
92 case 0:
93 pdg=411;
94 break;
95 case 1:
96 pdg=421;
97 break;
98 case 2:
99 pdg=413;
100 break;
101 case 3:
102 pdg=431;
103 break;
104 case 4:
105 pdg=421;
106 break;
107 case 5:
108 pdg=4122;
109 break;
110 }
111
112 SetMassLimits(0.1,pdg); //check range
113 fNPtBins=fRDCuts->GetNPtBins();
114
115 if(fDebug>1)fRDCuts->PrintAll();
116 // Output slot #1 writes into a TList container
117 DefineOutput(1,TList::Class()); //My private output
118 DefineOutput(2,TList::Class());
119}
120
121 //________________________________________________________________________
122AliAnalysisTaskSESignificance::~AliAnalysisTaskSESignificance()
123{
124 // Destructor
125 if (fOutput) {
126 delete fOutput;
127 fOutput = 0;
128 }
129 if (fCutList) {
130 delete fCutList;
131 fCutList = 0;
132 }
133 if(fHistNEvents){
134 delete fHistNEvents;
135 fHistNEvents=0;
136 }
137
138/*
139 if(fRDCuts) {
140 delete fRDCuts;
141 fRDCuts= 0;
142 }
143*/
144
145}
146//_________________________________________________________________
147void AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
148 Float_t mass=0;
149 Int_t abspdg=TMath::Abs(pdg);
150 if(abspdg==411||abspdg==421||abspdg==431)mass=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
151 fUpmasslimit = mass+range;
152 fLowmasslimit = mass-range;
153}
154//_________________________________________________________________
155void AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
156 if(uplimit>lowlimit)
157 {
158 fUpmasslimit = lowlimit;
159 fLowmasslimit = uplimit;
160 }
161}
162
163
164
165//________________________________________________________________________
166void AliAnalysisTaskSESignificance::LocalInit()
167{
168 // Initialization
169
170 if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
171
172 TList *mdvList = new TList();
173 mdvList->SetOwner();
174 mdvList = fCutList;
175
176 PostData(2,mdvList);
177
178 return;
179}
180//________________________________________________________________________
181void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
182{
183 // Create the output container
184
185 if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
186
187 // Several histograms are more conveniently managed in a TList
188 fOutput = new TList();
189 fOutput->SetOwner();
190 fOutput->SetName("OutputHistos");
191
192 //same number of steps in each multiDimVectorPtBin%d !
193 Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
194 nHist=nHist*fNPtBins;
195 for(Int_t i=0;i<nHist;i++){
196
197 TString hisname;
198 TString signame;
199 TString bkgname;
200 TString rflname;
201 TString title;
202
203 hisname.Form("hMass_%d",i);
204 signame.Form("hSig_%d",i);
205 bkgname.Form("hBkg_%d",i);
206 rflname.Form("hRfl_%d",i);
207
208 title.Form("Invariant mass;M[GeV/c^{2}];Entries");
209
210 fMassHist[i]=new TH1F(hisname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
211 fMassHist[i]->Sumw2();
212 fOutput->Add(fMassHist[i]);
213
214 if(fReadMC){
215 fSigHist[i]=new TH1F(signame.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
216 fSigHist[i]->Sumw2();
217 fOutput->Add(fSigHist[i]);
218
219 fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
220 fBkgHist[i]->Sumw2();
221 fOutput->Add(fBkgHist[i]);
222
223 fRflHist[i]=new TH1F(rflname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
224 fRflHist[i]->Sumw2();
225 fOutput->Add(fRflHist[i]);
226
227 }
228 }
229
230 fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",3,0,3.);
231 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
232 fHistNEvents->GetXaxis()->SetBinLabel(2,"nCandidatesSelected");
233 fHistNEvents->GetXaxis()->SetBinLabel(3,"nTotEntries Mass hists");
234 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
235 fOutput->Add(fHistNEvents);
236
237
238 return;
239}
240
241//________________________________________________________________________
242void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
243{
244 // Execute analysis for current event:
245 // heavy flavor candidates association to MC truth
246
247 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
248 if(fDebug>1) printf("Analysing decay %d\n",fDecChannel);
249 // Post the data already here
250 PostData(1,fOutput);
251 TClonesArray *arrayProng =0;
252
253 if(!aod && AODEvent() && IsStandardAOD()) {
254 // In case there is an AOD handler writing a standard AOD, use the AOD
255 // event in memory rather than the input (ESD) event.
256 aod = dynamic_cast<AliAODEvent*> (AODEvent());
257 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
258 // have to taken from the AOD event hold by the AliAODExtension
259 AliAODHandler* aodHandler = (AliAODHandler*)
260 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
261 if(aodHandler->GetExtensions()) {
262
263 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
264 AliAODEvent *aodFromExt = ext->GetAOD();
265
266
267 switch(fDecChannel){
268 case 0:
269 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
270 break;
271 case 1:
272 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
273 break;
274 case 2:
275 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
276 break;
277 case 3:
278 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
279 break;
280 case 4:
281 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
282 break;
283 case 5:
284 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
285 break;
286 }
287 }
288 } else {
289 switch(fDecChannel){
290 case 0:
291 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
292 break;
293 case 1:
294 arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
295 break;
296 case 2:
297 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
298 break;
299 case 3:
300 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
301 break;
302 case 4:
303 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
304 break;
305 case 5:
306 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
307 break;
308 }
309 }
310 if(!arrayProng) {
311 printf("AliAnalysisTaskSESignificance::UserExec: branch not found!\n");
312 return;
313 }
314
315 TClonesArray *arrayMC=0;
316 AliAODMCHeader *mcHeader=0;
317
318 // load MC particles
319 if(fReadMC){
320
321 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
322 if(!arrayMC) {
323 printf("AliAnalysisTaskSESignificance::UserExec: MC particles branch not found!\n");
324 // return;
325 }
326
327 // load MC header
328 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
329 if(!mcHeader) {
330 printf("AliAnalysisTaskSESignificance::UserExec: MC header branch not found!\n");
331 return;
332 }
333 }
334
335 fHistNEvents->Fill(0); // count event
336
337 AliAODRecoDecayHF *d=0;
338 Int_t nprongs=1;
339 Int_t *pdgdaughters=0x0;
340 Int_t prongpdg=211;
341
342
343 switch(fDecChannel){
344 case 0:
345 //D+
346 pdgdaughters =new Int_t[3];
347 pdgdaughters[0]=211;//pi
348 pdgdaughters[1]=321;//K
349 pdgdaughters[2]=211;//pi
350 nprongs=3;
351 prongpdg=411;
352 break;
353 case 1:
354 //D0
355 pdgdaughters =new Int_t[2];
356 pdgdaughters[0]=211;//pi
357 pdgdaughters[1]=321;//K
358 nprongs=2;
359 prongpdg=421;
360 break;
361 case 2:
362 //D*
363 pdgdaughters =new Int_t[3];
364 pdgdaughters[1]=211;//pi
365 pdgdaughters[0]=321;//K
366 pdgdaughters[2]=211;//pi (soft?)
367 nprongs=3;
368 prongpdg=413;
369 break;
370 case 3:
371 //Ds
372 pdgdaughters =new Int_t[3];
373 pdgdaughters[0]=211;//pi
374 pdgdaughters[1]=321;//K
375 pdgdaughters[2]=321;//K
376 nprongs=3;
377 prongpdg=431;
378 break;
379 case 4:
380 //D0
381 pdgdaughters =new Int_t[4];
382 pdgdaughters[0]=321;
383 pdgdaughters[1]=211;
384 pdgdaughters[2]=211;
385 pdgdaughters[3]=211;
386 nprongs=4;
387 prongpdg=421;
388 break;
389 case 5:
390 //Lambda_c
391 pdgdaughters =new Int_t[3];
392 //check order
393 pdgdaughters[0]=2212;//p
394 pdgdaughters[1]=321;//K
395 pdgdaughters[2]=211;//pi
396 nprongs=3;
397 prongpdg=4122;
398 break;
399 }
400
401 Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
402 Int_t nProng = arrayProng->GetEntriesFast();
403 if(fDebug>1) printf("Number of D2H: %d\n",nProng);
404
405 for (Int_t iProng = 0; iProng < nProng; iProng++) {
406 d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
407 Double_t invMass=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
408 Double_t invMassC=0;
409
410 Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel);
411
412 if(isSelected) {
413 fHistNEvents->Fill(1); // count selected with loosest cuts
414 if(fDebug>1) printf("Is Selected\n");
415
416 const Int_t nvars=fRDCuts->GetNVarsForOpt();
417 Float_t *vars = new Float_t[nvars];
418 Int_t nVals=0;
419
420 fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
421 Int_t ptbin=fRDCuts->PtBin(d->Pt());
422 if(ptbin==-1) continue;
423 TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
424 ULong64_t *addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
425 if(fDebug>1)printf("nvals = %d\n",nVals);
426
427 for(Int_t ivals=0;ivals<nVals;ivals++){
428 if(addresses[ivals]>=((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetNTotCells()){
429 if (fDebug>1) printf("Overflow!!\n");
430 return;
431 }
432 fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
433 fHistNEvents->Fill(2);
434 if(fReadMC){
435 Int_t lab=-1;
436 lab = d->MatchToMC(prongpdg,arrayMC,nprongs,pdgdaughters);
437 if(lab>=0){ //signal
438 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
439 Int_t pdgMC = dMC->GetPdgCode();
440
441 if(pdgMC==+prongpdg) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
442 else fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
443 }
444 else{ //background
445 fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
446 }
447
448 }
449 }
450
451
452 if(fDecChannel==AliAnalysisTaskSESignificance::kD0toKpi){ //repeat the cycle to checks cuts passed as D0bar
453 pdgdaughters[0]=321;
454 pdgdaughters[1]=211;
455 invMass=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
456 fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
457 addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
458 for(Int_t ivals=0;ivals<nVals;ivals++){
459 fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
460 if(fReadMC){
461 Int_t lab=-1;
462 lab = d->MatchToMC(-prongpdg,arrayMC,nprongs,pdgdaughters);
463 if(lab>=0){ //signal
464 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
465 Int_t pdgMC = dMC->GetPdgCode();
466 if(pdgMC==-prongpdg) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
467 else fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
468 }
469 else{ //background
470 fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
471 }
472
473 }
474 }
475 }
476
477 }// end if selected
478
479
480 }
481
482 PostData(1,fOutput);
483 return;
484}
485
486
487//________________________________________________________________________
488void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
489{
490 // Terminate analysis
491 //
492 if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
493
494
495 fOutput = dynamic_cast<TList*> (GetOutputData(1));
496 if (!fOutput) {
497 printf("ERROR: fOutput not available\n");
498 return;
499 }
500
501 fCutList = dynamic_cast<TList*> (GetOutputData(2));
502 if (!fCutList) {
503 printf("ERROR: fCutList not available\n");
504 return;
505 }
506
507 AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
508 if (!mdvtmp){
509 cout<<"multidimvec not found in TList"<<endl;
510 fCutList->ls();
511 return;
512 }
513 Int_t nHist=mdvtmp->GetNTotCells();
514 TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
515 Bool_t drawn=kFALSE;
516 for(Int_t i=0;i<nHist;i++){
517
518 TString hisname;
519 TString signame;
520 TString bkgname;
521 TString rflname;
522
523 hisname.Form("hMass_%d",i);
524 signame.Form("hSig_%d",i);
525 bkgname.Form("hBkg_%d",i);
526 rflname.Form("hRfl_%d",i);
527
528 fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
529
530 if (!drawn && fMassHist[i]->GetEntries() > 0){
531 c1->cd();
532 fMassHist[i]->Draw();
533 drawn=kTRUE;
534 }
535
536 if(fReadMC){
537 fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
538 fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
539 fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
540 }
541
542 }
543
544
545
546
547 return;
548}
549