]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGHF/vertexingHF/AliAnalysisTaskSESignificance.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Transition PWG3 --> PWGHF
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSESignificance.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
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
37#include <AliLog.h>
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"
53#include "AliRDHFCutsD0toKpipipi.h"
54#include "AliRDHFCutsDstoKKpi.h"
55#include "AliRDHFCutsDStartoKpipi.h"
56#include "AliRDHFCutsD0toKpi.h"
57#include "AliRDHFCutsLctopKpi.h"
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),
76 fUseSelBit(kFALSE),
77 fBFeedDown(kBoth),
78 fDecChannel(0),
79 fPDGmother(0),
80 fNProngs(0),
81 fBranchName(""),
82 fSelectionlevel(0),
83 fNVars(0),
84 fNBins(100),
85 fPartOrAndAntiPart(0),
86 fDsChannel(0)
87{
88 // Default constructor
89 SetPDGCodes();
90 SetDsChannel(kPhi);
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
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),
115 fUseSelBit(kFALSE),
116 fBFeedDown(kBoth),
117 fDecChannel(decaychannel),
118 fPDGmother(0),
119 fNProngs(0),
120 fBranchName(""),
121 fSelectionlevel(selectionlevel),
122 fNVars(0),
123 fNBins(100),
124 fPartOrAndAntiPart(0),
125 fDsChannel(0)
126{
127
128 SetPDGCodes();
129 SetDsChannel(kPhi);
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 }
136 fNPtBins=fRDCuts->GetNPtBins();
137
138 fNVars=fRDCuts->GetNVarsForOpt();
139 if(fNVars>kMaxCutVar) AliFatal(Form("Too large number of cut variables, maximum is %d",kMaxCutVar));
140
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());
145 DefineOutput(3,AliRDHFCuts::Class()); //class of the cuts
146 CheckConsistency();
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 }
165/*
166 if(fRDCuts) {
167 delete fRDCuts;
168 fRDCuts= 0;
169 }
170*/
171
172}
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
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){
261 AliFatal(Form("tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i));
262 return kFALSE;
263 }
264 Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
265 if(uppervars[ivar]&&lowermdv){
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;
268 }
269 if(!uppervars[ivar]&&!lowermdv){
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;
272 }
273 ic++;
274 }
275 }
276 }
277 return result;
278}
279//_________________________________________________________________
280void AliAnalysisTaskSESignificance::SetBFeedDown(FeedDownEnum flagB){
281 if(fReadMC)fBFeedDown=flagB;
282 else {
283 AliInfo("B feed down not allowed without MC info\n");
284 fBFeedDown=kBoth;
285 }
286}
287//_________________________________________________________________
288void AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
289 Float_t mass=0;
290 Int_t abspdg=TMath::Abs(pdg);
291 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
292 fUpmasslimit = mass+range;
293 fLowmasslimit = mass-range;
294}
295//_________________________________________________________________
296void AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
297 if(uplimit>lowlimit)
298 {
299 fUpmasslimit = uplimit;
300 fLowmasslimit = lowlimit;
301 }
302}
303
304
305
306//________________________________________________________________________
307void AliAnalysisTaskSESignificance::LocalInit()
308{
309 // Initialization
310
311 if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
312
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
361 TList *mdvList = new TList();
362 mdvList->SetOwner();
363 mdvList = fCutList;
364
365 PostData(2,mdvList);
366
367
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();
383 cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
384 nHist=nHist*fNPtBins;
385 cout<<"Total = "<<nHist<<endl;
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
401 fMassHist[i]=new TH1F(hisname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
402 fMassHist[i]->Sumw2();
403 fOutput->Add(fMassHist[i]);
404
405 if(fReadMC){
406 fSigHist[i]=new TH1F(signame.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
407 fSigHist[i]->Sumw2();
408 fOutput->Add(fSigHist[i]);
409
410 fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
411 fBkgHist[i]->Sumw2();
412 fOutput->Add(fBkgHist[i]);
413
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 }
419 }
420 }
421
422 fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",8,-0.5,7.5);
423 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
424 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvSelected (vtx)");
425 fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
426 fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
427 fHistNEvents->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
428 fHistNEvents->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
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 }
435 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
436 fOutput->Add(fHistNEvents);
437
438
439 PostData(1,fOutput);
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());
450 if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
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();
467 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
468
469 }
470 } else if(aod) {
471 arrayProng=(TClonesArray*)aod->GetList()->FindObject(fBranchName.Data());
472 }
473 if(!aod || !arrayProng) {
474 AliError("AliAnalysisTaskSESignificance::UserExec:Branch not found!\n");
475 return;
476 }
477
478 // fix for temporary bug in ESDfilter
479 // the AODs with null vertex pointer didn't pass the PhysSel
480 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
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) {
489 AliError("AliAnalysisTaskSESignificance::UserExec:MC particles branch not found!\n");
490 return;
491 }
492
493 // load MC header
494 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
495 if(!mcHeader) {
496 AliError("AliAnalysisTaskSESignificance::UserExec:MC header branch not found!\n");
497 return;
498 }
499 }
500
501 fHistNEvents->Fill(0); // count event
502
503 AliAODRecoDecayHF *d=0;
504
505
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
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
515 if(fRDCuts->IsEventSelected(aod)) {
516 fHistNEvents->Fill(1);
517 }else{
518 if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
519 fHistNEvents->Fill(4);
520 return;
521 }
522
523 for (Int_t iProng = 0; iProng < nProng; iProng++) {
524 fHistNEvents->Fill(6);
525 d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
526
527 AliAODRecoCascadeHF* DStarToD0pi = NULL;
528 AliAODRecoDecayHF2Prong* D0Particle = NULL;
529 fPDGDStarToD0pi[0] = 421; fPDGDStarToD0pi[1] = 211;
530 fPDGD0ToKpi[0] = 321; fPDGD0ToKpi[1] = 211;
531
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
554 if (fDecChannel==2) {
555 DStarToD0pi = (AliAODRecoCascadeHF*)arrayProng->At(iProng);
556 if (!DStarToD0pi->GetSecondaryVtx()) continue;
557 D0Particle = (AliAODRecoDecayHF2Prong*)DStarToD0pi->Get2Prong();
558 if (!D0Particle) continue;
559 }
560
561 Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(fPDGmother));
562 Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
563
564 if(fReadMC && fBFeedDown!=kBoth && isSelected){
565 Int_t labD = d->MatchToMC(fPDGmother,arrayMC,fNProngs,fPDGdaughters);
566 if(labD>=0){
567 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
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));
573 fHistNEvents->Fill(7);
574 if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
575 }
576 else {
577 //prompt particle
578 fHistNEvents->Fill(6);
579 if(fBFeedDown==kBeautyOnly)isSelected=kFALSE;
580 }
581 }
582 }
583
584 if(isSelected&&isFidAcc) {
585 fHistNEvents->Fill(2); // count selected with loosest cuts
586 if(fDebug>1) printf("+++++++Is Selected\n");
587
588 Int_t nVals=0;
589 if(fDecChannel==3) SetPDGdaughterDstoKKpi();
590 fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
591 Int_t ptbin=fRDCuts->PtBin(d->Pt());
592 if(ptbin==-1) continue;
593 TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
594 AliMultiDimVector* muvec=(AliMultiDimVector*)fCutList->FindObject(mdvname.Data());
595
596 ULong64_t *addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
597 if(fDebug>1)printf("nvals = %d\n",nVals);
598 for(Int_t ivals=0;ivals<nVals;ivals++){
599 if(addresses[ivals]>=muvec->GetNTotCells()){
600 if (fDebug>1) printf("Overflow!!\n");
601 delete [] addresses;
602 return;
603 }
604
605 fHistNEvents->Fill(3);
606
607 //fill the histograms with the appropriate method
608 switch (fDecChannel){
609 case 0:
610 FillDplus(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
611 break;
612 case 1:
613 FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
614 break;
615 case 2:
616 FillDstar(DStarToD0pi,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
617 break;
618 case 3:
619 if(isSelected&1){
620 FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,1);
621 }
622 break;
623 case 4:
624 FillD04p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
625 break;
626 case 5:
627 FillLambdac(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
628 break;
629 default:
630 break;
631 }
632
633 }
634
635 if (fDecChannel==3 && isSelected&2){
636 SetPDGdaughterDstopiKK();
637 nVals=0;
638 fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
639 delete [] addresses;
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
654 delete [] addresses;
655 }// end if selected
656
657 }
658
659
660 PostData(1,fOutput);
661 return;
662}
663
664//***************************************************************************
665
666// Methods used in the UserExec
667
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
674void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
675 //D+ channel
676 if(!isSel){
677 AliError("Candidate not selected\n");
678 return;
679 }
680
681 if(fPartOrAndAntiPart*d->GetCharge()<0)return;
682 Int_t pdgdaughters[3] = {211,321,211};
683 Double_t mass=d->InvMass(3,(UInt_t*)pdgdaughters);
684
685 fMassHist[index]->Fill(mass);
686
687
688 if(fReadMC){
689 Int_t lab=-1;
690 lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
691 if(lab>=0){ //signal
692 fSigHist[index]->Fill(mass);
693 } else{ //background
694 fBkgHist[index]->Fill(mass);
695 }
696 }
697}
698
699
700void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
701
702 //D0->Kpi channel
703
704 //mass histograms
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
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
715
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
725
726 matchtoMC = d->MatchToMC(fPDGmother,arrayMC,fNProngs,pdgdaughters);
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
752void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel){
753 //D* channel
754
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 }
802
803}
804
805
806void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay){
807
808 // Fill Ds histos
809
810
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
817 Int_t labDs=-1;
818 if(fReadMC){
819 labDs = d->MatchToMC(431,arrayMC,3,pdgDsKKpi);
820 }
821
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 }
850 }else{
851 fBkgHist[index]->Fill(masses[0]); // Background
852 }
853 }
854 }
855 }
856
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 }
874 }else{
875 fBkgHist[index]->Fill(masses[1]);
876 }
877 }
878 }
879 }
880}
881
882void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
883 //D0->Kpipipi channel
884 AliInfo("D0 in 4 prongs channel not implemented\n");
885
886}
887
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 }
956
957
958}
959
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()));
1013 if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
1014 }
1015
1016 }
1017
1018
1019
1020
1021 return;
1022}
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//_________________________________________________________________________________________________