From: fprino Date: Mon, 25 Nov 2013 14:47:45 +0000 (+0000) Subject: Usage of multi-vertexer to tag pileup events in D meson analysis (Zaida) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=5d3cf93bf935914db1d1bae1e2f046e11456e725;hp=87a7c4b0525d56a8f88526629171ae3ce77c2385;p=u%2Fmrichter%2FAliRoot.git Usage of multi-vertexer to tag pileup events in D meson analysis (Zaida) --- diff --git a/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx b/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx index 37a755b51b8..3492e31bb9a 100644 --- a/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx +++ b/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.cxx @@ -1,1547 +1,1547 @@ -/************************************************************************** - * Copyright(c) 1998-2009, 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. * - **************************************************************************/ -// -// -// Base class for DStar - Hadron Correlations Analysis -// -//----------------------------------------------------------------------- -// -// -// Author S.Bjelogrlic -// Utrecht University -// sandro.bjelogrlic@cern.ch -// -//----------------------------------------------------------------------- - -/* $Id$ */ - -//#include -#include -#include -#include -#include "TROOT.h" - -#include "AliAnalysisTaskDStarCorrelations.h" -#include "AliRDHFCutsDStartoKpipi.h" -#include "AliHFAssociatedTrackCuts.h" -#include "AliAODRecoDecay.h" -#include "AliAODRecoCascadeHF.h" -#include "AliAODRecoDecayHF2Prong.h" -#include "AliAODPidHF.h" -#include "AliVParticle.h" -#include "AliAnalysisManager.h" -#include "AliAODInputHandler.h" -#include "AliAODHandler.h" -#include "AliESDtrack.h" -#include "AliAODMCParticle.h" -#include "AliNormalizationCounter.h" -#include "AliReducedParticle.h" -#include "AliHFCorrelator.h" -#include "AliAODMCHeader.h" -#include "AliEventPoolManager.h" -#include "AliVertexingHFUtils.h" - -using std::cout; -using std::endl; - - -ClassImp(AliAnalysisTaskDStarCorrelations) - - -//__________________________________________________________________________ -AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() : -AliAnalysisTaskSE(), -fhandler(0x0), -fmcArray(0x0), -fCounter(0x0), -fCorrelator(0x0), -fselect(0), -fmontecarlo(kFALSE), -fmixing(kFALSE), -fFullmode(kFALSE), -fSystem(pp), -fEfficiencyVariable(kNone), -fReco(kTRUE), -fUseEfficiencyCorrection(kFALSE), -fUseDmesonEfficiencyCorrection(kFALSE), -fUseCentrality(kFALSE), -fUseHadronicChannelAtKineLevel(kFALSE), -fPhiBins(32), -fEvents(0), -fDebugLevel(0), -fDisplacement(0), -fDim(0), -fNofPtBins(0), -fMaxEtaDStar(0.9), -fDMesonSigmas(0), - -fD0Window(0), - -fOutput(0x0), -fOutputMC(0x0), -fCuts(0), -fCuts2(0), -fUtils(0), -fTracklets(0), -fDeffMapvsPt(0), -fDeffMapvsPtvsMult(0), -fDeffMapvsPtvsEta(0) -{ - SetDim(); - // default constructor -} - -//__________________________________________________________________________ -AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) : -AliAnalysisTaskSE(name), - -fhandler(0x0), -fmcArray(0x0), -fCounter(0x0), -fCorrelator(0x0), -fselect(0), -fmontecarlo(kFALSE), -fmixing(kFALSE), -fFullmode(mode), -fSystem(syst), -fEfficiencyVariable(kNone), -fReco(kTRUE), -fUseEfficiencyCorrection(kFALSE), -fUseDmesonEfficiencyCorrection(kFALSE), -fUseCentrality(kFALSE), -fUseHadronicChannelAtKineLevel(kFALSE), -fPhiBins(32), -fEvents(0), -fDebugLevel(0), -fDisplacement(0), -fDim(0), -fNofPtBins(0), -fMaxEtaDStar(0.9), -fDMesonSigmas(0), -fD0Window(0), - -fOutput(0x0), -fOutputMC(0x0), -fCuts(0), -fCuts2(AsscCuts), -fUtils(0), -fTracklets(0), -fDeffMapvsPt(0), -fDeffMapvsPtvsMult(0), -fDeffMapvsPtvsEta(0) -{ - Info("AliAnalysisTaskDStarCorrelations","Calling Constructor"); - - SetDim(); - if(fSystem == AA) fUseCentrality = kTRUE; else fUseCentrality = kFALSE; - - fCuts=cuts; - fNofPtBins= fCuts->GetNPtBins(); - //cout << "Enlarging the DZero window " << endl; - EnlargeDZeroMassWindow(); - // cout << "Done" << endl; - - - DefineInput(0, TChain::Class()); - - DefineOutput(1,TList::Class()); // histos from data and MC - DefineOutput(2,TList::Class()); // histos from MC - DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts - DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts - DefineOutput(5,AliNormalizationCounter::Class()); // normalization -} - -//__________________________________________________________________________ - -AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() { - // - // destructor - // - - Info("AliAnalysisTaskDStarCorrelations","Calling Destructor"); - - if(fhandler) {delete fhandler; fhandler = 0;} - //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;} - if(fmcArray) {delete fmcArray; fmcArray = 0;} - if(fCounter) {delete fCounter; fCounter = 0;} - if(fCorrelator) {delete fCorrelator; fCorrelator = 0;} - if(fOutput) {delete fOutput; fOutput = 0;} - if(fOutputMC) {delete fOutputMC; fOutputMC = 0;} - if(fCuts) {delete fCuts; fCuts = 0;} - if(fCuts2) {delete fCuts2; fCuts2=0;} - if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;} - if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;} - if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;} - -} - -//___________________________________________________________ -void AliAnalysisTaskDStarCorrelations::Init(){ - // - // Initialization - // - if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n"); - - AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts); - - - - // Post the D* cuts - PostData(3,copyfCuts); - - // Post the hadron cuts - PostData(4,fCuts2); - - return; -} - - -//_________________________________________________ -void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){ - Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); - - //slot #1 - //OpenFile(0); - fOutput = new TList(); - fOutput->SetOwner(); - - fOutputMC = new TList(); - fOutputMC->SetOwner(); - - // define histograms - DefineHistoForAnalysis(); - DefineThNSparseForAnalysis(); - - - fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName())); - fCounter->Init(); - - Double_t Pi = TMath::Pi(); - fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fUseCentrality); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb - fCorrelator->SetDeltaPhiInterval( -0.5*Pi, 1.5*Pi); // set correct phi interval - //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval - fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On - fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros - fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut - fCorrelator->SetUseMC(fmontecarlo); - fCorrelator->SetUseReco(fReco); - // fCorrelator->SetKinkRemoval(kTRUE); - Bool_t pooldef = fCorrelator->DefineEventPool(); - - if(!pooldef) AliInfo("Warning:: Event pool not defined properly"); - - fUtils = new AliAnalysisUtils(); - - - - PostData(1,fOutput); // set the outputs - PostData(2,fOutputMC); // set the outputs - PostData(5,fCounter); // set the outputs -} - -//_________________________________________________ -void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){ - - - if(fDebugLevel){ - - if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl; - if(!fReco) std::cout << "USING MC TRUTH" << std::endl; - std::cout << " " << std::endl; - std::cout << "=================================================================================" << std::endl; - if(!fmixing){ - if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl; - if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl; - if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl; - } - if(fmixing){ - if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl; - if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl; - if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl; - } - - }// end if debug - - - - if (!fInputEvent) { - Error("UserExec","NO EVENT FOUND!"); - return; - } - - AliAODEvent* aodEvent = dynamic_cast(fInputEvent); - if(!aodEvent){ - AliError("AOD event not found!"); - return; - } - - fTracklets = aodEvent->GetTracklets(); - - fEvents++; // event counter - ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0); - - fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo); - - // load MC array - fmcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); - if(fmontecarlo && !fmcArray){ - AliError("Array of MC particles not found"); - return; - } - - - - - // ********************************************** START EVENT SELECTION **************************************************** - - Bool_t isEvSel=fCuts->IsEventSelected(aodEvent); - - if(!isEvSel) { - - if(fCuts->IsEventRejectedDueToPileupSPD()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2); - if(fCuts->IsEventRejectedDueToCentrality()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3); - if(fCuts->IsEventRejectedDueToNotRecoVertex()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4); - if(fCuts->IsEventRejectedDueToVertexContributors()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5); - if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6); - if(fCuts->IsEventRejectedDueToTrigger()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7); - if(fCuts->IsEventRejectedDuePhysicsSelection()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8); - - return; - } - - // added event selection for pA - - if(fSystem == pA){ - - if(fUtils->IsFirstEventInChunk(aodEvent)) { - AliInfo("Rejecting the event - first in the chunk"); - ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9); - return; - } - if(!fUtils->IsVertexSelected2013pA(aodEvent)) { - ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10); - AliInfo("Rejecting the event - bad vertex"); - return; - } - } - // ******************************** END event selections ************************************************** - - AliCentrality *centralityObj = 0; - Double_t MultipOrCent = -1; - - if(fUseCentrality){ - /* if(fSystem == AA ){ */ centralityObj = aodEvent->GetHeader()->GetCentralityP(); - MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M"); - //AliInfo(Form("Centrality is %f", MultipOrCent)); - } - - if(!fUseCentrality) MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.); - - - fCorrelator->SetAODEvent(aodEvent); // set the event to be processed - - ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1); - - Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event is in pool settings - if(!correlatorON) { - AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event"); - return; - } - - if(fmontecarlo) fCorrelator->SetMCArray(fmcArray); - - // check the event type - // load MC header - if(fmontecarlo){ - AliAODMCHeader *mcHeader = dynamic_cast(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName())); - if (fmontecarlo && !mcHeader) { - AliError("Could not find MC Header in AOD"); - return; - } - - Bool_t isMCeventgood = kFALSE; - - - Int_t eventType = mcHeader->GetEventType(); - Int_t NMCevents = fCuts2->GetNofMCEventType(); - - for(Int_t k=0; kGetMCEventType(); - - if(eventType == MCEventType[k]) isMCeventgood= kTRUE; - ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType); - } - - if(NMCevents && !isMCeventgood){ - if(fDebugLevel) std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl; - return; - } - - } // end if montecarlo - - - // checks on vertex and multiplicity of the event - AliAODVertex *vtx = aodEvent->GetPrimaryVertex(); - Double_t zVtxPosition = vtx->GetZ(); // zvertex - - - if(fFullmode) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition); - - - - // D* reconstruction - TClonesArray *arrayDStartoD0pi=0; - if(!aodEvent && AODEvent() && IsStandardAOD()) { - // In case there is an AOD handler writing a standard AOD, use the AOD - // event in memory rather than the input (ESD) event. - aodEvent = dynamic_cast (AODEvent()); - // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) - // have to taken from the AOD event hold by the AliAODExtension - AliAODHandler* aodHandler = (AliAODHandler*) - ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); - if(aodHandler->GetExtensions()) { - AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); - AliAODEvent *aodFromExt = ext->GetAOD(); - arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar"); - } - } else { - arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar"); - } - - if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return; - - - // get the poolbin - - Int_t poolbin = fCuts2->GetPoolBin(MultipOrCent, zVtxPosition); - - - // initialize variables you will need for the D* - - Double_t ptDStar;// - Double_t phiDStar;// - Double_t etaDStar;// - Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag; - Double_t invMassDZero; - Double_t deltainvMDStar; - - - Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass(); - Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass(); - - - //MC tagging for DStar - //D* and D0 prongs needed to MatchToMC method - Int_t pdgDgDStartoD0pi[2]={421,211}; - Int_t pdgDgD0toKpi[2]={321,211}; - - Bool_t isDStarCand = kFALSE; - Bool_t isDfromB = kFALSE; - Bool_t isEventMixingFilledPeak = kFALSE; - Bool_t isEventMixingFilledSB = kFALSE; - Bool_t EventHasDStarCandidate = kFALSE; - Bool_t EventHasDZeroSideBandCandidate = kFALSE; - Bool_t EventHasDStarSideBandCandidate = kFALSE; - //loop on D* candidates - - Int_t looponDCands = 0; - if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast(); - if(!fReco) looponDCands = fmcArray->GetEntriesFast(); - - Int_t nOfDStarCandidates = 0; - Int_t nOfSBCandidates = 0; - - Double_t DmesonEfficiency = 1.; - Double_t DmesonWeight = 1.; - Double_t efficiencyvariable = -999; - - - - for (Int_t iDStartoD0pi = 0; iDStartoD0piAt(iDStartoD0pi); - if(!dstarD0pi->GetSecondaryVtx()) continue; - theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong(); - if (!theD0particle) continue; - - - // track quality cuts - Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks - // region of interest + topological cuts + PID - Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected - - //apply track selections - if(!isTkSelected) continue; - if(!isSelected) continue; - if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue; - - - ptDStar = dstarD0pi->Pt(); - phiDStar = dstarD0pi->Phi(); - etaDStar = dstarD0pi->Eta(); - if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; - if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr) efficiencyvariable = MultipOrCent; - if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar; - if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar(); - if(fEfficiencyVariable == kNone) efficiencyvariable = 0; - - // get the D meson efficiency - DmesonEfficiency = fCuts2->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable); - - - - if(fUseDmesonEfficiencyCorrection){ - if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency; - else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI - if(ptDStar>2.) DmesonWeight = 0.5; // at high pt a zero value in the efficiency can come only from stat fluctutations in MC -> 0.5 is an arbitrary asymptotic value - else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low - } - } - else DmesonWeight = 1.; - - // continue; - - Int_t mcLabelDStar = -999; - if(fmontecarlo){ - // find associated MC particle for D* ->D0toKpi - mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/); - if(mcLabelDStar>=0) isDStarMCtag = kTRUE; - if(!isDStarMCtag) continue; - AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar); - //check if DStar from B - Int_t labelMother = MCDStar->GetMother(); - AliAODMCParticle * mother = dynamic_cast(fmcArray->At(labelMother)); - if(!mother) continue; - Int_t motherPDG =TMath::Abs(mother->PdgCode()); - if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE; - } - - - phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); - - // set the phi of the D meson in the correct range - - Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt()); - - - Double_t dmDStarWindow = 0.0019/3;// 0.0019 = 3 sigma - // Double_t mD0Window=0.074/3; - - Double_t mD0Window= fD0Window[ptbin]/3; - //cout << "Check with new window " << fD0Window[ptbin]/3 << endl; - - - - invMassDZero = dstarD0pi->InvMassD0(); - if(!fmixing && fFullmode) ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero); - - deltainvMDStar = dstarD0pi->DeltaInvMass(); - - //good D0 candidates - if (TMath::Abs(invMassDZero-mPDGD0)FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight); - // good D*? - if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); - if(!fmixing) ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar); - isInPeak = kTRUE; - EventHasDStarCandidate = kTRUE; - nOfDStarCandidates++; - } // end Good D* - - // D* sideband? - if((deltainvMDStar-(mPDGDstar-mPDGD0)>fDMesonSigmas[2]*dmDStarWindow) && (deltainvMDStar-(mPDGDstar-mPDGD0)fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight); - if(!fmixing && fFullmode)((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero,DmesonWeight); - - if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); - isInDZeroSideBand = kTRUE; - EventHasDZeroSideBandCandidate = kTRUE; - nOfSBCandidates++; - if(!fmixing) ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar); - } - - }//end if sidebands - - - - - if(!isInPeak && !isInDStarSideBand && !isInDZeroSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME - - - // check properties of the events containing the D* - - - - isDStarCand = kTRUE; - - // charge of the daughter od the - daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge(); - - - trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID(); - trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID(); - trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID(); - - // end here the reco - - - }// end of if for applied reconstruction to D* - - Int_t DStarLabel = -1; - - if(!fReco){ // use pure MC information - - // get the DStar Particle - DStarMC = dynamic_cast(fmcArray->At(iDStartoD0pi)); - if (!DStarMC) { - AliWarning("Careful: DStar MC Particle not found in tree, skipping"); - continue; - } - DStarLabel = DStarMC->GetLabel(); - if(DStarLabel>0)isDStarMCtag = kTRUE; - - Int_t PDG =TMath::Abs(DStarMC->PdgCode()); - if(PDG !=413) continue; // skip if it is not a DStar - // check fiducial acceptance - if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue; - - //check if DStar from B - Int_t labelMother = DStarMC->GetMother(); - AliAODMCParticle * mother = dynamic_cast(fmcArray->At(labelMother)); - if(!mother) continue; - Int_t motherPDG =TMath::Abs(mother->PdgCode()); - if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE; - - Bool_t isDZero = kFALSE; - Bool_t isSoftPi = kFALSE; - - if(fUseHadronicChannelAtKineLevel){ - //check decay channel on MC ************************************************ - Int_t NDaugh = DStarMC->GetNDaughters(); - if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs - - for(Int_t i=0; iGetDaughter(i); - AliAODMCParticle* mcDaughter = dynamic_cast(fmcArray->At(daugh_label)); - if(!mcDaughter) continue; - Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode()); - if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl; - - if(daugh_pdg == 421) {isDZero = kTRUE; - Int_t NDaughD0 = mcDaughter->GetNDaughters(); - if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs - trackiddaugh0 = mcDaughter->GetDaughter(0); - trackiddaugh1 = mcDaughter->GetDaughter(1); - Bool_t isKaon = kFALSE; - Bool_t isPion = kFALSE; - - for(Int_t k=0;kGetDaughter(k); - AliAODMCParticle* mcGrandDaughter = dynamic_cast(fmcArray->At(labelD0daugh)); - if(!mcGrandDaughter) continue; - Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode()); - if(granddaugh_pdg==321) isKaon = kTRUE; - if(granddaugh_pdg==211) isPion = kTRUE; - } - if(!isKaon || !isKaon) continue; // skip if not correct decay channel of D0 - } - - if(daugh_pdg == 211) { - isSoftPi = kTRUE; - daughtercharge = mcDaughter->Charge(); - trackidsoftPi = daugh_label;} - } - if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel - } // end check decay channel - - ptDStar = DStarMC->Pt(); - phiDStar = DStarMC->Phi(); - etaDStar = DStarMC->Eta(); - - if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; - - } // end use pure MC information - - - // getting the number of triggers in the MCtag D* case - if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar); - if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar); - if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar); - - - fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters - fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge); - - - // ************************************************ CORRELATION ANALYSIS STARTS HERE - - - Bool_t execPool = fCorrelator->ProcessEventPool(); - - - if(fmixing && !execPool) { - AliInfo("Mixed event analysis: pool is not ready"); - if(!isEventMixingFilledPeak && isInPeak) { - ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(1); - isEventMixingFilledPeak = kTRUE; - } - if (!isEventMixingFilledSB && isInDZeroSideBand) { - ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(3); - isEventMixingFilledSB=kTRUE; - } - - continue; - } - - // check event topology - if(fmixing&&execPool){ - // pool is ready - run checks on bins filling - if(!isEventMixingFilledPeak && isInPeak) { - ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(0); - if(fFullmode) EventMixingChecks(aodEvent); - isEventMixingFilledPeak = kTRUE; - } - - if(!isEventMixingFilledSB && isInDZeroSideBand) { - ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(2); - isEventMixingFilledSB=kTRUE; - } - } - - Int_t NofEventsinPool = 1; - if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool(); - - - for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one - - Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix); - if(!analyzetracks) { - AliInfo("AliHFCorrelator::Cannot process the track array"); - continue; - } - - //initialization of variables for correlations with leading particles - Double_t DeltaPhiLeading = -999.; - Double_t DeltaEtaLeading = -999.; - - - Int_t NofTracks = fCorrelator->GetNofTracks(); - - - if(isInPeak && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(NofTracks); - if(isInDZeroSideBand && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(NofTracks); - - - - Double_t arraytofill[5]; - Double_t MCarraytofill[7]; - - - Double_t weight; - - for(Int_t iTrack = 0; iTrackCorrelate(iTrack); - if(!runcorrelation) continue; - - Double_t DeltaPhi = fCorrelator->GetDeltaPhi(); - Double_t DeltaEta = fCorrelator->GetDeltaEta(); - - AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle(); - if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;} - - Double_t ptHad = hadron->Pt(); - Double_t phiHad = hadron->Phi(); - Double_t etaHad = hadron->Eta(); - Int_t label = hadron->GetLabel(); - Int_t trackid = hadron->GetID(); - Double_t efficiency = hadron->GetWeight(); - - weight = 1; - if(fUseEfficiencyCorrection && efficiency){ - weight = DmesonWeight * (1./efficiency); - } - - phiHad = fCorrelator->SetCorrectPhiRange(phiHad); - - - if(fFullmode) ((TH2F*)fOutput->FindObject("WeightChecks"))->Fill(ptHad,efficiency); - - arraytofill[0] = DeltaPhi; - arraytofill[1] = DeltaEta; - arraytofill[2] = ptDStar; - arraytofill[3] = ptHad; - arraytofill[4] = poolbin; - - - MCarraytofill[0] = DeltaPhi; - MCarraytofill[1] = DeltaEta; - MCarraytofill[2] = ptDStar; - MCarraytofill[3] = ptHad; - MCarraytofill[4] = poolbin; - - - if(fmontecarlo){ - if(label<0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(0.,NofTracks); - if(label>=0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(1.,NofTracks); - if(label<0) continue; // skip track with wrong label - } - - Bool_t isDdaughter = kFALSE; - // skip the D daughters in the correlation - if(!fmixing && fReco){ - if(trackid == trackiddaugh0) continue; - if(trackid == trackiddaugh1) continue; - if(trackid == trackidsoftPi) continue; - } - - if(!fmixing && !fReco){ - AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label); - if(!part) continue; - if(IsDDaughter(DStarMC, part)) continue; - cout << "Skipping DStar daugheter " << endl; - } - if(!fmixing && !fReco && fmontecarlo){ // skip D* Daughetrs if it is Pure MCDStar - Int_t hadronlabel = label; - for(Int_t k=0; k<4;k++){ // go back 4 generations and check the mothers - if(DStarLabel<0){ break;} - if(hadronlabel<0) { break;} - AliAODMCParticle* mcParticle = dynamic_cast(fmcArray->At(hadronlabel)); - if(!mcParticle) {AliInfo("NO MC PARTICLE"); break;} - hadronlabel = mcParticle->GetMother(); - if(hadronlabel == DStarLabel) isDdaughter = kTRUE; - } - - if(isDdaughter && fDebugLevel){ - std::cout << "It is the D* daughter with label " << label << std::endl; - std::cout << "Daughter 0 label = " << trackiddaugh0 << std::endl; - std::cout << "Daughter 1 label = " << trackiddaugh1 << std::endl; - std::cout << "Soft pi label = " << trackidsoftPi << std::endl; - } - - if(isDdaughter) continue; // skip if track is from DStar - } - - // ================ FILL CORRELATION HISTOGRAMS =============================== - - // monte carlo case (mc tagged D*) - if((fmontecarlo && isDStarMCtag) || (fmontecarlo && !fReco)){ // check correlations of MC tagged DStars in MonteCarlo - - Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0) - - MCarraytofill[5] = 0; - if(PartSource[0]) MCarraytofill[5] = 1; - if(PartSource[1]) MCarraytofill[5] = 2; - if(PartSource[2]&&PartSource[0]) MCarraytofill[5] = 3; - if(PartSource[2]&&PartSource[1]) MCarraytofill[5] = 4; - if(PartSource[3]) MCarraytofill[5] = 5; - if(!isDfromB) MCarraytofill[6] = 0; - if(isDfromB) MCarraytofill[6] = 1; - if(!fReco && TMath::Abs(etaHad)>0.8) { - delete [] PartSource; - continue; // makes sure you study the correlation on MC truth only if particles are in acceptance - } - ((THnSparseF*)fOutputMC->FindObject("MCDStarCorrelationsDStarHadron"))->Fill(MCarraytofill); - - delete[] PartSource; - } - - // Good DStar canidates - if(isInPeak) { - - if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance - if(fselect==1) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarHadron"))->Fill(arraytofill,weight); - if(fselect==2) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKaon"))->Fill(arraytofill,weight); - if(fselect==3) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKZero"))->Fill(arraytofill,weight); - - ((TH3F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad,MultipOrCent); - if(fFullmode)((TH1D*)fOutput->FindObject("TracksInPeakSpectra"))->Fill(ptHad); - - } - - // Sidebands from D0 candidate - if(isInDZeroSideBand) { - - if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance - if(fselect==1) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight); - if(fselect==2) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight); - if(fselect==3) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight); - - if(fFullmode) ((TH1D*)fOutput->FindObject("TracksInSBSpectra"))->Fill(ptHad); - - } - - // Sidebands from D* candidate - if(isInDStarSideBand) { - - if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance - if(fselect==1 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight); - if(fselect==2 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight); - if(fselect==3 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight); - - } - - - } // end loop on track candidates - - - - // fill the leading particle histograms - - if(isInPeak && fFullmode) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading); - if(isInDZeroSideBand && fFullmode) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading); - - } // end loop on events in the pool - - }// end loop on D* candidates - - - - // check events with D* or SB canidates - if(fFullmode && EventHasDStarCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDStar"))->Fill(MultipOrCent,zVtxPosition); - if(fFullmode && EventHasDZeroSideBandCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDZeroSB"))->Fill(MultipOrCent,zVtxPosition); - - if(fFullmode && EventHasDStarCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDStarSB"))->Fill(MultipOrCent,zVtxPosition); - - - if(fFullmode) ((TH2F*)fOutput->FindObject("DStarCandidates"))->Fill(nOfDStarCandidates,MultipOrCent); - if(fFullmode) ((TH2F*)fOutput->FindObject("SBCandidates"))->Fill(nOfSBCandidates,MultipOrCent); - - // update event pool - Bool_t updated = fCorrelator->PoolUpdate(); - - // if(updated) EventMixingChecks(aodEvent); - if(!updated) AliInfo("Pool was not updated"); - - -} //end the exec - -//________________________________________ terminate ___________________________ -void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*) -{ - // The Terminate() function is the last function to be called during - // a query. It always runs on the client, it can be used to present - // the results graphically or save the results to file. - - AliAnalysisTaskSE::Terminate(); - - fOutput = dynamic_cast (GetOutputData(1)); - if (!fOutput) { - printf("ERROR: fOutput not available\n"); - return; - } - - return; -} -//_____________________________________________________ -Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const { - - //Daughter removal in MCKine case - Bool_t isDaughter = kFALSE; - Int_t labelD0 = d->GetLabel(); - - Int_t mother = track->GetMother(); - - //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!) - while (mother > 0){ - AliAODMCParticle* mcMoth = dynamic_cast(fmcArray->At(mother)); //it's the mother of the track! - if (mcMoth){ - if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE; - mother = mcMoth->GetMother(); //goes back by one - } else{ - AliError("Failed casting the mother particle!"); - break; - } - } - - return isDaughter; -} - -//_____________________________________________________ -void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){ - - //cout << "DEFINING THNSPARSES "<< endl; - - Double_t Pi = TMath::Pi(); - Int_t nbinscorr = fPhiBins; - Double_t lowcorrbin = -0.5*Pi; - Double_t upcorrbin = 1.5*Pi; - // define the THnSparseF - - //sparse bins - - //1 delta_phi - //2 delta_eta - //3 D* pt - //4 multiplicity - //5 track pt - //6 zVtx position - - Int_t nbinsPool = (fCuts2->GetNZvtxPoolBins())*(fCuts2->GetNCentPoolBins()); - - - Int_t nbinsSparse[5]= {nbinscorr, 32,30,250,nbinsPool}; - Double_t binLowLimitSparse[5]={lowcorrbin,-1.6, 0, 0,-0.5}; - Double_t binUpLimitSparse[5]= {upcorrbin, 1.6,30,25,nbinsPool-0.5}; - - Int_t MCnbinsSparse[7]= {nbinscorr, 32,30,250,nbinsPool,10,2}; - Double_t MCbinLowLimitSparse[7]={lowcorrbin,-1.6, 0, 0,-0.5,-0.5,-0.5}; // - Double_t MCbinUpLimitSparse[7]= {upcorrbin, 1.6,30,25,nbinsPool-0.5,9.5,1.5}; - - TString sparsename = "CorrelationsDStar"; - if(fselect==1) sparsename += "Hadron"; - if(fselect==2) sparsename += "Kaon"; - if(fselect==3) sparsename += "KZero"; - - TString D0Bkgsparsename = "DZeroBkg"; - D0Bkgsparsename += sparsename; - - TString DStarBkgsparsename = "DStarBkg"; - DStarBkgsparsename += sparsename; - - TString MCSparseName = "MCDStar"; - MCSparseName += sparsename; - // signal correlations - THnSparseF * Correlations = new THnSparseF(sparsename.Data(),"Correlations for signal",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse); - - // bkg correlations from D0 sidebands - THnSparseF * DZeroBkgCorrelations = new THnSparseF(D0Bkgsparsename.Data(),"Bkg Correlations estimated with D0 sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse); - - // bkg correlations from D* sidebands - THnSparseF * DStarBkgCorrelations = new THnSparseF(DStarBkgsparsename.Data(),"Bkg Correlations estimated with D* sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse); - - - THnSparseF * MCCorrelations = new THnSparseF(MCSparseName.Data(),"MC Correlations",7,MCnbinsSparse,MCbinLowLimitSparse,MCbinUpLimitSparse); - - MCCorrelations->GetAxis(5)->SetBinLabel(1," All "); - MCCorrelations->GetAxis(5)->SetBinLabel(2," from hadron Heavy flavour"); - MCCorrelations->GetAxis(5)->SetBinLabel(3," from c->D"); - MCCorrelations->GetAxis(5)->SetBinLabel(4," from b->D"); - MCCorrelations->GetAxis(5)->SetBinLabel(5," from b->B"); - MCCorrelations->GetAxis(5)->SetBinLabel(6," from quark Heavy flavour"); - MCCorrelations->GetAxis(5)->SetBinLabel(7," from c"); - MCCorrelations->GetAxis(5)->SetBinLabel(8," from b"); - - MCCorrelations->GetAxis(6)->SetBinLabel(1," if D* from c"); - MCCorrelations->GetAxis(6)->SetBinLabel(2," if D* from b"); - - Correlations->Sumw2(); - DZeroBkgCorrelations->Sumw2(); - DStarBkgCorrelations->Sumw2(); - - fOutput->Add(Correlations); - fOutput->Add(DZeroBkgCorrelations); - if(fFullmode) fOutput->Add(DStarBkgCorrelations); - if(fmontecarlo) fOutputMC->Add(MCCorrelations); - - -} -//__________________________________________________________________________________________________ -void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){ - - Double_t Pi = TMath::Pi(); - Int_t nbinscorr = fPhiBins; - Double_t lowcorrbin = -0.5*Pi ; // shift the bin by half the width so that at 0 is it the bin center - Double_t upcorrbin = 1.5*Pi ; - - // ========================= histograms for both Data and MonteCarlo - - - TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5); - NofEvents->GetXaxis()->SetBinLabel(1," All events"); - NofEvents->GetXaxis()->SetBinLabel(2," Selected events"); - NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup"); - NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality"); - NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx"); - NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr."); - NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc."); - NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger"); - NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel"); - NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk"); - NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx"); - fOutput->Add(NofEvents); - - - - - TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5); - if(!fmixing && fFullmode) fOutput->Add(D0InvMass); - - TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5); - if(!fmixing && fFullmode) fOutput->Add(D0InvMassinSB); - - //TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2); - //if(!fmixing) fOutput->Add(DeltaInvMass); - TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution; D* p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2); - if(!fmixing) fOutput->Add(DeltaInvMass); - - TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution; SB p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2); - if(!fmixing) fOutput->Add(bkgDeltaInvMass); - - DeltaInvMass->Sumw2(); - bkgDeltaInvMass->Sumw2(); - - TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50); - if(!fmixing) fOutput->Add(RecoPtDStar); - - TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50); - if(!fmixing) fOutput->Add(RecoPtBkg); - - TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50); - if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm); - - TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50); - if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty); - - TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50); - if(!fmixing) fOutput->Add(MCtagPtDStar); - - TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25); - if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectra); - - TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25); - if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectraifHF); - - TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","N of associated tracks per D trigger; Nof tracks; Entries",500,-0.5,499.5); - if(fFullmode) fOutput->Add(NofTracksInPeak); - - TH1D * NofTracksInSB = new TH1D("NofTracksInSB","N of associated tracks per SideBand trigger; Nof tracks; Entries",500,-0.5,499.5); - if(fFullmode) fOutput->Add(NofTracksInSB); - - TH1D * TracksInPeakSpectra = new TH1D("TracksInPeakSpectra","Pt Spectra tracks with D trigger; p_{T} GeV/c; Entries",500,-0.5,49.5); - if(fFullmode)fOutput->Add(TracksInPeakSpectra); - - TH1D * TracksInSBSpectra = new TH1D("TracksInSBSpectra","Pt Spectra tracks with SideBand trigger; p_{T} GeV/c; Entries",500,-0.5,49.5); - if(fFullmode)fOutput->Add(TracksInSBSpectra); - - - //TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5); - //if(fmixing) fOutput->Add(EventMixingCheck); - - - TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); - if(fFullmode)fOutput->Add(EventPropsCheck); - - TH2F * EventPropsCheckifDStar = new TH2F("EventPropsCheckifDStar","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); - if(fFullmode)fOutput->Add(EventPropsCheckifDStar); - - TH2F * EventPropsCheckifDZeroSB = new TH2F("EventPropsCheckifDZeroSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); - if(fFullmode)fOutput->Add(EventPropsCheckifDZeroSB); - - TH2F * EventPropsCheckifDStarSB = new TH2F("EventPropsCheckifDStarSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); - if(fFullmode)fOutput->Add(EventPropsCheckifDStarSB); - - - TH2F * WeightChecks = new TH2F("WeightChecks","Checks on efficiency correction",300,0,30,100,0.005,1.005); - if(fFullmode)fOutput->Add(WeightChecks); - - - - TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9); - fOutput->Add(PhiEtaTrigger); - - TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9); - fOutput->Add(PhiEtaSideBand); - - TH3F * PhiEtaPart = new TH3F("PhiEtaPart","#phi distribution of the associated particle; #phi; #eta; multiplicity",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9,100,0,1000); - fOutput->Add(PhiEtaPart); - - TH2F * DStarCandidates = new TH2F("DStarCandidates","# of D* candidates per event vs multiplicity",6,-0.5,5.5,50,0,500); - if(fFullmode)fOutput->Add(DStarCandidates); - - TH2F * SBCandidates = new TH2F("SBCandidates","# of SB candidates per event vs multiplicity",6,-0.5,5.5,50,0,500); - if(fFullmode)fOutput->Add(SBCandidates); - - - //correlations histograms - TString histoname1 = "DPhiDStar"; - if(fselect==1) histoname1 += "Hadron"; - if(fselect==2) histoname1 += "Kaon"; - if(fselect==3) histoname1 += "KZero"; - - /* - TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - //side band background histograms - TString histoname2 = "bkg"; - histoname2 += histoname1; - TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - - fOutput->Add(DPhiDStar); - - if(fselect==3){fOutput->Add(DPhiDStarKZero1);} - - fOutput->Add(bkgDPhiDStar); - - if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);} - - */ - // leading particle - TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - if(fFullmode)fOutput->Add(leadingcand); - if(fFullmode)fOutput->Add(leadingsidebands); - - // ========================= histos for analysis on MC only - - TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5); - if(fmontecarlo) fOutputMC->Add(EventTypeMC); - - TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5); - MCSources->GetXaxis()->SetBinLabel(1," All "); - MCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour"); - MCSources->GetXaxis()->SetBinLabel(3," from c->D"); - MCSources->GetXaxis()->SetBinLabel(4," from b->D"); - MCSources->GetXaxis()->SetBinLabel(5," from b->B"); - MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour"); - MCSources->GetXaxis()->SetBinLabel(7," from c"); - MCSources->GetXaxis()->SetBinLabel(8," from b"); - - if(fmontecarlo) fOutputMC->Add(MCSources); - - // leading particle from mc source - TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5); - LeadingMCSources->GetXaxis()->SetBinLabel(1," All "); - LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour"); - LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D"); - LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D"); - LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B"); - LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour"); - LeadingMCSources->GetXaxis()->SetBinLabel(7," from c"); - LeadingMCSources->GetXaxis()->SetBinLabel(8," from b"); - - if(fmontecarlo && fFullmode) fOutputMC->Add(LeadingMCSources); - - // all hadrons - TString histoname3 = "MCTag"; - histoname3 += histoname1; - TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString histoname44 = "CharmDOrigin"; - histoname44 += histoname1; - histoname44 += "MC"; - - TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - - TString histoname54 = "BeautyDOrigin"; - histoname54 += histoname1; - histoname54 += "MC"; - TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString histoname55 = "BeautyBOrigin"; - histoname55 += histoname1; - histoname55 += "MC"; - TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString histoname4 = "CharmQuarkOrigin"; - histoname4 += histoname1; - histoname4 += "MC"; - TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString histoname5 = "BeautyQuarkOrigin"; - histoname5 += histoname1; - histoname5 += "MC"; - TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString histoname6 = "NonHFOrigin"; - histoname6 += histoname1; - histoname6 += "MC"; - TH3D * NonHFOriginDPhiDStar = new TH3D(histoname6.Data(),histoname6.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - if(fmontecarlo && fFullmode){ - - fOutputMC->Add(MCTagDPhiDStar); - fOutputMC->Add(CharmDOriginDPhiDStar); - fOutputMC->Add(BeautyDOriginDPhiDStar); - fOutputMC->Add(BeautyBOriginDPhiDStar); - fOutputMC->Add(CharmQuarkOriginDPhiDStar); - fOutputMC->Add(BeautyQuarkOriginDPhiDStar); - fOutputMC->Add(NonHFOriginDPhiDStar); - - } - - // ========================= histos for analysis on MC - // all leading hadron - TString Leadinghistoname3 = "LeadingMCTag"; - Leadinghistoname3 += histoname1; - TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString Leadinghistoname44 = "LeadingCharmDOrigin"; - Leadinghistoname44 += histoname1; - Leadinghistoname44 += "MC"; - - TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - - TString Leadinghistoname54 = "LeadingBeautyDOrigin"; - Leadinghistoname54 += histoname1; - Leadinghistoname54 += "MC"; - TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString Leadinghistoname55 = "LeadingBeautyBOrigin"; - Leadinghistoname55 += histoname1; - Leadinghistoname55 += "MC"; - TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString Leadinghistoname4 = "LeadingCharmQuarkOrigin"; - Leadinghistoname4 += histoname1; - Leadinghistoname4 += "MC"; - TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin"; - Leadinghistoname5 += histoname1; - Leadinghistoname5 += "MC"; - TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); - - - - - if(fmontecarlo && fFullmode){ - - fOutputMC->Add(LeadingMCTagDPhiDStar); - fOutputMC->Add(LeadingCharmDOriginDPhiDStar); - fOutputMC->Add(LeadingBeautyDOriginDPhiDStar); - fOutputMC->Add(LeadingBeautyBOriginDPhiDStar); - fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar); - fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar); - - } - - TH3F * MCPhiEtaPart = new TH3F("MCPhiEtaPart","#phi distribution of the associated particle",nbinscorr,lowcorrbin,upcorrbin,50,-2.5,2.5,6,-0.5,6.5); - MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles"); - MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark"); - MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark"); - MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c"); - MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b"); - MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b"); - if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart); - - TH2D * TrackLabels = new TH2D("TrackLabels","NofEvents;track label; multiplicity",2,-0.5,1.5,500,-0.5,499.5); - if(fmontecarlo && fFullmode) fOutputMC->Add(TrackLabels); - - // ============================= EVENT MIXING CHECKS ====================================== - - Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool(); - Int_t MinNofTracks = fCuts2->GetMinNTracksInPool(); - Int_t NofCentBins = fCuts2->GetNCentPoolBins(); - Double_t * CentBins = fCuts2->GetCentPoolBins(); - Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins(); - Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins(); - - - - Int_t k =0; - - if(fSystem == AA) k = 100; // PbPb centrality - if(fSystem == pp || fSystem == pA) k = NofCentBins; // pp multiplicity - - - //Double_t minvalue = CentBins[0]; - //Double_t maxvalue = CentBins[NofCentBins+1]; - //Double_t Zminvalue = ZVrtxBins[0]; - //Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1]; - - Double_t minvalue, maxvalue; - Double_t Zminvalue, Zmaxvalue; - - Zminvalue = -15.; - Zmaxvalue = 15; - if(fSystem == AA) {minvalue = 0; maxvalue = 100;} // PbPb - if(fSystem == pp || fSystem == pA) {minvalue = 0; maxvalue = 500;} // multilpicity - - //Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents}; - // Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents}; - // Double_t * events = Nevents; - Double_t eventsv[] ={0,1000000}; - //Double_t * events = new Double_t[2]; - // events[0] = 0; -// events[1] = 1000000; - Double_t *events = eventsv; - Int_t Nevents = 1000000; - // TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,Nevents,events); - - TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,Nevents,events[0],events[1]); - - EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity "); - EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]"); - EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin"); - if(fmixing && fFullmode) fOutput->Add(EventsPerPoolBin); - - Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks; - //Int_t Diff = MaxNofTracks-MinNofTracks; - - //Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks}; - // Double_t * trackN = Ntracks; - - TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,MaxNofTracks,0,MaxNofTracks); - NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity "); - NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]"); - NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin"); - - if(fmixing && fFullmode) fOutput->Add(NofTracksPerPoolBin); - - TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Calls per pool bin",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins); - NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity "); - NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]"); - if(fmixing && fFullmode) fOutput->Add(NofPoolBinCalls); - - - - TH2D * EventProps = new TH2D("EventProps","Event properties",100,minvalue,maxvalue,100,Zminvalue,Zmaxvalue); - EventProps->GetXaxis()->SetTitle("Centrality/multiplicity "); - EventProps->GetYaxis()->SetTitle("Z vertex [cm]"); - if(fmixing && fFullmode) fOutput->Add(EventProps); - - TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5); - CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready"); - CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready"); - CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready"); - CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready"); - - if(fmixing) fOutput->Add(CheckPoolReadiness); - - -} - -//__________________________________________________________________________________________________ -void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){ - - - //Float_t* ptbins = fCuts->GetPtBinLimits(); - if(fD0Window) delete fD0Window; - fD0Window = new Float_t[fNofPtBins]; - - AliInfo("Enlarging the D0 mass windows from cut object\n"); - Int_t nvars = fCuts->GetNVars(); - - if(nvars<1){ - AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!"); - return; - } - Float_t** rdcutsvalmine=new Float_t*[nvars]; - for(Int_t iv=0;ivGetCutValue(0,j); - rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j); - cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl; - } - else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j); - - // set same windows - //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j); - } - } - - fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine); - - AliInfo("\n New windows set\n"); - fCuts->PrintAll(); - - - for(Int_t iv=0;ivGetNTracks(); - MultipOrCent = multiplicity; // convert from Int_t to Double_t - } - if(fSystem == AA){ // PbPb - - centralityObj = AOD->GetHeader()->GetCentralityP(); - MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M"); - AliInfo(Form("Centrality is %f", MultipOrCent)); - } - - AliAODVertex *vtx = AOD->GetPrimaryVertex(); - Double_t zvertex = vtx->GetZ(); // zvertex - - - - - AliEventPool * pool = fCorrelator->GetPool(); - - - - - ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool - ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties - - ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool - ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool -} - - - - - +/************************************************************************** + * Copyright(c) 1998-2009, 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. * + **************************************************************************/ +// +// +// Base class for DStar - Hadron Correlations Analysis +// +//----------------------------------------------------------------------- +// +// +// Author S.Bjelogrlic +// Utrecht University +// sandro.bjelogrlic@cern.ch +// +//----------------------------------------------------------------------- + +/* $Id$ */ + +//#include +#include +#include +#include +#include "TROOT.h" + +#include "AliAnalysisTaskDStarCorrelations.h" +#include "AliRDHFCutsDStartoKpipi.h" +#include "AliHFAssociatedTrackCuts.h" +#include "AliAODRecoDecay.h" +#include "AliAODRecoCascadeHF.h" +#include "AliAODRecoDecayHF2Prong.h" +#include "AliAODPidHF.h" +#include "AliVParticle.h" +#include "AliAnalysisManager.h" +#include "AliAODInputHandler.h" +#include "AliAODHandler.h" +#include "AliESDtrack.h" +#include "AliAODMCParticle.h" +#include "AliNormalizationCounter.h" +#include "AliReducedParticle.h" +#include "AliHFCorrelator.h" +#include "AliAODMCHeader.h" +#include "AliEventPoolManager.h" +#include "AliVertexingHFUtils.h" + +using std::cout; +using std::endl; + + +ClassImp(AliAnalysisTaskDStarCorrelations) + + +//__________________________________________________________________________ +AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() : +AliAnalysisTaskSE(), +fhandler(0x0), +fmcArray(0x0), +fCounter(0x0), +fCorrelator(0x0), +fselect(0), +fmontecarlo(kFALSE), +fmixing(kFALSE), +fFullmode(kFALSE), +fSystem(pp), +fEfficiencyVariable(kNone), +fReco(kTRUE), +fUseEfficiencyCorrection(kFALSE), +fUseDmesonEfficiencyCorrection(kFALSE), +fUseCentrality(kFALSE), +fUseHadronicChannelAtKineLevel(kFALSE), +fPhiBins(32), +fEvents(0), +fDebugLevel(0), +fDisplacement(0), +fDim(0), +fNofPtBins(0), +fMaxEtaDStar(0.9), +fDMesonSigmas(0), + +fD0Window(0), + +fOutput(0x0), +fOutputMC(0x0), +fCuts(0), +fCuts2(0), +fUtils(0), +fTracklets(0), +fDeffMapvsPt(0), +fDeffMapvsPtvsMult(0), +fDeffMapvsPtvsEta(0) +{ + SetDim(); + // default constructor +} + +//__________________________________________________________________________ +AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) : +AliAnalysisTaskSE(name), + +fhandler(0x0), +fmcArray(0x0), +fCounter(0x0), +fCorrelator(0x0), +fselect(0), +fmontecarlo(kFALSE), +fmixing(kFALSE), +fFullmode(mode), +fSystem(syst), +fEfficiencyVariable(kNone), +fReco(kTRUE), +fUseEfficiencyCorrection(kFALSE), +fUseDmesonEfficiencyCorrection(kFALSE), +fUseCentrality(kFALSE), +fUseHadronicChannelAtKineLevel(kFALSE), +fPhiBins(32), +fEvents(0), +fDebugLevel(0), +fDisplacement(0), +fDim(0), +fNofPtBins(0), +fMaxEtaDStar(0.9), +fDMesonSigmas(0), +fD0Window(0), + +fOutput(0x0), +fOutputMC(0x0), +fCuts(0), +fCuts2(AsscCuts), +fUtils(0), +fTracklets(0), +fDeffMapvsPt(0), +fDeffMapvsPtvsMult(0), +fDeffMapvsPtvsEta(0) +{ + Info("AliAnalysisTaskDStarCorrelations","Calling Constructor"); + + SetDim(); + if(fSystem == AA) fUseCentrality = kTRUE; else fUseCentrality = kFALSE; + + fCuts=cuts; + fNofPtBins= fCuts->GetNPtBins(); + //cout << "Enlarging the DZero window " << endl; + EnlargeDZeroMassWindow(); + // cout << "Done" << endl; + + + DefineInput(0, TChain::Class()); + + DefineOutput(1,TList::Class()); // histos from data and MC + DefineOutput(2,TList::Class()); // histos from MC + DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts + DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts + DefineOutput(5,AliNormalizationCounter::Class()); // normalization +} + +//__________________________________________________________________________ + +AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() { + // + // destructor + // + + Info("AliAnalysisTaskDStarCorrelations","Calling Destructor"); + + if(fhandler) {delete fhandler; fhandler = 0;} + //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;} + if(fmcArray) {delete fmcArray; fmcArray = 0;} + if(fCounter) {delete fCounter; fCounter = 0;} + if(fCorrelator) {delete fCorrelator; fCorrelator = 0;} + if(fOutput) {delete fOutput; fOutput = 0;} + if(fOutputMC) {delete fOutputMC; fOutputMC = 0;} + if(fCuts) {delete fCuts; fCuts = 0;} + if(fCuts2) {delete fCuts2; fCuts2=0;} + if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;} + if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;} + if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;} + +} + +//___________________________________________________________ +void AliAnalysisTaskDStarCorrelations::Init(){ + // + // Initialization + // + if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n"); + + AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts); + + + + // Post the D* cuts + PostData(3,copyfCuts); + + // Post the hadron cuts + PostData(4,fCuts2); + + return; +} + + +//_________________________________________________ +void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){ + Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); + + //slot #1 + //OpenFile(0); + fOutput = new TList(); + fOutput->SetOwner(); + + fOutputMC = new TList(); + fOutputMC->SetOwner(); + + // define histograms + DefineHistoForAnalysis(); + DefineThNSparseForAnalysis(); + + + fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName())); + fCounter->Init(); + + Double_t Pi = TMath::Pi(); + fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fUseCentrality); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb + fCorrelator->SetDeltaPhiInterval( -0.5*Pi, 1.5*Pi); // set correct phi interval + //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval + fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On + fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros + fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut + fCorrelator->SetUseMC(fmontecarlo); + fCorrelator->SetUseReco(fReco); + // fCorrelator->SetKinkRemoval(kTRUE); + Bool_t pooldef = fCorrelator->DefineEventPool(); + + if(!pooldef) AliInfo("Warning:: Event pool not defined properly"); + + fUtils = new AliAnalysisUtils(); + + + + PostData(1,fOutput); // set the outputs + PostData(2,fOutputMC); // set the outputs + PostData(5,fCounter); // set the outputs +} + +//_________________________________________________ +void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){ + + + if(fDebugLevel){ + + if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl; + if(!fReco) std::cout << "USING MC TRUTH" << std::endl; + std::cout << " " << std::endl; + std::cout << "=================================================================================" << std::endl; + if(!fmixing){ + if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl; + if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl; + if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl; + } + if(fmixing){ + if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl; + if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl; + if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl; + } + + }// end if debug + + + + if (!fInputEvent) { + Error("UserExec","NO EVENT FOUND!"); + return; + } + + AliAODEvent* aodEvent = dynamic_cast(fInputEvent); + if(!aodEvent){ + AliError("AOD event not found!"); + return; + } + + fTracklets = aodEvent->GetTracklets(); + + fEvents++; // event counter + ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0); + + fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo); + + // load MC array + fmcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); + if(fmontecarlo && !fmcArray){ + AliError("Array of MC particles not found"); + return; + } + + + + + // ********************************************** START EVENT SELECTION **************************************************** + + Bool_t isEvSel=fCuts->IsEventSelected(aodEvent); + + if(!isEvSel) { + + if(fCuts->IsEventRejectedDueToPileup()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2); + if(fCuts->IsEventRejectedDueToCentrality()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3); + if(fCuts->IsEventRejectedDueToNotRecoVertex()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4); + if(fCuts->IsEventRejectedDueToVertexContributors()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5); + if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6); + if(fCuts->IsEventRejectedDueToTrigger()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7); + if(fCuts->IsEventRejectedDuePhysicsSelection()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8); + + return; + } + + // added event selection for pA + + if(fSystem == pA){ + + if(fUtils->IsFirstEventInChunk(aodEvent)) { + AliInfo("Rejecting the event - first in the chunk"); + ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9); + return; + } + if(!fUtils->IsVertexSelected2013pA(aodEvent)) { + ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10); + AliInfo("Rejecting the event - bad vertex"); + return; + } + } + // ******************************** END event selections ************************************************** + + AliCentrality *centralityObj = 0; + Double_t MultipOrCent = -1; + + if(fUseCentrality){ + /* if(fSystem == AA ){ */ centralityObj = aodEvent->GetHeader()->GetCentralityP(); + MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M"); + //AliInfo(Form("Centrality is %f", MultipOrCent)); + } + + if(!fUseCentrality) MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.); + + + fCorrelator->SetAODEvent(aodEvent); // set the event to be processed + + ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1); + + Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event is in pool settings + if(!correlatorON) { + AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event"); + return; + } + + if(fmontecarlo) fCorrelator->SetMCArray(fmcArray); + + // check the event type + // load MC header + if(fmontecarlo){ + AliAODMCHeader *mcHeader = dynamic_cast(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName())); + if (fmontecarlo && !mcHeader) { + AliError("Could not find MC Header in AOD"); + return; + } + + Bool_t isMCeventgood = kFALSE; + + + Int_t eventType = mcHeader->GetEventType(); + Int_t NMCevents = fCuts2->GetNofMCEventType(); + + for(Int_t k=0; kGetMCEventType(); + + if(eventType == MCEventType[k]) isMCeventgood= kTRUE; + ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType); + } + + if(NMCevents && !isMCeventgood){ + if(fDebugLevel) std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl; + return; + } + + } // end if montecarlo + + + // checks on vertex and multiplicity of the event + AliAODVertex *vtx = aodEvent->GetPrimaryVertex(); + Double_t zVtxPosition = vtx->GetZ(); // zvertex + + + if(fFullmode) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition); + + + + // D* reconstruction + TClonesArray *arrayDStartoD0pi=0; + if(!aodEvent && AODEvent() && IsStandardAOD()) { + // In case there is an AOD handler writing a standard AOD, use the AOD + // event in memory rather than the input (ESD) event. + aodEvent = dynamic_cast (AODEvent()); + // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) + // have to taken from the AOD event hold by the AliAODExtension + AliAODHandler* aodHandler = (AliAODHandler*) + ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); + if(aodHandler->GetExtensions()) { + AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); + AliAODEvent *aodFromExt = ext->GetAOD(); + arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar"); + } + } else { + arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar"); + } + + if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return; + + + // get the poolbin + + Int_t poolbin = fCuts2->GetPoolBin(MultipOrCent, zVtxPosition); + + + // initialize variables you will need for the D* + + Double_t ptDStar;// + Double_t phiDStar;// + Double_t etaDStar;// + Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag; + Double_t invMassDZero; + Double_t deltainvMDStar; + + + Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass(); + Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass(); + + + //MC tagging for DStar + //D* and D0 prongs needed to MatchToMC method + Int_t pdgDgDStartoD0pi[2]={421,211}; + Int_t pdgDgD0toKpi[2]={321,211}; + + Bool_t isDStarCand = kFALSE; + Bool_t isDfromB = kFALSE; + Bool_t isEventMixingFilledPeak = kFALSE; + Bool_t isEventMixingFilledSB = kFALSE; + Bool_t EventHasDStarCandidate = kFALSE; + Bool_t EventHasDZeroSideBandCandidate = kFALSE; + Bool_t EventHasDStarSideBandCandidate = kFALSE; + //loop on D* candidates + + Int_t looponDCands = 0; + if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast(); + if(!fReco) looponDCands = fmcArray->GetEntriesFast(); + + Int_t nOfDStarCandidates = 0; + Int_t nOfSBCandidates = 0; + + Double_t DmesonEfficiency = 1.; + Double_t DmesonWeight = 1.; + Double_t efficiencyvariable = -999; + + + + for (Int_t iDStartoD0pi = 0; iDStartoD0piAt(iDStartoD0pi); + if(!dstarD0pi->GetSecondaryVtx()) continue; + theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong(); + if (!theD0particle) continue; + + + // track quality cuts + Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks + // region of interest + topological cuts + PID + Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected + + //apply track selections + if(!isTkSelected) continue; + if(!isSelected) continue; + if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue; + + + ptDStar = dstarD0pi->Pt(); + phiDStar = dstarD0pi->Phi(); + etaDStar = dstarD0pi->Eta(); + if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; + if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr) efficiencyvariable = MultipOrCent; + if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar; + if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar(); + if(fEfficiencyVariable == kNone) efficiencyvariable = 0; + + // get the D meson efficiency + DmesonEfficiency = fCuts2->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable); + + + + if(fUseDmesonEfficiencyCorrection){ + if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency; + else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI + if(ptDStar>2.) DmesonWeight = 0.5; // at high pt a zero value in the efficiency can come only from stat fluctutations in MC -> 0.5 is an arbitrary asymptotic value + else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low + } + } + else DmesonWeight = 1.; + + // continue; + + Int_t mcLabelDStar = -999; + if(fmontecarlo){ + // find associated MC particle for D* ->D0toKpi + mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/); + if(mcLabelDStar>=0) isDStarMCtag = kTRUE; + if(!isDStarMCtag) continue; + AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar); + //check if DStar from B + Int_t labelMother = MCDStar->GetMother(); + AliAODMCParticle * mother = dynamic_cast(fmcArray->At(labelMother)); + if(!mother) continue; + Int_t motherPDG =TMath::Abs(mother->PdgCode()); + if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE; + } + + + phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); + + // set the phi of the D meson in the correct range + + Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt()); + + + Double_t dmDStarWindow = 0.0019/3;// 0.0019 = 3 sigma + // Double_t mD0Window=0.074/3; + + Double_t mD0Window= fD0Window[ptbin]/3; + //cout << "Check with new window " << fD0Window[ptbin]/3 << endl; + + + + invMassDZero = dstarD0pi->InvMassD0(); + if(!fmixing && fFullmode) ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero); + + deltainvMDStar = dstarD0pi->DeltaInvMass(); + + //good D0 candidates + if (TMath::Abs(invMassDZero-mPDGD0)FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight); + // good D*? + if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); + if(!fmixing) ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar); + isInPeak = kTRUE; + EventHasDStarCandidate = kTRUE; + nOfDStarCandidates++; + } // end Good D* + + // D* sideband? + if((deltainvMDStar-(mPDGDstar-mPDGD0)>fDMesonSigmas[2]*dmDStarWindow) && (deltainvMDStar-(mPDGDstar-mPDGD0)fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight); + if(!fmixing && fFullmode)((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero,DmesonWeight); + + if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); + isInDZeroSideBand = kTRUE; + EventHasDZeroSideBandCandidate = kTRUE; + nOfSBCandidates++; + if(!fmixing) ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar); + } + + }//end if sidebands + + + + + if(!isInPeak && !isInDStarSideBand && !isInDZeroSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME + + + // check properties of the events containing the D* + + + + isDStarCand = kTRUE; + + // charge of the daughter od the + daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge(); + + + trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID(); + trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID(); + trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID(); + + // end here the reco + + + }// end of if for applied reconstruction to D* + + Int_t DStarLabel = -1; + + if(!fReco){ // use pure MC information + + // get the DStar Particle + DStarMC = dynamic_cast(fmcArray->At(iDStartoD0pi)); + if (!DStarMC) { + AliWarning("Careful: DStar MC Particle not found in tree, skipping"); + continue; + } + DStarLabel = DStarMC->GetLabel(); + if(DStarLabel>0)isDStarMCtag = kTRUE; + + Int_t PDG =TMath::Abs(DStarMC->PdgCode()); + if(PDG !=413) continue; // skip if it is not a DStar + // check fiducial acceptance + if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue; + + //check if DStar from B + Int_t labelMother = DStarMC->GetMother(); + AliAODMCParticle * mother = dynamic_cast(fmcArray->At(labelMother)); + if(!mother) continue; + Int_t motherPDG =TMath::Abs(mother->PdgCode()); + if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE; + + Bool_t isDZero = kFALSE; + Bool_t isSoftPi = kFALSE; + + if(fUseHadronicChannelAtKineLevel){ + //check decay channel on MC ************************************************ + Int_t NDaugh = DStarMC->GetNDaughters(); + if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs + + for(Int_t i=0; iGetDaughter(i); + AliAODMCParticle* mcDaughter = dynamic_cast(fmcArray->At(daugh_label)); + if(!mcDaughter) continue; + Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode()); + if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl; + + if(daugh_pdg == 421) {isDZero = kTRUE; + Int_t NDaughD0 = mcDaughter->GetNDaughters(); + if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs + trackiddaugh0 = mcDaughter->GetDaughter(0); + trackiddaugh1 = mcDaughter->GetDaughter(1); + Bool_t isKaon = kFALSE; + Bool_t isPion = kFALSE; + + for(Int_t k=0;kGetDaughter(k); + AliAODMCParticle* mcGrandDaughter = dynamic_cast(fmcArray->At(labelD0daugh)); + if(!mcGrandDaughter) continue; + Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode()); + if(granddaugh_pdg==321) isKaon = kTRUE; + if(granddaugh_pdg==211) isPion = kTRUE; + } + if(!isKaon || !isKaon) continue; // skip if not correct decay channel of D0 + } + + if(daugh_pdg == 211) { + isSoftPi = kTRUE; + daughtercharge = mcDaughter->Charge(); + trackidsoftPi = daugh_label;} + } + if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel + } // end check decay channel + + ptDStar = DStarMC->Pt(); + phiDStar = DStarMC->Phi(); + etaDStar = DStarMC->Eta(); + + if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; + + } // end use pure MC information + + + // getting the number of triggers in the MCtag D* case + if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar); + if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar); + if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar); + + + fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters + fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge); + + + // ************************************************ CORRELATION ANALYSIS STARTS HERE + + + Bool_t execPool = fCorrelator->ProcessEventPool(); + + + if(fmixing && !execPool) { + AliInfo("Mixed event analysis: pool is not ready"); + if(!isEventMixingFilledPeak && isInPeak) { + ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(1); + isEventMixingFilledPeak = kTRUE; + } + if (!isEventMixingFilledSB && isInDZeroSideBand) { + ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(3); + isEventMixingFilledSB=kTRUE; + } + + continue; + } + + // check event topology + if(fmixing&&execPool){ + // pool is ready - run checks on bins filling + if(!isEventMixingFilledPeak && isInPeak) { + ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(0); + if(fFullmode) EventMixingChecks(aodEvent); + isEventMixingFilledPeak = kTRUE; + } + + if(!isEventMixingFilledSB && isInDZeroSideBand) { + ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(2); + isEventMixingFilledSB=kTRUE; + } + } + + Int_t NofEventsinPool = 1; + if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool(); + + + for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one + + Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix); + if(!analyzetracks) { + AliInfo("AliHFCorrelator::Cannot process the track array"); + continue; + } + + //initialization of variables for correlations with leading particles + Double_t DeltaPhiLeading = -999.; + Double_t DeltaEtaLeading = -999.; + + + Int_t NofTracks = fCorrelator->GetNofTracks(); + + + if(isInPeak && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(NofTracks); + if(isInDZeroSideBand && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(NofTracks); + + + + Double_t arraytofill[5]; + Double_t MCarraytofill[7]; + + + Double_t weight; + + for(Int_t iTrack = 0; iTrackCorrelate(iTrack); + if(!runcorrelation) continue; + + Double_t DeltaPhi = fCorrelator->GetDeltaPhi(); + Double_t DeltaEta = fCorrelator->GetDeltaEta(); + + AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle(); + if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;} + + Double_t ptHad = hadron->Pt(); + Double_t phiHad = hadron->Phi(); + Double_t etaHad = hadron->Eta(); + Int_t label = hadron->GetLabel(); + Int_t trackid = hadron->GetID(); + Double_t efficiency = hadron->GetWeight(); + + weight = 1; + if(fUseEfficiencyCorrection && efficiency){ + weight = DmesonWeight * (1./efficiency); + } + + phiHad = fCorrelator->SetCorrectPhiRange(phiHad); + + + if(fFullmode) ((TH2F*)fOutput->FindObject("WeightChecks"))->Fill(ptHad,efficiency); + + arraytofill[0] = DeltaPhi; + arraytofill[1] = DeltaEta; + arraytofill[2] = ptDStar; + arraytofill[3] = ptHad; + arraytofill[4] = poolbin; + + + MCarraytofill[0] = DeltaPhi; + MCarraytofill[1] = DeltaEta; + MCarraytofill[2] = ptDStar; + MCarraytofill[3] = ptHad; + MCarraytofill[4] = poolbin; + + + if(fmontecarlo){ + if(label<0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(0.,NofTracks); + if(label>=0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(1.,NofTracks); + if(label<0) continue; // skip track with wrong label + } + + Bool_t isDdaughter = kFALSE; + // skip the D daughters in the correlation + if(!fmixing && fReco){ + if(trackid == trackiddaugh0) continue; + if(trackid == trackiddaugh1) continue; + if(trackid == trackidsoftPi) continue; + } + + if(!fmixing && !fReco){ + AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label); + if(!part) continue; + if(IsDDaughter(DStarMC, part)) continue; + cout << "Skipping DStar daugheter " << endl; + } + if(!fmixing && !fReco && fmontecarlo){ // skip D* Daughetrs if it is Pure MCDStar + Int_t hadronlabel = label; + for(Int_t k=0; k<4;k++){ // go back 4 generations and check the mothers + if(DStarLabel<0){ break;} + if(hadronlabel<0) { break;} + AliAODMCParticle* mcParticle = dynamic_cast(fmcArray->At(hadronlabel)); + if(!mcParticle) {AliInfo("NO MC PARTICLE"); break;} + hadronlabel = mcParticle->GetMother(); + if(hadronlabel == DStarLabel) isDdaughter = kTRUE; + } + + if(isDdaughter && fDebugLevel){ + std::cout << "It is the D* daughter with label " << label << std::endl; + std::cout << "Daughter 0 label = " << trackiddaugh0 << std::endl; + std::cout << "Daughter 1 label = " << trackiddaugh1 << std::endl; + std::cout << "Soft pi label = " << trackidsoftPi << std::endl; + } + + if(isDdaughter) continue; // skip if track is from DStar + } + + // ================ FILL CORRELATION HISTOGRAMS =============================== + + // monte carlo case (mc tagged D*) + if((fmontecarlo && isDStarMCtag) || (fmontecarlo && !fReco)){ // check correlations of MC tagged DStars in MonteCarlo + + Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0) + + MCarraytofill[5] = 0; + if(PartSource[0]) MCarraytofill[5] = 1; + if(PartSource[1]) MCarraytofill[5] = 2; + if(PartSource[2]&&PartSource[0]) MCarraytofill[5] = 3; + if(PartSource[2]&&PartSource[1]) MCarraytofill[5] = 4; + if(PartSource[3]) MCarraytofill[5] = 5; + if(!isDfromB) MCarraytofill[6] = 0; + if(isDfromB) MCarraytofill[6] = 1; + if(!fReco && TMath::Abs(etaHad)>0.8) { + delete [] PartSource; + continue; // makes sure you study the correlation on MC truth only if particles are in acceptance + } + ((THnSparseF*)fOutputMC->FindObject("MCDStarCorrelationsDStarHadron"))->Fill(MCarraytofill); + + delete[] PartSource; + } + + // Good DStar canidates + if(isInPeak) { + + if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance + if(fselect==1) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarHadron"))->Fill(arraytofill,weight); + if(fselect==2) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKaon"))->Fill(arraytofill,weight); + if(fselect==3) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKZero"))->Fill(arraytofill,weight); + + ((TH3F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad,MultipOrCent); + if(fFullmode)((TH1D*)fOutput->FindObject("TracksInPeakSpectra"))->Fill(ptHad); + + } + + // Sidebands from D0 candidate + if(isInDZeroSideBand) { + + if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance + if(fselect==1) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight); + if(fselect==2) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight); + if(fselect==3) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight); + + if(fFullmode) ((TH1D*)fOutput->FindObject("TracksInSBSpectra"))->Fill(ptHad); + + } + + // Sidebands from D* candidate + if(isInDStarSideBand) { + + if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance + if(fselect==1 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight); + if(fselect==2 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight); + if(fselect==3 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight); + + } + + + } // end loop on track candidates + + + + // fill the leading particle histograms + + if(isInPeak && fFullmode) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading); + if(isInDZeroSideBand && fFullmode) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading); + + } // end loop on events in the pool + + }// end loop on D* candidates + + + + // check events with D* or SB canidates + if(fFullmode && EventHasDStarCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDStar"))->Fill(MultipOrCent,zVtxPosition); + if(fFullmode && EventHasDZeroSideBandCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDZeroSB"))->Fill(MultipOrCent,zVtxPosition); + + if(fFullmode && EventHasDStarCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDStarSB"))->Fill(MultipOrCent,zVtxPosition); + + + if(fFullmode) ((TH2F*)fOutput->FindObject("DStarCandidates"))->Fill(nOfDStarCandidates,MultipOrCent); + if(fFullmode) ((TH2F*)fOutput->FindObject("SBCandidates"))->Fill(nOfSBCandidates,MultipOrCent); + + // update event pool + Bool_t updated = fCorrelator->PoolUpdate(); + + // if(updated) EventMixingChecks(aodEvent); + if(!updated) AliInfo("Pool was not updated"); + + +} //end the exec + +//________________________________________ terminate ___________________________ +void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*) +{ + // The Terminate() function is the last function to be called during + // a query. It always runs on the client, it can be used to present + // the results graphically or save the results to file. + + AliAnalysisTaskSE::Terminate(); + + fOutput = dynamic_cast (GetOutputData(1)); + if (!fOutput) { + printf("ERROR: fOutput not available\n"); + return; + } + + return; +} +//_____________________________________________________ +Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const { + + //Daughter removal in MCKine case + Bool_t isDaughter = kFALSE; + Int_t labelD0 = d->GetLabel(); + + Int_t mother = track->GetMother(); + + //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!) + while (mother > 0){ + AliAODMCParticle* mcMoth = dynamic_cast(fmcArray->At(mother)); //it's the mother of the track! + if (mcMoth){ + if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE; + mother = mcMoth->GetMother(); //goes back by one + } else{ + AliError("Failed casting the mother particle!"); + break; + } + } + + return isDaughter; +} + +//_____________________________________________________ +void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){ + + //cout << "DEFINING THNSPARSES "<< endl; + + Double_t Pi = TMath::Pi(); + Int_t nbinscorr = fPhiBins; + Double_t lowcorrbin = -0.5*Pi; + Double_t upcorrbin = 1.5*Pi; + // define the THnSparseF + + //sparse bins + + //1 delta_phi + //2 delta_eta + //3 D* pt + //4 multiplicity + //5 track pt + //6 zVtx position + + Int_t nbinsPool = (fCuts2->GetNZvtxPoolBins())*(fCuts2->GetNCentPoolBins()); + + + Int_t nbinsSparse[5]= {nbinscorr, 32,30,250,nbinsPool}; + Double_t binLowLimitSparse[5]={lowcorrbin,-1.6, 0, 0,-0.5}; + Double_t binUpLimitSparse[5]= {upcorrbin, 1.6,30,25,nbinsPool-0.5}; + + Int_t MCnbinsSparse[7]= {nbinscorr, 32,30,250,nbinsPool,10,2}; + Double_t MCbinLowLimitSparse[7]={lowcorrbin,-1.6, 0, 0,-0.5,-0.5,-0.5}; // + Double_t MCbinUpLimitSparse[7]= {upcorrbin, 1.6,30,25,nbinsPool-0.5,9.5,1.5}; + + TString sparsename = "CorrelationsDStar"; + if(fselect==1) sparsename += "Hadron"; + if(fselect==2) sparsename += "Kaon"; + if(fselect==3) sparsename += "KZero"; + + TString D0Bkgsparsename = "DZeroBkg"; + D0Bkgsparsename += sparsename; + + TString DStarBkgsparsename = "DStarBkg"; + DStarBkgsparsename += sparsename; + + TString MCSparseName = "MCDStar"; + MCSparseName += sparsename; + // signal correlations + THnSparseF * Correlations = new THnSparseF(sparsename.Data(),"Correlations for signal",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse); + + // bkg correlations from D0 sidebands + THnSparseF * DZeroBkgCorrelations = new THnSparseF(D0Bkgsparsename.Data(),"Bkg Correlations estimated with D0 sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse); + + // bkg correlations from D* sidebands + THnSparseF * DStarBkgCorrelations = new THnSparseF(DStarBkgsparsename.Data(),"Bkg Correlations estimated with D* sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse); + + + THnSparseF * MCCorrelations = new THnSparseF(MCSparseName.Data(),"MC Correlations",7,MCnbinsSparse,MCbinLowLimitSparse,MCbinUpLimitSparse); + + MCCorrelations->GetAxis(5)->SetBinLabel(1," All "); + MCCorrelations->GetAxis(5)->SetBinLabel(2," from hadron Heavy flavour"); + MCCorrelations->GetAxis(5)->SetBinLabel(3," from c->D"); + MCCorrelations->GetAxis(5)->SetBinLabel(4," from b->D"); + MCCorrelations->GetAxis(5)->SetBinLabel(5," from b->B"); + MCCorrelations->GetAxis(5)->SetBinLabel(6," from quark Heavy flavour"); + MCCorrelations->GetAxis(5)->SetBinLabel(7," from c"); + MCCorrelations->GetAxis(5)->SetBinLabel(8," from b"); + + MCCorrelations->GetAxis(6)->SetBinLabel(1," if D* from c"); + MCCorrelations->GetAxis(6)->SetBinLabel(2," if D* from b"); + + Correlations->Sumw2(); + DZeroBkgCorrelations->Sumw2(); + DStarBkgCorrelations->Sumw2(); + + fOutput->Add(Correlations); + fOutput->Add(DZeroBkgCorrelations); + if(fFullmode) fOutput->Add(DStarBkgCorrelations); + if(fmontecarlo) fOutputMC->Add(MCCorrelations); + + +} +//__________________________________________________________________________________________________ +void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){ + + Double_t Pi = TMath::Pi(); + Int_t nbinscorr = fPhiBins; + Double_t lowcorrbin = -0.5*Pi ; // shift the bin by half the width so that at 0 is it the bin center + Double_t upcorrbin = 1.5*Pi ; + + // ========================= histograms for both Data and MonteCarlo + + + TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5); + NofEvents->GetXaxis()->SetBinLabel(1," All events"); + NofEvents->GetXaxis()->SetBinLabel(2," Selected events"); + NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup"); + NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality"); + NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx"); + NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr."); + NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc."); + NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger"); + NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel"); + NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk"); + NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx"); + fOutput->Add(NofEvents); + + + + + TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5); + if(!fmixing && fFullmode) fOutput->Add(D0InvMass); + + TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5); + if(!fmixing && fFullmode) fOutput->Add(D0InvMassinSB); + + //TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2); + //if(!fmixing) fOutput->Add(DeltaInvMass); + TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution; D* p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2); + if(!fmixing) fOutput->Add(DeltaInvMass); + + TH2F *bkgDeltaInvMass = new TH2F("bkgDeltaInvMass","K#pi#pi - K#pi invariant mass distribution; SB p_{T}; #DeltaInvMass",30,0,30,750,0.1,0.2); + if(!fmixing) fOutput->Add(bkgDeltaInvMass); + + DeltaInvMass->Sumw2(); + bkgDeltaInvMass->Sumw2(); + + TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50); + if(!fmixing) fOutput->Add(RecoPtDStar); + + TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50); + if(!fmixing) fOutput->Add(RecoPtBkg); + + TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50); + if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm); + + TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50); + if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty); + + TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50); + if(!fmixing) fOutput->Add(MCtagPtDStar); + + TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25); + if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectra); + + TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25); + if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectraifHF); + + TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","N of associated tracks per D trigger; Nof tracks; Entries",500,-0.5,499.5); + if(fFullmode) fOutput->Add(NofTracksInPeak); + + TH1D * NofTracksInSB = new TH1D("NofTracksInSB","N of associated tracks per SideBand trigger; Nof tracks; Entries",500,-0.5,499.5); + if(fFullmode) fOutput->Add(NofTracksInSB); + + TH1D * TracksInPeakSpectra = new TH1D("TracksInPeakSpectra","Pt Spectra tracks with D trigger; p_{T} GeV/c; Entries",500,-0.5,49.5); + if(fFullmode)fOutput->Add(TracksInPeakSpectra); + + TH1D * TracksInSBSpectra = new TH1D("TracksInSBSpectra","Pt Spectra tracks with SideBand trigger; p_{T} GeV/c; Entries",500,-0.5,49.5); + if(fFullmode)fOutput->Add(TracksInSBSpectra); + + + //TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5); + //if(fmixing) fOutput->Add(EventMixingCheck); + + + TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); + if(fFullmode)fOutput->Add(EventPropsCheck); + + TH2F * EventPropsCheckifDStar = new TH2F("EventPropsCheckifDStar","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); + if(fFullmode)fOutput->Add(EventPropsCheckifDStar); + + TH2F * EventPropsCheckifDZeroSB = new TH2F("EventPropsCheckifDZeroSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); + if(fFullmode)fOutput->Add(EventPropsCheckifDZeroSB); + + TH2F * EventPropsCheckifDStarSB = new TH2F("EventPropsCheckifDStarSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10); + if(fFullmode)fOutput->Add(EventPropsCheckifDStarSB); + + + TH2F * WeightChecks = new TH2F("WeightChecks","Checks on efficiency correction",300,0,30,100,0.005,1.005); + if(fFullmode)fOutput->Add(WeightChecks); + + + + TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9); + fOutput->Add(PhiEtaTrigger); + + TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9); + fOutput->Add(PhiEtaSideBand); + + TH3F * PhiEtaPart = new TH3F("PhiEtaPart","#phi distribution of the associated particle; #phi; #eta; multiplicity",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9,100,0,1000); + fOutput->Add(PhiEtaPart); + + TH2F * DStarCandidates = new TH2F("DStarCandidates","# of D* candidates per event vs multiplicity",6,-0.5,5.5,50,0,500); + if(fFullmode)fOutput->Add(DStarCandidates); + + TH2F * SBCandidates = new TH2F("SBCandidates","# of SB candidates per event vs multiplicity",6,-0.5,5.5,50,0,500); + if(fFullmode)fOutput->Add(SBCandidates); + + + //correlations histograms + TString histoname1 = "DPhiDStar"; + if(fselect==1) histoname1 += "Hadron"; + if(fselect==2) histoname1 += "Kaon"; + if(fselect==3) histoname1 += "KZero"; + + /* + TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + //side band background histograms + TString histoname2 = "bkg"; + histoname2 += histoname1; + TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + + fOutput->Add(DPhiDStar); + + if(fselect==3){fOutput->Add(DPhiDStarKZero1);} + + fOutput->Add(bkgDPhiDStar); + + if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);} + + */ + // leading particle + TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + if(fFullmode)fOutput->Add(leadingcand); + if(fFullmode)fOutput->Add(leadingsidebands); + + // ========================= histos for analysis on MC only + + TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5); + if(fmontecarlo) fOutputMC->Add(EventTypeMC); + + TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5); + MCSources->GetXaxis()->SetBinLabel(1," All "); + MCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour"); + MCSources->GetXaxis()->SetBinLabel(3," from c->D"); + MCSources->GetXaxis()->SetBinLabel(4," from b->D"); + MCSources->GetXaxis()->SetBinLabel(5," from b->B"); + MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour"); + MCSources->GetXaxis()->SetBinLabel(7," from c"); + MCSources->GetXaxis()->SetBinLabel(8," from b"); + + if(fmontecarlo) fOutputMC->Add(MCSources); + + // leading particle from mc source + TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5); + LeadingMCSources->GetXaxis()->SetBinLabel(1," All "); + LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour"); + LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D"); + LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D"); + LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B"); + LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour"); + LeadingMCSources->GetXaxis()->SetBinLabel(7," from c"); + LeadingMCSources->GetXaxis()->SetBinLabel(8," from b"); + + if(fmontecarlo && fFullmode) fOutputMC->Add(LeadingMCSources); + + // all hadrons + TString histoname3 = "MCTag"; + histoname3 += histoname1; + TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString histoname44 = "CharmDOrigin"; + histoname44 += histoname1; + histoname44 += "MC"; + + TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + + TString histoname54 = "BeautyDOrigin"; + histoname54 += histoname1; + histoname54 += "MC"; + TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString histoname55 = "BeautyBOrigin"; + histoname55 += histoname1; + histoname55 += "MC"; + TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString histoname4 = "CharmQuarkOrigin"; + histoname4 += histoname1; + histoname4 += "MC"; + TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString histoname5 = "BeautyQuarkOrigin"; + histoname5 += histoname1; + histoname5 += "MC"; + TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString histoname6 = "NonHFOrigin"; + histoname6 += histoname1; + histoname6 += "MC"; + TH3D * NonHFOriginDPhiDStar = new TH3D(histoname6.Data(),histoname6.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + if(fmontecarlo && fFullmode){ + + fOutputMC->Add(MCTagDPhiDStar); + fOutputMC->Add(CharmDOriginDPhiDStar); + fOutputMC->Add(BeautyDOriginDPhiDStar); + fOutputMC->Add(BeautyBOriginDPhiDStar); + fOutputMC->Add(CharmQuarkOriginDPhiDStar); + fOutputMC->Add(BeautyQuarkOriginDPhiDStar); + fOutputMC->Add(NonHFOriginDPhiDStar); + + } + + // ========================= histos for analysis on MC + // all leading hadron + TString Leadinghistoname3 = "LeadingMCTag"; + Leadinghistoname3 += histoname1; + TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString Leadinghistoname44 = "LeadingCharmDOrigin"; + Leadinghistoname44 += histoname1; + Leadinghistoname44 += "MC"; + + TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + + TString Leadinghistoname54 = "LeadingBeautyDOrigin"; + Leadinghistoname54 += histoname1; + Leadinghistoname54 += "MC"; + TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString Leadinghistoname55 = "LeadingBeautyBOrigin"; + Leadinghistoname55 += histoname1; + Leadinghistoname55 += "MC"; + TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString Leadinghistoname4 = "LeadingCharmQuarkOrigin"; + Leadinghistoname4 += histoname1; + Leadinghistoname4 += "MC"; + TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin"; + Leadinghistoname5 += histoname1; + Leadinghistoname5 += "MC"; + TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2); + + + + + if(fmontecarlo && fFullmode){ + + fOutputMC->Add(LeadingMCTagDPhiDStar); + fOutputMC->Add(LeadingCharmDOriginDPhiDStar); + fOutputMC->Add(LeadingBeautyDOriginDPhiDStar); + fOutputMC->Add(LeadingBeautyBOriginDPhiDStar); + fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar); + fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar); + + } + + TH3F * MCPhiEtaPart = new TH3F("MCPhiEtaPart","#phi distribution of the associated particle",nbinscorr,lowcorrbin,upcorrbin,50,-2.5,2.5,6,-0.5,6.5); + MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles"); + MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark"); + MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark"); + MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c"); + MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b"); + MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b"); + if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart); + + TH2D * TrackLabels = new TH2D("TrackLabels","NofEvents;track label; multiplicity",2,-0.5,1.5,500,-0.5,499.5); + if(fmontecarlo && fFullmode) fOutputMC->Add(TrackLabels); + + // ============================= EVENT MIXING CHECKS ====================================== + + Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool(); + Int_t MinNofTracks = fCuts2->GetMinNTracksInPool(); + Int_t NofCentBins = fCuts2->GetNCentPoolBins(); + Double_t * CentBins = fCuts2->GetCentPoolBins(); + Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins(); + Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins(); + + + + Int_t k =0; + + if(fSystem == AA) k = 100; // PbPb centrality + if(fSystem == pp || fSystem == pA) k = NofCentBins; // pp multiplicity + + + //Double_t minvalue = CentBins[0]; + //Double_t maxvalue = CentBins[NofCentBins+1]; + //Double_t Zminvalue = ZVrtxBins[0]; + //Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1]; + + Double_t minvalue, maxvalue; + Double_t Zminvalue, Zmaxvalue; + + Zminvalue = -15.; + Zmaxvalue = 15; + if(fSystem == AA) {minvalue = 0; maxvalue = 100;} // PbPb + if(fSystem == pp || fSystem == pA) {minvalue = 0; maxvalue = 500;} // multilpicity + + //Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents}; + // Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents}; + // Double_t * events = Nevents; + Double_t eventsv[] ={0,1000000}; + //Double_t * events = new Double_t[2]; + // events[0] = 0; +// events[1] = 1000000; + Double_t *events = eventsv; + Int_t Nevents = 1000000; + // TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,Nevents,events); + + TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,Nevents,events[0],events[1]); + + EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity "); + EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]"); + EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin"); + if(fmixing && fFullmode) fOutput->Add(EventsPerPoolBin); + + Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks; + //Int_t Diff = MaxNofTracks-MinNofTracks; + + //Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks}; + // Double_t * trackN = Ntracks; + + TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,MaxNofTracks,0,MaxNofTracks); + NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity "); + NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]"); + NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin"); + + if(fmixing && fFullmode) fOutput->Add(NofTracksPerPoolBin); + + TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Calls per pool bin",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins); + NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity "); + NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]"); + if(fmixing && fFullmode) fOutput->Add(NofPoolBinCalls); + + + + TH2D * EventProps = new TH2D("EventProps","Event properties",100,minvalue,maxvalue,100,Zminvalue,Zmaxvalue); + EventProps->GetXaxis()->SetTitle("Centrality/multiplicity "); + EventProps->GetYaxis()->SetTitle("Z vertex [cm]"); + if(fmixing && fFullmode) fOutput->Add(EventProps); + + TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5); + CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready"); + CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready"); + CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready"); + CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready"); + + if(fmixing) fOutput->Add(CheckPoolReadiness); + + +} + +//__________________________________________________________________________________________________ +void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){ + + + //Float_t* ptbins = fCuts->GetPtBinLimits(); + if(fD0Window) delete fD0Window; + fD0Window = new Float_t[fNofPtBins]; + + AliInfo("Enlarging the D0 mass windows from cut object\n"); + Int_t nvars = fCuts->GetNVars(); + + if(nvars<1){ + AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!"); + return; + } + Float_t** rdcutsvalmine=new Float_t*[nvars]; + for(Int_t iv=0;ivGetCutValue(0,j); + rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j); + cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl; + } + else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j); + + // set same windows + //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j); + } + } + + fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine); + + AliInfo("\n New windows set\n"); + fCuts->PrintAll(); + + + for(Int_t iv=0;ivGetNTracks(); + MultipOrCent = multiplicity; // convert from Int_t to Double_t + } + if(fSystem == AA){ // PbPb + + centralityObj = AOD->GetHeader()->GetCentralityP(); + MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M"); + AliInfo(Form("Centrality is %f", MultipOrCent)); + } + + AliAODVertex *vtx = AOD->GetPrimaryVertex(); + Double_t zvertex = vtx->GetZ(); // zvertex + + + + + AliEventPool * pool = fCorrelator->GetPool(); + + + + + ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool + ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties + + ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool + ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool +} + + + + + diff --git a/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.h b/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.h index bf0cb3fc889..748d85c21a3 100644 --- a/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.h +++ b/PWGHF/correlationHF/AliAnalysisTaskDStarCorrelations.h @@ -1,166 +1,166 @@ -#ifndef AliAnalysisTaskDStarCorrelations_H -#define AliAnalysisTaskDStarCorrelations_H - -/************************************************************************** - * Copyright(c) 1998-2009, 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. * - **************************************************************************/ - -/* $Id$ */ - -//----------------------------------------------------------------------- -// -// -// Author S.Bjelogrlic -// Utrecht University -// sandro.bjelogrlic@cern.ch -// -//----------------------------------------------------------------------- - -#include -#include -#include -#include "AliAnalysisTaskSE.h" -#include "AliAODEvent.h" -#include "AliRDHFCutsDStartoKpipi.h" -#include "AliEventPoolManager.h" -#include "AliAODRecoCascadeHF.h" -#include "AliHFAssociatedTrackCuts.h" -#include "AliNormalizationCounter.h" -#include "AliHFCorrelator.h" -#include -#include "AliAnalysisUtils.h" -#include "AliVertexingHFUtils.h" -class TParticle ; -class TClonesArray ; -class AliAODMCParticle; -class AliAODEvent; -class AliVParticle; -class TObjArray; -class AliEventPoolManager; -class AliESDEvent; - - - -class AliAnalysisTaskDStarCorrelations : public AliAnalysisTaskSE -{ - - public : - - enum CollSyst {pp,pA,AA}; - enum DEffVariable{kNone,kMult,kCentr,kRapidity,kEta}; - - AliAnalysisTaskDStarCorrelations(); - AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts, AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode); - virtual ~AliAnalysisTaskDStarCorrelations(); - - - - // Implementation of interface methods - virtual void UserCreateOutputObjects(); - virtual void Init(); - virtual void LocalInit() {Init();} - virtual void UserExec(Option_t *option); - virtual void Terminate(Option_t *option); - - void DefineThNSparseForAnalysis(); - void DefineHistoForAnalysis(); - void EnlargeDZeroMassWindow(); - Bool_t IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const ; - - // checker for event mixing - void EventMixingChecks(AliAODEvent * AOD); - // setters - void SetMonteCarlo(Bool_t k) {fmontecarlo = k;} - void SetUseMixing (Bool_t j) {fmixing = j;} - void SetCorrelator(Int_t l) {fselect = l;} // select 1 for hadrons, 2 for Kaons, 3 for Kzeros - void SetUseDisplacement(Int_t m) {fDisplacement=m;} // select 0 for no displ, 1 for abs displ, 2 for d0/sigma_d0 - void SetCollSys(CollSyst system){fSystem=system;} // select between pp (kFALSE) or PbPb (kTRUE) - void SetEfficiencyVariable(DEffVariable var){fEfficiencyVariable = var;} // set the efficiency variable to use - void SetLevelOfDebug(Int_t debug){fDebugLevel=debug;} // set debug level - void SetUseReconstruction(Bool_t reco){fReco = reco;} - void SetDMesonSigmas(Float_t DStarWin, Float_t D0Win, Float_t SBmin, Float_t SBmax){ - fDMesonSigmas[0] = DStarWin; - fDMesonSigmas[1] = D0Win; - fDMesonSigmas[2] = SBmin; - fDMesonSigmas[3] = SBmax; - - } - - void SetUseEfficiencyCorrection(Bool_t correction){fUseEfficiencyCorrection = correction;} // setter for using the single track efficiency correction - void SetUseDmesonEfficiencyCorrection(Bool_t correction){fUseDmesonEfficiencyCorrection = correction;} // setter for using the single track efficiency correction - void SetUseHadronicChannelAtKineLevel (Bool_t use){fUseHadronicChannelAtKineLevel = use;} - - void SetDim(){fDim = 4; - fDMesonSigmas = new Float_t[4];} - void SetDeffMapvsPt(TH1D * map){fDeffMapvsPt = map;} - void SetDeffMapvsPtvsMult(TH2D * map){fDeffMapvsPtvsMult = (TH2D*)map;} - void SetDeffMapvsPtvsMultvsEta(TH2D * map){fDeffMapvsPtvsEta = map;} - void SetNofPhiBins(Int_t nbins){fPhiBins = nbins;} - void SetMaxDStarEta(Double_t eta){fMaxEtaDStar = eta;} - - - -private: - - AliAnalysisTaskDStarCorrelations(const AliAnalysisTaskDStarCorrelations &source); - AliAnalysisTaskDStarCorrelations& operator=(const AliAnalysisTaskDStarCorrelations& source); - - TObject* fhandler; //! Analysis Handler - TClonesArray* fmcArray; //mcarray - AliNormalizationCounter *fCounter; // counter - AliHFCorrelator * fCorrelator; // object for correlations - - - - Int_t fselect; // select what to correlate with a D* 1-chargedtracks,2-chargedkaons,3-k0s - Bool_t fmontecarlo;//switch for MC - Bool_t fmixing;// switch for event mixing - Bool_t fFullmode; - CollSyst fSystem; // pp, pPb or PbPb - DEffVariable fEfficiencyVariable; // set second variable to study efficiency (mult, centr, y, eta) - Bool_t fReco; // use reconstruction or MC truth - Bool_t fUseEfficiencyCorrection; // boolean variable to use or not the efficiency correction - Bool_t fUseDmesonEfficiencyCorrection; // boolean flag for the use of Dmeson efficiency correction - Bool_t fUseCentrality;// boolean to switch in between centrality or multiplicity - Bool_t fUseHadronicChannelAtKineLevel; // - - Int_t fPhiBins; - Int_t fEvents; //! number of event - Int_t fDebugLevel; //! debug level - Int_t fDisplacement; // set 0 for no displacement cut, 1 for absolute d0, 2 for d0/sigma_d0 - Int_t fDim;// - Int_t fNofPtBins; - Double_t fMaxEtaDStar; - Float_t *fDMesonSigmas;//[fDim] - Float_t * fD0Window; //[fNofPtBins] - - - - - TList *fOutput; //! user output data - TList *fOutputMC; //! outpu for MC - AliRDHFCutsDStartoKpipi *fCuts; // Cuts D* - AliHFAssociatedTrackCuts *fCuts2; // cuts for associated - AliAnalysisUtils *fUtils; - AliAODTracklets * fTracklets; // AliAODtracklets - - TH1D * fDeffMapvsPt; // histo for Deff mappin - TH2D * fDeffMapvsPtvsMult; // histo for Deff mappin - TH2D * fDeffMapvsPtvsEta; // histo for Deff mappin - - ClassDef(AliAnalysisTaskDStarCorrelations,6); // class for D meson correlations - -}; - -#endif +#ifndef AliAnalysisTaskDStarCorrelations_H +#define AliAnalysisTaskDStarCorrelations_H + +/************************************************************************** + * Copyright(c) 1998-2009, 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. * + **************************************************************************/ + +/* $Id$ */ + +//----------------------------------------------------------------------- +// +// +// Author S.Bjelogrlic +// Utrecht University +// sandro.bjelogrlic@cern.ch +// +//----------------------------------------------------------------------- + +#include +#include +#include +#include "AliAnalysisTaskSE.h" +#include "AliAODEvent.h" +#include "AliRDHFCutsDStartoKpipi.h" +#include "AliEventPoolManager.h" +#include "AliAODRecoCascadeHF.h" +#include "AliHFAssociatedTrackCuts.h" +#include "AliNormalizationCounter.h" +#include "AliHFCorrelator.h" +#include +#include "AliAnalysisUtils.h" +#include "AliVertexingHFUtils.h" +class TParticle ; +class TClonesArray ; +class AliAODMCParticle; +class AliAODEvent; +class AliVParticle; +class TObjArray; +class AliEventPoolManager; +class AliESDEvent; + + + +class AliAnalysisTaskDStarCorrelations : public AliAnalysisTaskSE +{ + + public : + + enum CollSyst {pp,pA,AA}; + enum DEffVariable{kNone,kMult,kCentr,kRapidity,kEta}; + + AliAnalysisTaskDStarCorrelations(); + AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts, AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode); + virtual ~AliAnalysisTaskDStarCorrelations(); + + + + // Implementation of interface methods + virtual void UserCreateOutputObjects(); + virtual void Init(); + virtual void LocalInit() {Init();} + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *option); + + void DefineThNSparseForAnalysis(); + void DefineHistoForAnalysis(); + void EnlargeDZeroMassWindow(); + Bool_t IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const ; + + // checker for event mixing + void EventMixingChecks(AliAODEvent * AOD); + // setters + void SetMonteCarlo(Bool_t k) {fmontecarlo = k;} + void SetUseMixing (Bool_t j) {fmixing = j;} + void SetCorrelator(Int_t l) {fselect = l;} // select 1 for hadrons, 2 for Kaons, 3 for Kzeros + void SetUseDisplacement(Int_t m) {fDisplacement=m;} // select 0 for no displ, 1 for abs displ, 2 for d0/sigma_d0 + void SetCollSys(CollSyst system){fSystem=system;} // select between pp (kFALSE) or PbPb (kTRUE) + void SetEfficiencyVariable(DEffVariable var){fEfficiencyVariable = var;} // set the efficiency variable to use + void SetLevelOfDebug(Int_t debug){fDebugLevel=debug;} // set debug level + void SetUseReconstruction(Bool_t reco){fReco = reco;} + void SetDMesonSigmas(Float_t DStarWin, Float_t D0Win, Float_t SBmin, Float_t SBmax){ + fDMesonSigmas[0] = DStarWin; + fDMesonSigmas[1] = D0Win; + fDMesonSigmas[2] = SBmin; + fDMesonSigmas[3] = SBmax; + + } + + void SetUseEfficiencyCorrection(Bool_t correction){fUseEfficiencyCorrection = correction;} // setter for using the single track efficiency correction + void SetUseDmesonEfficiencyCorrection(Bool_t correction){fUseDmesonEfficiencyCorrection = correction;} // setter for using the single track efficiency correction + void SetUseHadronicChannelAtKineLevel (Bool_t use){fUseHadronicChannelAtKineLevel = use;} + + void SetDim(){fDim = 4; + fDMesonSigmas = new Float_t[4];} + void SetDeffMapvsPt(TH1D * map){fDeffMapvsPt = map;} + void SetDeffMapvsPtvsMult(TH2D * map){fDeffMapvsPtvsMult = (TH2D*)map;} + void SetDeffMapvsPtvsMultvsEta(TH2D * map){fDeffMapvsPtvsEta = map;} + void SetNofPhiBins(Int_t nbins){fPhiBins = nbins;} + void SetMaxDStarEta(Double_t eta){fMaxEtaDStar = eta;} + + + +private: + + AliAnalysisTaskDStarCorrelations(const AliAnalysisTaskDStarCorrelations &source); + AliAnalysisTaskDStarCorrelations& operator=(const AliAnalysisTaskDStarCorrelations& source); + + TObject* fhandler; //! Analysis Handler + TClonesArray* fmcArray; //mcarray + AliNormalizationCounter *fCounter; // counter + AliHFCorrelator * fCorrelator; // object for correlations + + + + Int_t fselect; // select what to correlate with a D* 1-chargedtracks,2-chargedkaons,3-k0s + Bool_t fmontecarlo;//switch for MC + Bool_t fmixing;// switch for event mixing + Bool_t fFullmode; + CollSyst fSystem; // pp, pPb or PbPb + DEffVariable fEfficiencyVariable; // set second variable to study efficiency (mult, centr, y, eta) + Bool_t fReco; // use reconstruction or MC truth + Bool_t fUseEfficiencyCorrection; // boolean variable to use or not the efficiency correction + Bool_t fUseDmesonEfficiencyCorrection; // boolean flag for the use of Dmeson efficiency correction + Bool_t fUseCentrality;// boolean to switch in between centrality or multiplicity + Bool_t fUseHadronicChannelAtKineLevel; // + + Int_t fPhiBins; + Int_t fEvents; //! number of event + Int_t fDebugLevel; //! debug level + Int_t fDisplacement; // set 0 for no displacement cut, 1 for absolute d0, 2 for d0/sigma_d0 + Int_t fDim;// + Int_t fNofPtBins; + Double_t fMaxEtaDStar; + Float_t *fDMesonSigmas;//[fDim] + Float_t * fD0Window; //[fNofPtBins] + + + + + TList *fOutput; //! user output data + TList *fOutputMC; //! outpu for MC + AliRDHFCutsDStartoKpipi *fCuts; // Cuts D* + AliHFAssociatedTrackCuts *fCuts2; // cuts for associated + AliAnalysisUtils *fUtils; + AliAODTracklets * fTracklets; // AliAODtracklets + + TH1D * fDeffMapvsPt; // histo for Deff mappin + TH2D * fDeffMapvsPtvsMult; // histo for Deff mappin + TH2D * fDeffMapvsPtvsEta; // histo for Deff mappin + + ClassDef(AliAnalysisTaskDStarCorrelations,6); // class for D meson correlations + +}; + +#endif diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSEDs.cxx b/PWGHF/vertexingHF/AliAnalysisTaskSEDs.cxx index 78b81c1a64a..7c414964375 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskSEDs.cxx +++ b/PWGHF/vertexingHF/AliAnalysisTaskSEDs.cxx @@ -505,7 +505,7 @@ void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/) if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3); if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4); if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5); - if(fAnalysisCuts->IsEventRejectedDueToPileupSPD())fHistNEvents->Fill(6); + if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6); if(fAnalysisCuts->IsEventRejectedDueToCentrality()){ fHistNEvents->Fill(7); fHistCentrality[2]->Fill(evCentr); diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSEHFQA.cxx b/PWGHF/vertexingHF/AliAnalysisTaskSEHFQA.cxx index 15e80a75275..728c1b380f6 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskSEHFQA.cxx +++ b/PWGHF/vertexingHF/AliAnalysisTaskSEHFQA.cxx @@ -1537,7 +1537,7 @@ void AliAnalysisTaskSEHFQA::UserExec(Option_t */*option*/) //select event if(!fCuts->IsEventSelected(aod)) { evSelected=kFALSE; - if(fCuts->IsEventRejectedDueToPileupSPD()) { + if(fCuts->IsEventRejectedDueToPileup()) { if(hWhyEvRejected) hWhyEvRejected->Fill(0); evselByPileup=kFALSE; }// rejected for pileup diff --git a/PWGHF/vertexingHF/AliRDHFCuts.cxx b/PWGHF/vertexingHF/AliRDHFCuts.cxx index a96c6903c59..964ed7d1953 100644 --- a/PWGHF/vertexingHF/AliRDHFCuts.cxx +++ b/PWGHF/vertexingHF/AliRDHFCuts.cxx @@ -43,6 +43,7 @@ #include "AliAnalysisManager.h" #include "AliInputEventHandler.h" #include "AliPIDResponse.h" +#include "AliAnalysisUtils.h" #include "TRandom.h" #include @@ -572,7 +573,16 @@ Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) { Double_t cutz=(Double_t)fMinDzPileup; if(event->IsPileupFromSPD(cutc,cutz,3.,2.,10.)) { if(accept) fWhyRejection=1; - fEvRejectionBits+=1<