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