| 0a65e325 |
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 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame