]>
Commit | Line | Data |
---|---|---|
0ab29cfa | 1 | /* $Id$ */ |
2 | ||
3 | // | |
4 | // script to run the AliMultiplicityESDSelector | |
5 | // | |
6 | ||
7 | #include "../CreateESDChain.C" | |
0bd1f8a0 | 8 | #include "../PWG0Helper.C" |
0ab29cfa | 9 | |
03244034 | 10 | void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "") |
0ab29cfa | 11 | { |
0bd1f8a0 | 12 | if (aProof) |
6d81c2de | 13 | { |
14 | TProof::Mgr("lxb6046")->SetROOTVersion("vHEAD-07Jun2007b_dbg"); | |
0b4bfd98 | 15 | connectProof("lxb6046"); |
6d81c2de | 16 | } |
0bd1f8a0 | 17 | |
18 | TString libraries("libEG;libGeom;libESD;libPWG0base"); | |
19 | TString packages("PWG0base"); | |
10ebe68d | 20 | |
0ab29cfa | 21 | if (aMC != kFALSE) |
0bd1f8a0 | 22 | { |
23 | libraries += ";libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6"; | |
24 | packages += ";PWG0dep"; | |
25 | } | |
26 | ||
27 | if (!prepareQuery(libraries, packages, kTRUE)) | |
28 | return; | |
0ab29cfa | 29 | |
0ab29cfa | 30 | gROOT->ProcessLine(".L CreateCuts.C"); |
31 | gROOT->ProcessLine(".L drawPlots.C"); | |
32 | ||
0ab29cfa | 33 | // selection of esd tracks |
34 | AliESDtrackCuts* esdTrackCuts = CreateTrackCuts(); | |
35 | if (!esdTrackCuts) | |
36 | { | |
37 | printf("ERROR: esdTrackCuts could not be created\n"); | |
38 | return; | |
39 | } | |
40 | ||
0bd1f8a0 | 41 | TList inputList; |
42 | inputList.Add(esdTrackCuts); | |
43 | ||
6d81c2de | 44 | // pt study |
45 | if (TString(option).Contains("pt-spectrum-func")) | |
46 | { | |
47 | //TF1* func = new TF1("func", "0.7 + x", 0, 0.3); | |
48 | //TF1* func = new TF1("func", "1.3 - x", 0, 0.3); | |
49 | TF1* func = new TF1("func", "1", 0, 0.3); | |
50 | new TCanvas; func->Draw(); | |
51 | inputList.Add(func->GetHistogram()->Clone("pt-spectrum")); | |
52 | } | |
53 | ||
cfc19dd5 | 54 | TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE); |
0ab29cfa | 55 | |
56 | TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector"); | |
57 | AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo); | |
58 | ||
0bd1f8a0 | 59 | selectorName += ".cxx+"; |
0ab29cfa | 60 | |
61 | if (aDebug != kFALSE) | |
62 | selectorName += "g"; | |
63 | ||
447c325d | 64 | //Int_t result = chain->Process(selectorName, option); |
03244034 | 65 | Int_t result = executeQuery(chain, &inputList, selectorName, option); |
0ab29cfa | 66 | |
67 | if (result != 0) | |
68 | { | |
69 | printf("ERROR: Executing process failed with %d.\n", result); | |
70 | return; | |
71 | } | |
0a173978 | 72 | } |
0ab29cfa | 73 | |
0b4bfd98 | 74 | void SetTPC() |
0a173978 | 75 | { |
76 | gSystem->Load("libPWG0base"); | |
0b4bfd98 | 77 | AliMultiplicityCorrection::SetQualityRegions(kFALSE); |
78 | } | |
0a173978 | 79 | |
0b4bfd98 | 80 | void draw(const char* fileName = "multiplicityMC.root", const char* folder = "Multiplicity") |
81 | { | |
82 | gSystem->Load("libPWG0base"); | |
0a173978 | 83 | |
0b4bfd98 | 84 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder); |
0a173978 | 85 | |
0b4bfd98 | 86 | TFile::Open(fileName); |
87 | mult->LoadHistograms(); | |
0a173978 | 88 | mult->DrawHistograms(); |
447c325d | 89 | |
0b4bfd98 | 90 | return; |
91 | ||
447c325d | 92 | TH2* hist = (TH2*) gROOT->FindObject("fCorrelation3_zy"); |
93 | canvas = new TCanvas("c1", "c1", 600, 500); | |
94 | hist->SetStats(kFALSE); | |
95 | hist->Draw("COLZ"); | |
96 | hist->SetTitle(";true multiplicity in |#eta| < 2;measured multiplicity in |#eta| < 2"); | |
97 | hist->GetYaxis()->SetTitleOffset(1.1); | |
98 | gPad->SetRightMargin(0.15); | |
99 | gPad->SetLogz(); | |
100 | ||
101 | canvas->SaveAs("Plot_Correlation.pdf"); | |
0ab29cfa | 102 | } |
103 | ||
0b4bfd98 | 104 | void fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityMC_3M.root", Bool_t chi2 = kTRUE, Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE) |
9ca6be09 | 105 | { |
106 | gSystem->Load("libPWG0base"); | |
107 | ||
0b4bfd98 | 108 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder); |
9ca6be09 | 109 | |
110 | TFile::Open(fileNameMC); | |
0b4bfd98 | 111 | mult->LoadHistograms(); |
9ca6be09 | 112 | |
113 | TFile::Open(fileNameESD); | |
cfc19dd5 | 114 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); |
447c325d | 115 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID))); |
116 | //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID)); | |
9ca6be09 | 117 | |
cfc19dd5 | 118 | mult->SetMultiplicityESD(histID, hist); |
9ca6be09 | 119 | |
447c325d | 120 | // small hack to get around charge conservation for full phase space ;-) |
121 | if (fullPhaseSpace) | |
122 | { | |
123 | TH1* corr = mult->GetCorrelation(histID + 4); | |
dd701109 | 124 | |
447c325d | 125 | for (Int_t i=2; i<=corr->GetNbinsX(); i+=2) |
126 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) | |
127 | { | |
128 | corr->SetBinContent(i, j, corr->GetBinContent(i-1, j)); | |
129 | corr->SetBinError(i, j, corr->GetBinError(i-1, j)); | |
130 | } | |
131 | } | |
dd701109 | 132 | |
447c325d | 133 | /*mult->SetMultiplicityVtx(histID, hist2); |
134 | mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); | |
135 | return;*/ | |
cfc19dd5 | 136 | |
447c325d | 137 | if (chi2) |
138 | { | |
139 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4); | |
0b4bfd98 | 140 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e3); |
141 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5); | |
447c325d | 142 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); //, hist2->ProjectionY("mymchist")); |
143 | mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); | |
144 | } | |
145 | else | |
dd701109 | 146 | { |
6d81c2de | 147 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.2, 100); |
447c325d | 148 | mult->DrawComparison("Bayesian", histID, kFALSE, kTRUE, hist2->ProjectionY("mymchist2")); |
dd701109 | 149 | } |
150 | ||
151 | //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7); | |
152 | //mult->ApplyMinuitFit(histID, kFALSE); | |
153 | //mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY()); | |
154 | ||
447c325d | 155 | } |
156 | ||
157 | void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE) | |
158 | { | |
159 | gSystem->Load("libPWG0base"); | |
160 | ||
161 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
162 | ||
163 | TFile::Open(fileNameMC); | |
164 | mult->LoadHistograms("Multiplicity"); | |
165 | ||
166 | TFile::Open(fileNameESD); | |
167 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); | |
168 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID))); | |
169 | //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID)); | |
170 | ||
171 | mult->SetMultiplicityESD(histID, hist); | |
172 | ||
173 | // small hack to get around charge conservation for full phase space ;-) | |
174 | if (fullPhaseSpace) | |
175 | { | |
176 | TH1* corr = mult->GetCorrelation(histID + 4); | |
177 | ||
178 | for (Int_t i=2; i<=corr->GetNbinsX(); i+=2) | |
179 | for (Int_t j=1; j<=corr->GetNbinsY(); ++j) | |
180 | { | |
181 | corr->SetBinContent(i, j, corr->GetBinContent(i-1, j)); | |
182 | corr->SetBinError(i, j, corr->GetBinError(i-1, j)); | |
183 | } | |
184 | } | |
185 | ||
186 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); | |
187 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE); | |
188 | mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); | |
189 | ||
190 | TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult"); | |
191 | ||
192 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000); | |
193 | mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result); | |
194 | mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist")); | |
0a173978 | 195 | |
196 | return mult; | |
197 | } | |
cfc19dd5 | 198 | |
199 | const char* GetRegName(Int_t type) | |
200 | { | |
201 | switch (type) | |
202 | { | |
203 | case AliMultiplicityCorrection::kNone: return "None"; break; | |
204 | case AliMultiplicityCorrection::kPol0: return "Pol0"; break; | |
205 | case AliMultiplicityCorrection::kPol1: return "Pol1"; break; | |
206 | case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break; | |
207 | case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break; | |
0b4bfd98 | 208 | case AliMultiplicityCorrection::kLog : return "Log"; break; |
dd701109 | 209 | } |
210 | return 0; | |
211 | } | |
212 | ||
213 | const char* GetEventTypeName(Int_t type) | |
214 | { | |
215 | switch (type) | |
216 | { | |
217 | case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break; | |
218 | case AliMultiplicityCorrection::kMB: return "minimum bias"; break; | |
219 | case AliMultiplicityCorrection::kINEL: return "inelastic"; break; | |
cfc19dd5 | 220 | } |
221 | return 0; | |
222 | } | |
223 | ||
0b4bfd98 | 224 | /*void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
cfc19dd5 | 225 | { |
dd701109 | 226 | gSystem->mkdir(targetDir); |
227 | ||
228 | gSystem->Load("libPWG0base"); | |
229 | ||
230 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
231 | TFile::Open(fileNameMC); | |
232 | mult->LoadHistograms("Multiplicity"); | |
233 | ||
234 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); | |
235 | TFile::Open(fileNameESD); | |
236 | multESD->LoadHistograms("Multiplicity"); | |
237 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); | |
238 | ||
239 | TCanvas* canvas = new TCanvas("EvaluateBayesianMethod", "EvaluateBayesianMethod", 800, 600); | |
240 | TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); | |
241 | legend->SetFillColor(0); | |
242 | ||
243 | Float_t min = 1e10; | |
244 | Float_t max = 0; | |
245 | ||
246 | TGraph* first = 0; | |
247 | Int_t count = 0; // just to order the saved images... | |
248 | ||
249 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) | |
250 | { | |
251 | TGraph* fitResultsMC = new TGraph; | |
252 | fitResultsMC->SetTitle(";Weight Parameter"); | |
253 | TGraph* fitResultsRes = new TGraph; | |
254 | fitResultsRes->SetTitle(";Weight Parameter"); | |
255 | ||
256 | fitResultsMC->SetFillColor(0); | |
257 | fitResultsRes->SetFillColor(0); | |
258 | fitResultsMC->SetMarkerStyle(20+type); | |
259 | fitResultsRes->SetMarkerStyle(24+type); | |
260 | fitResultsRes->SetMarkerColor(kRed); | |
261 | fitResultsRes->SetLineColor(kRed); | |
262 | ||
263 | legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); | |
264 | legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); | |
265 | ||
447c325d | 266 | for (Float_t weight = 0; weight < 1.01; weight += 0.1) |
dd701109 | 267 | { |
268 | Float_t chi2MC = 0; | |
269 | Float_t residuals = 0; | |
270 | ||
0b4bfd98 | 271 | mult->ApplyBayesianMethod(histID, kFALSE, type, weight, 100, 0, kFALSE); |
dd701109 | 272 | mult->DrawComparison(Form("%s/Bayesian_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
273 | mult->GetComparisonResults(&chi2MC, 0, &residuals); | |
274 | ||
275 | fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC); | |
276 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); | |
277 | ||
278 | min = TMath::Min(min, TMath::Min(chi2MC, residuals)); | |
279 | max = TMath::Max(max, TMath::Max(chi2MC, residuals)); | |
280 | } | |
281 | ||
282 | fitResultsMC->Print(); | |
283 | fitResultsRes->Print(); | |
284 | ||
285 | canvas->cd(); | |
286 | fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); | |
287 | fitResultsRes->Draw("SAME CP"); | |
288 | ||
289 | if (first == 0) | |
290 | first = fitResultsMC; | |
291 | } | |
292 | ||
293 | gPad->SetLogy(); | |
294 | printf("min = %f, max = %f\n", min, max); | |
295 | if (min <= 0) | |
296 | min = 1e-5; | |
297 | first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); | |
298 | ||
299 | legend->Draw(); | |
300 | ||
301 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); | |
302 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); | |
303 | } | |
304 | ||
305 | void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2) | |
306 | { | |
307 | gSystem->mkdir(targetDir); | |
308 | ||
309 | gSystem->Load("libPWG0base"); | |
310 | ||
311 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
312 | TFile::Open(fileNameMC); | |
313 | mult->LoadHistograms("Multiplicity"); | |
314 | ||
315 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); | |
316 | TFile::Open(fileNameESD); | |
317 | multESD->LoadHistograms("Multiplicity"); | |
318 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); | |
319 | ||
320 | TCanvas* canvas = new TCanvas("EvaluateBayesianMethodIterations", "EvaluateBayesianMethodIterations", 800, 600); | |
321 | TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98); | |
322 | legend->SetFillColor(0); | |
323 | ||
324 | Float_t min = 1e10; | |
325 | Float_t max = 0; | |
326 | ||
327 | TGraph* first = 0; | |
328 | Int_t count = 0; // just to order the saved images... | |
329 | ||
330 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) | |
331 | { | |
332 | TGraph* fitResultsMC = new TGraph; | |
333 | fitResultsMC->SetTitle(";Iterations"); | |
334 | TGraph* fitResultsRes = new TGraph; | |
335 | fitResultsRes->SetTitle(";Iterations"); | |
336 | ||
337 | fitResultsMC->SetFillColor(0); | |
338 | fitResultsRes->SetFillColor(0); | |
339 | fitResultsMC->SetMarkerStyle(20+type); | |
340 | fitResultsRes->SetMarkerStyle(24+type); | |
341 | fitResultsRes->SetMarkerColor(kRed); | |
342 | fitResultsRes->SetLineColor(kRed); | |
343 | ||
344 | legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type))); | |
345 | legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type))); | |
346 | ||
347 | for (Int_t iter = 5; iter <= 50; iter += 5) | |
348 | { | |
349 | Float_t chi2MC = 0; | |
350 | Float_t residuals = 0; | |
351 | ||
352 | mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1, iter); | |
353 | mult->DrawComparison(Form("%s/BayesianIter_%02d_%d_%d", targetDir, count++, type, iter), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); | |
354 | mult->GetComparisonResults(&chi2MC, 0, &residuals); | |
355 | ||
356 | fitResultsMC->SetPoint(fitResultsMC->GetN(), iter, chi2MC); | |
357 | fitResultsRes->SetPoint(fitResultsRes->GetN(), iter, residuals); | |
358 | ||
359 | min = TMath::Min(min, TMath::Min(chi2MC, residuals)); | |
360 | max = TMath::Max(max, TMath::Max(chi2MC, residuals)); | |
361 | } | |
362 | ||
363 | fitResultsMC->Print(); | |
364 | fitResultsRes->Print(); | |
365 | ||
366 | canvas->cd(); | |
367 | fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME")); | |
368 | fitResultsRes->Draw("SAME CP"); | |
369 | ||
370 | if (first == 0) | |
371 | first = fitResultsMC; | |
372 | } | |
373 | ||
374 | gPad->SetLogy(); | |
375 | printf("min = %f, max = %f\n", min, max); | |
376 | if (min <= 0) | |
377 | min = 1e-5; | |
378 | first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5); | |
379 | ||
380 | legend->Draw(); | |
381 | ||
382 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); | |
383 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); | |
0b4bfd98 | 384 | }*/ |
dd701109 | 385 | |
447c325d | 386 | void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
387 | { | |
388 | gSystem->mkdir(targetDir); | |
389 | ||
390 | gSystem->Load("libPWG0base"); | |
391 | ||
392 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
393 | TFile::Open(fileNameMC); | |
394 | mult->LoadHistograms("Multiplicity"); | |
395 | ||
396 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); | |
397 | TFile::Open(fileNameESD); | |
398 | multESD->LoadHistograms("Multiplicity"); | |
399 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); | |
400 | ||
401 | Int_t count = 0; // just to order the saved images... | |
402 | ||
403 | TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE"); | |
404 | ||
0b4bfd98 | 405 | Int_t colors[3] = {1, 2, 4}; |
406 | Int_t markers[12] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3}; | |
407 | ||
408 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type) | |
409 | //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type) | |
447c325d | 410 | { |
411 | TString tmp; | |
412 | tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type)); | |
413 | ||
414 | TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600); | |
415 | ||
0b4bfd98 | 416 | for (Int_t i = 1; i <= 4; i++) |
447c325d | 417 | { |
0b4bfd98 | 418 | Int_t iterArray[4] = {5, 20, 50, 100}; |
419 | //Int_t iter = i * 40 - 20; | |
420 | Int_t iter = iterArray[i-1]; | |
421 | ||
422 | TGraph* fitResultsMC[3]; | |
423 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) | |
424 | { | |
425 | fitResultsMC[region] = new TGraph; | |
426 | fitResultsMC[region]->SetTitle(Form("%d iter. - reg. %d", iter, region+1)); | |
427 | fitResultsMC[region]->GetXaxis()->SetTitle("smoothing parameter #alpha"); | |
428 | fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region)); | |
429 | fitResultsMC[region]->SetName(Form("%s_MC_%d", tmp.Data(), i * AliMultiplicityCorrection::kQualityRegions + region - 2)); | |
430 | fitResultsMC[region]->SetFillColor(0); | |
431 | fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]); | |
432 | fitResultsMC[region]->SetLineColor(colors[region]); | |
433 | } | |
434 | ||
447c325d | 435 | TGraph* fitResultsRes = new TGraph; |
436 | fitResultsRes->SetTitle(Form("%d iterations", iter)); | |
437 | fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i)); | |
438 | fitResultsRes->GetXaxis()->SetTitle("smoothing parameter"); | |
439 | fitResultsRes->GetYaxis()->SetTitle("P_{2}"); | |
440 | ||
447c325d | 441 | fitResultsRes->SetFillColor(0); |
447c325d | 442 | fitResultsRes->SetMarkerStyle(19+i); |
0b4bfd98 | 443 | fitResultsRes->SetMarkerColor(1); |
444 | fitResultsRes->SetLineColor(1); | |
447c325d | 445 | |
446 | for (Float_t weight = 0.0; weight < 1.01; weight += 0.2) | |
447 | { | |
448 | Float_t chi2MC = 0; | |
449 | Float_t residuals = 0; | |
450 | ||
0b4bfd98 | 451 | mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter, 0, kFALSE); |
447c325d | 452 | mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY()); |
453 | mult->GetComparisonResults(&chi2MC, 0, &residuals); | |
454 | ||
0b4bfd98 | 455 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
456 | fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region)); | |
457 | ||
447c325d | 458 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
459 | } | |
460 | ||
0b4bfd98 | 461 | graphFile->cd(); |
462 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) | |
463 | fitResultsMC[region]->Write(); | |
464 | ||
447c325d | 465 | fitResultsRes->Write(); |
466 | } | |
467 | } | |
0b4bfd98 | 468 | |
469 | graphFile->Close(); | |
447c325d | 470 | } |
471 | ||
472 | void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE) | |
473 | { | |
474 | gSystem->Load("libPWG0base"); | |
475 | ||
476 | TString name; | |
477 | TFile* graphFile = 0; | |
478 | if (type == -1) | |
479 | { | |
480 | name = "EvaluateChi2Method"; | |
481 | graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir)); | |
482 | } | |
483 | else | |
484 | { | |
485 | name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type)); | |
486 | graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir)); | |
487 | } | |
488 | ||
489 | TCanvas* canvas = new TCanvas(name, name, 800, 500); | |
490 | if (type == -1) | |
0b4bfd98 | 491 | { |
447c325d | 492 | canvas->SetLogx(); |
0b4bfd98 | 493 | canvas->SetLogy(); |
494 | } | |
495 | canvas->SetTopMargin(0.05); | |
496 | canvas->SetGridx(); | |
497 | canvas->SetGridy(); | |
447c325d | 498 | |
0b4bfd98 | 499 | TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98); |
447c325d | 500 | legend->SetFillColor(0); |
501 | ||
502 | Int_t count = 1; | |
503 | ||
504 | Float_t xMin = 1e20; | |
505 | Float_t xMax = 0; | |
506 | ||
507 | Float_t yMin = 1e20; | |
508 | Float_t yMax = 0; | |
509 | ||
0b4bfd98 | 510 | Float_t yMinRegion[3]; |
511 | for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i) | |
512 | yMinRegion[i] = 1e20; | |
513 | ||
447c325d | 514 | TString xaxis, yaxis; |
515 | ||
516 | while (1) | |
517 | { | |
518 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); | |
519 | TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count)); | |
520 | ||
0b4bfd98 | 521 | if (!mc) |
447c325d | 522 | break; |
523 | ||
524 | xaxis = mc->GetXaxis()->GetTitle(); | |
525 | yaxis = mc->GetYaxis()->GetTitle(); | |
526 | ||
527 | mc->Print(); | |
447c325d | 528 | |
0b4bfd98 | 529 | if (res) |
530 | res->Print(); | |
447c325d | 531 | |
0b4bfd98 | 532 | xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin()); |
533 | yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin()); | |
447c325d | 534 | |
0b4bfd98 | 535 | xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax()); |
536 | yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax()); | |
537 | ||
538 | if (plotRes && res) | |
447c325d | 539 | { |
0b4bfd98 | 540 | xMin = TMath::Min(xMin, res->GetXaxis()->GetXmin()); |
541 | yMin = TMath::Min(yMin, res->GetYaxis()->GetXmin()); | |
447c325d | 542 | |
0b4bfd98 | 543 | xMax = TMath::Max(xMax, res->GetXaxis()->GetXmax()); |
544 | yMax = TMath::Max(yMax, res->GetYaxis()->GetXmax()); | |
447c325d | 545 | } |
0b4bfd98 | 546 | |
547 | for (Int_t i=0; i<mc->GetN(); ++i) | |
548 | yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]); | |
549 | ||
550 | count++; | |
447c325d | 551 | } |
552 | ||
0b4bfd98 | 553 | for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i) |
554 | Printf("Minimum for region %d is %f", i, yMinRegion[i]); | |
555 | ||
447c325d | 556 | if (type >= 0) |
557 | { | |
558 | xaxis = "smoothing parameter"; | |
559 | } | |
560 | else if (type == -1) | |
561 | { | |
562 | xaxis = "weight parameter"; | |
0b4bfd98 | 563 | xMax *= 5; |
447c325d | 564 | } |
0b4bfd98 | 565 | //yaxis = "P_{1} (2 <= t <= 150)"; |
447c325d | 566 | |
567 | printf("%f %f %f %f\n", xMin, xMax, yMin, yMax); | |
568 | ||
569 | TGraph* dummy = new TGraph; | |
570 | dummy->SetPoint(0, xMin, yMin); | |
571 | dummy->SetPoint(1, xMax, yMax); | |
572 | dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data())); | |
573 | ||
574 | dummy->SetMarkerColor(0); | |
575 | dummy->Draw("AP"); | |
0b4bfd98 | 576 | dummy->GetYaxis()->SetMoreLogLabels(1); |
447c325d | 577 | |
578 | count = 1; | |
579 | ||
580 | while (1) | |
581 | { | |
582 | TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count)); | |
583 | TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count)); | |
584 | ||
0b4bfd98 | 585 | //printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res); |
447c325d | 586 | |
0b4bfd98 | 587 | if (!mc) |
447c325d | 588 | break; |
589 | ||
590 | printf("Loaded %d sucessful.\n", count); | |
591 | ||
0b4bfd98 | 592 | /*if (type == -1) |
447c325d | 593 | { |
0b4bfd98 | 594 | legend->AddEntry(mc, Form("Eq. (%d)", (count-1)/3+11)); |
447c325d | 595 | } |
0b4bfd98 | 596 | else*/ |
447c325d | 597 | legend->AddEntry(mc); |
598 | ||
599 | mc->Draw("SAME PC"); | |
600 | ||
0b4bfd98 | 601 | if (plotRes && res) |
447c325d | 602 | { |
603 | legend->AddEntry(res); | |
604 | res->Draw("SAME PC"); | |
605 | } | |
606 | ||
607 | count++; | |
608 | } | |
609 | ||
610 | legend->Draw(); | |
611 | ||
612 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); | |
613 | canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName())); | |
614 | } | |
615 | ||
0b4bfd98 | 616 | void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", const char* targetDir, Int_t histID = 3) |
dd701109 | 617 | { |
618 | gSystem->mkdir(targetDir); | |
619 | ||
cfc19dd5 | 620 | gSystem->Load("libPWG0base"); |
621 | ||
622 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
623 | ||
624 | TFile::Open(fileNameMC); | |
625 | mult->LoadHistograms("Multiplicity"); | |
626 | ||
627 | TFile::Open(fileNameESD); | |
628 | TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID)); | |
629 | TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID)); | |
630 | ||
631 | mult->SetMultiplicityESD(histID, hist); | |
632 | ||
dd701109 | 633 | Int_t count = 0; // just to order the saved images... |
0b4bfd98 | 634 | Int_t colors[3] = {1, 2, 4}; |
635 | Int_t markers[12] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3}; | |
dd701109 | 636 | |
447c325d | 637 | TGraph* fitResultsRes = 0; |
638 | ||
639 | TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE"); | |
cfc19dd5 | 640 | |
0b4bfd98 | 641 | for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kLog; ++type) |
642 | // for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kEntropy; type <= AliMultiplicityCorrection::kEntropy; ++type) | |
643 | // for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kPol0; ++type) | |
cfc19dd5 | 644 | { |
0b4bfd98 | 645 | TGraph* fitResultsMC[3]; |
646 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) | |
647 | { | |
648 | fitResultsMC[region] = new TGraph; | |
649 | fitResultsMC[region]->SetTitle(Form("Eq. (%d) - reg. %d", type+10, region+1)); | |
650 | fitResultsMC[region]->GetXaxis()->SetTitle("weight parameter #alpha"); | |
651 | fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region)); | |
652 | fitResultsMC[region]->SetName(Form("EvaluateChi2Method_MC_%d", type * AliMultiplicityCorrection::kQualityRegions + region - 2)); | |
653 | fitResultsMC[region]->SetFillColor(0); | |
654 | fitResultsMC[region]->SetMarkerStyle(markers[(type-1) * AliMultiplicityCorrection::kQualityRegions + region]); | |
655 | fitResultsMC[region]->SetLineColor(colors[region]); | |
656 | } | |
657 | ||
447c325d | 658 | fitResultsRes = new TGraph; |
659 | fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type))); | |
660 | fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type)); | |
661 | fitResultsRes->GetXaxis()->SetTitle("Weight Parameter"); | |
cfc19dd5 | 662 | |
cfc19dd5 | 663 | fitResultsRes->SetFillColor(0); |
dd701109 | 664 | fitResultsRes->SetMarkerStyle(23+type); |
cfc19dd5 | 665 | fitResultsRes->SetMarkerColor(kRed); |
666 | fitResultsRes->SetLineColor(kRed); | |
667 | ||
0b4bfd98 | 668 | for (Int_t i=0; i<7; ++i) |
447c325d | 669 | { |
670 | Float_t weight = TMath::Power(TMath::Sqrt(10), i+6); | |
671 | //Float_t weight = TMath::Power(10, i+2); | |
cfc19dd5 | 672 | |
447c325d | 673 | //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5; |
cfc19dd5 | 674 | |
cfc19dd5 | 675 | Float_t chi2MC = 0; |
676 | Float_t residuals = 0; | |
447c325d | 677 | Float_t chi2Limit = 0; |
678 | ||
679 | TString runName; | |
680 | runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight); | |
cfc19dd5 | 681 | |
682 | mult->SetRegularizationParameters(type, weight); | |
447c325d | 683 | Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
684 | mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY()); | |
685 | if (status != 0) | |
686 | { | |
687 | printf("MINUIT did not succeed. Skipping...\n"); | |
688 | continue; | |
689 | } | |
690 | ||
dd701109 | 691 | mult->GetComparisonResults(&chi2MC, 0, &residuals); |
447c325d | 692 | TH1* result = mult->GetMultiplicityESDCorrected(histID); |
693 | result->SetName(runName); | |
694 | result->Write(); | |
cfc19dd5 | 695 | |
0b4bfd98 | 696 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
697 | fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region)); | |
698 | ||
cfc19dd5 | 699 | fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals); |
cfc19dd5 | 700 | } |
701 | ||
447c325d | 702 | graphFile->cd(); |
0b4bfd98 | 703 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
704 | fitResultsMC[region]->Write(); | |
447c325d | 705 | fitResultsRes->Write(); |
cfc19dd5 | 706 | } |
707 | ||
447c325d | 708 | graphFile->Close(); |
dd701109 | 709 | } |
710 | ||
711 | void EvaluateChi2MethodAll() | |
712 | { | |
713 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); | |
714 | EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); | |
715 | EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); | |
716 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); | |
717 | EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); | |
718 | } | |
719 | ||
720 | void EvaluateBayesianMethodAll() | |
721 | { | |
722 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M"); | |
723 | EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M"); | |
724 | EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD"); | |
725 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M"); | |
726 | EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD"); | |
727 | } | |
728 | ||
447c325d | 729 | void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3) |
dd701109 | 730 | { |
731 | gSystem->mkdir(targetDir); | |
732 | ||
733 | gSystem->Load("libPWG0base"); | |
734 | ||
735 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
736 | ||
737 | TFile::Open(fileNameMC); | |
738 | mult->LoadHistograms("Multiplicity"); | |
739 | ||
740 | TFile::Open(fileNameESD); | |
741 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); | |
742 | multESD->LoadHistograms("Multiplicity"); | |
743 | ||
744 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); | |
745 | ||
447c325d | 746 | TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200); |
747 | canvas->Divide(3, 3); | |
dd701109 | 748 | |
749 | Int_t count = 0; | |
750 | ||
447c325d | 751 | for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type) |
dd701109 | 752 | { |
447c325d | 753 | TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc"); |
754 | mc->Sumw2(); | |
dd701109 | 755 | |
447c325d | 756 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 | 757 | mult->ApplyMinuitFit(histID, kFALSE, type); |
758 | mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc); | |
759 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result"); | |
760 | ||
761 | mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1); | |
762 | mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc); | |
763 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult"); | |
764 | ||
765 | mc->GetXaxis()->SetRangeUser(0, 150); | |
766 | chi2Result->GetXaxis()->SetRangeUser(0, 150); | |
767 | ||
447c325d | 768 | /* // skip errors for now |
dd701109 | 769 | for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i) |
770 | { | |
771 | chi2Result->SetBinError(i, 0); | |
772 | bayesResult->SetBinError(i, 0); | |
447c325d | 773 | }*/ |
dd701109 | 774 | |
775 | canvas->cd(++count); | |
776 | mc->SetFillColor(kYellow); | |
777 | mc->DrawCopy(); | |
778 | chi2Result->SetLineColor(kRed); | |
779 | chi2Result->DrawCopy("SAME"); | |
780 | bayesResult->SetLineColor(kBlue); | |
781 | bayesResult->DrawCopy("SAME"); | |
782 | gPad->SetLogy(); | |
783 | ||
784 | canvas->cd(count + 3); | |
447c325d | 785 | chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio"); |
786 | bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio"); | |
787 | chi2ResultRatio->Divide(chi2Result, mc, 1, 1, ""); | |
788 | bayesResultRatio->Divide(bayesResult, mc, 1, 1, ""); | |
dd701109 | 789 | |
447c325d | 790 | chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5); |
dd701109 | 791 | |
447c325d | 792 | chi2ResultRatio->DrawCopy("HIST"); |
793 | bayesResultRatio->DrawCopy("SAME HIST"); | |
dd701109 | 794 | |
447c325d | 795 | canvas->cd(count + 6); |
796 | chi2Result->Divide(chi2Result, bayesResult, 1, 1, ""); | |
797 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); | |
798 | chi2Result->DrawCopy("HIST"); | |
dd701109 | 799 | } |
800 | ||
801 | canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName())); | |
802 | canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName())); | |
803 | } | |
804 | ||
447c325d | 805 | void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
dd701109 | 806 | { |
807 | gSystem->Load("libPWG0base"); | |
808 | ||
809 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
810 | ||
811 | TFile::Open(fileNameMC); | |
812 | mult->LoadHistograms("Multiplicity"); | |
813 | ||
447c325d | 814 | const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" }; |
dd701109 | 815 | |
0b4bfd98 | 816 | TGraph* fitResultsChi2[3]; |
817 | TGraph* fitResultsBayes[3]; | |
818 | ||
819 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) | |
820 | { | |
821 | fitResultsChi2[region] = new TGraph; | |
822 | fitResultsChi2[region]->SetTitle(";Nevents;Chi2"); | |
823 | fitResultsChi2[region]->SetName(Form("fitResultsChi2_%d", region)); | |
824 | fitResultsChi2[region]->SetMarkerStyle(20+region); | |
825 | ||
826 | fitResultsBayes[region] = new TGraph; | |
827 | fitResultsBayes[region]->SetTitle(";Nevents;Chi2"); | |
828 | fitResultsBayes[region]->SetName(Form("fitResultsBayes_%d", region)); | |
829 | fitResultsBayes[region]->SetMarkerStyle(20+region); | |
830 | fitResultsBayes[region]->SetMarkerColor(2); | |
831 | } | |
832 | ||
447c325d | 833 | TGraph* fitResultsChi2Limit = new TGraph; fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach"); |
834 | TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach"); | |
0b4bfd98 | 835 | TGraph* fitResultsChi2Res = new TGraph; fitResultsChi2Res->SetTitle(";Nevents;Chi2"); |
836 | TGraph* fitResultsBayesRes = new TGraph; fitResultsBayesRes->SetTitle(";Nevents;Chi2"); | |
447c325d | 837 | |
447c325d | 838 | fitResultsChi2Limit->SetName("fitResultsChi2Limit"); |
839 | fitResultsBayesLimit->SetName("fitResultsBayesLimit"); | |
0b4bfd98 | 840 | fitResultsChi2Res->SetName("fitResultsChi2Res"); |
841 | fitResultsBayesRes->SetName("fitResultsBayesRes"); | |
dd701109 | 842 | |
843 | TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600); | |
844 | canvas->Divide(5, 2); | |
845 | ||
846 | Float_t min = 1e10; | |
847 | Float_t max = 0; | |
848 | ||
447c325d | 849 | TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE"); |
850 | file->Close(); | |
851 | ||
dd701109 | 852 | for (Int_t i=0; i<5; ++i) |
853 | { | |
854 | TFile::Open(files[i]); | |
855 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); | |
856 | multESD->LoadHistograms("Multiplicity"); | |
857 | ||
858 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); | |
859 | Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries(); | |
447c325d | 860 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i)); |
dd701109 | 861 | |
447c325d | 862 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 | 863 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx); |
864 | mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); | |
865 | ||
dd701109 | 866 | Int_t chi2MCLimit = 0; |
0b4bfd98 | 867 | Float_t chi2Residuals = 0; |
868 | mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals); | |
869 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) | |
870 | { | |
871 | fitResultsChi2[region]->SetPoint(fitResultsChi2[region]->GetN(), nEntries, mult->GetQuality(region)); | |
872 | min = TMath::Min(min, mult->GetQuality(region)); | |
873 | max = TMath::Max(max, mult->GetQuality(region)); | |
874 | } | |
dd701109 | 875 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit); |
0b4bfd98 | 876 | fitResultsChi2Res->SetPoint(fitResultsChi2Res->GetN(), nEntries, chi2Residuals); |
dd701109 | 877 | |
447c325d | 878 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
dd701109 | 879 | |
0b4bfd98 | 880 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0, kFALSE); |
dd701109 | 881 | mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
447c325d | 882 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
0b4bfd98 | 883 | mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals); |
884 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) | |
885 | { | |
886 | fitResultsBayes[region]->SetPoint(fitResultsBayes[region]->GetN(), nEntries, mult->GetQuality(region)); | |
887 | min = TMath::Min(min, mult->GetQuality(region)); | |
888 | max = TMath::Max(max, mult->GetQuality(region)); | |
889 | } | |
dd701109 | 890 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit); |
0b4bfd98 | 891 | fitResultsBayesRes->SetPoint(fitResultsBayesRes->GetN(), nEntries, chi2Residuals); |
dd701109 | 892 | |
dd701109 | 893 | mc->GetXaxis()->SetRangeUser(0, 150); |
894 | chi2Result->GetXaxis()->SetRangeUser(0, 150); | |
895 | ||
896 | // skip errors for now | |
897 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) | |
898 | { | |
899 | chi2Result->SetBinError(j, 0); | |
900 | bayesResult->SetBinError(j, 0); | |
901 | } | |
902 | ||
903 | canvas->cd(i+1); | |
904 | mc->SetFillColor(kYellow); | |
905 | mc->DrawCopy(); | |
906 | chi2Result->SetLineColor(kRed); | |
907 | chi2Result->DrawCopy("SAME"); | |
908 | bayesResult->SetLineColor(kBlue); | |
909 | bayesResult->DrawCopy("SAME"); | |
910 | gPad->SetLogy(); | |
911 | ||
912 | canvas->cd(i+6); | |
913 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); | |
914 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); | |
915 | ||
916 | // skip errors for now | |
917 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) | |
918 | { | |
919 | chi2Result->SetBinError(j, 0); | |
920 | bayesResult->SetBinError(j, 0); | |
921 | } | |
922 | ||
923 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); | |
924 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); | |
925 | ||
926 | chi2Result->DrawCopy(""); | |
927 | bayesResult->DrawCopy("SAME"); | |
447c325d | 928 | |
929 | TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE"); | |
930 | mc->Write(); | |
931 | chi2Result->Write(); | |
932 | bayesResult->Write(); | |
933 | file->Close(); | |
dd701109 | 934 | } |
935 | ||
936 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); | |
937 | canvas->SaveAs(Form("%s.C", canvas->GetName())); | |
938 | ||
939 | TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400); | |
940 | canvas2->Divide(2, 1); | |
941 | ||
942 | canvas2->cd(1); | |
dd701109 | 943 | |
0b4bfd98 | 944 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) |
945 | { | |
946 | fitResultsChi2[region]->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); | |
947 | fitResultsChi2[region]->Draw(((region == 0) ? "AP" : "P SAME")); | |
948 | ||
949 | fitResultsBayes[region]->Draw("P SAME"); | |
950 | } | |
dd701109 | 951 | |
952 | gPad->SetLogy(); | |
953 | ||
954 | canvas2->cd(2); | |
955 | fitResultsChi2Limit->SetMarkerStyle(20); | |
956 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); | |
957 | fitResultsChi2Limit->Draw("AP"); | |
958 | ||
959 | fitResultsBayesLimit->SetMarkerStyle(3); | |
960 | fitResultsBayesLimit->SetMarkerColor(2); | |
961 | fitResultsBayesLimit->Draw("P SAME"); | |
962 | ||
963 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); | |
964 | canvas2->SaveAs(Form("%s.C", canvas2->GetName())); | |
447c325d | 965 | |
966 | TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE"); | |
0b4bfd98 | 967 | |
968 | for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region) | |
969 | { | |
970 | fitResultsChi2[region]->Write(); | |
971 | fitResultsBayes[region]->Write(); | |
972 | } | |
447c325d | 973 | fitResultsChi2Limit->Write(); |
974 | fitResultsBayesLimit->Write(); | |
0b4bfd98 | 975 | fitResultsChi2Res->Write(); |
976 | fitResultsBayesRes->Write(); | |
447c325d | 977 | file->Close(); |
dd701109 | 978 | } |
979 | ||
447c325d | 980 | void StartingConditions(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
dd701109 | 981 | { |
982 | gSystem->Load("libPWG0base"); | |
983 | ||
984 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
985 | ||
986 | TFile::Open(fileNameMC); | |
987 | mult->LoadHistograms("Multiplicity"); | |
988 | ||
447c325d | 989 | const char* files[] = { "multiplicityMC_1M_3.root", "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root" } |
dd701109 | 990 | |
991 | // this one we try to unfold | |
992 | TFile::Open(files[0]); | |
993 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD"); | |
994 | multESD->LoadHistograms("Multiplicity"); | |
995 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); | |
447c325d | 996 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc"); |
dd701109 | 997 | |
998 | TGraph* fitResultsChi2 = new TGraph; | |
999 | fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); | |
1000 | TGraph* fitResultsBayes = new TGraph; | |
1001 | fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); | |
1002 | TGraph* fitResultsChi2Limit = new TGraph; | |
1003 | fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); | |
1004 | TGraph* fitResultsBayesLimit = new TGraph; | |
1005 | fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); | |
1006 | ||
1007 | TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600); | |
1008 | canvas->Divide(8, 2); | |
1009 | ||
1010 | TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400); | |
1011 | canvas3->Divide(2, 1); | |
1012 | ||
1013 | Float_t min = 1e10; | |
1014 | Float_t max = 0; | |
1015 | ||
1016 | TH1* firstChi = 0; | |
1017 | TH1* firstBayesian = 0; | |
1018 | TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond"); | |
1019 | ||
1020 | TLegend* legend = new TLegend(0.7, 0.7, 1, 1); | |
1021 | ||
447c325d | 1022 | TFile* file = TFile::Open("StartingConditions.root", "RECREATE"); |
1023 | mc->Write(); | |
1024 | file->Close(); | |
1025 | ||
dd701109 | 1026 | for (Int_t i=0; i<8; ++i) |
1027 | { | |
1028 | if (i == 0) | |
1029 | { | |
1030 | startCond = (TH1*) mc->Clone("startCond2"); | |
1031 | } | |
1032 | else if (i < 6) | |
1033 | { | |
1034 | TFile::Open(files[i-1]); | |
1035 | AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); | |
1036 | multESD2->LoadHistograms("Multiplicity"); | |
1037 | startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond"); | |
1038 | } | |
1039 | else if (i == 6) | |
1040 | { | |
1041 | func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50); | |
1042 | func->SetParNames("scaling", "averagen", "k"); | |
1043 | func->SetParLimits(0, 1e-3, 1e10); | |
1044 | func->SetParLimits(1, 0.001, 1000); | |
1045 | func->SetParLimits(2, 0.001, 1000); | |
1046 | ||
1047 | func->SetParameters(1, 10, 2); | |
1048 | for (Int_t j=2; j<=startCond->GetNbinsX(); j++) | |
1049 | startCond->SetBinContent(j, func->Eval(j-1)); | |
1050 | } | |
1051 | else if (i == 7) | |
1052 | { | |
1053 | for (Int_t j=1; j<=startCond->GetNbinsX(); j++) | |
1054 | startCond->SetBinContent(j, 1); | |
1055 | } | |
1056 | ||
447c325d | 1057 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); |
dd701109 | 1058 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond); |
1059 | mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); | |
1060 | ||
1061 | Float_t chi2MC = 0; | |
1062 | Int_t chi2MCLimit = 0; | |
1063 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); | |
1064 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); | |
1065 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); | |
1066 | min = TMath::Min(min, chi2MC); | |
1067 | max = TMath::Max(max, chi2MC); | |
1068 | ||
447c325d | 1069 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); |
dd701109 | 1070 | if (!firstChi) |
1071 | firstChi = (TH1*) chi2Result->Clone("firstChi"); | |
1072 | ||
447c325d | 1073 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, startCond); |
dd701109 | 1074 | mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); |
447c325d | 1075 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); |
dd701109 | 1076 | if (!firstBayesian) |
1077 | firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); | |
1078 | ||
1079 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); | |
1080 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); | |
1081 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); | |
1082 | ||
447c325d | 1083 | TFile* file = TFile::Open("StartingConditions.root", "UPDATE"); |
1084 | chi2Result->Write(); | |
1085 | bayesResult->Write(); | |
1086 | file->Close(); | |
1087 | ||
dd701109 | 1088 | min = TMath::Min(min, chi2MC); |
1089 | max = TMath::Max(max, chi2MC); | |
1090 | mc->GetXaxis()->SetRangeUser(0, 150); | |
1091 | chi2Result->GetXaxis()->SetRangeUser(0, 150); | |
1092 | ||
1093 | // skip errors for now | |
1094 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) | |
1095 | { | |
1096 | chi2Result->SetBinError(j, 0); | |
1097 | bayesResult->SetBinError(j, 0); | |
1098 | } | |
1099 | ||
1100 | canvas3->cd(1); | |
1101 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); | |
1102 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); | |
1103 | tmp->Divide(firstChi); | |
1104 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1105 | tmp->GetXaxis()->SetRangeUser(0, 200); | |
1106 | tmp->SetLineColor(i+1); | |
1107 | legend->AddEntry(tmp, Form("%d", i)); | |
1108 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); | |
1109 | ||
1110 | canvas3->cd(2); | |
1111 | tmp = (TH1*) bayesResult->Clone("tmp"); | |
1112 | tmp->SetTitle("Difference to best initial conditions;Npart;Ratio"); | |
1113 | tmp->Divide(firstBayesian); | |
1114 | tmp->SetLineColor(i+1); | |
1115 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1116 | tmp->GetXaxis()->SetRangeUser(0, 200); | |
1117 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); | |
1118 | ||
1119 | canvas->cd(i+1); | |
1120 | mc->SetFillColor(kYellow); | |
1121 | mc->DrawCopy(); | |
1122 | chi2Result->SetLineColor(kRed); | |
1123 | chi2Result->DrawCopy("SAME"); | |
1124 | bayesResult->SetLineColor(kBlue); | |
1125 | bayesResult->DrawCopy("SAME"); | |
1126 | gPad->SetLogy(); | |
1127 | ||
1128 | canvas->cd(i+9); | |
1129 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); | |
1130 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); | |
1131 | ||
1132 | // skip errors for now | |
1133 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) | |
1134 | { | |
1135 | chi2Result->SetBinError(j, 0); | |
1136 | bayesResult->SetBinError(j, 0); | |
1137 | } | |
1138 | ||
1139 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); | |
1140 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1141 | ||
1142 | chi2Result->DrawCopy(""); | |
1143 | bayesResult->DrawCopy("SAME"); | |
1144 | } | |
1145 | ||
1146 | canvas3->cd(1); | |
1147 | legend->Draw(); | |
1148 | ||
cfc19dd5 | 1149 | canvas->SaveAs(Form("%s.gif", canvas->GetName())); |
dd701109 | 1150 | |
1151 | TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400); | |
1152 | canvas2->Divide(2, 1); | |
1153 | ||
1154 | canvas2->cd(1); | |
1155 | fitResultsChi2->SetMarkerStyle(20); | |
1156 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); | |
1157 | fitResultsChi2->Draw("AP"); | |
1158 | ||
1159 | fitResultsBayes->SetMarkerStyle(3); | |
1160 | fitResultsBayes->SetMarkerColor(2); | |
1161 | fitResultsBayes->Draw("P SAME"); | |
1162 | ||
1163 | gPad->SetLogy(); | |
1164 | ||
1165 | canvas2->cd(2); | |
1166 | fitResultsChi2Limit->SetMarkerStyle(20); | |
1167 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); | |
1168 | fitResultsChi2Limit->Draw("AP"); | |
1169 | ||
1170 | fitResultsBayesLimit->SetMarkerStyle(3); | |
1171 | fitResultsBayesLimit->SetMarkerColor(2); | |
1172 | fitResultsBayesLimit->Draw("P SAME"); | |
1173 | ||
1174 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); | |
1175 | canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); | |
cfc19dd5 | 1176 | } |
1177 | ||
447c325d | 1178 | void DifferentSamples(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3) |
1179 | { | |
1180 | gSystem->Load("libPWG0base"); | |
1181 | ||
1182 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
1183 | ||
1184 | TFile::Open(fileNameMC); | |
1185 | mult->LoadHistograms("Multiplicity"); | |
1186 | ||
1187 | const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" }; | |
1188 | ||
1189 | TGraph* fitResultsChi2 = new TGraph; | |
1190 | fitResultsChi2->SetTitle(";Input Dist ID;Chi2"); | |
1191 | TGraph* fitResultsBayes = new TGraph; | |
1192 | fitResultsBayes->SetTitle(";Input Dist ID;Chi2"); | |
1193 | TGraph* fitResultsChi2Limit = new TGraph; | |
1194 | fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach"); | |
1195 | TGraph* fitResultsBayesLimit = new TGraph; | |
1196 | fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach"); | |
1197 | ||
1198 | TCanvas* canvasA = new TCanvas("DifferentSamplesA", "DifferentSamplesA", 1200, 600); | |
1199 | canvasA->Divide(4, 2); | |
1200 | ||
1201 | TCanvas* canvasB = new TCanvas("DifferentSamplesB", "DifferentSamplesB", 1200, 600); | |
1202 | canvasB->Divide(4, 2); | |
1203 | ||
1204 | TCanvas* canvas4 = new TCanvas("DifferentSamples4", "DifferentSamples4", 1000, 400); | |
1205 | canvas4->Divide(2, 1); | |
1206 | ||
1207 | TCanvas* canvas3 = new TCanvas("DifferentSamples3", "DifferentSamples3", 1000, 400); | |
1208 | canvas3->Divide(2, 1); | |
1209 | ||
1210 | Float_t min = 1e10; | |
1211 | Float_t max = 0; | |
1212 | ||
1213 | TH1* firstChi = 0; | |
1214 | TH1* firstBayesian = 0; | |
1215 | ||
1216 | TLegend* legend = new TLegend(0.7, 0.7, 1, 1); | |
1217 | ||
1218 | TFile* file = TFile::Open("DifferentSamples.root", "RECREATE"); | |
1219 | file->Close(); | |
1220 | ||
1221 | for (Int_t i=0; i<8; ++i) | |
1222 | { | |
1223 | TFile::Open(files[i]); | |
1224 | AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2"); | |
1225 | multESD->LoadHistograms("Multiplicity"); | |
1226 | mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID)); | |
1227 | TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i)); | |
1228 | mc->Sumw2(); | |
1229 | ||
1230 | mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000); | |
1231 | mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE); | |
1232 | mult->DrawComparison(Form("DifferentSamples_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc); | |
1233 | ||
1234 | Float_t chi2MC = 0; | |
1235 | Int_t chi2MCLimit = 0; | |
1236 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); | |
1237 | fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC); | |
1238 | fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit); | |
1239 | min = TMath::Min(min, chi2MC); | |
1240 | max = TMath::Max(max, chi2MC); | |
1241 | ||
1242 | TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i)); | |
1243 | if (!firstChi) | |
1244 | firstChi = (TH1*) chi2Result->Clone("firstChi"); | |
1245 | ||
1246 | mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100); | |
1247 | mult->DrawComparison(Form("DifferentSamples_%d_Bayesian", i), histID, kFALSE, kTRUE, mc); | |
1248 | TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i)); | |
1249 | if (!firstBayesian) | |
1250 | firstBayesian = (TH1*) bayesResult->Clone("firstBayesian"); | |
1251 | ||
1252 | TFile* file = TFile::Open("DifferentSamples.root", "UPDATE"); | |
1253 | mc->Write(); | |
1254 | chi2Result->Write(); | |
1255 | bayesResult->Write(); | |
1256 | file->Close(); | |
1257 | ||
1258 | mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0); | |
1259 | fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC); | |
1260 | fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit); | |
1261 | ||
1262 | min = TMath::Min(min, chi2MC); | |
1263 | max = TMath::Max(max, chi2MC); | |
1264 | mc->GetXaxis()->SetRangeUser(0, 150); | |
1265 | chi2Result->GetXaxis()->SetRangeUser(0, 150); | |
1266 | ||
1267 | // skip errors for now | |
1268 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) | |
1269 | { | |
1270 | chi2Result->SetBinError(j, 0); | |
1271 | bayesResult->SetBinError(j, 0); | |
1272 | } | |
1273 | ||
1274 | canvas4->cd(1); | |
1275 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); | |
1276 | tmp->SetTitle("Unfolded/MC;Npart;Ratio"); | |
1277 | tmp->Divide(mc); | |
1278 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1279 | tmp->GetXaxis()->SetRangeUser(0, 200); | |
1280 | tmp->SetLineColor(i+1); | |
1281 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); | |
1282 | ||
1283 | canvas4->cd(2); | |
1284 | tmp = (TH1*) bayesResult->Clone("tmp"); | |
1285 | tmp->SetTitle("Unfolded/MC;Npart;Ratio"); | |
1286 | tmp->Divide(mc); | |
1287 | tmp->SetLineColor(i+1); | |
1288 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1289 | tmp->GetXaxis()->SetRangeUser(0, 200); | |
1290 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); | |
1291 | ||
1292 | canvas3->cd(1); | |
1293 | TH1* tmp = (TH1*) chi2Result->Clone("tmp"); | |
1294 | tmp->SetTitle("Ratio to first result;Npart;Ratio"); | |
1295 | tmp->Divide(firstChi); | |
1296 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1297 | tmp->GetXaxis()->SetRangeUser(0, 200); | |
1298 | tmp->SetLineColor(i+1); | |
1299 | legend->AddEntry(tmp, Form("%d", i)); | |
1300 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); | |
1301 | ||
1302 | canvas3->cd(2); | |
1303 | tmp = (TH1*) bayesResult->Clone("tmp"); | |
1304 | tmp->SetTitle("Ratio to first result;Npart;Ratio"); | |
1305 | tmp->Divide(firstBayesian); | |
1306 | tmp->SetLineColor(i+1); | |
1307 | tmp->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1308 | tmp->GetXaxis()->SetRangeUser(0, 200); | |
1309 | tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST"); | |
1310 | ||
1311 | if (i < 4) | |
1312 | { | |
1313 | canvasA->cd(i+1); | |
1314 | } | |
1315 | else | |
1316 | canvasB->cd(i+1-4); | |
1317 | ||
1318 | mc->SetFillColor(kYellow); | |
1319 | mc->DrawCopy(); | |
1320 | chi2Result->SetLineColor(kRed); | |
1321 | chi2Result->DrawCopy("SAME"); | |
1322 | bayesResult->SetLineColor(kBlue); | |
1323 | bayesResult->DrawCopy("SAME"); | |
1324 | gPad->SetLogy(); | |
1325 | ||
1326 | if (i < 4) | |
1327 | { | |
1328 | canvasA->cd(i+5); | |
1329 | } | |
1330 | else | |
1331 | canvasB->cd(i+5-4); | |
1332 | ||
1333 | chi2Result->Divide(chi2Result, mc, 1, 1, "B"); | |
1334 | bayesResult->Divide(bayesResult, mc, 1, 1, "B"); | |
1335 | ||
1336 | // skip errors for now | |
1337 | for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j) | |
1338 | { | |
1339 | chi2Result->SetBinError(j, 0); | |
1340 | bayesResult->SetBinError(j, 0); | |
1341 | } | |
1342 | ||
1343 | chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC"); | |
1344 | chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5); | |
1345 | ||
1346 | chi2Result->DrawCopy(""); | |
1347 | bayesResult->DrawCopy("SAME"); | |
1348 | } | |
1349 | ||
1350 | canvas3->cd(1); | |
1351 | legend->Draw(); | |
1352 | ||
1353 | canvasA->SaveAs(Form("%s.gif", canvasA->GetName())); | |
1354 | canvasB->SaveAs(Form("%s.gif", canvasB->GetName())); | |
1355 | ||
1356 | TCanvas* canvas2 = new TCanvas("DifferentSamples2", "DifferentSamples2", 800, 400); | |
1357 | canvas2->Divide(2, 1); | |
1358 | ||
1359 | canvas2->cd(1); | |
1360 | fitResultsChi2->SetMarkerStyle(20); | |
1361 | fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max); | |
1362 | fitResultsChi2->Draw("AP"); | |
1363 | ||
1364 | fitResultsBayes->SetMarkerStyle(3); | |
1365 | fitResultsBayes->SetMarkerColor(2); | |
1366 | fitResultsBayes->Draw("P SAME"); | |
1367 | ||
1368 | gPad->SetLogy(); | |
1369 | ||
1370 | canvas2->cd(2); | |
1371 | fitResultsChi2Limit->SetMarkerStyle(20); | |
1372 | fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax())); | |
1373 | fitResultsChi2Limit->Draw("AP"); | |
1374 | ||
1375 | fitResultsBayesLimit->SetMarkerStyle(3); | |
1376 | fitResultsBayesLimit->SetMarkerColor(2); | |
1377 | fitResultsBayesLimit->Draw("P SAME"); | |
1378 | ||
1379 | canvas2->SaveAs(Form("%s.gif", canvas2->GetName())); | |
1380 | canvas3->SaveAs(Form("%s.gif", canvas3->GetName())); | |
1381 | canvas4->SaveAs(Form("%s.gif", canvas4->GetName())); | |
1382 | } | |
1383 | ||
cfc19dd5 | 1384 | void Merge(Int_t n, const char** files, const char* output) |
1385 | { | |
447c325d | 1386 | // const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_100k_2.root", "multiplicityMC_100k_3.root", "multiplicityMC_100k_4.root", "multiplicityMC_100k_5.root", "multiplicityMC_100k_6.root", "multiplicityMC_100k_7.root", "multiplicityMC_100k_8.root" }; |
1387 | ||
1388 | ||
cfc19dd5 | 1389 | gSystem->Load("libPWG0base"); |
1390 | ||
1391 | AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n]; | |
1392 | TList list; | |
1393 | for (Int_t i=0; i<n; ++i) | |
1394 | { | |
1395 | TString name("Multiplicity"); | |
1396 | if (i > 0) | |
1397 | name.Form("Multiplicity%d", i); | |
1398 | ||
1399 | TFile::Open(files[i]); | |
1400 | data[i] = new AliMultiplicityCorrection(name, name); | |
1401 | data[i]->LoadHistograms("Multiplicity"); | |
1402 | if (i > 0) | |
1403 | list.Add(data[i]); | |
1404 | } | |
1405 | ||
1406 | data[0]->Merge(&list); | |
1407 | ||
447c325d | 1408 | //data[0]->DrawHistograms(); |
cfc19dd5 | 1409 | |
1410 | TFile::Open(output, "RECREATE"); | |
1411 | data[0]->SaveHistograms(); | |
1412 | gFile->Close(); | |
1413 | } | |
1414 | ||
1415 | void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root") | |
1416 | { | |
1417 | gSystem->Load("libPWG0base"); | |
1418 | ||
1419 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
1420 | ||
1421 | TFile::Open(fileName); | |
1422 | mult->LoadHistograms("Multiplicity"); | |
1423 | ||
1424 | TF1* func = 0; | |
1425 | ||
1426 | if (caseNo >= 4) | |
1427 | { | |
0b4bfd98 | 1428 | func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 500); |
cfc19dd5 | 1429 | func->SetParNames("scaling", "averagen", "k"); |
1430 | } | |
1431 | ||
1432 | switch (caseNo) | |
1433 | { | |
0b4bfd98 | 1434 | case 0: func = new TF1("flat", "1000"); break; |
cfc19dd5 | 1435 | case 1: func = new TF1("flat", "501-x"); break; |
1436 | case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break; | |
1437 | case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break; | |
dd701109 | 1438 | case 4: func->SetParameters(1e7, 10, 2); break; |
447c325d | 1439 | case 5: func->SetParameters(1, 13, 7); break; |
dd701109 | 1440 | case 6: func->SetParameters(1e7, 30, 4); break; |
0b4bfd98 | 1441 | case 7: func->SetParameters(1e7, 30, 2); break; // *** |
dd701109 | 1442 | case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break; |
cfc19dd5 | 1443 | |
1444 | default: return; | |
1445 | } | |
1446 | ||
447c325d | 1447 | new TCanvas; |
1448 | func->Draw(); | |
1449 | ||
1450 | mult->SetGenMeasFromFunc(func, 3); | |
cfc19dd5 | 1451 | |
dd701109 | 1452 | TFile::Open("out.root", "RECREATE"); |
1453 | mult->SaveHistograms(); | |
1454 | ||
0b4bfd98 | 1455 | new TCanvas; mult->GetMultiplicityESD(3)->ProjectionY()->DrawCopy(); |
1456 | new TCanvas; mult->GetMultiplicityVtx(3)->ProjectionY()->DrawCopy(); | |
1457 | ||
dd701109 | 1458 | //mult->ApplyBayesianMethod(2, kFALSE); |
1459 | //mult->ApplyMinuitFit(2, kFALSE); | |
cfc19dd5 | 1460 | //mult->ApplyGaussianMethod(2, kFALSE); |
447c325d | 1461 | //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx); |
dd701109 | 1462 | } |
1463 | ||
1464 | void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2) | |
1465 | { | |
1466 | gSystem->Load("libPWG0base"); | |
1467 | ||
1468 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
1469 | ||
1470 | TFile::Open(fileName); | |
1471 | mult->LoadHistograms("Multiplicity"); | |
1472 | ||
1473 | // empty under/overflow bins in x, otherwise Project3D takes them into account | |
1474 | TH3* corr = mult->GetCorrelation(corrMatrix); | |
1475 | for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j) | |
1476 | { | |
1477 | for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k) | |
1478 | { | |
1479 | corr->SetBinContent(0, j, k, 0); | |
1480 | corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0); | |
1481 | } | |
1482 | } | |
1483 | ||
1484 | TH2* proj = (TH2*) corr->Project3D("zy"); | |
1485 | ||
1486 | // normalize correction for given nPart | |
1487 | for (Int_t i=1; i<=proj->GetNbinsX(); ++i) | |
1488 | { | |
1489 | Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY()); | |
1490 | if (sum <= 0) | |
1491 | continue; | |
1492 | ||
1493 | for (Int_t j=1; j<=proj->GetNbinsY(); ++j) | |
1494 | { | |
1495 | // npart sum to 1 | |
1496 | proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum); | |
1497 | proj->SetBinError(i, j, proj->GetBinError(i, j) / sum); | |
1498 | } | |
1499 | } | |
1500 | ||
1501 | new TCanvas; | |
447c325d | 1502 | proj->DrawCopy("COLZ"); |
dd701109 | 1503 | |
1504 | TH1* scaling = proj->ProjectionY("scaling", 1, 1); | |
1505 | scaling->Reset(); | |
1506 | scaling->SetMarkerStyle(3); | |
1507 | //scaling->GetXaxis()->SetRangeUser(0, 50); | |
1508 | TH1* mean = (TH1F*) scaling->Clone("mean"); | |
1509 | TH1* width = (TH1F*) scaling->Clone("width"); | |
1510 | ||
1511 | TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500); | |
1512 | lognormal->SetParNames("scaling", "mean", "sigma"); | |
1513 | lognormal->SetParameters(1, 1, 1); | |
1514 | lognormal->SetParLimits(0, 1, 1); | |
1515 | lognormal->SetParLimits(1, 0, 100); | |
1516 | lognormal->SetParLimits(2, 1e-3, 1); | |
1517 | ||
1518 | TF1* nbd = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50); | |
1519 | nbd->SetParNames("scaling", "averagen", "k"); | |
1520 | nbd->SetParameters(1, 13, 5); | |
1521 | nbd->SetParLimits(0, 1, 1); | |
1522 | nbd->SetParLimits(1, 1, 100); | |
1523 | nbd->SetParLimits(2, 1, 1e8); | |
1524 | ||
1525 | TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50); | |
1526 | poisson->SetParNames("scaling", "k", "deltax"); | |
1527 | poisson->SetParameters(1, 1, 0); | |
1528 | poisson->SetParLimits(0, 0, 10); | |
1529 | poisson->SetParLimits(1, 0.01, 100); | |
1530 | poisson->SetParLimits(2, 0, 10); | |
1531 | ||
1532 | TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50); | |
1533 | mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth"); | |
1534 | mygaus->SetParameters(1, 0, 1, 1, 0, 1); | |
1535 | mygaus->SetParLimits(2, 1e-5, 10); | |
1536 | mygaus->SetParLimits(4, 1, 1); | |
1537 | mygaus->SetParLimits(5, 1e-5, 10); | |
1538 | ||
1539 | //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50); | |
1540 | TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50); | |
1541 | sqrt->SetParNames("ydelta", "exp", "xdelta"); | |
1542 | sqrt->SetParameters(0, 0, 1); | |
1543 | sqrt->SetParLimits(1, 0, 10); | |
1544 | ||
1545 | const char* fitWith = "gaus"; | |
1546 | ||
1547 | for (Int_t i=1; i<=150; ++i) | |
1548 | { | |
1549 | printf("Fitting %d...\n", i); | |
1550 | ||
1551 | TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e"); | |
447c325d | 1552 | |
dd701109 | 1553 | //hist->GetXaxis()->SetRangeUser(0, 50); |
1554 | //lognormal->SetParameter(0, hist->GetMaximum()); | |
1555 | hist->Fit(fitWith, "0 M", ""); | |
1556 | ||
1557 | TF1* func = hist->GetListOfFunctions()->FindObject(fitWith); | |
1558 | ||
447c325d | 1559 | if (0 && (i % 5 == 0)) |
dd701109 | 1560 | { |
447c325d | 1561 | pad = new TCanvas; |
dd701109 | 1562 | hist->Draw(); |
1563 | func->Clone()->Draw("SAME"); | |
447c325d | 1564 | pad->SetLogy(); |
dd701109 | 1565 | } |
1566 | ||
1567 | scaling->Fill(i, func->GetParameter(0)); | |
1568 | mean->Fill(i, func->GetParameter(1)); | |
1569 | width->Fill(i, func->GetParameter(2)); | |
1570 | } | |
1571 | ||
1572 | TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500); | |
1573 | log->SetParameters(0, 1, 1); | |
1574 | log->SetParLimits(1, 0, 100); | |
1575 | log->SetParLimits(2, 1e-3, 10); | |
1576 | ||
1577 | TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500); | |
1578 | over->SetParameters(0, 1, 0); | |
1579 | //over->SetParLimits(0, 0, 100); | |
1580 | over->SetParLimits(1, 1e-3, 10); | |
1581 | over->SetParLimits(2, 0, 100); | |
1582 | ||
1583 | c1 = new TCanvas("fitparams", "fitparams", 1200, 400); | |
1584 | c1->Divide(3, 1); | |
1585 | ||
1586 | c1->cd(1); | |
1587 | scaling->Draw("P"); | |
1588 | ||
1589 | //TF1* scalingFit = new TF1("mypol0", "[0]"); | |
1590 | TF1* scalingFit = over; | |
447c325d | 1591 | scaling->Fit(scalingFit, "", "", 3, 140); |
1592 | scalingFit->SetRange(0, 200); | |
1593 | scalingFit->Draw("SAME"); | |
dd701109 | 1594 | |
1595 | c1->cd(2); | |
1596 | mean->Draw("P"); | |
1597 | ||
1598 | //TF1* meanFit = log; | |
1599 | TF1* meanFit = new TF1("mypol1", "[0]+[1]*x"); | |
447c325d | 1600 | mean->Fit(meanFit, "", "", 3, 140); |
1601 | meanFit->SetRange(0, 200); | |
1602 | meanFit->Draw("SAME"); | |
dd701109 | 1603 | |
1604 | c1->cd(3); | |
1605 | width->Draw("P"); | |
1606 | ||
1607 | //TF1* widthFit = over; | |
447c325d | 1608 | TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)"); |
1609 | widthFit->SetParLimits(2, 1e-5, 1e5); | |
1610 | width->Fit(widthFit, "", "", 5, 140); | |
1611 | widthFit->SetRange(0, 200); | |
1612 | widthFit->Draw("SAME"); | |
dd701109 | 1613 | |
1614 | // build new correction matrix | |
1615 | TH2* new = (TH2*) proj->Clone("new"); | |
1616 | new->Reset(); | |
1617 | Float_t x, y; | |
1618 | for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) | |
1619 | { | |
1620 | TF1* func = (TF1*) gROOT->FindObject(fitWith); | |
1621 | x = new->GetXaxis()->GetBinCenter(i); | |
1622 | //if (i == 1) | |
1623 | // x = 0.1; | |
1624 | x++; | |
1625 | func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); | |
1626 | printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x)); | |
1627 | ||
1628 | for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) | |
1629 | { | |
447c325d | 1630 | if (i < 11) |
dd701109 | 1631 | { |
1632 | // leave bins 1..20 untouched | |
1633 | new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j)); | |
1634 | } | |
1635 | else | |
1636 | { | |
1637 | y = new->GetYaxis()->GetBinCenter(j); | |
1638 | if (j == 1) | |
1639 | y = 0.1; | |
1640 | if (func->Eval(y) > 1e-4) | |
1641 | new->SetBinContent(i, j, func->Eval(y)); | |
1642 | } | |
1643 | } | |
1644 | } | |
1645 | ||
1646 | // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0 | |
1647 | // we take the values from the old response matrix | |
1648 | //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1) | |
1649 | // new->SetBinContent(i, 1, proj->GetBinContent(i, 1)); | |
1650 | ||
1651 | //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1) | |
1652 | // new->SetBinContent(1, j, proj->GetBinContent(1, j)); | |
1653 | ||
1654 | // normalize correction for given nPart | |
1655 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) | |
1656 | { | |
1657 | Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY()); | |
1658 | if (sum <= 0) | |
1659 | continue; | |
1660 | ||
1661 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) | |
1662 | { | |
1663 | // npart sum to 1 | |
1664 | new->SetBinContent(i, j, new->GetBinContent(i, j) / sum); | |
1665 | new->SetBinError(i, j, new->GetBinError(i, j) / sum); | |
1666 | } | |
1667 | } | |
1668 | ||
1669 | new TCanvas; | |
1670 | new->Draw("COLZ"); | |
1671 | ||
1672 | TH2* diff = (TH2*) new->Clone("diff"); | |
1673 | diff->Add(proj, -1); | |
1674 | ||
1675 | new TCanvas; | |
1676 | diff->Draw("COLZ"); | |
1677 | diff->SetMinimum(-0.05); | |
1678 | diff->SetMaximum(0.05); | |
1679 | ||
1680 | corr->Reset(); | |
1681 | ||
1682 | for (Int_t i=1; i<=new->GetNbinsX(); ++i) | |
1683 | for (Int_t j=1; j<=new->GetNbinsY(); ++j) | |
1684 | corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j)); | |
1685 | ||
1686 | new TCanvas; | |
1687 | corr->Project3D("zy")->Draw("COLZ"); | |
1688 | ||
1689 | TFile::Open("out.root", "RECREATE"); | |
1690 | mult->SaveHistograms(); | |
1691 | ||
1692 | TH1* proj1 = proj->ProjectionY("proj1", 36, 36); | |
1693 | TH1* proj2 = new->ProjectionY("proj2", 36, 36); | |
1694 | proj2->SetLineColor(2); | |
1695 | ||
1696 | new TCanvas; | |
1697 | proj1->Draw(); | |
1698 | proj2->Draw("SAME"); | |
cfc19dd5 | 1699 | } |
447c325d | 1700 | |
0b4bfd98 | 1701 | void buildCorrelationMap(const char* fileName = "multiplicityMC_2M.root", Int_t corrMatrix = 3) |
1702 | { | |
1703 | gSystem->Load("libPWG0base"); | |
1704 | ||
1705 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
1706 | ||
1707 | TFile::Open(fileName); | |
1708 | mult->LoadHistograms("Multiplicity"); | |
1709 | ||
1710 | TH3F* new = mult->GetCorrelation(corrMatrix); | |
1711 | new->Reset(); | |
1712 | ||
1713 | TF1* func = new TF1("func", "gaus(0)"); | |
1714 | ||
1715 | Int_t vtxBin = new->GetNbinsX() / 2; | |
1716 | if (vtxBin == 0) | |
1717 | vtxBin = 1; | |
1718 | ||
1719 | Float_t sigma = 2; | |
1720 | for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1) | |
1721 | { | |
1722 | Float_t x = new->GetYaxis()->GetBinCenter(i); | |
1723 | func->SetParameters(1, x * 0.8, sigma); | |
1724 | //func->SetParameters(1, x, sigma); | |
1725 | ||
1726 | for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1) | |
1727 | { | |
1728 | Float_t y = new->GetYaxis()->GetBinCenter(j); | |
1729 | ||
1730 | // cut at 1 sigma | |
1731 | if (TMath::Abs(y-x*0.8) < sigma) | |
1732 | new->SetBinContent(vtxBin, i, j, func->Eval(y)); | |
1733 | ||
1734 | // test only bin 40 has smearing | |
1735 | //if (x != 40) | |
1736 | // new->SetBinContent(vtxBin, i, j, (i == j)); | |
1737 | } | |
1738 | } | |
1739 | ||
1740 | new TCanvas; | |
1741 | new->Project3D("zy")->DrawCopy("COLZ"); | |
1742 | ||
1743 | TFile* file = TFile::Open("out.root", "RECREATE"); | |
1744 | mult->SetCorrelation(corrMatrix, new); | |
1745 | mult->SaveHistograms(); | |
1746 | file->Close(); | |
1747 | } | |
1748 | ||
447c325d | 1749 | void GetCrossSections(const char* fileName) |
1750 | { | |
1751 | gSystem->Load("libPWG0base"); | |
1752 | ||
1753 | AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); | |
1754 | ||
1755 | TFile::Open(fileName); | |
1756 | mult->LoadHistograms("Multiplicity"); | |
1757 | ||
1758 | TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2"); | |
1759 | xSection2->Sumw2(); | |
1760 | xSection2->Scale(1.0 / xSection2->Integral()); | |
1761 | ||
1762 | TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15"); | |
1763 | xSection15->Sumw2(); | |
1764 | xSection15->Scale(1.0 / xSection15->Integral()); | |
1765 | ||
1766 | TFile::Open("crosssection.root", "RECREATE"); | |
1767 | xSection2->Write(); | |
1768 | xSection15->Write(); | |
1769 | gFile->Close(); | |
1770 | } | |
0b4bfd98 | 1771 | |
1772 | void BuildResponseFromTree(const char* fileName, const char* target) | |
1773 | { | |
1774 | // | |
1775 | // builds several response matrices with different particle ratios (systematic study) | |
1776 | // | |
1777 | ||
1778 | gSystem->Load("libPWG0base"); | |
1779 | ||
1780 | TFile::Open(fileName); | |
1781 | TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies"); | |
1782 | ||
1783 | TFile* file = TFile::Open(target, "RECREATE"); | |
1784 | file->Close(); | |
1785 | ||
1786 | Int_t tracks = 0; // control variables | |
1787 | Int_t noLabel = 0; | |
6d81c2de | 1788 | Int_t secondaries = 0; |
0b4bfd98 | 1789 | Int_t doubleCount = 0; |
1790 | ||
6d81c2de | 1791 | for (Int_t num = 0; num < 1; num++) |
0b4bfd98 | 1792 | { |
1793 | AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(Form("Multiplicity_%d", num), Form("Multiplicity_%d", num)); | |
1794 | ||
1795 | Float_t ratio[4]; // pi, K, p, other | |
1796 | for (Int_t i = 0; i < 4; i++) | |
1797 | ratio[i] = 1; | |
1798 | ||
1799 | switch (num) | |
1800 | { | |
1801 | case 1 : ratio[1] = 0.5; break; | |
1802 | case 2 : ratio[2] = 0.5; break; | |
1803 | case 3 : ratio[1] = 1.5; break; | |
1804 | case 4 : ratio[2] = 1.5; break; | |
1805 | case 5 : ratio[1] = 0.5; ratio[2] = 0.5; break; | |
1806 | case 6 : ratio[1] = 1.5; ratio[2] = 1.5; break; | |
1807 | } | |
1808 | ||
1809 | for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++) | |
1810 | { | |
1811 | fParticleSpecies->GetEvent(i); | |
1812 | ||
1813 | Float_t* f = fParticleSpecies->GetArgs(); | |
1814 | ||
1815 | Float_t gene = 0; | |
1816 | Float_t meas = 0; | |
1817 | ||
1818 | for (Int_t j = 0; j < 4; j++) | |
1819 | { | |
1820 | gene += ratio[j] * f[j+1]; | |
1821 | meas += ratio[j] * f[j+1+4]; | |
1822 | tracks += f[j+1+4]; | |
1823 | } | |
1824 | ||
1825 | // add the ones w/o label | |
1826 | tracks += f[9]; | |
1827 | noLabel += f[9]; | |
1828 | ||
6d81c2de | 1829 | // secondaries are already part of meas! |
1830 | secondaries += f[10]; | |
1831 | ||
1832 | // double counted are already part of meas! | |
1833 | doubleCount += f[11]; | |
0b4bfd98 | 1834 | |
6d81c2de | 1835 | // ones w/o label are added without weighting to allow comparison to default analysis. however this is only valid when their fraction is low enough! |
0b4bfd98 | 1836 | meas += f[9]; |
0b4bfd98 | 1837 | |
1838 | //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]); | |
1839 | ||
6d81c2de | 1840 | fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, 0, meas, meas, meas, meas); |
1841 | fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, gene, gene, gene, gene, 0); | |
1842 | fMultiplicity->FillMeasured(f[0], meas, meas, meas, meas); | |
0b4bfd98 | 1843 | } |
1844 | ||
1845 | //fMultiplicity->DrawHistograms(); | |
1846 | ||
1847 | TFile* file = TFile::Open(target, "UPDATE"); | |
1848 | fMultiplicity->SaveHistograms(); | |
1849 | file->Close(); | |
1850 | ||
1851 | if (num == 0) | |
6d81c2de | 1852 | { |
1853 | Printf("%d total tracks, %d w/o label = %.2f %%, %d double counted = %.2f %%, secondaries = %.2f %%", tracks, noLabel, 100.0 * noLabel / tracks, doubleCount, 100.0 * doubleCount / tracks, 100.0 * secondaries / tracks); | |
1854 | if ((Float_t) noLabel / tracks > 0.02) | |
1855 | Printf("WARNING: More than 2%% of tracks without label, this might bias the study!"); | |
1856 | } | |
0b4bfd98 | 1857 | } |
1858 | } | |
1859 | ||
1860 | void MergeModifyCrossSection(const char* output) | |
1861 | { | |
1862 | const char* files[] = { "multiplicityMC_400k_syst_nd.root", "multiplicityMC_400k_syst_sd.root", "multiplicityMC_400k_syst_dd.root" }; | |
1863 | ||
1864 | gSystem->Load("libPWG0base"); | |
1865 | ||
1866 | TFile::Open(output, "RECREATE"); | |
1867 | gFile->Close(); | |
1868 | ||
1869 | for (Int_t num=0; num<7; ++num) | |
1870 | { | |
1871 | AliMultiplicityCorrection* data[3]; | |
1872 | TList list; | |
1873 | ||
1874 | Float_t ratio[3]; | |
1875 | switch (num) | |
1876 | { | |
1877 | case 0: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.0; break; | |
1878 | case 1: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.0; break; | |
1879 | case 2: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 1.0; break; | |
1880 | case 3: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 1.5; break; | |
1881 | case 4: ratio[0] = 1.0; ratio[1] = 1.0; ratio[2] = 0.5; break; | |
1882 | case 5: ratio[0] = 1.0; ratio[1] = 1.5; ratio[2] = 1.5; break; | |
1883 | case 6: ratio[0] = 1.0; ratio[1] = 0.5; ratio[2] = 0.5; break; | |
1884 | default: return; | |
1885 | } | |
1886 | ||
1887 | for (Int_t i=0; i<3; ++i) | |
1888 | { | |
1889 | TString name; | |
1890 | name.Form("Multiplicity_%d", num); | |
1891 | if (i > 0) | |
1892 | name.Form("Multiplicity_%d_%d", num, i); | |
1893 | ||
1894 | TFile::Open(files[i]); | |
1895 | data[i] = new AliMultiplicityCorrection(name, name); | |
1896 | data[i]->LoadHistograms("Multiplicity"); | |
1897 | ||
1898 | // modify x-section | |
1899 | for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++) | |
1900 | { | |
1901 | data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]); | |
1902 | data[i]->GetMultiplicityMB(j)->Scale(ratio[i]); | |
1903 | data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]); | |
1904 | } | |
1905 | ||
1906 | for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++) | |
1907 | data[i]->GetMultiplicityESD(j)->Scale(ratio[i]); | |
1908 | ||
1909 | for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++) | |
1910 | data[i]->GetCorrelation(j)->Scale(ratio[i]); | |
1911 | ||
1912 | if (i > 0) | |
1913 | list.Add(data[i]); | |
1914 | } | |
1915 | ||
1916 | printf("Case %d, %s: Entries in response matrix 3: ND: %.2f SD: %.2f DD: %.2f", num, data[0]->GetName(), data[0]->GetCorrelation(3)->Integral(), data[1]->GetCorrelation(3)->Integral(), data[2]->GetCorrelation(3)->Integral()); | |
1917 | ||
1918 | data[0]->Merge(&list); | |
1919 | ||
1920 | Printf(" Total: %.2f", data[0]->GetCorrelation(3)->Integral()); | |
1921 | ||
1922 | TFile::Open(output, "UPDATE"); | |
1923 | data[0]->SaveHistograms(); | |
1924 | gFile->Close(); | |
1925 | ||
1926 | list.Clear(); | |
1927 | ||
1928 | for (Int_t i=0; i<3; ++i) | |
1929 | delete data[i]; | |
1930 | } | |
1931 | } |