Task single muon and jpsi as a fonction of the tracklet multiplicity (Matthieu)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 24 Feb 2011 18:51:40 +0000 (18:51 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 24 Feb 2011 18:51:40 +0000 (18:51 +0000)
PWG3/CMakelibPWG3muon.pkg
PWG3/PWG3muonLinkDef.h
PWG3/muon/AddTaskMuonCollisionMultiplicity.C [new file with mode: 0644]
PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.cxx [new file with mode: 0644]
PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.h [new file with mode: 0644]
PWG3/muon/AnalysisFunctionOfMultiplicity.C [new file with mode: 0644]

index 9301291..2aad588 100644 (file)
@@ -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}" )
 
index 6cc9000..3656395 100644 (file)
@@ -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 (file)
index 0000000..ce69cc4
--- /dev/null
@@ -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 (file)
index 0000000..61da122
--- /dev/null
@@ -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 <TH2D.h>
+#include <THnSparse.h>
+#include <TChain.h>
+#include <TList.h>
+#include <TArrayD.h>
+#include <Riostream.h>
+#include <TParticle.h>
+#include <TLorentzVector.h>
+
+//___________________________________________________
+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<AliAODEvent *> (InputEvent());
+  if (!fAOD)
+    fESD = dynamic_cast<AliESDEvent *> (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 (file)
index 0000000..68f63e7
--- /dev/null
@@ -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 (file)
index 0000000..6085e16
--- /dev/null
@@ -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 <TH1D.h>
+#include <TH2D.h>
+#include <TH3D.h>
+#include <THnSparse.h>
+#include <TGraphErrors.h>
+#include <TGraphAsymmErrors.h>
+#include <TFile.h>
+#include <TList.h>
+#include <Riostream.h>
+#include <TMath.h>
+
+// RooFit includes
+#include <RooRealVar.h>
+#include <RooAbsReal.h>
+#include <RooArgSet.h>
+#include <RooCBShape.h>
+#include <RooGaussian.h>
+#include <RooExponential.h>
+#include <RooAddPdf.h>
+#include <RooDataHist.h>
+#include <RooFitResult.h>
+#include <RooPlot.h>
+
+// std includes
+#include <vector>
+
+
+
+
+// These functions are in charge of doing all the analysis
+void ProduceTriggerGraph(TFile *inputFile, TFile *outputFile, std::vector<Double_t> multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut);
+
+void ProduceSingleMuonRawGraph(TFile *inputFile, TFile *outputFile, std::vector<Double_t> 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<Double_t> 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<Double_t> 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 <Double_t> 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<Double_t> 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<Double_t> 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 <Double_t> pTRanges;
+  while(pTFile >> pTRange)
+    pTRanges.push_back(pTRange);
+
+  ifstream etaFile ("etaRanges.txt");
+  Double_t etaRange = 0.0;
+  std::vector <Double_t> 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<Double_t> 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 <Double_t> pTRanges;
+  while(pTFile >> pTRange)
+    pTRanges.push_back(pTRange);
+
+  ifstream etaFile ("etaRanges.txt");
+  Double_t etaRange = 0.0;
+  std::vector <Double_t> 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<Int_t>(multiplicityRanges[iMult]), static_cast<Int_t>(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<Double_t> 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<Int_t>(multiplicityRanges[iMult]), static_cast<Int_t>(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);
+}