-/**************************************************************************\r
- * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
- * *\r
- * Author: The ALICE Off-line Project. *\r
- * Contributors are mentioned in the code where appropriate. *\r
- * *\r
- * Permission to use, copy, modify and distribute this software and its *\r
- * documentation strictly for non-commercial purposes is hereby granted *\r
- * without fee, provided that the above copyright notice appears in all *\r
- * copies and that both the copyright notice and this permission notice *\r
- * appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is *\r
- * provided "as is" without express or implied warranty. *\r
- **************************************************************************/\r
-//\r
-//\r
-// Base class for DStar - Hadron Correlations Analysis\r
-//\r
-//-----------------------------------------------------------------------\r
-// \r
-//\r
-// Author S.Bjelogrlic\r
-// Utrecht University \r
-// sandro.bjelogrlic@cern.ch\r
-//\r
-//-----------------------------------------------------------------------\r
-\r
-/* $Id$ */\r
-\r
-//#include <TDatabasePDG.h>\r
-#include <TParticle.h>\r
-#include <TVector3.h>\r
-#include <TChain.h>\r
-#include "TROOT.h"\r
-\r
-#include "AliAnalysisTaskDStarCorrelations.h"\r
-#include "AliRDHFCutsDStartoKpipi.h"\r
-#include "AliHFAssociatedTrackCuts.h"\r
-#include "AliAODRecoDecay.h"\r
-#include "AliAODRecoCascadeHF.h"\r
-#include "AliAODRecoDecayHF2Prong.h"\r
-#include "AliAODPidHF.h"\r
-#include "AliVParticle.h"\r
-#include "AliAnalysisManager.h"\r
-#include "AliAODInputHandler.h"\r
-#include "AliAODHandler.h"\r
-#include "AliESDtrack.h"\r
-#include "AliAODMCParticle.h"\r
-#include "AliNormalizationCounter.h"\r
-#include "AliReducedParticle.h"\r
-#include "AliHFCorrelator.h"\r
-#include "AliAODMCHeader.h"\r
-#include "AliEventPoolManager.h"\r
-#include "AliVertexingHFUtils.h"\r
-\r
-using std::cout;\r
-using std::endl;\r
-\r
-\r
-ClassImp(AliAnalysisTaskDStarCorrelations)\r
-\r
-\r
-//__________________________________________________________________________\r
-AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations() :\r
-AliAnalysisTaskSE(),\r
-fhandler(0x0),\r
-fmcArray(0x0),\r
-fCounter(0x0),\r
-fCorrelator(0x0),\r
-fselect(0),\r
-fmontecarlo(kFALSE),\r
-fmixing(kFALSE),\r
-fFullmode(kFALSE),\r
-fSystem(pp),\r
-fEfficiencyVariable(kNone),\r
-fReco(kTRUE),\r
-fUseEfficiencyCorrection(kFALSE),\r
-fUseDmesonEfficiencyCorrection(kFALSE),\r
-fUseCentrality(kFALSE),\r
-fUseHadronicChannelAtKineLevel(kFALSE),\r
-fPhiBins(32),\r
-fEvents(0),\r
-fDebugLevel(0),\r
-fDisplacement(0),\r
-fDim(0),\r
-fNofPtBins(0),\r
-fMaxEtaDStar(0.9),\r
-fDMesonSigmas(0),\r
-\r
-fD0Window(0),\r
-\r
-fOutput(0x0),\r
-fOutputMC(0x0),\r
-fCuts(0),\r
-fCuts2(0),\r
-fUtils(0),\r
-fTracklets(0),\r
-fDeffMapvsPt(0),\r
-fDeffMapvsPtvsMult(0),\r
-fDeffMapvsPtvsEta(0)\r
-{\r
- SetDim();\r
- // default constructor\r
-}\r
-\r
-//__________________________________________________________________________\r
-AliAnalysisTaskDStarCorrelations::AliAnalysisTaskDStarCorrelations(const Char_t* name,AliRDHFCutsDStartoKpipi* cuts, AliHFAssociatedTrackCuts *AsscCuts,AliAnalysisTaskDStarCorrelations::CollSyst syst,Bool_t mode) :\r
-AliAnalysisTaskSE(name),\r
-\r
-fhandler(0x0),\r
-fmcArray(0x0),\r
-fCounter(0x0),\r
-fCorrelator(0x0),\r
-fselect(0),\r
-fmontecarlo(kFALSE),\r
-fmixing(kFALSE),\r
-fFullmode(mode),\r
-fSystem(syst),\r
-fEfficiencyVariable(kNone),\r
-fReco(kTRUE),\r
-fUseEfficiencyCorrection(kFALSE),\r
-fUseDmesonEfficiencyCorrection(kFALSE),\r
-fUseCentrality(kFALSE),\r
-fUseHadronicChannelAtKineLevel(kFALSE),\r
-fPhiBins(32),\r
-fEvents(0),\r
-fDebugLevel(0),\r
-fDisplacement(0),\r
-fDim(0),\r
-fNofPtBins(0),\r
-fMaxEtaDStar(0.9),\r
-fDMesonSigmas(0),\r
-fD0Window(0),\r
-\r
-fOutput(0x0),\r
-fOutputMC(0x0),\r
-fCuts(0),\r
-fCuts2(AsscCuts),\r
-fUtils(0),\r
-fTracklets(0),\r
-fDeffMapvsPt(0),\r
-fDeffMapvsPtvsMult(0),\r
-fDeffMapvsPtvsEta(0)\r
-{\r
- Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");\r
- \r
- SetDim();\r
- if(fSystem == AA) fUseCentrality = kTRUE; else fUseCentrality = kFALSE;\r
- \r
- fCuts=cuts;\r
- fNofPtBins= fCuts->GetNPtBins();\r
- //cout << "Enlarging the DZero window " << endl;\r
- EnlargeDZeroMassWindow();\r
- // cout << "Done" << endl;\r
- \r
- \r
- DefineInput(0, TChain::Class());\r
- \r
- DefineOutput(1,TList::Class()); // histos from data and MC\r
- DefineOutput(2,TList::Class()); // histos from MC\r
- DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts\r
- DefineOutput(4,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts\r
- DefineOutput(5,AliNormalizationCounter::Class()); // normalization\r
-}\r
-\r
-//__________________________________________________________________________\r
-\r
-AliAnalysisTaskDStarCorrelations::~AliAnalysisTaskDStarCorrelations() {\r
- //\r
- // destructor\r
- //\r
- \r
- Info("AliAnalysisTaskDStarCorrelations","Calling Destructor"); \r
- \r
- if(fhandler) {delete fhandler; fhandler = 0;}\r
- //if(fPoolMgr) {delete fPoolMgr; fPoolMgr = 0;} \r
- if(fmcArray) {delete fmcArray; fmcArray = 0;}\r
- if(fCounter) {delete fCounter; fCounter = 0;}\r
- if(fCorrelator) {delete fCorrelator; fCorrelator = 0;}\r
- if(fOutput) {delete fOutput; fOutput = 0;}\r
- if(fOutputMC) {delete fOutputMC; fOutputMC = 0;}\r
- if(fCuts) {delete fCuts; fCuts = 0;}\r
- if(fCuts2) {delete fCuts2; fCuts2=0;}\r
- if(fDeffMapvsPt){delete fDeffMapvsPt; fDeffMapvsPt=0;}\r
- if(fDeffMapvsPtvsMult){delete fDeffMapvsPtvsMult; fDeffMapvsPtvsMult=0;}\r
- if(fDeffMapvsPtvsEta){delete fDeffMapvsPtvsEta; fDeffMapvsPtvsEta=0;}\r
-\r
-}\r
-\r
-//___________________________________________________________\r
-void AliAnalysisTaskDStarCorrelations::Init(){\r
- //\r
- // Initialization\r
- //\r
- if(fDebugLevel > 1) printf("AliAnalysisTaskDStarCorrelations::Init() \n");\r
- \r
- AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);\r
- \r
- \r
- \r
- // Post the D* cuts\r
- PostData(3,copyfCuts);\r
- \r
- // Post the hadron cuts\r
- PostData(4,fCuts2);\r
- \r
- return;\r
-}\r
-\r
-\r
-//_________________________________________________\r
-void AliAnalysisTaskDStarCorrelations::UserCreateOutputObjects(){\r
- Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r
- \r
- //slot #1\r
- //OpenFile(0);\r
- fOutput = new TList();\r
- fOutput->SetOwner();\r
- \r
- fOutputMC = new TList();\r
- fOutputMC->SetOwner();\r
- \r
- // define histograms\r
- DefineHistoForAnalysis();\r
- DefineThNSparseForAnalysis();\r
- \r
-\r
- fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));\r
- fCounter->Init();\r
- \r
- Double_t Pi = TMath::Pi();\r
- fCorrelator = new AliHFCorrelator("Correlator",fCuts2,fUseCentrality); // fCuts2 is the hadron cut object, fSystem to switch between pp or PbPb\r
- fCorrelator->SetDeltaPhiInterval( -0.5*Pi, 1.5*Pi); // set correct phi interval\r
- //fCorrelator->SetDeltaPhiInterval((-0.5)*Pi,(1.5)*Pi); // set correct phi interval\r
- fCorrelator->SetEventMixing(fmixing); //set kFALSE/kTRUE for mixing Off/On\r
- fCorrelator->SetAssociatedParticleType(fselect); // set 1/2/3 for hadron/kaons/kzeros\r
- fCorrelator->SetApplyDisplacementCut(fDisplacement); //set kFALSE/kTRUE for using the displacement cut\r
- fCorrelator->SetUseMC(fmontecarlo);\r
- fCorrelator->SetUseReco(fReco);\r
- // fCorrelator->SetKinkRemoval(kTRUE);\r
- Bool_t pooldef = fCorrelator->DefineEventPool();\r
- \r
- if(!pooldef) AliInfo("Warning:: Event pool not defined properly");\r
- \r
- fUtils = new AliAnalysisUtils();\r
- \r
- \r
- \r
- PostData(1,fOutput); // set the outputs\r
- PostData(2,fOutputMC); // set the outputs\r
- PostData(5,fCounter); // set the outputs\r
-}\r
-\r
-//_________________________________________________\r
-void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){\r
- \r
- \r
- if(fDebugLevel){\r
- \r
- if(fReco) std::cout << "USING RECONSTRUCTION" << std::endl;\r
- if(!fReco) std::cout << "USING MC TRUTH" << std::endl;\r
- std::cout << " " << std::endl;\r
- std::cout << "=================================================================================" << std::endl;\r
- if(!fmixing){\r
- if(fselect==1) std::cout << "TASK::Correlation with hadrons on SE "<< std::endl;\r
- if(fselect==2) std::cout << "TASK::Correlation with kaons on SE "<< std::endl;\r
- if(fselect==3) std::cout << "TASK::Correlation with kzeros on SE "<< std::endl;\r
- }\r
- if(fmixing){\r
- if(fselect==1) std::cout << "TASK::Correlation with hadrons on ME "<< std::endl;\r
- if(fselect==2) std::cout << "TASK::Correlation with kaons on ME "<< std::endl;\r
- if(fselect==3) std::cout << "TASK::Correlation with kzeros on ME "<< std::endl;\r
- }\r
- \r
- }// end if debug\r
- \r
- \r
- \r
- if (!fInputEvent) {\r
- Error("UserExec","NO EVENT FOUND!");\r
- return;\r
- }\r
- \r
- AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r
- if(!aodEvent){\r
- AliError("AOD event not found!");\r
- return;\r
- }\r
- \r
- fTracklets = aodEvent->GetTracklets();\r
- \r
- fEvents++; // event counter\r
- ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);\r
- \r
- fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo);\r
- \r
- // load MC array\r
- fmcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r
- if(fmontecarlo && !fmcArray){\r
- AliError("Array of MC particles not found");\r
- return;\r
- }\r
- \r
- \r
- \r
- \r
- // ********************************************** START EVENT SELECTION ****************************************************\r
- \r
- Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);\r
- \r
- if(!isEvSel) {\r
- \r
- if(fCuts->IsEventRejectedDueToPileupSPD()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(2);\r
- if(fCuts->IsEventRejectedDueToCentrality()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3);\r
- if(fCuts->IsEventRejectedDueToNotRecoVertex()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4);\r
- if(fCuts->IsEventRejectedDueToVertexContributors()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5);\r
- if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6);\r
- if(fCuts->IsEventRejectedDueToTrigger()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7);\r
- if(fCuts->IsEventRejectedDuePhysicsSelection()) ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8);\r
- \r
- return;\r
- }\r
- \r
- // added event selection for pA\r
- \r
- if(fSystem == pA){\r
- \r
- if(fUtils->IsFirstEventInChunk(aodEvent)) {\r
- AliInfo("Rejecting the event - first in the chunk");\r
- ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(9);\r
- return;\r
- }\r
- if(!fUtils->IsVertexSelected2013pA(aodEvent)) {\r
- ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(10);\r
- AliInfo("Rejecting the event - bad vertex");\r
- return;\r
- }\r
- }\r
- // ******************************** END event selections **************************************************\r
- \r
- AliCentrality *centralityObj = 0;\r
- Double_t MultipOrCent = -1;\r
- \r
- if(fUseCentrality){\r
- /* if(fSystem == AA ){ */ centralityObj = aodEvent->GetHeader()->GetCentralityP();\r
- MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
- //AliInfo(Form("Centrality is %f", MultipOrCent));\r
- }\r
- \r
- if(!fUseCentrality) MultipOrCent = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.);\r
- \r
- \r
- fCorrelator->SetAODEvent(aodEvent); // set the event to be processed\r
- \r
- ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1);\r
- \r
- Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event is in pool settings\r
- if(!correlatorON) {\r
- AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");\r
- return;\r
- }\r
- \r
- if(fmontecarlo) fCorrelator->SetMCArray(fmcArray);\r
- \r
- // check the event type\r
- // load MC header\r
- if(fmontecarlo){\r
- AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));\r
- if (fmontecarlo && !mcHeader) {\r
- AliError("Could not find MC Header in AOD");\r
- return;\r
- }\r
- \r
- Bool_t isMCeventgood = kFALSE;\r
- \r
- \r
- Int_t eventType = mcHeader->GetEventType();\r
- Int_t NMCevents = fCuts2->GetNofMCEventType();\r
- \r
- for(Int_t k=0; k<NMCevents; k++){\r
- Int_t * MCEventType = fCuts2->GetMCEventType();\r
- \r
- if(eventType == MCEventType[k]) isMCeventgood= kTRUE;\r
- ((TH1D*)fOutputMC->FindObject("EventTypeMC"))->Fill(eventType);\r
- }\r
- \r
- if(NMCevents && !isMCeventgood){\r
- if(fDebugLevel) std::cout << "The MC event " << eventType << " not interesting for this analysis: skipping" << std::endl;\r
- return;\r
- }\r
- \r
- } // end if montecarlo\r
- \r
- \r
- // checks on vertex and multiplicity of the event\r
- AliAODVertex *vtx = aodEvent->GetPrimaryVertex();\r
- Double_t zVtxPosition = vtx->GetZ(); // zvertex\r
- \r
- \r
- if(fFullmode) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);\r
- \r
- \r
- \r
- // D* reconstruction\r
- TClonesArray *arrayDStartoD0pi=0;\r
- if(!aodEvent && AODEvent() && IsStandardAOD()) {\r
- // In case there is an AOD handler writing a standard AOD, use the AOD\r
- // event in memory rather than the input (ESD) event.\r
- aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());\r
- // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
- // have to taken from the AOD event hold by the AliAODExtension\r
- AliAODHandler* aodHandler = (AliAODHandler*)\r
- ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
- if(aodHandler->GetExtensions()) {\r
- AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
- AliAODEvent *aodFromExt = ext->GetAOD();\r
- arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");\r
- }\r
- } else {\r
- arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");\r
- }\r
- \r
- if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r
- \r
- \r
- // get the poolbin\r
- \r
- Int_t poolbin = fCuts2->GetPoolBin(MultipOrCent, zVtxPosition);\r
- \r
- \r
- // initialize variables you will need for the D*\r
- \r
- Double_t ptDStar;//\r
- Double_t phiDStar;//\r
- Double_t etaDStar;//\r
- Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag;\r
- Double_t invMassDZero;\r
- Double_t deltainvMDStar;\r
- \r
- \r
- Double_t mPDGD0=1.8648;//TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
- Double_t mPDGDstar=2.01022;//TDatabasePDG::Instance()->GetParticle(413)->Mass();\r
- \r
- \r
- //MC tagging for DStar\r
- //D* and D0 prongs needed to MatchToMC method\r
- Int_t pdgDgDStartoD0pi[2]={421,211};\r
- Int_t pdgDgD0toKpi[2]={321,211};\r
- \r
- Bool_t isDStarCand = kFALSE;\r
- Bool_t isDfromB = kFALSE;\r
- Bool_t isEventMixingFilledPeak = kFALSE;\r
- Bool_t isEventMixingFilledSB = kFALSE;\r
- Bool_t EventHasDStarCandidate = kFALSE;\r
- Bool_t EventHasDZeroSideBandCandidate = kFALSE;\r
- Bool_t EventHasDStarSideBandCandidate = kFALSE;\r
- //loop on D* candidates\r
- \r
- Int_t looponDCands = 0;\r
- if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast();\r
- if(!fReco) looponDCands = fmcArray->GetEntriesFast();\r
- \r
- Int_t nOfDStarCandidates = 0;\r
- Int_t nOfSBCandidates = 0;\r
- \r
- Double_t DmesonEfficiency = 1.;\r
- Double_t DmesonWeight = 1.;\r
- Double_t efficiencyvariable = -999;\r
- \r
- \r
- \r
- for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {\r
- isInPeak = kFALSE;\r
- isInDStarSideBand = kFALSE;\r
- isInDZeroSideBand = kFALSE;\r
- isDStarMCtag = kFALSE;\r
- isDfromB = kFALSE;\r
- ptDStar = -123.4;\r
- phiDStar = -999;\r
- etaDStar = -56.;\r
- invMassDZero = - 999;\r
- deltainvMDStar = -998;\r
- AliAODRecoCascadeHF* dstarD0pi;\r
- AliAODRecoDecayHF2Prong* theD0particle;\r
- AliAODMCParticle* DStarMC=0x0;\r
- Short_t daughtercharge = -2;\r
- Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info\r
- Int_t trackiddaugh1 = -1;\r
- Int_t trackidsoftPi = -1;\r
- \r
- // start the if reconstructed candidates from here ************************************************\r
- \r
- if(fReco){//// if reconstruction is applied\r
- dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);\r
- if(!dstarD0pi->GetSecondaryVtx()) continue;\r
- theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();\r
- if (!theD0particle) continue;\r
- \r
- \r
- // track quality cuts\r
- Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks\r
- // region of interest + topological cuts + PID\r
- Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected\r
- \r
- //apply track selections\r
- if(!isTkSelected) continue;\r
- if(!isSelected) continue;\r
- if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;\r
- \r
- \r
- ptDStar = dstarD0pi->Pt();\r
- phiDStar = dstarD0pi->Phi();\r
- etaDStar = dstarD0pi->Eta();\r
- if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;\r
- if(fEfficiencyVariable == kMult || fEfficiencyVariable == kCentr) efficiencyvariable = MultipOrCent;\r
- if(fEfficiencyVariable == kEta) efficiencyvariable = etaDStar;\r
- if(fEfficiencyVariable == kRapidity) efficiencyvariable = dstarD0pi->YDstar();\r
- if(fEfficiencyVariable == kNone) efficiencyvariable = 0;\r
- \r
- // get the D meson efficiency\r
- DmesonEfficiency = fCuts2->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);\r
- \r
- \r
-\r
- if(fUseDmesonEfficiencyCorrection){\r
- if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;\r
- else {// THIS ELSE STATEMENT MUST BE REFINED: THE EFFICIENCY MAP HAS TO BE REPLACED WITH A WEIGHT MAP COOKED A PRIORI\r
- 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\r
- else DmesonWeight = 1.e+5; // at low pt it can be that the efficiency is really low\r
- }\r
- }\r
- else DmesonWeight = 1.; \r
- \r
- // continue;\r
- \r
- Int_t mcLabelDStar = -999;\r
- if(fmontecarlo){\r
- // find associated MC particle for D* ->D0toKpi\r
- mcLabelDStar = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,fmcArray/*,kFALSE*/);\r
- if(mcLabelDStar>=0) isDStarMCtag = kTRUE;\r
- if(!isDStarMCtag) continue;\r
- AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);\r
- //check if DStar from B\r
- Int_t labelMother = MCDStar->GetMother();\r
- AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));\r
- if(!mother) continue;\r
- Int_t motherPDG =TMath::Abs(mother->PdgCode());\r
- if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;\r
- }\r
- \r
- \r
- phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar);\r
- \r
- // set the phi of the D meson in the correct range\r
- \r
- Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());\r
- \r
- \r
- Double_t dmDStarWindow = 0.0019/3;// 0.0019 = 3 sigma\r
- // Double_t mD0Window=0.074/3;\r
- \r
- Double_t mD0Window= fD0Window[ptbin]/3;\r
- //cout << "Check with new window " << fD0Window[ptbin]/3 << endl;\r
- \r
- \r
- \r
- invMassDZero = dstarD0pi->InvMassD0();\r
- if(!fmixing && fFullmode) ((TH2F*)fOutput->FindObject("D0InvMass"))->Fill(ptDStar,invMassDZero);\r
- \r
- deltainvMDStar = dstarD0pi->DeltaInvMass();\r
- \r
- //good D0 candidates\r
- if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){\r
- \r
- if(!fmixing) ((TH2F*)fOutput->FindObject("DeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);\r
- // good D*?\r
- if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow){\r
- \r
- if(!fmixing) ((TH1F*)fOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight);\r
- if(!fmixing) ((TH2F*)fOutput->FindObject("PhiEtaTrigger"))->Fill(phiDStar,etaDStar);\r
- isInPeak = kTRUE;\r
- EventHasDStarCandidate = kTRUE;\r
- nOfDStarCandidates++;\r
- } // end Good D*\r
- \r
- // D* sideband?\r
- if((deltainvMDStar-(mPDGDstar-mPDGD0)>fDMesonSigmas[2]*dmDStarWindow) && (deltainvMDStar-(mPDGDstar-mPDGD0)<fDMesonSigmas[3]*dmDStarWindow)){\r
- isInDStarSideBand = kTRUE;\r
- EventHasDStarSideBandCandidate = kTRUE;\r
- } // end D* sideband\r
- \r
- }// end good D0 candidates\r
- \r
- //D0 sidebands\r
- if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){\r
- if(!fmixing)((TH2F*)fOutput->FindObject("bkgDeltaInvMass"))->Fill(ptDStar,deltainvMDStar,DmesonWeight);\r
- if(!fmixing && fFullmode)((TH2F*)fOutput->FindObject("D0InvMassinSB"))->Fill(ptDStar,invMassDZero,DmesonWeight);\r
- \r
- if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0] *dmDStarWindow){ // is in DStar peak region?\r
- if(!fmixing) ((TH1F*)fOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight);\r
- isInDZeroSideBand = kTRUE;\r
- EventHasDZeroSideBandCandidate = kTRUE;\r
- nOfSBCandidates++;\r
- if(!fmixing) ((TH2F*)fOutput->FindObject("PhiEtaSideBand"))->Fill(phiDStar,etaDStar);\r
- }\r
- \r
- }//end if sidebands\r
- \r
- \r
- \r
- \r
- if(!isInPeak && !isInDStarSideBand && !isInDZeroSideBand) continue; // skip if it is not side band or peak event - SAVE CPU TIME\r
- \r
- \r
- // check properties of the events containing the D*\r
-\r
- \r
- \r
- isDStarCand = kTRUE;\r
- \r
- // charge of the daughter od the\r
- daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();\r
- \r
- \r
- trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();\r
- trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();\r
- trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();\r
- \r
- // end here the reco\r
- \r
- \r
- }// end of if for applied reconstruction to D*\r
- \r
- Int_t DStarLabel = -1;\r
- \r
- if(!fReco){ // use pure MC information\r
- \r
- // get the DStar Particle\r
- DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));\r
- if (!DStarMC) {\r
- AliWarning("Careful: DStar MC Particle not found in tree, skipping");\r
- continue;\r
- }\r
- DStarLabel = DStarMC->GetLabel();\r
- if(DStarLabel>0)isDStarMCtag = kTRUE;\r
- \r
- Int_t PDG =TMath::Abs(DStarMC->PdgCode());\r
- if(PDG !=413) continue; // skip if it is not a DStar\r
- // check fiducial acceptance\r
- if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;\r
- \r
- //check if DStar from B\r
- Int_t labelMother = DStarMC->GetMother();\r
- AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));\r
- if(!mother) continue;\r
- Int_t motherPDG =TMath::Abs(mother->PdgCode());\r
- if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;\r
- \r
- Bool_t isDZero = kFALSE;\r
- Bool_t isSoftPi = kFALSE;\r
- \r
- if(fUseHadronicChannelAtKineLevel){\r
- //check decay channel on MC ************************************************\r
- Int_t NDaugh = DStarMC->GetNDaughters();\r
- if(NDaugh != 2) continue; // skip decay channels w/0 2 prongs\r
- \r
- for(Int_t i=0; i<NDaugh;i++){ // loop on daughters\r
- Int_t daugh_label = DStarMC->GetDaughter(i);\r
- AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(daugh_label));\r
- if(!mcDaughter) continue;\r
- Int_t daugh_pdg = TMath::Abs(mcDaughter->GetPdgCode());\r
- if(fDebugLevel) std::cout << "Daughter " << i << " pdg code is " << daugh_pdg << std::endl;\r
- \r
- if(daugh_pdg == 421) {isDZero = kTRUE;\r
- Int_t NDaughD0 = mcDaughter->GetNDaughters();\r
- if(NDaughD0 != 2) continue; // skip decay channels w/0 2 prongs\r
- trackiddaugh0 = mcDaughter->GetDaughter(0);\r
- trackiddaugh1 = mcDaughter->GetDaughter(1);\r
- Bool_t isKaon = kFALSE;\r
- Bool_t isPion = kFALSE;\r
- \r
- for(Int_t k=0;k<NDaughD0;k++){\r
- Int_t labelD0daugh = mcDaughter->GetDaughter(k);\r
- AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelD0daugh));\r
- if(!mcGrandDaughter) continue;\r
- Int_t granddaugh_pdg = TMath::Abs(mcGrandDaughter->GetPdgCode());\r
- if(granddaugh_pdg==321) isKaon = kTRUE;\r
- if(granddaugh_pdg==211) isPion = kTRUE;\r
- }\r
- if(!isKaon || !isKaon) continue; // skip if not correct decay channel of D0\r
- }\r
- \r
- if(daugh_pdg == 211) {\r
- isSoftPi = kTRUE;\r
- daughtercharge = mcDaughter->Charge();\r
- trackidsoftPi = daugh_label;}\r
- }\r
- if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel\r
- } // end check decay channel\r
- \r
- ptDStar = DStarMC->Pt();\r
- phiDStar = DStarMC->Phi();\r
- etaDStar = DStarMC->Eta();\r
- \r
- if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;\r
- \r
- } // end use pure MC information\r
- \r
- \r
- // getting the number of triggers in the MCtag D* case\r
- if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutput->FindObject("MCtagPtDStar"))->Fill(ptDStar);\r
- if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar);\r
- if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar);\r
- \r
- \r
- fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters\r
- fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);\r
- \r
- \r
- // ************************************************ CORRELATION ANALYSIS STARTS HERE\r
- \r
- \r
- Bool_t execPool = fCorrelator->ProcessEventPool();\r
- \r
- \r
- if(fmixing && !execPool) {\r
- AliInfo("Mixed event analysis: pool is not ready");\r
- if(!isEventMixingFilledPeak && isInPeak) {\r
- ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(1);\r
- isEventMixingFilledPeak = kTRUE;\r
- }\r
- if (!isEventMixingFilledSB && isInDZeroSideBand) {\r
- ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(3);\r
- isEventMixingFilledSB=kTRUE;\r
- }\r
- \r
- continue;\r
- }\r
- \r
- // check event topology\r
- if(fmixing&&execPool){\r
- // pool is ready - run checks on bins filling\r
- if(!isEventMixingFilledPeak && isInPeak) {\r
- ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(0);\r
- if(fFullmode) EventMixingChecks(aodEvent);\r
- isEventMixingFilledPeak = kTRUE;\r
- }\r
- \r
- if(!isEventMixingFilledSB && isInDZeroSideBand) {\r
- ((TH1D*)fOutput->FindObject("CheckPoolReadiness"))->Fill(2);\r
- isEventMixingFilledSB=kTRUE;\r
- }\r
- }\r
- \r
- Int_t NofEventsinPool = 1;\r
- if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();\r
- \r
- \r
- for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){// loop on events in the pool; if it is SE analysis, stops at one\r
- \r
- Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix);\r
- if(!analyzetracks) {\r
- AliInfo("AliHFCorrelator::Cannot process the track array");\r
- continue;\r
- }\r
- \r
- //initialization of variables for correlations with leading particles\r
- Double_t DeltaPhiLeading = -999.;\r
- Double_t DeltaEtaLeading = -999.;\r
- \r
- \r
- Int_t NofTracks = fCorrelator->GetNofTracks();\r
- \r
- \r
- if(isInPeak && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInPeak"))->Fill(NofTracks);\r
- if(isInDZeroSideBand && fFullmode) ((TH1D*)fOutput->FindObject("NofTracksInSB"))->Fill(NofTracks);\r
- \r
- \r
- \r
- Double_t arraytofill[5];\r
- Double_t MCarraytofill[7];\r
- \r
- \r
- Double_t weight;\r
- \r
- for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates\r
- Bool_t runcorrelation = fCorrelator->Correlate(iTrack);\r
- if(!runcorrelation) continue;\r
- \r
- Double_t DeltaPhi = fCorrelator->GetDeltaPhi();\r
- Double_t DeltaEta = fCorrelator->GetDeltaEta();\r
- \r
- AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();\r
- if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}\r
- \r
- Double_t ptHad = hadron->Pt();\r
- Double_t phiHad = hadron->Phi();\r
- Double_t etaHad = hadron->Eta();\r
- Int_t label = hadron->GetLabel();\r
- Int_t trackid = hadron->GetID();\r
- Double_t efficiency = hadron->GetWeight();\r
- \r
- weight = 1;\r
- if(fUseEfficiencyCorrection && efficiency){\r
- weight = DmesonWeight * (1./efficiency);\r
- }\r
- \r
- phiHad = fCorrelator->SetCorrectPhiRange(phiHad);\r
- \r
- \r
- if(fFullmode) ((TH2F*)fOutput->FindObject("WeightChecks"))->Fill(ptHad,efficiency);\r
- \r
- arraytofill[0] = DeltaPhi;\r
- arraytofill[1] = DeltaEta;\r
- arraytofill[2] = ptDStar;\r
- arraytofill[3] = ptHad;\r
- arraytofill[4] = poolbin;\r
- \r
- \r
- MCarraytofill[0] = DeltaPhi;\r
- MCarraytofill[1] = DeltaEta;\r
- MCarraytofill[2] = ptDStar;\r
- MCarraytofill[3] = ptHad;\r
- MCarraytofill[4] = poolbin;\r
- \r
- \r
- if(fmontecarlo){\r
- if(label<0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(0.,NofTracks);\r
- if(label>=0 && fFullmode) ((TH2D*)fOutputMC->FindObject("TrackLabels"))->Fill(1.,NofTracks);\r
- if(label<0) continue; // skip track with wrong label\r
- }\r
- \r
- Bool_t isDdaughter = kFALSE;\r
- // skip the D daughters in the correlation\r
- if(!fmixing && fReco){\r
- if(trackid == trackiddaugh0) continue;\r
- if(trackid == trackiddaugh1) continue;\r
- if(trackid == trackidsoftPi) continue;\r
- }\r
- \r
- if(!fmixing && !fReco){\r
- AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);\r
- if(!part) continue;\r
- if(IsDDaughter(DStarMC, part)) continue;\r
- cout << "Skipping DStar daugheter " << endl;\r
- }\r
- if(!fmixing && !fReco && fmontecarlo){ // skip D* Daughetrs if it is Pure MCDStar\r
- Int_t hadronlabel = label;\r
- for(Int_t k=0; k<4;k++){ // go back 4 generations and check the mothers\r
- if(DStarLabel<0){ break;}\r
- if(hadronlabel<0) { break;}\r
- AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fmcArray->At(hadronlabel));\r
- if(!mcParticle) {AliInfo("NO MC PARTICLE"); break;}\r
- hadronlabel = mcParticle->GetMother();\r
- if(hadronlabel == DStarLabel) isDdaughter = kTRUE;\r
- }\r
- \r
- if(isDdaughter && fDebugLevel){\r
- std::cout << "It is the D* daughter with label " << label << std::endl;\r
- std::cout << "Daughter 0 label = " << trackiddaugh0 << std::endl;\r
- std::cout << "Daughter 1 label = " << trackiddaugh1 << std::endl;\r
- std::cout << "Soft pi label = " << trackidsoftPi << std::endl;\r
- }\r
- \r
- if(isDdaughter) continue; // skip if track is from DStar\r
- }\r
- \r
- // ================ FILL CORRELATION HISTOGRAMS ===============================\r
- \r
- // monte carlo case (mc tagged D*)\r
- if((fmontecarlo && isDStarMCtag) || (fmontecarlo && !fReco)){ // check correlations of MC tagged DStars in MonteCarlo\r
- \r
- Bool_t* PartSource = fCuts2->IsMCpartFromHF(label,fmcArray); // check source of associated particle (hadron/kaon/K0)\r
- \r
- MCarraytofill[5] = 0;\r
- if(PartSource[0]) MCarraytofill[5] = 1;\r
- if(PartSource[1]) MCarraytofill[5] = 2;\r
- if(PartSource[2]&&PartSource[0]) MCarraytofill[5] = 3;\r
- if(PartSource[2]&&PartSource[1]) MCarraytofill[5] = 4;\r
- if(PartSource[3]) MCarraytofill[5] = 5;\r
- if(!isDfromB) MCarraytofill[6] = 0;\r
- if(isDfromB) MCarraytofill[6] = 1;\r
- if(!fReco && TMath::Abs(etaHad)>0.8) {\r
- delete [] PartSource;\r
- continue; // makes sure you study the correlation on MC truth only if particles are in acceptance\r
- }\r
- ((THnSparseF*)fOutputMC->FindObject("MCDStarCorrelationsDStarHadron"))->Fill(MCarraytofill);\r
- \r
- delete[] PartSource;\r
- }\r
- \r
- // Good DStar canidates\r
- if(isInPeak) {\r
- \r
- if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance\r
- if(fselect==1) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
- if(fselect==2) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
- if(fselect==3) ((THnSparseF*)fOutput->FindObject("CorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
- \r
- ((TH3F*)fOutput->FindObject("PhiEtaPart"))->Fill(phiHad,etaHad,MultipOrCent);\r
- if(fFullmode)((TH1D*)fOutput->FindObject("TracksInPeakSpectra"))->Fill(ptHad);\r
- \r
- }\r
- \r
- // Sidebands from D0 candidate\r
- if(isInDZeroSideBand) {\r
- \r
- if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance\r
- if(fselect==1) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
- if(fselect==2) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
- if(fselect==3) ((THnSparseF*)fOutput->FindObject("DZeroBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
- \r
- if(fFullmode) ((TH1D*)fOutput->FindObject("TracksInSBSpectra"))->Fill(ptHad);\r
- \r
- }\r
- \r
- // Sidebands from D* candidate\r
- if(isInDStarSideBand) {\r
- \r
- if(!fReco && TMath::Abs(etaHad)>0.8) continue; // makes sure you study the correlation on MC truth only if particles are in acceptance\r
- if(fselect==1 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarHadron"))->Fill(arraytofill,weight);\r
- if(fselect==2 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKaon"))->Fill(arraytofill,weight);\r
- if(fselect==3 && fFullmode) ((THnSparseF*)fOutput->FindObject("DStarBkgCorrelationsDStarKZero"))->Fill(arraytofill,weight);\r
-\r
- }\r
- \r
- \r
- } // end loop on track candidates\r
- \r
- \r
- \r
- // fill the leading particle histograms\r
- \r
- if(isInPeak && fFullmode) ((TH3D*)fOutput->FindObject("LeadingCand"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);\r
- if(isInDZeroSideBand && fFullmode) ((TH3D*)fOutput->FindObject("LeadingSB"))->Fill(DeltaPhiLeading,ptDStar,DeltaEtaLeading);\r
- \r
- } // end loop on events in the pool\r
- \r
- }// end loop on D* candidates\r
- \r
- \r
- \r
- // check events with D* or SB canidates\r
- if(fFullmode && EventHasDStarCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDStar"))->Fill(MultipOrCent,zVtxPosition);\r
- if(fFullmode && EventHasDZeroSideBandCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDZeroSB"))->Fill(MultipOrCent,zVtxPosition);\r
- \r
- if(fFullmode && EventHasDStarCandidate) ((TH2F*)fOutput->FindObject("EventPropsCheckifDStarSB"))->Fill(MultipOrCent,zVtxPosition);\r
- \r
- \r
- if(fFullmode) ((TH2F*)fOutput->FindObject("DStarCandidates"))->Fill(nOfDStarCandidates,MultipOrCent);\r
- if(fFullmode) ((TH2F*)fOutput->FindObject("SBCandidates"))->Fill(nOfSBCandidates,MultipOrCent);\r
- \r
- // update event pool\r
- Bool_t updated = fCorrelator->PoolUpdate();\r
- \r
- // if(updated) EventMixingChecks(aodEvent);\r
- if(!updated) AliInfo("Pool was not updated");\r
- \r
- \r
-} //end the exec\r
-\r
-//________________________________________ terminate ___________________________\r
-void AliAnalysisTaskDStarCorrelations::Terminate(Option_t*)\r
-{ \r
- // The Terminate() function is the last function to be called during\r
- // a query. It always runs on the client, it can be used to present\r
- // the results graphically or save the results to file.\r
- \r
- AliAnalysisTaskSE::Terminate();\r
- \r
- fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
- if (!fOutput) { \r
- printf("ERROR: fOutput not available\n");\r
- return;\r
- }\r
-\r
- return;\r
-}\r
-//_____________________________________________________\r
-Bool_t AliAnalysisTaskDStarCorrelations::IsDDaughter(AliAODMCParticle* d, AliAODMCParticle* track) const {\r
- \r
- //Daughter removal in MCKine case\r
- Bool_t isDaughter = kFALSE;\r
- Int_t labelD0 = d->GetLabel();\r
- \r
- Int_t mother = track->GetMother();\r
- \r
- //Loop on the mothers to find the D0 label (it must be the trigger D0, not a generic D0!)\r
- while (mother > 0){\r
- AliAODMCParticle* mcMoth = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!\r
- if (mcMoth){\r
- if (mcMoth->GetLabel() == labelD0) isDaughter = kTRUE;\r
- mother = mcMoth->GetMother(); //goes back by one\r
- } else{\r
- AliError("Failed casting the mother particle!");\r
- break;\r
- }\r
- }\r
- \r
- return isDaughter;\r
-}\r
-\r
-//_____________________________________________________\r
-void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){\r
- \r
- //cout << "DEFINING THNSPARSES "<< endl;\r
- \r
- Double_t Pi = TMath::Pi();\r
- Int_t nbinscorr = fPhiBins;\r
- Double_t lowcorrbin = -0.5*Pi;\r
- Double_t upcorrbin = 1.5*Pi;\r
- // define the THnSparseF\r
- \r
- //sparse bins\r
- \r
- //1 delta_phi\r
- //2 delta_eta\r
- //3 D* pt\r
- //4 multiplicity\r
- //5 track pt\r
- //6 zVtx position\r
- \r
- Int_t nbinsPool = (fCuts2->GetNZvtxPoolBins())*(fCuts2->GetNCentPoolBins());\r
- \r
- \r
- Int_t nbinsSparse[5]= {nbinscorr, 32,30,250,nbinsPool};\r
- Double_t binLowLimitSparse[5]={lowcorrbin,-1.6, 0, 0,-0.5};\r
- Double_t binUpLimitSparse[5]= {upcorrbin, 1.6,30,25,nbinsPool-0.5};\r
- \r
- Int_t MCnbinsSparse[7]= {nbinscorr, 32,30,250,nbinsPool,10,2}; \r
- Double_t MCbinLowLimitSparse[7]={lowcorrbin,-1.6, 0, 0,-0.5,-0.5,-0.5}; // \r
- Double_t MCbinUpLimitSparse[7]= {upcorrbin, 1.6,30,25,nbinsPool-0.5,9.5,1.5};\r
- \r
- TString sparsename = "CorrelationsDStar";\r
- if(fselect==1) sparsename += "Hadron";\r
- if(fselect==2) sparsename += "Kaon";\r
- if(fselect==3) sparsename += "KZero";\r
- \r
- TString D0Bkgsparsename = "DZeroBkg";\r
- D0Bkgsparsename += sparsename;\r
- \r
- TString DStarBkgsparsename = "DStarBkg";\r
- DStarBkgsparsename += sparsename;\r
- \r
- TString MCSparseName = "MCDStar";\r
- MCSparseName += sparsename;\r
- // signal correlations\r
- THnSparseF * Correlations = new THnSparseF(sparsename.Data(),"Correlations for signal",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
- \r
- // bkg correlations from D0 sidebands\r
- THnSparseF * DZeroBkgCorrelations = new THnSparseF(D0Bkgsparsename.Data(),"Bkg Correlations estimated with D0 sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
- \r
- // bkg correlations from D* sidebands\r
- THnSparseF * DStarBkgCorrelations = new THnSparseF(DStarBkgsparsename.Data(),"Bkg Correlations estimated with D* sidebands",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);\r
- \r
- \r
- THnSparseF * MCCorrelations = new THnSparseF(MCSparseName.Data(),"MC Correlations",7,MCnbinsSparse,MCbinLowLimitSparse,MCbinUpLimitSparse);\r
- \r
- MCCorrelations->GetAxis(5)->SetBinLabel(1," All ");\r
- MCCorrelations->GetAxis(5)->SetBinLabel(2," from hadron Heavy flavour");\r
- MCCorrelations->GetAxis(5)->SetBinLabel(3," from c->D");\r
- MCCorrelations->GetAxis(5)->SetBinLabel(4," from b->D");\r
- MCCorrelations->GetAxis(5)->SetBinLabel(5," from b->B");\r
- MCCorrelations->GetAxis(5)->SetBinLabel(6," from quark Heavy flavour");\r
- MCCorrelations->GetAxis(5)->SetBinLabel(7," from c");\r
- MCCorrelations->GetAxis(5)->SetBinLabel(8," from b");\r
- \r
- MCCorrelations->GetAxis(6)->SetBinLabel(1," if D* from c");\r
- MCCorrelations->GetAxis(6)->SetBinLabel(2," if D* from b");\r
- \r
- Correlations->Sumw2();\r
- DZeroBkgCorrelations->Sumw2();\r
- DStarBkgCorrelations->Sumw2();\r
- \r
- fOutput->Add(Correlations);\r
- fOutput->Add(DZeroBkgCorrelations);\r
- if(fFullmode) fOutput->Add(DStarBkgCorrelations);\r
- if(fmontecarlo) fOutputMC->Add(MCCorrelations);\r
- \r
- \r
-}\r
-//__________________________________________________________________________________________________\r
-void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){\r
- \r
- Double_t Pi = TMath::Pi();\r
- Int_t nbinscorr = fPhiBins;\r
- Double_t lowcorrbin = -0.5*Pi ; // shift the bin by half the width so that at 0 is it the bin center\r
- Double_t upcorrbin = 1.5*Pi ;\r
- \r
- // ========================= histograms for both Data and MonteCarlo\r
- \r
- \r
- TH1D * NofEvents = new TH1D("NofEvents","NofEvents",12,-0.5,11.5);\r
- NofEvents->GetXaxis()->SetBinLabel(1," All events");\r
- NofEvents->GetXaxis()->SetBinLabel(2," Selected events");\r
- NofEvents->GetXaxis()->SetBinLabel(3," Rejected - SPD Pileup");\r
- NofEvents->GetXaxis()->SetBinLabel(4," Rejected - Centrality");\r
- NofEvents->GetXaxis()->SetBinLabel(5," Rejected - No Reco Vtx");\r
- NofEvents->GetXaxis()->SetBinLabel(6," Rejected - Vtx Contr.");\r
- NofEvents->GetXaxis()->SetBinLabel(7," Rejected - Vtx outside fid.acc.");\r
- NofEvents->GetXaxis()->SetBinLabel(8," Rejected - Trigger");\r
- NofEvents->GetXaxis()->SetBinLabel(9," Rejected - Phys.Sel");\r
- NofEvents->GetXaxis()->SetBinLabel(10," Rejected - pA - 1st in chunk");\r
- NofEvents->GetXaxis()->SetBinLabel(11," Rejected - pA - bad vtx");\r
- fOutput->Add(NofEvents);\r
- \r
- \r
- \r
- \r
- TH2F *D0InvMass = new TH2F("D0InvMass","K#pi invariant mass distribution",300,0,30,1500,0.5,3.5);\r
- if(!fmixing && fFullmode) fOutput->Add(D0InvMass);\r
- \r
- TH2F *D0InvMassinSB = new TH2F("D0InvMassinSB","K#pi invariant mass distribution in sb",300,0,30,1500,0.5,3.5);\r
- if(!fmixing && fFullmode) fOutput->Add(D0InvMassinSB);\r
- \r
- //TH2F *DeltaInvMass = new TH2F("DeltaInvMass","K#pi#pi - K#pi invariant mass distribution",300,0,30,750,0.1,0.2);\r
- //if(!fmixing) fOutput->Add(DeltaInvMass);\r
- 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);\r
- if(!fmixing) fOutput->Add(DeltaInvMass);\r
- \r
- 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);\r
- if(!fmixing) fOutput->Add(bkgDeltaInvMass);\r
- \r
- DeltaInvMass->Sumw2();\r
- bkgDeltaInvMass->Sumw2();\r
- \r
- TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",50,0,50);\r
- if(!fmixing) fOutput->Add(RecoPtDStar);\r
- \r
- TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",50,0,50);\r
- if(!fmixing) fOutput->Add(RecoPtBkg);\r
- \r
- TH1D *MCtagPtDStarfromCharm = new TH1D("MCtagPtDStarfromCharm","RECO pt of MCtagged DStars from charm",50,0,50);\r
- if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromCharm);\r
- \r
- TH1D *MCtagPtDStarfromBeauty = new TH1D("MCtagPtDStarfromBeauty","RECO pt of MCtagged DStars from beauty",50,0,50);\r
- if(fmontecarlo) fOutputMC->Add(MCtagPtDStarfromBeauty);\r
- \r
- TH1F *MCtagPtDStar = new TH1F("MCtagPtDStar","RECO pt of MCtagged DStars side bands",50,0,50);\r
- if(!fmixing) fOutput->Add(MCtagPtDStar);\r
- \r
- TH2F *KZeroSpectra = new TH2F("KZeroSpectra","Spectra of K0s",500,0.3,0.8,250,0,25);\r
- if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectra);\r
- \r
- TH2F *KZeroSpectraifHF = new TH2F("KZeroSpectraifHF","Spectra of K0s in association with a D*",500,0.3,0.8,250,0,25);\r
- if(fselect==3 && fFullmode) fOutput->Add(KZeroSpectraifHF);\r
- \r
- TH1D * NofTracksInPeak = new TH1D("NofTracksInPeak","N of associated tracks per D trigger; Nof tracks; Entries",500,-0.5,499.5);\r
- if(fFullmode) fOutput->Add(NofTracksInPeak);\r
- \r
- TH1D * NofTracksInSB = new TH1D("NofTracksInSB","N of associated tracks per SideBand trigger; Nof tracks; Entries",500,-0.5,499.5);\r
- if(fFullmode) fOutput->Add(NofTracksInSB);\r
- \r
- TH1D * TracksInPeakSpectra = new TH1D("TracksInPeakSpectra","Pt Spectra tracks with D trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);\r
- if(fFullmode)fOutput->Add(TracksInPeakSpectra);\r
- \r
- TH1D * TracksInSBSpectra = new TH1D("TracksInSBSpectra","Pt Spectra tracks with SideBand trigger; p_{T} GeV/c; Entries",500,-0.5,49.5);\r
- if(fFullmode)fOutput->Add(TracksInSBSpectra);\r
- \r
- \r
- //TH2I * EventMixingCheck = new TH2I("EventMixingCheck","EventMixingCheck",5,-0.5,4.5,7,-0.5,6.5);\r
- //if(fmixing) fOutput->Add(EventMixingCheck);\r
- \r
- \r
- TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
- if(fFullmode)fOutput->Add(EventPropsCheck);\r
- \r
- TH2F * EventPropsCheckifDStar = new TH2F("EventPropsCheckifDStar","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
- if(fFullmode)fOutput->Add(EventPropsCheckifDStar);\r
- \r
- TH2F * EventPropsCheckifDZeroSB = new TH2F("EventPropsCheckifDZeroSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
- if(fFullmode)fOutput->Add(EventPropsCheckifDZeroSB);\r
- \r
- TH2F * EventPropsCheckifDStarSB = new TH2F("EventPropsCheckifDStarSB","Properties of the event with D* Cand; Multiplicity; ZVtx Position [cm]",1000,0,1000,40,-10,10);\r
- if(fFullmode)fOutput->Add(EventPropsCheckifDStarSB);\r
- \r
- \r
- TH2F * WeightChecks = new TH2F("WeightChecks","Checks on efficiency correction",300,0,30,100,0.005,1.005);\r
- if(fFullmode)fOutput->Add(WeightChecks);\r
- \r
- \r
- \r
- TH2F * PhiEtaTrigger = new TH2F("PhiEtaTrigger","#phi distribution of the trigger particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);\r
- fOutput->Add(PhiEtaTrigger);\r
- \r
- TH2F * PhiEtaSideBand = new TH2F("PhiEtaSideBand","#phi distribution of the sideband particle",nbinscorr,lowcorrbin,upcorrbin,18,-0.9,0.9);\r
- fOutput->Add(PhiEtaSideBand);\r
- \r
- 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);\r
- fOutput->Add(PhiEtaPart);\r
- \r
- TH2F * DStarCandidates = new TH2F("DStarCandidates","# of D* candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);\r
- if(fFullmode)fOutput->Add(DStarCandidates);\r
- \r
- TH2F * SBCandidates = new TH2F("SBCandidates","# of SB candidates per event vs multiplicity",6,-0.5,5.5,50,0,500);\r
- if(fFullmode)fOutput->Add(SBCandidates);\r
- \r
- \r
- //correlations histograms\r
- TString histoname1 = "DPhiDStar";\r
- if(fselect==1) histoname1 += "Hadron";\r
- if(fselect==2) histoname1 += "Kaon";\r
- if(fselect==3) histoname1 += "KZero";\r
- \r
- /*\r
- TH3D * DPhiDStar = new TH3D(histoname1.Data(),histoname1.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TH3D * DPhiDStarKZero1 = new TH3D("DPhiDStarKZero1","DPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- //side band background histograms\r
- TString histoname2 = "bkg";\r
- histoname2 += histoname1;\r
- TH3D * bkgDPhiDStar = new TH3D(histoname2.Data(),histoname2.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- TH3D * bkgDPhiDStarKZero1 = new TH3D("bkgDPhiDStarKZero1","bkgDPhiDStarKZero1",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- \r
- fOutput->Add(DPhiDStar);\r
- \r
- if(fselect==3){fOutput->Add(DPhiDStarKZero1);}\r
- \r
- fOutput->Add(bkgDPhiDStar);\r
- \r
- if(fselect==3){fOutput->Add(bkgDPhiDStarKZero1);}\r
- \r
- */\r
- // leading particle\r
- TH3D * leadingcand = new TH3D("LeadingCand","LeadingCand",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- TH3D * leadingsidebands = new TH3D("LeadingSB","LeadingSB",nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- if(fFullmode)fOutput->Add(leadingcand);\r
- if(fFullmode)fOutput->Add(leadingsidebands);\r
- \r
- // ========================= histos for analysis on MC only\r
- \r
- TH1D * EventTypeMC = new TH1D("EventTypeMC","EventTypeMC",100,-0.5,99.5);\r
- if(fmontecarlo) fOutputMC->Add(EventTypeMC);\r
- \r
- TH1F * MCSources = new TH1F("MCSources","Origin of associated particles in MC", 10, -0.5, 9.5);\r
- MCSources->GetXaxis()->SetBinLabel(1," All ");\r
- MCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");\r
- MCSources->GetXaxis()->SetBinLabel(3," from c->D");\r
- MCSources->GetXaxis()->SetBinLabel(4," from b->D");\r
- MCSources->GetXaxis()->SetBinLabel(5," from b->B");\r
- MCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");\r
- MCSources->GetXaxis()->SetBinLabel(7," from c");\r
- MCSources->GetXaxis()->SetBinLabel(8," from b");\r
- \r
- if(fmontecarlo) fOutputMC->Add(MCSources);\r
- \r
- // leading particle from mc source\r
- TH1F * LeadingMCSources = new TH1F("LeadingMCSources","Origin of associated leading particles in MC", 10, -0.5, 9.5);\r
- LeadingMCSources->GetXaxis()->SetBinLabel(1," All ");\r
- LeadingMCSources->GetXaxis()->SetBinLabel(2," from hadron Heavy flavour");\r
- LeadingMCSources->GetXaxis()->SetBinLabel(3," from c->D");\r
- LeadingMCSources->GetXaxis()->SetBinLabel(4," from b->D");\r
- LeadingMCSources->GetXaxis()->SetBinLabel(5," from b->B");\r
- LeadingMCSources->GetXaxis()->SetBinLabel(6," from quark Heavy flavour");\r
- LeadingMCSources->GetXaxis()->SetBinLabel(7," from c");\r
- LeadingMCSources->GetXaxis()->SetBinLabel(8," from b");\r
- \r
- if(fmontecarlo && fFullmode) fOutputMC->Add(LeadingMCSources);\r
- \r
- // all hadrons\r
- TString histoname3 = "MCTag";\r
- histoname3 += histoname1;\r
- TH3D * MCTagDPhiDStar = new TH3D(histoname3.Data(),histoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString histoname44 = "CharmDOrigin";\r
- histoname44 += histoname1;\r
- histoname44 += "MC";\r
- \r
- TH3D * CharmDOriginDPhiDStar = new TH3D(histoname44.Data(),histoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- \r
- TString histoname54 = "BeautyDOrigin";\r
- histoname54 += histoname1;\r
- histoname54 += "MC";\r
- TH3D * BeautyDOriginDPhiDStar = new TH3D(histoname54.Data(),histoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString histoname55 = "BeautyBOrigin";\r
- histoname55 += histoname1;\r
- histoname55 += "MC";\r
- TH3D * BeautyBOriginDPhiDStar = new TH3D(histoname55.Data(),histoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString histoname4 = "CharmQuarkOrigin";\r
- histoname4 += histoname1;\r
- histoname4 += "MC";\r
- TH3D * CharmQuarkOriginDPhiDStar = new TH3D(histoname4.Data(),histoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString histoname5 = "BeautyQuarkOrigin";\r
- histoname5 += histoname1;\r
- histoname5 += "MC";\r
- TH3D * BeautyQuarkOriginDPhiDStar = new TH3D(histoname5.Data(),histoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString histoname6 = "NonHFOrigin";\r
- histoname6 += histoname1;\r
- histoname6 += "MC";\r
- TH3D * NonHFOriginDPhiDStar = new TH3D(histoname6.Data(),histoname6.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- if(fmontecarlo && fFullmode){\r
- \r
- fOutputMC->Add(MCTagDPhiDStar);\r
- fOutputMC->Add(CharmDOriginDPhiDStar);\r
- fOutputMC->Add(BeautyDOriginDPhiDStar);\r
- fOutputMC->Add(BeautyBOriginDPhiDStar);\r
- fOutputMC->Add(CharmQuarkOriginDPhiDStar);\r
- fOutputMC->Add(BeautyQuarkOriginDPhiDStar);\r
- fOutputMC->Add(NonHFOriginDPhiDStar);\r
- \r
- }\r
- \r
- // ========================= histos for analysis on MC\r
- // all leading hadron\r
- TString Leadinghistoname3 = "LeadingMCTag";\r
- Leadinghistoname3 += histoname1;\r
- TH3D * LeadingMCTagDPhiDStar = new TH3D(Leadinghistoname3.Data(),Leadinghistoname3.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString Leadinghistoname44 = "LeadingCharmDOrigin";\r
- Leadinghistoname44 += histoname1;\r
- Leadinghistoname44 += "MC";\r
- \r
- TH3D * LeadingCharmDOriginDPhiDStar = new TH3D(Leadinghistoname44.Data(),Leadinghistoname44.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- \r
- TString Leadinghistoname54 = "LeadingBeautyDOrigin";\r
- Leadinghistoname54 += histoname1;\r
- Leadinghistoname54 += "MC";\r
- TH3D * LeadingBeautyDOriginDPhiDStar = new TH3D(Leadinghistoname54.Data(),Leadinghistoname54.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString Leadinghistoname55 = "LeadingBeautyBOrigin";\r
- Leadinghistoname55 += histoname1;\r
- Leadinghistoname55 += "MC";\r
- TH3D * LeadingBeautyBOriginDPhiDStar = new TH3D(Leadinghistoname55.Data(),Leadinghistoname55.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString Leadinghistoname4 = "LeadingCharmQuarkOrigin";\r
- Leadinghistoname4 += histoname1;\r
- Leadinghistoname4 += "MC";\r
- TH3D * LeadingCharmQuarkOriginDPhiDStar = new TH3D(Leadinghistoname4.Data(),Leadinghistoname4.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- TString Leadinghistoname5 = "LeadingBeautyQuarkOrigin";\r
- Leadinghistoname5 += histoname1;\r
- Leadinghistoname5 += "MC";\r
- TH3D * LeadingBeautyQuarkOriginDPhiDStar = new TH3D(Leadinghistoname5.Data(),Leadinghistoname5.Data(),nbinscorr,lowcorrbin,upcorrbin,50,0,50,39,-2,2);\r
- \r
- \r
- \r
- \r
- if(fmontecarlo && fFullmode){\r
- \r
- fOutputMC->Add(LeadingMCTagDPhiDStar);\r
- fOutputMC->Add(LeadingCharmDOriginDPhiDStar);\r
- fOutputMC->Add(LeadingBeautyDOriginDPhiDStar);\r
- fOutputMC->Add(LeadingBeautyBOriginDPhiDStar);\r
- fOutputMC->Add(LeadingCharmQuarkOriginDPhiDStar);\r
- fOutputMC->Add(LeadingBeautyQuarkOriginDPhiDStar);\r
- \r
- }\r
- \r
- TH3F * MCPhiEtaPart = new TH3F("MCPhiEtaPart","#phi distribution of the associated particle",nbinscorr,lowcorrbin,upcorrbin,50,-2.5,2.5,6,-0.5,6.5);\r
- MCPhiEtaPart->GetZaxis()->SetBinLabel(1,"All particles");\r
- MCPhiEtaPart->GetZaxis()->SetBinLabel(2,"from c quark");\r
- MCPhiEtaPart->GetZaxis()->SetBinLabel(3,"from b quark");\r
- MCPhiEtaPart->GetZaxis()->SetBinLabel(4,"from D from c");\r
- MCPhiEtaPart->GetZaxis()->SetBinLabel(5,"from D from b");\r
- MCPhiEtaPart->GetZaxis()->SetBinLabel(6,"from B from b");\r
- if(fmontecarlo) fOutputMC->Add(MCPhiEtaPart);\r
- \r
- TH2D * TrackLabels = new TH2D("TrackLabels","NofEvents;track label; multiplicity",2,-0.5,1.5,500,-0.5,499.5);\r
- if(fmontecarlo && fFullmode) fOutputMC->Add(TrackLabels);\r
- \r
- // ============================= EVENT MIXING CHECKS ======================================\r
- \r
- Int_t MaxNofEvents = fCuts2->GetMaxNEventsInPool();\r
- Int_t MinNofTracks = fCuts2->GetMinNTracksInPool();\r
- Int_t NofCentBins = fCuts2->GetNCentPoolBins();\r
- Double_t * CentBins = fCuts2->GetCentPoolBins();\r
- Int_t NofZVrtxBins = fCuts2->GetNZvtxPoolBins();\r
- Double_t *ZVrtxBins = fCuts2->GetZvtxPoolBins();\r
- \r
- \r
- \r
- Int_t k =0;\r
- \r
- if(fSystem == AA) k = 100; // PbPb centrality\r
- if(fSystem == pp || fSystem == pA) k = NofCentBins; // pp multiplicity\r
- \r
- \r
- //Double_t minvalue = CentBins[0];\r
- //Double_t maxvalue = CentBins[NofCentBins+1];\r
- //Double_t Zminvalue = ZVrtxBins[0];\r
- //Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1];\r
- \r
- Double_t minvalue, maxvalue;\r
- Double_t Zminvalue, Zmaxvalue;\r
- \r
- Zminvalue = -15.;\r
- Zmaxvalue = 15;\r
- if(fSystem == AA) {minvalue = 0; maxvalue = 100;} // PbPb\r
- if(fSystem == pp || fSystem == pA) {minvalue = 0; maxvalue = 500;} // multilpicity\r
- \r
- //Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};\r
- // Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};\r
- // Double_t * events = Nevents;\r
- Double_t eventsv[] ={0,1000000};\r
- //Double_t * events = new Double_t[2];\r
- // events[0] = 0;\r
-// events[1] = 1000000;\r
- Double_t *events = eventsv;\r
- Int_t Nevents = 1000000;\r
- // TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,Nevents,events);\r
- \r
- TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,Nevents,events[0],events[1]);\r
- \r
- EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
- EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");\r
- EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin");\r
- if(fmixing && fFullmode) fOutput->Add(EventsPerPoolBin);\r
- \r
- Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks;\r
- //Int_t Diff = MaxNofTracks-MinNofTracks;\r
- \r
- //Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks};\r
- // Double_t * trackN = Ntracks;\r
- \r
- TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,minvalue,maxvalue,NofZVrtxBins,-15,15,MaxNofTracks,0,MaxNofTracks);\r
- NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
- NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");\r
- NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin");\r
- \r
- if(fmixing && fFullmode) fOutput->Add(NofTracksPerPoolBin);\r
- \r
- TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Calls per pool bin",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins);\r
- NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
- NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");\r
- if(fmixing && fFullmode) fOutput->Add(NofPoolBinCalls);\r
- \r
- \r
- \r
- TH2D * EventProps = new TH2D("EventProps","Event properties",100,minvalue,maxvalue,100,Zminvalue,Zmaxvalue);\r
- EventProps->GetXaxis()->SetTitle("Centrality/multiplicity ");\r
- EventProps->GetYaxis()->SetTitle("Z vertex [cm]");\r
- if(fmixing && fFullmode) fOutput->Add(EventProps);\r
- \r
- TH1D * CheckPoolReadiness = new TH1D("CheckPoolReadiness","Pool readiness",5,-0.5,4.5);\r
- CheckPoolReadiness->GetXaxis()->SetBinLabel(1,"Have a D cand, pool is ready");\r
- CheckPoolReadiness->GetXaxis()->SetBinLabel(2,"Have a D cand, pool is not ready");\r
- CheckPoolReadiness->GetXaxis()->SetBinLabel(3,"Have a SB cand, pool is ready");\r
- CheckPoolReadiness->GetXaxis()->SetBinLabel(4,"Have a SB cand, pool is not ready");\r
- \r
- if(fmixing) fOutput->Add(CheckPoolReadiness);\r
- \r
- \r
-}\r
-\r
-//__________________________________________________________________________________________________\r
-void AliAnalysisTaskDStarCorrelations::EnlargeDZeroMassWindow(){\r
- \r
-\r
- //Float_t* ptbins = fCuts->GetPtBinLimits();\r
- if(fD0Window) delete fD0Window;\r
- fD0Window = new Float_t[fNofPtBins];\r
- \r
- AliInfo("Enlarging the D0 mass windows from cut object\n"); \r
- Int_t nvars = fCuts->GetNVars();\r
-\r
- if(nvars<1){\r
- AliWarning("EnlargeDZeroMassWindow: 0 variables in cut object... check!");\r
- return;\r
- }\r
- Float_t** rdcutsvalmine=new Float_t*[nvars];\r
- for(Int_t iv=0;iv<nvars;iv++){\r
- rdcutsvalmine[iv]=new Float_t[fNofPtBins];\r
- }\r
- \r
- \r
- for (Int_t k=0;k<nvars;k++){\r
- for (Int_t j=0;j<fNofPtBins;j++){\r
- \r
- // enlarge D0 window\r
- if(k==0) {\r
- fD0Window[j] =fCuts->GetCutValue(0,j);\r
- rdcutsvalmine[k][j] = 5.* fCuts->GetCutValue(0,j);\r
- cout << "the set window = " << fD0Window[j] << " for ptbin " << j << endl;\r
- }\r
- else rdcutsvalmine[k][j] =fCuts->GetCutValue(k,j);\r
- \r
- // set same windows\r
- //rdcutsvalmine[k][j] =oldCuts->GetCutValue(k,j);\r
- }\r
- }\r
- \r
- fCuts->SetCuts(nvars,fNofPtBins,rdcutsvalmine);\r
- \r
- AliInfo("\n New windows set\n"); \r
- fCuts->PrintAll();\r
- \r
- \r
- for(Int_t iv=0;iv<nvars;iv++){\r
- delete rdcutsvalmine[iv];\r
- }\r
- delete [] rdcutsvalmine;\r
- \r
-}\r
-\r
-\r
-//____________________________ Run checks on event mixing ___________________________________________________\r
-void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){\r
- \r
- \r
- AliCentrality *centralityObj = 0;\r
- Int_t multiplicity = -1;\r
- Double_t MultipOrCent = -1;\r
- \r
- // get the pool for event mixing\r
- if(fSystem != AA){ // pp\r
- multiplicity = AOD->GetNTracks();\r
- MultipOrCent = multiplicity; // convert from Int_t to Double_t\r
- }\r
- if(fSystem == AA){ // PbPb\r
- \r
- centralityObj = AOD->GetHeader()->GetCentralityP();\r
- MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");\r
- AliInfo(Form("Centrality is %f", MultipOrCent));\r
- }\r
- \r
- AliAODVertex *vtx = AOD->GetPrimaryVertex();\r
- Double_t zvertex = vtx->GetZ(); // zvertex\r
- \r
- \r
- \r
- \r
- AliEventPool * pool = fCorrelator->GetPool();\r
- \r
- \r
- \r
- \r
- ((TH2D*)fOutput->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool\r
- ((TH2D*)fOutput->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties\r
- \r
- ((TH3D*)fOutput->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of events in the pool\r
- ((TH3D*)fOutput->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of calls of pool\r
-}\r
- \r
-\r
-\r
-\r
-\r
+/**************************************************************************
+ * 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: AliAnalysisTaskDStarCorrelations.cxx 65139 2013-11-25 14:47:45Z fprino $ */
+
+//#include <TDatabasePDG.h>
+#include <TParticle.h>
+#include <TVector3.h>
+#include <TChain.h>
+#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),
+fBkgMethod(kDZeroSB),
+fReco(kTRUE),
+fUseEfficiencyCorrection(kFALSE),
+fUseDmesonEfficiencyCorrection(kFALSE),
+fUseCentrality(kFALSE),
+fUseHadronicChannelAtKineLevel(kFALSE),
+fRemoveMoreThanOneDmesonCandidate(kFALSE),
+fPhiBins(32),
+fEvents(0),
+fDebugLevel(0),
+fDisplacement(0),
+fDim(0),
+fNofPtBins(0),
+fMaxEtaDStar(0.9),
+fDMesonSigmas(0),
+
+fD0Window(0),
+
+fOutput(0x0),
+fDmesonOutput(0x0),
+fTracksOutput(0x0),
+fEMOutput(0x0),
+fCorrelationOutput(0x0),
+fOutputMC(0x0),
+
+fCuts(0),
+fAssocCuts(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),
+fBkgMethod(kDZeroSB),
+fReco(kTRUE),
+fUseEfficiencyCorrection(kFALSE),
+fUseDmesonEfficiencyCorrection(kFALSE),
+fUseCentrality(kFALSE),
+fUseHadronicChannelAtKineLevel(kFALSE),
+fRemoveMoreThanOneDmesonCandidate(kFALSE),
+fPhiBins(32),
+fEvents(0),
+fDebugLevel(0),
+fDisplacement(0),
+fDim(0),
+fNofPtBins(0),
+fMaxEtaDStar(0.9),
+fDMesonSigmas(0),
+fD0Window(0),
+
+fOutput(0x0),
+fDmesonOutput(0x0),
+fTracksOutput(0x0),
+fEMOutput(0x0),
+fCorrelationOutput(0x0),
+fOutputMC(0x0),
+
+fCuts(0),
+fAssocCuts(AsscCuts),
+fUtils(0),
+fTracklets(0),
+fDeffMapvsPt(0),
+fDeffMapvsPtvsMult(0),
+fDeffMapvsPtvsEta(0)
+{
+ Info("AliAnalysisTaskDStarCorrelations","Calling Constructor");
+
+ SetDim();
+ if(fSystem == AA) fUseCentrality = kTRUE; else fUseCentrality = kFALSE;
+
+ // if(fSystem == AA) fBkgMethod = kDStarSB; else fBkgMethod = kDZeroSB;
+
+ 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,TList::Class()); // histos from data and MC
+ DefineOutput(4,TList::Class()); // histos from MC
+ DefineOutput(5,TList::Class()); // histos from data and MC
+ DefineOutput(6,TList::Class()); // histos from data and MC
+ DefineOutput(7,AliNormalizationCounter::Class()); // normalization
+
+ DefineOutput(8,AliRDHFCutsDStartoKpipi::Class()); // my D meson cuts
+ DefineOutput(9,AliHFAssociatedTrackCuts::Class()); // my associated tracks cuts
+}
+
+//__________________________________________________________________________
+
+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(fAssocCuts) {delete fAssocCuts; fAssocCuts=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(8,copyfCuts);
+
+ // Post the hadron cuts
+ PostData(9,fAssocCuts);
+
+ 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();
+
+ fDmesonOutput = new TList();
+ fDmesonOutput->SetOwner();
+
+ fTracksOutput = new TList();
+ fTracksOutput->SetOwner();
+
+ fEMOutput = new TList();
+ fEMOutput->SetOwner();
+
+ fCorrelationOutput = new TList();
+ fCorrelationOutput->SetOwner();
+
+ // define histograms
+ DefineHistoForAnalysis();
+ DefineThNSparseForAnalysis();
+
+
+
+
+ fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(7)->GetContainer()->GetName()));
+ fCounter->Init();
+
+ Double_t Pi = TMath::Pi();
+ //AliHFCorrelator(const Char_t* name, AliHFAssociatedTrackCuts *cuts, Bool_t useCentrality, AliRDHFCuts *cutObject)
+ fCorrelator = new AliHFCorrelator("Correlator",fAssocCuts,fUseCentrality,fCuts); // fAssocCuts 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,fDmesonOutput); // set the outputs
+ PostData(3,fTracksOutput); // set the outputs
+ PostData(4,fEMOutput); // set the outputs
+ PostData(5,fCorrelationOutput); // set the outputs
+ PostData(6,fOutputMC); // set the outputs
+ PostData(7,fCounter); // set the outputs
+}
+
+//________________________________________ user exec ____________
+void AliAnalysisTaskDStarCorrelations::UserExec(Option_t *){
+ // cout << "Task debug check 1 " << endl;
+ // ********************************************** LOAD THE EVENT ****************************************************
+
+ if (!fInputEvent) {
+ Error("UserExec","NO EVENT FOUND!");
+ return;
+ }
+
+ AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
+ if(!aodEvent){
+ AliError("AOD event not found!");
+ return;
+ }
+
+ fEvents++; // event counter
+ ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(0);
+
+ fCounter->StoreEvent(aodEvent,fCuts,fmontecarlo); // store event in AliNormalizationCounter
+
+ // load MC array
+ fmcArray = dynamic_cast<TClonesArray*>(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); /*cout << "Reject PILEUP" << endl;*/}
+ if(fCuts->IsEventRejectedDueToCentrality()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(3); /*cout << "Reject CENTRALITY" << endl;*/}
+ if(fCuts->IsEventRejectedDueToNotRecoVertex()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(4); /*cout << "Reject NOT RECO VTX" << endl;*/}
+ if(fCuts->IsEventRejectedDueToVertexContributors()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(5); /*cout << "Reject VTX CONTRIB" << endl;*/}
+ if(fCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(6); /*cout << "Reject VTX no fid reg " << endl;*/}
+ if(fCuts->IsEventRejectedDueToTrigger()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(7); /*cout << "Reject TRIGGER" << endl;*/}
+ if(fCuts->IsEventRejectedDuePhysicsSelection()) {((TH1D*)fOutput->FindObject("NofEvents"))->Fill(8); /*cout << "Reject PHYS SEL" << endl;*/}
+
+ 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;
+ }
+ }
+
+ fCorrelator->SetAODEvent(aodEvent); // set the event in the correlator class
+
+ Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing and check if event within the pool definition
+ if(!correlatorON) {
+ AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
+ return; // if not the case, skips the event
+ }
+
+ ((TH1D*)fOutput->FindObject("NofEvents"))->Fill(1); // count events that are passing the event selection
+
+
+
+ // cout << "Task debug check 2 " << endl;
+ // ********************************************** CENTRALITY/MULTIPLICITY ****************************************************
+
+
+ Double_t MultipOrCent = fCorrelator->GetCentrality(); // returns centrality or multiplicity
+
+
+
+ // ********************************************** MC SETUP ****************************************************
+
+ if(fmontecarlo) {
+
+ fCorrelator->SetMCArray(fmcArray); // export MC array into correlatior class
+
+ AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+ if (fmontecarlo && !mcHeader) {
+ AliError("Could not find MC Header in AOD");
+ return;
+ }
+
+ // select MC event type (PP, GS etc) - those are defined in the associated tracks cut object
+ Bool_t isMCeventgood = kFALSE;
+
+
+ Int_t eventType = mcHeader->GetEventType();
+ Int_t NMCevents = fAssocCuts->GetNofMCEventType();
+
+ for(Int_t k=0; k<NMCevents; k++){
+ Int_t * MCEventType = fAssocCuts->GetMCEventType();
+
+ 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 fmontecarlo
+
+
+ // ********************************************** EVENT PROPERTIES ****************************************************
+ // checks on vertex and multiplicity of the event
+ AliAODVertex *vtx = aodEvent->GetPrimaryVertex();
+ Double_t zVtxPosition = vtx->GetZ(); // zvertex
+
+
+
+ if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
+
+ if(!fmixing) ((TH2F*)fOutput->FindObject("EventPropsCheck"))->Fill(MultipOrCent,zVtxPosition);
+ if(!fmixing) ((TH1D*)fOutput->FindObject("MultiplicityOnlyCheck"))->Fill(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.));
+
+
+ Int_t poolbin = fAssocCuts->GetPoolBin(MultipOrCent, zVtxPosition); // get the event pool bin - commented for the moment - to be checked
+
+
+
+ // ********************************************** D * selection ****************************************************
+
+
+ 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<AliAODEvent*> (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");
+ }
+
+
+
+ // D* related variables
+
+ Double_t ptDStar; // D* pt
+ Double_t phiDStar; // D* phi
+ Double_t etaDStar; // D* eta
+ Bool_t isInPeak, isInDZeroSideBand, isInDStarSideBand, isDStarMCtag; // flags for selection
+ Double_t invMassDZero; // D0 inv mass
+ Double_t deltainvMDStar; // D* - D0 invarian mass
+
+ 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;
+ Bool_t isTrackForPeakFilled = kFALSE;
+ Bool_t isTrackForSBFilled = kFALSE;
+ //loop on D* candidates
+
+ Int_t looponDCands = 0;
+ if(fReco) looponDCands = arrayDStartoD0pi->GetEntriesFast(); // load N of D* candidates from AOD array
+ if(!fReco) looponDCands = fmcArray->GetEntriesFast(); // load N of D* candidates from MC array
+
+ Int_t nOfDStarCandidates = 0; // D candidates counter
+ Int_t nOfSBCandidates = 0; // SB candidates counter
+
+ Double_t DmesonEfficiency = 1.; // Efficiency of D meson
+ Double_t DmesonWeight = 1.; // Applyed weight
+ Double_t efficiencyvariable = -999; // Variable to be used (defined by the AddTask)
+ Double_t dmDStarWindow = 0.0019/3; // DStar window
+
+ Int_t ncandidates = 0;
+
+ // cout << "Task debug check 3 " << endl;
+ // loop over D meson candidates
+ for (Int_t iDStartoD0pi = 0; iDStartoD0pi<looponDCands; iDStartoD0pi++) {
+
+ if(fRemoveMoreThanOneDmesonCandidate && ncandidates) continue; // if there is more than one D candidate, skip it
+
+ // initialize variables
+ isInPeak = kFALSE;
+ isInDStarSideBand = kFALSE;
+ isInDZeroSideBand = kFALSE;
+
+ EventHasDStarCandidate = kFALSE;
+ EventHasDZeroSideBandCandidate = kFALSE;
+ EventHasDStarSideBandCandidate = kFALSE;
+
+
+ isDStarMCtag = kFALSE;
+ isDfromB = kFALSE;
+ ptDStar = -123.4;
+ phiDStar = -999;
+ etaDStar = -56.;
+ invMassDZero = - 999;
+ deltainvMDStar = -998;
+ AliAODRecoCascadeHF* dstarD0pi;
+ AliAODRecoDecayHF2Prong* theD0particle;
+ AliAODMCParticle* DStarMC=0x0;
+ Short_t daughtercharge = -2;
+ Int_t trackiddaugh0 = -1; // track id if it is reconstruction - label if it is montecarlo info
+ Int_t trackiddaugh1 = -1;
+ Int_t trackidsoftPi = -1;
+ Int_t ptbin = -1;
+
+ // ============================================ using reconstruction on Data or MC
+ if(fReco){
+ // cout << "Task debug check 4 " << endl;
+ // get the D objects
+ dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
+ if(!dstarD0pi->GetSecondaryVtx()) continue;
+ theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
+ if (!theD0particle) continue;
+
+
+ // apply topological cuts
+
+ // 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 topological cuts
+ if(!isTkSelected) continue;
+ if(!isSelected) continue;
+ if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
+
+ // get D candidate variables
+ ptDStar = dstarD0pi->Pt();
+ if(ptDStar<fCuts->GetMinPtCandidate()||ptDStar>fCuts->GetMaxPtCandidate()) continue;
+
+ 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;
+
+
+ // if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue; // skip candidates outside the defined eta range
+
+ phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar); // set the Phi of the D* in the range defined a priori (-0.5 Pi - 1.5 Pi)
+ ptbin=fCuts->PtBin(dstarD0pi->Pt()); // get the pt bin of the D*
+
+ // cout << "DStar pt = " << ptDStar << endl;
+ // cout << "pt bin = " << ptbin << endl;
+ // if(ptbin<3) continue;
+
+ Double_t mD0Window= fD0Window[ptbin]/3;
+ invMassDZero = dstarD0pi->InvMassD0();
+ deltainvMDStar = dstarD0pi->DeltaInvMass();
+
+
+ // get the D meson efficiency
+ DmesonEfficiency = fAssocCuts->GetTrigWeight(dstarD0pi->Pt(),efficiencyvariable);
+
+ // check this!
+ if(fUseDmesonEfficiencyCorrection){
+ if(DmesonEfficiency>1.e-5) DmesonWeight = 1./DmesonEfficiency;
+ else DmesonWeight = 0; // do not consider te entry if the efficiency is too low
+ // else DmesonWeight = 1.;
+ }
+
+ // using montecarlo on reconstruction
+ 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;
+ if(isDStarMCtag){
+ AliAODMCParticle *MCDStar = (AliAODMCParticle*)fmcArray->At(mcLabelDStar);
+ //check if DStar from B
+ Int_t labelMother = MCDStar->GetMother();
+ AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelMother));
+ if(!mother) continue;
+ Int_t motherPDG =TMath::Abs(mother->PdgCode());
+ if((motherPDG>=500 && motherPDG <600) || (motherPDG>=5000 && motherPDG<6000 )) isDfromB = kTRUE;
+ }
+ }
+
+
+ // fill mass histograms
+ // cout << "crash here 1" << endl;
+ // plot D0 vs DStar mass
+ if(!fmixing){
+ // cout << Form("histDZerovsDStarMass_%d",ptbin) <<endl;
+ ((TH2D*)fDmesonOutput->FindObject(Form("histDZerovsDStarMass_%d",ptbin)))->Fill(invMassDZero,deltainvMDStar);
+ if(fUseDmesonEfficiencyCorrection) ((TH2D*)fDmesonOutput->FindObject(Form("histDZerovsDStarMassWeight_%d",ptbin)))->Fill(invMassDZero,deltainvMDStar,DmesonWeight);
+ }
+ // cout << "crash here 2" << endl;
+
+
+ // fill D0 invariant mass
+ if(!fmixing) {
+ ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMass_%d",ptbin)))->Fill(invMassDZero);
+ // cout << "Task debug check 5.1 " << endl;
+ if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDZeroMassWeight_%d",ptbin)))->Fill(invMassDZero,DmesonWeight);
+ } // end if !mixing
+
+
+ // Check for D* and D0 invariant candidates for phi and eta
+ if(fFullmode){
+ Double_t arrayDStar[5];
+
+ arrayDStar[0] = deltainvMDStar;
+ arrayDStar[1] = invMassDZero;
+ arrayDStar[2] = phiDStar;
+ arrayDStar[3] = etaDStar;
+ arrayDStar[4] = ptDStar;
+ if(!fmixing)((THnSparseF*)fDmesonOutput->FindObject("DmesonMassCheck"))->Fill(arrayDStar);
+ }
+
+
+ // good D0 canidates
+ if (TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[1]*mD0Window){
+ // fill D* invariant mass
+ if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMass_%d",ptbin)))->Fill(deltainvMDStar);
+ // fill D* invariant mass if weighting
+ if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);} // end if no mixing
+ isInPeak = kTRUE;
+ // fill other good candidate distributions - pt, phi eta
+ if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow) {
+ ((TH1F*)fDmesonOutput->FindObject("RecoPtDStar"))->Fill(ptDStar,DmesonWeight); // fill pt
+ if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiInclusiveDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
+ if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaInclusiveDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
+ nOfDStarCandidates++;
+ ncandidates ++;
+ EventHasDStarCandidate=kTRUE;
+ } // end if in good D* mass window
+
+ // count D* sideband candidates
+ if(fBkgMethod == kDStarSB ){
+ if(deltainvMDStar>0.1473 && deltainvMDStar<0.1644){
+ // if ((deltainvMDStar-(mPDGDstar-mPDGD0))>fDMesonSigmas[2]*dmDStarWindow && (deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[3]*dmDStarWindow ){
+ ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
+ if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
+ if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
+ nOfSBCandidates++;
+ ncandidates ++;
+ EventHasDStarSideBandCandidate = kTRUE;
+ }
+
+ } // end if using DStar sidebans
+
+
+ }// end good D0
+
+ // cout << "Task debug check 6 " << endl;
+ // Sideband D0
+ if (TMath::Abs(invMassDZero-mPDGD0)>fDMesonSigmas[2]*mD0Window && TMath::Abs(invMassDZero-mPDGD0)<fDMesonSigmas[3]*mD0Window ){
+ isInDZeroSideBand = kTRUE;
+ if(!fmixing){ ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMass_%d",ptbin)))->Fill(deltainvMDStar);
+ if(fUseDmesonEfficiencyCorrection) ((TH1F*)fDmesonOutput->FindObject(Form("histDStarFromSBMassWeight_%d",ptbin)))->Fill(deltainvMDStar,DmesonWeight);
+
+ if(fBkgMethod == kDZeroSB){
+ if(TMath::Abs(deltainvMDStar-(mPDGDstar-mPDGD0))<fDMesonSigmas[0]*dmDStarWindow) {
+
+ ((TH1F*)fDmesonOutput->FindObject("RecoPtBkg"))->Fill(ptDStar,DmesonWeight); // fill pt
+ if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("PhiSidebandDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
+ if(!fmixing) ((TH2F*)fDmesonOutput->FindObject("EtaSidebandDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
+ nOfSBCandidates++;
+ ncandidates ++;
+ EventHasDZeroSideBandCandidate = kTRUE;
+ }
+ }
+ }
+ }// end SideBand D0
+ // cout << "Task debug check 7 " << endl;
+
+ if(! isInPeak && !isInDZeroSideBand) continue; // skip candidates that are not interesting
+ if(TMath::Abs(deltainvMDStar)<0.142 || TMath::Abs(deltainvMDStar)>0.1644) continue; // skip all D* canidates outside the interesting mass ranges
+ // cout << "Good D* candidate" << endl;
+
+ // get charge of soft pion
+ daughtercharge = ((AliAODTrack*)dstarD0pi->GetBachelor())->Charge();
+
+ // get ID of daughters used for reconstruction
+ trackiddaugh0 = ((AliAODTrack*)theD0particle->GetDaughter(0))->GetID();
+ trackiddaugh1 = ((AliAODTrack*)theD0particle->GetDaughter(1))->GetID();
+ trackidsoftPi = ((AliAODTrack*)dstarD0pi->GetBachelor())->GetID();
+ }// end if reconstruction
+
+ // cout << "crash here 3" << endl;
+
+ // ============================================ using MC kinematics only
+ Int_t DStarLabel = -1;
+
+ if(!fReco){ // use pure MC information
+ // cout << "0 - Running on kine " << endl;
+ // get the DStar Particle
+ DStarMC = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iDStartoD0pi));
+ if (!DStarMC) {
+ AliWarning("Careful: DStar MC Particle not found in tree, skipping");
+ continue;
+ }
+ // cout << "1 - Have D* " << endl;
+ 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
+ // cout << "PDG = " << PDG << endl;
+ // check fiducial acceptance
+ if(!fCuts->IsInFiducialAcceptance(DStarMC->Pt(),DStarMC->Y())) continue;
+ // cout << "2 - Have D* in fiducial acceptance " << endl;
+
+ if(DStarMC->Pt()<fCuts->GetMinPtCandidate()||DStarMC->Pt()>fCuts->GetMaxPtCandidate()) continue;
+ ptbin=fCuts->PtBin(DStarMC->Pt());
+ // cout << "3 - Have D* in ptrange acceptance " << endl;
+ //check if DStar from B
+ Int_t labelMother = DStarMC->GetMother();
+ AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(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; i<NDaugh;i++){ // loop on DStar daughters
+ Int_t daugh_label = DStarMC->GetDaughter(i);
+ AliAODMCParticle* mcDaughter = dynamic_cast<AliAODMCParticle*>(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) {
+ 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;k<NDaughD0;k++){
+ Int_t labelD0daugh = mcDaughter->GetDaughter(k);
+ AliAODMCParticle* mcGrandDaughter = dynamic_cast<AliAODMCParticle*>(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 || !isPion) break; // skip if not correct decay channel of D0
+ isDZero = kTRUE;
+ } // end check on Dzero
+
+ if(daugh_pdg == 211) {
+ isSoftPi = kTRUE;
+ daughtercharge = mcDaughter->Charge();
+ trackidsoftPi = daugh_label;}
+ } // end loop on D* daughters
+ if(!isDZero || !isSoftPi) continue; // skip if not correct decay channel
+ } // end check decay channel
+
+ ptDStar = DStarMC->Pt();
+ phiDStar = DStarMC->Phi();
+ etaDStar = DStarMC->Eta();
+
+ phiDStar = fCorrelator->SetCorrectPhiRange(phiDStar);
+
+ if(TMath::Abs(etaDStar) > fMaxEtaDStar) continue;
+ // cout << "Dstars are selected" << endl;
+
+ }// end if pure MC information
+
+ // check kinematics for tagged D*
+ if(fmontecarlo && !fmixing) ((TH2F*)fOutputMC->FindObject("MCTagEtaInclusiveDStar"))->Fill(etaDStar,ptDStar); // fill phi, eta
+ if(fmontecarlo && !fmixing) ((TH2F*)fOutputMC->FindObject("MCTagPhiInclusiveDStar"))->Fill(phiDStar,ptDStar); // fill phi, eta
+
+
+ // getting the number of triggers in the MCtag D* case
+ if(fmontecarlo && isDStarMCtag) ((TH1F*)fOutputMC->FindObject("MCtagPtDStar"))->Fill(ptDStar,DmesonWeight);
+ if(fmontecarlo && isDStarMCtag && !isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromCharm"))->Fill(ptDStar,DmesonWeight);
+ if(fmontecarlo && isDStarMCtag && isDfromB) ((TH1D*)fOutputMC->FindObject("MCtagPtDStarfromBeauty"))->Fill(ptDStar,DmesonWeight);
+
+
+ fCorrelator->SetTriggerParticleProperties(ptDStar,phiDStar,etaDStar); // pass to the object the necessary trigger part parameters
+ fCorrelator->SetTriggerParticleDaughterCharge(daughtercharge);
+
+
+ // cout << "crash here 4" << endl;
+
+ // ************************************************ CORRELATION ANALYSIS STARTS HERE *****************************************************
+
+ // cout << "Correlating " << endl;
+
+ Bool_t execPool = fCorrelator->ProcessEventPool(); // checks pool readiness for mixed events
+
+ // if analysis is on mixed event and pool settings are not satisfied, fill relevant histograms and skip
+ if(fmixing && !execPool) {
+ AliInfo("Mixed event analysis: pool is not ready");
+ if(!isEventMixingFilledPeak && isInPeak) {
+ ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(1);
+ isEventMixingFilledPeak = kTRUE;
+ }
+ if (!isEventMixingFilledSB && isInDZeroSideBand) {
+ ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(3);
+ isEventMixingFilledSB=kTRUE;
+ }
+ continue;
+ } // end if pool not ok
+ // if analysis is on mixed event and pool settings are satisfied, fill relevant histograms and continue
+ if(fmixing&&execPool){
+ // pool is ready - run checks on bins filling
+ if(!isEventMixingFilledPeak && isInPeak) {
+ ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(0);
+ if(fFullmode) EventMixingChecks(aodEvent);
+ isEventMixingFilledPeak = kTRUE;
+ }
+
+ if(!isEventMixingFilledSB && isInDZeroSideBand) {
+ ((TH1D*)fEMOutput->FindObject("CheckPoolReadiness"))->Fill(2);
+ isEventMixingFilledSB=kTRUE;
+ }
+ } // end if pool ok
+
+
+
+
+ Int_t NofEventsinPool = 1;
+ if(fmixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
+
+ // cout << "crash here 5" << endl;
+ //************************************************** LOOP ON EVENTS IN EVENT POOL *****************************************************
+
+ 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); // process the associated tracks
+ if(!analyzetracks) {
+ AliInfo("AliHFCorrelator::Cannot process the track array");
+ continue;
+ }
+
+ Double_t arraytofill[5];
+ Double_t MCarraytofill[6]; // think about this
+ Double_t weight;
+
+ Int_t NofTracks = fCorrelator->GetNofTracks(); // number of tracks in event
+
+
+ if(!fmixing){
+ if(EventHasDStarCandidate) ((TH1D*)fTracksOutput->FindObject("TracksPerDcandidate"))->Fill(NofTracks);
+ if(fBkgMethod == kDStarSB && EventHasDStarSideBandCandidate) ((TH1D*)fTracksOutput->FindObject("TracksPerSBcandidate"))->Fill(NofTracks);
+ if(fBkgMethod == kDZeroSB && EventHasDZeroSideBandCandidate) ((TH1D*)fTracksOutput->FindObject("TracksPerSBcandidate"))->Fill(NofTracks);
+
+ if(isDStarMCtag && fmontecarlo) ((TH1D*)fTracksOutput->FindObject("TracksPerDMC"))->Fill(NofTracks);
+ }
+
+ //************************************************** LOOP ON TRACKS *****************************************************
+ // cout << "crash here 6" << endl;
+ for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ // looping on track candidates
+ // cout << "crash correlation 1" << endl;
+ Bool_t runcorrelation = fCorrelator->Correlate(iTrack); // calculate correlations
+ if(!runcorrelation) continue;
+
+ Double_t DeltaPhi = fCorrelator->GetDeltaPhi();
+ Double_t DeltaEta = fCorrelator->GetDeltaEta();
+
+ AliReducedParticle * hadron = fCorrelator->GetAssociatedParticle();
+ if(!hadron) {/*cout << "No Hadron" << endl;*/ continue;}
+
+
+ Int_t trackid = hadron->GetID();
+
+ // check D if it is a D meson daughter
+ if(!fmixing && fReco){ // for reconstruction
+ if(trackid == trackiddaugh0) continue;
+ if(trackid == trackiddaugh1) continue;
+ if(trackid == trackidsoftPi) continue;
+ }
+
+
+ Int_t label = hadron->GetLabel();
+
+
+ if(fmontecarlo){
+ AliAODMCParticle *part = (AliAODMCParticle*)fmcArray->At(label);
+ if(!part) continue;
+ if (!part->IsPhysicalPrimary()) continue; // removing secondary tracks
+ if(!fmixing && !fReco){
+ if(IsDDaughter(DStarMC, part)) continue; // skipping D* daughter
+ }
+ }
+
+
+ // if it is ok, then do the rest
+
+ Double_t ptHad = hadron->Pt();
+ Double_t phiHad = hadron->Phi();
+ phiHad = fCorrelator->SetCorrectPhiRange(phiHad); // set phi in correct range
+ Double_t etaHad = hadron->Eta();
+
+
+ Double_t efficiency = hadron->GetWeight();
+
+
+
+ Bool_t *trackOrigin = NULL;
+
+ if(fmontecarlo) {trackOrigin = fAssocCuts->IsMCpartFromHF(label,fmcArray);
+ if(trackOrigin[4]) {/*cout << "Something is wrong with hadron in MC - skipping" << endl; */continue;}
+ }
+ // cout << "crash correlation 3" << endl;
+
+ if(!isTrackForPeakFilled && !fmixing && EventHasDStarCandidate){
+
+ ((TH2F*)fTracksOutput->FindObject("PhiInclusiveTracks"))->Fill(phiHad,ptHad); // fill phi, eta
+ ((TH2F*)fTracksOutput->FindObject("EtaInclusiveTracks"))->Fill(etaHad,ptHad); // fill phi, eta
+ isTrackForPeakFilled = kTRUE; // makes sure you do not fill twice in case of more candidates
+ }
+
+ if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDZeroSB) && EventHasDZeroSideBandCandidate){
+ ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
+ ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
+ isTrackForSBFilled = kTRUE;
+ }
+
+ if(!isTrackForSBFilled && !fmixing && (fBkgMethod == kDStarSB) && EventHasDStarSideBandCandidate){
+ ((TH2F*)fTracksOutput->FindObject("PhiSidebandTracks"))->Fill(phiHad,ptHad); // fill phi, eta
+ ((TH2F*)fTracksOutput->FindObject("EtaSidebandTracks"))->Fill(etaHad,ptHad); // fill phi, eta
+ isTrackForSBFilled = kTRUE;
+ }
+ // cout << "crash correlation 4" << endl;
+
+ weight = 1;
+ if(fUseEfficiencyCorrection && efficiency){
+ weight = DmesonWeight * (1./efficiency);
+ }
+
+ // cout << "crash correlation 5" << endl;
+ arraytofill[0] = DeltaPhi;
+ arraytofill[1] = deltainvMDStar;
+ arraytofill[2] = DeltaEta;
+ arraytofill[3] = ptHad;
+ arraytofill[4] = poolbin;
+
+ if(fmontecarlo){
+ MCarraytofill[0] = DeltaPhi;
+ MCarraytofill[1] = deltainvMDStar;
+ MCarraytofill[2] = DeltaEta;
+ MCarraytofill[3] = ptHad;
+ MCarraytofill[4] = poolbin;
+
+ if(trackOrigin[0]) MCarraytofill[5] = 1 ; // is from charm
+ else if(trackOrigin[1]) MCarraytofill[5] = 2 ; // is from beauty
+ else MCarraytofill[5] = 0; // non HF track
+ }
+
+
+ // ============================================= FILL CORRELATION THNSparses ===============================
+
+ // filling signal
+ // if(isInPeak){
+ if(isInPeak){
+ // cout << "Filling signal " << endl;
+ // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
+ //cout ("CorrelationsDStarHadron_%d",ptbin)
+ if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarHadron_%d",ptbin)))->Fill(arraytofill,weight);
+ if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKaon_%d",ptbin)))->Fill(arraytofill,weight);
+ if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsDStarKZero_%d",ptbin)))->Fill(arraytofill,weight);
+ }
+
+ // filling bkg
+ // if(fBkgMethod == kDStarSB && isInPeak){ // bkg from DStar
+ if(fBkgMethod == kDStarSB && EventHasDStarSideBandCandidate){ // bkg from DStar
+ // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
+ // cout << "Filling bkg from D* sidebands " << endl;
+ if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
+ if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
+ if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDStarSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
+
+ } // end if bkg from DStar
+
+ // if(fBkgMethod == kDZeroSB && isInDZeroSideBand){ // bkg from DStar
+ if(fBkgMethod == kDZeroSB && EventHasDZeroSideBandCandidate){
+ // if(!fReco && TMath::Abs(etaHad)>0.8) continue;
+ // cout << "Filling bkg from Dzero sidebands " << endl;
+ if(fselect==1)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBHadron_%d",ptbin)))->Fill(arraytofill,weight);
+ if(fselect==2)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKaon_%d",ptbin)))->Fill(arraytofill,weight);
+ if(fselect==3)((THnSparseF*)fCorrelationOutput->FindObject(Form("CorrelationsBkgFromDZeroSBKZero_%d",ptbin)))->Fill(arraytofill,weight);
+
+ } // end if bkg from DZero
+
+
+
+ if(fmontecarlo){
+ // add the montecarlos
+ if(!fReco && TMath::Abs(etaHad)>0.8) continue;
+ if(!isDfromB){
+ //cout << "Filling correlations from charm " << endl;
+ // cout << "Ik zoek op " << Form("CorrelationsMCfromCharmHadron_%d",ptbin) << endl;
+ // cout << "de lijst bevat : " << endl;
+ fOutputMC->ls();
+ if(fselect==1)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmHadron_%d",ptbin)))->Fill(MCarraytofill,weight);
+ if(fselect==2)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmKaon_%d",ptbin)))->Fill(MCarraytofill,weight);
+ if(fselect==3)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromCharmKZero_%d",ptbin)))->Fill(MCarraytofill,weight);
+ }
+ if(isDfromB){
+ //cout << "Filling correlations from beauty " << endl;
+ if(fselect==1)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyHadron_%d",ptbin)))->Fill(MCarraytofill,weight);
+ if(fselect==2)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyKaon_%d",ptbin)))->Fill(MCarraytofill,weight);
+ if(fselect==3)((THnSparseF*)fOutputMC->FindObject(Form("CorrelationsMCfromBeautyKZero_%d",ptbin)))->Fill(MCarraytofill,weight);
+ }
+ }
+
+ delete [] trackOrigin;
+ } // end loop on associated tracks
+
+ } // end loop on events in event pool
+
+ if(fmixing){
+ ((TH1D*)fEMOutput->FindObject("PoolBinDistribution"))->Fill(poolbin);
+ ((TH2D*)fEMOutput->FindObject("EventDistributionPerPoolBin"))->Fill(poolbin,NofEventsinPool);
+ }
+ } // end loop on D mesons
+
+ if(!fmixing){
+ if(nOfDStarCandidates) ((TH1D*)fDmesonOutput->FindObject("NumberOfDCandidates"))->Fill(nOfDStarCandidates);
+ if(nOfSBCandidates) ((TH1D*)fDmesonOutput->FindObject("NumberOfSBCandidates"))->Fill(nOfSBCandidates);
+ }
+
+
+ Bool_t updated = fCorrelator->PoolUpdate(); // update event pool
+
+
+
+
+ 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<TList*> (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<AliAODMCParticle*>(fmcArray->At(mother)); //it's the mother of the track!
+ if (mcMoth){
+ if (mcMoth->GetLabel() == labelD0) {isDaughter = kTRUE; return isDaughter;}
+ mother = mcMoth->GetMother(); //goes back by one
+ } else{
+ AliError("Failed casting the mother particle!");
+ break;
+ }
+ }
+
+ return isDaughter;
+}
+
+//_____________________________________________________
+void AliAnalysisTaskDStarCorrelations::DefineThNSparseForAnalysis(){
+
+ Double_t Pi = TMath::Pi();
+ Int_t nbinscorr = fPhiBins;
+ Double_t lowcorrbin = -0.5*Pi ;
+ Double_t upcorrbin = 1.5*Pi ;
+
+
+ // create ThNSparses
+
+ Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
+
+
+ //sparse bins
+
+ //1 delta_phi
+ //2 invariant mass D *
+ //3 delta eta
+ //4 track pt
+
+
+ Int_t nbinsPool = (fAssocCuts->GetNZvtxPoolBins())*(fAssocCuts->GetNCentPoolBins());
+
+
+ // for reconstruction on Data and MC
+ Int_t nbinsSparse[5]= {nbinscorr , 32 , 32, 10,nbinsPool};
+ Double_t binLowLimitSparse[5]={lowcorrbin,0.14314 ,-1.6, 0,-0.5};
+ Double_t binUpLimitSparse[5]= {upcorrbin ,0.14794 , 1.6, 5,nbinsPool-0.5};
+
+ // Int_t nbinsSparseDStarSB[5]= {nbinscorr , 40 , 32, 10,nbinsPool};
+ // Double_t binLowLimitSparseDStarSB[5]={lowcorrbin,0.14788 ,-1.6, 0,-0.5};
+ // Double_t binUpLimitSparseDStarSB[5]= {upcorrbin ,0.1504 , 1.6, 5,nbinsPool-0.5};
+
+ Int_t nbinsSparseDStarSB[5]= {nbinscorr , 27 , 32, 10,nbinsPool};
+ Double_t binLowLimitSparseDStarSB[5]={lowcorrbin,0.1473 ,-1.6, 0,-0.5};
+ Double_t binUpLimitSparseDStarSB[5]= {upcorrbin ,0.1644 , 1.6, 5,nbinsPool-0.5};
+
+ Int_t nbinsSparseMC[6]= {nbinscorr , 32 , 32, 10,nbinsPool,3};
+ Double_t binLowLimitSparseMC[6]={lowcorrbin,0.14314 ,-1.6, 0,-0.5,-0.5};
+ Double_t binUpLimitSparseMC[6]= {upcorrbin ,0.14794 , 1.6, 5,nbinsPool-0.5,2.5};
+
+ TString signalSparseName = "";
+ TString bkgSparseName = "";
+ TString MCSparseNameCharm = "";
+ TString MCSparseNameBeauty = "";
+
+ THnSparseF * CorrelationsSignal = NULL;
+ THnSparseF * CorrelationsBkg = NULL;
+
+ THnSparseF * MCCorrelationsCharm = NULL;
+ THnSparseF * MCCorrelationsBeauty = NULL;
+
+
+ Float_t * ptbinlims = fCuts->GetPtBinLimits();
+
+
+
+
+ for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
+
+
+
+ if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
+
+
+
+ signalSparseName = "CorrelationsDStar";
+ if(fselect==1) signalSparseName += "Hadron_";
+ if(fselect==2) signalSparseName += "Kaon_";
+ if(fselect==3) signalSparseName += "KZero_";
+
+ bkgSparseName = "CorrelationsBkg";
+ if(fBkgMethod == kDStarSB ) bkgSparseName+="FromDStarSB";
+ if(fBkgMethod == kDZeroSB ) bkgSparseName+="FromDZeroSB";
+ if(fselect==1) bkgSparseName += "Hadron_";
+ if(fselect==2) bkgSparseName += "Kaon_";
+ if(fselect==3) bkgSparseName += "KZero_";
+
+ MCSparseNameCharm = "CorrelationsMCfromCharm";
+ if(fselect==1) MCSparseNameCharm += "Hadron_";
+ if(fselect==2) MCSparseNameCharm += "Kaon_";
+ if(fselect==3) MCSparseNameCharm += "KZero_";
+
+ MCSparseNameBeauty = "CorrelationsMCfromBeauty";
+ if(fselect==1) MCSparseNameBeauty += "Hadron_";
+ if(fselect==2) MCSparseNameBeauty += "Kaon_";
+ if(fselect==3) MCSparseNameBeauty += "KZero_";
+
+ signalSparseName+=Form("%d",iBin);
+ bkgSparseName+=Form("%d",iBin);
+ MCSparseNameCharm+=Form("%d",iBin);
+ MCSparseNameBeauty+=Form("%d",iBin);
+ //cout << "ThNSparses name = " << signalSparseName << endl;
+
+ // define thnsparses for signal candidates
+ CorrelationsSignal = new THnSparseF(signalSparseName.Data(),"Correlations for signal; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
+ CorrelationsSignal->Sumw2();
+ fCorrelationOutput->Add(CorrelationsSignal);
+
+ // define thnsparses for bkg candidates from DStar
+ if(fBkgMethod == kDStarSB ){
+ CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DStar; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",5,nbinsSparseDStarSB,binLowLimitSparseDStarSB,binUpLimitSparseDStarSB);
+ CorrelationsBkg->Sumw2();
+ fCorrelationOutput->Add(CorrelationsBkg);
+ }
+
+ // define thnsparses for bkg candidates from DZero
+ if(fBkgMethod == kDZeroSB ){
+ CorrelationsBkg = new THnSparseF(bkgSparseName.Data(),"Correlations for bkg from DZero; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",5,nbinsSparse,binLowLimitSparse,binUpLimitSparse);
+ CorrelationsBkg->Sumw2();
+ fCorrelationOutput->Add(CorrelationsBkg);
+ }
+
+ if(fmontecarlo){
+ MCCorrelationsCharm = new THnSparseF(MCSparseNameCharm.Data(),"Correlations for DStar from charm; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",6,nbinsSparseMC,binLowLimitSparseMC,binUpLimitSparseMC);
+ MCCorrelationsCharm->Sumw2();
+ fOutputMC->Add(MCCorrelationsCharm);
+
+ MCCorrelationsBeauty = new THnSparseF(MCSparseNameBeauty.Data(),"Correlations for DStar from beauty; #Delta#Phi; invariant mass; #Delta #eta;AssocTrk p_{T}",6,nbinsSparseMC,binLowLimitSparseMC,binUpLimitSparseMC);
+ MCCorrelationsBeauty->Sumw2();
+ fOutputMC->Add(MCCorrelationsBeauty);
+ }
+
+ } // end loop on bins
+
+ // make Dstar mass ThnSparse
+
+ Int_t NDstarSparse[5]= {30 , 15 , 32, 18,14};
+ Double_t binLowDstarSparse[5]={0.14,1.7 ,-0.5*TMath::Pi(), -0.9,2};
+ Double_t binUpDstarSparse[5]= {0.16,2.0 , 1.5*TMath::Pi(), 0.9,16};
+
+ THnSparseF * DmesonMassCheck = new THnSparseF("DmesonMassCheck","Check D meson mass; #Delta Inv Mass GeV/c^{2}; M(K#pi) GeV/c^{2}; #varphi; #eta; p_{T} GeV/c",5,NDstarSparse,binLowDstarSparse,binUpDstarSparse);
+ if(!fmixing) fDmesonOutput->Add(DmesonMassCheck);
+
+}
+
+//__________________________________________________________________________________________________
+void AliAnalysisTaskDStarCorrelations::DefineHistoForAnalysis(){
+
+ Double_t Pi = TMath::Pi();
+ Int_t nbinscorr = fPhiBins;
+ Double_t lowcorrbin = -0.5*Pi ;
+ Double_t upcorrbin = 1.5*Pi ;
+
+ // event counter
+ 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);
+
+ //event properties
+ TH2F * EventPropsCheck = new TH2F("EventPropsCheck","Properties of the event; Multiplicity/Centrality; ZVtx Position [cm]",1000,0,1000,40,-10,10);
+ fOutput->Add(EventPropsCheck);
+
+ //event properties
+ TH1D * SPDMultiplicty = new TH1D("MultiplicityOnlyCheck","Properties of the event; SPD Multiplicity",20000,0,20000);
+ fOutput->Add(SPDMultiplicty);
+
+
+ // =================================================== D* histograms
+ TH1F * D0mass = NULL;
+ TH1F * DStarMass = NULL;
+ TH1F * DStarFromSBMass = NULL;
+ TH2D * DZerovsDStarMass = NULL;
+
+ TH1F * D0massWeighted = NULL;
+ TH1F * DStarMassWeighted = NULL;
+ TH1F * DStarFromSBMassWeighted = NULL;
+ TH2D * DZerovsDStarMassWeighted = NULL;
+
+
+ TString nameDZeroMass = "", nameDStarMass = "", nameDStarFromSBMass = "", nameDZerovsDStarMass = "";
+
+ Int_t nofPtBins = fCuts->GetNPtBins();// number of ptbins
+ Float_t * ptbinlims = fCuts->GetPtBinLimits();
+
+ //GetMinPtCandidate()
+
+
+ for(Int_t iBin =0; iBin < nofPtBins; iBin++){ // create a mass histogram for each ptbin
+
+
+
+ if(ptbinlims[iBin]<fCuts->GetMinPtCandidate() || ptbinlims[iBin]>fCuts->GetMaxPtCandidate())continue;
+
+
+ // std::cout << ">>> Ptbin = " << iBin << " limit = " << ptbinlims[iBin] << std::endl;
+
+ nameDZeroMass = "histDZeroMass_";
+ nameDStarMass = "histDStarMass_";
+ nameDStarFromSBMass = "histDStarFromSBMass_";
+ nameDZerovsDStarMass = "histDZerovsDStarMass_";
+
+ nameDZeroMass+=Form("%d",iBin);
+ nameDStarMass+=Form("%d",iBin);
+ nameDStarFromSBMass+=Form("%d",iBin);
+ nameDZerovsDStarMass+=Form("%d",iBin);
+ // cout << "D vs D histogram: " << nameDZerovsDStarMass << endl;
+
+ D0mass = new TH1F(nameDZeroMass.Data(), Form("D^{0} invariant mass in bin %d; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
+ DStarMass = new TH1F(nameDStarMass.Data(), Form("Delta invariant mass for candidates in bin %d; M(K#pi#pi)- M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
+ DStarFromSBMass = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invariant mass for sideband in bin %d; M(K#pi#pi)- M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
+ DZerovsDStarMass = new TH2D(nameDZerovsDStarMass.Data(),Form("Delta invariant mass for sideband in bin %d; M(K#pi) GeV/c^{2};M(K#pi#pi)- M(K#pi) GeV/c^{2}",iBin),200,1.75,1.95,200,0.1,0.2);
+
+ if(!fmixing){
+ fDmesonOutput->Add(D0mass);
+ fDmesonOutput->Add(DStarMass);
+ fDmesonOutput->Add(DStarFromSBMass);
+ fDmesonOutput->Add(DZerovsDStarMass);
+ }
+
+ // if using D meson efficiency, define weighted histos
+ if(fUseDmesonEfficiencyCorrection){
+
+ nameDZeroMass = "histDZeroMassWeight_";
+ nameDStarMass = "histDStarMassWeight_";
+ nameDStarFromSBMass = "histDStarFromSBMassWeight_";
+ nameDZerovsDStarMass = "histDZerovsDStarMassWeight_";
+
+ nameDZeroMass+=Form("%d",iBin);
+ nameDStarMass+=Form("%d",iBin);
+ nameDStarFromSBMass+=Form("%d",iBin);
+ nameDZerovsDStarMass+=Form("%d",iBin);
+
+ D0massWeighted = new TH1F(nameDZeroMass.Data(), Form("D^{0} invariant mass in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,1.75,1.95);
+ DStarMassWeighted = new TH1F(nameDStarMass.Data(), Form("Delta invariant mass for candidates in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
+ DStarFromSBMassWeighted = new TH1F(nameDStarFromSBMass.Data(), Form("Delta invariant mass for sideband in bin %d eff weight; M(K#pi) GeV/c^{2};",iBin),200,0.1,0.2);
+ DZerovsDStarMassWeighted = new TH2D(nameDZerovsDStarMass.Data(),Form("Delta invariant mass for sideband in bin %d; M(K#pi) GeV/c^{2};M(K#pi#pi)- M(K#pi) GeV/c^{2}",iBin),200,1.75,1.95,200,0.1,0.2);
+
+ if(!fmixing){
+ fDmesonOutput->Add(D0massWeighted);
+ fDmesonOutput->Add(DStarMassWeighted);
+ fDmesonOutput->Add(DStarFromSBMassWeighted);
+ fDmesonOutput->Add(DZerovsDStarMassWeighted);
+ }
+ }
+ }// end loop on pt bins
+
+
+ // pt distributions
+ TH1F *RecoPtDStar = new TH1F("RecoPtDStar","RECO DStar pt distribution",60,0,60);
+ fDmesonOutput->Add(RecoPtDStar);
+ TH1F *RecoPtBkg = new TH1F("RecoPtBkg","RECO pt distribution side bands",60,0,60);
+ fDmesonOutput->Add(RecoPtBkg);
+
+ TH1D *NumberOfDCandidates = new TH1D("NumberOfDCandidates","Number of D candidates",10,-0.5,9.5);
+ TH1D *NumberOfSBCandidates = new TH1D("NumberOfSBCandidates","Number of SB candidates",10,-0.5,9.5);
+ if(!fmixing) fDmesonOutput->Add(NumberOfDCandidates);
+ if(!fmixing) fDmesonOutput->Add(NumberOfSBCandidates);
+
+ // phi distribution
+ TH2F * PhiInclusiveDStar = new TH2F("PhiInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,50,0,50);
+ TH2F * PhiSidebandDStar = new TH2F("PhiSidebandDStar","Azimuthal distributions of Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,50,0,50);
+
+ // eta distribution
+ TH2F * EtaInclusiveDStar = new TH2F("EtaInclusiveDStar","Azimuthal distributions of Inclusive Dmesons; #eta; pT;Entries",20,-1,1,50,0,50);
+ TH2F * EtaSidebandDStar = new TH2F("EtaSidebandDStar","Azimuthal distributions of Sideband Dmesons; #eta; pT;Entries",20,-1,1,50,0,50);
+
+ if(!fmixing) fDmesonOutput->Add(PhiInclusiveDStar);
+ if(!fmixing) fDmesonOutput->Add(PhiSidebandDStar);
+ if(!fmixing) fDmesonOutput->Add(EtaInclusiveDStar);
+ if(!fmixing) fDmesonOutput->Add(EtaSidebandDStar);
+
+
+ // single track related histograms
+ // phi distribution
+ TH2F * PhiInclusiveTracks = new TH2F("PhiInclusiveTracks","Azimuthal distributions tracks if Inclusive Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
+ TH2F * PhiSidebandTracks = new TH2F("PhiSidebandTracks","Azimuthal distributions tracks if Sideband Dmesons; #phi; pT;Entries",nbinscorr,lowcorrbin,upcorrbin,100,0,10);
+
+ // eta distribution
+ TH2F * EtaInclusiveTracks = new TH2F("EtaInclusiveTracks","Azimuthal distributions of tracks if Inclusive Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
+ TH2F * EtaSidebandTracks = new TH2F("EtaSidebandTracks","Azimuthal distributions of tracks if Sideband Dmesons; #eta; pT;Entries",20,-1,1,100,0,10);
+
+ TH1D * TracksPerDcandidate = new TH1D("TracksPerDcandidate","Distribution of number of tracks per D meson candidate; N tracks; counts",20000,-0.5,19999.5);
+ TH1D * TracksPerSBcandidate = new TH1D("TracksPerSBcandidate","Distribution of number of tracks per sideband candidate; N tracks; counts",20000,-0.5,19999.5);
+ TH1D * TracksPerDMC = new TH1D("TracksPerDMC","Distribution of number of tracks per tagged D ; N tracks; counts",20000,-0.5,19999.5);
+
+ if(!fmixing) fTracksOutput->Add(PhiInclusiveTracks);
+ if(!fmixing) fTracksOutput->Add(PhiSidebandTracks);
+ if(!fmixing) fTracksOutput->Add(EtaInclusiveTracks);
+ if(!fmixing) fTracksOutput->Add(EtaSidebandTracks);
+
+ if(!fmixing) fTracksOutput->Add(TracksPerDcandidate);
+ if(!fmixing) fTracksOutput->Add(TracksPerSBcandidate);
+ if(!fmixing && fmontecarlo) fTracksOutput->Add(TracksPerDMC);
+
+
+ // Montecarlo for D*
+ 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(fmontecarlo) fOutputMC->Add(MCtagPtDStar);
+
+ TH2F * MCTagEtaInclusiveDStar = new TH2F("MCTagEtaInclusiveDStar","MC Tag eta distributions of Inclusive Dmesons; #eta; pT;Entries",20,-1,1,50,0,50);
+ if(fmontecarlo && !fmixing) fOutputMC->Add(MCTagEtaInclusiveDStar);
+
+ TH2F * MCTagPhiInclusiveDStar = new TH2F("MCTagPhiInclusiveDStar","MC Tag Azimuthal distributions of Inclusive Dmesons; #eta; pT;Entries",20,-0.5*TMath::Pi(),1.5*TMath::Pi(),50,0,50);
+ if(fmontecarlo && !fmixing) fOutputMC->Add(MCTagPhiInclusiveDStar);
+
+
+
+
+
+ // event mixing histograms
+ 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) fEMOutput->Add(CheckPoolReadiness);
+
+
+ Int_t NofCentBins = fAssocCuts->GetNCentPoolBins();
+ Int_t NofZVrtxBins = fAssocCuts->GetNZvtxPoolBins();
+ Int_t nPoolBins = NofCentBins*NofZVrtxBins;
+
+ Int_t maxevents = fAssocCuts->GetMaxNEventsInPool();
+
+
+ TH1D * PoolBinDistribution = new TH1D("PoolBinDistribution","Pool Bin Checks; PoolBin; Entry",nPoolBins,-0.5,nPoolBins-0.5);
+ if(fmixing) fEMOutput->Add(PoolBinDistribution);
+
+ TH2D * EventDistributionPerPoolBin = new TH2D("EventDistributionPerPoolBin","Pool Bin Checks; PoolBin; Entry",nPoolBins,-0.5,nPoolBins-0.5,maxevents+2,0,maxevents+2);
+ if(fmixing) fEMOutput->Add(EventDistributionPerPoolBin);
+
+}
+
+
+//__________________________________________________________________________________________________
+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;iv<nvars;iv++){
+ rdcutsvalmine[iv]=new Float_t[fNofPtBins];
+ }
+
+
+ for (Int_t k=0;k<nvars;k++){
+ for (Int_t j=0;j<fNofPtBins;j++){
+
+ // enlarge D0 window
+ if(k==0) {
+ fD0Window[j] =fCuts->GetCutValue(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;iv<nvars;iv++){
+ delete rdcutsvalmine[iv];
+ }
+ delete [] rdcutsvalmine;
+
+}
+
+
+//____________________________ Run checks on event mixing ___________________________________________________
+void AliAnalysisTaskDStarCorrelations::EventMixingChecks(AliAODEvent* AOD){
+ // check this function
+
+ AliCentrality *centralityObj = 0;
+ Int_t multiplicity = -1;
+ Double_t MultipOrCent = -1;
+
+ // get the pool for event mixing
+ if(fSystem != AA){ // pp
+ multiplicity = AOD->GetNumberOfTracks();
+ MultipOrCent = multiplicity; // convert from Int_t to Double_t
+ }
+ if(fSystem == AA){ // PbPb
+
+ centralityObj = ((AliVAODHeader*)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
+}
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff