PWG/muon/AnalysisFunctionOfMultiplicity.C

   1 /*************************************************************************************************************
   2
   3 This macro analyses the production of SingleMuon and J/Psi as a function of the collision multiplicity in the SPD.
   4 It reads and analyses the output of the Analysis Task PWG3/muon/AliAnalysisTaskMuonCollisionMultiplicity.
   5 Refer to the corresponding files to know what the output of this task is.
   6
   7 Thismacro use a number of input files whose names are hard-coded
   8 These are :
   9 pTRanges.txt : different ranges in pT in which the study is done.
  10 etaRnages.txt : different ranges in eta in which the study is done.
  11
  12
  13 *************************************************************************************************************/
  14
  15
  16
  17
  18 // ROOT includes
  19 #include <TH1D.h>
  20 #include <TH2D.h>
  21 #include <TH3D.h>
  22 #include <THnSparse.h>
  23 #include <TGraphErrors.h>
  24 #include <TGraphAsymmErrors.h>
  25 #include <TFile.h>
  26 #include <TList.h>
  27 #include <Riostream.h>
  28 #include <TMath.h>
  29
  30 // RooFit includes
  31 #include <RooRealVar.h>
  32 #include <RooAbsReal.h>
  33 #include <RooArgSet.h>
  34 #include <RooCBShape.h>
  35 #include <RooGaussian.h>
  36 #include <RooExponential.h>
  37 #include <RooAddPdf.h>
  38 #include <RooDataHist.h>
  39 #include <RooFitResult.h>
  40 #include <RooPlot.h>
  41
  42 // std includes
  43 #include <vector>
  44
  45
  46
  47
  48 // These functions are in charge of doing all the analysis
  49 void ProduceTriggerGraph(TFile *inputFile, TFile *outputFile, std::vector<Double_t> multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut);
  50
  51 void ProduceSingleMuonRawGraph(TFile *inputFile, TFile *outputFile, std::vector<Double_t> multiplicityRanges,
  52                                Bool_t applyZCut, Bool_t applyPileUpCut, Bool_t applyMatchTrigCut, Bool_t applyThetaAbsCut, Bool_t applyPDCACut);
  53
  54 void AnalysisSingleMuon(TFile *outputFile);
  55
  56 void SingleMuonYieldGraph(TFile *outputFile, TGraphErrors *CINT1B);
  57
  58 void SingleMuonYieldOverMeanMultGraph(TFile *outputFile, TGraphErrors *CINT1B);
  59
  60 void SingleMuonYieldNormalisedGraph(TFile *outputFile, TGraphErrors *CMUS1B);
  61
  62 void ProduceFitResults(TFile *inputFile, TFile *outputFile, std::vector<Double_t> multiplicityRanges,
  63                        Bool_t applyZCut, Bool_t applyPileUpCut, Double_t numberMatchTrig, Bool_t applyThetaAbsCut, Bool_t applyPDCACut);
  64
  65 void FitInvariantMass(TList *fitList, TH1D *invariantMass, Bool_t areParametersFixed, TH1D *referenceParameters);
  66
  67 void ProduceDimuonGraph(TFile *outputFile, std::vector<Double_t> multiplicityRanges);
  68
  69
  70 ///////////////////////////////////////////////////////////////////////////////////
  71 ///////////////////////////////////////////////////////////////////////////////////
  72 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")
  73 {
  74   // The main function of the macro
  75
  76   // Open the input file and create the output file
  77   TFile *inputFile = TFile::Open(inputName, "read");
  78   TFile *outputFile = TFile::Open(outputName, "recreate");
  79
  80   // First of all, we need to figure which conditions we are using
  81   // For the event : the cut on z_vertex and if we have pile up from the SPD
  82   // For the muons : the usual cuts for SMR
  83   Bool_t applyZCut = kTRUE;
  84   Bool_t applyPileUpCut = kTRUE;
  85   Bool_t applyMatchTrigCut = kTRUE;
  86   Bool_t applyThetaAbsCut = kTRUE;
  87   Bool_t applyPDCACut = kTRUE;
  88   Double_t nMatchTrig = 1.0;
  89
  90
  91   // Now, we need the multiplicity ranges
  92   // We get it from an input file
  93   ifstream multiplicityFile(multiplicityFileName);
  94   Double_t multiplicityRange = 0.0;
  95   std::vector <Double_t> multiplicityRanges;
  96   while(multiplicityFile >> multiplicityRange)
  97     multiplicityRanges.push_back(multiplicityRange);
  98
  99   // With that, we can produce the number of CINT1B in each bin
 100   // Specifficaly, we'll save two TGrpahErrors in the output :
 101   // - the graph containing the number of minimum bias event in the bin
 102   // - the graph containing the mean multiplicity in the bin
 103   ProduceTriggerGraph(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut);
 104
 105   // Now, the muon analysis
 106   if (doSingleMuonAnalysis)
 107     {
 108       // First, we produce the single muon raw graph
 109       // it is the graph giving the number of muons detected in different the multiplicity bins, and in different eta and pT domains
 110       ProduceSingleMuonRawGraph(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut, applyMatchTrigCut, applyThetaAbsCut, applyPDCACut);
 111
 112       // Raw graph are created, it is time to analyse them
 113       AnalysisSingleMuon(outputFile);
 114     }
 115
 116
 117   // Then, the J/Psi analysis
 118   if (doJPsiAnalysis)
 119     {
 120       ProduceFitResults(inputFile, outputFile, multiplicityRanges, applyZCut, applyPileUpCut, nMatchTrig, applyThetaAbsCut, applyPDCACut);
 121      ProduceDimuonGraph(outputFile, multiplicityRanges);
 122     }
 123
 124   inputFile->Close();
 125   outputFile->Close();
 126   cout << "End" << endl;
 127 }
 128
 129
 130
 131 ///////////////////////////////////////////////////////////////////////////////////
 132 ///////////////////////////////////////////////////////////////////////////////////
 133 void ProduceTriggerGraph(TFile *inputFile, TFile *outputFile, std::vector<Double_t> multiplicityRanges, Bool_t applyZCut, Bool_t applyPileUpCut)
 134 {
 135   // Compute the number of CINT1B and CMUS1B in each bin
 136   // 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
 137
 138
 139   // First of all, we need to get the number of CINT1B as a function of the multiplicity in a TH1D
 140   // Retrieving the THnSparse in the input
 141   THnSparseD *inputCINT1B = (THnSparseD *) ((TList *) (inputFile->Get("Trigger;1"))->FindObject("CINT1B") );
 142   THnSparseD *inputCMUS1B = (THnSparseD *) ((TList *) (inputFile->Get("Trigger;1"))->FindObject("CMUS1B") );
 143
 144
 145   // Reminder of the architecture of this THnSparse
 146   // dimension 0 : multiplicity of the event
 147   // dimension 1 : is the vertex in the z range (0 for no, 1 for yes)?
 148   // dimension 2 : is it an event without pile up (0 for no, 1 for yes)?
 149
 150   // Now, applying the cuts on z vertex and pile-up
 151   if (applyZCut)
 152     {
 153       inputCINT1B->GetAxis(1)->SetRangeUser(1.0, 2.0);
 154       inputCMUS1B->GetAxis(1)->SetRangeUser(1.0, 2.0);
 155     }
 156   if (applyPileUpCut)
 157     {
 158       inputCINT1B->GetAxis(2)->SetRangeUser(1.0, 2.0);
 159       inputCMUS1B->GetAxis(2)->SetRangeUser(1.0, 2.0);
 160     }
 161
 162   // cuts applied, we can project the THnSparse on TH1D, because it is easier to handle
 163   TH1D *histoCINT1B = inputCINT1B->Projection(0, "E");
 164   TH1D *histoCMUS1B = inputCMUS1B->Projection(0, "E");
 165
 166
 167   // We can now create and fill the two TGraphErrors
 168   TGraphErrors *graphCINT1B = new TGraphErrors(multiplicityRanges.size() - 1);
 169   graphCINT1B->SetName("graphCINT1B");
 170   TGraphErrors *graphCMUS1B = new TGraphErrors(multiplicityRanges.size() - 1);
 171   graphCMUS1B->SetName("graphCMUS1B");
 172
 173   // Graphs for plotting purposes
 174   TGraphAsymmErrors *plotCINT1B = new TGraphAsymmErrors(multiplicityRanges.size() - 1);
 175   plotCINT1B->SetName("plotCINT1B");
 176   TGraphAsymmErrors *plotCMUS1B = new TGraphAsymmErrors(multiplicityRanges.size() - 1);
 177   plotCMUS1B->SetName("plotCMUS1B");
 178
 179   for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++)
 180     {
 181       // CINT1B first
 182       Int_t firstBin = histoCINT1B->GetXaxis()->FindBin(multiplicityRanges[iRange]);
 183       Int_t lastBin = histoCINT1B->GetXaxis()->FindBin(multiplicityRanges[iRange+1]);
 184       histoCINT1B->GetXaxis()->SetRange(firstBin, lastBin);
 185
 186       Double_t total = 0.0;
 187       Double_t totalErr = 0.0;
 188       total = histoCINT1B->IntegralAndError(firstBin, lastBin, totalErr);
 189
 190       Double_t mean = histoCINT1B->GetMean();
 191       Double_t meanErr = histoCINT1B->GetMeanError();
 192
 193       graphCINT1B->SetPoint(iRange, mean, total);
 194       graphCINT1B->SetPointError(iRange, meanErr, totalErr);
 195
 196       plotCINT1B->SetPoint(iRange, mean, total);
 197       plotCINT1B->SetPointError(iRange,
 198                                 mean - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - mean,
 199                                 totalErr, totalErr);
 200
 201       // CMUS1B second
 202       firstBin = histoCMUS1B->GetXaxis()->FindBin(multiplicityRanges[iRange]);
 203       lastBin = histoCMUS1B->GetXaxis()->FindBin(multiplicityRanges[iRange+1]);
 204       histoCMUS1B->GetXaxis()->SetRange(firstBin, lastBin);
 205
 206       total = histoCMUS1B->IntegralAndError(firstBin, lastBin, totalErr);
 207
 208       mean = histoCMUS1B->GetMean();
 209       meanErr = histoCMUS1B->GetMeanError();
 210
 211       graphCMUS1B->SetPoint(iRange, mean, total);
 212       graphCMUS1B->SetPointError(iRange, meanErr, totalErr);
 213
 214       plotCMUS1B->SetPoint(iRange, mean, total);
 215       plotCMUS1B->SetPointError(iRange,
 216                                 mean - multiplicityRanges[iRange], multiplicityRanges[iRange+1] - mean,
 217                                 totalErr, totalErr);
 218     }
 219
 220   // Saving the graphs
 221   outputFile->WriteTObject(plotCINT1B);
 222   outputFile->WriteTObject(plotCMUS1B);
 223   outputFile->WriteTObject(graphCINT1B);
 224   outputFile->WriteTObject(graphCMUS1B);
 225
 226   delete graphCINT1B;
 227   delete histoCINT1B;
 228   delete inputCINT1B;
 229   delete graphCMUS1B;
 230   delete histoCMUS1B;
 231   delete inputCMUS1B;
 232 }
 233
 234
 235 void ProduceSingleMuonRawGraph(TFile *inputFile, TFile *outputFile, std::vector<Double_t> multiplicityRanges,
 236                                Bool_t applyZCut, Bool_t applyPileUpCut, Bool_t applyMatchTrigCut, Bool_t applyThetaAbsCut, Bool_t applyPDCACut)
 237 {
 238   // Produce the raw graph for Single Muon that will be used all along the analysis
 239   // There are several graph :
 240   //  - over all the range in both eta and pT
 241   //  - Single Muons Reference : pT > 1 GeV
 242   //  - Heavy Flavor Muons : pT > 4 GeV
 243   //  - bin by bin in pT and over all the range in eta
 244   //  - bin by bin in eta and with single muon reference
 245   //  and right now, there is not enough stat to do it bin by bin in both pT and eta
 246
 247
 248   // Retrieve the THnSparse and the Minimum bias histo
 249   THnSparseD *inputSingleMuon = (THnSparseD *) ((TList *) (inputFile->Get("SingleMuon;1"))->FindObject("muonCMUS1B") );
 250   TGraphErrors *CINT1B = (TGraphErrors *) (outputFile->Get("graphCINT1B;1") );
 251
 252   // Reminder of the architecture of this THnSparse
 253   // dimension 0  : multiplicity of the event
 254   // dimension 1  : is the vertex in the z range (0 for no, 1 for yes)?
 255   // dimension 2  : is it an event without pile up (0 for no, 1 for yes)?
 256   // dimension 3  : does the muon match the trigger (0 for no, 1 for yes)?
 257   // dimension 4  : theta_abs of the muon
 258   // dimension 5  : eta of the muon
 259   // dimension 6  : p DCA_x of the muon
 260   // dimension 7  : p DCA_y of the muon
 261   // dimension 8  : p DCA of the muon
 262   // dimension 9  : p of the muon
 263   // dimension 10 : pT of the muon
 264
 265
 266
 267   // get the pT and eta ranges from files.
 268   // Beware, the names of the files are hard-coded, so make sure to have them in your folder.
 269   // I might change this in the future, but I felt the main function had enough parameters already.
 270   ifstream pTFile ("pTRanges.txt");
 271   Double_t pTRange = 0.0;
 272   std::vector <Double_t> pTRanges;
 273   while(pTFile >> pTRange)
 274     pTRanges.push_back(pTRange);
 275
 276   ifstream etaFile ("etaRanges.txt");
 277   Double_t etaRange = 0.0;
 278   std::vector <Double_t> etaRanges;
 279   while(etaFile >> etaRange)
 280     etaRanges.push_back(etaRange);
 281
 282
 283   // Apply all the cuts
 284   // Beware, theta_abs and pDCA cuts are hard-coded
 285   if (applyZCut)
 286     inputSingleMuon->GetAxis(1)->SetRangeUser(1.0, 2.0);
 287   if (applyPileUpCut)
 288     inputSingleMuon->GetAxis(2)->SetRangeUser(1.0, 2.0);
 289   if (applyMatchTrigCut)
 290     inputSingleMuon->GetAxis(3)->SetRangeUser(1.0, 2.0);
 291   if (applyThetaAbsCut)
 292     inputSingleMuon->GetAxis(4)->SetRangeUser(2.0, 9.0);
 293   if (applyPDCACut)
 294     inputSingleMuon->GetAxis(8)->SetRangeUser(0.0, 450.0); // no cut yet, I need to check what are the usual values
 295
 296
 297   // First, we get all the raw histos
 298   TList *rawSingleMuon = new TList();
 299   rawSingleMuon->SetName("rawSingleMuon");
 300
 301   // All integrated
 302   // There are some hard cuts
 303   // - pT : 0.5 -> 8 GeV, because we are not confident on what we are seeing outside of these limits
 304   // - eta : -4.0, -> -2.5, acceptance of the spectrometer
 305
 306   Double_t etaMinAbsolute = -4.0;
 307   Double_t etaMaxAbsolute = -2.5;
 308   inputSingleMuon->GetAxis(5)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute);
 309
 310   Double_t pTMinAbsolute = 0.5;
 311   Double_t pTMaxAbsolute = 8.0;
 312   inputSingleMuon->GetAxis(10)->SetRangeUser(pTMinAbsolute, pTMaxAbsolute);
 313
 314   // Putting this in a TH3D, since it is easier to use and allow for computation of errors on the integral
 315   // TH3D is multiplicity, eta and pT
 316   TH3D *histoSingleMuon = inputSingleMuon->Projection(0, 5, 10, "E");
 317
 318   // Declare all the histos
 319   TGraphAsymmErrors *allSingleMuon = new TGraphAsymmErrors(multiplicityRanges.size() - 1);
 320   allSingleMuon->SetName("allSingleMuon");
 321   TGraphAsymmErrors *rawSMR = new TGraphAsymmErrors(multiplicityRanges.size() - 1); // Single Muon Reference (pT > 1.0 GeV)
 322   rawSMR->SetName("rawSMR");
 323   TGraphAsymmErrors *rawHFM = new TGraphAsymmErrors(multiplicityRanges.size() - 1); // Heavy Flavor Muon (pT > 4.0 GeV)
 324   rawHFM->SetName("rawHFM");
 325
 326
 327   for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++)
 328     {
 329       Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]);
 330       Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]);
 331       Int_t firstBinY = histoSingleMuon->GetYaxis()->GetFirst();
 332       Int_t lastBinY = histoSingleMuon->GetYaxis()->GetLast();
 333       Int_t firstBinZ = histoSingleMuon->GetZaxis()->GetFirst();
 334       Int_t lastBinZ = histoSingleMuon->GetZaxis()->GetLast();
 335
 336       Double_t nSingleMuon = 0.0;
 337       Double_t errSingleMuon = 0.0;
 338       nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon);
 339
 340       allSingleMuon->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon);
 341       allSingleMuon->SetPointError(iRange,
 342                                    CINT1B->GetX()[iRange] - multiplicityRanges[iRange],
 343                                    multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange],
 344                                    errSingleMuon,
 345                                    errSingleMuon);
 346
 347       // SMR
 348       firstBinZ = histoSingleMuon->GetZaxis()->FindBin(1.0);
 349       lastBinZ = histoSingleMuon->GetZaxis()->GetLast();
 350       nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon);
 351
 352       rawSMR->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon);
 353       rawSMR->SetPointError(iRange,
 354                             CINT1B->GetX()[iRange] - multiplicityRanges[iRange],
 355                             multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange],
 356                             errSingleMuon,
 357                             errSingleMuon);
 358
 359       // HFM
 360       firstBinZ = histoSingleMuon->GetZaxis()->FindBin(4.0);
 361       lastBinZ = histoSingleMuon->GetZaxis()->GetLast();
 362       nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon);
 363
 364       rawHFM->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon);
 365       rawHFM->SetPointError(iRange,
 366                             CINT1B->GetX()[iRange] - multiplicityRanges[iRange],
 367                             multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange],
 368                             errSingleMuon,
 369                             errSingleMuon);
 370     }
 371
 372   TList *integratedSingleMuon = new TList();
 373   integratedSingleMuon->SetName("rawSingleMuonIntegrated");
 374
 375   integratedSingleMuon->AddAt(allSingleMuon, 0);
 376   integratedSingleMuon->AddAt(rawSMR, 1);
 377   integratedSingleMuon->AddAt(rawHFM, 2);
 378
 379   rawSingleMuon->AddAt(integratedSingleMuon, 0);
 380   //delete allSingleMuon;
 381   //delete rawSMR;
 382   //delete rawHFM;
 383
 384   // Now, for the study in pT
 385   TList *rawSingleMuonPt = new TList();
 386   rawSingleMuonPt->SetName("rawSingleMuonPt");
 387
 388   for (UInt_t ipT = 0; ipT < pTRanges.size()-1; ipT++)
 389     {
 390       TString nameGraph;
 391       nameGraph.Form("pTRange_%1f-%1f", pTRanges[ipT], pTRanges[ipT+1]);
 392       TGraphAsymmErrors *singleMuonPtRange = new TGraphAsymmErrors();
 393       singleMuonPtRange->SetName(nameGraph);
 394
 395       for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++)
 396         {
 397           Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]);
 398           Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]);
 399           Int_t firstBinY = histoSingleMuon->GetYaxis()->GetFirst();
 400           Int_t lastBinY = histoSingleMuon->GetYaxis()->GetLast();
 401           Int_t firstBinZ = histoSingleMuon->GetZaxis()->FindBin(pTRanges[ipT]);
 402           Int_t lastBinZ = histoSingleMuon->GetZaxis()->FindBin(pTRanges[ipT]+1);
 403
 404           Double_t nSingleMuon = 0.0;
 405           Double_t errSingleMuon = 0.0;
 406           nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon);
 407
 408           singleMuonPtRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon);
 409           singleMuonPtRange->SetPointError(iRange,
 410                                            CINT1B->GetX()[iRange] - multiplicityRanges[iRange],
 411                                            multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange],
 412                                            errSingleMuon,
 413                                            errSingleMuon);
 414
 415         }
 416
 417       rawSingleMuonPt->AddAt(singleMuonPtRange, ipT);
 418       //delete singleMuonPtRange;
 419     }
 420
 421   rawSingleMuon->AddAt(rawSingleMuonPt, 1);
 422
 423   // At last, the study in eta
 424   // Both on SMR and HFM
 425   TList *rawSingleMuonEta = new TList();
 426   rawSingleMuonEta->SetName("rawSingleMuonEta");
 427
 428   for (UInt_t iEta = 0; iEta < etaRanges.size()-1; iEta++)
 429     {
 430       TString nameGraph;
 431       nameGraph.Form("SMRetaRange_%1f-%1f", etaRanges[iEta], etaRanges[iEta+1]);
 432       TGraphAsymmErrors *SMREtaRange = new TGraphAsymmErrors();
 433       SMREtaRange->SetName(nameGraph);
 434
 435       nameGraph.Form("HFMetaRange_%1f-%1f", etaRanges[iEta], etaRanges[iEta+1]);
 436       TGraphAsymmErrors *HFMEtaRange = new TGraphAsymmErrors();
 437       HFMEtaRange->SetName(nameGraph);
 438
 439       for (UInt_t iRange = 0; iRange < multiplicityRanges.size() - 1; iRange++)
 440         {
 441           // First, SMR
 442           Int_t firstBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange]);
 443           Int_t lastBinX = histoSingleMuon->GetXaxis()->FindBin(multiplicityRanges[iRange+1]);
 444           Int_t firstBinY = histoSingleMuon->GetYaxis()->FindBin(etaRanges[iEta]);
 445           Int_t lastBinY = histoSingleMuon->GetYaxis()->FindBin(etaRanges[iEta]+1);
 446           Int_t firstBinZ = histoSingleMuon->GetZaxis()->FindBin(1.0);
 447           Int_t lastBinZ = histoSingleMuon->GetZaxis()->GetLast();
 448
 449           Double_t nSingleMuon = 0.0;
 450           Double_t errSingleMuon = 0.0;
 451           nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon);
 452
 453           SMREtaRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon);
 454           SMREtaRange->SetPointError(iRange,
 455                                      CINT1B->GetX()[iRange] - multiplicityRanges[iRange],
 456                                      multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange],
 457                                      errSingleMuon,
 458                                      errSingleMuon);
 459
 460           // Then, HFM
 461           firstBinZ = histoSingleMuon->GetZaxis()->FindBin(4.0);
 462           lastBinZ = histoSingleMuon->GetZaxis()->GetLast();
 463
 464           nSingleMuon = histoSingleMuon->IntegralAndError(firstBinX, lastBinX, firstBinY, lastBinY, firstBinZ, lastBinZ, errSingleMuon);
 465
 466           HFMEtaRange->SetPoint(iRange, CINT1B->GetX()[iRange], nSingleMuon);
 467           HFMEtaRange->SetPointError(iRange,
 468                                      CINT1B->GetX()[iRange] - multiplicityRanges[iRange],
 469                                      multiplicityRanges[iRange+1] - CINT1B->GetX()[iRange],
 470                                      errSingleMuon,
 471                                      errSingleMuon);
 472
 473
 474         }
 475
 476       rawSingleMuonEta->AddAt(SMREtaRange, iEta);
 477       rawSingleMuonEta->AddAt(HFMEtaRange, iEta + etaRanges.size()-1);
 478       //delete SMREtaRange;
 479       //delete HFMEtaRange;
 480     }
 481
 482   rawSingleMuon->AddAt(rawSingleMuonEta, 2);
 483
 484   outputFile->WriteTObject(rawSingleMuon);
 485   delete integratedSingleMuon;
 486   delete rawSingleMuonPt;
 487   delete rawSingleMuonEta;
 488 }
 489
 490
 491 void AnalysisSingleMuon(TFile *outputFile)
 492 {
 493   // Analyse the raw data for single muons
 494
 495   // First, we retrieve the CINT1B graph, as usual
 496   TGraphErrors *CINT1B = (TGraphErrors *) outputFile->Get("graphCINT1B;1");
 497   TGraphErrors *CMUS1B = (TGraphErrors *) outputFile->Get("graphCMUS1B;1");
 498
 499   // Creation of all the yield graphs
 500   SingleMuonYieldGraph(outputFile, CINT1B);
 501
 502   // Creation of the yield over mean mult graph
 503   SingleMuonYieldOverMeanMultGraph(outputFile, CINT1B);
 504
 505   // Creation of the yield with the reference bin normalised to 1, for comparison
 506   SingleMuonYieldNormalisedGraph(outputFile, CMUS1B);
 507 }
 508
 509
 510
 511
 512 //////////////////////////////////////////////////////////////////////////////////////////////////
 513 //////////////////////////////////////////////////////////////////////////////////////////////////
 514 void SingleMuonYieldGraph(TFile *outputFile, TGraphErrors *CINT1B)
 515 {
 516   // This function use the raw graphs of outputfile and the CINT1B graph to create the yield graph
 517   // yield are raw divided by CINT1B
 518   TList *yieldSingleMuon = new TList();
 519   yieldSingleMuon->SetName("yieldSingleMuon");
 520
 521   for (Int_t iList = 0; iList < ((TList *) outputFile->Get("rawSingleMuon;1"))->GetEntries(); iList++)
 522     {
 523       TList *currentList = (TList *) (( TList *) outputFile->Get("rawSingleMuon;1"))->At(iList);
 524       TList *newList = new TList();
 525       TString listName = currentList->GetName();
 526       listName.ReplaceAll("raw", 3, "yield", 5);
 527       newList->SetName(listName);
 528
 529       for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++)
 530         {
 531           TGraphAsymmErrors *rawGraph = (TGraphAsymmErrors *) currentList->At(iGraph);
 532
 533           TString yieldName = rawGraph->GetName();
 534           yieldName.ReplaceAll("raw", 3, "yield", 5);
 535           TGraphAsymmErrors *yieldGraph = new TGraphAsymmErrors(rawGraph->GetN());
 536           yieldGraph->SetName(yieldName);
 537
 538           // fill the yield graph
 539           for (Int_t iBin = 0; iBin < yieldGraph->GetN(); iBin++ )
 540             {
 541               if (CINT1B->GetY()[iBin] != 0.0)
 542                 {
 543                   yieldGraph->SetPoint(iBin,
 544                                        rawGraph->GetX()[iBin],
 545                                        rawGraph->GetY()[iBin]/CINT1B->GetY()[iBin]);
 546                   if (rawGraph->GetY()[iBin] != 0.0)
 547                     {
 548                       Double_t error = yieldGraph->GetY()[iBin] *
 549                         TMath::Sqrt((rawGraph->GetEYlow()[iBin]/rawGraph->GetY()[iBin])*(rawGraph->GetEYlow()[iBin]/rawGraph->GetY()[iBin]) +
 550                                     (CINT1B->GetEY()[iBin]/CINT1B->GetY()[iBin])*(CINT1B->GetEY()[iBin]/CINT1B->GetY()[iBin]));
 551                       yieldGraph->SetPointError(iBin,
 552                                               rawGraph->GetEXlow()[iBin],
 553                                               rawGraph->GetEXhigh()[iBin],
 554                                               error, error);
 555                     }
 556                 }
 557
 558               else
 559                 yieldGraph->SetPoint(iBin, 0, 0);
 560             }
 561
 562           newList->AddAt(yieldGraph, iGraph);
 563         }
 564
 565       yieldSingleMuon->AddAt(newList, iList);
 566     }
 567
 568   outputFile->WriteTObject(yieldSingleMuon);
 569 }
 570
 571
 572
 573 //////////////////////////////////////////////////////////////////////////////////////////////////
 574 //////////////////////////////////////////////////////////////////////////////////////////////////
 575 void SingleMuonYieldOverMeanMultGraph(TFile *outputFile, TGraphErrors *CINT1B)
 576 {
 577   // This function use the yield graphs of outputfile to create the yield graph divided by the mean multiplicity in each bin
 578   // The CINT1B graph is necessary to get the error on the mean multiplicity
 579   TList *yieldOverMeanMultSingleMuon = new TList();
 580   yieldOverMeanMultSingleMuon->SetName("yieldOverMeanMultSingleMuon");
 581
 582   for (Int_t iList = 0; iList < ((TList *) outputFile->Get("yieldSingleMuon;1"))->GetEntries(); iList++)
 583     {
 584       TList *currentList = (TList *) (( TList *) outputFile->Get("yieldSingleMuon;1"))->At(iList);
 585       TList *newList = new TList();
 586       TString listName = currentList->GetName();
 587       listName.ReplaceAll("yield", 5, "yieldOverMeanMult", 17);
 588       newList->SetName(listName);
 589
 590       for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++)
 591         {
 592           TGraphAsymmErrors *yieldGraph = (TGraphAsymmErrors *) currentList->At(iGraph);
 593
 594           TString yieldOverMeanMultName = yieldGraph->GetName();
 595           yieldOverMeanMultName.ReplaceAll("yield", 5, "yieldOverMeanMult", 17);
 596           TGraphAsymmErrors *yieldOverMeanMultGraph = new TGraphAsymmErrors(yieldGraph->GetN());
 597           yieldOverMeanMultGraph->SetName(yieldOverMeanMultName);
 598
 599           // fill the yield graph
 600           for (Int_t iBin = 0; iBin < yieldOverMeanMultGraph->GetN(); iBin++ )
 601             {
 602               if (CINT1B->GetX()[iBin] != 0.0)
 603                 {
 604                   yieldOverMeanMultGraph->SetPoint(iBin,
 605                                                    yieldGraph->GetX()[iBin],
 606                                                    yieldGraph->GetY()[iBin]/CINT1B->GetX()[iBin]);
 607                   if (yieldGraph->GetY()[iBin] != 0.0)
 608                     {
 609                       Double_t error = yieldOverMeanMultGraph->GetY()[iBin]*
 610                         TMath::Sqrt((yieldGraph->GetEYlow()[iBin]/yieldGraph->GetY()[iBin])*(yieldGraph->GetEYlow()[iBin]/yieldGraph->GetY()[iBin]) +
 611                                     (CINT1B->GetEX()[iBin]/CINT1B->GetX()[iBin])*(CINT1B->GetEX()[iBin]/CINT1B->GetX()[iBin]));
 612                       yieldOverMeanMultGraph->SetPointError(iBin,
 613                                                             yieldGraph->GetEXlow()[iBin],
 614                                                             yieldGraph->GetEXhigh()[iBin],
 615                                                             error, error);
 616                     }
 617                 }
 618
 619               else
 620                 yieldOverMeanMultGraph->SetPoint(iBin, 0, 0);
 621             }
 622
 623           newList->AddAt(yieldOverMeanMultGraph, iGraph);
 624         }
 625
 626       yieldOverMeanMultSingleMuon->AddAt(newList, iList);
 627     }
 628
 629   outputFile->WriteTObject(yieldOverMeanMultSingleMuon);
 630 }
 631
 632
 633 //////////////////////////////////////////////////////////////////////////////////////////////////
 634 //////////////////////////////////////////////////////////////////////////////////////////////////
 635 void SingleMuonYieldNormalisedGraph(TFile *outputFile, TGraphErrors *CMUS1B)
 636 {
 637   // This function normalise the yield over mean multiplicity so that the reference bin in multiplicity is equal to 1.
 638   // The reference bin is defined it as the bin with the maximum number of CMUS1B
 639   TList *yieldNormalisedSingleMuon = new TList();
 640   yieldNormalisedSingleMuon->SetName("yieldNormalisedSingleMuon");
 641
 642
 643   // The first step is to find the reference bin in multiplicity
 644   // we define it as the bin with the maximum number of CMUS1B
 645
 646   Double_t maxCMUS1B = 0.0;
 647   Int_t referenceBin = 0;
 648   for (Int_t iBin = 0; iBin < CMUS1B->GetN(); iBin++)
 649     if (CMUS1B->GetX()[iBin] > 10.0)
 650       if (CMUS1B->GetY()[iBin] > maxCMUS1B)
 651         {
 652           maxCMUS1B = CMUS1B->GetY()[iBin];
 653           referenceBin = iBin;
 654         }
 655
 656   // Now the loop, as usual
 657   // We don't propagate the error on the normalisation factor, since this is artificial
 658   for (Int_t iList = 0; iList < ((TList *) outputFile->Get("yieldOverMeanMultSingleMuon;1"))->GetEntries(); iList++)
 659     {
 660       TList *currentList = (TList *) (( TList *) outputFile->Get("yieldOverMeanMultSingleMuon;1"))->At(iList);
 661       TList *newList = new TList();
 662       TString listName = currentList->GetName();
 663       listName.ReplaceAll("yieldOverMeanMult", 17, "yieldNormalised", 15);
 664       newList->SetName(listName);
 665
 666       for (Int_t iGraph = 0; iGraph < currentList->GetEntries(); iGraph++)
 667         {
 668           TGraphAsymmErrors *yieldOverMeanMultGraph = (TGraphAsymmErrors *) currentList->At(iGraph);
 669
 670           TString yieldNormalisedName = yieldOverMeanMultGraph->GetName();
 671           yieldNormalisedName.ReplaceAll("yieldOverMeanMult", 17, "yieldNormalised", 15);
 672           TGraphAsymmErrors *yieldNormalisedGraph = new TGraphAsymmErrors(yieldOverMeanMultGraph->GetN());
 673           yieldNormalisedGraph->SetName(yieldNormalisedName);
 674
 675           // fill the yield graph
 676           for (Int_t iBin = 0; iBin < yieldNormalisedGraph->GetN(); iBin++ )
 677             {
 678               if (yieldOverMeanMultGraph->GetY()[iBin] != 0.0)
 679                 {
 680                   yieldNormalisedGraph->SetPoint(iBin,
 681                                                  yieldOverMeanMultGraph->GetX()[iBin],
 682                                                  yieldOverMeanMultGraph->GetY()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin]);
 683                   if (yieldOverMeanMultGraph->GetY()[iBin] != 0.0)
 684                     yieldNormalisedGraph->SetPointError(iBin,
 685                                                         yieldOverMeanMultGraph->GetEXlow()[iBin],
 686                                                         yieldOverMeanMultGraph->GetEXhigh()[iBin],
 687                                                         yieldOverMeanMultGraph->GetEYlow()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin],
 688                                                         yieldOverMeanMultGraph->GetEYhigh()[iBin]/yieldOverMeanMultGraph->GetY()[referenceBin]);
 689                 }
 690
 691               else
 692                 yieldNormalisedGraph->SetPoint(iBin, 0, 0);
 693             }
 694
 695           newList->AddAt(yieldNormalisedGraph, iGraph);
 696         }
 697
 698       yieldNormalisedSingleMuon->AddAt(newList, iList);
 699     }
 700
 701   outputFile->WriteTObject(yieldNormalisedSingleMuon);
 702 }
 703
 704
 705 /////////////////////////////////////////////////////////////////////////////////////////////////
 706 /////////////////////////////////////////////////////////////////////////////////////////////////
 707 void ProduceFitResults(TFile *inputFile, TFile *outputFile, std::vector<Double_t> multiplicityRanges,
 708                        Bool_t applyZCut, Bool_t applyPileUpCut, Double_t numberMatchTrig, Bool_t applyThetaAbsCut, Bool_t applyPDCACut)
 709 {
 710   // Produce the raw graph for Dimuons that will be used all along the analysis
 711   // There are two graphs :
 712   //  - J/Psi
 713   //  - Background in the J/Psi range
 714
 715   // Retrieve the THnSparse and the Minimum bias histo
 716   THnSparseD *inputDimuon = (THnSparseD *) ((TList *) (inputFile->Get("Dimuon;1"))->FindObject("dimuonCMUS1B") );
 717
 718   // Reminder of the structure of the dimuon histo
 719   // dimension 0  : multiplicity of the event
 720   // dimension 1  : is the vertex in the z range (0 for no, 1 for yes)?
 721   // dimension 2  : is it an event without pile up (0 for no, 1 for yes)?
 722   // dimension 3  : does the first muon match the trigger (0 for no, 1 for yes)?
 723   // dimension 4  : does the second muon match the trigger (0 for no, 1 for yes)?
 724   // dimension 5  : number of muons matching the trigger in the dimuon
 725   // dimension 6  : theta_abs of the first muon
 726   // dimension 7  : theta_abs of the second muon
 727   // dimension 8  : eta of the first muon
 728   // dimension 9  : eta of the second muon
 729   // dimension 10 : p DCA_x of the first muon
 730   // dimension 11 : p DCA_y of the first muon
 731   // dimension 12 : p DCA_x of the second muon
 732   // dimension 13 : p DCA_y of the second muon
 733   // dimension 14 : p of the first muon
 734   // dimension 15 : p of the second muon
 735   // dimension 16 : p of the dimuon
 736   // dimension 17 : pT of the dimuon
 737   // dimension 18 : invariant mass of the dimuon
 738
 739
 740
 741   // get the pT and eta ranges from files.
 742   // Beware, the names of the files are hard-coded, so make sure to have them in your folder.
 743   // I might change this in the future, but I felt the main function had enough parameters already.
 744   ifstream pTFile ("pTRanges.txt");
 745   Double_t pTRange = 0.0;
 746   std::vector <Double_t> pTRanges;
 747   while(pTFile >> pTRange)
 748     pTRanges.push_back(pTRange);
 749
 750   ifstream etaFile ("etaRanges.txt");
 751   Double_t etaRange = 0.0;
 752   std::vector <Double_t> etaRanges;
 753   while(etaFile >> etaRange)
 754     etaRanges.push_back(etaRange);
 755
 756
 757   // Apply all the cuts
 758   // Beware, theta_abs and pDCA cuts are hard-coded
 759   if (applyZCut)
 760     inputDimuon->GetAxis(1)->SetRangeUser(1.0, 2.0);
 761   if (applyPileUpCut)
 762     inputDimuon->GetAxis(2)->SetRangeUser(1.0, 2.0);
 763   inputDimuon->GetAxis(5)->SetRangeUser(numberMatchTrig, 3.0);
 764   if (applyThetaAbsCut)
 765     {
 766       inputDimuon->GetAxis(6)->SetRangeUser(2.0, 9.0);
 767       inputDimuon->GetAxis(7)->SetRangeUser(2.0, 9.0);
 768     }
 769   if (applyPDCACut)
 770     {// no cut yet, I need to check what are the usual values
 771       inputDimuon->GetAxis(12)->SetRangeUser(0.0, 450.0);
 772       inputDimuon->GetAxis(15)->SetRangeUser(0.0, 450.0);
 773     }
 774
 775
 776   // First, we get the invariant mass histos
 777   TList *dimuonFitResults = new TList();
 778   dimuonFitResults->SetName("dimuonFitResults");
 779
 780   // There are some hard cuts
 781   // - eta : -4.0, -> -2.5, acceptance of the spectrometer
 782
 783   Double_t etaMinAbsolute = -4.0;
 784   Double_t etaMaxAbsolute = -2.5;
 785   inputDimuon->GetAxis(8)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute);
 786   inputDimuon->GetAxis(9)->SetRangeUser(etaMinAbsolute, etaMaxAbsolute);
 787
 788   // First, we get the invariant mass histo for the whole range in multiplicity
 789   inputDimuon->GetAxis(0)->SetRangeUser(multiplicityRanges[0], multiplicityRanges[multiplicityRanges.size()-1]);
 790   TH1D *invariantMassIntegrated = inputDimuon->Projection(18, "E");
 791   invariantMassIntegrated->SetName("invariantMassIntegrated");
 792
 793
 794   // Now for the fit on all the multiplicity
 795   // This is used to get a value for some parameters
 796   TList *fitAll = new TList();
 797   fitAll->SetName("AllMultiplicity");
 798
 799
 800   // Create the container for the reference parameters
 801   TH1D *referenceParameters = new TH1D("referenceParameters", "parameters", 8, 1.0, 9.0);
 802   referenceParameters->Sumw2();
 803   referenceParameters->GetXaxis()->SetBinLabel(1, "meanJPsi");
 804   referenceParameters->GetXaxis()->SetBinLabel(2, "sigmaJPsi");
 805   referenceParameters->GetXaxis()->SetBinLabel(3, "alphaJPsi");
 806   referenceParameters->GetXaxis()->SetBinLabel(4, "nJPsi");
 807   referenceParameters->GetXaxis()->SetBinLabel(5, "meanPsiPrime");
 808   referenceParameters->GetXaxis()->SetBinLabel(6, "sigmaPsiPrime");
 809   referenceParameters->GetXaxis()->SetBinLabel(7, "expBkg1");
 810   referenceParameters->GetXaxis()->SetBinLabel(8, "expBkg2");
 811
 812
 813   FitInvariantMass(fitAll, invariantMassIntegrated, kFALSE, referenceParameters);
 814   dimuonFitResults->AddAt(fitAll, 0);
 815
 816   // Now, we make a fit for each range in multiplicity
 817   for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++)
 818     {
 819       TString name;
 820       name.Form("Multiplicity_%d-%d", static_cast<Int_t>(multiplicityRanges[iMult]), static_cast<Int_t>(multiplicityRanges[iMult+1]));
 821       TList *fitRange = new TList();
 822       fitRange->SetName(name);
 823
 824       inputDimuon->GetAxis(0)->SetRangeUser(multiplicityRanges[iMult], multiplicityRanges[iMult+1]);
 825       TH1D *invariantMassRange = inputDimuon->Projection(18, "E");
 826       invariantMassRange->SetName("invariantMass");
 827
 828       FitInvariantMass(fitRange, invariantMassRange, kTRUE, referenceParameters);
 829       dimuonFitResults->AddAt(fitRange, iMult + 1);
 830     }
 831
 832   outputFile->WriteTObject(dimuonFitResults);
 833 }
 834
 835
 836
 837
 838 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 839 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 840 void FitInvariantMass(TList *fitList, TH1D *invariantMass, Bool_t areParametersFixed, TH1D *referenceParameters)
 841 {
 842   // Fit the given histo and put the results in outputFile
 843   // The fit is : Crystal Ball (J/Psi) + gaussian (Psi') + two exponentials (background)
 844
 845   // The parameters histo architecture is the following :
 846   // bin 1 : mean J/Psi
 847   // bin 2 : sigma J/Psi
 848   // bin 3 : alpha J/Psi
 849   // bin 4 : n J/Psi
 850   // bin 5 : mean Psi'
 851   // bin 6 : sigma Psi'
 852   // bin 7 : exponential factor background 1
 853   // bin 8 : exponential factor background 2
 854
 855   // The results histo architecture is the following :
 856   // bin 1 : J/Psi
 857   // bin 2 : Background J/Psi at 2 sigma
 858   // bin 3 : signal over background
 859   // bin 4 : chi2
 860
 861   // There are some parameters that are not always free :
 862   // mean of J/Psi peak
 863   // sigma of J/Psi peak
 864   // alpha of J/Psi function (Crystall Ball function)
 865   // n of J/Psi function
 866   // mean of Psi' peak
 867   // sigma of Psi' peak
 868   // If there are fixed, their values are taken from the parameters histo
 869
 870
 871   // First, we have to create two histos, one that will contain the parameters of the fit, and the other the results
 872   TH1D *parameters = new TH1D("parameters", "parameters", 8, 1.0, 9.0);
 873   parameters->Sumw2();
 874   parameters->GetXaxis()->SetBinLabel(1, "meanJPsi");
 875   parameters->GetXaxis()->SetBinLabel(2, "sigmaJPsi");
 876   parameters->GetXaxis()->SetBinLabel(3, "alphaJPsi");
 877   parameters->GetXaxis()->SetBinLabel(4, "nJPsi");
 878   parameters->GetXaxis()->SetBinLabel(5, "meanPsiPrime");
 879   parameters->GetXaxis()->SetBinLabel(6, "sigmaPsiPrime");
 880   parameters->GetXaxis()->SetBinLabel(7, "expBkg1");
 881   parameters->GetXaxis()->SetBinLabel(8, "expBkg2");
 882   TH1D *results = new TH1D("results", "results", 4, 1.0, 5.0);
 883   results->Sumw2();
 884   results->GetXaxis()->SetBinLabel(1, "JPsi");
 885   results->GetXaxis()->SetBinLabel(2, "bkg");
 886   results->GetXaxis()->SetBinLabel(3, "SB");
 887   results->GetXaxis()->SetBinLabel(4, "chi2");
 888
 889
 890   // Declare observable M
 891   RooRealVar M("M","Dimuon invariant mass [GeV/c^2]", 2.0, 5.0);
 892
 893
 894   // Declare Cristal Ball for J/Psi
 895   RooRealVar mean_jpsi("mean_jpsi","mean_jpsi", 3.096, 2.8, 3.2);
 896   RooRealVar sigma_jpsi("sigma_jpsi","sigma_jpsi", 0.070, 0, 0.2);
 897   RooRealVar alpha_jpsi("alpha","alpha", 5, 0, 10);
 898   RooRealVar n_jpsi("n","n",5,0,10);
 899   RooCBShape jPsi("jpsi", "crystal ball PDF", M, mean_jpsi, sigma_jpsi, alpha_jpsi, n_jpsi);
 900
 901
 902   //Declare gaussian for Psi prime
 903   RooRealVar mean_psip("mean_psip", "mean_psip", 3.686, 3.6, 3.75);
 904   RooRealVar sigma_psip("sigma_psip","sigma_psip", 0.086);
 905   RooGaussian psiPrime("psip","Gaussian", M, mean_psip, sigma_psip);
 906
 907   //Declare exponentials for background
 908   RooRealVar c1("c1", "c1", -10, -0.1);
 909   RooExponential bkg_exp1("bkg_exp1", "background", M, c1);
 910
 911   RooRealVar c2("c2", "c2", -5.0, -1);
 912   RooExponential bkg_exp2("bkg_exp2", "background", M, c2);
 913
 914   // Give a value and fix sigma, alpha and n in the Cristal Ball
 915   // The values correspond to a previous fit (ideally, with all the statistic)
 916   // I should automatize this one of these day
 917   if (areParametersFixed)
 918     {
 919       Double_t fixMeanJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("meanJPsi"));
 920       mean_jpsi.setVal(fixMeanJPsi);
 921       mean_jpsi.setError(0.0);
 922       mean_jpsi.setRange(fixMeanJPsi, fixMeanJPsi);
 923
 924       Double_t fixSigmaJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi"));
 925       sigma_jpsi.setVal(fixSigmaJPsi);
 926       sigma_jpsi.setError(0.0);
 927       sigma_jpsi.setRange(fixSigmaJPsi, fixSigmaJPsi);
 928
 929       Double_t fixAlphaJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi"));
 930       alpha_jpsi.setVal(fixAlphaJPsi);
 931       alpha_jpsi.setError(0.0);
 932       alpha_jpsi.setRange(fixAlphaJPsi, fixAlphaJPsi);
 933
 934       Double_t fixNJPsi = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("nJPsi"));
 935       n_jpsi.setVal(fixNJPsi);
 936       n_jpsi.setError(0.0);
 937       n_jpsi.setRange(fixNJPsi, fixNJPsi);
 938
 939       Double_t fixMeanPsiPrime = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime"));
 940       mean_psip.setVal(fixMeanPsiPrime);
 941       mean_psip.setError(0.0);
 942       mean_psip.setRange(fixMeanPsiPrime, fixMeanPsiPrime);
 943
 944       Double_t fixSigmaPsiPrime = referenceParameters->GetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime"));
 945       sigma_psip.setVal(fixSigmaPsiPrime);
 946       sigma_psip.setError(0.0);
 947       sigma_psip.setRange(fixSigmaPsiPrime, fixSigmaPsiPrime);
 948     }
 949
 950
 951   // Sum the composite signal and background into an extended pdf nsig*sig+nbkg*bkg
 952   RooRealVar fitJPsi("fitJPsi", "number of J/Psi events", 1600 ,0.0, 15000);
 953   RooRealVar fitPsiPrime("fitPsiPrime", "number of Psi Prime events", 50, 0.0, 200);
 954   RooRealVar fitBkg1("fitBkg1", "number of background events", 5000, 0.0, 15000);
 955   RooRealVar fitBkg2("fitBkg2", "number of background events", 5000, 0.0, 15000);
 956
 957
 958   RooAddPdf *fitFunction = new RooAddPdf("model", "(CB+Gauss+2exp)", RooArgList(bkg_exp1, bkg_exp2, jPsi, psiPrime), RooArgList(fitBkg1, fitBkg2, fitJPsi, fitPsiPrime));
 959
 960   // Define the histo to be fitted
 961   RooDataHist *invariantMassFit = new RooDataHist("invariantMassFit", "invariantMassFit", M, invariantMass);
 962
 963   // Fit the invariant mass
 964   // This is the important part
 965   RooFitResult *fitResult = fitFunction->fitTo(*invariantMassFit, RooFit::Save());
 966   //RooFitResult *fitResult = fitFunction->chi2FitTo(*invariantMassFit, RooFit::Save());
 967
 968   RooPlot *plot = M.frame();
 969   TString title = "fittedInvariantMass";
 970   plot->SetTitle(title);
 971   plot->SetName(title);
 972
 973   invariantMassFit->plotOn(plot, RooFit::Name("invariantMassFit"));
 974   fitFunction->plotOn(plot, RooFit::Name("fitFunction"), RooFit::Range(2.0, 5.0));
 975
 976
 977   // Fill the histo with the final values of the parameters
 978   parameters->SetBinContent(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getVal());
 979   parameters->SetBinError(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getError());
 980
 981   parameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getVal());
 982   parameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getError());
 983
 984   parameters->SetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getVal());
 985   parameters->SetBinError(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getError());
 986
 987   parameters->SetBinContent(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getVal());
 988   parameters->SetBinError(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getError());
 989
 990   parameters->SetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getVal());
 991   parameters->SetBinError(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getError());
 992
 993   parameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getVal());
 994   parameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getError());
 995
 996   // Fill the histo with the reference parameters
 997   if (!areParametersFixed)
 998     {
 999       referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getVal());
