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