/************************************************************************** * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ ///////////////////////////////////////////////////////////// // // AliAnalysisTaskSE for the comparison of heavy flavor // decay candidates with the MC truth. // // Author: A.Dainese, andrea.dainese@lnl.infn.it ///////////////////////////////////////////////////////////// #include #include #include #include #include "AliAODEvent.h" #include "AliAODVertex.h" #include "AliAODTrack.h" #include "AliAODMCHeader.h" #include "AliAODMCParticle.h" #include "AliAODRecoDecayHF2Prong.h" #include "AliAnalysisTaskSE.h" #include "AliAnalysisTaskSECompareHF.h" ClassImp(AliAnalysisTaskSECompareHF) //________________________________________________________________________ AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(): AliAnalysisTaskSE(), fOutput(0), fNtupleD0Cmp(0), fHistMass(0) { // Default constructor SetD0toKpiCuts(); // Output slot #1 writes into a TList container DefineOutput(1,TList::Class()); //My private output } //________________________________________________________________________ AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name): AliAnalysisTaskSE(name), fOutput(0), fNtupleD0Cmp(0), fHistMass(0) { // Default constructor SetD0toKpiCuts(); // Output slot #1 writes into a TList container DefineOutput(1,TList::Class()); //My private output } //________________________________________________________________________ AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF() { // Destructor if (fOutput) { delete fOutput; fOutput = 0; } } //________________________________________________________________________ void AliAnalysisTaskSECompareHF::Init() { // Initialization if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n"); return; } //________________________________________________________________________ void AliAnalysisTaskSECompareHF::UserCreateOutputObjects() { // Create the output container // if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n"); // Several histograms are more conveniently managed in a TList fOutput = new TList(); fOutput->SetOwner(); fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965); fHistMass->Sumw2(); fHistMass->SetMinimum(0); fOutput->Add(fHistMass); fNtupleD0Cmp = new TNtuple("fNtupleD0Cmp","D0 comparison","pdg:VxRec:VxTrue:PtRec:PtTrue"); fOutput->Add(fNtupleD0Cmp); return; } //________________________________________________________________________ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/) { // Execute analysis for current event: // heavy flavor candidates association to MC truth AliAODEvent *aod = dynamic_cast (InputEvent()); // load D0->Kpi candidates TClonesArray *inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi"); if(!inputArrayD0toKpi) { printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n"); return; } // AOD primary vertex AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); //vtx1->Print(); // load MC particles TClonesArray *mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()); if(!mcArray) { printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n"); return; } // load MC header AliAODMCHeader *mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()); if(!mcHeader) { printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n"); return; } // loop over D0->Kpi candidates Int_t nInD0toKpi = inputArrayD0toKpi->GetEntriesFast(); printf("Number of D0->Kpi: %d\n",nInD0toKpi); Int_t lab0,lab1,labMother,labD0daugh0,labD0daugh1,pdgMother,pdgD0; for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) { AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayD0toKpi->UncheckedAt(iD0toKpi); labD0daugh0=-1; labD0daugh1=-1; Bool_t unsetvtx=kFALSE; if(!d->GetOwnPrimaryVtx()) { d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables unsetvtx=kTRUE; } Int_t okD0=0,okD0bar=0; if(d->SelectD0(fD0toKpiCuts,okD0,okD0bar)) { // get daughter AOD tracks AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0); AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1); lab0 = trk0->GetLabel(); lab1 = trk1->GetLabel(); AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(lab0); if(!part0) { printf("no MC particle\n"); continue; } while(part0->GetMother()>=0) { labMother=part0->GetMother(); part0 = (AliAODMCParticle*)mcArray->At(labMother); if(!part0) { printf("no MC mother particle\n"); break; } pdgMother = TMath::Abs(part0->GetPdgCode()); if(pdgMother==421) { labD0daugh0=labMother; break; } } AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(lab1); if(!part1) { printf("no MC particle\n"); continue; } while(part1->GetMother()>=0) { labMother=part1->GetMother(); part1 = (AliAODMCParticle*)mcArray->At(labMother); if(!part1) { printf("no MC mother particle\n"); break; } pdgMother = TMath::Abs(part1->GetPdgCode()); if(pdgMother==421) { labD0daugh1=labMother; break; } } // check if the candidate is signal if(labD0daugh0>=0 && labD0daugh1>=0 && labD0daugh0==labD0daugh1) { AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(labD0daugh0); pdgD0 = partD0->GetPdgCode(); Double_t invmass = (pdgD0==421 ? d->InvMassD0() : d->InvMassD0bar()); fHistMass->Fill(invmass); // Post the data already here PostData(1,fOutput); fNtupleD0Cmp->Fill(pdgD0,d->Xv(),partD0->Xv(),d->Pt(),partD0->Pt()); } } if(unsetvtx) d->UnsetOwnPrimaryVtx(); } // end loop on D0->Kpi return; } //________________________________________________________________________ void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/) { // Terminate analysis // if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n"); fOutput = dynamic_cast (GetOutputData(1)); if (!fOutput) { printf("ERROR: fOutput not available\n"); return; } fNtupleD0Cmp = dynamic_cast(fOutput->FindObject("fNtupleD0Cmp")); fHistMass = dynamic_cast(fOutput->FindObject("fHistMass")); return; } //________________________________________________________________________ void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(Double_t cut0,Double_t cut1, Double_t cut2,Double_t cut3,Double_t cut4, Double_t cut5,Double_t cut6, Double_t cut7,Double_t cut8) { // Set the cuts for D0 selection // cuts[0] = inv. mass half width [GeV] // cuts[1] = dca [cm] // cuts[2] = cosThetaStar // cuts[3] = pTK [GeV/c] // cuts[4] = pTPi [GeV/c] // cuts[5] = d0K [cm] upper limit! // cuts[6] = d0Pi [cm] upper limit! // cuts[7] = d0d0 [cm^2] // cuts[8] = cosThetaPoint fD0toKpiCuts[0] = cut0; fD0toKpiCuts[1] = cut1; fD0toKpiCuts[2] = cut2; fD0toKpiCuts[3] = cut3; fD0toKpiCuts[4] = cut4; fD0toKpiCuts[5] = cut5; fD0toKpiCuts[6] = cut6; fD0toKpiCuts[7] = cut7; fD0toKpiCuts[8] = cut8; return; } //________________________________________________________________________ void AliAnalysisTaskSECompareHF::SetD0toKpiCuts(const Double_t cuts[9]) { // Set the cuts for D0 selection // cuts[0] = inv. mass half width [GeV] // cuts[1] = dca [cm] // cuts[2] = cosThetaStar // cuts[3] = pTK [GeV/c] // cuts[4] = pTPi [GeV/c] // cuts[5] = d0K [cm] upper limit! // cuts[6] = d0Pi [cm] upper limit! // cuts[7] = d0d0 [cm^2] // cuts[8] = cosThetaPoint for(Int_t i=0; i<9; i++) fD0toKpiCuts[i] = cuts[i]; return; }