1000       referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("meanJPsi"), mean_jpsi.getError());
1001
1002       referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getVal());
1003       referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaJPsi"), sigma_jpsi.getError());
1004
1005       referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getVal());
1006       referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("alphaJPsi"), alpha_jpsi.getError());
1007
1008       referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getVal());
1009       referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("nJPsi"), n_jpsi.getError());
1010
1011       referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getVal());
1012       referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("meanPsiPrime"), mean_psip.getError());
1013
1014       referenceParameters->SetBinContent(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getVal());
1015       referenceParameters->SetBinError(parameters->GetXaxis()->FindBin("sigmaPsiPrime"), sigma_psip.getError());
1016     }
1017
1018   parameters->SetBinContent(parameters->GetXaxis()->FindBin("expBkg1"), c1.getVal());
1019   parameters->SetBinError(parameters->GetXaxis()->FindBin("expBkg1"), c1.getError());
1020
1021   parameters->SetBinContent(parameters->GetXaxis()->FindBin("expBkg2"), c2.getVal());
1022   parameters->SetBinError(parameters->GetXaxis()->FindBin("expBkg2"), c2.getError());
1023
1024
1025   // Define the range to compute the results
1026   // We choose two sigma
1027   // We also need two new variables to compute the background
1028   Double_t lowerBoundJPsi = mean_jpsi.getVal() - 2.0*sigma_jpsi.getVal();
1029   Double_t upperBoundJPsi = mean_jpsi.getVal() + 2.0*sigma_jpsi.getVal();
1030   M.setRange("JPsiRange" ,lowerBoundJPsi, upperBoundJPsi);
1031
1032   RooAbsReal *fracBkgRange1 = bkg_exp1.createIntegral(M, M, "JPsiRange");
1033   RooAbsReal *fracBkgRange2 = bkg_exp2.createIntegral(M, M, "JPsiRange");
1034
1035   // Get the results, and fill the results histo
1036   Double_t nJPsi = fitJPsi.getVal();
1037   Double_t errJPsi = fitJPsi.getError();
1038   results->SetBinContent(results->GetXaxis()->FindBin("JPsi"), nJPsi);
1039   results->SetBinError(results->GetXaxis()->FindBin("JPsi"), errJPsi);
1040
1041   Double_t nBkg = (fitBkg1.getVal() * fracBkgRange1->getVal() + fitBkg2.getVal() * fracBkgRange2->getVal());
1042   Double_t errBkg = (fitBkg1.getError() * fracBkgRange1->getVal() + fitBkg2.getError() * fracBkgRange2->getVal());
1043   results->SetBinContent(results->GetXaxis()->FindBin("bkg"), nBkg);
1044   results->SetBinError(results->GetXaxis()->FindBin("bkg"), errBkg);
1045
1046   Double_t SB = 0.0;
1047   Double_t errSB = 0.0;
1048   if (nJPsi != 0.0 && nBkg != 0.0)
1049     {
1050       SB = nJPsi/nBkg;
1051       errSB = SB * (errJPsi/nJPsi + errBkg/nBkg);
1052     }
1053   results->SetBinContent(results->GetXaxis()->FindBin("SB"), SB);
1054   results->SetBinError(results->GetXaxis()->FindBin("SB"), errSB);
1055
1056   Int_t nDF = fitResult->floatParsFinal().getSize();
1057   Double_t chi2 = plot->chiSquare("fitFunction", "invariantMassFit", nDF);
1058   results->SetBinContent(results->GetXaxis()->FindBin("chi2"), chi2);
1059   results->SetBinError(results->GetXaxis()->FindBin("chi2"), 0.0);
1060
1061   fitList->AddAt(invariantMass, 0);
1062   fitList->AddAt(plot, 1);
1063   fitList->AddAt(parameters, 2);
1064   fitList->AddAt(results, 3);
1065 }
1066
1067
1068 void ProduceDimuonGraph(TFile *outputFile, std::vector<Double_t> multiplicityRanges)
1069 {
1070   // This function will create all the graphs for dimuons (J/Psi and Background)
1071   // - raw
1072   // - yield
1073   // - yield over mean mult
1074   // - yield over mean mult normalised
1075
1076   // First, get the CINT1B and CMUS1B graphs
1077   TGraphErrors *CINT1B = (TGraphErrors *) outputFile->Get("graphCINT1B");
1078   TGraphErrors *CMUS1B = (TGraphErrors *) outputFile->Get("graphCMUS1B");
1079
1080   TList *JPsiGraph = new TList();
1081   JPsiGraph->SetName("JPsiGraph");
1082   TList *bkgGraph = new TList();
1083   bkgGraph->SetName("bkgGraph");
1084
1085   // Raw Graphs
1086   TGraphAsymmErrors *rawJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1087   rawJPsiGraph->SetName("rawJPsiGraph");
1088   TGraphAsymmErrors *rawBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1089   rawBkgGraph->SetName("rawBkgGraph");
1090
1091   for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++)
1092     {
1093       TString name;
1094       name.Form("Multiplicity_%d-%d", static_cast<Int_t>(multiplicityRanges[iMult]), static_cast<Int_t>(multiplicityRanges[iMult+1]));
1095       TH1D *results = (TH1D *) ((TList *) ((TList *) outputFile->Get("dimuonFitResults;1"))->FindObject(name))->FindObject("results");
1096
1097       // J/Psi
1098       Double_t nJPsi = results->GetBinContent(results->GetXaxis()->FindBin("JPsi"));
1099       Double_t errJPsi = results->GetBinError(results->GetXaxis()->FindBin("JPsi"));
1100
1101       rawJPsiGraph->SetPoint(iMult, CINT1B->GetX()[iMult], nJPsi);
1102       rawJPsiGraph->SetPointError(iMult,
1103                                   CINT1B->GetX()[iMult] - multiplicityRanges[iMult],
1104                                   multiplicityRanges[iMult+1] - CINT1B->GetX()[iMult],
1105                                   errJPsi, errJPsi);
1106
1107       // Background
1108       Double_t nBkg = results->GetBinContent(results->GetXaxis()->FindBin("bkg"));
1109       Double_t errBkg = results->GetBinError(results->GetXaxis()->FindBin("bkg"));
1110
1111       rawJPsiGraph->SetPoint(iMult, CINT1B->GetX()[iMult], nBkg);
1112       rawJPsiGraph->SetPointError(iMult,
1113                                   CINT1B->GetX()[iMult] - multiplicityRanges[iMult],
1114                                   multiplicityRanges[iMult+1] - CINT1B->GetX()[iMult],
1115                                   errBkg, errBkg);
1116     }
1117
1118   JPsiGraph->AddAt(rawJPsiGraph, 0);
1119   bkgGraph->AddAt(rawBkgGraph, 0);
1120
1121   // Yield Graphs
1122   TGraphAsymmErrors *yieldJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1123   yieldJPsiGraph->SetName("yieldJPsiGraph");
1124   TGraphAsymmErrors *yieldBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1125   yieldBkgGraph->SetName("yieldBkgGraph");
1126
1127   for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++)
1128     {
1129       if (CINT1B->GetY()[iMult] != 0.0)
1130         {
1131           // J/Psi
1132           yieldJPsiGraph->SetPoint(iMult, rawJPsiGraph->GetX()[iMult], rawJPsiGraph->GetY()[iMult]/CINT1B->GetY()[iMult]);
1133           if (rawJPsiGraph->GetY()[iMult] != 0.0)
1134             {
1135               Double_t error = yieldJPsiGraph->GetY()[iMult]*
1136                 TMath::Sqrt((rawJPsiGraph->GetEYlow()[iMult]/rawJPsiGraph->GetY()[iMult])*(rawJPsiGraph->GetEYlow()[iMult]/rawJPsiGraph->GetY()[iMult]) +
1137                             (CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])*(CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult]));
1138               yieldJPsiGraph->SetPointError(iMult,
1139                                             rawJPsiGraph->GetEXlow()[iMult],
1140                                             rawJPsiGraph->GetEXhigh()[iMult],
1141                                             error, error);
1142             }
1143           // Background
1144           yieldBkgGraph->SetPoint(iMult, rawBkgGraph->GetX()[iMult], rawBkgGraph->GetY()[iMult]/CINT1B->GetY()[iMult]);
1145           if (rawBkgGraph->GetY()[iMult] != 0.0)
1146             {
1147               Double_t error = yieldBkgGraph->GetY()[iMult]*
1148                 TMath::Sqrt((rawBkgGraph->GetEYlow()[iMult]/rawBkgGraph->GetY()[iMult])*(rawBkgGraph->GetEYlow()[iMult]/rawBkgGraph->GetY()[iMult]) +
1149                             (CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult])*(CINT1B->GetEY()[iMult]/CINT1B->GetY()[iMult]));
1150
1151               yieldBkgGraph->SetPointError(iMult,
1152                                            rawBkgGraph->GetEXlow()[iMult],
1153                                            rawBkgGraph->GetEXhigh()[iMult],
1154                                            error, error);
1155             }
1156         }
1157
1158       else
1159         {
1160           yieldJPsiGraph->SetPoint(iMult, rawJPsiGraph->GetX()[iMult], 0);
1161           yieldBkgGraph->SetPoint(iMult, rawBkgGraph->GetX()[iMult], 0);
1162         }
1163     }
1164
1165   JPsiGraph->AddAt(yieldJPsiGraph, 1);
1166   bkgGraph->AddAt(yieldBkgGraph, 1);
1167
1168   // Yield over Mean Mult graph
1169   TGraphAsymmErrors *yieldOverMeanMultJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1170   yieldOverMeanMultJPsiGraph->SetName("yieldOverMeanMultJPsiGraph");
1171   TGraphAsymmErrors *yieldOverMeanMultBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1172   yieldOverMeanMultBkgGraph->SetName("yieldOverMeanMultBkgGraph");
1173
1174   for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++)
1175     {
1176       if (CINT1B->GetX()[iMult] != 0.0)
1177         {
1178           // J/Psi
1179           yieldOverMeanMultJPsiGraph->SetPoint(iMult, yieldJPsiGraph->GetX()[iMult], yieldJPsiGraph->GetY()[iMult]/CINT1B->GetX()[iMult]);
1180           if (yieldJPsiGraph->GetY()[iMult] != 0.0)
1181             {
1182               Double_t error = yieldOverMeanMultJPsiGraph->GetY()[iMult]*
1183                 TMath::Sqrt((yieldJPsiGraph->GetEYlow()[iMult]/yieldJPsiGraph->GetY()[iMult])*(yieldJPsiGraph->GetEYlow()[iMult]/yieldJPsiGraph->GetY()[iMult]) +
1184                             (CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])*(CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult]));
1185               yieldOverMeanMultJPsiGraph->SetPointError(iMult,
1186                                                         yieldJPsiGraph->GetEXlow()[iMult],
1187                                                         yieldJPsiGraph->GetEXhigh()[iMult],
1188                                                         error, error);
1189             }
1190           // Background
1191           yieldOverMeanMultBkgGraph->SetPoint(iMult, yieldBkgGraph->GetX()[iMult], yieldBkgGraph->GetY()[iMult]/CINT1B->GetX()[iMult]);
1192           if (yieldBkgGraph->GetY()[iMult] != 0.0)
1193             {
1194               Double_t error = yieldOverMeanMultBkgGraph->GetY()[iMult]*
1195                 TMath::Sqrt((yieldBkgGraph->GetEYlow()[iMult]/yieldBkgGraph->GetY()[iMult])*(yieldBkgGraph->GetEYlow()[iMult]/yieldBkgGraph->GetY()[iMult]) +
1196                             (CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult])*(CINT1B->GetEX()[iMult]/CINT1B->GetX()[iMult]));
1197               yieldOverMeanMultBkgGraph->SetPointError(iMult,
1198                                                        yieldBkgGraph->GetEXlow()[iMult],
1199                                                        yieldBkgGraph->GetEXhigh()[iMult],
1200                                                        error, error);
1201             }
1202         }
1203
1204       else
1205         {
1206           yieldOverMeanMultJPsiGraph->SetPoint(iMult, yieldJPsiGraph->GetX()[iMult], 0);
1207           yieldOverMeanMultBkgGraph->SetPoint(iMult, yieldBkgGraph->GetX()[iMult], 0);
1208         }
1209     }
1210
1211   JPsiGraph->AddAt(yieldOverMeanMultJPsiGraph, 2);
1212   bkgGraph->AddAt(yieldOverMeanMultBkgGraph, 2);
1213
1214   // Yield over Mean Mult normalised to get the bin with the highest number of CMUS1B equal to 1
1215   Double_t maxCMUS1B = 0.0;
1216   Int_t referenceBin = 0;
1217   for (Int_t iBin = 0; iBin < CMUS1B->GetN(); iBin++)
1218     if (CMUS1B->GetX()[iBin] > 10.0)
1219       if (CMUS1B->GetY()[iBin] > maxCMUS1B)
1220         {
1221           maxCMUS1B = CMUS1B->GetY()[iBin];
1222           referenceBin = iBin;
1223         }
1224
1225   TGraphAsymmErrors *yieldNormalisedJPsiGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1226   yieldNormalisedJPsiGraph->SetName("yieldNormalisedJPsiGraph");
1227   TGraphAsymmErrors *yieldNormalisedBkgGraph = new TGraphAsymmErrors(multiplicityRanges.size()-1);
1228   yieldNormalisedBkgGraph->SetName("yieldNormalisedBkgGraph");
1229
1230   for (UInt_t iMult = 0; iMult < multiplicityRanges.size()-1; iMult++)
1231     {
1232       // JPsi
1233       if (yieldOverMeanMultJPsiGraph->GetY()[referenceBin] != 0.0)
1234         {
1235           yieldNormalisedJPsiGraph->SetPoint(iMult,
1236                                              yieldOverMeanMultJPsiGraph->GetX()[iMult],
1237                                              yieldOverMeanMultJPsiGraph->GetY()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin]);
1238           yieldNormalisedJPsiGraph->SetPointError(iMult,
1239                                                   yieldOverMeanMultJPsiGraph->GetEXlow()[iMult],
1240                                                   yieldOverMeanMultJPsiGraph->GetEXhigh()[iMult],
1241                                                   yieldOverMeanMultJPsiGraph->GetEYlow()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin],
1242                                                   yieldOverMeanMultJPsiGraph->GetEYhigh()[iMult]/yieldOverMeanMultJPsiGraph->GetY()[referenceBin]);
1243         }
1244
1245       // Bkg
1246       if (yieldOverMeanMultBkgGraph->GetY()[referenceBin] != 0.0)
1247         {
1248           yieldNormalisedBkgGraph->SetPoint(iMult,
1249                                             yieldOverMeanMultBkgGraph->GetX()[iMult],
1250                                             yieldOverMeanMultBkgGraph->GetY()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin]);
1251           yieldNormalisedBkgGraph->SetPointError(iMult,
1252                                                  yieldOverMeanMultBkgGraph->GetEXlow()[iMult],
1253                                                  yieldOverMeanMultBkgGraph->GetEXhigh()[iMult],
1254                                                  yieldOverMeanMultBkgGraph->GetEYlow()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin],
1255                                                  yieldOverMeanMultBkgGraph->GetEYhigh()[iMult]/yieldOverMeanMultBkgGraph->GetY()[referenceBin]);
1256         }
1257
1258     }
1259
1260   JPsiGraph->AddAt(yieldNormalisedJPsiGraph, 3);
1261   bkgGraph->AddAt(yieldNormalisedBkgGraph, 3);
1262
1263   outputFile->WriteTObject(JPsiGraph);
1264   outputFile->WriteTObject(bkgGraph);
1265 }