From 0a65e3255c8a59c9506572b192cb552c3c175d3d Mon Sep 17 00:00:00 2001 From: martinez Date: Thu, 24 Feb 2011 18:51:40 +0000 Subject: [PATCH] Task single muon and jpsi as a fonction of the tracklet multiplicity (Matthieu) --- PWG3/CMakelibPWG3muon.pkg | 2 +- PWG3/PWG3muonLinkDef.h | 1 + PWG3/muon/AddTaskMuonCollisionMultiplicity.C | 67 + ...iAnalysisTaskMuonCollisionMultiplicity.cxx | 671 +++++++++ ...AliAnalysisTaskMuonCollisionMultiplicity.h | 75 + PWG3/muon/AnalysisFunctionOfMultiplicity.C | 1265 +++++++++++++++++ 6 files changed, 2080 insertions(+), 1 deletion(-) create mode 100644 PWG3/muon/AddTaskMuonCollisionMultiplicity.C create mode 100644 PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.cxx create mode 100644 PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.h create mode 100644 PWG3/muon/AnalysisFunctionOfMultiplicity.C diff --git a/PWG3/CMakelibPWG3muon.pkg b/PWG3/CMakelibPWG3muon.pkg index 9301291a465..2aad5889dde 100644 --- a/PWG3/CMakelibPWG3muon.pkg +++ b/PWG3/CMakelibPWG3muon.pkg @@ -25,7 +25,7 @@ # SHLIBS - Shared Libraries and objects for linking (Executables only) # #--------------------------------------------------------------------------------# -set ( SRCS muon/AliHistogramCollection.cxx muon/AliAnalysisTaskESDMuonFilter.cxx muon/AliAnalysisTaskMuonAODfromGeneral.cxx muon/AliAnalysisTaskFromStandardToMuonAOD.cxx muon/AliAnalysisTaskSingleMu.cxx muon/AliAnalysisTaskLUT.cxx muon/AliAnalysisTaskTrigChEff.cxx muon/AliAnalysisTaskLinkToMC.cxx muon/AliAODEventInfo.cxx muon/AliESDMuonTrackCuts.cxx muon/AliAnalysisTaskSingleMuESD.cxx muon/AliCFMuonResTask1.cxx muon/AliCFMuonSingleTask1.cxx muon/AliEventPoolMuon.cxx muon/AliAnalysisTaskCreateMixedDimuons.cxx muon/AliAnalysisTaskMuonAODCreation.cxx muon/AliAnalysisTaskMuonDistributions.cxx muon/AliMuonInfoStoreRD.cxx muon/AliDimuInfoStoreRD.cxx muon/AliMuonInfoStoreMC.cxx muon/AliDimuInfoStoreMC.cxx muon/AliMuonsHFHeader.cxx muon/AliAnalysisTaskSEMuonsHF.cxx muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx muon/AliAnalysisTaskMuonTreeBuilder.cxx muon/AliAnalysisTaskMuonQA.cxx muon/AliAODMuonReplicator.cxx) +set ( SRCS muon/AliHistogramCollection.cxx muon/AliAnalysisTaskESDMuonFilter.cxx muon/AliAnalysisTaskMuonAODfromGeneral.cxx muon/AliAnalysisTaskFromStandardToMuonAOD.cxx muon/AliAnalysisTaskSingleMu.cxx muon/AliAnalysisTaskLUT.cxx muon/AliAnalysisTaskTrigChEff.cxx muon/AliAnalysisTaskLinkToMC.cxx muon/AliAODEventInfo.cxx muon/AliESDMuonTrackCuts.cxx muon/AliAnalysisTaskSingleMuESD.cxx muon/AliCFMuonResTask1.cxx muon/AliCFMuonSingleTask1.cxx muon/AliEventPoolMuon.cxx muon/AliAnalysisTaskCreateMixedDimuons.cxx muon/AliAnalysisTaskMuonAODCreation.cxx muon/AliAnalysisTaskMuonDistributions.cxx muon/AliMuonInfoStoreRD.cxx muon/AliDimuInfoStoreRD.cxx muon/AliMuonInfoStoreMC.cxx muon/AliDimuInfoStoreMC.cxx muon/AliMuonsHFHeader.cxx muon/AliAnalysisTaskSEMuonsHF.cxx muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx muon/AliAnalysisTaskMuonTreeBuilder.cxx muon/AliAnalysisTaskMuonQA.cxx muon/AliAODMuonReplicator.cxx muon/AliAnalysisTaskMuonCollisionMultiplicity.cxx) string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" ) diff --git a/PWG3/PWG3muonLinkDef.h b/PWG3/PWG3muonLinkDef.h index 6cc9000ec0d..36563955237 100644 --- a/PWG3/PWG3muonLinkDef.h +++ b/PWG3/PWG3muonLinkDef.h @@ -34,5 +34,6 @@ #pragma link C++ class AliAnalysisTaskMuonQA+; #pragma link C++ class AliHistogramCollection+; #pragma link C++ class AliHistogramCollectionIterator+; +#pragma link C++ class AliAnalysisTaskMuonCollisionMultiplicity+; #endif diff --git a/PWG3/muon/AddTaskMuonCollisionMultiplicity.C b/PWG3/muon/AddTaskMuonCollisionMultiplicity.C new file mode 100644 index 00000000000..ce69cc48619 --- /dev/null +++ b/PWG3/muon/AddTaskMuonCollisionMultiplicity.C @@ -0,0 +1,67 @@ +AliAnalysisTaskMuonCollisionMultiplicity *AddTaskMuonCollisionMultiplicity() +{ + // + // Task for the determination of the MUON trigger chamber efficiency + // + // lenhardt@subatech.in2p3.fr + // + + // Get the pointer to the existing analysis manager via the static access method. + //============================================================================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTask", "No analysis manager to connect to."); + return NULL; + } + + const Char_t* fileName = "MUON.Multiplicity.THnSparse.root"; + + // Create the task + AliAnalysisTaskMuonCollisionMultiplicity* taskMuonMultiplicity = new AliAnalysisTaskMuonCollisionMultiplicity("MuonMultiplicity"); + + taskMuonMultiplicity->SetZCut(10.0); + taskMuonMultiplicity->SetEtaCut(1.6); + + // Add to the manager + mgr->AddTask(taskMuonMultiplicity); + /* + AliAODInputHandler* aodH = new AliAODInputHandler(); + mgr->SetInputEventHandler(aodH); + */ + + AliESDInputHandler* esdH = new AliESDInputHandler(); + mgr->SetInputEventHandler(esdH); + + /* + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + mcHandler->SetReadTR(kTRUE); + //mcHandler->SetInputPath("alien:///alice/sim/LHC10d3/117222/107/"); + */ + + // + // Create containers for input/output + AliAnalysisDataContainer* cinput0 = mgr->GetCommonInputContainer(); + + AliAnalysisDataContainer *coutput0 = + mgr->CreateContainer("Test", TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + AliAnalysisDataContainer *coutput1 = + mgr->CreateContainer("Trigger", TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + AliAnalysisDataContainer *coutput2 = + mgr->CreateContainer("SingleMuon", TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + AliAnalysisDataContainer *coutput3 = + mgr->CreateContainer("Dimuon", TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + AliAnalysisDataContainer *coutput4 = + mgr->CreateContainer("MonteCarlo", TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + + // Attach input + mgr->ConnectInput (taskMuonMultiplicity, 0, cinput0 ); + // Attach output + mgr->ConnectOutput(taskMuonMultiplicity, 0, coutput0); + mgr->ConnectOutput(taskMuonMultiplicity, 1, coutput1); + mgr->ConnectOutput(taskMuonMultiplicity, 2, coutput2); + mgr->ConnectOutput(taskMuonMultiplicity, 3, coutput3); + mgr->ConnectOutput(taskMuonMultiplicity, 4, coutput4); + + return taskMuonMultiplicity; +} diff --git a/PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.cxx b/PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.cxx new file mode 100644 index 00000000000..61da122d14b --- /dev/null +++ b/PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.cxx @@ -0,0 +1,671 @@ +/************************************************************************** + * Copyright(c) 1998-1999, 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. * + **************************************************************************/ + +/// Compute the number of Muons tracks as a function of the SPD tracklets multiplicity +/// Compare with Monte Carlo tracks +/// Author Matthieu LENHARDT - SUBATECH, Nantes + + +//PWG3/muon includes +#include "AliAnalysisTaskMuonCollisionMultiplicity.h" + +//STEER includes +#include "AliAODEvent.h" +#include "AliAODVertex.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" +#include "AliAODTrack.h" +#include "AliESDMuonTrack.h" +#include "AliAODTracklets.h" +#include "AliAODInputHandler.h" +#include "AliAnalysisManager.h" +#include "AliAODDimuon.h" +#include "AliMCEvent.h" +#include "AliMCParticle.h" +#include "AliMultiplicity.h" + +//ROOT includes +#include +#include +#include +#include +#include +#include +#include +#include + +//___________________________________________________ +AliAnalysisTaskMuonCollisionMultiplicity::AliAnalysisTaskMuonCollisionMultiplicity() + : + AliAnalysisTaskSE(), + fIsInit(kFALSE), + fAOD(0x0), + fESD(0x0), + fZCut(0), + fEtaCut(0), + fTrackletMultiplicity(0), + fTriggerList(0), + fSingleMuonList(0), + fDimuonList(0), + fMonteCarloList(0) +{ + ///Default Constructor +} + + +//___________________________________________________ +AliAnalysisTaskMuonCollisionMultiplicity::AliAnalysisTaskMuonCollisionMultiplicity(const AliAnalysisTaskMuonCollisionMultiplicity& src) + : + AliAnalysisTaskSE(), + fIsInit(kFALSE), + fAOD(0x0), + fESD(0x0), + fZCut(0), + fEtaCut(0), + fTrackletMultiplicity(0), + fTriggerList(0), + fSingleMuonList(0), + fDimuonList(0), + fMonteCarloList(0) +{ + /// copy ctor + src.Copy(*this); + +} + +//___________________________________________________ +AliAnalysisTaskMuonCollisionMultiplicity& AliAnalysisTaskMuonCollisionMultiplicity::operator=(const AliAnalysisTaskMuonCollisionMultiplicity& src) +{ + /// assignement operator + if ( this != &src ) + { + src.Copy(*this); + } + return *this; +} + + + + +//___________________________________________________ +AliAnalysisTaskMuonCollisionMultiplicity::AliAnalysisTaskMuonCollisionMultiplicity(const Char_t *name) + : + AliAnalysisTaskSE(name), + fIsInit(kFALSE), + fAOD(0x0), + fESD(0x0), + fZCut(0), + fEtaCut(0), + fTrackletMultiplicity(0), + fTriggerList(0), + fSingleMuonList(0), + fDimuonList(0), + fMonteCarloList(0) +{ + // Define Inputs and outputs + //DefineInput(0, TChain::Class()); + //DefineInput(1, TChain::Class()); + DefineOutput(0, TList::Class()); + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); + DefineOutput(3, TList::Class()); + DefineOutput(4, TList::Class()); +} + + +//______________________________________________________________________________ +AliAnalysisTaskMuonCollisionMultiplicity::~AliAnalysisTaskMuonCollisionMultiplicity() +{ +// Destructor. + delete fAOD; + delete fESD; + delete fTriggerList; + delete fSingleMuonList; + delete fDimuonList; + delete fMonteCarloList; +} + + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::UserCreateOutputObjects() +{ + if (!fIsInit) + Init(); + OpenFile(0); +} + + + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::UserExec(Option_t */*option*/) +{ + fAOD = 0x0; + fESD = 0x0; + + if (!fIsInit) + Init(); + + fAOD = dynamic_cast (InputEvent()); + if (!fAOD) + fESD = dynamic_cast (InputEvent()); + + + if (fAOD) + CheckEventAOD(); + + if (fESD) + CheckEventESD(); + + PostData(1, fTriggerList); + PostData(2, fSingleMuonList); + PostData(3, fDimuonList); + PostData(4, fMonteCarloList); +} + + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::NotifyRun() +{ + +} + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::FinishTaskOutput() +{ + +} + + +//__________________________________________________________________________ +Bool_t AliAnalysisTaskMuonCollisionMultiplicity::CheckEventAOD() +{ + AliAODVertex *vertex = fAOD->GetPrimaryVertex(); + + if (!vertex) + return kFALSE; + + if (vertex->GetNContributors() < 1) + return kFALSE; + + ComputeMultiplicity(); + if (fTrackletMultiplicity < 1) + return kFALSE; + + // Variables use to determine the type of trigger : + // 0 for minimum bias : CINT1B, CINT1-B or MB1 + // 1 for muon events : CMUS1B, CMUS1-B or MULow + // -1 for everything else + TString trigger = fAOD->GetFiredTriggerClasses(); + Int_t triggerClass = -1; + + if (trigger.Contains("CINT1B") || trigger.Contains("CINT1-B") || trigger.Contains("MB1")) + triggerClass = 0; + + if (trigger.Contains("CMUS1B") || trigger.Contains("CMUS1-B") || trigger.Contains("MULow")) + triggerClass = 1; + + if (triggerClass >= 0) + FillHistosAOD(triggerClass); + + return kTRUE; +} + + + +//__________________________________________________________________________ +Bool_t AliAnalysisTaskMuonCollisionMultiplicity::CheckEventESD() +{ + const AliESDVertex *vertex = fESD->GetPrimaryVertex(); + + if (!vertex) + return kFALSE; + + if (vertex->GetNContributors() < 1) + return kFALSE; + + ComputeMultiplicity(); + if (fTrackletMultiplicity < 1) + return kFALSE; + + // Variables use to determine the type of trigger : + // 0 for minimum bias : CINT1B, CINT1-B or MB1 + // 1 for muon events : CMUS1B, CMUS1-B or MULow + // -1 for everything else + TString trigger = fESD->GetFiredTriggerClasses(); + Int_t triggerClass = -1; + + if (trigger.Contains("CINT1B") || trigger.Contains("CINT1-B") || trigger.Contains("MB1") || trigger.Contains("CMBAC-B") || trigger.Contains("CMBACS2-B")) + triggerClass = 0; + + if (trigger.Contains("CMUS1B") || trigger.Contains("CMUS1-B") || trigger.Contains("MULow")) + triggerClass = 1; + + if (triggerClass >= 0) + FillHistosESD(triggerClass); + + if (fMCEvent) + FillHistosMC(); + + return kTRUE; +} + + + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::FillHistosAOD(Int_t triggerClass) +{ + // Fill histos for AOD events + Int_t nTracks = fAOD->GetNTracks(); + Int_t nDimuons = fAOD->GetNDimuons(); + + // Fill histos + Double_t vertexCut = (TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) < fZCut); + Double_t pileUp = !(fAOD->IsPileupFromSPD()); + + Double_t valuesTrigger[3] = {fTrackletMultiplicity, vertexCut, pileUp}; + ((THnSparseD *)fTriggerList->At(triggerClass))->Fill(valuesTrigger); + + // Loop on the muons tracks + for (Int_t ii = 0; ii < nTracks; ii++) + if (IsUsableMuon(fAOD->GetTrack(ii))) + { + Double_t matchTrigger = fAOD->GetTrack(ii)->GetMatchTrigger(); + if (matchTrigger > 1.0) + matchTrigger = 1.0; // We don't care what type of trigger it is + + Double_t thetaAbs = (180.0 / TMath::Pi()) * TMath::ATan(fAOD->GetTrack(ii)->GetRAtAbsorberEnd()/505.0); + Double_t eta = fAOD->GetTrack(ii)->Eta(); + Double_t dcaX = fAOD->GetTrack(ii)->XAtDCA(); + Double_t dcaY = fAOD->GetTrack(ii)->YAtDCA(); + Double_t p = fAOD->GetTrack(ii)->P(); + Double_t pT = fAOD->GetTrack(ii)->Pt(); + Double_t pDCA = p*TMath::Sqrt(dcaX*dcaX + dcaY*dcaY); + + Double_t valuesMuon[11] = {fTrackletMultiplicity, vertexCut, pileUp, matchTrigger, thetaAbs, eta, p*dcaX, p*dcaY, pDCA, p, pT}; + ((THnSparseD *)fSingleMuonList->At(triggerClass))->Fill(valuesMuon); + } + + // Loop on Dimuons + for (Int_t ii = 0; ii < nDimuons; ii++) + if (fAOD->GetDimuon(ii)->Charge() == 0.0) + if (IsUsableMuon(fAOD->GetDimuon(ii)->GetMu(0))) + if (IsUsableMuon(fAOD->GetDimuon(ii)->GetMu(1))) + { + Double_t matchTrigger1 = fAOD->GetDimuon(ii)->GetMu(0)->GetMatchTrigger(); + if (matchTrigger1 > 0.0) + matchTrigger1 = 1.0; + Double_t matchTrigger2 = fAOD->GetDimuon(ii)->GetMu(1)->GetMatchTrigger(); + if (matchTrigger2 > 0.0) + matchTrigger2 = 1.0; + + Double_t nMatchTrigger = matchTrigger1 + matchTrigger2; + + Double_t thetaAbs1 = (180.0 / TMath::Pi()) * TMath::ATan(fAOD->GetDimuon(ii)->GetMu(0)->GetRAtAbsorberEnd()/505.0); + Double_t thetaAbs2 = (180.0 / TMath::Pi()) * TMath::ATan(fAOD->GetDimuon(ii)->GetMu(1)->GetRAtAbsorberEnd()/505.0); + Double_t eta1 = fAOD->GetDimuon(ii)->GetMu(0)->Eta(); + Double_t eta2 = fAOD->GetDimuon(ii)->GetMu(1)->Eta(); + Double_t dcaX1 = fAOD->GetDimuon(ii)->GetMu(0)->XAtDCA(); + Double_t dcaY1 = fAOD->GetDimuon(ii)->GetMu(0)->YAtDCA(); + Double_t dcaX2 = fAOD->GetDimuon(ii)->GetMu(1)->XAtDCA(); + Double_t dcaY2 = fAOD->GetDimuon(ii)->GetMu(1)->YAtDCA(); + + Double_t p1 = fAOD->GetDimuon(ii)->GetMu(0)->P(); + Double_t p2 = fAOD->GetDimuon(ii)->GetMu(1)->P(); + + Double_t pDCA1 = p1*TMath::Sqrt(dcaX1*dcaX1 + dcaY1*dcaY1); + Double_t pDCA2 = p2 * TMath::Sqrt(dcaX2*dcaX2 + dcaY2*dcaY2); + + Double_t p = fAOD->GetDimuon(ii)->P(); + Double_t pT = fAOD->GetDimuon(ii)->Pt(); + Double_t M = fAOD->GetDimuon(ii)->M(); + + Double_t valuesDimuon[19] = {fTrackletMultiplicity, vertexCut, pileUp, matchTrigger1, matchTrigger2, nMatchTrigger, thetaAbs1, thetaAbs2, + eta1, eta2, p1*dcaX1, p1*dcaY1, p2*dcaX2, p2*dcaY2, p1, p2, p, pT, M}; + ((THnSparseD *)fDimuonList->At(triggerClass))->Fill(valuesDimuon); + } +} + + + +void AliAnalysisTaskMuonCollisionMultiplicity::FillHistosESD(Int_t triggerClass) +{ + // Fill the histos for ESD events + // This is not yet tested (in particular the Dimuon part) + + Int_t nTracks = fESD->GetNumberOfMuonTracks(); + + Double_t vertexCut = (TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) < fZCut); + Double_t pileUp = !(fESD->IsPileupFromSPD()); + + Double_t valuesTrigger[3] = {fTrackletMultiplicity, vertexCut, pileUp}; + ((THnSparseD *)fTriggerList->At(triggerClass))->Fill(valuesTrigger); + + // Loop on the muons tracks + for (Int_t ii = 0; ii < nTracks; ii++) + if (IsUsableMuon(fESD->GetMuonTrack(ii))) + { + Double_t matchTrigger1 = fESD->GetMuonTrack(ii)->GetMatchTrigger(); + if (matchTrigger1 > 1.0) + matchTrigger1 = 1.0; // We don't care what type of trigger it is + + Double_t thetaAbs1 = (180.0 / TMath::Pi()) * TMath::ATan(fESD->GetMuonTrack(ii)->GetRAtAbsorberEnd()/505.0); + Double_t eta1 = fESD->GetMuonTrack(ii)->Eta(); + Double_t dcaX1 = fESD->GetMuonTrack(ii)->GetNonBendingCoorAtDCA(); + Double_t dcaY1 = fESD->GetMuonTrack(ii)->GetBendingCoorAtDCA(); + Double_t p1 = fESD->GetMuonTrack(ii)->P(); + Double_t pUncor1 = fESD->GetMuonTrack(ii)->PUncorrected(); + Double_t pT1 = fESD->GetMuonTrack(ii)->Pt(); + + // The pDCA is computed with the average of p and pUncor + Double_t pDCAX1 = (p1+pUncor1) * dcaX1 / 2.0; + Double_t pDCAY1 = (p1+pUncor1) * dcaY1 / 2.0; + Double_t pDCA1 = (p1+pUncor1) * TMath::Sqrt(dcaX1*dcaX1 + dcaY1*dcaY1) / 2.0; + + Double_t valuesMuon[11] = {fTrackletMultiplicity, vertexCut, pileUp, matchTrigger1, thetaAbs1, eta1, pDCAX1, pDCAY1, pDCA1, p1, pT1}; + ((THnSparseD *)fSingleMuonList->At(triggerClass))->Fill(valuesMuon); + + // Second loop on muons, to fill the dimuons histos + for (Int_t jj = ii+1; jj < nTracks; jj++) + if (IsUsableMuon(fESD->GetMuonTrack(jj))) + if (fESD->GetMuonTrack(ii)->Charge() + fESD->GetMuonTrack(jj)->Charge() == 0.0) + { + Double_t matchTrigger2 = fESD->GetMuonTrack(jj)->GetMatchTrigger(); + if (matchTrigger2 > 0.0) + matchTrigger2 = 1.0; + + Double_t nMatchTrigger = matchTrigger1 + matchTrigger2; + + Double_t thetaAbs2 = (180.0 / TMath::Pi()) * TMath::ATan(fESD->GetMuonTrack(jj)->GetRAtAbsorberEnd()/505.0); + Double_t eta2 = fESD->GetMuonTrack(jj)->Eta(); + Double_t dcaX2 = fESD->GetMuonTrack(jj)->GetNonBendingCoorAtDCA(); + Double_t dcaY2 = fESD->GetMuonTrack(jj)->GetBendingCoorAtDCA(); + Double_t p2 = fESD->GetMuonTrack(jj)->P(); + Double_t pUncor2 = fESD->GetMuonTrack(jj)->PUncorrected(); + + // The pDCA is computed with the average of p and pUncor + Double_t pDCAX2 = (p2+pUncor2) * dcaX2/2.0; + Double_t pDCAY2 = (p2+pUncor2) * dcaY2/2.0; + Double_t pDCA2 = (p2+pUncor2) * TMath::Sqrt(dcaX2*dcaX2 + dcaY2*dcaY2) / 2.0; + + // To compute the p, pT and M of the dimuon, we need a TLorentz vector of the dimuon + Double_t E = fESD->GetMuonTrack(ii)->E() + fESD->GetMuonTrack(jj)->E(); + Double_t pX = fESD->GetMuonTrack(ii)->Px() + fESD->GetMuonTrack(jj)->Px(); + Double_t pY = fESD->GetMuonTrack(ii)->Py() + fESD->GetMuonTrack(jj)->Py(); + Double_t pZ = fESD->GetMuonTrack(ii)->Pz() + fESD->GetMuonTrack(jj)->Pz(); + TLorentzVector *dimuonVector = new TLorentzVector(pX, pY, pZ, E); + dimuonVector->SetPxPyPzE(pX, pY, pZ, E); + + Double_t p = dimuonVector->P(); + Double_t pT = TMath::Sqrt(pX*pX + pY*pY); + Double_t M = dimuonVector->M(); + + Double_t valuesDimuon[19] = {fTrackletMultiplicity, vertexCut, pileUp, matchTrigger1, matchTrigger2, nMatchTrigger, thetaAbs1, thetaAbs2, + eta1, eta2, pDCAX1, pDCAY1, pDCAX2, pDCAY2, p1, p2, p, pT, M}; + ((THnSparseD *)fDimuonList->At(triggerClass))->Fill(valuesDimuon); + delete dimuonVector; + } + } + +} + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::FillHistosMC() +{ + // Fill the histo of the correlation between MC tracks and ESD tracks + + Int_t multiplicityGenerated = 0; + const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks(); + + for (Int_t nn = 0; nn < fMCEvent->GetNumberOfTracks(); nn++) + { + AliMCParticle *particle = (AliMCParticle *) fMCEvent->GetTrack(nn); + Bool_t isGoodMult = kTRUE; + + if (particle->Particle()->GetStatusCode() != 1) + isGoodMult = kFALSE; + + if (particle->Charge() == 0) + isGoodMult = kFALSE; + + if (TMath::Abs(particle->Eta()) > 1.0) + isGoodMult = kFALSE; + + // Check if the particle is a pion, kaon, proton, electron or muon + if (TMath::Abs(particle->PdgCode()) != 211 && TMath::Abs(particle->PdgCode()) != 321 && TMath::Abs(particle->PdgCode()) != 2212 && + TMath::Abs(particle->PdgCode()) != 11 && TMath::Abs(particle->PdgCode()) != 13) + isGoodMult = kFALSE; + + // Check if the distance to vertex is inferior to 1 cm + Double_t distanceToVertex = TMath::Sqrt((particle->Xv() - vertex->GetXv())*(particle->Xv() - vertex->GetXv()) + + (particle->Yv() - vertex->GetYv())*(particle->Yv() - vertex->GetYv()) + + (particle->Zv() - vertex->GetZv())*(particle->Zv() - vertex->GetZv())); + if (distanceToVertex > 1.0) + isGoodMult = kFALSE; + + if (isGoodMult) + multiplicityGenerated += 1; + } + + ((TH2D *)fMonteCarloList->At(0))->Fill(multiplicityGenerated, fTrackletMultiplicity); +} + + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity() +{ + // Compute the collision multiplicity based on AOD or ESD tracklets + + Int_t multiplicity = 0; + + if (fAOD) + { + AliAODTracklets *tracklets = fAOD->GetTracklets(); + Int_t nTracklets = tracklets->GetNumberOfTracklets(); + for (Int_t nn = 0; nn < nTracklets; nn++) + { + Double_t theta = tracklets->GetTheta(nn); + Double_t eta = -TMath::Log(TMath::Tan(theta/2.0)); + + if (TMath::Abs(eta) < fEtaCut) + multiplicity += 1; + } + } + + + else + for (Int_t nn = 0; nn < fESD->GetMultiplicity()->GetNumberOfTracklets(); nn++) + if (TMath::Abs(fESD->GetMultiplicity()->GetEta(nn)) < fEtaCut) + multiplicity += 1; + + + + fTrackletMultiplicity = multiplicity; +} + + + +//________________________________________________________________________ +Bool_t AliAnalysisTaskMuonCollisionMultiplicity::IsUsableMuon(AliAODTrack *track) +{ + // Check if the track is a usable muon track + // Cuts applied : + // - is it a muon track? + // - does it have a pT > 0.0? + + Bool_t isGood = kFALSE; + + if (!track->IsMuonTrack()) + return isGood; + + if (!(track->Pt() > 0.0)) + return isGood; + + isGood = kTRUE; + return isGood; +} + + + +//________________________________________________________________________ +Bool_t AliAnalysisTaskMuonCollisionMultiplicity::IsUsableMuon(AliESDMuonTrack *track) +{ + // Check if the track is a usable muon track + // Cuts applied : + // - is it a muon track? + // - does it have a pT > 0.0? + + Bool_t isGood = kFALSE; + + if (!track->ContainTrackerData()) + return isGood; + + if (!(track->Pt() > 0.0)) + return isGood; + + isGood = kTRUE; + return isGood; +} + + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::Init() +{ + // Initialize the object + + fTriggerList = new TList(); + fSingleMuonList = new TList(); + fDimuonList = new TList(); + fMonteCarloList = new TList(); + + fTriggerList->SetOwner(); + fSingleMuonList->SetOwner(); + fDimuonList->SetOwner(); + fMonteCarloList->SetOwner(); + + + + + // Trigger histos + // dimension 0 : multiplicity of the event + // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? + // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? + Int_t nBinsTrigger[3] = { 150, 2, 2}; + Double_t minRangeTrigger[3] = { 0.0, 0.0, 0.0}; + Double_t maxRangeTrigger[3] = {150.0, 2.0, 2.0}; + THnSparseD *CINT1B = new THnSparseD ("CINT1B", "CINT1B", 3, nBinsTrigger, minRangeTrigger, maxRangeTrigger); + THnSparseD *CMUS1B = new THnSparseD ("CMUS1B", "CMUS1B", 3, nBinsTrigger, minRangeTrigger, maxRangeTrigger); + CINT1B->Sumw2(); + CMUS1B->Sumw2(); + + fTriggerList->AddAt(CINT1B, 0); + fTriggerList->AddAt(CMUS1B, 1); + + + + + + // Muons histos + // dimension 0 : multiplicity of the event + // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? + // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? + // dimension 3 : does the muon match the trigger (0 for no, 1 for yes)? + // dimension 4 : theta_abs of the muon + // dimension 5 : eta of the muon + // dimension 6 : p DCA_x of the muon + // dimension 7 : p DCA_y of the muon + // dimension 8 : p DCA of the muon + // dimension 9 : p of the muon + // dimension 10 : pT of the muon + + Int_t nBinsMuon[11] = { 150, 2, 2, 2, 110, 35, 150, 150, 150, 500, 300}; + Double_t minRangeMuon[11] = { 0.0, 0.0, 0.0, 0.0, 0.0, -5.0, -300.0, -300.0, 0.0, 0.0, 0.0}; + Double_t maxRangeMuon[11] = {150.0, 2.0, 2.0, 2.0, 11.0, -1.5, 300.0, 300.0, 450.0, 100.0, 30.0}; + + THnSparseD *muonCINT1B = new THnSparseD("muonCINT1B", "muonCINT1B", 11, nBinsMuon, minRangeMuon, maxRangeMuon); + THnSparseD *muonCMUS1B = new THnSparseD("muonCMUS1B", "muonCMUS1B", 11, nBinsMuon, minRangeMuon, maxRangeMuon); + muonCINT1B->Sumw2(); + muonCMUS1B->Sumw2(); + + fSingleMuonList->AddAt(muonCINT1B, 0); + fSingleMuonList->AddAt(muonCMUS1B, 1); + + + // Dimuons histos + // dimension 0 : multiplicity of the event + // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? + // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? + // dimension 3 : does the first muon match the trigger (0 for no, 1 for yes)? + // dimension 4 : does the second muon match the trigger (0 for no, 1 for yes)? + // dimension 5 : number of muons matching the trigger in the dimuon + // dimension 6 : theta_abs of the first muon + // dimension 7 : theta_abs of the second muon + // dimension 8 : eta of the first muon + // dimension 9 : eta of the second muon + // dimension 10 : p DCA_x of the first muon + // dimension 11 : p DCA_y of the first muon + // dimension 12 : p DCA_x of the second muon + // dimension 13 : p DCA_y of the second muon + // dimension 14 : p of the first muon + // dimension 15 : p of the second muon + // dimension 16 : p of the dimuon + // dimension 17 : pT of the dimuon + // dimension 18 : invariant mass of the dimuon + + Int_t nBinsDimuon[19] = { 150, 2, 2, 2, 2, 3, 110, 110, 35, 35, 150, 150, 150, 150, 500, 500, 500, 300, 375}; + Double_t minRangeDimuon[19] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -5.0, -5.0, -300.0, -300.0, -300.0, -300.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + Double_t maxRangeDimuon[19] = {150.0, 2.0, 2.0, 2.0, 2.0, 3.0, 11.0, 11.0, -1.5, -1.5, 300.0, 300.0, 300.0, 300.0, 100.0, 100.0, 100.0, 30.0, 15.0}; + + THnSparseD *dimuonCINT1B = new THnSparseD("dimuonCINT1B", "dimuonCINT1B", 19, nBinsDimuon, minRangeDimuon, maxRangeDimuon); + THnSparseD *dimuonCMUS1B = new THnSparseD("dimuonCMUS1B", "dimuonCMUS1B", 19, nBinsDimuon, minRangeDimuon, maxRangeDimuon); + + dimuonCINT1B->Sumw2(); + dimuonCMUS1B->Sumw2(); + + fDimuonList->AddAt(dimuonCINT1B, 0); + fDimuonList->AddAt(dimuonCMUS1B, 1); + + // MonteCarlo Histo + TH2D *correlGenerReco = new TH2D("correlGenerReco", "correlGenerReco", 250, 0.0, 250.0, 250, 0.0, 250.0); + correlGenerReco->GetXaxis()->SetTitle("N ch gener"); + correlGenerReco->GetYaxis()->SetTitle("N reco tracklets"); + + correlGenerReco->Sumw2(); + + fMonteCarloList->AddAt(correlGenerReco, 0); + + fIsInit = kTRUE; +} + + + + + +//________________________________________________________________________ +void AliAnalysisTaskMuonCollisionMultiplicity::Terminate(Option_t */*option*/) +{ +//Terminate analysis + + fTriggerList = (TList *) GetOutputData(0); + fSingleMuonList = (TList *) GetOutputData(1); + fDimuonList = (TList *) GetOutputData(2); +} diff --git a/PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.h b/PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.h new file mode 100644 index 00000000000..68f63e771df --- /dev/null +++ b/PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.h @@ -0,0 +1,75 @@ +#ifndef ALIANALYSISTASKMUONCOLLISIONMULTIPLICITY_H +#define ALIANALYSISTASKMUONCOLLISIONMULTIPLICITY_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/// \ingroup base +/// \class AliAnalysisTaskMultiplicity +/// \compute the number of Muons tracks as a function of the SPD tracklets multiplicity +/// Author Matthieu LENHARDT - SUBATECH, Nantes + +#include "AliAnalysisTaskSE.h" + +class AliAODEvent; +class AliESDEvent; +class AliAODTrack; +class AliESDMuonTrack; +class AliAODDimuon; +class TList; +class TArrayD; + +class AliAnalysisTaskMuonCollisionMultiplicity : public AliAnalysisTaskSE +{ + public: + AliAnalysisTaskMuonCollisionMultiplicity(); + AliAnalysisTaskMuonCollisionMultiplicity(const AliAnalysisTaskMuonCollisionMultiplicity& rhs); + AliAnalysisTaskMuonCollisionMultiplicity& operator=(const AliAnalysisTaskMuonCollisionMultiplicity&); + AliAnalysisTaskMuonCollisionMultiplicity(const Char_t* name); + virtual ~AliAnalysisTaskMuonCollisionMultiplicity(); + + // Implementation of interface methods + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *option); + virtual void NotifyRun(); + virtual void FinishTaskOutput(); + + + Double_t GetZCut() {return fZCut;}; + Double_t GetEtaCut() {return fEtaCut;}; + void SetZCut(Double_t zCut) {fZCut = zCut;}; + void SetEtaCut(Double_t etaCut) {fEtaCut = etaCut;}; + + private: + + Bool_t fIsInit; //< Is the class initialized? + + AliAODEvent *fAOD; //!< AOD Event + AliESDEvent *fESD; //!< ESD Event + + Double_t fZCut; //< Cut on the |z| of the primary vertex + Double_t fEtaCut; //< Cut on the eta cut of the SPD tracklets + + Int_t fTrackletMultiplicity; //< SPD tracklets multiplicity in the current event + + TList *fTriggerList; //< list of all trigger histos + TList *fSingleMuonList; //< List of all single muons histos + TList *fDimuonList; //< List of all dimuons histos + TList *fMonteCarloList; //< List of all histos containing MC info + + + void Init(); + Bool_t CheckEventAOD(); + Bool_t CheckEventESD(); + void ComputeMultiplicity(); + Bool_t IsUsableMuon(AliAODTrack *track); + Bool_t IsUsableMuon(AliESDMuonTrack *track); + void FillHistosAOD(Int_t triggerClass); + void FillHistosESD(Int_t triggerClass); + void FillHistosMC(); + + ClassDef(AliAnalysisTaskMuonCollisionMultiplicity, 1); +}; + +#endif + diff --git a/PWG3/muon/AnalysisFunctionOfMultiplicity.C b/PWG3/muon/AnalysisFunctionOfMultiplicity.C new file mode 100644 index 00000000000..6085e164c6a --- /dev/null +++ b/PWG3/muon/AnalysisFunctionOfMultiplicity.C @@ -0,0 +1,1265 @@ +/************************************************************************************************************* + +This macro analyses the production of SingleMuon and J/Psi as a function of the collision multiplicity in the SPD. +It reads and analyses the output of the Analysis Task PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity. +Refer to the corresponding files to know what the output of this task is. + +Thismacro use a number of input files whose names are hard-coded +These are : +pTRanges.txt : different ranges in pT in which the study is done. +etaRnages.txt : different ranges in eta in which the study is done. + + +*************************************************************************************************************/ + + + + +// ROOT includes +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// RooFit includes +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// std includes +#include + + + + +// These functions are in charge of doing all the analysis +void ProduceTriggerGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut); + +void ProduceSingleMuonRawGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, + Bool_t applyZCut, Bool_t applyPileUpCut, Bool_t applyMatchTrigCut, Bool_t applyThetaAbsCut, Bool_t applyPDCACut); + +void AnalysisSingleMuon(TFile *outputFile); + +void SingleMuonYieldGraph(TFile *outputFile, TGraphErrors *CINT1B); + +void SingleMuonYieldOverMeanMultGraph(TFile *outputFile, TGraphErrors *CINT1B); + +void SingleMuonYieldNormalisedGraph(TFile *outputFile, TGraphErrors *CMUS1B); + +void ProduceFitResults(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, + Bool_t applyZCut, Bool_t applyPileUpCut, Double_t numberMatchTrig, Bool_t applyThetaAbsCut, Bool_t applyPDCACut); + +void FitInvariantMass(TList *fitList, TH1D *invariantMass, Bool_t areParametersFixed, TH1D *referenceParameters); + +void ProduceDimuonGraph(TFile *outputFile, std::vector multiplicityRanges); + + +/////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////// +void AnalysisFunctionOfMultiplicity(Bool_t doSingleMuonAnalysis = kTRUE, Bool_t doJPsiAnalysis = kTRUE, TString inputName = "./Input/MUON.Multiplicity.THnSparse.ESDs.root", TString multiplicityFileName = "multiplicityJPsi.txt", TString outputName = "./Output/MultiplicityResults.LHC10e.JPsi.ESDs.root") +{ + // The main function of the macro + + // Open the input file and create the output file + TFile *inputFile = TFile::Open(inputName, "read"); + TFile *outputFile = TFile::Open(outputName, "recreate"); + + // First of all, we need to figure which conditions we are using + // For the event : the cut on z_vertex and if we have pile up from the SPD + // For the muons : the usual cuts for SMR + Bool_t applyZCut = kTRUE; + Bool_t applyPileUpCut = kTRUE; + Bool_t applyMatchTrigCut = kTRUE; + Bool_t applyThetaAbsCut = kTRUE; + Bool_t applyPDCACut = kTRUE; + Double_t nMatchTrig = 1.0; + + + // Now, we need the multiplicity ranges + // We get it from an input file + ifstream multiplicityFile(multiplicityFileName); + Double_t multiplicityRange = 0.0; + std::vector multiplicityRanges; + while(multiplicityFile >> multiplicityRange) + multiplicityRanges.push_back(multiplicityRange); + + // With that, we can produce the number of CINT1B in each bin + // Specifficaly, we'll save two TGrpahErrors in the output : + // - the graph containing the number of minimum bias event in the bin + // - the graph containing the mean multiplicity in the bin + ProduceTriggerGraph(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut); + + // Now, the muon analysis + if (doSingleMuonAnalysis) + { + // First, we produce the single muon raw graph + // it is the graph giving the number of muons detected in different the multiplicity bins, and in different eta and pT domains + ProduceSingleMuonRawGraph(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut, applyMatchTrigCut, applyThetaAbsCut, applyPDCACut); + + // Raw graph are created, it is time to analyse them + AnalysisSingleMuon(outputFile); + } + + + // Then, the J/Psi analysis + if (doJPsiAnalysis) + { + ProduceFitResults(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut, nMatchTrig, applyThetaAbsCut, applyPDCACut); + ProduceDimuonGraph(outputFile, multiplicityRanges); + } + + inputFile->Close(); + outputFile->Close(); + cout << "End" << endl; +} + + + +/////////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////////// +void ProduceTriggerGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut) +{ + // Compute the number of CINT1B and CMUS1B in each bin + // Along the x axis, the points are placed at the mean multiplicity of the bin, with an error equal to the error on the mean multiplicity + + + // First of all, we need to get the number of CINT1B as a function of the multiplicity in a TH1D + // Retrieving the THnSparse in the input + THnSparseD *inputCINT1B = (THnSparseD *) ((TList *) (inputFile->Get("Trigger;1"))->FindObject("CINT1B") ); + THnSparseD *inputCMUS1B = (THnSparseD *) ((TList *) (inputFile->Get("Trigger;1"))->FindObject("CMUS1B") ); + + + // Reminder of the architecture of this THnSparse + // dimension 0 : multiplicity of the event + // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? + // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? + + // Now, applying the cuts on z vertex and pile-up + if (applyZCut) + { + inputCINT1B->GetAxis(1)->SetRangeUser(1.0, 2.0); + inputCMUS1B->GetAxis(1)->SetRangeUser(1.0, 2.0); + } + if (applyPileUpCut) + { + inputCINT1B->GetAxis(2)->SetRangeUser(1.0, 2.0); + inputCMUS1B->GetAxis(2)->SetRangeUser(1.0, 2.0); + } + + // cuts applied, we can project the THnSparse on TH1D, because it is easier to handle + TH1D *histoCINT1B = inputCINT1B->Projection(0, "E"); + TH1D *histoCMUS1B = inputCMUS1B->Projection(0, "E"); + + + // We can now create and fill the two TGraphErrors + TGraphErrors *graphCINT1B = new TGraphErrors(multiplicityRanges.size() - 1); + graphCINT1B->SetName("graphCINT1B"); + TGraphErrors *graphCMUS1B = new TGraphErrors(multiplicityRanges.size() - 1); + graphCMUS1B->SetName("graphCMUS1B"); + + // Graphs for plotting purposes + TGraphAsymmErrors *plotCINT1B = new TGraphAsymmErrors(multiplicityRanges.size() - 1); + plotCINT1B->SetName("plotCINT1B"); + TGraphAsymmErrors *plotCMUS1B = new TGraphAsymmErrors(multiplicityRanges.size() - 1); + plotCMUS1B->SetName("plotCMUS1B"); + + for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) + { + // CINT1B first + Int_t firstBin = histoCINT1B->GetXaxis()->FindBin(multiplicityRanges[iRange]); + Int_t lastBin = histoCINT1B->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); + histoCINT1B->GetXaxis()->SetRange(firstBin, lastBin); + + Double_t total = 0.0; + Double_t totalErr = 0.0; + total = histoCINT1B->IntegralAndError(firstBin, lastBin, totalErr); + + Double_t mean = histoCINT1B->GetMean(); + Double_t meanErr = histoCINT1B->GetMeanError(); + + graphCINT1B->SetPoint(iRange, mean, total); + graphCINT1B->SetPointError(iRange, meanErr, totalErr); + + plotCINT1B->SetPoint(iRange, mean, total); + plotCINT1B->SetPointError(iRange, + mean - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - mean, + totalErr, totalErr); + + // CMUS1B second + firstBin = histoCMUS1B->GetXaxis()->FindBin(multiplicityRanges[iRange]); + lastBin = histoCMUS1B->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); + histoCMUS1B->GetXaxis()->SetRange(firstBin, lastBin); + + total = histoCMUS1B->IntegralAndError(firstBin, lastBin, totalErr); + + mean = histoCMUS1B->GetMean(); + meanErr = histoCMUS1B->GetMeanError(); + + graphCMUS1B->SetPoint(iRange, mean, total); + graphCMUS1B->SetPointError(iRange, meanErr, totalErr); + + plotCMUS1B->SetPoint(iRange, mean, total); + plotCMUS1B->SetPointError(iRange, + mean - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - mean, + totalErr, totalErr); + } + + // Saving the graphs + outputFile->WriteTObject(plotCINT1B); + outputFile->WriteTObject(plotCMUS1B); + outputFile->WriteTObject(graphCINT1B); + outputFile->WriteTObject(graphCMUS1B); + + delete graphCINT1B; + delete histoCINT1B; + delete inputCINT1B; + delete graphCMUS1B; + delete histoCMUS1B; + delete inputCMUS1B; +} + + +void ProduceSingleMuonRawGraph(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, + Bool_t applyZCut, Bool_t applyPileUpCut, Bool_t applyMatchTrigCut, Bool_t applyThetaAbsCut, Bool_t applyPDCACut) +{ + // Produce the raw graph for Single Muon that will be used all along the analysis + // There are several graph : + // - over all the range in both eta and pT + // - Single Muons Reference : pT > 1 GeV + // - Heavy Flavor Muons : pT > 4 GeV + // - bin by bin in pT and over all the range in eta + // - bin by bin in eta and with single muon reference + // and right now, there is not enough stat to do it bin by bin in both pT and eta + + + // Retrieve the THnSparse and the Minimum bias histo + THnSparseD *inputSingleMuon = (THnSparseD *) ((TList *) (inputFile->Get("SingleMuon;1"))->FindObject("muonCMUS1B") ); + TGraphErrors *CINT1B = (TGraphErrors *) (outputFile->Get("graphCINT1B;1") ); + + // Reminder of the architecture of this THnSparse + // dimension 0 : multiplicity of the event + // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? + // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? + // dimension 3 : does the muon match the trigger (0 for no, 1 for yes)? + // dimension 4 : theta_abs of the muon + // dimension 5 : eta of the muon + // dimension 6 : p DCA_x of the muon + // dimension 7 : p DCA_y of the muon + // dimension 8 : p DCA of the muon + // dimension 9 : p of the muon + // dimension 10 : pT of the muon + + + + // get the pT and eta ranges from files. + // Beware, the names of the files are hard-coded, so make sure to have them in your folder. + // I might change this in the future, but I felt the main function had enough parameters already. + ifstream pTFile ("pTRanges.txt"); + Double_t pTRange = 0.0; + std::vector pTRanges; + while(pTFile >> pTRange) + pTRanges.push_back(pTRange); + + ifstream etaFile ("etaRanges.txt"); + Double_t etaRange = 0.0; + std::vector etaRanges; + while(etaFile >> etaRange) + etaRanges.push_back(etaRange); + + + // Apply all the cuts + // Beware, theta_abs and pDCA cuts are hard-coded + if (applyZCut) + inputSingleMuon->GetAxis(1)->SetRangeUser(1.0, 2.0); + if (applyPileUpCut) + inputSingleMuon->GetAxis(2)->SetRangeUser(1.0, 2.0); + if (applyMatchTrigCut) + inputSingleMuon->GetAxis(3)->SetRangeUser(1.0, 2.0); + if (applyThetaAbsCut) + inputSingleMuon->GetAxis(4)->SetRangeUser(2.0, 9.0); + if (applyPDCACut) + inputSingleMuon->GetAxis(8)->SetRangeUser(0.0, 450.0); // no cut yet, I need to check what are the usual values + + + // First, we get all the raw histos + TList *rawSingleMuon = new TList(); + rawSingleMuon->SetName("rawSingleMuon"); + + // All integrated + // There are some hard cuts + // - pT : 0.5 -> 8 GeV, because we are not confident on what we are seeing outside of these limits + // - eta : -4.0, -> -2.5, acceptance of the spectrometer + + Double_t etaMinAbsolute = -4.0; + Double_t etaMaxAbsolute = -2.5; + inputSingleMuon->GetAxis(5)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute); + + Double_t pTMinAbsolute = 0.5; + Double_t pTMaxAbsolute = 8.0; + inputSingleMuon->GetAxis(10)->SetRangeUser(pTMinAbsolute, pTMaxAbsolute); + + // Putting this in a TH3D, since it is easier to use and allow for computation of errors on the integral + // TH3D is multiplicity, eta and pT + TH3D *histoSingleMuon = inputSingleMuon->Projection(0, 5, 10, "E"); + + // Declare all the histos + TGraphAsymmErrors *allSingleMuon = new TGraphAsymmErrors(multiplicityRanges.size() - 1); + allSingleMuon->SetName("allSingleMuon"); + TGraphAsymmErrors *rawSMR = new TGraphAsymmErrors(multiplicityRanges.size() - 1); // Single Muon Reference (pT > 1.0 GeV) + rawSMR->SetName("rawSMR"); + TGraphAsymmErrors *rawHFM = new TGraphAsymmErrors(multiplicityRanges.size() - 1); // Heavy Flavor Muon (pT > 4.0 GeV) + rawHFM->SetName("rawHFM"); + + + for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) + { + Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]); + Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); + Int_t firstBinY = histoSingleMuon->GetYaxis()->GetFirst(); + Int_t lastBinY = histoSingleMuon->GetYaxis()->GetLast(); + Int_t firstBinZ = histoSingleMuon->GetZaxis()->GetFirst(); + Int_t lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); + + Double_t nSingleMuon = 0.0; + Double_t errSingleMuon = 0.0; + nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); + + allSingleMuon->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); + allSingleMuon->SetPointError(iRange, + CINT1B->GetX()[iRange] - multiplicityRanges[iRange], + multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], + errSingleMuon, + errSingleMuon); + + // SMR + firstBinZ = histoSingleMuon->GetZaxis()->FindBin(1.0); + lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); + nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); + + rawSMR->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); + rawSMR->SetPointError(iRange, + CINT1B->GetX()[iRange] - multiplicityRanges[iRange], + multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], + errSingleMuon, + errSingleMuon); + + // HFM + firstBinZ = histoSingleMuon->GetZaxis()->FindBin(4.0); + lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); + nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); + + rawHFM->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); + rawHFM->SetPointError(iRange, + CINT1B->GetX()[iRange] - multiplicityRanges[iRange], + multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], + errSingleMuon, + errSingleMuon); + } + + TList *integratedSingleMuon = new TList(); + integratedSingleMuon->SetName("rawSingleMuonIntegrated"); + + integratedSingleMuon->AddAt(allSingleMuon, 0); + integratedSingleMuon->AddAt(rawSMR, 1); + integratedSingleMuon->AddAt(rawHFM, 2); + + rawSingleMuon->AddAt(integratedSingleMuon, 0); + //delete allSingleMuon; + //delete rawSMR; + //delete rawHFM; + + // Now, for the study in pT + TList *rawSingleMuonPt = new TList(); + rawSingleMuonPt->SetName("rawSingleMuonPt"); + + for (UInt_t ipT = 0; ipT < pTRanges.size()-1; ipT++) + { + TString nameGraph; + nameGraph.Form("pTRange_%1f-%1f", pTRanges[ipT], pTRanges[ipT+1]); + TGraphAsymmErrors *singleMuonPtRange = new TGraphAsymmErrors(); + singleMuonPtRange->SetName(nameGraph); + + for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) + { + Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]); + Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); + Int_t firstBinY = histoSingleMuon->GetYaxis()->GetFirst(); + Int_t lastBinY = histoSingleMuon->GetYaxis()->GetLast(); + Int_t firstBinZ = histoSingleMuon->GetZaxis()->FindBin(pTRanges[ipT]); + Int_t lastBinZ = histoSingleMuon->GetZaxis()->FindBin(pTRanges[ipT]+1); + + Double_t nSingleMuon = 0.0; + Double_t errSingleMuon = 0.0; + nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); + + singleMuonPtRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); + singleMuonPtRange->SetPointError(iRange, + CINT1B->GetX()[iRange] - multiplicityRanges[iRange], + multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], + errSingleMuon, + errSingleMuon); + + } + + rawSingleMuonPt->AddAt(singleMuonPtRange, ipT); + //delete singleMuonPtRange; + } + + rawSingleMuon->AddAt(rawSingleMuonPt, 1); + + // At last, the study in eta + // Both on SMR and HFM + TList *rawSingleMuonEta = new TList(); + rawSingleMuonEta->SetName("rawSingleMuonEta"); + + for (UInt_t iEta = 0; iEta < etaRanges.size()-1; iEta++) + { + TString nameGraph; + nameGraph.Form("SMRetaRange_%1f-%1f", etaRanges[iEta], etaRanges[iEta+1]); + TGraphAsymmErrors *SMREtaRange = new TGraphAsymmErrors(); + SMREtaRange->SetName(nameGraph); + + nameGraph.Form("HFMetaRange_%1f-%1f", etaRanges[iEta], etaRanges[iEta+1]); + TGraphAsymmErrors *HFMEtaRange = new TGraphAsymmErrors(); + HFMEtaRange->SetName(nameGraph); + + for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++) + { + // First, SMR + Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]); + Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]); + Int_t firstBinY = histoSingleMuon->GetYaxis()->FindBin(etaRanges[iEta]); + Int_t lastBinY = histoSingleMuon->GetYaxis()->FindBin(etaRanges[iEta]+1); + Int_t firstBinZ = histoSingleMuon->GetZaxis()->FindBin(1.0); + Int_t lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); + + Double_t nSingleMuon = 0.0; + Double_t errSingleMuon = 0.0; + nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); + + SMREtaRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); + SMREtaRange->SetPointError(iRange, + CINT1B->GetX()[iRange] - multiplicityRanges[iRange], + multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], + errSingleMuon, + errSingleMuon); + + // Then, HFM + firstBinZ = histoSingleMuon->GetZaxis()->FindBin(4.0); + lastBinZ = histoSingleMuon->GetZaxis()->GetLast(); + + nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon); + + HFMEtaRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon); + HFMEtaRange->SetPointError(iRange, + CINT1B->GetX()[iRange] - multiplicityRanges[iRange], + multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange], + errSingleMuon, + errSingleMuon); + + + } + + rawSingleMuonEta->AddAt(SMREtaRange, iEta); + rawSingleMuonEta->AddAt(HFMEtaRange, iEta + etaRanges.size()-1); + //delete SMREtaRange; + //delete HFMEtaRange; + } + + rawSingleMuon->AddAt(rawSingleMuonEta, 2); + + outputFile->WriteTObject(rawSingleMuon); + delete integratedSingleMuon; + delete rawSingleMuonPt; + delete rawSingleMuonEta; +} + + +void AnalysisSingleMuon(TFile *outputFile) +{ + // Analyse the raw data for single muons + + // First, we retrieve the CINT1B graph, as usual + TGraphErrors *CINT1B = (TGraphErrors *) outputFile->Get("graphCINT1B;1"); + TGraphErrors *CMUS1B = (TGraphErrors *) outputFile->Get("graphCMUS1B;1"); + + // Creation of all the yield graphs + SingleMuonYieldGraph(outputFile, CINT1B); + + // Creation of the yield over mean mult graph + SingleMuonYieldOverMeanMultGraph(outputFile, CINT1B); + + // Creation of the yield with the reference bin normalised to 1, for comparison + SingleMuonYieldNormalisedGraph(outputFile, CMUS1B); +} + + + + +////////////////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////////////////// +void SingleMuonYieldGraph(TFile *outputFile, TGraphErrors *CINT1B) +{ + // This function use the raw graphs of outputfile and the CINT1B graph to create the yield graph + // yield are raw divided by CINT1B + TList *yieldSingleMuon = new TList(); + yieldSingleMuon->SetName("yieldSingleMuon"); + + for (Int_t iList = 0; iList < ((TList *) outputFile->Get("rawSingleMuon;1"))->GetEntries(); iList++) + { + TList *currentList = (TList *) (( TList *) outputFile->Get("rawSingleMuon;1"))->At(iList); + TList *newList = new TList(); + TString listName = currentList->GetName(); + listName.ReplaceAll("raw", 3, "yield", 5); + newList->SetName(listName); + + for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++) + { + TGraphAsymmErrors *rawGraph = (TGraphAsymmErrors *) currentList->At(iGraph); + + TString yieldName = rawGraph->GetName(); + yieldName.ReplaceAll("raw", 3, "yield", 5); + TGraphAsymmErrors *yieldGraph = new TGraphAsymmErrors(rawGraph->GetN()); + yieldGraph->SetName(yieldName); + + // fill the yield graph + for (Int_t iBin = 0; iBin < yieldGraph->GetN(); iBin++ ) + { + if (CINT1B->GetY()[iBin] != 0.0) + { + yieldGraph->SetPoint(iBin, + rawGraph->GetX()[iBin], + rawGraph->GetY()[iBin]/CINT1B->GetY()[iBin]); + if (rawGraph->GetY()[iBin] != 0.0) + { + Double_t error = yieldGraph->GetY()[iBin] * + TMath::Sqrt((rawGraph->GetEYlow()[iBin]/rawGraph->GetY()[iBin])*(rawGraph->GetEYlow()[iBin]/rawGraph->GetY()[iBin]) + + (CINT1B->GetEY()[iBin]/CINT1B->GetY()[iBin])*(CINT1B->GetEY()[iBin]/CINT1B->GetY()[iBin])); + yieldGraph->SetPointError(iBin, + rawGraph->GetEXlow()[iBin], + rawGraph->GetEXhigh()[iBin], + error, error); + } + } + + else + yieldGraph->SetPoint(iBin, 0, 0); + } + + newList->AddAt(yieldGraph, iGraph); + } + + yieldSingleMuon->AddAt(newList, iList); + } + + outputFile->WriteTObject(yieldSingleMuon); +} + + + +////////////////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////////////////// +void SingleMuonYieldOverMeanMultGraph(TFile *outputFile, TGraphErrors *CINT1B) +{ + // This function use the yield graphs of outputfile to create the yield graph divided by the mean multiplicity in each bin + // The CINT1B graph is necessary to get the error on the mean multiplicity + TList *yieldOverMeanMultSingleMuon = new TList(); + yieldOverMeanMultSingleMuon->SetName("yieldOverMeanMultSingleMuon"); + + for (Int_t iList = 0; iList < ((TList *) outputFile->Get("yieldSingleMuon;1"))->GetEntries(); iList++) + { + TList *currentList = (TList *) (( TList *) outputFile->Get("yieldSingleMuon;1"))->At(iList); + TList *newList = new TList(); + TString listName = currentList->GetName(); + listName.ReplaceAll("yield", 5, "yieldOverMeanMult", 17); + newList->SetName(listName); + + for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++) + { + TGraphAsymmErrors *yieldGraph = (TGraphAsymmErrors *) currentList->At(iGraph); + + TString yieldOverMeanMultName = yieldGraph->GetName(); + yieldOverMeanMultName.ReplaceAll("yield", 5, "yieldOverMeanMult", 17); + TGraphAsymmErrors *yieldOverMeanMultGraph = new TGraphAsymmErrors(yieldGraph->GetN()); + yieldOverMeanMultGraph->SetName(yieldOverMeanMultName); + + // fill the yield graph + for (Int_t iBin = 0; iBin < yieldOverMeanMultGraph->GetN(); iBin++ ) + { + if (CINT1B->GetX()[iBin] != 0.0) + { + yieldOverMeanMultGraph->SetPoint(iBin, + yieldGraph->GetX()[iBin], + yieldGraph->GetY()[iBin]/CINT1B->GetX()[iBin]); + if (yieldGraph->GetY()[iBin] != 0.0) + { + Double_t error = yieldOverMeanMultGraph->GetY()[iBin]* + TMath::Sqrt((yieldGraph->GetEYlow()[iBin]/yieldGraph->GetY()[iBin])*(yieldGraph->GetEYlow()[iBin]/yieldGraph->GetY()[iBin]) + + (CINT1B->GetEX()[iBin]/CINT1B->GetX()[iBin])*(CINT1B->GetEX()[iBin]/CINT1B->GetX()[iBin])); + yieldOverMeanMultGraph->SetPointError(iBin, + yieldGraph->GetEXlow()[iBin], + yieldGraph->GetEXhigh()[iBin], + error, error); + } + } + + else + yieldOverMeanMultGraph->SetPoint(iBin, 0, 0); + } + + newList->AddAt(yieldOverMeanMultGraph, iGraph); + } + + yieldOverMeanMultSingleMuon->AddAt(newList, iList); + } + + outputFile->WriteTObject(yieldOverMeanMultSingleMuon); +} + + +////////////////////////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////////////////////// +void SingleMuonYieldNormalisedGraph(TFile *outputFile, TGraphErrors *CMUS1B) +{ + // This function normalise the yield over mean multiplicity so that the reference bin in multiplicity is equal to 1. + // The reference bin is defined it as the bin with the maximum number of CMUS1B + TList *yieldNormalisedSingleMuon = new TList(); + yieldNormalisedSingleMuon->SetName("yieldNormalisedSingleMuon"); + + + // The first step is to find the reference bin in multiplicity + // we define it as the bin with the maximum number of CMUS1B + + Double_t maxCMUS1B = 0.0; + Int_t referenceBin = 0; + for (Int_t iBin = 0; iBin < CMUS1B->GetN(); iBin++) + if (CMUS1B->GetX()[iBin] > 10.0) + if (CMUS1B->GetY()[iBin] > maxCMUS1B) + { + maxCMUS1B = CMUS1B->GetY()[iBin]; + referenceBin = iBin; + } + + // Now the loop, as usual + // We don't propagate the error on the normalisation factor, since this is artificial + for (Int_t iList = 0; iList < ((TList *) outputFile->Get("yieldOverMeanMultSingleMuon;1"))->GetEntries(); iList++) + { + TList *currentList = (TList *) (( TList *) outputFile->Get("yieldOverMeanMultSingleMuon;1"))->At(iList); + TList *newList = new TList(); + TString listName = currentList->GetName(); + listName.ReplaceAll("yieldOverMeanMult", 17, "yieldNormalised", 15); + newList->SetName(listName); + + for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++) + { + TGraphAsymmErrors *yieldOverMeanMultGraph = (TGraphAsymmErrors *) currentList->At(iGraph); + + TString yieldNormalisedName = yieldOverMeanMultGraph->GetName(); + yieldNormalisedName.ReplaceAll("yieldOverMeanMult", 17, "yieldNormalised", 15); + TGraphAsymmErrors *yieldNormalisedGraph = new TGraphAsymmErrors(yieldOverMeanMultGraph->GetN()); + yieldNormalisedGraph->SetName(yieldNormalisedName); + + // fill the yield graph + for (Int_t iBin = 0; iBin < yieldNormalisedGraph->GetN(); iBin++ ) + { + if (yieldOverMeanMultGraph->GetY()[iBin] != 0.0) + { + yieldNormalisedGraph->SetPoint(iBin, + yieldOverMeanMultGraph->GetX()[iBin], + yieldOverMeanMultGraph->GetY()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin]); + if (yieldOverMeanMultGraph->GetY()[iBin] != 0.0) + yieldNormalisedGraph->SetPointError(iBin, + yieldOverMeanMultGraph->GetEXlow()[iBin], + yieldOverMeanMultGraph->GetEXhigh()[iBin], + yieldOverMeanMultGraph->GetEYlow()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin], + yieldOverMeanMultGraph->GetEYhigh()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin]); + } + + else + yieldNormalisedGraph->SetPoint(iBin, 0, 0); + } + + newList->AddAt(yieldNormalisedGraph, iGraph); + } + + yieldNormalisedSingleMuon->AddAt(newList, iList); + } + + outputFile->WriteTObject(yieldNormalisedSingleMuon); +} + + +///////////////////////////////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////////////////////////////// +void ProduceFitResults(TFile *inputFile, TFile *outputFile, std::vector multiplicityRanges, + Bool_t applyZCut, Bool_t applyPileUpCut, Double_t numberMatchTrig, Bool_t applyThetaAbsCut, Bool_t applyPDCACut) +{ + // Produce the raw graph for Dimuons that will be used all along the analysis + // There are two graphs : + // - J/Psi + // - Background in the J/Psi range + + // Retrieve the THnSparse and the Minimum bias histo + THnSparseD *inputDimuon = (THnSparseD *) ((TList *) (inputFile->Get("Dimuon;1"))->FindObject("dimuonCMUS1B") ); + + // Reminder of the structure of the dimuon histo + // dimension 0 : multiplicity of the event + // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)? + // dimension 2 : is it an event without pile up (0 for no, 1 for yes)? + // dimension 3 : does the first muon match the trigger (0 for no, 1 for yes)? + // dimension 4 : does the second muon match the trigger (0 for no, 1 for yes)? + // dimension 5 : number of muons matching the trigger in the dimuon + // dimension 6 : theta_abs of the first muon + // dimension 7 : theta_abs of the second muon + // dimension 8 : eta of the first muon + // dimension 9 : eta of the second muon + // dimension 10 : p DCA_x of the first muon + // dimension 11 : p DCA_y of the first muon + // dimension 12 : p DCA_x of the second muon + // dimension 13 : p DCA_y of the second muon + // dimension 14 : p of the first muon + // dimension 15 : p of the second muon + // dimension 16 : p of the dimuon + // dimension 17 : pT of the dimuon + // dimension 18 : invariant mass of the dimuon + + + + // get the pT and eta ranges from files. + // Beware, the names of the files are hard-coded, so make sure to have them in your folder. + // I might change this in the future, but I felt the main function had enough parameters already. + ifstream pTFile ("pTRanges.txt"); + Double_t pTRange = 0.0; + std::vector pTRanges; + while(pTFile >> pTRange) + pTRanges.push_back(pTRange); + + ifstream etaFile ("etaRanges.txt"); + Double_t etaRange = 0.0; + std::vector etaRanges; + while(etaFile >> etaRange) + etaRanges.push_back(etaRange); + + + // Apply all the cuts + // Beware, theta_abs and pDCA cuts are hard-coded + if (applyZCut) + inputDimuon->GetAxis(1)->SetRangeUser(1.0, 2.0); + if (applyPileUpCut) + inputDimuon->GetAxis(2)->SetRangeUser(1.0, 2.0); + inputDimuon->GetAxis(5)->SetRangeUser(numberMatchTrig, 3.0); + if (applyThetaAbsCut) + { + inputDimuon->GetAxis(6)->SetRangeUser(2.0, 9.0); + inputDimuon->GetAxis(7)->SetRangeUser(2.0, 9.0); + } + if (applyPDCACut) + {// no cut yet, I need to check what are the usual values + inputDimuon->GetAxis(12)->SetRangeUser(0.0, 450.0); + inputDimuon->GetAxis(15)->SetRangeUser(0.0, 450.0); + } + + + // First, we get the invariant mass histos + TList *dimuonFitResults = new TList(); + dimuonFitResults->SetName("dimuonFitResults"); + + // There are some hard cuts + // - eta : -4.0, -> -2.5, acceptance of the spectrometer + + Double_t etaMinAbsolute = -4.0; + Double_t etaMaxAbsolute = -2.5; + inputDimuon->GetAxis(8)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute); + inputDimuon->GetAxis(9)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute); + + // First, we get the invariant mass histo for the whole range in multiplicity + inputDimuon->GetAxis(0)->SetRangeUser(multiplicityRanges[0], multiplicityRanges[multiplicityRanges.size()-1]); + TH1D *invariantMassIntegrated = inputDimuon->Projection(18, "E"); + invariantMassIntegrated->SetName("invariantMassIntegrated"); + + + // Now for the fit on all the multiplicity + // This is used to get a value for some parameters + TList *fitAll = new TList(); + fitAll->SetName("AllMultiplicity"); + + + // Create the container for the reference parameters + TH1D *referenceParameters = new TH1D("referenceParameters", "parameters", 8, 1.0, 9.0); + referenceParameters->Sumw2(); + referenceParameters->GetXaxis()->SetBinLabel(1, "meanJPsi"); + referenceParameters->GetXaxis()->SetBinLabel(2, "sigmaJPsi"); + referenceParameters->GetXaxis()->SetBinLabel(3, "alphaJPsi"); + referenceParameters->GetXaxis()->SetBinLabel(4, "nJPsi"); + referenceParameters->GetXaxis()->SetBinLabel(5, "meanPsiPrime"); + referenceParameters->GetXaxis()->SetBinLabel(6, "sigmaPsiPrime"); + referenceParameters->GetXaxis()->SetBinLabel(7, "expBkg1"); + referenceParameters->GetXaxis()->SetBinLabel(8, "expBkg2"); + + + FitInvariantMass(fitAll, invariantMassIntegrated, kFALSE, referenceParameters); + dimuonFitResults->AddAt(fitAll, 0); + + // Now, we make a fit for each range in multiplicity + for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) + { + TString name; + name.Form("Multiplicity_%d-%d", static_cast(multiplicityRanges[iMult]), static_cast(multiplicityRanges[iMult+1])); + TList *fitRange = new TList(); + fitRange->SetName(name); + + inputDimuon->GetAxis(0)->SetRangeUser(multiplicityRanges[iMult], multiplicityRanges[iMult+1]); + TH1D *invariantMassRange = inputDimuon->Projection(18, "E"); + invariantMassRange->SetName("invariantMass"); + + FitInvariantMass(fitRange, invariantMassRange, kTRUE, referenceParameters); + dimuonFitResults->AddAt(fitRange, iMult + 1); + } + + outputFile->WriteTObject(dimuonFitResults); +} + + + + +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// +void FitInvariantMass(TList *fitList, TH1D *invariantMass, Bool_t areParametersFixed, TH1D *referenceParameters) +{ + // Fit the given histo and put the results in outputFile + // The fit is : Crystal Ball (J/Psi) + gaussian (Psi') + two exponentials (background) + + // The parameters histo architecture is the following : + // bin 1 : mean J/Psi + // bin 2 : sigma J/Psi + // bin 3 : alpha J/Psi + // bin 4 : n J/Psi + // bin 5 : mean Psi' + // bin 6 : sigma Psi' + // bin 7 : exponential factor background 1 + // bin 8 : exponential factor background 2 + + // The results histo architecture is the following : + // bin 1 : J/Psi + // bin 2 : Background J/Psi at 2 sigma + // bin 3 : signal over background + // bin 4 : chi2 + + // There are some parameters that are not always free : + // mean of J/Psi peak + // sigma of J/Psi peak + // alpha of J/Psi function (Crystall Ball function) + // n of J/Psi function + // mean of Psi' peak + // sigma of Psi' peak + // If there are fixed, their values are taken from the parameters histo + + + // First, we have to create two histos, one that will contain the parameters of the fit, and the other the results + TH1D *parameters = new TH1D("parameters", "parameters", 8, 1.0, 9.0); + parameters->Sumw2(); + parameters->GetXaxis()->SetBinLabel(1, "meanJPsi"); + parameters->GetXaxis()->SetBinLabel(2, "sigmaJPsi"); + parameters->GetXaxis()->SetBinLabel(3, "alphaJPsi"); + parameters->GetXaxis()->SetBinLabel(4, "nJPsi"); + parameters->GetXaxis()->SetBinLabel(5, "meanPsiPrime"); + parameters->GetXaxis()->SetBinLabel(6, "sigmaPsiPrime"); + parameters->GetXaxis()->SetBinLabel(7, "expBkg1"); + parameters->GetXaxis()->SetBinLabel(8, "expBkg2"); + TH1D *results = new TH1D("results", "results", 4, 1.0, 5.0); + results->Sumw2(); + results->GetXaxis()->SetBinLabel(1, "JPsi"); + results->GetXaxis()->SetBinLabel(2, "bkg"); + results->GetXaxis()->SetBinLabel(3, "SB"); + results->GetXaxis()->SetBinLabel(4, "chi2"); + + + // Declare observable M + RooRealVar M("M","Dimuon invariant mass [GeV/c^2]", 2.0, 5.0); + + + // Declare Cristal Ball for J/Psi + RooRealVar mean_jpsi("mean_jpsi","mean_jpsi", 3.096, 2.8, 3.2); + RooRealVar sigma_jpsi("sigma_jpsi","sigma_jpsi", 0.070, 0, 0.2); + RooRealVar alpha_jpsi("alpha","alpha", 5, 0, 10); + RooRealVar n_jpsi("n","n",5,0,10); + RooCBShape jPsi("jpsi", "crystal ball PDF", M, mean_jpsi, sigma_jpsi, alpha_jpsi, n_jpsi); + + + //Declare gaussian for Psi prime + RooRealVar mean_psip("mean_psip", "mean_psip", 3.686, 3.6, 3.75); + RooRealVar sigma_psip("sigma_psip","sigma_psip", 0.086); + RooGaussian psiPrime("psip","Gaussian", M, mean_psip, sigma_psip); + + //Declare exponentials for background + RooRealVar c1("c1", "c1", -10, -0.1); + RooExponential bkg_exp1("bkg_exp1", "background", M, c1); + + RooRealVar c2("c2", "c2", -5.0, -1); + RooExponential bkg_exp2("bkg_exp2", "background", M, c2); + + // Give a value and fix sigma, alpha and n in the Cristal Ball + // The values correspond to a previous fit (ideally, with all the statistic) + // I should automatize this one of these day + if (areParametersFixed) + { + Double_t fixMeanJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("meanJPsi")); + mean_jpsi.setVal(fixMeanJPsi); + mean_jpsi.setError(0.0); + mean_jpsi.setRange(fixMeanJPsi, fixMeanJPsi); + + Double_t fixSigmaJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi")); + sigma_jpsi.setVal(fixSigmaJPsi); + sigma_jpsi.setError(0.0); + sigma_jpsi.setRange(fixSigmaJPsi, fixSigmaJPsi); + + Double_t fixAlphaJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi")); + alpha_jpsi.setVal(fixAlphaJPsi); + alpha_jpsi.setError(0.0); + alpha_jpsi.setRange(fixAlphaJPsi, fixAlphaJPsi); + + Double_t fixNJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("nJPsi")); + n_jpsi.setVal(fixNJPsi); + n_jpsi.setError(0.0); + n_jpsi.setRange(fixNJPsi, fixNJPsi); + + Double_t fixMeanPsiPrime = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime")); + mean_psip.setVal(fixMeanPsiPrime); + mean_psip.setError(0.0); + mean_psip.setRange(fixMeanPsiPrime, fixMeanPsiPrime); + + Double_t fixSigmaPsiPrime = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime")); + sigma_psip.setVal(fixSigmaPsiPrime); + sigma_psip.setError(0.0); + sigma_psip.setRange(fixSigmaPsiPrime, fixSigmaPsiPrime); + } + + + // Sum the composite signal and background into an extended pdf nsig*sig+nbkg*bkg + RooRealVar fitJPsi("fitJPsi", "number of J/Psi events", 1600 ,0.0, 15000); + RooRealVar fitPsiPrime("fitPsiPrime", "number of Psi Prime events", 50, 0.0, 200); + RooRealVar fitBkg1("fitBkg1", "number of background events", 5000, 0.0, 15000); + RooRealVar fitBkg2("fitBkg2", "number of background events", 5000, 0.0, 15000); + + + RooAddPdf *fitFunction = new RooAddPdf("model", "(CB+Gauss+2exp)", RooArgList(bkg_exp1, bkg_exp2, jPsi, psiPrime), RooArgList(fitBkg1, fitBkg2, fitJPsi, fitPsiPrime)); + + // Define the histo to be fitted + RooDataHist *invariantMassFit = new RooDataHist("invariantMassFit", "invariantMassFit", M, invariantMass); + + // Fit the invariant mass + // This is the important part + RooFitResult *fitResult = fitFunction->fitTo(*invariantMassFit, RooFit::Save()); + //RooFitResult *fitResult = fitFunction->chi2FitTo(*invariantMassFit, RooFit::Save()); + + RooPlot *plot = M.frame(); + TString title = "fittedInvariantMass"; + plot->SetTitle(title); + plot->SetName(title); + + invariantMassFit->plotOn(plot, RooFit::Name("invariantMassFit")); + fitFunction->plotOn(plot, RooFit::Name("fitFunction"), RooFit::Range(2.0, 5.0)); + + + // Fill the histo with the final values of the parameters + parameters->SetBinContent(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getError()); + + parameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getError()); + + parameters->SetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getError()); + + parameters->SetBinContent(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getError()); + + parameters->SetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getError()); + + parameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getError()); + + // Fill the histo with the reference parameters + if (!areParametersFixed) + { + referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getVal()); + referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getError()); + + referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getVal()); + referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getError()); + + referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getVal()); + referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getError()); + + referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getVal()); + referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getError()); + + referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getVal()); + referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getError()); + + referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getVal()); + referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getError()); + } + + parameters->SetBinContent(parameters->GetXaxis()->FindBin("expBkg1"), c1.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("expBkg1"), c1.getError()); + + parameters->SetBinContent(parameters->GetXaxis()->FindBin("expBkg2"), c2.getVal()); + parameters->SetBinError(parameters->GetXaxis()->FindBin("expBkg2"), c2.getError()); + + + // Define the range to compute the results + // We choose two sigma + // We also need two new variables to compute the background + Double_t lowerBoundJPsi = mean_jpsi.getVal() - 2.0*sigma_jpsi.getVal(); + Double_t upperBoundJPsi = mean_jpsi.getVal() + 2.0*sigma_jpsi.getVal(); + M.setRange("JPsiRange" ,lowerBoundJPsi, upperBoundJPsi); + + RooAbsReal *fracBkgRange1 = bkg_exp1.createIntegral(M, M, "JPsiRange"); + RooAbsReal *fracBkgRange2 = bkg_exp2.createIntegral(M, M, "JPsiRange"); + + // Get the results, and fill the results histo + Double_t nJPsi = fitJPsi.getVal(); + Double_t errJPsi = fitJPsi.getError(); + results->SetBinContent(results->GetXaxis()->FindBin("JPsi"), nJPsi); + results->SetBinError(results->GetXaxis()->FindBin("JPsi"), errJPsi); + + Double_t nBkg = (fitBkg1.getVal() * fracBkgRange1->getVal() + fitBkg2.getVal() * fracBkgRange2->getVal()); + Double_t errBkg = (fitBkg1.getError() * fracBkgRange1->getVal() + fitBkg2.getError() * fracBkgRange2->getVal()); + results->SetBinContent(results->GetXaxis()->FindBin("bkg"), nBkg); + results->SetBinError(results->GetXaxis()->FindBin("bkg"), errBkg); + + Double_t SB = 0.0; + Double_t errSB = 0.0; + if (nJPsi != 0.0 && nBkg != 0.0) + { + SB = nJPsi/nBkg; + errSB = SB * (errJPsi/nJPsi + errBkg/nBkg); + } + results->SetBinContent(results->GetXaxis()->FindBin("SB"), SB); + results->SetBinError(results->GetXaxis()->FindBin("SB"), errSB); + + Int_t nDF = fitResult->floatParsFinal().getSize(); + Double_t chi2 = plot->chiSquare("fitFunction", "invariantMassFit", nDF); + results->SetBinContent(results->GetXaxis()->FindBin("chi2"), chi2); + results->SetBinError(results->GetXaxis()->FindBin("chi2"), 0.0); + + fitList->AddAt(invariantMass, 0); + fitList->AddAt(plot, 1); + fitList->AddAt(parameters, 2); + fitList->AddAt(results, 3); +} + + +void ProduceDimuonGraph(TFile *outputFile, std::vector multiplicityRanges) +{ + // This function will create all the graphs for dimuons (J/Psi and Background) + // - raw + // - yield + // - yield over mean mult + // - yield over mean mult normalised + + // First, get the CINT1B and CMUS1B graphs + TGraphErrors *CINT1B = (TGraphErrors *) outputFile->Get("graphCINT1B"); + TGraphErrors *CMUS1B = (TGraphErrors *) outputFile->Get("graphCMUS1B"); + + TList *JPsiGraph = new TList(); + JPsiGraph->SetName("JPsiGraph"); + TList *bkgGraph = new TList(); + bkgGraph->SetName("bkgGraph"); + + // Raw Graphs + TGraphAsymmErrors *rawJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + rawJPsiGraph->SetName("rawJPsiGraph"); + TGraphAsymmErrors *rawBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + rawBkgGraph->SetName("rawBkgGraph"); + + for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) + { + TString name; + name.Form("Multiplicity_%d-%d", static_cast(multiplicityRanges[iMult]), static_cast(multiplicityRanges[iMult+1])); + TH1D *results = (TH1D *) ((TList *) ((TList *) outputFile->Get("dimuonFitResults;1"))->FindObject(name))->FindObject("results"); + + // J/Psi + Double_t nJPsi = results->GetBinContent(results->GetXaxis()->FindBin("JPsi")); + Double_t errJPsi = results->GetBinError(results->GetXaxis()->FindBin("JPsi")); + + rawJPsiGraph->SetPoint(iMult, CINT1B->GetX()[iMult], nJPsi); + rawJPsiGraph->SetPointError(iMult, + CINT1B->GetX()[iMult] - multiplicityRanges[iMult], + multiplicityRanges[iMult+1] - CINT1B->GetX()[iMult], + errJPsi, errJPsi); + + // Background + Double_t nBkg = results->GetBinContent(results->GetXaxis()->FindBin("bkg")); + Double_t errBkg = results->GetBinError(results->GetXaxis()->FindBin("bkg")); + + rawJPsiGraph->SetPoint(iMult, CINT1B->GetX()[iMult], nBkg); + rawJPsiGraph->SetPointError(iMult, + CINT1B->GetX()[iMult] - multiplicityRanges[iMult], + multiplicityRanges[iMult+1] - CINT1B->GetX()[iMult], + errBkg, errBkg); + } + + JPsiGraph->AddAt(rawJPsiGraph, 0); + bkgGraph->AddAt(rawBkgGraph, 0); + + // Yield Graphs + TGraphAsymmErrors *yieldJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + yieldJPsiGraph->SetName("yieldJPsiGraph"); + TGraphAsymmErrors *yieldBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + yieldBkgGraph->SetName("yieldBkgGraph"); + + for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) + { + if (CINT1B->GetY()[iMult] != 0.0) + { + // J/Psi + yieldJPsiGraph->SetPoint(iMult, rawJPsiGraph->GetX()[iMult], rawJPsiGraph->GetY()[iMult]/CINT1B->GetY()[iMult]); + if (rawJPsiGraph->GetY()[iMult] != 0.0) + { + Double_t error = yieldJPsiGraph->GetY()[iMult]* + TMath::Sqrt((rawJPsiGraph->GetEYlow()[iMult]/rawJPsiGraph->GetY()[iMult])*(rawJPsiGraph->GetEYlow()[iMult]/rawJPsiGraph->GetY()[iMult]) + + (CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])*(CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])); + yieldJPsiGraph->SetPointError(iMult, + rawJPsiGraph->GetEXlow()[iMult], + rawJPsiGraph->GetEXhigh()[iMult], + error, error); + } + // Background + yieldBkgGraph->SetPoint(iMult, rawBkgGraph->GetX()[iMult], rawBkgGraph->GetY()[iMult]/CINT1B->GetY()[iMult]); + if (rawBkgGraph->GetY()[iMult] != 0.0) + { + Double_t error = yieldBkgGraph->GetY()[iMult]* + TMath::Sqrt((rawBkgGraph->GetEYlow()[iMult]/rawBkgGraph->GetY()[iMult])*(rawBkgGraph->GetEYlow()[iMult]/rawBkgGraph->GetY()[iMult]) + + (CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])*(CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])); + + yieldBkgGraph->SetPointError(iMult, + rawBkgGraph->GetEXlow()[iMult], + rawBkgGraph->GetEXhigh()[iMult], + error, error); + } + } + + else + { + yieldJPsiGraph->SetPoint(iMult, rawJPsiGraph->GetX()[iMult], 0); + yieldBkgGraph->SetPoint(iMult, rawBkgGraph->GetX()[iMult], 0); + } + } + + JPsiGraph->AddAt(yieldJPsiGraph, 1); + bkgGraph->AddAt(yieldBkgGraph, 1); + + // Yield over Mean Mult graph + TGraphAsymmErrors *yieldOverMeanMultJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + yieldOverMeanMultJPsiGraph->SetName("yieldOverMeanMultJPsiGraph"); + TGraphAsymmErrors *yieldOverMeanMultBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + yieldOverMeanMultBkgGraph->SetName("yieldOverMeanMultBkgGraph"); + + for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) + { + if (CINT1B->GetX()[iMult] != 0.0) + { + // J/Psi + yieldOverMeanMultJPsiGraph->SetPoint(iMult, yieldJPsiGraph->GetX()[iMult], yieldJPsiGraph->GetY()[iMult]/CINT1B->GetX()[iMult]); + if (yieldJPsiGraph->GetY()[iMult] != 0.0) + { + Double_t error = yieldOverMeanMultJPsiGraph->GetY()[iMult]* + TMath::Sqrt((yieldJPsiGraph->GetEYlow()[iMult]/yieldJPsiGraph->GetY()[iMult])*(yieldJPsiGraph->GetEYlow()[iMult]/yieldJPsiGraph->GetY()[iMult]) + + (CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])*(CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])); + yieldOverMeanMultJPsiGraph->SetPointError(iMult, + yieldJPsiGraph->GetEXlow()[iMult], + yieldJPsiGraph->GetEXhigh()[iMult], + error, error); + } + // Background + yieldOverMeanMultBkgGraph->SetPoint(iMult, yieldBkgGraph->GetX()[iMult], yieldBkgGraph->GetY()[iMult]/CINT1B->GetX()[iMult]); + if (yieldBkgGraph->GetY()[iMult] != 0.0) + { + Double_t error = yieldOverMeanMultBkgGraph->GetY()[iMult]* + TMath::Sqrt((yieldBkgGraph->GetEYlow()[iMult]/yieldBkgGraph->GetY()[iMult])*(yieldBkgGraph->GetEYlow()[iMult]/yieldBkgGraph->GetY()[iMult]) + + (CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])*(CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])); + yieldOverMeanMultBkgGraph->SetPointError(iMult, + yieldBkgGraph->GetEXlow()[iMult], + yieldBkgGraph->GetEXhigh()[iMult], + error, error); + } + } + + else + { + yieldOverMeanMultJPsiGraph->SetPoint(iMult, yieldJPsiGraph->GetX()[iMult], 0); + yieldOverMeanMultBkgGraph->SetPoint(iMult, yieldBkgGraph->GetX()[iMult], 0); + } + } + + JPsiGraph->AddAt(yieldOverMeanMultJPsiGraph, 2); + bkgGraph->AddAt(yieldOverMeanMultBkgGraph, 2); + + // Yield over Mean Mult normalised to get the bin with the highest number of CMUS1B equal to 1 + Double_t maxCMUS1B = 0.0; + Int_t referenceBin = 0; + for (Int_t iBin = 0; iBin < CMUS1B->GetN(); iBin++) + if (CMUS1B->GetX()[iBin] > 10.0) + if (CMUS1B->GetY()[iBin] > maxCMUS1B) + { + maxCMUS1B = CMUS1B->GetY()[iBin]; + referenceBin = iBin; + } + + TGraphAsymmErrors *yieldNormalisedJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + yieldNormalisedJPsiGraph->SetName("yieldNormalisedJPsiGraph"); + TGraphAsymmErrors *yieldNormalisedBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1); + yieldNormalisedBkgGraph->SetName("yieldNormalisedBkgGraph"); + + for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++) + { + // JPsi + if (yieldOverMeanMultJPsiGraph->GetY()[referenceBin] != 0.0) + { + yieldNormalisedJPsiGraph->SetPoint(iMult, + yieldOverMeanMultJPsiGraph->GetX()[iMult], + yieldOverMeanMultJPsiGraph->GetY()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin]); + yieldNormalisedJPsiGraph->SetPointError(iMult, + yieldOverMeanMultJPsiGraph->GetEXlow()[iMult], + yieldOverMeanMultJPsiGraph->GetEXhigh()[iMult], + yieldOverMeanMultJPsiGraph->GetEYlow()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin], + yieldOverMeanMultJPsiGraph->GetEYhigh()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin]); + } + + // Bkg + if (yieldOverMeanMultBkgGraph->GetY()[referenceBin] != 0.0) + { + yieldNormalisedBkgGraph->SetPoint(iMult, + yieldOverMeanMultBkgGraph->GetX()[iMult], + yieldOverMeanMultBkgGraph->GetY()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin]); + yieldNormalisedBkgGraph->SetPointError(iMult, + yieldOverMeanMultBkgGraph->GetEXlow()[iMult], + yieldOverMeanMultBkgGraph->GetEXhigh()[iMult], + yieldOverMeanMultBkgGraph->GetEYlow()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin], + yieldOverMeanMultBkgGraph->GetEYhigh()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin]); + } + + } + + JPsiGraph->AddAt(yieldNormalisedJPsiGraph, 3); + bkgGraph->AddAt(yieldNormalisedBkgGraph, 3); + + outputFile->WriteTObject(JPsiGraph); + outputFile->WriteTObject(bkgGraph); +} -- 2.43.0