git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Modifications to take into account new muon esd structure (P. Cortese)
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / runMultiplicitySelector.C
CommitLineData
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 10void 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)
cfc19dd5 13 connectProof("jgrosseo@lxb6046");
0bd1f8a0 14
15 TString libraries("libEG;libGeom;libESD;libPWG0base");
16 TString packages("PWG0base");
10ebe68d 17
0ab29cfa 18 if (aMC != kFALSE)
0bd1f8a0 19 {
20 libraries += ";libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6";
21 packages += ";PWG0dep";
22 }
23
24 if (!prepareQuery(libraries, packages, kTRUE))
25 return;
0ab29cfa 26
0ab29cfa 27 gROOT->ProcessLine(".L CreateCuts.C");
28 gROOT->ProcessLine(".L drawPlots.C");
29
0ab29cfa 30 // selection of esd tracks
31 AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
32 if (!esdTrackCuts)
33 {
34 printf("ERROR: esdTrackCuts could not be created\n");
35 return;
36 }
37
0bd1f8a0 38 TList inputList;
39 inputList.Add(esdTrackCuts);
40
cfc19dd5 41 TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE);
0ab29cfa 42
43 TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector");
44 AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
45
0bd1f8a0 46 selectorName += ".cxx+";
0ab29cfa 47
48 if (aDebug != kFALSE)
49 selectorName += "g";
50
03244034 51 Int_t result = executeQuery(chain, &inputList, selectorName, option);
0ab29cfa 52
53 if (result != 0)
54 {
55 printf("ERROR: Executing process failed with %d.\n", result);
56 return;
57 }
0a173978 58}
0ab29cfa 59
0a173978 60void draw(const char* fileName = "multiplicityMC.root")
61{
62 gSystem->Load("libPWG0base");
63
64 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
65
66 TFile::Open(fileName);
67 mult->LoadHistograms("Multiplicity");
68
69 mult->DrawHistograms();
0ab29cfa 70}
71
cfc19dd5 72void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
0a173978 73{
74 gSystem->Load("libPWG0base");
75
76 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
77
78 TFile::Open(fileName);
79 mult->LoadHistograms("Multiplicity");
80
dd701109 81 //mult->ApplyLaszloMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
82
83 //return;
84
85
86 mult->ApplyBayesianMethod(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
87 mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
88
89 return;
90
91 TStopwatch timer;
cfc19dd5 92
dd701109 93 timer.Start();
94
95 mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1);
96 mult->ApplyMinuitFit(hist, kFALSE, AliMultiplicityCorrection::kTrVtx);
97 mult->DrawComparison("MinuitChi2", hist, kFALSE, kFALSE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
98
99 timer.Stop();
100 timer.Print();
101
102 return 0;
cfc19dd5 103
104 //mult->ApplyGaussianMethod(hist, kFALSE);
9ca6be09 105
dd701109 106 mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0);
107 mult->ApplyNBDFit(hist, kFALSE);
108 mult->DrawComparison("NBDChi2Fit", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY());
109
9ca6be09 110 return mult;
111}
112
dd701109 113void* fitOther(const char* fileNameMC = "multiplicityMC_3M.root", const char* fileNameESD = "multiplicityMC_3M.root", Int_t histID = 2)
9ca6be09 114{
115 gSystem->Load("libPWG0base");
116
117 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
118
119 TFile::Open(fileNameMC);
120 mult->LoadHistograms("Multiplicity");
121
122 TFile::Open(fileNameESD);
cfc19dd5 123 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
124 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
dd701109 125 //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityMB%d", histID));
9ca6be09 126
cfc19dd5 127 mult->SetMultiplicityESD(histID, hist);
9ca6be09 128
dd701109 129 mult->SetMultiplicityVtx(histID, hist2);
130 mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
131 return;
132
133 mult->SetRegularizationParameters(AliMultiplicityCorrection::kTest, 1.1);
134 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
135 mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
136
137 return;
138
cfc19dd5 139 //mult->ApplyGaussianMethod(histID, kFALSE);
cfc19dd5 140
dd701109 141 for (Float_t f=0.1; f<=0.11; f+=0.05)
142 {
143 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, f);
144 mult->DrawComparison(Form("Bayesian_%f", f), histID, kFALSE, kTRUE, hist2->ProjectionY());
145 }
146
147 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 1e7);
148 //mult->ApplyMinuitFit(histID, kFALSE);
149 //mult->DrawComparison("MinuitChi2", histID, kFALSE, kTRUE, hist2->ProjectionY());
150
0a173978 151
152 return mult;
153}
cfc19dd5 154
155const char* GetRegName(Int_t type)
156{
157 switch (type)
158 {
159 case AliMultiplicityCorrection::kNone: return "None"; break;
160 case AliMultiplicityCorrection::kPol0: return "Pol0"; break;
161 case AliMultiplicityCorrection::kPol1: return "Pol1"; break;
162 case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
163 case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break;
dd701109 164 case AliMultiplicityCorrection::kTest : return "Test"; break;
165 }
166 return 0;
167}
168
169const char* GetEventTypeName(Int_t type)
170{
171 switch (type)
172 {
173 case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break;
174 case AliMultiplicityCorrection::kMB: return "minimum bias"; break;
175 case AliMultiplicityCorrection::kINEL: return "inelastic"; break;
cfc19dd5 176 }
177 return 0;
178}
179
dd701109 180void EvaluateBayesianMethod(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
cfc19dd5 181{
dd701109 182 gSystem->mkdir(targetDir);
183
184 gSystem->Load("libPWG0base");
185
186 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
187 TFile::Open(fileNameMC);
188 mult->LoadHistograms("Multiplicity");
189
190 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
191 TFile::Open(fileNameESD);
192 multESD->LoadHistograms("Multiplicity");
193 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
194
195 TCanvas* canvas = new TCanvas("EvaluateBayesianMethod", "EvaluateBayesianMethod", 800, 600);
196 TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
197 legend->SetFillColor(0);
198
199 Float_t min = 1e10;
200 Float_t max = 0;
201
202 TGraph* first = 0;
203 Int_t count = 0; // just to order the saved images...
204
205 for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
206 {
207 TGraph* fitResultsMC = new TGraph;
208 fitResultsMC->SetTitle(";Weight Parameter");
209 TGraph* fitResultsRes = new TGraph;
210 fitResultsRes->SetTitle(";Weight Parameter");
211
212 fitResultsMC->SetFillColor(0);
213 fitResultsRes->SetFillColor(0);
214 fitResultsMC->SetMarkerStyle(20+type);
215 fitResultsRes->SetMarkerStyle(24+type);
216 fitResultsRes->SetMarkerColor(kRed);
217 fitResultsRes->SetLineColor(kRed);
218
219 legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
220 legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
221
222 for (Float_t weight = 0; weight < 0.301; weight += 0.02)
223 {
224 Float_t chi2MC = 0;
225 Float_t residuals = 0;
226
227 mult->ApplyBayesianMethod(histID, kFALSE, type, weight);
228 mult->DrawComparison(Form("%s/Bayesian_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
229 mult->GetComparisonResults(&chi2MC, 0, &residuals);
230
231 fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
232 fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
233
234 min = TMath::Min(min, TMath::Min(chi2MC, residuals));
235 max = TMath::Max(max, TMath::Max(chi2MC, residuals));
236 }
237
238 fitResultsMC->Print();
239 fitResultsRes->Print();
240
241 canvas->cd();
242 fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME"));
243 fitResultsRes->Draw("SAME CP");
244
245 if (first == 0)
246 first = fitResultsMC;
247 }
248
249 gPad->SetLogy();
250 printf("min = %f, max = %f\n", min, max);
251 if (min <= 0)
252 min = 1e-5;
253 first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
254
255 legend->Draw();
256
257 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
258 canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
259}
260
261void EvaluateBayesianMethodIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
262{
263 gSystem->mkdir(targetDir);
264
265 gSystem->Load("libPWG0base");
266
267 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
268 TFile::Open(fileNameMC);
269 mult->LoadHistograms("Multiplicity");
270
271 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
272 TFile::Open(fileNameESD);
273 multESD->LoadHistograms("Multiplicity");
274 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
275
276 TCanvas* canvas = new TCanvas("EvaluateBayesianMethodIterations", "EvaluateBayesianMethodIterations", 800, 600);
277 TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
278 legend->SetFillColor(0);
279
280 Float_t min = 1e10;
281 Float_t max = 0;
282
283 TGraph* first = 0;
284 Int_t count = 0; // just to order the saved images...
285
286 for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
287 {
288 TGraph* fitResultsMC = new TGraph;
289 fitResultsMC->SetTitle(";Iterations");
290 TGraph* fitResultsRes = new TGraph;
291 fitResultsRes->SetTitle(";Iterations");
292
293 fitResultsMC->SetFillColor(0);
294 fitResultsRes->SetFillColor(0);
295 fitResultsMC->SetMarkerStyle(20+type);
296 fitResultsRes->SetMarkerStyle(24+type);
297 fitResultsRes->SetMarkerColor(kRed);
298 fitResultsRes->SetLineColor(kRed);
299
300 legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetEventTypeName(type)));
301 legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetEventTypeName(type)));
302
303 for (Int_t iter = 5; iter <= 50; iter += 5)
304 {
305 Float_t chi2MC = 0;
306 Float_t residuals = 0;
307
308 mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1, iter);
309 mult->DrawComparison(Form("%s/BayesianIter_%02d_%d_%d", targetDir, count++, type, iter), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
310 mult->GetComparisonResults(&chi2MC, 0, &residuals);
311
312 fitResultsMC->SetPoint(fitResultsMC->GetN(), iter, chi2MC);
313 fitResultsRes->SetPoint(fitResultsRes->GetN(), iter, residuals);
314
315 min = TMath::Min(min, TMath::Min(chi2MC, residuals));
316 max = TMath::Max(max, TMath::Max(chi2MC, residuals));
317 }
318
319 fitResultsMC->Print();
320 fitResultsRes->Print();
321
322 canvas->cd();
323 fitResultsMC->Draw(Form("%s CP", (first == 0) ? "A" : "SAME"));
324 fitResultsRes->Draw("SAME CP");
325
326 if (first == 0)
327 first = fitResultsMC;
328 }
329
330 gPad->SetLogy();
331 printf("min = %f, max = %f\n", min, max);
332 if (min <= 0)
333 min = 1e-5;
334 first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
335
336 legend->Draw();
337
338 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
339 canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
340}
341
342void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
343{
344 gSystem->mkdir(targetDir);
345
cfc19dd5 346 gSystem->Load("libPWG0base");
347
348 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
349
350 TFile::Open(fileNameMC);
351 mult->LoadHistograms("Multiplicity");
352
353 TFile::Open(fileNameESD);
354 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
355 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
356
357 mult->SetMultiplicityESD(histID, hist);
358
359 TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600);
dd701109 360 TLegend* legend = new TLegend(0.2, 0.7, 0.58, 0.98);
cfc19dd5 361 legend->SetFillColor(0);
362
363 Float_t min = 1e10;
364 Float_t max = 0;
365
366 TGraph* first = 0;
dd701109 367 Int_t count = 0; // just to order the saved images...
368
369 Bool_t firstLoop = kTRUE;
cfc19dd5 370
371 for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
dd701109 372 //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
cfc19dd5 373 {
374 TGraph* fitResultsMC = new TGraph;
375 fitResultsMC->SetTitle(";Weight Parameter");
376 TGraph* fitResultsRes = new TGraph;
377 fitResultsRes->SetTitle(";Weight Parameter");
378
379 fitResultsMC->SetFillColor(0);
380 fitResultsRes->SetFillColor(0);
381 fitResultsMC->SetMarkerStyle(19+type);
dd701109 382 fitResultsRes->SetMarkerStyle(23+type);
cfc19dd5 383 fitResultsRes->SetMarkerColor(kRed);
384 fitResultsRes->SetLineColor(kRed);
385
386 legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetRegName(type)));
387 legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetRegName(type)));
388
389 if (first == 0)
390 first = fitResultsMC;
391
dd701109 392 for (Float_t weight = 1e-4; weight < 2e4; weight *= 10)
393 //for (Float_t weight = 0.1; weight < 10; weight *= TMath::Sqrt(TMath::Sqrt(10)))
cfc19dd5 394 {
395 Float_t chi2MC = 0;
396 Float_t residuals = 0;
397
398 mult->SetRegularizationParameters(type, weight);
dd701109 399 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
400 mult->DrawComparison(Form("%s/MinuitChi2_%02d_%d_%f", targetDir, count++, type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
401 mult->GetComparisonResults(&chi2MC, 0, &residuals);
cfc19dd5 402
403 fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
404 fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
405
406 min = TMath::Min(min, TMath::Min(chi2MC, residuals));
407 max = TMath::Max(max, TMath::Max(chi2MC, residuals));
408 }
409
410 fitResultsMC->Print();
411 fitResultsRes->Print();
412
413 canvas->cd();
dd701109 414 fitResultsMC->Draw(Form("%s CP", (firstLoop) ? "A" : "SAME"));
cfc19dd5 415 fitResultsRes->Draw("SAME CP");
dd701109 416
417 firstLoop = kFALSE;
cfc19dd5 418 }
419
420 gPad->SetLogx();
421 gPad->SetLogy();
422 printf("min = %f, max = %f\n", min, max);
423 if (min <= 0)
424 min = 1e-5;
425 first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
426
427 legend->Draw();
428
dd701109 429 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
430 canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
431}
432
433void EvaluateChi2MethodAll()
434{
435 EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
436 EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
437 EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
438 EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
439 EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
440}
441
442void EvaluateBayesianMethodAll()
443{
444 EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
445 EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
446 EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
447 EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
448 EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
449}
450
451void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 2)
452{
453 gSystem->mkdir(targetDir);
454
455 gSystem->Load("libPWG0base");
456
457 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
458
459 TFile::Open(fileNameMC);
460 mult->LoadHistograms("Multiplicity");
461
462 TFile::Open(fileNameESD);
463 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
464 multESD->LoadHistograms("Multiplicity");
465
466 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
467
468 TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 800);
469 canvas->Divide(3, 2);
470
471 Int_t count = 0;
472
473 for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
474 {
475 TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY();
476
477 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
478 mult->ApplyMinuitFit(histID, kFALSE, type);
479 mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc);
480 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
481
482 mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1);
483 mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc);
484 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
485
486 mc->GetXaxis()->SetRangeUser(0, 150);
487 chi2Result->GetXaxis()->SetRangeUser(0, 150);
488
489 // skip errors for now
490 for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
491 {
492 chi2Result->SetBinError(i, 0);
493 bayesResult->SetBinError(i, 0);
494 }
495
496 canvas->cd(++count);
497 mc->SetFillColor(kYellow);
498 mc->DrawCopy();
499 chi2Result->SetLineColor(kRed);
500 chi2Result->DrawCopy("SAME");
501 bayesResult->SetLineColor(kBlue);
502 bayesResult->DrawCopy("SAME");
503 gPad->SetLogy();
504
505 canvas->cd(count + 3);
506 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
507 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
508
509 // skip errors for now
510 for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
511 {
512 chi2Result->SetBinError(i, 0);
513 bayesResult->SetBinError(i, 0);
514 }
515
516 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
517
518 chi2Result->DrawCopy("");
519 bayesResult->DrawCopy("SAME");
520 }
521
522 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
523 canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
524}
525
526void StatisticsPlot(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
527{
528 gSystem->Load("libPWG0base");
529
530 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
531
532 TFile::Open(fileNameMC);
533 mult->LoadHistograms("Multiplicity");
534
535 const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
536
537 TGraph* fitResultsChi2 = new TGraph;
538 fitResultsChi2->SetTitle(";Nevents;Chi2");
539 TGraph* fitResultsBayes = new TGraph;
540 fitResultsBayes->SetTitle(";Nevents;Chi2");
541 TGraph* fitResultsChi2Limit = new TGraph;
542 fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
543 TGraph* fitResultsBayesLimit = new TGraph;
544 fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
545
546 TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
547 canvas->Divide(5, 2);
548
549 Float_t min = 1e10;
550 Float_t max = 0;
551
552 for (Int_t i=0; i<5; ++i)
553 {
554 TFile::Open(files[i]);
555 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
556 multESD->LoadHistograms("Multiplicity");
557
558 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
559 Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries();
560 TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
561
562 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
563 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
564 mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
565
566 Float_t chi2MC = 0;
567 Int_t chi2MCLimit = 0;
568 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
569 fitResultsChi2->SetPoint(fitResultsChi2->GetN(), nEntries, chi2MC);
570 fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit);
571 min = TMath::Min(min, chi2MC);
572 max = TMath::Max(max, chi2MC);
573
574 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
575
576 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1);
577 mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
578 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
579 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
580 fitResultsBayes->SetPoint(fitResultsBayes->GetN(), nEntries, chi2MC);
581 fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
582
583 min = TMath::Min(min, chi2MC);
584 max = TMath::Max(max, chi2MC);
585 mc->GetXaxis()->SetRangeUser(0, 150);
586 chi2Result->GetXaxis()->SetRangeUser(0, 150);
587
588 // skip errors for now
589 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
590 {
591 chi2Result->SetBinError(j, 0);
592 bayesResult->SetBinError(j, 0);
593 }
594
595 canvas->cd(i+1);
596 mc->SetFillColor(kYellow);
597 mc->DrawCopy();
598 chi2Result->SetLineColor(kRed);
599 chi2Result->DrawCopy("SAME");
600 bayesResult->SetLineColor(kBlue);
601 bayesResult->DrawCopy("SAME");
602 gPad->SetLogy();
603
604 canvas->cd(i+6);
605 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
606 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
607
608 // skip errors for now
609 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
610 {
611 chi2Result->SetBinError(j, 0);
612 bayesResult->SetBinError(j, 0);
613 }
614
615 chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
616 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
617
618 chi2Result->DrawCopy("");
619 bayesResult->DrawCopy("SAME");
620 }
621
622 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
623 canvas->SaveAs(Form("%s.C", canvas->GetName()));
624
625 TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400);
626 canvas2->Divide(2, 1);
627
628 canvas2->cd(1);
629 fitResultsChi2->SetMarkerStyle(20);
630 fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
631 fitResultsChi2->Draw("AP");
632
633 fitResultsBayes->SetMarkerStyle(3);
634 fitResultsBayes->SetMarkerColor(2);
635 fitResultsBayes->Draw("P SAME");
636
637 gPad->SetLogy();
638
639 canvas2->cd(2);
640 fitResultsChi2Limit->SetMarkerStyle(20);
641 fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
642 fitResultsChi2Limit->Draw("AP");
643
644 fitResultsBayesLimit->SetMarkerStyle(3);
645 fitResultsBayesLimit->SetMarkerColor(2);
646 fitResultsBayesLimit->Draw("P SAME");
647
648 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
649 canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
650}
651
652void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 2)
653{
654 gSystem->Load("libPWG0base");
655
656 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
657
658 TFile::Open(fileNameMC);
659 mult->LoadHistograms("Multiplicity");
660
661 const char* files[] = { "multiplicityMC_0.5M.root", "multiplicityMC_0.75M.root", "multiplicityMC_1M_3.root", "multiplicityMC_1.25M.root", "multiplicityMC_1.5M.root" };
662
663 // this one we try to unfold
664 TFile::Open(files[0]);
665 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
666 multESD->LoadHistograms("Multiplicity");
667 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
668 TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY();
669
670 TGraph* fitResultsChi2 = new TGraph;
671 fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
672 TGraph* fitResultsBayes = new TGraph;
673 fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
674 TGraph* fitResultsChi2Limit = new TGraph;
675 fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
676 TGraph* fitResultsBayesLimit = new TGraph;
677 fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
678
679 TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600);
680 canvas->Divide(8, 2);
681
682 TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400);
683 canvas3->Divide(2, 1);
684
685 Float_t min = 1e10;
686 Float_t max = 0;
687
688 TH1* firstChi = 0;
689 TH1* firstBayesian = 0;
690 TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond");
691
692 TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
693
694 for (Int_t i=0; i<8; ++i)
695 {
696 if (i == 0)
697 {
698 startCond = (TH1*) mc->Clone("startCond2");
699 }
700 else if (i < 6)
701 {
702 TFile::Open(files[i-1]);
703 AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
704 multESD2->LoadHistograms("Multiplicity");
705 startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond");
706 }
707 else if (i == 6)
708 {
709 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);
710 func->SetParNames("scaling", "averagen", "k");
711 func->SetParLimits(0, 1e-3, 1e10);
712 func->SetParLimits(1, 0.001, 1000);
713 func->SetParLimits(2, 0.001, 1000);
714
715 func->SetParameters(1, 10, 2);
716 for (Int_t j=2; j<=startCond->GetNbinsX(); j++)
717 startCond->SetBinContent(j, func->Eval(j-1));
718 }
719 else if (i == 7)
720 {
721 for (Int_t j=1; j<=startCond->GetNbinsX(); j++)
722 startCond->SetBinContent(j, 1);
723 }
724
725 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 0.3);
726 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
727 mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
728
729 Float_t chi2MC = 0;
730 Int_t chi2MCLimit = 0;
731 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
732 fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
733 fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
734 min = TMath::Min(min, chi2MC);
735 max = TMath::Max(max, chi2MC);
736
737 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
738 if (!firstChi)
739 firstChi = (TH1*) chi2Result->Clone("firstChi");
740
741 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 0.1, 30, startCond);
742 mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
743 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
744 if (!firstBayesian)
745 firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
746
747 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
748 fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
749 fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
750
751 min = TMath::Min(min, chi2MC);
752 max = TMath::Max(max, chi2MC);
753 mc->GetXaxis()->SetRangeUser(0, 150);
754 chi2Result->GetXaxis()->SetRangeUser(0, 150);
755
756 // skip errors for now
757 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
758 {
759 chi2Result->SetBinError(j, 0);
760 bayesResult->SetBinError(j, 0);
761 }
762
763 canvas3->cd(1);
764 TH1* tmp = (TH1*) chi2Result->Clone("tmp");
765 tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
766 tmp->Divide(firstChi);
767 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
768 tmp->GetXaxis()->SetRangeUser(0, 200);
769 tmp->SetLineColor(i+1);
770 legend->AddEntry(tmp, Form("%d", i));
771 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
772
773 canvas3->cd(2);
774 tmp = (TH1*) bayesResult->Clone("tmp");
775 tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
776 tmp->Divide(firstBayesian);
777 tmp->SetLineColor(i+1);
778 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
779 tmp->GetXaxis()->SetRangeUser(0, 200);
780 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
781
782 canvas->cd(i+1);
783 mc->SetFillColor(kYellow);
784 mc->DrawCopy();
785 chi2Result->SetLineColor(kRed);
786 chi2Result->DrawCopy("SAME");
787 bayesResult->SetLineColor(kBlue);
788 bayesResult->DrawCopy("SAME");
789 gPad->SetLogy();
790
791 canvas->cd(i+9);
792 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
793 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
794
795 // skip errors for now
796 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
797 {
798 chi2Result->SetBinError(j, 0);
799 bayesResult->SetBinError(j, 0);
800 }
801
802 chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
803 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
804
805 chi2Result->DrawCopy("");
806 bayesResult->DrawCopy("SAME");
807 }
808
809 canvas3->cd(1);
810 legend->Draw();
811
cfc19dd5 812 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
dd701109 813
814 TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400);
815 canvas2->Divide(2, 1);
816
817 canvas2->cd(1);
818 fitResultsChi2->SetMarkerStyle(20);
819 fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
820 fitResultsChi2->Draw("AP");
821
822 fitResultsBayes->SetMarkerStyle(3);
823 fitResultsBayes->SetMarkerColor(2);
824 fitResultsBayes->Draw("P SAME");
825
826 gPad->SetLogy();
827
828 canvas2->cd(2);
829 fitResultsChi2Limit->SetMarkerStyle(20);
830 fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
831 fitResultsChi2Limit->Draw("AP");
832
833 fitResultsBayesLimit->SetMarkerStyle(3);
834 fitResultsBayesLimit->SetMarkerColor(2);
835 fitResultsBayesLimit->Draw("P SAME");
836
837 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
838 canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
cfc19dd5 839}
840
841void Merge(Int_t n, const char** files, const char* output)
842{
843 gSystem->Load("libPWG0base");
844
845 AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
846 TList list;
847 for (Int_t i=0; i<n; ++i)
848 {
849 TString name("Multiplicity");
850 if (i > 0)
851 name.Form("Multiplicity%d", i);
852
853 TFile::Open(files[i]);
854 data[i] = new AliMultiplicityCorrection(name, name);
855 data[i]->LoadHistograms("Multiplicity");
856 if (i > 0)
857 list.Add(data[i]);
858 }
859
860 data[0]->Merge(&list);
861
862 data[0]->DrawHistograms();
863
864 TFile::Open(output, "RECREATE");
865 data[0]->SaveHistograms();
866 gFile->Close();
867}
868
869void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
870{
871 gSystem->Load("libPWG0base");
872
873 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
874
875 TFile::Open(fileName);
876 mult->LoadHistograms("Multiplicity");
877
878 TF1* func = 0;
879
880 if (caseNo >= 4)
881 {
882 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);
883 func->SetParNames("scaling", "averagen", "k");
884 }
885
886 switch (caseNo)
887 {
888 case 0: func = new TF1("flat", "1"); break;
889 case 1: func = new TF1("flat", "501-x"); break;
890 case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
891 case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
dd701109 892 case 4: func->SetParameters(1e7, 10, 2); break;
893 case 5: func->SetParameters(1e7, 20, 3); break;
894 case 6: func->SetParameters(1e7, 30, 4); break;
895 case 7: func->SetParameters(1e7, 70, 2); break;
896 case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
cfc19dd5 897
898 default: return;
899 }
900
901 mult->SetGenMeasFromFunc(func, 2);
902
dd701109 903 TFile::Open("out.root", "RECREATE");
904 mult->SaveHistograms();
905
906 //mult->ApplyBayesianMethod(2, kFALSE);
907 //mult->ApplyMinuitFit(2, kFALSE);
cfc19dd5 908 //mult->ApplyGaussianMethod(2, kFALSE);
dd701109 909 mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
910}
911
912void smoothCorrelationMap(const char* fileName = "multiplicityMC.root", Int_t corrMatrix = 2)
913{
914 gSystem->Load("libPWG0base");
915
916 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
917
918 TFile::Open(fileName);
919 mult->LoadHistograms("Multiplicity");
920
921 // empty under/overflow bins in x, otherwise Project3D takes them into account
922 TH3* corr = mult->GetCorrelation(corrMatrix);
923 for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j)
924 {
925 for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k)
926 {
927 corr->SetBinContent(0, j, k, 0);
928 corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0);
929 }
930 }
931
932 TH2* proj = (TH2*) corr->Project3D("zy");
933
934 // normalize correction for given nPart
935 for (Int_t i=1; i<=proj->GetNbinsX(); ++i)
936 {
937 Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY());
938 if (sum <= 0)
939 continue;
940
941 for (Int_t j=1; j<=proj->GetNbinsY(); ++j)
942 {
943 // npart sum to 1
944 proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum);
945 proj->SetBinError(i, j, proj->GetBinError(i, j) / sum);
946 }
947 }
948
949 new TCanvas;
950 proj->Draw("COLZ");
951
952 TH1* scaling = proj->ProjectionY("scaling", 1, 1);
953 scaling->Reset();
954 scaling->SetMarkerStyle(3);
955 //scaling->GetXaxis()->SetRangeUser(0, 50);
956 TH1* mean = (TH1F*) scaling->Clone("mean");
957 TH1* width = (TH1F*) scaling->Clone("width");
958
959 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
960 lognormal->SetParNames("scaling", "mean", "sigma");
961 lognormal->SetParameters(1, 1, 1);
962 lognormal->SetParLimits(0, 1, 1);
963 lognormal->SetParLimits(1, 0, 100);
964 lognormal->SetParLimits(2, 1e-3, 1);
965
966 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);
967 nbd->SetParNames("scaling", "averagen", "k");
968 nbd->SetParameters(1, 13, 5);
969 nbd->SetParLimits(0, 1, 1);
970 nbd->SetParLimits(1, 1, 100);
971 nbd->SetParLimits(2, 1, 1e8);
972
973 TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50);
974 poisson->SetParNames("scaling", "k", "deltax");
975 poisson->SetParameters(1, 1, 0);
976 poisson->SetParLimits(0, 0, 10);
977 poisson->SetParLimits(1, 0.01, 100);
978 poisson->SetParLimits(2, 0, 10);
979
980 TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50);
981 mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth");
982 mygaus->SetParameters(1, 0, 1, 1, 0, 1);
983 mygaus->SetParLimits(2, 1e-5, 10);
984 mygaus->SetParLimits(4, 1, 1);
985 mygaus->SetParLimits(5, 1e-5, 10);
986
987 //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50);
988 TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50);
989 sqrt->SetParNames("ydelta", "exp", "xdelta");
990 sqrt->SetParameters(0, 0, 1);
991 sqrt->SetParLimits(1, 0, 10);
992
993 const char* fitWith = "gaus";
994
995 for (Int_t i=1; i<=150; ++i)
996 {
997 printf("Fitting %d...\n", i);
998
999 TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e");
1000 //hist->GetXaxis()->SetRangeUser(0, 50);
1001 //lognormal->SetParameter(0, hist->GetMaximum());
1002 hist->Fit(fitWith, "0 M", "");
1003
1004 TF1* func = hist->GetListOfFunctions()->FindObject(fitWith);
1005
1006 if (((i-1) % 15 == 0) || ((i % 5 == 0) && i < 30))
1007 {
1008 new TCanvas;
1009 hist->Draw();
1010 func->Clone()->Draw("SAME");
1011 gPad->SetLogy();
1012 }
1013
1014 scaling->Fill(i, func->GetParameter(0));
1015 mean->Fill(i, func->GetParameter(1));
1016 width->Fill(i, func->GetParameter(2));
1017 }
1018
1019 TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500);
1020 log->SetParameters(0, 1, 1);
1021 log->SetParLimits(1, 0, 100);
1022 log->SetParLimits(2, 1e-3, 10);
1023
1024 TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500);
1025 over->SetParameters(0, 1, 0);
1026 //over->SetParLimits(0, 0, 100);
1027 over->SetParLimits(1, 1e-3, 10);
1028 over->SetParLimits(2, 0, 100);
1029
1030 c1 = new TCanvas("fitparams", "fitparams", 1200, 400);
1031 c1->Divide(3, 1);
1032
1033 c1->cd(1);
1034 scaling->Draw("P");
1035
1036 //TF1* scalingFit = new TF1("mypol0", "[0]");
1037 TF1* scalingFit = over;
1038 scaling->Fit(scalingFit, "", "", 3, 100);
1039
1040 c1->cd(2);
1041 mean->Draw("P");
1042
1043 //TF1* meanFit = log;
1044 TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
1045 mean->Fit(meanFit, "", "", 3, 100);
1046
1047 c1->cd(3);
1048 width->Draw("P");
1049
1050 //TF1* widthFit = over;
1051 TF1* widthFit = new TF1("mypol2", "[0]+[1]*x+[2]*x*x");
1052 width->Fit(widthFit, "", "", 5, 100);
1053
1054 // build new correction matrix
1055 TH2* new = (TH2*) proj->Clone("new");
1056 new->Reset();
1057 Float_t x, y;
1058 for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
1059 {
1060 TF1* func = (TF1*) gROOT->FindObject(fitWith);
1061 x = new->GetXaxis()->GetBinCenter(i);
1062 //if (i == 1)
1063 // x = 0.1;
1064 x++;
1065 func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
1066 printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
1067
1068 for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
1069 {
1070 if (i < 21)
1071 {
1072 // leave bins 1..20 untouched
1073 new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j));
1074 }
1075 else
1076 {
1077 y = new->GetYaxis()->GetBinCenter(j);
1078 if (j == 1)
1079 y = 0.1;
1080 if (func->Eval(y) > 1e-4)
1081 new->SetBinContent(i, j, func->Eval(y));
1082 }
1083 }
1084 }
1085
1086 // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0
1087 // we take the values from the old response matrix
1088 //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
1089 // new->SetBinContent(i, 1, proj->GetBinContent(i, 1));
1090
1091 //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
1092 // new->SetBinContent(1, j, proj->GetBinContent(1, j));
1093
1094 // normalize correction for given nPart
1095 for (Int_t i=1; i<=new->GetNbinsX(); ++i)
1096 {
1097 Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY());
1098 if (sum <= 0)
1099 continue;
1100
1101 for (Int_t j=1; j<=new->GetNbinsY(); ++j)
1102 {
1103 // npart sum to 1
1104 new->SetBinContent(i, j, new->GetBinContent(i, j) / sum);
1105 new->SetBinError(i, j, new->GetBinError(i, j) / sum);
1106 }
1107 }
1108
1109 new TCanvas;
1110 new->Draw("COLZ");
1111
1112 TH2* diff = (TH2*) new->Clone("diff");
1113 diff->Add(proj, -1);
1114
1115 new TCanvas;
1116 diff->Draw("COLZ");
1117 diff->SetMinimum(-0.05);
1118 diff->SetMaximum(0.05);
1119
1120 corr->Reset();
1121
1122 for (Int_t i=1; i<=new->GetNbinsX(); ++i)
1123 for (Int_t j=1; j<=new->GetNbinsY(); ++j)
1124 corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2, i, j, new->GetBinContent(i, j));
1125
1126 new TCanvas;
1127 corr->Project3D("zy")->Draw("COLZ");
1128
1129 TFile::Open("out.root", "RECREATE");
1130 mult->SaveHistograms();
1131
1132 TH1* proj1 = proj->ProjectionY("proj1", 36, 36);
1133 TH1* proj2 = new->ProjectionY("proj2", 36, 36);
1134 proj2->SetLineColor(2);
1135
1136 new TCanvas;
1137 proj1->Draw();
1138 proj2->Draw("SAME");
cfc19dd5 1139}