]> 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
Added event-level selection (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());
3aa6ce53 119 CheckConsistency();
cbedddce 120}
121
122 //________________________________________________________________________
123AliAnalysisTaskSESignificance::~AliAnalysisTaskSESignificance()
124{
125 // Destructor
126 if (fOutput) {
127 delete fOutput;
128 fOutput = 0;
129 }
130 if (fCutList) {
131 delete fCutList;
132 fCutList = 0;
133 }
134 if(fHistNEvents){
135 delete fHistNEvents;
136 fHistNEvents=0;
137 }
138
139/*
140 if(fRDCuts) {
141 delete fRDCuts;
142 fRDCuts= 0;
143 }
144*/
145
3aa6ce53 146}
147//_________________________________________________________________
148Bool_t AliAnalysisTaskSESignificance::CheckConsistency(){
149
150 Bool_t result = kTRUE;
151
152 const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
153 //Float_t *vars = new Float_t[nvars];
154 Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
155 Bool_t *uppervars = fRDCuts->GetIsUpperCut();
156 TString *names = fRDCuts->GetVarNames();
157
158 for(Int_t i=0;i<fNPtBins;i++){
159 TString mdvname=Form("multiDimVectorPtBin%d",i);
160 Int_t ic=0;
161
162 for(Int_t ivar=0;ivar<nvars;ivar++){
163 if(varsforopt[ivar]){
164 Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
165 Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
166 if(min==max){
167 printf("AliAnalysisTaskSESignificance::CheckConsistency: ERROR! \n tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i);
168 result = kFALSE;
169 }
170 Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
171 if(uppervars[ivar]&&lowermdv){
172 printf("AliAnalysisTaskSESignificance::CheckConsistency: WARNING! %s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---Task will use swapped Tight/Loose cuts \n ",names[ivar].Data(),i);
173 ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->SwapLimits(ic);
174 result = kTRUE;
175 }
176 if(!uppervars[ivar]&&!lowermdv){
177 printf("AliAnalysisTaskSESignificance::CheckConsistency: WARNING! %s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---Task will use swapped Tight/Loose cuts \n",names[ivar].Data(),i);
178 ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->SwapLimits(ic);
179 result = kTRUE;
180 }
181 ic++;
182 }
183 }
184 }
185 return result;
186}
cbedddce 187//_________________________________________________________________
188void AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
189 Float_t mass=0;
190 Int_t abspdg=TMath::Abs(pdg);
5ecc6e9b 191 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
cbedddce 192 fUpmasslimit = mass+range;
193 fLowmasslimit = mass-range;
194}
195//_________________________________________________________________
196void AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
197 if(uplimit>lowlimit)
198 {
199 fUpmasslimit = lowlimit;
200 fLowmasslimit = uplimit;
201 }
202}
203
204
205
206//________________________________________________________________________
207void AliAnalysisTaskSESignificance::LocalInit()
208{
209 // Initialization
210
211 if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
212
213 TList *mdvList = new TList();
214 mdvList->SetOwner();
215 mdvList = fCutList;
216
217 PostData(2,mdvList);
218
219 return;
220}
221//________________________________________________________________________
222void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
223{
224 // Create the output container
225
226 if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
227
228 // Several histograms are more conveniently managed in a TList
229 fOutput = new TList();
230 fOutput->SetOwner();
231 fOutput->SetName("OutputHistos");
232
233 //same number of steps in each multiDimVectorPtBin%d !
234 Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
235 nHist=nHist*fNPtBins;
236 for(Int_t i=0;i<nHist;i++){
237
238 TString hisname;
239 TString signame;
240 TString bkgname;
241 TString rflname;
242 TString title;
243
244 hisname.Form("hMass_%d",i);
245 signame.Form("hSig_%d",i);
246 bkgname.Form("hBkg_%d",i);
247 rflname.Form("hRfl_%d",i);
248
249 title.Form("Invariant mass;M[GeV/c^{2}];Entries");
250
251 fMassHist[i]=new TH1F(hisname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
252 fMassHist[i]->Sumw2();
253 fOutput->Add(fMassHist[i]);
254
255 if(fReadMC){
256 fSigHist[i]=new TH1F(signame.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
257 fSigHist[i]->Sumw2();
258 fOutput->Add(fSigHist[i]);
259
260 fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
261 fBkgHist[i]->Sumw2();
262 fOutput->Add(fBkgHist[i]);
263
264 fRflHist[i]=new TH1F(rflname.Data(),title.Data(),100,fLowmasslimit,fUpmasslimit);
265 fRflHist[i]->Sumw2();
266 fOutput->Add(fRflHist[i]);
267
268 }
269 }
270
5ecc6e9b 271 fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",4,0,4.);
cbedddce 272 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
5ecc6e9b 273 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvNoSelected");
274 fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
275 fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
cbedddce 276 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
277 fOutput->Add(fHistNEvents);
278
279
280 return;
281}
282
283//________________________________________________________________________
284void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
285{
286 // Execute analysis for current event:
287 // heavy flavor candidates association to MC truth
288
289 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
290 if(fDebug>1) printf("Analysing decay %d\n",fDecChannel);
291 // Post the data already here
292 PostData(1,fOutput);
293 TClonesArray *arrayProng =0;
294
295 if(!aod && AODEvent() && IsStandardAOD()) {
296 // In case there is an AOD handler writing a standard AOD, use the AOD
297 // event in memory rather than the input (ESD) event.
298 aod = dynamic_cast<AliAODEvent*> (AODEvent());
299 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
300 // have to taken from the AOD event hold by the AliAODExtension
301 AliAODHandler* aodHandler = (AliAODHandler*)
302 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
303 if(aodHandler->GetExtensions()) {
304
305 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
306 AliAODEvent *aodFromExt = ext->GetAOD();
307
308
309 switch(fDecChannel){
310 case 0:
311 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
312 break;
313 case 1:
314 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
315 break;
316 case 2:
317 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
318 break;
319 case 3:
320 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
321 break;
322 case 4:
323 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
324 break;
325 case 5:
326 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
327 break;
328 }
329 }
330 } else {
331 switch(fDecChannel){
332 case 0:
333 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
334 break;
335 case 1:
336 arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
337 break;
338 case 2:
339 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
340 break;
341 case 3:
342 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
343 break;
344 case 4:
345 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
346 break;
347 case 5:
348 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
349 break;
350 }
351 }
352 if(!arrayProng) {
353 printf("AliAnalysisTaskSESignificance::UserExec: branch not found!\n");
354 return;
355 }
356
357 TClonesArray *arrayMC=0;
358 AliAODMCHeader *mcHeader=0;
359
360 // load MC particles
361 if(fReadMC){
362
363 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
364 if(!arrayMC) {
365 printf("AliAnalysisTaskSESignificance::UserExec: MC particles branch not found!\n");
366 // return;
367 }
368
369 // load MC header
370 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
371 if(!mcHeader) {
372 printf("AliAnalysisTaskSESignificance::UserExec: MC header branch not found!\n");
373 return;
374 }
375 }
376
377 fHistNEvents->Fill(0); // count event
378
379 AliAODRecoDecayHF *d=0;
380 Int_t nprongs=1;
381 Int_t *pdgdaughters=0x0;
382 Int_t prongpdg=211;
383
384
385 switch(fDecChannel){
386 case 0:
387 //D+
388 pdgdaughters =new Int_t[3];
389 pdgdaughters[0]=211;//pi
390 pdgdaughters[1]=321;//K
391 pdgdaughters[2]=211;//pi
392 nprongs=3;
393 prongpdg=411;
394 break;
395 case 1:
396 //D0
397 pdgdaughters =new Int_t[2];
398 pdgdaughters[0]=211;//pi
399 pdgdaughters[1]=321;//K
400 nprongs=2;
401 prongpdg=421;
402 break;
403 case 2:
404 //D*
405 pdgdaughters =new Int_t[3];
406 pdgdaughters[1]=211;//pi
407 pdgdaughters[0]=321;//K
408 pdgdaughters[2]=211;//pi (soft?)
409 nprongs=3;
410 prongpdg=413;
411 break;
412 case 3:
413 //Ds
414 pdgdaughters =new Int_t[3];
415 pdgdaughters[0]=211;//pi
416 pdgdaughters[1]=321;//K
417 pdgdaughters[2]=321;//K
418 nprongs=3;
419 prongpdg=431;
420 break;
421 case 4:
422 //D0
423 pdgdaughters =new Int_t[4];
424 pdgdaughters[0]=321;
425 pdgdaughters[1]=211;
426 pdgdaughters[2]=211;
427 pdgdaughters[3]=211;
428 nprongs=4;
429 prongpdg=421;
430 break;
431 case 5:
432 //Lambda_c
433 pdgdaughters =new Int_t[3];
434 //check order
435 pdgdaughters[0]=2212;//p
436 pdgdaughters[1]=321;//K
437 pdgdaughters[2]=211;//pi
438 nprongs=3;
439 prongpdg=4122;
440 break;
441 }
442
443 Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
444 Int_t nProng = arrayProng->GetEntriesFast();
445 if(fDebug>1) printf("Number of D2H: %d\n",nProng);
446
5ecc6e9b 447 if(!fRDCuts->IsEventSelected(aod)){
448 fHistNEvents->Fill(1);
449 return;
450 }
451
cbedddce 452 for (Int_t iProng = 0; iProng < nProng; iProng++) {
453 d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
454 Double_t invMass=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
455 Double_t invMassC=0;
456
457 Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel);
458
459 if(isSelected) {
5ecc6e9b 460 fHistNEvents->Fill(2); // count selected with loosest cuts
cbedddce 461 if(fDebug>1) printf("Is Selected\n");
462
463 const Int_t nvars=fRDCuts->GetNVarsForOpt();
464 Float_t *vars = new Float_t[nvars];
465 Int_t nVals=0;
466
467 fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
468 Int_t ptbin=fRDCuts->PtBin(d->Pt());
469 if(ptbin==-1) continue;
470 TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
471 ULong64_t *addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
472 if(fDebug>1)printf("nvals = %d\n",nVals);
cbedddce 473 for(Int_t ivals=0;ivals<nVals;ivals++){
474 if(addresses[ivals]>=((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetNTotCells()){
475 if (fDebug>1) printf("Overflow!!\n");
476 return;
477 }
478 fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
5ecc6e9b 479 fHistNEvents->Fill(3);
cbedddce 480 if(fReadMC){
481 Int_t lab=-1;
482 lab = d->MatchToMC(prongpdg,arrayMC,nprongs,pdgdaughters);
483 if(lab>=0){ //signal
484 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
485 Int_t pdgMC = dMC->GetPdgCode();
486
487 if(pdgMC==+prongpdg) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
488 else fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
489 }
490 else{ //background
491 fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMass);
492 }
493
494 }
495 }
496
497
498 if(fDecChannel==AliAnalysisTaskSESignificance::kD0toKpi){ //repeat the cycle to checks cuts passed as D0bar
499 pdgdaughters[0]=321;
500 pdgdaughters[1]=211;
501 invMass=d->InvMass(nprongs,(UInt_t*)pdgdaughters);
502 fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
503 addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
504 for(Int_t ivals=0;ivals<nVals;ivals++){
505 fMassHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
506 if(fReadMC){
507 Int_t lab=-1;
508 lab = d->MatchToMC(-prongpdg,arrayMC,nprongs,pdgdaughters);
509 if(lab>=0){ //signal
510 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(lab);
511 Int_t pdgMC = dMC->GetPdgCode();
512 if(pdgMC==-prongpdg) fSigHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
513 else fRflHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
514 }
515 else{ //background
516 fBkgHist[ptbin*nHistpermv+addresses[ivals]]->Fill(invMassC);
517 }
518
519 }
520 }
521 }
522
523 }// end if selected
524
525
526 }
527
528 PostData(1,fOutput);
529 return;
530}
531
532
533//________________________________________________________________________
534void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
535{
536 // Terminate analysis
537 //
538 if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
539
540
541 fOutput = dynamic_cast<TList*> (GetOutputData(1));
542 if (!fOutput) {
543 printf("ERROR: fOutput not available\n");
544 return;
545 }
546
547 fCutList = dynamic_cast<TList*> (GetOutputData(2));
548 if (!fCutList) {
549 printf("ERROR: fCutList not available\n");
550 return;
551 }
552
553 AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
554 if (!mdvtmp){
555 cout<<"multidimvec not found in TList"<<endl;
556 fCutList->ls();
557 return;
558 }
559 Int_t nHist=mdvtmp->GetNTotCells();
560 TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
561 Bool_t drawn=kFALSE;
562 for(Int_t i=0;i<nHist;i++){
563
564 TString hisname;
565 TString signame;
566 TString bkgname;
567 TString rflname;
568
569 hisname.Form("hMass_%d",i);
570 signame.Form("hSig_%d",i);
571 bkgname.Form("hBkg_%d",i);
572 rflname.Form("hRfl_%d",i);
573
574 fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
575
576 if (!drawn && fMassHist[i]->GetEntries() > 0){
577 c1->cd();
578 fMassHist[i]->Draw();
579 drawn=kTRUE;
580 }
581
582 if(fReadMC){
583 fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
584 fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
585 fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
586 }
587
588 }
589
590
591
592
593 return;
594}
595