1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////
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 /////////////////////////////////////////////////////////////
28 #include <Riostream.h>
29 #include <TClonesArray.h>
35 #include <TDatabasePDG.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"
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"
60 #include "AliAnalysisTaskSESignificance.h"
65 ClassImp(AliAnalysisTaskSESignificance)
68 //________________________________________________________________________
69 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance():
88 fPartOrAndAntiPart(0),
91 // Default constructor
95 for(Int_t i=0;i<4;i++) fPDGdaughters[i]=0;
96 for(Int_t i=0;i<2;i++) {fPDGD0ToKpi[i]=0;fPDGDStarToD0pi[i]=0;}
97 for(Int_t i=0;i<kMaxCutVar;i++) fVars[i]=0.;
98 for(Int_t i=0;i<kMaxNHist;i++) {
107 //________________________________________________________________________
108 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
109 AliAnalysisTaskSE(name),
120 fDecChannel(decaychannel),
124 fSelectionlevel(selectionlevel),
127 fPartOrAndAntiPart(0),
133 if (fDecChannel!=2) SetMassLimits(0.15,fPDGmother); //check range
137 SetMassLimits(min, max);
139 fNPtBins=fRDCuts->GetNPtBins();
141 fNVars=fRDCuts->GetNVarsForOpt();
142 if(fNVars>kMaxCutVar) AliFatal(Form("Too large number of cut variables, maximum is %d",kMaxCutVar));
144 if(fDebug>1)fRDCuts->PrintAll();
145 // Output slot #1 writes into a TList container
146 DefineOutput(1,TList::Class()); //My private output
147 DefineOutput(2,TList::Class());
148 DefineOutput(3,AliRDHFCuts::Class()); //class of the cuts
152 //________________________________________________________________________
153 AliAnalysisTaskSESignificance::~AliAnalysisTaskSESignificance()
176 //_________________________________________________________________
177 void AliAnalysisTaskSESignificance::SetPDGCodes(){
178 // sets channel dependent quantities
185 fPDGdaughters[0]=211;//pi
186 fPDGdaughters[1]=321;//K
187 fPDGdaughters[2]=211;//pi
188 fPDGdaughters[3]=0; //empty
189 fBranchName="Charm3Prong";
195 fPDGdaughters[0]=211;//pi
196 fPDGdaughters[1]=321;//K
197 fPDGdaughters[2]=0; //empty
198 fPDGdaughters[3]=0; //empty
199 fBranchName="D0toKpi";
205 fPDGdaughters[1]=211;//pi
206 fPDGdaughters[0]=321;//K
207 fPDGdaughters[2]=211;//pi (soft?)
208 fPDGdaughters[3]=0; //empty
215 fPDGdaughters[0]=321;//K
216 fPDGdaughters[1]=321;//K
217 fPDGdaughters[2]=211;//pi
218 fPDGdaughters[3]=0; //empty
219 fBranchName="Charm3Prong";
225 fPDGdaughters[0]=321;
226 fPDGdaughters[1]=211;
227 fPDGdaughters[2]=211;
228 fPDGdaughters[3]=211;
229 fBranchName="Charm4Prong";
235 fPDGdaughters[0]=2212;//p
236 fPDGdaughters[1]=321;//K
237 fPDGdaughters[2]=211;//pi
238 fPDGdaughters[3]=0; //empty
239 fBranchName="Charm3Prong";
244 //_________________________________________________________________
245 Bool_t AliAnalysisTaskSESignificance::CheckConsistency(){
247 Bool_t result = kTRUE;
249 const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
250 //Float_t *vars = new Float_t[nvars];
251 Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
252 Bool_t *uppervars = fRDCuts->GetIsUpperCut();
253 TString *names = fRDCuts->GetVarNames();
255 for(Int_t i=0;i<fNPtBins;i++){
256 TString mdvname=Form("multiDimVectorPtBin%d",i);
259 for(Int_t ivar=0;ivar<nvars;ivar++){
260 if(varsforopt[ivar]){
261 Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
262 Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
264 AliFatal(Form("tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i));
267 Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
268 if(uppervars[ivar]&&lowermdv){
269 AliFatal(Form("%s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---please check your cuts \n ",names[ivar].Data(),i));
272 if(!uppervars[ivar]&&!lowermdv){
273 AliFatal(Form("%s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---please check your cuts \n",names[ivar].Data(),i));
282 //_________________________________________________________________
283 void AliAnalysisTaskSESignificance::SetBFeedDown(FeedDownEnum flagB){
284 if(fReadMC)fBFeedDown=flagB;
286 AliInfo("B feed down not allowed without MC info\n");
290 //_________________________________________________________________
291 void AliAnalysisTaskSESignificance::SetMassLimits(Float_t range, Int_t pdg){
293 Int_t abspdg=TMath::Abs(pdg);
294 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
295 fUpmasslimit = mass+range;
296 fLowmasslimit = mass-range;
298 //_________________________________________________________________
299 void AliAnalysisTaskSESignificance::SetMassLimits(Float_t lowlimit, Float_t uplimit){
302 fUpmasslimit = uplimit;
303 fLowmasslimit = lowlimit;
309 //________________________________________________________________________
310 void AliAnalysisTaskSESignificance::LocalInit()
314 if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
319 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
326 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
333 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
340 AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fRDCuts)));
347 AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fRDCuts)));
354 AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fRDCuts)));
364 TList *mdvList = new TList();
372 //________________________________________________________________________
373 void AliAnalysisTaskSESignificance::UserCreateOutputObjects()
375 // Create the output container
377 if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
379 // Several histograms are more conveniently managed in a TList
380 fOutput = new TList();
382 fOutput->SetName("OutputHistos");
384 //same number of steps in each multiDimVectorPtBin%d !
385 Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
386 cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
387 nHist=nHist*fNPtBins;
388 cout<<"Total = "<<nHist<<endl;
389 for(Int_t i=0;i<nHist;i++){
397 hisname.Form("hMass_%d",i);
398 signame.Form("hSig_%d",i);
399 bkgname.Form("hBkg_%d",i);
400 rflname.Form("hRfl_%d",i);
402 title.Form("Invariant mass;M[GeV/c^{2}];Entries");
404 fMassHist[i]=new TH1F(hisname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
405 fMassHist[i]->Sumw2();
406 fOutput->Add(fMassHist[i]);
409 fSigHist[i]=new TH1F(signame.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
410 fSigHist[i]->Sumw2();
411 fOutput->Add(fSigHist[i]);
413 fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
414 fBkgHist[i]->Sumw2();
415 fOutput->Add(fBkgHist[i]);
417 if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi){
418 fRflHist[i]=new TH1F(rflname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
419 fRflHist[i]->Sumw2();
420 fOutput->Add(fRflHist[i]);
425 fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",8,-0.5,7.5);
426 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
427 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvSelected (vtx)");
428 fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
429 fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
430 fHistNEvents->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
431 fHistNEvents->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
433 fHistNEvents->GetXaxis()->SetBinLabel(7,"MC Cand from c");
434 fHistNEvents->GetXaxis()->SetBinLabel(8,"MC Cand from b");
436 fHistNEvents->GetXaxis()->SetBinLabel(7,"N candidates");
438 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
439 fOutput->Add(fHistNEvents);
446 //________________________________________________________________________
447 void AliAnalysisTaskSESignificance::UserExec(Option_t */*option*/)
449 // Execute analysis for current event:
450 // heavy flavor candidates association to MC truth
452 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
453 if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
454 // Post the data already here
456 TClonesArray *arrayProng =0;
458 if(!aod && AODEvent() && IsStandardAOD()) {
459 // In case there is an AOD handler writing a standard AOD, use the AOD
460 // event in memory rather than the input (ESD) event.
461 aod = dynamic_cast<AliAODEvent*> (AODEvent());
462 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
463 // have to taken from the AOD event hold by the AliAODExtension
464 AliAODHandler* aodHandler = (AliAODHandler*)
465 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
466 if(aodHandler->GetExtensions()) {
468 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
469 AliAODEvent *aodFromExt = ext->GetAOD();
470 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
474 arrayProng=(TClonesArray*)aod->GetList()->FindObject(fBranchName.Data());
476 if(!aod || !arrayProng) {
477 AliError("AliAnalysisTaskSESignificance::UserExec:Branch not found!\n");
481 // fix for temporary bug in ESDfilter
482 // the AODs with null vertex pointer didn't pass the PhysSel
483 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
484 TClonesArray *arrayMC=0;
485 AliAODMCHeader *mcHeader=0;
490 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
492 AliError("AliAnalysisTaskSESignificance::UserExec:MC particles branch not found!\n");
497 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
499 AliError("AliAnalysisTaskSESignificance::UserExec:MC header branch not found!\n");
504 fHistNEvents->Fill(0); // count event
506 AliAODRecoDecayHF *d=0;
510 Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
511 Int_t nProng = arrayProng->GetEntriesFast();
512 if(fDebug>1) printf("Number of D2H: %d\n",nProng);
514 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
515 TString trigclass=aod->GetFiredTriggerClasses();
516 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(5);
518 if(fRDCuts->IsEventSelected(aod)) {
519 fHistNEvents->Fill(1);
521 if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
522 fHistNEvents->Fill(4);
526 for (Int_t iProng = 0; iProng < nProng; iProng++) {
527 fHistNEvents->Fill(6);
528 d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
530 AliAODRecoCascadeHF* DStarToD0pi = NULL;
531 AliAODRecoDecayHF2Prong* D0Particle = NULL;
532 fPDGDStarToD0pi[0] = 421; fPDGDStarToD0pi[1] = 211;
533 fPDGD0ToKpi[0] = 321; fPDGD0ToKpi[1] = 211;
535 Bool_t isSelBit=kTRUE;
538 isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
541 isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
544 isSelBit=d->HasSelectionBit(AliRDHFCuts::kDstarCuts);
547 isSelBit=d->HasSelectionBit(AliRDHFCuts::kDsCuts);
549 if(fDecChannel==5) isSelBit=d->HasSelectionBit(AliRDHFCuts::kLcCuts);
555 if(!isSelBit) continue;
557 if (fDecChannel==2) {
558 DStarToD0pi = (AliAODRecoCascadeHF*)arrayProng->At(iProng);
559 if (!DStarToD0pi->GetSecondaryVtx()) continue;
560 D0Particle = (AliAODRecoDecayHF2Prong*)DStarToD0pi->Get2Prong();
561 if (!D0Particle) continue;
564 Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(fPDGmother));
565 Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
567 if(fReadMC && fBFeedDown!=kBoth && isSelected){
568 Int_t labD = d->MatchToMC(fPDGmother,arrayMC,fNProngs,fPDGdaughters);
570 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
571 Int_t pdgGranma = CheckOrigin(partD, arrayMC);
572 Int_t abspdgGranma = TMath::Abs(pdgGranma);
573 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
575 AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
576 fHistNEvents->Fill(7);
577 if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
581 fHistNEvents->Fill(6);
582 if(fBFeedDown==kBeautyOnly)isSelected=kFALSE;
587 if(isSelected&&isFidAcc) {
588 fHistNEvents->Fill(2); // count selected with loosest cuts
589 if(fDebug>1) printf("+++++++Is Selected\n");
592 if(fDecChannel==3) SetPDGdaughterDstoKKpi();
593 fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
594 Int_t ptbin=fRDCuts->PtBin(d->Pt());
595 if(ptbin==-1) continue;
596 TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
597 AliMultiDimVector* muvec=(AliMultiDimVector*)fCutList->FindObject(mdvname.Data());
599 ULong64_t *addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
600 if(fDebug>1)printf("nvals = %d\n",nVals);
601 for(Int_t ivals=0;ivals<nVals;ivals++){
602 if(addresses[ivals]>=muvec->GetNTotCells()){
603 if (fDebug>1) printf("Overflow!!\n");
608 fHistNEvents->Fill(3);
610 //fill the histograms with the appropriate method
611 switch (fDecChannel){
613 FillDplus(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
616 FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
619 FillDstar(DStarToD0pi,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
623 FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,1);
627 FillD04p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
630 FillLambdac(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
638 if (fDecChannel==3 && isSelected&2){
639 SetPDGdaughterDstopiKK();
641 fRDCuts->GetCutVarsForOpt(d,fVars,fNVars,fPDGdaughters,aod);
643 addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
644 if(fDebug>1)printf("nvals = %d\n",nVals);
645 for(Int_t ivals=0;ivals<nVals;ivals++){
646 if(addresses[ivals]>=muvec->GetNTotCells()){
647 if (fDebug>1) printf("Overflow!!\n");
651 FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,0);
667 //***************************************************************************
669 // Methods used in the UserExec
672 //********************************************************************************************
674 //Methods to fill the istograms with MC information, one for each candidate
675 //NB: the implementation for each candidate is responsibility of the corresponding developer
677 void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
680 AliError("Candidate not selected\n");
684 if(fPartOrAndAntiPart*d->GetCharge()<0)return;
685 Int_t pdgdaughters[3] = {211,321,211};
686 Double_t mass=d->InvMass(3,(UInt_t*)pdgdaughters);
688 fMassHist[index]->Fill(mass);
693 lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
695 fSigHist[index]->Fill(mass);
697 fBkgHist[index]->Fill(mass);
703 void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
708 Int_t pdgdaughtersD0[2]={211,321};//pi,K
709 Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
711 masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
712 masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
714 if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(masses[0]);
715 if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(masses[1]);
725 Int_t pdgdaughters[2];
726 pdgdaughters[0]=211;//pi
727 pdgdaughters[1]=321;//K
729 matchtoMC = d->MatchToMC(fPDGmother,arrayMC,fNProngs,pdgdaughters);
731 Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
732 if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0){ //D0
734 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
735 Int_t pdgMC = dMC->GetPdgCode();
737 if(pdgMC==prongPdgPlus) fSigHist[index]->Fill(masses[0]);
738 else fRflHist[index]->Fill(masses[0]);
740 } else fBkgHist[index]->Fill(masses[0]);
743 if(isSel>=2 && fPartOrAndAntiPart<=0){ //D0bar
745 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
746 Int_t pdgMC = dMC->GetPdgCode();
748 if(pdgMC==prongPdgMinus) fSigHist[index]->Fill(masses[1]);
749 else fRflHist[index]->Fill(masses[1]);
750 } else fBkgHist[index]->Fill(masses[1]);
755 void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel){
758 AliInfo("Dstar selected\n");
760 Double_t mass = dstarD0pi->DeltaInvMass();
762 if((isSel>0) && TMath::Abs(fPartOrAndAntiPart)>=0) fMassHist[index]->Fill(mass);
765 Int_t matchtoMC = -1;
766 matchtoMC = dstarD0pi->MatchToMC(413,421,fPDGDStarToD0pi, fPDGD0ToKpi, arrayMC);
768 Int_t prongPdgDStarPlus=413;//,prongPdgDStarMinus=(-1)*413;
770 if ((isSel>1) && TMath::Abs(fPartOrAndAntiPart)>=0) {
773 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
774 Int_t pdgMC = dMC->GetPdgCode();
776 if (pdgMC==prongPdgDStarPlus) fSigHist[index]->Fill(mass);
778 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
779 mass = dstarD0pi->DeltaInvMass();
780 fRflHist[index]->Fill(mass);
781 dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
784 else fBkgHist[index]->Fill(mass);
790 void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay){
795 Int_t pdgDsKKpi[3]={321,321,211};//K,K,pi
796 Int_t pdgDspiKK[3]={211,321,321};//pi,K,K
798 masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgDsKKpi); //Ds
799 masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgDspiKK); //Dsbar
803 labDs = d->MatchToMC(431,arrayMC,3,pdgDsKKpi);
806 Int_t isKKpi=isSel&1;
807 Int_t ispiKK=isSel&2;
808 Int_t isPhiKKpi=isSel&4;
809 Int_t isPhipiKK=isSel&8;
810 Int_t isK0starKKpi=isSel&16;
811 Int_t isK0starpiKK=isSel&32;
814 if(fDsChannel==kPhi && (isPhiKKpi==0 && isPhipiKK==0)) return;
815 if(fDsChannel==kK0star && (isK0starKKpi==0 && isK0starpiKK==0)) return;
818 if(isKKpi && fPartOrAndAntiPart*d->GetCharge()>=0) {
819 if(fDsChannel==kPhi && isPhiKKpi==0) return;
820 if(fDsChannel==kK0star && isK0starKKpi==0) return;
822 fMassHist[index]->Fill(masses[0]);
826 Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
827 AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
828 Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
830 fSigHist[index]->Fill(masses[0]); //signal
832 fRflHist[index]->Fill(masses[0]); //Reflected signal
835 fBkgHist[index]->Fill(masses[0]); // Background
842 if(ispiKK && fPartOrAndAntiPart*d->GetCharge()>=0){
843 if(fDsChannel==kPhi && isPhipiKK==0) return;
844 if(fDsChannel==kK0star && isK0starpiKK==0) return;
846 fMassHist[index]->Fill(masses[1]);
850 Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
851 AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
852 Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
854 fSigHist[index]->Fill(masses[1]);
856 fRflHist[index]->Fill(masses[1]);
859 fBkgHist[index]->Fill(masses[1]);
866 void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
867 //D0->Kpipipi channel
868 AliInfo("D0 in 4 prongs channel not implemented\n");
872 void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* d,TClonesArray* arrayMC,Int_t index,Int_t isSel){
876 Int_t pdgdaughtersLc[3]={2212,321,211}; //p,K,pi
877 masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc);
878 Int_t pdgdaughtersLc2[3]={211,321,2212}; //pi,K,p
879 masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc2);
882 if(fPartOrAndAntiPart==0 || fPartOrAndAntiPart==d->GetCharge()) {
884 // isSel=1 : p K pi ; isSel=2 : pi K p ;
885 if(isSel==1 || isSel==3) fMassHist[index]->Fill(masses[0]);
886 if(isSel>=2) fMassHist[index]->Fill(masses[1]);
888 // Check the MC truth
895 Int_t matchtoMC = -1;
897 Int_t pPDG = 2212; // p
898 Int_t kPDG = 321; // K
899 Int_t piPDG = 211; // pi
900 Int_t absPdgMom = 4122;
901 matchtoMC = d->MatchToMC(absPdgMom,arrayMC,fNProngs,pdgdaughtersLc);
905 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
906 Int_t pdgMC = dMC->GetPdgCode();
907 if (TMath::Abs(pdgMC)!=absPdgMom) AliInfo("What's up, isn't it a lambdac ?!");
908 Int_t labDau0 = ((AliAODTrack*)d->GetDaughter(0))->GetLabel();
909 AliAODMCParticle* p0 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
910 Int_t pdgCode0 = TMath::Abs(p0->GetPdgCode());
911 Int_t labDau1 = ((AliAODTrack*)d->GetDaughter(1))->GetLabel();
912 AliAODMCParticle* p1 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau1);
913 Int_t pdgCode1 = TMath::Abs(p1->GetPdgCode());
914 Int_t labDau2 = ((AliAODTrack*)d->GetDaughter(2))->GetLabel();
915 AliAODMCParticle* p2 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau2);
916 Int_t pdgCode2 = TMath::Abs(p2->GetPdgCode());
918 // Fill in the histograms in case of p K pi decays
920 if(pdgCode0==pPDG && pdgCode1==kPDG && pdgCode2==piPDG){
921 fSigHist[index]->Fill(masses[0]);
923 fRflHist[index]->Fill(masses[0]);
926 // Fill in the histograms in case of pi K p decays
928 if(pdgCode0==piPDG && pdgCode1==kPDG && pdgCode2==pPDG){
929 fSigHist[index]->Fill(masses[1]);
931 fRflHist[index]->Fill(masses[1]);
935 if(ispKpi==1) fBkgHist[index]->Fill(masses[0]);
936 if(ispiKp==2) fBkgHist[index]->Fill(masses[1]);
945 //________________________________________________________________________
946 void AliAnalysisTaskSESignificance::Terminate(Option_t */*option*/)
948 // Terminate analysis
950 if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
953 fOutput = dynamic_cast<TList*> (GetOutputData(1));
955 printf("ERROR: fOutput not available\n");
959 fCutList = dynamic_cast<TList*> (GetOutputData(2));
961 printf("ERROR: fCutList not available\n");
965 AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
967 cout<<"multidimvec not found in TList"<<endl;
971 Int_t nHist=mdvtmp->GetNTotCells();
972 TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
974 for(Int_t i=0;i<nHist;i++){
981 hisname.Form("hMass_%d",i);
982 signame.Form("hSig_%d",i);
983 bkgname.Form("hBkg_%d",i);
984 rflname.Form("hRfl_%d",i);
986 fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
988 if (!drawn && fMassHist[i]->GetEntries() > 0){
990 fMassHist[i]->Draw();
995 fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
996 fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
997 if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
1007 //_________________________________________________________________________________________________
1008 Int_t AliAnalysisTaskSESignificance::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1011 // checking whether the very mother of the D0 is a charm or a bottom quark
1014 Int_t pdgGranma = 0;
1016 mother = mcPart->GetMother();
1020 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1021 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1022 if(!mcGranma) break;
1023 pdgGranma = mcGranma->GetPdgCode();
1024 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1025 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1026 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1029 mother = mcGranma->GetMother();
1033 //_________________________________________________________________________________________________