]>
Commit | Line | Data |
---|---|---|
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 | } |