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