]> 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
Coverity
[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
83e7d8d8 35#include <AliLog.h>
cbedddce 36#include "AliAnalysisManager.h"
37#include "AliAODHandler.h"
38#include "AliAODEvent.h"
39#include "AliAODVertex.h"
40#include "AliAODTrack.h"
41#include "AliAODMCHeader.h"
42#include "AliAODMCParticle.h"
43#include "AliAODRecoDecayHF3Prong.h"
44#include "AliAODRecoDecayHF.h"
45#include "AliAODRecoDecayHF2Prong.h"
46#include "AliAODRecoDecayHF4Prong.h"
47#include "AliAODRecoCascadeHF.h"
48
49#include "AliAnalysisTaskSE.h"
50#include "AliRDHFCutsDplustoKpipi.h"
83e7d8d8 51#include "AliRDHFCutsD0toKpipipi.h"
52#include "AliRDHFCutsDstoKKpi.h"
53#include "AliRDHFCutsDStartoKpipi.h"
cbedddce 54#include "AliRDHFCutsD0toKpi.h"
83e7d8d8 55#include "AliRDHFCutsLctopKpi.h"
cbedddce 56#include "AliMultiDimVector.h"
57
58#include "AliAnalysisTaskSESignificance.h"
59
60ClassImp(AliAnalysisTaskSESignificance)
61
62
63//________________________________________________________________________
64AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance():
65 AliAnalysisTaskSE(),
66 fOutput(0),
67 fCutList(0),
68 fHistNEvents(0),
69 fUpmasslimit(1.965),
70 fLowmasslimit(1.765),
71 fRDCuts(0),
72 fNPtBins(0),
73 fReadMC(kFALSE),
33c5731e 74 fBFeedDown(kBoth),
cbedddce 75 fDecChannel(0),
f80a6da9 76 fSelectionlevel(0),
77 fNBins(100),
78 fPartOrAndAntiPart(0)
cbedddce 79{
80 // Default constructor
81}
82
83//________________________________________________________________________
84AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
85 AliAnalysisTaskSE(name),
86 fOutput(0),
87 fCutList(listMDV),
88 fHistNEvents(0),
89 fUpmasslimit(0),
90 fLowmasslimit(0),
91 fRDCuts(rdCuts),
92 fNPtBins(0),
93 fReadMC(kFALSE),
33c5731e 94 fBFeedDown(kBoth),
cbedddce 95 fDecChannel(decaychannel),
f80a6da9 96 fSelectionlevel(selectionlevel),
97 fNBins(100),
98 fPartOrAndAntiPart(0)
cbedddce 99{
100
101 Int_t pdg=421;
102 switch(fDecChannel){
103 case 0:
104 pdg=411;
105 break;
106 case 1:
107 pdg=421;
108 break;
109 case 2:
110 pdg=413;
111 break;
112 case 3:
113 pdg=431;
114 break;
115 case 4:
116 pdg=421;
117 break;
118 case 5:
119 pdg=4122;
120 break;
121 }
122
f80a6da9 123 SetMassLimits(0.15,pdg); //check range
cbedddce 124 fNPtBins=fRDCuts->GetNPtBins();
83e7d8d8 125
cbedddce 126 if(fDebug>1)fRDCuts->PrintAll();
127 // Output slot #1 writes into a TList container
128 DefineOutput(1,TList::Class()); //My private output
129 DefineOutput(2,TList::Class());
795316ee 130 DefineOutput(3,AliRDHFCuts::Class()); //class of the cuts
3aa6ce53 131 CheckConsistency();
cbedddce 132}
133
134 //________________________________________________________________________
135AliAnalysisTaskSESignificance::~AliAnalysisTaskSESignificance()
136{
137 // Destructor
138 if (fOutput) {
139 delete fOutput;
140 fOutput = 0;
141 }
142 if (fCutList) {
143 delete fCutList;
144 fCutList = 0;
145 }
146 if(fHistNEvents){
147 delete fHistNEvents;
148 fHistNEvents=0;
149 }
150
151/*
152 if(fRDCuts) {
153 delete fRDCuts;
154 fRDCuts= 0;
155 }
156*/
157
3aa6ce53 158}
159//_________________________________________________________________
160Bool_t AliAnalysisTaskSESignificance::CheckConsistency(){
161
162 Bool_t result = kTRUE;
163
164 const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
165 //Float_t *vars = new Float_t[nvars];
166 Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
167 Bool_t *uppervars = fRDCuts->GetIsUpperCut();
168 TString *names = fRDCuts->GetVarNames();
169
170 for(Int_t i=0;i<fNPtBins;i++){
171 TString mdvname=Form("multiDimVectorPtBin%d",i);
172 Int_t ic=0;
173
174 for(Int_t ivar=0;ivar<nvars;ivar++){
175 if(varsforopt[ivar]){
176 Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
177 Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
178 if(min==max){
e4c00a7b 179 AliFatal(Form("tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i));
180 return kFALSE;
3aa6ce53 181 }
182 Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
183 if(uppervars[ivar]&&lowermdv){
e4c00a7b 184 AliFatal(Form("%s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---please check your cuts \n ",names[ivar].Data(),i));
185 return kFALSE;
3aa6ce53 186 }
187 if(!uppervars[ivar]&&!lowermdv){
e4c00a7b 188 AliFatal(Form("%s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---please check your cuts \n",names[ivar].Data(),i));
189 return kFALSE;
3aa6ce53 190 }
191 ic++;
192 }
193 }
194 }
195 return result;
196}
cbedddce 197//_________________________________________________________________
33c5731e 198void AliAnalysisTaskSESignificance::SetBFeedDown(FeedDownEnum flagB){
199 if(fReadMC)fBFeedDown=flagB;
200 else {
795316ee 201 AliInfo("B feed down not allowed without MC info\n");
202 fBFeedDown=kBoth;
33c5731e 203 }
204}
205//_________________________________________________________________
cbedddce 206void AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
207 Float_t mass=0;
208 Int_t abspdg=TMath::Abs(pdg);
5ecc6e9b 209 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
cbedddce 210 fUpmasslimit = mass+range;
211 fLowmasslimit = mass-range;
212}
213//_________________________________________________________________
214void AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
215 if(uplimit>lowlimit)
216 {
f80a6da9 217 fUpmasslimit = uplimit;
218 fLowmasslimit = lowlimit;
cbedddce 219 }
220}
221
222
223
224//________________________________________________________________________
225void AliAnalysisTaskSESignificance::LocalInit()
226{
227 // Initialization
228
229 if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
230
795316ee 231 switch(fDecChannel){
232 case 0:
233 {
234 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
235 // Post the data
236 PostData(3,copycut);
237 }
238 break;
239 case 1:
240 {
241 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
242 // Post the data
243 PostData(3,copycut);
244 }
245 break;
246 case 2:
247 {
248 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
249 // Post the data
250 PostData(3,copycut);
251 }
252 break;
253 case 3:
254 {
255 AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fRDCuts)));
256 // Post the data
257 PostData(3,copycut);
258 }
259 break;
260 case 4:
261 {
262 AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fRDCuts)));
263 // Post the data
264 PostData(3,copycut);
265 }
266 break;
267 case 5:
268 {
269 AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fRDCuts)));
270 // Post the data
271 PostData(3,copycut);
272 }
273 break;
274
275 default:
276 return;
277 }
278
cbedddce 279 TList *mdvList = new TList();
280 mdvList->SetOwner();
281 mdvList = fCutList;
795316ee 282
cbedddce 283 PostData(2,mdvList);
284
795316ee 285
cbedddce 286}
287//________________________________________________________________________
288void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
289{
290 // Create the output container
291
292 if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
293
294 // Several histograms are more conveniently managed in a TList
295 fOutput = new TList();
296 fOutput->SetOwner();
297 fOutput->SetName("OutputHistos");
298
299 //same number of steps in each multiDimVectorPtBin%d !
300 Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
83e7d8d8 301 cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
cbedddce 302 nHist=nHist*fNPtBins;
83e7d8d8 303 cout<<"Total = "<<nHist<<endl;
cbedddce 304 for(Int_t i=0;i<nHist;i++){
305
306 TString hisname;
307 TString signame;
308 TString bkgname;
309 TString rflname;
310 TString title;
311
312 hisname.Form("hMass_%d",i);
313 signame.Form("hSig_%d",i);
314 bkgname.Form("hBkg_%d",i);
315 rflname.Form("hRfl_%d",i);
316
317 title.Form("Invariant mass;M[GeV/c^{2}];Entries");
318
f80a6da9 319 fMassHist[i]=new TH1F(hisname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
cbedddce 320 fMassHist[i]->Sumw2();
321 fOutput->Add(fMassHist[i]);
322
323 if(fReadMC){
f80a6da9 324 fSigHist[i]=new TH1F(signame.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
cbedddce 325 fSigHist[i]->Sumw2();
326 fOutput->Add(fSigHist[i]);
327
f80a6da9 328 fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
cbedddce 329 fBkgHist[i]->Sumw2();
330 fOutput->Add(fBkgHist[i]);
331
f80a6da9 332 if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi){
333 fRflHist[i]=new TH1F(rflname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
334 fRflHist[i]->Sumw2();
335 fOutput->Add(fRflHist[i]);
336 }
cbedddce 337 }
338 }
339
795316ee 340 fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",8,-0.5,7.5);
cbedddce 341 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
dc850ba8 342 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvSelected (vtx)");
5ecc6e9b 343 fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
344 fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
d52f7b50 345 fHistNEvents->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
dc850ba8 346 fHistNEvents->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
95af3e50 347 if(fReadMC){
348 fHistNEvents->GetXaxis()->SetBinLabel(7,"MC Cand from c");
349 fHistNEvents->GetXaxis()->SetBinLabel(8,"MC Cand from b");
350 } else{
351 fHistNEvents->GetXaxis()->SetBinLabel(7,"N candidates");
352 }
cbedddce 353 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
354 fOutput->Add(fHistNEvents);
355
356
357 return;
358}
359
360//________________________________________________________________________
361void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
362{
363 // Execute analysis for current event:
364 // heavy flavor candidates association to MC truth
365
366 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
83e7d8d8 367 if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
cbedddce 368 // Post the data already here
369 PostData(1,fOutput);
370 TClonesArray *arrayProng =0;
371
372 if(!aod && AODEvent() && IsStandardAOD()) {
373 // In case there is an AOD handler writing a standard AOD, use the AOD
374 // event in memory rather than the input (ESD) event.
375 aod = dynamic_cast<AliAODEvent*> (AODEvent());
376 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
377 // have to taken from the AOD event hold by the AliAODExtension
378 AliAODHandler* aodHandler = (AliAODHandler*)
379 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
380 if(aodHandler->GetExtensions()) {
381
382 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
383 AliAODEvent *aodFromExt = ext->GetAOD();
384
385
386 switch(fDecChannel){
387 case 0:
388 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
389 break;
390 case 1:
391 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
392 break;
393 case 2:
394 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
395 break;
396 case 3:
397 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
398 break;
399 case 4:
400 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
401 break;
402 case 5:
403 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
404 break;
405 }
406 }
dc222f77 407 } else if(aod) {
cbedddce 408 switch(fDecChannel){
409 case 0:
410 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
411 break;
412 case 1:
413 arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
414 break;
415 case 2:
416 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
417 break;
418 case 3:
419 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
420 break;
421 case 4:
422 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
423 break;
424 case 5:
425 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
426 break;
427 }
428 }
dc222f77 429 if(!aod || !arrayProng) {
f80a6da9 430 AliError("AliAnalysisTaskSESignificance::UserExec:Branch not found!\n");
cbedddce 431 return;
432 }
433
7c23877d 434 // fix for temporary bug in ESDfilter
435 // the AODs with null vertex pointer didn't pass the PhysSel
f80a6da9 436 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
cbedddce 437 TClonesArray *arrayMC=0;
438 AliAODMCHeader *mcHeader=0;
439
440 // load MC particles
441 if(fReadMC){
442
443 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
444 if(!arrayMC) {
dc222f77 445 AliError("AliAnalysisTaskSESignificance::UserExec:MC particles branch not found!\n");
446 return;
cbedddce 447 }
448
449 // load MC header
450 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
451 if(!mcHeader) {
f80a6da9 452 AliError("AliAnalysisTaskSESignificance::UserExec:MC header branch not found!\n");
cbedddce 453 return;
454 }
455 }
456
457 fHistNEvents->Fill(0); // count event
458
459 AliAODRecoDecayHF *d=0;
460 Int_t nprongs=1;
461 Int_t *pdgdaughters=0x0;
f80a6da9 462 Int_t absPdgMom=411;
cbedddce 463
464
465 switch(fDecChannel){
466 case 0:
467 //D+
468 pdgdaughters =new Int_t[3];
469 pdgdaughters[0]=211;//pi
470 pdgdaughters[1]=321;//K
471 pdgdaughters[2]=211;//pi
472 nprongs=3;
f80a6da9 473 absPdgMom=411;
cbedddce 474 break;
475 case 1:
476 //D0
477 pdgdaughters =new Int_t[2];
478 pdgdaughters[0]=211;//pi
479 pdgdaughters[1]=321;//K
480 nprongs=2;
f80a6da9 481 absPdgMom=421;
cbedddce 482 break;
483 case 2:
484 //D*
485 pdgdaughters =new Int_t[3];
486 pdgdaughters[1]=211;//pi
487 pdgdaughters[0]=321;//K
488 pdgdaughters[2]=211;//pi (soft?)
489 nprongs=3;
f80a6da9 490 absPdgMom=413;
cbedddce 491 break;
492 case 3:
493 //Ds
494 pdgdaughters =new Int_t[3];
4d829de7 495 pdgdaughters[0]=321;//K
cbedddce 496 pdgdaughters[1]=321;//K
4d829de7 497 pdgdaughters[2]=211;//pi
cbedddce 498 nprongs=3;
f80a6da9 499 absPdgMom=431;
cbedddce 500 break;
501 case 4:
f80a6da9 502 //D0 in 4 prongs
cbedddce 503 pdgdaughters =new Int_t[4];
504 pdgdaughters[0]=321;
505 pdgdaughters[1]=211;
506 pdgdaughters[2]=211;
507 pdgdaughters[3]=211;
508 nprongs=4;
f80a6da9 509 absPdgMom=421;
cbedddce 510 break;
511 case 5:
512 //Lambda_c
513 pdgdaughters =new Int_t[3];
cbedddce 514 pdgdaughters[0]=2212;//p
515 pdgdaughters[1]=321;//K
516 pdgdaughters[2]=211;//pi
517 nprongs=3;
f80a6da9 518 absPdgMom=4122;
cbedddce 519 break;
520 }
521
522 Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
523 Int_t nProng = arrayProng->GetEntriesFast();
524 if(fDebug>1) printf("Number of D2H: %d\n",nProng);
525
dc850ba8 526 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
527 TString trigclass=aod->GetFiredTriggerClasses();
528 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(5);
529
fb9311e8 530 if(fRDCuts->IsEventSelected(aod)) {
531 fHistNEvents->Fill(1);
532 }else{
d52f7b50 533 if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
534 fHistNEvents->Fill(4);
fb9311e8 535 delete [] pdgdaughters;
5ecc6e9b 536 return;
537 }
dc850ba8 538
cbedddce 539 for (Int_t iProng = 0; iProng < nProng; iProng++) {
95af3e50 540 fHistNEvents->Fill(6);
cbedddce 541 d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
cbedddce 542
f80a6da9 543 Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
a4bb4427 544 Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
33c5731e 545
795316ee 546 if(fReadMC && fBFeedDown!=kBoth && isSelected){
33c5731e 547 Int_t labD = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
548 if(labD>=0){
549 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
550 Int_t label=partD->GetMother();
551 AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(label);
552 while(label>=0){//get first mother
553 mot = (AliAODMCParticle*)arrayMC->At(label);
554 label=mot->GetMother();
555 }
556 Int_t pdgMotCode = mot->GetPdgCode();
795316ee 557
33c5731e 558 if(TMath::Abs(pdgMotCode)<=4){
795316ee 559 fHistNEvents->Fill(6);
33c5731e 560 if(fBFeedDown==kBeautyOnly)isSelected=kFALSE; //from primary charm
561 }else{
795316ee 562 fHistNEvents->Fill(7);
563 if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
33c5731e 564 }
795316ee 565
33c5731e 566 }
567 }
cbedddce 568
33c5731e 569
f80a6da9 570 if(isSelected&&isFidAcc) {
5ecc6e9b 571 fHistNEvents->Fill(2); // count selected with loosest cuts
83e7d8d8 572 if(fDebug>1) printf("+++++++Is Selected\n");
cbedddce 573
f80a6da9 574 Double_t* invMass=0x0;
575 Int_t nmasses;
576 CalculateInvMasses(d,invMass,nmasses);
577
cbedddce 578 const Int_t nvars=fRDCuts->GetNVarsForOpt();
579 Float_t *vars = new Float_t[nvars];
580 Int_t nVals=0;
581
582 fRDCuts->GetCutVarsForOpt(d,vars,nvars,pdgdaughters);
583 Int_t ptbin=fRDCuts->PtBin(d->Pt());
584 if(ptbin==-1) continue;
585 TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
586 ULong64_t *addresses = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGlobalAddressesAboveCuts(vars,(Float_t)d->Pt(),nVals);
587 if(fDebug>1)printf("nvals = %d\n",nVals);
cbedddce 588 for(Int_t ivals=0;ivals<nVals;ivals++){
589 if(addresses[ivals]>=((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetNTotCells()){
590 if (fDebug>1) printf("Overflow!!\n");
dc222f77 591 delete [] addresses;
fb9311e8 592 delete [] pdgdaughters;
cbedddce 593 return;
594 }
cbedddce 595
f80a6da9 596 fHistNEvents->Fill(3);
cbedddce 597
f80a6da9 598 //fill the histograms with the appropriate method
599 switch (fDecChannel){
600 case 0:
601 FillDplus(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),invMass,isSelected);
602 break;
603 case 1:
604 FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),invMass,isSelected);
605 break;
606 case 2:
607 FillDstar(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),invMass,isSelected);
608 break;
609 case 3:
610 FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),invMass,isSelected);
611 break;
612 case 4:
613 FillD04p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),invMass,isSelected);
614 break;
615 case 5:
616 FillLambdac(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),invMass,isSelected);
617 break;
618 default:
619 break;
cbedddce 620 }
f80a6da9 621
cbedddce 622 }
dc222f77 623 delete [] addresses;
cbedddce 624 }// end if selected
cbedddce 625
626 }
f80a6da9 627
e11ae259 628 if(pdgdaughters) {delete [] pdgdaughters; pdgdaughters=NULL;}
629
cbedddce 630 PostData(1,fOutput);
631 return;
632}
633
f80a6da9 634//***************************************************************************
635
636// Methods used in the UserExec
637
638void AliAnalysisTaskSESignificance::CalculateInvMasses(AliAODRecoDecayHF* d,Double_t*& masses,Int_t& nmasses){
639 //Calculates all the possible invariant masses for each candidate
640 //NB: the implementation for each candidate is responsibility of the corresponding developer
641
642 switch(fDecChannel){
643 case 0:
644 //D+ -- Giacomo, Renu ?
645 {
646 nmasses=1;
647 masses=new Double_t[nmasses];
648 Int_t pdgdaughters[3] = {211,321,211};
649 masses[0]=d->InvMass(3,(UInt_t*)pdgdaughters);
650 }
651 break;
652 case 1:
653 //D0 (Kpi) -- Chiara
654 {
655 const Int_t ndght=2;
656 nmasses=2;
657 masses=new Double_t[nmasses];
658 Int_t pdgdaughtersD0[ndght]={211,321};//pi,K
659 masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0); //D0
660 Int_t pdgdaughtersD0bar[ndght]={321,211};//K,pi
661 masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0bar); //D0bar
662 }
663 break;
664 case 2:
665 //D* -- ? Yifei, Alessandro
666 {
667 //.....
668 }
669 break;
670 case 3:
671 //Ds -- Sergey, Sahdana
672 {
4d829de7 673 const Int_t ndght=3;
674 nmasses=2;
675 masses=new Double_t[nmasses];
676 Int_t pdgdaughtersDsKKpi[ndght]={321,321,211};//K,K,pi
677 masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersDsKKpi); //Ds
678 Int_t pdgdaughtersDspiKK[ndght]={211,321,321};//pi,K,K
679 masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersDspiKK); //Dsbar
680
f80a6da9 681 //.....
682 }
683 break;
684 case 4:
685 //D0 (Kpipipi) -- ? Rossella, Fabio ???
686 {
687 //.....
688 }
689 break;
690 case 5:
691 //Lambda_c -- Rossella
692 {
693 //.....
694 }
695 break;
696 default:
697 break;
698 }
699}
700
701//********************************************************************************************
702
703//Methods to fill the istograms with MC information, one for each candidate
704//NB: the implementation for each candidate is responsibility of the corresponding developer
705
706void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Double_t* masses,Int_t isSel){
707 //D+ channel
708 if(!isSel){
709 AliError("Candidate not selected\n");
710 return;
711 }
712
713 if(fPartOrAndAntiPart*d->GetCharge()<0)return;
714
715 fMassHist[index]->Fill(masses[0]);
716
717 Int_t pdgdaughters[3] = {211,321,211};
718
719 if(fReadMC){
720 Int_t lab=-1;
721 lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
722 if(lab>=0){ //signal
723 fSigHist[index]->Fill(masses[0]);
724 } else{ //background
725 fBkgHist[index]->Fill(masses[0]);
726 }
727 }
728}
729
730
731void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Double_t* masses,Int_t isSel){
732
733 //D0->Kpi channel
734
735 //mass histograms
736 if(!masses){
737 AliError("Masses not calculated\n");
738 return;
739 }
740 if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(masses[0]);
741 if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(masses[1]);
742
743
4d829de7 744
f80a6da9 745 //MC histograms
746 if(fReadMC){
747
748 Int_t matchtoMC=-1;
749
750 //D0
751 Int_t pdgdaughters[2];
752 pdgdaughters[0]=211;//pi
753 pdgdaughters[1]=321;//K
754 Int_t nprongs=2;
755 Int_t absPdgMom=421;
756
757 matchtoMC = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
758
759 Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
760 if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0){ //D0
761 if(matchtoMC>=0){
762 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
763 Int_t pdgMC = dMC->GetPdgCode();
764
765 if(pdgMC==prongPdgPlus) fSigHist[index]->Fill(masses[0]);
766 else fRflHist[index]->Fill(masses[0]);
767
768 } else fBkgHist[index]->Fill(masses[0]);
769
770 }
771 if(isSel>=2 && fPartOrAndAntiPart<=0){ //D0bar
772 if(matchtoMC>=0){
773 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
774 Int_t pdgMC = dMC->GetPdgCode();
775
776 if(pdgMC==prongPdgMinus) fSigHist[index]->Fill(masses[1]);
777 else fRflHist[index]->Fill(masses[1]);
778 } else fBkgHist[index]->Fill(masses[1]);
779 }
780 }
781}
782
783void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Double_t* /*masses*/,Int_t /*matchtoMC*/){
784 //D* channel
785
786 AliInfo("Dstar channel not implemented\n");
787
788}
789
790
4d829de7 791void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Double_t* masses,Int_t isSel){
792
793 //AliInfo("Ds channel not implemented\n");
794
795
796 if(!masses){
797 AliError("Masses not calculated\n");
798 return;
799 }
800
801 Int_t pdgDstoKKpi[3]={321,321,211};
802 Int_t labDs=-1;
803 if(fReadMC){
804 labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
805 }
806
807
808 if(isSel&1 && fPartOrAndAntiPart*d->GetCharge()>=0) {
809
810 fMassHist[index]->Fill(masses[0]);
811
812 if(fReadMC){
f80a6da9 813
4d829de7 814 if(labDs>=0){
815 Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
816 AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
817 Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
818
819 if(pdgCode0==321) {
820
821 fSigHist[index]->Fill(masses[0]); //signal
822 }else{
823 fRflHist[index]->Fill(masses[0]); //Reflected signal
824 }
825 }else{
826 fBkgHist[index]->Fill(masses[0]); // Background
827 }
828 }
829 }
830
831 if(isSel&2 && fPartOrAndAntiPart*d->GetCharge()>=0){
832 fMassHist[index]->Fill(masses[1]);
833 if(fReadMC){
834 if(labDs>=0){
835 Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
836 AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
837 Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
838
839 if(pdgCode0==211) {
840
841 fSigHist[index]->Fill(masses[1]);
842 }else{
843 fRflHist[index]->Fill(masses[1]);
844 }
845 }else{
846 fBkgHist[index]->Fill(masses[1]);
847 }
848 }
849 }
850
851
f80a6da9 852
853}
854
855void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Double_t* /*masses*/,Int_t /*matchtoMC*/){
856 //D0->Kpipipi channel
857 AliInfo("D0 in 4 prongs channel not implemented\n");
858
859}
860
861void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Double_t* /*masses*/,Int_t /*matchtoMC*/){
862 AliInfo("Lambdac channel not implemented\n");
863
95af3e50 864
f80a6da9 865
866}
867
cbedddce 868
869//________________________________________________________________________
870void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
871{
872 // Terminate analysis
873 //
874 if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
875
876
877 fOutput = dynamic_cast<TList*> (GetOutputData(1));
878 if (!fOutput) {
879 printf("ERROR: fOutput not available\n");
880 return;
881 }
882
883 fCutList = dynamic_cast<TList*> (GetOutputData(2));
884 if (!fCutList) {
885 printf("ERROR: fCutList not available\n");
886 return;
887 }
888
889 AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
890 if (!mdvtmp){
891 cout<<"multidimvec not found in TList"<<endl;
892 fCutList->ls();
893 return;
894 }
895 Int_t nHist=mdvtmp->GetNTotCells();
896 TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
897 Bool_t drawn=kFALSE;
898 for(Int_t i=0;i<nHist;i++){
899
900 TString hisname;
901 TString signame;
902 TString bkgname;
903 TString rflname;
904
905 hisname.Form("hMass_%d",i);
906 signame.Form("hSig_%d",i);
907 bkgname.Form("hBkg_%d",i);
908 rflname.Form("hRfl_%d",i);
909
910 fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
911
912 if (!drawn && fMassHist[i]->GetEntries() > 0){
913 c1->cd();
914 fMassHist[i]->Draw();
915 drawn=kTRUE;
916 }
917
918 if(fReadMC){
919 fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
920 fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
f80a6da9 921 if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
cbedddce 922 }
923
924 }
925
926
927
928
929 return;
930}
ab3e5f67 931//-------------------------------------------
cbedddce 932