]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/selectors/multiplicity/correct.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGUD / selectors / multiplicity / correct.C
CommitLineData
0ab29cfa 1/* $Id$ */
2
3//
dca331bb 4// script to correct the multiplicity spectrum + helpers
0ab29cfa 5//
6
2701e5ae 7Bool_t is900GeV = 1;
eb9356d5 8Bool_t is2360TeV = 0;
9
10// 900 GeV, MC
2701e5ae 11const Int_t kBinLimits[] = { 42, 57, 60 };
12const Int_t kTrustLimits[] = { 26, 42, 54 };
eb9356d5 13
14// 2.36 TeV
15//const Int_t kBinLimits[] = { 43, 67, 83 };
16//const Int_t kTrustLimits[] = { 33, 57, 73 };
17
18// 7 TeV
2701e5ae 19//const Int_t kBinLimits[] = { 60, 120, 60 };
20//const Int_t kTrustLimits[] = { 26, 70, 54 };
eb9356d5 21
0b4bfd98 22void SetTPC()
0a173978 23{
24 gSystem->Load("libPWG0base");
0b4bfd98 25 AliMultiplicityCorrection::SetQualityRegions(kFALSE);
26}
0a173978 27
69b09e3b 28const char* GetMultLabel(Int_t etaR = -1, Bool_t trueM = kTRUE)
29{
30 if (etaR == -1)
31 etaR = etaRange;
32
33 TString tmpStr((trueM) ? "True " : "Measured ");
34
35 if (etaR == 4)
36 {
37 tmpStr += "multiplicity (full phase space)";
38 }
39 else
40 tmpStr += Form("multiplicity in |#eta| < %.1f", (etaR+1)* 0.5);
41 return Form("%s", tmpStr.Data());
42}
43
2440928d 44void draw(const char* fileName = "multiplicity.root", const char* folder = "Multiplicity")
0b4bfd98 45{
2440928d 46 loadlibs();
0a173978 47
0b4bfd98 48 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
0a173978 49
0b4bfd98 50 TFile::Open(fileName);
51 mult->LoadHistograms();
0a173978 52 mult->DrawHistograms();
447c325d 53
f3eb27f6 54 TH2* hist = (TH2*) gROOT->FindObject("fCorrelation2_zy");
447c325d 55 canvas = new TCanvas("c1", "c1", 600, 500);
f3eb27f6 56 canvas->SetTopMargin(0.05);
447c325d 57 hist->SetStats(kFALSE);
58 hist->Draw("COLZ");
f3eb27f6 59 hist->SetTitle(";true multiplicity in |#eta| < 1.5;measured multiplicity in |#eta| < 1.5");
447c325d 60 hist->GetYaxis()->SetTitleOffset(1.1);
f3eb27f6 61 gPad->SetRightMargin(0.12);
447c325d 62 gPad->SetLogz();
63
f3eb27f6 64 canvas->SaveAs("responsematrix.eps");
0ab29cfa 65}
66
0f67a57c 67void loadlibs()
9ca6be09 68{
0f67a57c 69 gSystem->Load("libTree");
70 gSystem->Load("libVMC");
71
72 gSystem->Load("libSTEERBase");
4070f709 73 gSystem->Load("libESD");
74 gSystem->Load("libAOD");
0f67a57c 75 gSystem->Load("libANALYSIS");
2440928d 76 gSystem->Load("libANALYSISalice");
9ca6be09 77 gSystem->Load("libPWG0base");
0f67a57c 78}
9ca6be09 79
eb9356d5 80void LoadAndInitialize(void* multVoid, void* esdVoid, void* multTriggerVoid, Int_t histID, Bool_t fullPhaseSpace, Int_t* geneLimits)
0f67a57c 81{
eb9356d5 82 AliMultiplicityCorrection* mult = (AliMultiplicityCorrection*) multVoid;
83 AliMultiplicityCorrection* esd = (AliMultiplicityCorrection*) esdVoid;
84 AliMultiplicityCorrection* multTrigger = (AliMultiplicityCorrection*) multTriggerVoid;
b3b260d1 85
eb9356d5 86 for (Int_t i=0; i<3; i++)
87 geneLimits[i] = kBinLimits[i];
88
89 // REBINNING
90 if (1)
447c325d 91 {
eb9356d5 92 // 900 GeV
93 if (is900GeV)
94 {
95 if (1)
96 {
97 Int_t bins05 = 31;
98 Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
99 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
100 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5,
101 100.5 };
102 }
103
104 if (0)
105 {
106 Int_t bins05 = 29;
107 Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
108 10.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
109 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 34.5, 38.5, 42.5, 50.5,
110 100.5 };
111 }
112
113 if (0)
114 {
115 Int_t bins05 = 25;
116 Double_t* newBinsEta05 = new Double_t[bins05+1];
117
118 //newBinsEta05[0] = -0.5;
119 //newBinsEta05[1] = 0.5;
120 //newBinsEta05[2] = 1.5;
121 //newBinsEta05[3] = 2.5;
122
123 for (Int_t i=0; i<bins05+1; i++)
124 newBinsEta05[i] = -0.5 + i*2;
125 //newBinsEta05[i] = 4.5 + (i-4)*2;
126 }
127
128 Int_t bins10 = 54;
129 Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
130 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
131 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
132 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
133 40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5,
134 60.5, 65.5, 70.5, 100.5 };
135
136 geneLimits[0] = bins05;
137 geneLimits[1] = bins10;
138 geneLimits[2] = bins10;
139 }
140
141 // 2.36 TeV
142 if (is2360TeV)
70fdd197 143 {
eb9356d5 144 Int_t bins05 = 36;
145 Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
146 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
147 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5,
148 40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
149 Int_t bins10 = 64;
150 Double_t newBinsEta10[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
151 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
152 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
153 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
154 40.5, 42.5, 44.5, 46.5, 48.5, 50.5, 52.5, 54.5, 56.5, 58.5,
155 60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5,
156 80.5, 85.5, 90.5, 100.5 };
157
158 geneLimits[0] = bins05;
159 geneLimits[1] = bins10;
160 geneLimits[2] = bins10;
70fdd197 161 }
eb9356d5 162
163 // 7 TeV
164 if (!is900GeV && !is2360TeV)
165 {
166 Int_t bins05 = 36;
167 Double_t newBinsEta05[] = { -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
168 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
169 20.5, 22.5, 24.5, 26.5, 28.5, 30.5, 32.5, 34.5, 36.5, 38.5,
170 40.5, 45.5, 50.5, 55.5, 60.5, 100.5 };
171
172 Int_t bins10 = 85;
173 Double_t newBinsEta10[] = {
174 -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
175 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5,
176 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
177 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
178 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5,
179 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5,
180 60.5, 62.5, 64.5, 66.5, 68.5, 70.5, 72.5, 74.5, 76.5, 78.5,
181 80.5, 82.5, 84.5, 86.5, 88.5, 90.5, 92.5, 94.5, 96.5, 98.5,
182 100.5, 105.5, 110.5, 115.5, 120.5 };
183
184 geneLimits[0] = bins05;
185 geneLimits[1] = bins10;
186 geneLimits[2] = bins10;
187 }
188
189 esd->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
190 mult->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
191 multTrigger->RebinGenerated(bins05, newBinsEta05, bins10, newBinsEta10, bins10, newBinsEta10);
192 }
70fdd197 193
eb9356d5 194 for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
195 {
b3b260d1 196 TH2F* hist = esd->GetMultiplicityESD(hID);
b3b260d1 197
7dcf959e 198 mult->SetMultiplicityESD(hID, hist);
eb9356d5 199 mult->SetTriggeredEvents(hID, esd->GetTriggeredEvents(hID));
200
201 // insert trigger efficiency in flat response matrix
202 for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
203 mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
204
205 mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
206 mult->FixTriggerEfficiencies(10);
207
b3b260d1 208 // small hack to get around charge conservation for full phase space ;-)
209 if (fullPhaseSpace)
210 {
211 TH1* corr = mult->GetCorrelation(histID + 4);
212
213 for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
214 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
215 {
216 corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
217 corr->SetBinError(i, j, corr->GetBinError(i-1, j));
218 }
219 }
eb9356d5 220 }
221}
222
223void correct(const char* fileNameMC = "multiplicityMC.root", const char* folder = "Multiplicity", const char* fileNameESD = "multiplicityESD.root", Bool_t chi2 = 1, Int_t histID = 1, Bool_t fullPhaseSpace = kFALSE, Float_t beta = -25, Int_t eventType = 2 /* AliMultiplicityCorrection::kTrVtx */, const char* targetFile = "unfolded.root", const char* folderESD = "Multiplicity", Bool_t calcBias = kFALSE)
224{
225 loadlibs();
226
227 Int_t geneLimits[] = { 0, 0, 0 };
b3b260d1 228
eb9356d5 229 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
230 AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
231 AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
232
233 LoadAndInitialize(mult, esd, multTrigger, histID, fullPhaseSpace, geneLimits);
234
235 //AliUnfolding::SetDebug(1);
236
237 for (Int_t hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
238 {
239 AliUnfolding::SetNbins(kBinLimits[hID], geneLimits[hID]);
240
241 TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
242
b3b260d1 243 if (chi2)
244 {
eb9356d5 245 if (is900GeV)
246 {
247 //Float_t regParamPol0[] = { 2, 2, 10 }; // TPCITS
248 Float_t regParamPol0[] = { 5, 15, 25 }; // SPD
249 Float_t regParamPol1[] = { 0.15, 0.25, 0.5 };
250 }
251 else if (is2360TeV)
252 {
253 Float_t regParamPol0[] = { 10, 12, 40 };
254 Float_t regParamPol1[] = { 0.25, 0.25, 2 };
255 }
256 else
257 {
258 Float_t regParamPol0[] = { 1, 25, 10 };
259 Float_t regParamPol1[] = { 0.15, 0.5, 0.5 };
260 AliUnfolding::SetSkipBinsBegin(3);
261 }
262
263 Int_t reg = AliUnfolding::kPol0;
264 if (beta > 0)
265 reg = AliUnfolding::kPol1;
266
267 Float_t regParam = TMath::Abs(beta);
268 if (histID == -1)
269 {
270 if (beta > 0)
271 regParam = regParamPol1[hID];
272 else
273 regParam = regParamPol0[hID];
274 }
275 AliUnfolding::SetChi2Regularization(reg, regParam);
276
277 //AliUnfolding::SetChi2Regularization(AliUnfolding::kLog, 1000000);
278 //AliUnfolding::SetChi2Regularization(AliUnfolding::kRatio, 10);
279 //TVirtualFitter::SetDefaultFitter("Minuit2");
280
70fdd197 281 //AliUnfolding::SetCreateOverflowBin(kFALSE);
b3b260d1 282 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0); //mult->SetCreateBigBin(kFALSE);
283 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125); mult->SetCreateBigBin(kFALSE);
70fdd197 284 // AliUnfolding::SetChi2Regularization(AliUnfolding::kEntropy, beta);
b3b260d1 285 //mult->SetRegularizationParameters(AliMultiplicityCorrection::kLog, 1e5);
286
eb9356d5 287 if (0)
288 {
289 // part for checking
290 TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
291 mcCompare->Sumw2();
292 mult->ApplyMinuitFit(hID, fullPhaseSpace, AliMultiplicityCorrection::kMB, 0, kTRUE, mcCompare);
293 mult->SetMultiplicityESDCorrected(hID, (TH1F*) mcCompare);
294 }
295 else
296 {
297 // Zero Bin
298 Int_t zeroBin = 0;
299 if (is900GeV) // from data
300 {
301 // background subtraction
302 Int_t background = 0;
303
304 //background = 1091 + 4398; // MB1 for run 104824 - 52 (SPD)
305 //background = 1087 + 4398; // MB1 for run 104824 - 52 (TPCITS)
306
2701e5ae 307 background = 417 + 1758; // MB1 for run 104867 - 92 (SPD)
eb9356d5 308 //background = 1162+422; // MB1 for run 104892 (TPCITS)
309 //background = 5830 + 1384; // MB1 for run 104824,25,45,52,67,92 (TPC runs) (TPCITS)
310
2701e5ae 311 //background = 20 + 0; // V0AND for run 104824 - 52
eb9356d5 312 //background = 10 + 0; // V0AND for run 104867 - 92
313
314 Printf("NOTE: Subtracting %d background events", background);
315 gSystem->Sleep(1000);
316
317 zeroBin = mult->GetTriggeredEvents(hID)->GetBinContent(1) - background - mult->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1);
318 }
319 else if (is2360TeV)
320 {
321 // from MC
322 Float_t fractionZeroBin = (multTrigger->GetTriggeredEvents(hID)->GetBinContent(1) - multTrigger->GetMultiplicityESD(hID)->Integral(0, 2, 1, 1)) / multTrigger->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
323 zeroBin = fractionZeroBin * mult->GetMultiplicityESD(hID)->Integral(1, 1, 1, mult->GetMultiplicityESD(hID)->GetNbinsY());
324
325 Printf("Zero bin from MC: Estimating %d events with trigger but without vertex", zeroBin);
326 gSystem->Sleep(1000);
327 }
328 else
329 {
330 AliUnfolding::SetSkip0BinInChi2(kTRUE);
331 }
b3b260d1 332
eb9356d5 333 //mult->SetVertexRange(3, 4);
334 mult->ApplyMinuitFit(hID, fullPhaseSpace, eventType, zeroBin, kFALSE, 0, calcBias);
335 }
b3b260d1 336 }
337 else
338 {
eb9356d5 339 // HACK
340 //mult->GetMultiplicityESD(hID)->SetBinContent(1, 1, 0);
341 //for (Int_t bin=1; bin<=mult->GetCorrelation(hID)->GetNbinsY(); bin++)
342 // mult->GetCorrelation(hID)->SetBinContent(1, bin, 1, 0);
343 AliUnfolding::SetChi2Regularization(AliUnfolding::kNone, 0);
344 mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, beta, 0, kTRUE);
70fdd197 345 //mult->ApplyBayesianMethod(hID, fullPhaseSpace, eventType, 1, 5, 0, kFALSE);
b3b260d1 346 }
dd701109 347 }
69b09e3b 348
349 Printf("Writing result to %s", targetFile);
350 TFile* file = TFile::Open(targetFile, "RECREATE");
eb9356d5 351 mult->SaveHistograms("Multiplicity");
144ff489 352 file->Write();
353 file->Close();
0a0f4adb 354
eb9356d5 355 if (histID == -1)
356 return;
357
b3b260d1 358 for (hID = ((histID == -1) ? 0 : histID); hID <= ((histID == -1) ? 2 : histID); hID++)
359 {
360 TH2F* hist2 = esd->GetMultiplicityMC(hID, eventType);
361 TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
eb9356d5 362 mult->DrawComparison(Form("%s_%d_%f", (chi2) ? "MinuitChi2" : "Bayesian", hID, beta), hID, fullPhaseSpace, kTRUE, mcCompare, kFALSE, eventType);
363
364 Printf("<n> = %f", mult->GetMultiplicityESDCorrected(hID)->GetMean());
b3b260d1 365 }
69b09e3b 366}
367
eb9356d5 368void correctAll(Int_t eventType)
369{
370 const char* suffix = "";
371 switch (eventType)
372 {
373 case 0: suffix = "trvtx"; break;
374 case 1: suffix = "mb"; break;
375 case 2: suffix = "inel"; break;
376 case 3: suffix = "nsd"; break;
377 }
378
379 correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0, -1, eventType, TString(Form("chi2a_%s.root", suffix)));
380 correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, -1, 0, 1, eventType, TString(Form("chi2b_%s.root", suffix)));
381 correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 0, -1, 0, 40, eventType, TString(Form("bayes_%s.root", suffix)));
382
383 if (eventType == 3)
384 drawAll(1);
385 else if (eventType == 2)
386 drawAllINEL();
387}
388
389void drawAll(Bool_t showUA5 = kFALSE)
390{
391 const char* files[] = { "chi2a_nsd.root", "chi2b_nsd.root", "bayes_nsd.root" };
392 drawAll(files, showUA5);
393}
394
395void drawAllINEL()
396{
397 const char* files[] = { "chi2a_inel.root", "chi2b_inel.root", "bayes_inel.root" };
398 drawAll(files);
399}
400
401void drawAllMB()
402{
403 const char* files[] = { "chi2a_mb.root", "chi2b_mb.root", "bayes_mb.root" };
404 drawAll(files);
405}
406
407void drawAll(const char** files, Bool_t showUA5 = kFALSE, Bool_t normalize = kFALSE)
70fdd197 408{
409 loadlibs();
410
eb9356d5 411 Int_t colors[] = { 1, 2, 4 };
412
413 c = new TCanvas(Form("c%d", gRandom->Uniform(100)), "c", 1800, 600);
414 c->Divide(3, 1);
415
416 l = new TLegend(0.6, 0.6, 0.99, 0.99);
417 l->SetFillColor(0);
418
419 TH1* hist0 = 0;
420
421 for (Int_t hID=0; hID<3; hID++)
422 {
423 c->cd(hID+1)->SetLogy();
424 gPad->SetGridx();
425 gPad->SetGridy();
426 for (Int_t i=0; i<3; i++)
427 {
428 TFile::Open(files[i]);
429
430 hist = (TH1*) gFile->Get(Form("Multiplicity/fMultiplicityESDCorrected%d", hID));
70fdd197 431
eb9356d5 432 Float_t average0 = hist->GetMean();
433 hist2 = (TH1*) hist->Clone("temp");
434 hist2->SetBinContent(1, 0);
435 Float_t average1 = hist2->GetMean();
436 Printf("%d: <N> = %.2f <N>/(eta) = %.2f | without 0: <N> = %.2f <N>/(eta) = %.2f", hID, average0, average0 / ((hID+1) - 0.4 * (hID / 2)), average1, average1 / ((hID+1) - 0.4 * (hID / 2)));
437
438 hist->SetLineColor(colors[i]);
439
440 hist->SetStats(0);
441 hist->GetXaxis()->SetRangeUser(0, kBinLimits[hID]);
442
443 Float_t min = 0.1;
444 Float_t max = hist->GetMaximum() * 1.5;
445
446 if (normalize)
447 min = 1e-6;
448
449 hist->GetYaxis()->SetRangeUser(min, max);
450 hist->SetTitle(";unfolded multiplicity;events");
451
452 if (hID == 0)
453 {
454 l->AddEntry(hist, files[i]);
455 }
456
457 if (hist->Integral() <= 0)
458 continue;
459
460 if (normalize)
461 hist->Scale(1.0 / hist->Integral());
462
463 AliPWG0Helper::NormalizeToBinWidth(hist);
464
465 hist->DrawCopy((i == 0) ? "" : "SAME");
466
467 if (!hist0)
468 hist0 = (TH1*) hist->Clone("hist0");
469 }
470
471 if (hist0)
472 {
473 line = new TLine(kTrustLimits[hID], hist0->GetMinimum(), kTrustLimits[hID], hist0->GetMaximum());
474 line->SetLineWidth(3);
475 line->Draw();
476 }
477 }
70fdd197 478
eb9356d5 479 c->cd(1);
480 l->Draw();
70fdd197 481
eb9356d5 482 if (showUA5)
483 {
484 TGraphErrors *gre = new TGraphErrors(23);
485 gre->SetName("Graph");
486 gre->SetTitle("UA5");
487 gre->SetFillColor(1);
488 gre->SetMarkerStyle(24);
489 gre->SetPoint(0,0,0.15);
490 gre->SetPointError(0,0.5,0.01);
491 gre->SetPoint(1,1,0.171);
492 gre->SetPointError(1,0.5,0.007);
493 gre->SetPoint(2,2,0.153);
494 gre->SetPointError(2,0.5,0.007);
495 gre->SetPoint(3,3,0.124);
496 gre->SetPointError(3,0.5,0.006);
497 gre->SetPoint(4,4,0.099);
498 gre->SetPointError(4,0.5,0.005);
499 gre->SetPoint(5,5,0.076);
500 gre->SetPointError(5,0.5,0.005);
501 gre->SetPoint(6,6,0.057);
502 gre->SetPointError(6,0.5,0.004);
503 gre->SetPoint(7,7,0.043);
504 gre->SetPointError(7,0.5,0.004);
505 gre->SetPoint(8,8,0.032);
506 gre->SetPointError(8,0.5,0.003);
507 gre->SetPoint(9,9,0.024);
508 gre->SetPointError(9,0.5,0.003);
509 gre->SetPoint(10,10,0.018);
510 gre->SetPointError(10,0.5,0.002);
511 gre->SetPoint(11,11,0.013);
512 gre->SetPointError(11,0.5,0.002);
513 gre->SetPoint(12,12,0.01);
514 gre->SetPointError(12,0.5,0.002);
515 gre->SetPoint(13,13,0.007);
516 gre->SetPointError(13,0.5,0.001);
517 gre->SetPoint(14,14,0.005);
518 gre->SetPointError(14,0.5,0.001);
519 gre->SetPoint(15,15,0.004);
520 gre->SetPointError(15,0.5,0.001);
521 gre->SetPoint(16,16,0.0027);
522 gre->SetPointError(16,0.5,0.0009);
523 gre->SetPoint(17,17,0.0021);
524 gre->SetPointError(17,0.5,0.0008);
525 gre->SetPoint(18,18,0.0015);
526 gre->SetPointError(18,0.5,0.0007);
527 gre->SetPoint(19,19,0.0011);
528 gre->SetPointError(19,0.5,0.0006);
529 gre->SetPoint(20,20,0.0008);
530 gre->SetPointError(20,0.5,0.0005);
531 gre->SetPoint(21,22,0.00065);
532 gre->SetPointError(21,1,0.0003);
533 gre->SetPoint(22,27.5,0.00013);
534 gre->SetPointError(22,3.5,7.1e-05);
535
536 for (Int_t i=0; i<gre->GetN(); i++)
537 {
538 gre->GetY()[i] *= hist0->Integral();
539 gre->GetEY()[i] *= hist0->Integral();
540 }
541
542 gre->Draw("P");
543
544 ratio = (TGraphErrors*) gre->Clone("ratio");
545
546 for (Int_t i=0; i<gre->GetN(); i++)
547 {
548 //Printf("%f %d", hist0->GetBinContent(hist0->FindBin(ratio->GetX()[i])), hist0->FindBin(ratio->GetX()[i]));
549 Int_t bin = hist0->FindBin(gre->GetX()[i]);
550 if (hist0->GetBinContent(bin) > 0)
551 {
552 ratio->GetEY()[i] = TMath::Sqrt((ratio->GetEY()[i] / ratio->GetY()[i])**2 + (hist0->GetBinError(bin) / hist0->GetBinContent(bin))**2);
553 ratio->GetY()[i] /= hist0->GetBinContent(bin);
554 ratio->GetEY()[i] *= ratio->GetY()[i];
555 }
556 }
557 new TCanvas;
558 gPad->SetGridx();
559 gPad->SetGridy();
560 //ratio->Print();
561 ratio->Draw("AP");
562 ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
563 }
564
565 c->SaveAs("draw_all.png");
566}
567
568TH1* GetChi2Bias(Float_t alpha = -5)
569{
570 loadlibs();
571
572 AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kMB;
573 Int_t hID = 1;
70fdd197 574
eb9356d5 575 correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "nobias.root", "Multiplicity", kFALSE);
576
577 correct("multiplicityMC.root", "Multiplicity", "multiplicityESD.root", 1, hID, 0, alpha, eventType, "withbias.root", "Multiplicity", kTRUE);
578
70fdd197 579 // without bias calculation
eb9356d5 580 mult = AliMultiplicityCorrection::Open("nobias.root");
70fdd197 581 baseold = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("baseold");
582
583 // with bias calculation
eb9356d5 584 mult = AliMultiplicityCorrection::Open("withbias.root");
70fdd197 585 base = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("base");
586
587 // relative error plots
588 ratio1 = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("hist1");
589 ratio1->SetStats(0);
590 ratio1->SetTitle(";unfolded multiplicity;relative error");
591 ratio1->Reset();
592 ratio2 = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("hist1");
593 ratio2->Reset();
594
eb9356d5 595 for (Int_t t = 0; t<baseold->GetNbinsX(); t++)
70fdd197 596 {
597 Printf("Bin %d; Content: %f; Chi2 Error: %f; Bias: %f; In percent: %.2f %.2f", t+1, base->GetBinContent(t+1), baseold->GetBinError(t+1), base->GetBinError(t+1), (base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1, (base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
598 if (base->GetBinContent(t+1) <= 0)
599 continue;
600 ratio1->SetBinContent(t+1, baseold->GetBinError(t+1) / base->GetBinContent(t+1));
601 ratio2->SetBinContent(t+1, base->GetBinError(t+1) / base->GetBinContent(t+1));
602 }
603
604 //new TCanvas; base->Draw(); gPad->SetLogy();
605
eb9356d5 606 c = new TCanvas;
70fdd197 607 ratio1->GetYaxis()->SetRangeUser(0, 1);
eb9356d5 608 ratio1->GetXaxis()->SetRangeUser(0, 50);
70fdd197 609 ratio1->Draw();
610 ratio2->SetLineColor(2);
611 ratio2->Draw("SAME");
612
613 return base;
614}
615
616void CheckBayesianBias()
617{
618 loadlibs();
619
620 const char* fileNameMC = "multiplicityMC.root";
621 const char* folder = "Multiplicity";
622 const char* fileNameESD = "multiplicityESD.root";
623 const char* folderESD = "Multiplicity";
624
625 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
626 AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
627
628 const Int_t kMaxBins = 35;
629 AliUnfolding::SetNbins(kMaxBins, kMaxBins);
630 //AliUnfolding::SetDebug(1);
631
632 Int_t hID = 0;
633
634 TH2F* hist = esd->GetMultiplicityESD(hID);
635 mult->SetMultiplicityESD(hID, hist);
636
637 // without bias calculation
638 mult->ApplyBayesianMethod(hID, kFALSE, 0, 1, 10, 0, 1);
639 baseold = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("baseold");
640
641 // with bias calculation
642 mult->ApplyBayesianMethod(hID, kFALSE, 0, 1, 10, 0, 2);
643 base = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("base");
644
645 TH1* hist1 = new TH1F("hist1", "", 100, 0, 20);
646 TH1* hist2 = new TH1F("hist2", "", 100, 0, 20);
647
648 for (Int_t t = 0; t<kMaxBins; t++)
649 {
650 Printf("Bin %d; Content: %f; Randomization Error: %f; Bias: %f; In percent: %.2f %.2f", t+1, base->GetBinContent(t+1), baseold->GetBinError(t+1), base->GetBinError(t+1), (base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1, (base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
651
652 hist1->Fill((base->GetBinContent(t+1) > 0) ? 100.0 * baseold->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
653 hist2->Fill((base->GetBinContent(t+1) > 0) ? 100.0 * base->GetBinError(t+1) / base->GetBinContent(t+1) : -1);
654 }
655
656 new TCanvas;
657 hist1->Draw();
658 hist2->SetLineColor(2);
659 hist2->Draw("SAME");
660}
661
eb9356d5 662void DrawUnfoldingLimit(Int_t hID, Float_t min, Float_t max)
663{
664 line = new TLine(kTrustLimits[hID], min, kTrustLimits[hID], max);
665 line->SetLineWidth(2);
666 line->Draw();
667}
668
669void DataScan(Int_t hID, Bool_t redo = kTRUE)
70fdd197 670{
671 // makes a set of unfolded spectra and compares
672 // don't forget FindUnfoldedLimit in plots.C
673
674 loadlibs();
675
676 // files...
eb9356d5 677 Bool_t mc = 1;
70fdd197 678 const char* fileNameMC = "multiplicityMC.root";
679 const char* folder = "Multiplicity";
680 const char* fileNameESD = "multiplicityESD.root";
681 const char* folderESD = "Multiplicity";
eb9356d5 682 const Int_t kMaxBins = kBinLimits[hID];
683
70fdd197 684 AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx;
685 Bool_t evaluteBias = kFALSE;
686
687 Int_t referenceCase = 2; // all results are shown as ratio to this case
688
689 // chi2 range
690 AliUnfolding::RegularizationType regBegin = AliUnfolding::kPol0;
eb9356d5 691 AliUnfolding::RegularizationType regEnd = AliUnfolding::kPol1;
692 Float_t regParamBegin[] = { 0, 1, 0.2, 3 };
693 Float_t regParamEnd[] = { 0, 11, 0.5, 31 };
694 Float_t regParamStep[] = { 0, 2, 2, TMath::Sqrt(10) };
70fdd197 695
696 // bayesian range
eb9356d5 697 Int_t iterBegin = 40;
698 Int_t iterEnd = 41;
699 Int_t iterStep = 20;
70fdd197 700
701 TList labels;
702
703 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open(fileNameMC, folder);
704 AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open(fileNameESD, folderESD);
705
706 mult->SetMultiplicityESD(hID, esd->GetMultiplicityESD(hID));
707
eb9356d5 708 // insert trigger efficiency in flat response matrix
709 AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
710 for (Int_t evType = AliMultiplicityCorrection::kTrVtx; evType <= AliMultiplicityCorrection::kNSD; evType++)
711 mult->SetMultiplicityMC(hID, evType, multTrigger->GetMultiplicityMC(hID, evType));
712
713 mult->SetNoVertexEvents(hID, multTrigger->GetNoVertexEvents(hID));
714 mult->FixTriggerEfficiencies(10);
715
716 Float_t fraction0Generated = multTrigger->GetFraction0Generated(hID);
717
70fdd197 718 if (mc)
719 mult->SetMultiplicityMC(hID, eventType, esd->GetMultiplicityMC(hID, eventType));
720
721 AliUnfolding::SetNbins(kMaxBins, kMaxBins);
722 //AliUnfolding::SetDebug(1);
723
724 Int_t count = -1;
725
726 for (AliUnfolding::RegularizationType reg = regBegin; reg <= regEnd; reg++)
727 {
728 for (Float_t regParam = regParamBegin[reg]; regParam < regParamEnd[reg]; regParam *= regParamStep[reg])
729 {
730 count++;
731 labels.Add(new TObjString(Form("#chi^{2} Reg %d #beta = %.2f", (Int_t) reg, regParam)));
732
733 if (!redo)
734 continue;
735
736 AliUnfolding::SetChi2Regularization(reg, regParam);
737
eb9356d5 738 mult->ApplyMinuitFit(hID, kFALSE, eventType, fraction0Generated, kFALSE, 0, evaluteBias);
70fdd197 739
740 file = TFile::Open(Form("datascan_%d.root", count), "RECREATE");
741 mult->SaveHistograms();
742 file->Close();
743 }
744 }
745
746 for (Int_t iter = iterBegin; iter <= iterEnd; iter += iterStep)
747 {
748 count++;
749 labels.Add(new TObjString(Form("Bayesian Iter = %d", iter)));
750
751 if (!redo)
752 continue;
753
754 mult->ApplyBayesianMethod(hID, kFALSE, eventType, 1, iter, 0, kTRUE);
755
756 file = TFile::Open(Form("datascan_%d.root", count), "RECREATE");
757 mult->SaveHistograms();
758 file->Close();
759 }
760
761 // 1) all distributions
762 // 2) ratio to MC
763 // 3) ratio to ref point
764 // 4) relative error
765 // 5) residuals
766 c = new TCanvas("c", "c", 1200, 800);
767 c->Divide(3, 2);
768
769 c->cd(1)->SetLogy();
770
771 // get reference point
772 mult = AliMultiplicityCorrection::Open(Form("datascan_%d.root", referenceCase));
773 if (!mult)
774 return;
775 refPoint = (TH1*) mult->GetMultiplicityESDCorrected(hID)->Clone("refPoint");
776
777 legend = new TLegend(0.5, 0.5, 0.99, 0.99);
778 legend->SetFillColor(0);
779
780 TH1* residualSum = new TH1F("residualSum", ";;residual squared sum", count + 1, 0.5, count + 1.5);
781
eb9356d5 782 Int_t allowedCount = count;
70fdd197 783 count = 0;
784 while (1)
785 {
786 mult = AliMultiplicityCorrection::Open(Form("datascan_%d.root", count++));
787 if (!mult)
788 break;
eb9356d5 789 if (count > allowedCount+1)
790 break;
70fdd197 791
792 hist = (TH1*) mult->GetMultiplicityESDCorrected(hID);
793 hist->SetLineColor((count-1) % 4 + 1);
794 hist->SetLineStyle((count-1) / 4 + 1);
795 hist->GetXaxis()->SetRangeUser(0, kMaxBins);
eb9356d5 796 hist->GetYaxis()->SetRangeUser(0.1, hist->GetMaximum() * 1.5);
70fdd197 797 hist->SetStats(kFALSE);
798 hist->SetTitle("");
799
800 legend->AddEntry(hist->Clone(), labels.At(count-1)->GetName());
801
802 c->cd(1);
803 hist->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
804
eb9356d5 805 DrawUnfoldingLimit(hID, 0.1, hist->GetMaximum());
806
70fdd197 807 if (mc)
808 {
809 TH2* mcHist2d = mult->GetMultiplicityMC(hID, eventType);
810 mcHist = mcHist2d->ProjectionY("mcmchist", 1, mcHist2d->GetNbinsX());
811
812 c->cd(1);
813 mcHist->SetMarkerStyle(24);
814 mcHist->Draw("P SAME");
eb9356d5 815
70fdd197 816 c->cd(2);
817 // calculate ratio using only the error on the mc bin
818 ratio = (TH1*) hist->Clone("ratio");
819 for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
820 {
821 if (mcHist->GetBinContent(bin) <= 0)
822 continue;
823 ratio->SetBinContent(bin, hist->GetBinContent(bin) / mcHist->GetBinContent(bin));
824 ratio->SetBinError(bin, mcHist->GetBinError(bin) / mcHist->GetBinContent(bin) * ratio->GetBinContent(bin));
825 }
826
827 ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
828 ratio->GetYaxis()->SetTitle("Ratio unfolded / MC");
829 ratio->DrawCopy((count == 1) ? "HIST" : "HISTSAME");
eb9356d5 830
831 DrawUnfoldingLimit(hID, 0.5, 1.5);
70fdd197 832 }
833
834 c->cd(3);
835 // calculate ratio using no errors for now
836 ratio = (TH1*) hist->Clone("ratio");
837 for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
838 {
839 if (refPoint->GetBinContent(bin) <= 0)
840 continue;
841 ratio->SetBinContent(bin, hist->GetBinContent(bin) / refPoint->GetBinContent(bin));
842 ratio->SetBinError(bin, 0);
843 }
844
845 ratio->GetYaxis()->SetRangeUser(0.5, 1.5);
846 ratio->GetYaxis()->SetTitle("Ratio unfolded / unfolded reference case");
847 ratio->DrawCopy((count == 1) ? "" : "SAME");
848
eb9356d5 849 DrawUnfoldingLimit(hID, 0.5, 1.5);
850
70fdd197 851 c->cd(4);
852 ratio = (TH1*) hist->Clone("ratio");
853 for (Int_t bin=1; bin<=ratio->GetNbinsX(); bin++)
854 {
855 if (ratio->GetBinContent(bin) <= 0)
856 continue;
857 ratio->SetBinContent(bin, ratio->GetBinError(bin) / ratio->GetBinContent(bin));
858 ratio->SetBinError(bin, 0);
859 }
860 ratio->GetYaxis()->SetRangeUser(0, 1);
861 ratio->GetYaxis()->SetTitle("Relative error");
862 ratio->DrawCopy((count == 1) ? "" : "SAME");
863
eb9356d5 864 DrawUnfoldingLimit(hID, 0, 1);
865
70fdd197 866 c->cd(5);
867 Float_t sum;
868 residuals = mult->GetResiduals(hID, AliMultiplicityCorrection::kTrVtx, sum);
869 residuals->SetStats(0);
870 residuals->GetXaxis()->SetRangeUser(0, kMaxBins);
871 residuals->SetStats(kFALSE);
872 residuals->SetLineColor((count-1) % 4 + 1);
873 residuals->SetLineStyle((count-1) / 4 + 1);
874 residuals->GetYaxis()->SetRangeUser(-5, 5);
875 residuals->DrawCopy((count == 1) ? "" : "SAME");
876
877 residualSum->Fill(count, sum);
878 residualSum->GetXaxis()->SetBinLabel(count, labels.At(count-1)->GetName());
879 }
880
881 c->cd(6);
882 residualSum->Draw();
883 legend->Draw();
884}
885
69b09e3b 886void DrawBayesianIterations(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 1, Int_t eventType = 0 /* AliMultiplicityCorrection::kTrVtx */)
887{
888 loadlibs();
889
890 const char* folder = "Multiplicity";
891
892 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection(folder, folder);
893 TFile::Open(fileNameMC);
894 mult->LoadHistograms();
895
896 AliMultiplicityCorrection* esd = new AliMultiplicityCorrection(folder, folder);
897 TFile::Open(fileNameESD);
898 esd->LoadHistograms();
899
900 TH2F* hist = esd->GetMultiplicityESD(histID);
901 TH2F* hist2 = esd->GetMultiplicityMC(histID, eventType);
902
903 mult->SetMultiplicityESD(histID, hist);
904 mult->SetMultiplicityMC(histID, eventType, hist2);
905
906 TH1* mcCompare = hist2->ProjectionY("mcmchist", 1, hist2->GetNbinsX());
907 //mcCompare->Scale(1.0 / mcCompare->Integral());
908
909 TH1* esdHist = (TH1*) hist->ProjectionY("myesd", 1, 1)->Clone("myesd");
910 //esdHist->Scale(1.0 / esdHist->Integral());
911
912 c = new TCanvas("c", "c", 1200, 600);
913 c->Divide(2, 1);
914
915 c->cd(1);
916 gPad->SetLeftMargin(0.12);
917 gPad->SetTopMargin(0.05);
918 gPad->SetRightMargin(0.05);
919 gPad->SetLogy();
920 gPad->SetGridx();
921 gPad->SetGridy();
922
923 mcCompare->GetXaxis()->SetRangeUser(0, 80);
924 mcCompare->SetStats(0);
925 mcCompare->SetFillColor(5);
926 mcCompare->SetLineColor(5);
927 mcCompare->SetTitle(Form(";%s;Entries", GetMultLabel(1)));
928 mcCompare->GetYaxis()->SetRangeUser(2, esdHist->GetMaximum() * 2);
929 mcCompare->GetYaxis()->SetTitleOffset(1.3);
930 mcCompare->Draw("HIST");
931 esdHist->SetMarkerStyle(5);
932 esdHist->Draw("P HIST SAME");
933
934 c->cd(2);
935 gPad->SetTopMargin(0.05);
936 gPad->SetRightMargin(0.05);
937 gPad->SetGridx();
938 gPad->SetGridy();
939
940 trueMeasuredRatio = (TH1*) mcCompare->Clone("trueMeasuredRatio");
941 trueMeasuredRatio->Divide(esdHist);
942 trueMeasuredRatio->SetStats(0);
943 trueMeasuredRatio->SetTitle(Form(";%s;Ratio", GetMultLabel(1)));
944 trueMeasuredRatio->GetYaxis()->SetTitleOffset(1.3);
945 trueMeasuredRatio->GetYaxis()->SetRangeUser(0, 2);
946 // ROOT displays all values > 2 at 2 which looks weird
947 for (Int_t i=1; i<=trueMeasuredRatio->GetNbinsX(); i++)
948 if (trueMeasuredRatio->GetBinContent(i) > 2)
949 trueMeasuredRatio->SetBinContent(i, 0);
950 trueMeasuredRatio->SetMarkerStyle(5);
951 trueMeasuredRatio->Draw("P");
952
953 Int_t colors[] = {1, 2, 4, 6};
954
955 legend = new TLegend(0.15, 0.13, 0.5, 0.5);
956 legend->SetFillColor(0);
957 legend->SetTextSize(0.04);
958
959 legend->AddEntry(mcCompare, "True", "F");
960 legend->AddEntry(esdHist, "Measured", "P");
961
962 legend2 = new TLegend(0.15, 0.13, 0.5, 0.4);
963 legend2->SetFillColor(0);
964 legend2->SetTextSize(0.04);
965
966 legend2->AddEntry(trueMeasuredRatio, "Measured", "P");
967
968 Int_t iters[] = {1, 3, 10, -1};
969 for (Int_t i = 0; i<4; i++)
970 {
971 Int_t iter = iters[i];
972 mult->ApplyBayesianMethod(histID, kFALSE, eventType, 1, iter, 0, 0);
973 corr = mult->GetMultiplicityESDCorrected(histID);
974 corr->Scale(1.0 / corr->Integral());
975 corr->Scale(mcCompare->Integral());
976 corr->GetXaxis()->SetRangeUser(0, 80);
977 //corr->SetMarkerStyle(20+iter);
978 corr->SetLineColor(colors[i]);
979 corr->SetLineStyle(i+1);
980 corr->SetLineWidth(2);
981
982 c->cd(1);
983 corr->DrawCopy("SAME HIST");
984
985 c->cd(2);
986 clone = (TH1*) corr->Clone("clone");
987 clone->Divide(mcCompare, corr);
988 clone->GetYaxis()->SetRangeUser(0, 2);
989 clone->DrawCopy("SAME HIST");
990
991 legend->AddEntry((TH1*) corr->Clone(), (iter > -1) ? Form("%d iteration%s", iter, (iter == 1) ? "" : "s") : "Convergence", "L");
992 legend2->AddEntry((TH1*) corr->Clone(), (iter > -1) ? Form("%d iteration%s", iter, (iter == 1) ? "" : "s") : "Convergence", "L");
993 }
994
995 c->cd(1);
996 legend->Draw();
997
998 c->cd(2);
999 legend2->Draw();
1000
1001 c->SaveAs("bayesian_iterations.eps");
0a0f4adb 1002}
1003
f3eb27f6 1004void CompareChi2Bayesian(Int_t histID = 2, const char* chi2File = "chi2.root", const char* bayesianFile = "bayesian.root", const char* label1 = "Chi2", const char* label2 = "Bayesian", const char* mcFile = 0, Float_t simpleCorrect = 0)
0a0f4adb 1005{
1006 const char* folder = "Multiplicity";
1007
1008 loadlibs();
1009
1010 AliMultiplicityCorrection* chi2 = new AliMultiplicityCorrection(folder, folder);
1011 TFile::Open(chi2File);
1012 chi2->LoadHistograms();
1013
1014 AliMultiplicityCorrection* bayesian = new AliMultiplicityCorrection(folder, folder);
1015 TFile::Open(bayesianFile);
1016 bayesian->LoadHistograms();
1017
f3eb27f6 1018 histRAW = chi2->GetMultiplicityESD(histID)->ProjectionY("raw", 1, chi2->GetMultiplicityESD(histID)->GetNbinsX());
1019 histRAW->Scale(1.0 / histRAW->Integral());
1020
0a0f4adb 1021 histC = chi2->GetMultiplicityESDCorrected(histID);
1022 histB = bayesian->GetMultiplicityESDCorrected(histID);
1023
1024 c = new TCanvas("CompareChi2Bayesian", "CompareChi2Bayesian", 800, 600);
1025 c->SetRightMargin(0.05);
1026 c->SetTopMargin(0.05);
1027 c->SetLogy();
f3eb27f6 1028 c->SetGridx();
1029 c->SetGridy();
0a0f4adb 1030
1031 histC->SetTitle(";N;P(N)");
1032 histC->SetStats(kFALSE);
1033 histC->GetXaxis()->SetRangeUser(0, 100);
1034
1035 histC->SetLineColor(1);
1036 histB->SetLineColor(2);
f3eb27f6 1037 histRAW->SetLineColor(3);
0a0f4adb 1038
51f6de65 1039 histC->DrawCopy("HISTE");
1040 histB->DrawCopy("HISTE SAME");
f3eb27f6 1041 histRAW->DrawCopy("SAME");
0a0f4adb 1042
1043 legend = new TLegend(0.2, 0.2, 0.4, 0.4);
1044 legend->SetFillColor(0);
1045
f3eb27f6 1046 legend->AddEntry(histC, label1);
1047 legend->AddEntry(histB, label2);
1048 legend->AddEntry(histRAW, "raw ESD");
0a0f4adb 1049
51f6de65 1050 histGraph = 0;
f3eb27f6 1051 if (simpleCorrect > 0)
1052 {
1053 graph = new TGraph;
1054 graph->SetMarkerStyle(25);
1055 graph->SetFillColor(0);
1056 for (Int_t bin=1; bin<=histRAW->GetNbinsX(); bin++)
1057 graph->SetPoint(graph->GetN(), histRAW->GetXaxis()->GetBinCenter(bin) * simpleCorrect, histRAW->GetBinContent(bin));
1058
1059 graph->Draw("PSAME");
1060 legend->AddEntry(graph, "weighting");
1061
1062 // now create histogram from graph and normalize
1063 histGraph = (TH1*) histRAW->Clone();
1064 histGraph->Reset();
1065 for (Int_t bin=1; bin<=histGraph->GetNbinsX(); bin++)
1066 {
1067 Int_t j=1;
1068 for (j=1; j<graph->GetN(); j++)
1069 if (graph->GetX()[j] > histGraph->GetXaxis()->GetBinCenter(bin))
1070 break;
1071 if (j == graph->GetN())
1072 continue;
1073 if (histGraph->GetXaxis()->GetBinCenter(bin) - graph->GetX()[j] < graph->GetX()[j-1] - histGraph->GetXaxis()->GetBinCenter(bin))
1074 j--;
1075 histGraph->SetBinContent(bin, graph->GetY()[j]);
1076 }
1077
1078 Printf("Integral = %f", histGraph->Integral());
1079 histGraph->Scale(1.0 / histGraph->Integral());
1080
1081 histGraph->SetLineColor(6);
1082 histGraph->DrawCopy("SAME");
1083 legend->AddEntry(histGraph, "weighting normalized");
1084 }
1085
1086 if (mcFile)
0a0f4adb 1087 {
1088 AliMultiplicityCorrection* mc = new AliMultiplicityCorrection(folder, folder);
1089 TFile::Open(mcFile);
1090 mc->LoadHistograms();
1091
1092 histMC = mc->GetMultiplicityVtx(histID)->ProjectionY("mc", 1, mc->GetMultiplicityVtx(histID)->GetNbinsX());
f3eb27f6 1093 histMC->Sumw2();
0a0f4adb 1094 histMC->Scale(1.0 / histMC->Integral());
1095
51f6de65 1096 histMC->Draw("HISTE SAME");
0a0f4adb 1097 histMC->SetLineColor(4);
1098 legend->AddEntry(histMC, "MC");
1099 }
1100
1101 legend->Draw();
1102
f3eb27f6 1103 c->SaveAs(Form("%s.png", c->GetName()));
1104 c->SaveAs(Form("%s.eps", c->GetName()));
1105
1106 if (!mcFile)
1107 return;
1108
1109 // build ratios
1110
1111 c = new TCanvas("CompareChi2BayesianRatio", "CompareChi2BayesianRatio", 800, 600);
1112 c->SetRightMargin(0.05);
1113 c->SetTopMargin(0.05);
1114 c->SetGridx();
1115 c->SetGridy();
1116
1117 for (Int_t bin=1; bin<=histC->GetNbinsX(); bin++)
1118 {
1119 if (histMC->GetBinContent(bin) > 0)
1120 {
1121 histC->SetBinContent(bin, histC->GetBinContent(bin) / histMC->GetBinContent(bin));
1122 histB->SetBinContent(bin, histB->GetBinContent(bin) / histMC->GetBinContent(bin));
51f6de65 1123
1124 /*
1125 if (histGraph)
1126 {
1127 histGraph->SetBinContent(bin, histGraph->GetBinContent(bin) / histMC->GetBinContent(bin));
1128 histGraph->SetBinError(bin, 0);
1129 }
1130 */
f3eb27f6 1131
1132 // TODO errors?
1133 histC->SetBinError(bin, 0);
1134 histB->SetBinError(bin, 0);
f3eb27f6 1135 }
1136 }
1137
1138 histC->GetYaxis()->SetRangeUser(0.5, 2);
1139 histC->GetYaxis()->SetTitle("Unfolded / MC");
1140
1141 histC->Draw("HIST");
1142 histB->Draw("HIST SAME");
51f6de65 1143 /*
1144 if (histGraph)
1145 histGraph->Draw("HIST SAME");
1146 */
f3eb27f6 1147
1148 /*
1149 if (simpleCorrect > 0)
1150 {
1151 graph2 = new TGraph;
1152 graph2->SetMarkerStyle(25);
1153 graph2->SetFillColor(0);
1154 for (Int_t i=0; i<graph->GetN(); i++)
1155 {
1156 Float_t mcValue = histMC->GetBinContent(histMC->FindBin(graph->GetX()[i]));
1157 Float_t mcError = histMC->GetBinError(histMC->FindBin(graph->GetX()[i]));
1158 if (mcValue > 0)
1159 graph2->SetPoint(graph2->GetN(), graph->GetX()[i], graph->GetY()[i] / mcValue);
1160 }
1161 graph2->Draw("PSAME");
1162 }
1163 */
1164
1165 c->SaveAs(Form("%s.png", c->GetName()));
1166 c->SaveAs(Form("%s.eps", c->GetName()));
1167}
1168
1169
1170void CompareMC(Int_t histID, Int_t eventType, const char* file1, const char* file2, const char* label1, const char* label2)
1171{
1172 const char* folder = "Multiplicity";
1173
1174 loadlibs();
1175
1176 AliMultiplicityCorrection* mc1 = new AliMultiplicityCorrection(folder, folder);
1177 TFile::Open(file1);
1178 mc1->LoadHistograms();
1179 histMC1 = mc1->GetMultiplicityMC(histID, eventType)->ProjectionY("mc1", 1, mc1->GetMultiplicityMC(histID, eventType)->GetNbinsX());
1180 histMC1->Scale(1.0 / histMC1->Integral());
1181
1182 AliMultiplicityCorrection* mc2 = new AliMultiplicityCorrection(folder, folder);
1183 TFile::Open(file2);
1184 mc2->LoadHistograms();
1185 histMC2 = mc2->GetMultiplicityMC(histID, eventType)->ProjectionY("mc2", 1, mc2->GetMultiplicityMC(histID, eventType)->GetNbinsX());
1186 histMC2->Scale(1.0 / histMC2->Integral());
1187
1188 c = new TCanvas("CompareMC", "CompareMC", 800, 600);
1189 c->SetRightMargin(0.05);
1190 c->SetTopMargin(0.05);
1191 c->SetLogy();
1192
1193 histMC1->SetTitle(";N;P(N)");
1194 histMC1->SetStats(kFALSE);
1195 histMC1->GetXaxis()->SetRangeUser(0, 100);
1196
1197 histMC1->SetLineColor(1);
1198 histMC2->SetLineColor(2);
1199
1200 histMC1->Draw();
1201 histMC2->Draw("SAME");
1202
1203 legend = new TLegend(0.2, 0.2, 0.4, 0.4);
1204 legend->SetFillColor(0);
1205
1206 legend->AddEntry(histMC1, label1);
1207 legend->AddEntry(histMC2, label2);
1208
1209 legend->Draw();
1210
0a0f4adb 1211 c->SaveAs(Form("%s.gif", c->GetName()));
1212 c->SaveAs(Form("%s.eps", c->GetName()));
447c325d 1213}
1214
1215void* fit2Step(const char* fileNameMC = "multiplicityMC_2M.root", const char* fileNameESD = "multiplicityMC_1M_3.root", Int_t histID = 3, Bool_t fullPhaseSpace = kFALSE)
1216{
1217 gSystem->Load("libPWG0base");
1218
1219 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1220
1221 TFile::Open(fileNameMC);
1222 mult->LoadHistograms("Multiplicity");
1223
1224 TFile::Open(fileNameESD);
1225 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
1226 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", ((fullPhaseSpace) ? 4 : histID)));
1227 //hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityINEL%d", histID));
1228
1229 mult->SetMultiplicityESD(histID, hist);
1230
1231 // small hack to get around charge conservation for full phase space ;-)
1232 if (fullPhaseSpace)
1233 {
1234 TH1* corr = mult->GetCorrelation(histID + 4);
1235
1236 for (Int_t i=2; i<=corr->GetNbinsX(); i+=2)
1237 for (Int_t j=1; j<=corr->GetNbinsY(); ++j)
1238 {
1239 corr->SetBinContent(i, j, corr->GetBinContent(i-1, j));
1240 corr->SetBinError(i, j, corr->GetBinError(i-1, j));
1241 }
1242 }
1243
1244 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
1245 mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
1246 mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
1247
1248 TH1* result = (TH1*) mult->GetMultiplicityESDCorrected((fullPhaseSpace) ? 4 : histID))->Clone("firstresult");
1249
1250 mult->SetRegularizationParameters(AliMultiplicityCorrection::kEntropy, 100000);
1251 mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE, result);
1252 mult->DrawComparison("MinuitChi2_Step2", histID, fullPhaseSpace, kTRUE, hist2->ProjectionY("mymchist"));
0a173978 1253
1254 return mult;
1255}
cfc19dd5 1256
1257const char* GetRegName(Int_t type)
1258{
1259 switch (type)
1260 {
1261 case AliMultiplicityCorrection::kNone: return "None"; break;
1262 case AliMultiplicityCorrection::kPol0: return "Pol0"; break;
1263 case AliMultiplicityCorrection::kPol1: return "Pol1"; break;
1264 case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
1265 case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break;
0b4bfd98 1266 case AliMultiplicityCorrection::kLog : return "Log"; break;
dd701109 1267 }
1268 return 0;
1269}
1270
1271const char* GetEventTypeName(Int_t type)
1272{
1273 switch (type)
1274 {
1275 case AliMultiplicityCorrection::kTrVtx: return "trigger, vertex"; break;
1276 case AliMultiplicityCorrection::kMB: return "minimum bias"; break;
1277 case AliMultiplicityCorrection::kINEL: return "inelastic"; break;
cfc19dd5 1278 }
1279 return 0;
1280}
1281
69b09e3b 1282void EvaluateBayesianMethodIterationsSmoothing(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", const char* targetDir = "bayesian", Int_t histID = 1)
447c325d 1283{
1284 gSystem->mkdir(targetDir);
1285
1286 gSystem->Load("libPWG0base");
1287
1288 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1289 TFile::Open(fileNameMC);
1290 mult->LoadHistograms("Multiplicity");
1291
1292 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
1293 TFile::Open(fileNameESD);
1294 multESD->LoadHistograms("Multiplicity");
1295 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1296
1297 Int_t count = 0; // just to order the saved images...
1298
1299 TFile* graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir), "RECREATE");
1300
69b09e3b 1301 Int_t colors[] = {1, 2, 3, 4, 6};
1302 Int_t markers[] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3, 4, 5, 6};
0b4bfd98 1303
1304 for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
1305 //for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kINEL; ++type)
447c325d 1306 {
1307 TString tmp;
1308 tmp.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
1309
1310 TCanvas* canvas = new TCanvas(tmp, tmp, 800, 600);
1311
0f67a57c 1312 for (Int_t i = 1; i <= 5; i++)
447c325d 1313 {
69b09e3b 1314 Int_t iterArray[5] = {2, 5, 10, 20, -1};
0b4bfd98 1315 Int_t iter = iterArray[i-1];
1316
1317 TGraph* fitResultsMC[3];
1318 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1319 {
1320 fitResultsMC[region] = new TGraph;
69b09e3b 1321 fitResultsMC[region]->SetTitle(Form("%d iterations - reg. %d", iter, region+1));
0b4bfd98 1322 fitResultsMC[region]->GetXaxis()->SetTitle("smoothing parameter #alpha");
1323 fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
1324 fitResultsMC[region]->SetName(Form("%s_MC_%d", tmp.Data(), i * AliMultiplicityCorrection::kQualityRegions + region - 2));
1325 fitResultsMC[region]->SetFillColor(0);
0f67a57c 1326 //fitResultsMC[region]->SetMarkerStyle(markers[(i-1) * AliMultiplicityCorrection::kQualityRegions + region]);
69b09e3b 1327 fitResultsMC[region]->SetMarkerStyle(markers[i-1]);
1328 fitResultsMC[region]->SetLineColor(colors[i-1]);
1329 fitResultsMC[region]->SetMarkerColor(colors[i-1]);
0b4bfd98 1330 }
1331
447c325d 1332 TGraph* fitResultsRes = new TGraph;
1333 fitResultsRes->SetTitle(Form("%d iterations", iter));
1334 fitResultsRes->SetName(Form("%s_Res_%d", tmp.Data(), i));
1335 fitResultsRes->GetXaxis()->SetTitle("smoothing parameter");
1336 fitResultsRes->GetYaxis()->SetTitle("P_{2}");
1337
447c325d 1338 fitResultsRes->SetFillColor(0);
69b09e3b 1339 fitResultsRes->SetMarkerStyle(markers[i-1]);
1340 fitResultsRes->SetMarkerColor(colors[i-1]);
1341 fitResultsRes->SetLineColor(colors[i-1]);
447c325d 1342
1343 for (Float_t weight = 0.0; weight < 1.01; weight += 0.2)
1344 {
1345 Float_t chi2MC = 0;
1346 Float_t residuals = 0;
1347
0b4bfd98 1348 mult->ApplyBayesianMethod(histID, kFALSE, type, weight, iter, 0, kFALSE);
447c325d 1349 mult->DrawComparison(Form("%s/BayesianIterSmooth_%03d_%d_%d_%f", targetDir, count++, type, iter, weight), histID, kFALSE, kTRUE, multESD->GetMultiplicityMC(histID, type)->ProjectionY());
1350 mult->GetComparisonResults(&chi2MC, 0, &residuals);
1351
0b4bfd98 1352 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1353 fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
1354
447c325d 1355 fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
1356 }
1357
0b4bfd98 1358 graphFile->cd();
1359 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1360 fitResultsMC[region]->Write();
1361
447c325d 1362 fitResultsRes->Write();
1363 }
1364 }
0b4bfd98 1365
1366 graphFile->Close();
447c325d 1367}
1368
1369void EvaluateDrawResult(const char* targetDir, Int_t type = 0, Bool_t plotRes = kTRUE)
1370{
1371 gSystem->Load("libPWG0base");
1372
1373 TString name;
1374 TFile* graphFile = 0;
1375 if (type == -1)
1376 {
1377 name = "EvaluateChi2Method";
1378 graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
1379 }
1380 else
1381 {
1382 name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
1383 graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
1384 }
1385
1386 TCanvas* canvas = new TCanvas(name, name, 800, 500);
1387 if (type == -1)
0b4bfd98 1388 {
447c325d 1389 canvas->SetLogx();
0b4bfd98 1390 canvas->SetLogy();
1391 }
1392 canvas->SetTopMargin(0.05);
1393 canvas->SetGridx();
1394 canvas->SetGridy();
447c325d 1395
0b4bfd98 1396 TLegend* legend = new TLegend(0.8, 0.15, 0.98, 0.98);
447c325d 1397 legend->SetFillColor(0);
1398
1399 Int_t count = 1;
1400
1401 Float_t xMin = 1e20;
1402 Float_t xMax = 0;
1403
1404 Float_t yMin = 1e20;
1405 Float_t yMax = 0;
1406
0b4bfd98 1407 Float_t yMinRegion[3];
1408 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1409 yMinRegion[i] = 1e20;
1410
447c325d 1411 TString xaxis, yaxis;
1412
1413 while (1)
1414 {
1415 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1416 TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1417
0b4bfd98 1418 if (!mc)
447c325d 1419 break;
1420
1421 xaxis = mc->GetXaxis()->GetTitle();
1422 yaxis = mc->GetYaxis()->GetTitle();
1423
1424 mc->Print();
447c325d 1425
0b4bfd98 1426 if (res)
1427 res->Print();
447c325d 1428
0b4bfd98 1429 xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
1430 yMin = TMath::Min(yMin, mc->GetYaxis()->GetXmin());
447c325d 1431
0b4bfd98 1432 xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
1433 yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
1434
1435 if (plotRes && res)
447c325d 1436 {
0b4bfd98 1437 xMin = TMath::Min(xMin, res->GetXaxis()->GetXmin());
1438 yMin = TMath::Min(yMin, res->GetYaxis()->GetXmin());
447c325d 1439
0b4bfd98 1440 xMax = TMath::Max(xMax, res->GetXaxis()->GetXmax());
1441 yMax = TMath::Max(yMax, res->GetYaxis()->GetXmax());
447c325d 1442 }
0b4bfd98 1443
1444 for (Int_t i=0; i<mc->GetN(); ++i)
1445 yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
1446
1447 count++;
447c325d 1448 }
1449
0b4bfd98 1450 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1451 Printf("Minimum for region %d is %f", i, yMinRegion[i]);
1452
447c325d 1453 if (type >= 0)
1454 {
1455 xaxis = "smoothing parameter";
1456 }
1457 else if (type == -1)
1458 {
1459 xaxis = "weight parameter";
0b4bfd98 1460 xMax *= 5;
447c325d 1461 }
0b4bfd98 1462 //yaxis = "P_{1} (2 <= t <= 150)";
447c325d 1463
1464 printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
1465
1466 TGraph* dummy = new TGraph;
1467 dummy->SetPoint(0, xMin, yMin);
1468 dummy->SetPoint(1, xMax, yMax);
1469 dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
1470
1471 dummy->SetMarkerColor(0);
1472 dummy->Draw("AP");
0b4bfd98 1473 dummy->GetYaxis()->SetMoreLogLabels(1);
447c325d 1474
1475 count = 1;
1476
1477 while (1)
1478 {
1479 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1480 TGraph* res = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1481
0b4bfd98 1482 //printf("%s_MC_%d %p %p\n", name.Data(), count, mc, res);
447c325d 1483
0b4bfd98 1484 if (!mc)
447c325d 1485 break;
1486
1487 printf("Loaded %d sucessful.\n", count);
1488
44466df2 1489 if (type == -1)
447c325d 1490 {
44466df2 1491 legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
447c325d 1492 }
44466df2 1493 else
447c325d 1494 legend->AddEntry(mc);
1495
1496 mc->Draw("SAME PC");
1497
0b4bfd98 1498 if (plotRes && res)
447c325d 1499 {
1500 legend->AddEntry(res);
1501 res->Draw("SAME PC");
1502 }
1503
1504 count++;
1505 }
1506
1507 legend->Draw();
1508
1509 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
1510 canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
1511}
1512
0f67a57c 1513void EvaluateDrawResultRegions(const char* targetDir, Int_t type = 0)
1514{
69b09e3b 1515 loadlibs();
0f67a57c 1516
1517 TString name;
1518 TFile* graphFile = 0;
1519 if (type == -1)
1520 {
1521 name = "EvaluateChi2Method";
1522 graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir));
1523 }
1524 else
1525 {
1526 name.Form("EvaluateBayesianMethodIterationsSmoothing_%s", GetEventTypeName(type));
1527 graphFile = TFile::Open(Form("%s/EvaluateBayesianMethodIterationsSmoothing.root", targetDir));
1528 }
1529
69b09e3b 1530 TCanvas* canvas = new TCanvas(name, name, 1200, 800);
0f67a57c 1531 //canvas->Divide(1, AliMultiplicityCorrection::kQualityRegions, 0, 0);
69b09e3b 1532 canvas->SetTopMargin(0);
1533 canvas->SetBottomMargin(0);
1534 canvas->SetRightMargin(0.05);
0f67a57c 1535 canvas->Range(0, 0, 1, 1);
1536
69b09e3b 1537 TPad* pad[4];
1538 pad[3] = new TPad(Form("%s_pad1", name.Data()), "", 0, 0, 0.5, 0.5); pad[3]->SetTopMargin(0); pad[3]->SetBottomMargin(0.15);
1539 pad[1] = new TPad(Form("%s_pad2", name.Data()), "", 0, 0.5, 0.5, 1); pad[1]->SetBottomMargin(0);
1540 pad[0] = new TPad(Form("%s_pad3", name.Data()), "", 0.5, 0, 1, 0.5); pad[0]->SetTopMargin(0); pad[0]->SetBottomMargin(0.15);
1541 pad[2] = new TPad(Form("%s_pad4", name.Data()), "", 0.5, 0.5, 1, 1); pad[2]->SetBottomMargin(0);
0f67a57c 1542
1543 Float_t yMinRegion[3];
1544 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1545 yMinRegion[i] = 1e20;
1546
69b09e3b 1547 for (Int_t region = 0; region <= AliMultiplicityCorrection::kQualityRegions; region++)
0f67a57c 1548 {
1549 canvas->cd();
69b09e3b 1550 pad[region]->Draw();
1551 pad[region]->cd();
1552 pad[region]->SetRightMargin(0.05);
1553
0f67a57c 1554 if (type == -1)
1555 {
69b09e3b 1556 pad[region]->SetLogx();
0f67a57c 1557 }
69b09e3b 1558 pad[region]->SetLogy();
1559
1560 pad[region]->SetGridx();
1561 pad[region]->SetGridy();
0f67a57c 1562
69b09e3b 1563 TLegend* legend = new TLegend(0.5, 0.4, 0.85, 0.85);
0f67a57c 1564 legend->SetFillColor(0);
69b09e3b 1565 //legend->SetNColumns(3);
1566 legend->SetTextSize(0.06);
0f67a57c 1567 Int_t count = 0;
1568
1569 Float_t xMin = 1e20;
1570 Float_t xMax = 0;
1571
1572 Float_t yMin = 1e20;
1573 Float_t yMax = 0;
1574
1575 TString xaxis, yaxis;
1576
1577 while (1)
1578 {
1579 count++;
1580
69b09e3b 1581 if (region > 0)
1582 {
1583 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1584 if (!mc)
1585 break;
1586
1587 if (TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
1588 continue;
1589
1590 for (Int_t i=0; i<mc->GetN(); ++i)
1591 yMinRegion[(count-1) % 3] = TMath::Min(yMinRegion[(count-1) % 3], mc->GetY()[i]);
1592 }
1593 else
1594 {
1595 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1596 if (!mc)
1597 break;
1598 }
0f67a57c 1599
1600 xaxis = mc->GetXaxis()->GetTitle();
1601 yaxis = mc->GetYaxis()->GetTitle();
1602
69b09e3b 1603 //mc->Print();
0f67a57c 1604
1605 xMin = TMath::Min(xMin, mc->GetXaxis()->GetXmin());
69b09e3b 1606 yMin = TMath::Min(yMin, TMath::MinElement(mc->GetN(), mc->GetY()));
0f67a57c 1607
1608 xMax = TMath::Max(xMax, mc->GetXaxis()->GetXmax());
1609 yMax = TMath::Max(yMax, mc->GetYaxis()->GetXmax());
69b09e3b 1610
1611 //Printf("%f %f %f %f", xMin, xMax, yMin, yMax);
0f67a57c 1612 }
1613
1614 if (type >= 0)
1615 {
69b09e3b 1616 xaxis = "Smoothing parameter #alpha";
0f67a57c 1617 }
1618 else if (type == -1)
1619 {
69b09e3b 1620 xaxis = "Weight parameter #beta";
1621 //xMax *= 5;
1622 }
1623 if (region > 0)
1624 {
1625 yaxis = Form("Q_{1} in region %d ", region); // (2 <= t <= 150)";
0f67a57c 1626 }
69b09e3b 1627 else
1628 yaxis = "Q_{2} ";
0f67a57c 1629
1630 printf("%f %f %f %f\n", xMin, xMax, yMin, yMax);
1631
69b09e3b 1632 if (region > 0)
1633 {
1634 if (type == -1)
1635 {
1636 yMin = 0.534622;
1637 yMax = 19.228941;
1638 }
1639 else
1640 {
1641 yMin = 0.759109;
1642 yMax = 10;
1643 }
1644 }
1645
1646 if (type != -1)
1647 xMax = 1.0;
1648
1649 TH1* dummy = new TH2F("dummy", "dummy", 100, xMin, xMax, 100, yMin * 0.9, yMax * 1.1);
0f67a57c 1650 dummy->SetTitle(Form(";%s;%s", xaxis.Data(), yaxis.Data()));
1651
69b09e3b 1652 if (region == 0 && type != -1)
1653 dummy->GetYaxis()->SetMoreLogLabels(1);
1654 dummy->SetLabelSize(0.06, "xy");
1655 dummy->SetTitleSize(0.06, "xy");
1656 dummy->GetXaxis()->SetTitleOffset(1.2);
1657 dummy->GetYaxis()->SetTitleOffset(0.8);
1658 dummy->SetStats(0);
1659
1660 dummy->DrawCopy();
1661
0f67a57c 1662
1663 count = 0;
1664
1665 while (1)
1666 {
1667 count++;
1668
69b09e3b 1669 if (region > 0)
1670 {
1671 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_MC_%d", name.Data(), count));
1672 if (mc && TString(mc->GetTitle()).Contains(Form("reg. %d", region)) == kFALSE)
1673 continue;
1674 }
1675 else
1676 TGraph* mc = (TGraph*) graphFile->Get(Form("%s_Res_%d", name.Data(), count));
1677
0f67a57c 1678 if (!mc)
1679 break;
69b09e3b 1680 //printf("Loaded %d sucessful.\n", count);
0f67a57c 1681
1682 if (type == -1)
1683 {
69b09e3b 1684 //legend->AddEntry(mc, Form("Eq. (%d) - reg. %d", 10 + (count-1) / 3, 1+ (count-1) % 3));
1685 //legend->AddEntry(mc, Form("Eq. (7.%d)", 11 + (count-1) / 3));
1686 const char* names[] = { "Const", "Linear", "Log" };
1687 legend->AddEntry(mc, names[(count-1) / 3], "PL");
0f67a57c 1688 }
1689 else
69b09e3b 1690 {
1691 TString label(mc->GetTitle());
1692 TString newLabel(label(0, label.Index(" ")));
1693 //Printf("%s", newLabel.Data());
1694 if (newLabel.Atoi() == -1)
1695 {
1696 newLabel = "Convergence";
1697 }
1698 else
1699 newLabel = label(0, label.Index("iterations") + 10);
1700
1701 legend->AddEntry(mc, newLabel, "PL");
1702 }
0f67a57c 1703
1704 mc->Draw("SAME PC");
1705
69b09e3b 1706 }
0f67a57c 1707
69b09e3b 1708 if (region == 2)
1709 legend->Draw();
0f67a57c 1710 }
1711
1712 for (Int_t i=0; i<AliMultiplicityCorrection::kQualityRegions; ++i)
1713 Printf("Minimum for region %d is %f", i, yMinRegion[i]);
1714
1715 canvas->Modified();
1716 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
1717 canvas->SaveAs(Form("%s/%s.eps", targetDir, canvas->GetName()));
1718}
1719
69b09e3b 1720void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", const char* targetDir = "chi2compare", Int_t histID = 1)
dd701109 1721{
1722 gSystem->mkdir(targetDir);
1723
cfc19dd5 1724 gSystem->Load("libPWG0base");
1725
1726 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1727
1728 TFile::Open(fileNameMC);
1729 mult->LoadHistograms("Multiplicity");
1730
1731 TFile::Open(fileNameESD);
1732 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
1733 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
1734
1735 mult->SetMultiplicityESD(histID, hist);
1736
dd701109 1737 Int_t count = 0; // just to order the saved images...
69b09e3b 1738 Int_t colors[] = {1, 2, 4, 3};
1739 Int_t markers[] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 3};
dd701109 1740
447c325d 1741 TGraph* fitResultsRes = 0;
1742
1743 TFile* graphFile = TFile::Open(Form("%s/EvaluateChi2Method.root", targetDir), "RECREATE");
cfc19dd5 1744
69b09e3b 1745 //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kLog; ++type)
1746 //for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kEntropy; type <= AliMultiplicityCorrection::kEntropy; ++type)
1747 for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol1; type <= AliMultiplicityCorrection::kPol1; ++type)
cfc19dd5 1748 {
0b4bfd98 1749 TGraph* fitResultsMC[3];
1750 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1751 {
1752 fitResultsMC[region] = new TGraph;
44466df2 1753 fitResultsMC[region]->SetTitle(Form("Eq. (%d) - reg. %d", type+9, region+1));
0b4bfd98 1754 fitResultsMC[region]->GetXaxis()->SetTitle("weight parameter #alpha");
1755 fitResultsMC[region]->GetYaxis()->SetTitle(Form("P_{1} in region %d", region));
1756 fitResultsMC[region]->SetName(Form("EvaluateChi2Method_MC_%d", type * AliMultiplicityCorrection::kQualityRegions + region - 2));
1757 fitResultsMC[region]->SetFillColor(0);
69b09e3b 1758 //fitResultsMC[region]->SetMarkerStyle(markers[(type-1) * AliMultiplicityCorrection::kQualityRegions + region]);
1759 //fitResultsMC[region]->SetLineColor(colors[region]);
1760 fitResultsMC[region]->SetMarkerStyle(markers[type-1]);
1761 fitResultsMC[region]->SetMarkerColor(colors[type-1]);
1762 fitResultsMC[region]->SetLineColor(colors[type-1]);
0b4bfd98 1763 }
1764
447c325d 1765 fitResultsRes = new TGraph;
1766 fitResultsRes->SetTitle(Form("%s residual chi2", GetRegName(type)));
1767 fitResultsRes->SetName(Form("EvaluateChi2Method_Res_%d", type));
1768 fitResultsRes->GetXaxis()->SetTitle("Weight Parameter");
cfc19dd5 1769
cfc19dd5 1770 fitResultsRes->SetFillColor(0);
69b09e3b 1771 fitResultsRes->SetMarkerStyle(markers[type-1]);
1772 fitResultsRes->SetMarkerColor(colors[type-1]);
1773 fitResultsRes->SetLineColor(colors[type-1]);
cfc19dd5 1774
69b09e3b 1775 for (Int_t i=0; i<15; ++i)
447c325d 1776 {
69b09e3b 1777 Float_t weight = TMath::Power(TMath::Sqrt(10), i+1);
447c325d 1778 //Float_t weight = TMath::Power(10, i+2);
cfc19dd5 1779
447c325d 1780 //if (type == AliMultiplicityCorrection::kEntropy) weight = 1e4 * (i+1) * 1.5;
cfc19dd5 1781
cfc19dd5 1782 Float_t chi2MC = 0;
1783 Float_t residuals = 0;
447c325d 1784 Float_t chi2Limit = 0;
1785
1786 TString runName;
1787 runName.Form("MinuitChi2_%02d_%d_%f", count++, type, weight);
cfc19dd5 1788
1789 mult->SetRegularizationParameters(type, weight);
447c325d 1790 Int_t status = mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
69b09e3b 1791 mult->DrawComparison(Form("%s/%s", targetDir, runName.Data()), histID, kFALSE, kTRUE, hist2->ProjectionY("mymc", 1, hist2->GetNbinsX()));
447c325d 1792 if (status != 0)
1793 {
1794 printf("MINUIT did not succeed. Skipping...\n");
1795 continue;
1796 }
1797
dd701109 1798 mult->GetComparisonResults(&chi2MC, 0, &residuals);
447c325d 1799 TH1* result = mult->GetMultiplicityESDCorrected(histID);
1800 result->SetName(runName);
1801 result->Write();
cfc19dd5 1802
0b4bfd98 1803 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1804 fitResultsMC[region]->SetPoint(fitResultsMC[region]->GetN(), weight, mult->GetQuality(region));
1805
cfc19dd5 1806 fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
cfc19dd5 1807 }
1808
447c325d 1809 graphFile->cd();
0b4bfd98 1810 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1811 fitResultsMC[region]->Write();
447c325d 1812 fitResultsRes->Write();
cfc19dd5 1813 }
1814
447c325d 1815 graphFile->Close();
dd701109 1816}
1817
1818void EvaluateChi2MethodAll()
1819{
1820 EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
1821 EvaluateChi2Method("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
1822 EvaluateChi2Method("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
1823 EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
1824 EvaluateChi2Method("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
1825}
1826
1827void EvaluateBayesianMethodAll()
1828{
1829 EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M.root", "eval-3M-3M");
1830 EvaluateBayesianMethod("multiplicityMC_2M.root", "multiplicityMC_1M_3.root", "eval-2M-1M");
1831 EvaluateBayesianMethod("multiplicityMC_3M.root", "multiplicityMC_3M_NBD.root", "eval-3M-NBD");
1832 EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_1M_3.root", "eval-2MS-1M");
1833 EvaluateBayesianMethod("multiplicityMC_2M_smoothed.root", "multiplicityMC_3M_NBD.root", "eval-2MS-NBD");
1834}
1835
447c325d 1836void CompareMethods(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", const char* targetDir, Int_t histID = 3)
dd701109 1837{
1838 gSystem->mkdir(targetDir);
1839
1840 gSystem->Load("libPWG0base");
1841
1842 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1843
1844 TFile::Open(fileNameMC);
1845 mult->LoadHistograms("Multiplicity");
1846
1847 TFile::Open(fileNameESD);
1848 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
1849 multESD->LoadHistograms("Multiplicity");
1850
1851 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1852
447c325d 1853 TCanvas* canvas = new TCanvas("CompareMethods", "CompareMethods", 1200, 1200);
1854 canvas->Divide(3, 3);
dd701109 1855
1856 Int_t count = 0;
1857
447c325d 1858 for (AliMultiplicityCorrection::EventType type = AliMultiplicityCorrection::kTrVtx; type <= AliMultiplicityCorrection::kTrVtx; ++type)
dd701109 1859 {
447c325d 1860 TH1* mc = multESD->GetMultiplicityMC(histID, type)->ProjectionY("mymc");
1861 mc->Sumw2();
dd701109 1862
447c325d 1863 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
dd701109 1864 mult->ApplyMinuitFit(histID, kFALSE, type);
1865 mult->DrawComparison(Form("%s/CompareMethods_%d_MinuitChi2", targetDir, type), histID, kFALSE, kTRUE, mc);
1866 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("chi2Result");
1867
1868 mult->ApplyBayesianMethod(histID, kFALSE, type, 0.1);
1869 mult->DrawComparison(Form("%s/CompareMethods_%d_Bayesian", targetDir, type), histID, kFALSE, kTRUE, mc);
1870 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone("bayesResult");
1871
1872 mc->GetXaxis()->SetRangeUser(0, 150);
1873 chi2Result->GetXaxis()->SetRangeUser(0, 150);
1874
447c325d 1875/* // skip errors for now
dd701109 1876 for (Int_t i=1; i<=chi2Result->GetNbinsX(); ++i)
1877 {
1878 chi2Result->SetBinError(i, 0);
1879 bayesResult->SetBinError(i, 0);
447c325d 1880 }*/
dd701109 1881
1882 canvas->cd(++count);
1883 mc->SetFillColor(kYellow);
1884 mc->DrawCopy();
1885 chi2Result->SetLineColor(kRed);
1886 chi2Result->DrawCopy("SAME");
1887 bayesResult->SetLineColor(kBlue);
1888 bayesResult->DrawCopy("SAME");
1889 gPad->SetLogy();
1890
1891 canvas->cd(count + 3);
447c325d 1892 chi2ResultRatio = (TH1*) chi2Result->Clone("chi2ResultRatio");
1893 bayesResultRatio = (TH1*) bayesResult->Clone("bayesResultRatio");
1894 chi2ResultRatio->Divide(chi2Result, mc, 1, 1, "");
1895 bayesResultRatio->Divide(bayesResult, mc, 1, 1, "");
dd701109 1896
447c325d 1897 chi2ResultRatio->GetYaxis()->SetRangeUser(0.5, 1.5);
dd701109 1898
447c325d 1899 chi2ResultRatio->DrawCopy("HIST");
1900 bayesResultRatio->DrawCopy("SAME HIST");
dd701109 1901
447c325d 1902 canvas->cd(count + 6);
1903 chi2Result->Divide(chi2Result, bayesResult, 1, 1, "");
1904 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
1905 chi2Result->DrawCopy("HIST");
dd701109 1906 }
1907
1908 canvas->SaveAs(Form("%s/%s.gif", targetDir, canvas->GetName()));
1909 canvas->SaveAs(Form("%s/%s.C", targetDir, canvas->GetName()));
1910}
1911
447c325d 1912void StatisticsPlot(const char* fileNameMC = "multiplicityMC_2M.root", Int_t histID = 3)
dd701109 1913{
1914 gSystem->Load("libPWG0base");
1915
1916 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
1917
1918 TFile::Open(fileNameMC);
1919 mult->LoadHistograms("Multiplicity");
1920
447c325d 1921 const char* files[] = { "multiplicityMC_100k_1.root", "multiplicityMC_200k.root", "multiplicityMC_400k.root", "multiplicityMC_600k.root", "multiplicityMC_800k.root" };
dd701109 1922
0b4bfd98 1923 TGraph* fitResultsChi2[3];
1924 TGraph* fitResultsBayes[3];
1925
1926 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1927 {
1928 fitResultsChi2[region] = new TGraph;
1929 fitResultsChi2[region]->SetTitle(";Nevents;Chi2");
1930 fitResultsChi2[region]->SetName(Form("fitResultsChi2_%d", region));
1931 fitResultsChi2[region]->SetMarkerStyle(20+region);
1932
1933 fitResultsBayes[region] = new TGraph;
1934 fitResultsBayes[region]->SetTitle(";Nevents;Chi2");
1935 fitResultsBayes[region]->SetName(Form("fitResultsBayes_%d", region));
1936 fitResultsBayes[region]->SetMarkerStyle(20+region);
1937 fitResultsBayes[region]->SetMarkerColor(2);
1938 }
1939
447c325d 1940 TGraph* fitResultsChi2Limit = new TGraph; fitResultsChi2Limit->SetTitle(";Nevents;Multiplicity reach");
1941 TGraph* fitResultsBayesLimit = new TGraph; fitResultsBayesLimit->SetTitle(";Nevents;Multiplicity reach");
0b4bfd98 1942 TGraph* fitResultsChi2Res = new TGraph; fitResultsChi2Res->SetTitle(";Nevents;Chi2");
1943 TGraph* fitResultsBayesRes = new TGraph; fitResultsBayesRes->SetTitle(";Nevents;Chi2");
447c325d 1944
447c325d 1945 fitResultsChi2Limit->SetName("fitResultsChi2Limit");
1946 fitResultsBayesLimit->SetName("fitResultsBayesLimit");
0b4bfd98 1947 fitResultsChi2Res->SetName("fitResultsChi2Res");
1948 fitResultsBayesRes->SetName("fitResultsBayesRes");
dd701109 1949
1950 TCanvas* canvas = new TCanvas("StatisticsPlot", "StatisticsPlot", 1200, 600);
1951 canvas->Divide(5, 2);
1952
1953 Float_t min = 1e10;
1954 Float_t max = 0;
1955
447c325d 1956 TFile* file = TFile::Open("StatisticsPlot.root", "RECREATE");
1957 file->Close();
1958
dd701109 1959 for (Int_t i=0; i<5; ++i)
1960 {
1961 TFile::Open(files[i]);
1962 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
1963 multESD->LoadHistograms("Multiplicity");
1964
1965 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
1966 Int_t nEntries = multESD->GetMultiplicityESD(histID)->GetEntries();
447c325d 1967 TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY(Form("mc_%d", i));
dd701109 1968
447c325d 1969 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 10000);
dd701109 1970 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
1971 mult->DrawComparison(Form("StatisticsPlot_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
1972
dd701109 1973 Int_t chi2MCLimit = 0;
0b4bfd98 1974 Float_t chi2Residuals = 0;
1975 mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
1976 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1977 {
1978 fitResultsChi2[region]->SetPoint(fitResultsChi2[region]->GetN(), nEntries, mult->GetQuality(region));
1979 min = TMath::Min(min, mult->GetQuality(region));
1980 max = TMath::Max(max, mult->GetQuality(region));
1981 }
dd701109 1982 fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), nEntries, chi2MCLimit);
0b4bfd98 1983 fitResultsChi2Res->SetPoint(fitResultsChi2Res->GetN(), nEntries, chi2Residuals);
dd701109 1984
447c325d 1985 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
dd701109 1986
0b4bfd98 1987 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 100, 0, kFALSE);
dd701109 1988 mult->DrawComparison(Form("StatisticsPlot_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
447c325d 1989 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
0b4bfd98 1990 mult->GetComparisonResults(0, &chi2MCLimit, &chi2Residuals);
1991 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
1992 {
1993 fitResultsBayes[region]->SetPoint(fitResultsBayes[region]->GetN(), nEntries, mult->GetQuality(region));
1994 min = TMath::Min(min, mult->GetQuality(region));
1995 max = TMath::Max(max, mult->GetQuality(region));
1996 }
dd701109 1997 fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), nEntries, chi2MCLimit);
0b4bfd98 1998 fitResultsBayesRes->SetPoint(fitResultsBayesRes->GetN(), nEntries, chi2Residuals);
dd701109 1999
dd701109 2000 mc->GetXaxis()->SetRangeUser(0, 150);
2001 chi2Result->GetXaxis()->SetRangeUser(0, 150);
2002
2003 // skip errors for now
2004 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2005 {
2006 chi2Result->SetBinError(j, 0);
2007 bayesResult->SetBinError(j, 0);
2008 }
2009
2010 canvas->cd(i+1);
2011 mc->SetFillColor(kYellow);
2012 mc->DrawCopy();
2013 chi2Result->SetLineColor(kRed);
2014 chi2Result->DrawCopy("SAME");
2015 bayesResult->SetLineColor(kBlue);
2016 bayesResult->DrawCopy("SAME");
2017 gPad->SetLogy();
2018
2019 canvas->cd(i+6);
2020 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
2021 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
2022
2023 // skip errors for now
2024 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2025 {
2026 chi2Result->SetBinError(j, 0);
2027 bayesResult->SetBinError(j, 0);
2028 }
2029
2030 chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
2031 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
2032
2033 chi2Result->DrawCopy("");
2034 bayesResult->DrawCopy("SAME");
447c325d 2035
2036 TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
2037 mc->Write();
2038 chi2Result->Write();
2039 bayesResult->Write();
2040 file->Close();
dd701109 2041 }
2042
2043 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
2044 canvas->SaveAs(Form("%s.C", canvas->GetName()));
2045
2046 TCanvas* canvas2 = new TCanvas("StatisticsPlot2", "StatisticsPlot2", 800, 400);
2047 canvas2->Divide(2, 1);
2048
2049 canvas2->cd(1);
dd701109 2050
0b4bfd98 2051 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
2052 {
2053 fitResultsChi2[region]->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
2054 fitResultsChi2[region]->Draw(((region == 0) ? "AP" : "P SAME"));
2055
2056 fitResultsBayes[region]->Draw("P SAME");
2057 }
dd701109 2058
2059 gPad->SetLogy();
2060
2061 canvas2->cd(2);
2062 fitResultsChi2Limit->SetMarkerStyle(20);
2063 fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
2064 fitResultsChi2Limit->Draw("AP");
2065
2066 fitResultsBayesLimit->SetMarkerStyle(3);
2067 fitResultsBayesLimit->SetMarkerColor(2);
2068 fitResultsBayesLimit->Draw("P SAME");
2069
2070 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
2071 canvas2->SaveAs(Form("%s.C", canvas2->GetName()));
447c325d 2072
2073 TFile* file = TFile::Open("StatisticsPlot.root", "UPDATE");
0b4bfd98 2074
2075 for (Int_t region=0; region<AliMultiplicityCorrection::kQualityRegions; ++region)
2076 {
2077 fitResultsChi2[region]->Write();
2078 fitResultsBayes[region]->Write();
2079 }
447c325d 2080 fitResultsChi2Limit->Write();
2081 fitResultsBayesLimit->Write();
0b4bfd98 2082 fitResultsChi2Res->Write();
2083 fitResultsBayesRes->Write();
447c325d 2084 file->Close();
dd701109 2085}
2086
69b09e3b 2087void StartingConditions(const char* fileNameMC = "multiplicityMC.root", Int_t histID = 1)
dd701109 2088{
69b09e3b 2089 loadlibs();
dd701109 2090
2091 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2092
2093 TFile::Open(fileNameMC);
2094 mult->LoadHistograms("Multiplicity");
2095
69b09e3b 2096 const char* files[] = { "multiplicityESD.root", "multiplicityESD_100k_1.root", "multiplicityESD_100k_2.root", "multiplicityESD_100k_1.root", "multiplicityESD_100k_2.root" }
dd701109 2097
2098 // this one we try to unfold
2099 TFile::Open(files[0]);
2100 AliMultiplicityCorrection* multESD = new AliMultiplicityCorrection("MultiplicityESD", "MultiplicityESD");
2101 multESD->LoadHistograms("Multiplicity");
2102 mult->SetMultiplicityESD(histID, multESD->GetMultiplicityESD(histID));
69b09e3b 2103 TH1* mc = multESD->GetMultiplicityMC(histID, AliMultiplicityCorrection::kTrVtx)->ProjectionY("mc", 1, 1);
dd701109 2104
2105 TGraph* fitResultsChi2 = new TGraph;
2106 fitResultsChi2->SetTitle(";Input Dist ID;Chi2");
2107 TGraph* fitResultsBayes = new TGraph;
2108 fitResultsBayes->SetTitle(";Input Dist ID;Chi2");
2109 TGraph* fitResultsChi2Limit = new TGraph;
2110 fitResultsChi2Limit->SetTitle(";Input Dist ID;Multiplicity reach");
2111 TGraph* fitResultsBayesLimit = new TGraph;
2112 fitResultsBayesLimit->SetTitle(";Input Dist ID;Multiplicity reach");
2113
2114 TCanvas* canvas = new TCanvas("StartingConditions", "StartingConditions", 1200, 600);
2115 canvas->Divide(8, 2);
2116
2117 TCanvas* canvas3 = new TCanvas("StartingConditions3", "StartingConditions3", 1000, 400);
2118 canvas3->Divide(2, 1);
2119
2120 Float_t min = 1e10;
2121 Float_t max = 0;
2122
2123 TH1* firstChi = 0;
2124 TH1* firstBayesian = 0;
2125 TH1* startCond = multESD->GetMultiplicityESD(histID)->ProjectionY("startCond");
2126
2127 TLegend* legend = new TLegend(0.7, 0.7, 1, 1);
2128
447c325d 2129 TFile* file = TFile::Open("StartingConditions.root", "RECREATE");
2130 mc->Write();
2131 file->Close();
2132
dd701109 2133 for (Int_t i=0; i<8; ++i)
2134 {
2135 if (i == 0)
2136 {
2137 startCond = (TH1*) mc->Clone("startCond2");
2138 }
2139 else if (i < 6)
2140 {
2141 TFile::Open(files[i-1]);
2142 AliMultiplicityCorrection* multESD2 = new AliMultiplicityCorrection("MultiplicityESD2", "MultiplicityESD2");
2143 multESD2->LoadHistograms("Multiplicity");
2144 startCond = multESD2->GetMultiplicityESD(histID)->ProjectionY("startCond");
2145 }
2146 else if (i == 6)
2147 {
69b09e3b 2148 func = new TF1("nbd", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 50);
2149
dd701109 2150 func->SetParNames("scaling", "averagen", "k");
2151 func->SetParLimits(0, 1e-3, 1e10);
2152 func->SetParLimits(1, 0.001, 1000);
2153 func->SetParLimits(2, 0.001, 1000);
2154
2155 func->SetParameters(1, 10, 2);
2156 for (Int_t j=2; j<=startCond->GetNbinsX(); j++)
2157 startCond->SetBinContent(j, func->Eval(j-1));
2158 }
2159 else if (i == 7)
2160 {
2161 for (Int_t j=1; j<=startCond->GetNbinsX(); j++)
2162 startCond->SetBinContent(j, 1);
2163 }
2164
69b09e3b 2165 mult->SetRegularizationParameters(AliMultiplicityCorrection::kPol0, 1e5);
dd701109 2166 mult->ApplyMinuitFit(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, kFALSE, startCond);
69b09e3b 2167 //mult->DrawComparison(Form("StartingConditions_%d_MinuitChi2", i), histID, kFALSE, kTRUE, mc);
dd701109 2168
2169 Float_t chi2MC = 0;
2170 Int_t chi2MCLimit = 0;
2171 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
2172 fitResultsChi2->SetPoint(fitResultsChi2->GetN(), i, chi2MC);
2173 fitResultsChi2Limit->SetPoint(fitResultsChi2Limit->GetN(), i, chi2MCLimit);
2174 min = TMath::Min(min, chi2MC);
2175 max = TMath::Max(max, chi2MC);
2176
447c325d 2177 TH1* chi2Result = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("chi2Result_%d", i));
dd701109 2178 if (!firstChi)
2179 firstChi = (TH1*) chi2Result->Clone("firstChi");
2180
69b09e3b 2181 mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx, 1, 10, startCond);
2182 //mult->DrawComparison(Form("StartingConditions_%d_Bayesian", i), histID, kFALSE, kTRUE, mc);
447c325d 2183 TH1* bayesResult = (TH1*) mult->GetMultiplicityESDCorrected(histID)->Clone(Form("bayesResult_%d", i));
dd701109 2184 if (!firstBayesian)
2185 firstBayesian = (TH1*) bayesResult->Clone("firstBayesian");
2186
2187 mult->GetComparisonResults(&chi2MC, &chi2MCLimit, 0);
2188 fitResultsBayes->SetPoint(fitResultsBayes->GetN(), i, chi2MC);
2189 fitResultsBayesLimit->SetPoint(fitResultsBayesLimit->GetN(), i, chi2MCLimit);
2190
447c325d 2191 TFile* file = TFile::Open("StartingConditions.root", "UPDATE");
2192 chi2Result->Write();
2193 bayesResult->Write();
2194 file->Close();
2195
dd701109 2196 min = TMath::Min(min, chi2MC);
2197 max = TMath::Max(max, chi2MC);
2198 mc->GetXaxis()->SetRangeUser(0, 150);
2199 chi2Result->GetXaxis()->SetRangeUser(0, 150);
2200
2201 // skip errors for now
2202 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2203 {
2204 chi2Result->SetBinError(j, 0);
2205 bayesResult->SetBinError(j, 0);
2206 }
2207
2208 canvas3->cd(1);
2209 TH1* tmp = (TH1*) chi2Result->Clone("tmp");
2210 tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
2211 tmp->Divide(firstChi);
2212 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
2213 tmp->GetXaxis()->SetRangeUser(0, 200);
2214 tmp->SetLineColor(i+1);
2215 legend->AddEntry(tmp, Form("%d", i));
2216 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
2217
2218 canvas3->cd(2);
2219 tmp = (TH1*) bayesResult->Clone("tmp");
2220 tmp->SetTitle("Difference to best initial conditions;Npart;Ratio");
2221 tmp->Divide(firstBayesian);
2222 tmp->SetLineColor(i+1);
2223 tmp->GetYaxis()->SetRangeUser(0.5, 1.5);
2224 tmp->GetXaxis()->SetRangeUser(0, 200);
2225 tmp->DrawCopy((i > 0) ? "SAME HIST" : "HIST");
2226
2227 canvas->cd(i+1);
2228 mc->SetFillColor(kYellow);
2229 mc->DrawCopy();
2230 chi2Result->SetLineColor(kRed);
2231 chi2Result->DrawCopy("SAME");
2232 bayesResult->SetLineColor(kBlue);
2233 bayesResult->DrawCopy("SAME");
2234 gPad->SetLogy();
2235
2236 canvas->cd(i+9);
2237 chi2Result->Divide(chi2Result, mc, 1, 1, "B");
2238 bayesResult->Divide(bayesResult, mc, 1, 1, "B");
2239
2240 // skip errors for now
2241 for (Int_t j=0; j<=chi2Result->GetNbinsX(); ++j)
2242 {
2243 chi2Result->SetBinError(j, 0);
2244 bayesResult->SetBinError(j, 0);
2245 }
2246
2247 chi2Result->SetTitle("Ratios;Npart;unfolded measured/MC");
2248 chi2Result->GetYaxis()->SetRangeUser(0.5, 1.5);
2249
2250 chi2Result->DrawCopy("");
2251 bayesResult->DrawCopy("SAME");
2252 }
2253
2254 canvas3->cd(1);
2255 legend->Draw();
2256
cfc19dd5 2257 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
dd701109 2258
2259 TCanvas* canvas2 = new TCanvas("StartingConditions2", "StartingConditions2", 800, 400);
2260 canvas2->Divide(2, 1);
2261
2262 canvas2->cd(1);
2263 fitResultsChi2->SetMarkerStyle(20);
2264 fitResultsChi2->GetYaxis()->SetRangeUser(0.5 * min, 1.5 * max);
2265 fitResultsChi2->Draw("AP");
2266
2267 fitResultsBayes->SetMarkerStyle(3);
2268 fitResultsBayes->SetMarkerColor(2);
2269 fitResultsBayes->Draw("P SAME");
2270
2271 gPad->SetLogy();
2272
2273 canvas2->cd(2);
2274 fitResultsChi2Limit->SetMarkerStyle(20);
2275 fitResultsChi2Limit->GetYaxis()->SetRangeUser(0.9 * TMath::Min(fitResultsChi2Limit->GetYaxis()->GetXmin(), fitResultsBayesLimit->GetYaxis()->GetXmin()), 1.1 * TMath::Max(fitResultsChi2Limit->GetYaxis()->GetXmax(), fitResultsBayesLimit->GetYaxis()->GetXmax()));
2276 fitResultsChi2Limit->Draw("AP");
2277
2278 fitResultsBayesLimit->SetMarkerStyle(3);
2279 fitResultsBayesLimit->SetMarkerColor(2);
2280 fitResultsBayesLimit->Draw("P SAME");
2281
2282 canvas2->SaveAs(Form("%s.gif", canvas2->GetName()));
2283 canvas3->SaveAs(Form("%s.gif", canvas3->GetName()));
cfc19dd5 2284}
2285
2286void Merge(Int_t n, const char** files, const char* output)
2287{
447c325d 2288 // 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" };
2289
2290
cfc19dd5 2291 gSystem->Load("libPWG0base");
2292
2293 AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
2294 TList list;
2295 for (Int_t i=0; i<n; ++i)
2296 {
2297 TString name("Multiplicity");
2298 if (i > 0)
2299 name.Form("Multiplicity%d", i);
2300
2301 TFile::Open(files[i]);
eb9356d5 2302 Printf("Loading from file %s", files[i]);
cfc19dd5 2303 data[i] = new AliMultiplicityCorrection(name, name);
2304 data[i]->LoadHistograms("Multiplicity");
2305 if (i > 0)
2306 list.Add(data[i]);
2307 }
2308
2309 data[0]->Merge(&list);
2310
447c325d 2311 //data[0]->DrawHistograms();
cfc19dd5 2312
2313 TFile::Open(output, "RECREATE");
2314 data[0]->SaveHistograms();
2315 gFile->Close();
2316}
2317
2318void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
2319{
69b09e3b 2320 loadlibs();
cfc19dd5 2321
2322 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2323
2324 TFile::Open(fileName);
2325 mult->LoadHistograms("Multiplicity");
2326
2327 TF1* func = 0;
2328
2329 if (caseNo >= 4)
2330 {
69b09e3b 2331 TF1* func = new TF1("nbd", "exp(log([0]) + TMath::LnGamma([2]+x) - TMath::LnGamma([2]) - TMath::LnGamma(x+1) + log([1] / ([1]+[2])) * x + log(1.0 + [1]/[2]) * -[2])", 0, 200);
2332 //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, 200);
cfc19dd5 2333 func->SetParNames("scaling", "averagen", "k");
2334 }
2335
2336 switch (caseNo)
2337 {
0b4bfd98 2338 case 0: func = new TF1("flat", "1000"); break;
cfc19dd5 2339 case 1: func = new TF1("flat", "501-x"); break;
2340 case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
2341 case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
69b09e3b 2342 case 4: func->SetParameters(2e5, 15, 2); break;
447c325d 2343 case 5: func->SetParameters(1, 13, 7); break;
dd701109 2344 case 6: func->SetParameters(1e7, 30, 4); break;
0b4bfd98 2345 case 7: func->SetParameters(1e7, 30, 2); break; // ***
dd701109 2346 case 8: func = new TF1("testlaszlo", "10*1000*x*exp(-0.1*x)"); break;
cfc19dd5 2347
2348 default: return;
2349 }
2350
447c325d 2351 new TCanvas;
2352 func->Draw();
2353
69b09e3b 2354 mult->SetGenMeasFromFunc(func, 1);
cfc19dd5 2355
dd701109 2356 TFile::Open("out.root", "RECREATE");
2357 mult->SaveHistograms();
2358
69b09e3b 2359 new TCanvas; mult->GetMultiplicityESD(1)->ProjectionY()->DrawCopy();
2360 new TCanvas; mult->GetMultiplicityVtx(1)->ProjectionY("mc", 1, mult->GetMultiplicityVtx(1)->GetNbinsX())->DrawCopy();
0b4bfd98 2361
dd701109 2362 //mult->ApplyBayesianMethod(2, kFALSE);
2363 //mult->ApplyMinuitFit(2, kFALSE);
cfc19dd5 2364 //mult->ApplyGaussianMethod(2, kFALSE);
447c325d 2365 //mult->ApplyLaszloMethod(2, kFALSE, AliMultiplicityCorrection::kTrVtx);
dd701109 2366}
2367
70fdd197 2368void smoothCorrelationMap(const char* fileName = "multiplicity.root", Int_t corrMatrix = 1)
dd701109 2369{
2370 gSystem->Load("libPWG0base");
2371
2372 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2373
2374 TFile::Open(fileName);
2375 mult->LoadHistograms("Multiplicity");
2376
2377 // empty under/overflow bins in x, otherwise Project3D takes them into account
2378 TH3* corr = mult->GetCorrelation(corrMatrix);
2379 for (Int_t j=1; j<=corr->GetYaxis()->GetNbins(); ++j)
2380 {
2381 for (Int_t k=1; k<=corr->GetZaxis()->GetNbins(); ++k)
2382 {
2383 corr->SetBinContent(0, j, k, 0);
2384 corr->SetBinContent(corr->GetXaxis()->GetNbins()+1, j, k, 0);
2385 }
2386 }
2387
2388 TH2* proj = (TH2*) corr->Project3D("zy");
2389
2390 // normalize correction for given nPart
2391 for (Int_t i=1; i<=proj->GetNbinsX(); ++i)
2392 {
2393 Double_t sum = proj->Integral(i, i, 1, proj->GetNbinsY());
2394 if (sum <= 0)
2395 continue;
2396
2397 for (Int_t j=1; j<=proj->GetNbinsY(); ++j)
2398 {
2399 // npart sum to 1
2400 proj->SetBinContent(i, j, proj->GetBinContent(i, j) / sum);
2401 proj->SetBinError(i, j, proj->GetBinError(i, j) / sum);
2402 }
2403 }
2404
2405 new TCanvas;
447c325d 2406 proj->DrawCopy("COLZ");
dd701109 2407
2408 TH1* scaling = proj->ProjectionY("scaling", 1, 1);
2409 scaling->Reset();
2410 scaling->SetMarkerStyle(3);
2411 //scaling->GetXaxis()->SetRangeUser(0, 50);
2412 TH1* mean = (TH1F*) scaling->Clone("mean");
2413 TH1* width = (TH1F*) scaling->Clone("width");
2414
2415 TF1* lognormal = new TF1("lognormal", "[0]*exp(-(log(x)-[1])^2/(2*[2]^2))/(x*[2]*TMath::Sqrt(2*TMath::Pi()))", 0.01, 500);
2416 lognormal->SetParNames("scaling", "mean", "sigma");
2417 lognormal->SetParameters(1, 1, 1);
2418 lognormal->SetParLimits(0, 1, 1);
2419 lognormal->SetParLimits(1, 0, 100);
2420 lognormal->SetParLimits(2, 1e-3, 1);
2421
2422 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);
2423 nbd->SetParNames("scaling", "averagen", "k");
2424 nbd->SetParameters(1, 13, 5);
2425 nbd->SetParLimits(0, 1, 1);
2426 nbd->SetParLimits(1, 1, 100);
2427 nbd->SetParLimits(2, 1, 1e8);
2428
2429 TF1* poisson = new TF1("poisson", "[0] * exp(-(x+[2])) * (x+[2])**[1] / TMath::Factorial([1])", 0.01, 50);
2430 poisson->SetParNames("scaling", "k", "deltax");
2431 poisson->SetParameters(1, 1, 0);
2432 poisson->SetParLimits(0, 0, 10);
2433 poisson->SetParLimits(1, 0.01, 100);
2434 poisson->SetParLimits(2, 0, 10);
2435
2436 TF1* mygaus = new TF1("mygaus", "[0] * exp(-(x-[1])**2 / 2 / [2] - [3] * log(x + [4]) / [5])", 0.01, 50);
2437 mygaus->SetParNames("scaling", "mean", "width", "scale2log", "logmean", "logwidth");
2438 mygaus->SetParameters(1, 0, 1, 1, 0, 1);
2439 mygaus->SetParLimits(2, 1e-5, 10);
2440 mygaus->SetParLimits(4, 1, 1);
2441 mygaus->SetParLimits(5, 1e-5, 10);
2442
2443 //TF1* sqrt = new TF1("sqrt", "[0] + [1] * sqrt((x + [3]) * [2])", 0, 50);
2444 TF1* sqrt = new TF1("sqrt", "[0] + (x + [1])**[2]", 0, 50);
2445 sqrt->SetParNames("ydelta", "exp", "xdelta");
2446 sqrt->SetParameters(0, 0, 1);
2447 sqrt->SetParLimits(1, 0, 10);
2448
2449 const char* fitWith = "gaus";
2450
70fdd197 2451 for (Int_t i=1; i<=80; ++i)
dd701109 2452 {
2453 printf("Fitting %d...\n", i);
2454
2455 TH1* hist = proj->ProjectionY(Form("proj%d", i), i, i, "e");
447c325d 2456
dd701109 2457 //hist->GetXaxis()->SetRangeUser(0, 50);
2458 //lognormal->SetParameter(0, hist->GetMaximum());
2459 hist->Fit(fitWith, "0 M", "");
2460
2461 TF1* func = hist->GetListOfFunctions()->FindObject(fitWith);
2462
447c325d 2463 if (0 && (i % 5 == 0))
dd701109 2464 {
447c325d 2465 pad = new TCanvas;
dd701109 2466 hist->Draw();
2467 func->Clone()->Draw("SAME");
447c325d 2468 pad->SetLogy();
dd701109 2469 }
2470
70fdd197 2471 scaling->SetBinContent(i, func->GetParameter(0));
2472 mean->SetBinContent(i, func->GetParameter(1));
2473 width->SetBinContent(i, func->GetParameter(2));
2474
2475 scaling->SetBinError(i, func->GetParError(0));
2476 mean->SetBinError(i, func->GetParError(1));
2477 width->SetBinError(i, func->GetParError(2));
dd701109 2478 }
2479
2480 TF1* log = new TF1("log", "[0] + [1] * log([2] * x)", 0.01, 500);
2481 log->SetParameters(0, 1, 1);
2482 log->SetParLimits(1, 0, 100);
2483 log->SetParLimits(2, 1e-3, 10);
2484
2485 TF1* over = new TF1("over", "[0] + [1] / (x+[2])", 0.01, 500);
2486 over->SetParameters(0, 1, 0);
2487 //over->SetParLimits(0, 0, 100);
2488 over->SetParLimits(1, 1e-3, 10);
2489 over->SetParLimits(2, 0, 100);
2490
2491 c1 = new TCanvas("fitparams", "fitparams", 1200, 400);
2492 c1->Divide(3, 1);
2493
2494 c1->cd(1);
2495 scaling->Draw("P");
2496
2497 //TF1* scalingFit = new TF1("mypol0", "[0]");
2498 TF1* scalingFit = over;
447c325d 2499 scaling->Fit(scalingFit, "", "", 3, 140);
2500 scalingFit->SetRange(0, 200);
2501 scalingFit->Draw("SAME");
70fdd197 2502
dd701109 2503 c1->cd(2);
2504 mean->Draw("P");
2505
2506 //TF1* meanFit = log;
2507 TF1* meanFit = new TF1("mypol1", "[0]+[1]*x");
447c325d 2508 mean->Fit(meanFit, "", "", 3, 140);
2509 meanFit->SetRange(0, 200);
2510 meanFit->Draw("SAME");
dd701109 2511
2512 c1->cd(3);
2513 width->Draw("P");
2514
2515 //TF1* widthFit = over;
447c325d 2516 TF1* widthFit = new TF1("mypol", "[0]+[1]*TMath::Sqrt([2]*x)");
2517 widthFit->SetParLimits(2, 1e-5, 1e5);
2518 width->Fit(widthFit, "", "", 5, 140);
2519 widthFit->SetRange(0, 200);
2520 widthFit->Draw("SAME");
70fdd197 2521
dd701109 2522 // build new correction matrix
2523 TH2* new = (TH2*) proj->Clone("new");
2524 new->Reset();
2525 Float_t x, y;
2526 for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
2527 {
2528 TF1* func = (TF1*) gROOT->FindObject(fitWith);
2529 x = new->GetXaxis()->GetBinCenter(i);
2530 //if (i == 1)
2531 // x = 0.1;
2532 x++;
2533 func->SetParameters(scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
2534 printf("%f %f %f %f\n", x, scalingFit->Eval(x), meanFit->Eval(x), widthFit->Eval(x));
2535
2536 for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
2537 {
447c325d 2538 if (i < 11)
dd701109 2539 {
2540 // leave bins 1..20 untouched
2541 new->SetBinContent(i, j, corr->Integral(1, corr->GetNbinsX(), i, i, j, j));
2542 }
2543 else
2544 {
2545 y = new->GetYaxis()->GetBinCenter(j);
2546 if (j == 1)
2547 y = 0.1;
2548 if (func->Eval(y) > 1e-4)
2549 new->SetBinContent(i, j, func->Eval(y));
2550 }
2551 }
2552 }
2553
2554 // fill 0 multiplicity bins, this cannot be done with the function because it does not accept 0
2555 // we take the values from the old response matrix
2556 //for (Int_t i=1; i<=new->GetXaxis()->GetNbins(); i+=1)
2557 // new->SetBinContent(i, 1, proj->GetBinContent(i, 1));
2558
2559 //for (Int_t j=1; j<=new->GetYaxis()->GetNbins(); j+=1)
2560 // new->SetBinContent(1, j, proj->GetBinContent(1, j));
2561
2562 // normalize correction for given nPart
2563 for (Int_t i=1; i<=new->GetNbinsX(); ++i)
2564 {
2565 Double_t sum = new->Integral(i, i, 1, proj->GetNbinsY());
2566 if (sum <= 0)
2567 continue;
2568
2569 for (Int_t j=1; j<=new->GetNbinsY(); ++j)
2570 {
2571 // npart sum to 1
2572 new->SetBinContent(i, j, new->GetBinContent(i, j) / sum);
2573 new->SetBinError(i, j, new->GetBinError(i, j) / sum);
2574 }
2575 }
2576
2577 new TCanvas;
2578 new->Draw("COLZ");
2579
2580 TH2* diff = (TH2*) new->Clone("diff");
2581 diff->Add(proj, -1);
2582
2583 new TCanvas;
2584 diff->Draw("COLZ");
2585 diff->SetMinimum(-0.05);
2586 diff->SetMaximum(0.05);
2587
2588 corr->Reset();
2589
2590 for (Int_t i=1; i<=new->GetNbinsX(); ++i)
2591 for (Int_t j=1; j<=new->GetNbinsY(); ++j)
70fdd197 2592 corr->SetBinContent(corr->GetXaxis()->GetNbins() / 2 + 1, i, j, new->GetBinContent(i, j));
dd701109 2593
2594 new TCanvas;
2595 corr->Project3D("zy")->Draw("COLZ");
2596
2597 TFile::Open("out.root", "RECREATE");
2598 mult->SaveHistograms();
2599
2600 TH1* proj1 = proj->ProjectionY("proj1", 36, 36);
2601 TH1* proj2 = new->ProjectionY("proj2", 36, 36);
2602 proj2->SetLineColor(2);
2603
2604 new TCanvas;
2605 proj1->Draw();
2606 proj2->Draw("SAME");
cfc19dd5 2607}
447c325d 2608
0b4bfd98 2609void buildCorrelationMap(const char* fileName = "multiplicityMC_2M.root", Int_t corrMatrix = 3)
2610{
2611 gSystem->Load("libPWG0base");
2612
2613 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2614
2615 TFile::Open(fileName);
2616 mult->LoadHistograms("Multiplicity");
2617
2618 TH3F* new = mult->GetCorrelation(corrMatrix);
2619 new->Reset();
2620
2621 TF1* func = new TF1("func", "gaus(0)");
2622
2623 Int_t vtxBin = new->GetNbinsX() / 2;
2624 if (vtxBin == 0)
2625 vtxBin = 1;
2626
2627 Float_t sigma = 2;
2628 for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1)
2629 {
2630 Float_t x = new->GetYaxis()->GetBinCenter(i);
2631 func->SetParameters(1, x * 0.8, sigma);
2632 //func->SetParameters(1, x, sigma);
2633
2634 for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1)
2635 {
2636 Float_t y = new->GetYaxis()->GetBinCenter(j);
2637
2638 // cut at 1 sigma
2639 if (TMath::Abs(y-x*0.8) < sigma)
2640 new->SetBinContent(vtxBin, i, j, func->Eval(y));
2641
2642 // test only bin 40 has smearing
2643 //if (x != 40)
2644 // new->SetBinContent(vtxBin, i, j, (i == j));
2645 }
2646 }
2647
2648 new TCanvas;
2649 new->Project3D("zy")->DrawCopy("COLZ");
2650
2651 TFile* file = TFile::Open("out.root", "RECREATE");
2652 mult->SetCorrelation(corrMatrix, new);
2653 mult->SaveHistograms();
2654 file->Close();
2655}
2656
447c325d 2657void GetCrossSections(const char* fileName)
2658{
2659 gSystem->Load("libPWG0base");
2660
2661 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2662
2663 TFile::Open(fileName);
2664 mult->LoadHistograms("Multiplicity");
2665
2666 TH1* xSection2 = mult->GetCorrelation(3)->Project3D("y")->Clone("xSection2");
2667 xSection2->Sumw2();
2668 xSection2->Scale(1.0 / xSection2->Integral());
2669
2670 TH1* xSection15 = mult->GetCorrelation(2)->Project3D("y")->Clone("xSection15");
2671 xSection15->Sumw2();
2672 xSection15->Scale(1.0 / xSection15->Integral());
2673
2674 TFile::Open("crosssection.root", "RECREATE");
2675 xSection2->Write();
2676 xSection15->Write();
2677 gFile->Close();
2678}
0b4bfd98 2679
44466df2 2680void AnalyzeSpeciesTree(const char* fileName)
2681{
2682 //
2683 // prints statistics about fParticleSpecies
2684 //
2685
2686 gSystem->Load("libPWG0base");
2687
2688 TFile::Open(fileName);
2689 TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies");
2690
2691 const Int_t nFields = 8;
2692 Long_t totals[8];
2693 for (Int_t i=0; i<nFields; i++)
2694 totals[i] = 0;
2695
2696 for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
2697 {
2698 fParticleSpecies->GetEvent(i);
2699
2700 Float_t* f = fParticleSpecies->GetArgs();
2701
2702 for (Int_t j=0; j<nFields; j++)
2703 totals[j] += f[j+1];
2704 }
2705
2706 for (Int_t i=0; i<nFields; i++)
2707 Printf("%d --> %ld", i, totals[i]);
2708}
2709
0b4bfd98 2710void BuildResponseFromTree(const char* fileName, const char* target)
2711{
2712 //
2713 // builds several response matrices with different particle ratios (systematic study)
eb9356d5 2714 // particle-species study
69b09e3b 2715 //
2716 // WARNING doesn't work uncompiled, see test.C
0b4bfd98 2717
69b09e3b 2718 loadlibs();
0b4bfd98 2719
2720 TFile::Open(fileName);
2721 TNtuple* fParticleSpecies = (TNtuple*) gFile->Get("fParticleSpecies");
2722
2723 TFile* file = TFile::Open(target, "RECREATE");
2724 file->Close();
2725
2726 Int_t tracks = 0; // control variables
2727 Int_t noLabel = 0;
6d81c2de 2728 Int_t secondaries = 0;
0b4bfd98 2729 Int_t doubleCount = 0;
2730
eb9356d5 2731 for (Int_t num = 0; num < 7; num++)
0b4bfd98 2732 {
eb9356d5 2733 Printf("%d", num);
2734
69b09e3b 2735 TString name;
2736 name.Form("Multiplicity_%d", num);
2737 AliMultiplicityCorrection* fMultiplicity = new AliMultiplicityCorrection(name, name);
0b4bfd98 2738
2739 Float_t ratio[4]; // pi, K, p, other
2740 for (Int_t i = 0; i < 4; i++)
2741 ratio[i] = 1;
2742
eb9356d5 2743 if (num == 1)
2744 ratio[1] = 0.7;
2745 if (num == 2)
2746 ratio[1] = 1.3;
2747
2748 if (num == 3)
2749 ratio[2] = 0.7;
2750 if (num == 4)
2751 ratio[2] = 1.3;
2752
2753 if (num == 5)
2754 ratio[3] = 0.7;
2755 if (num == 6)
2756 ratio[3] = 1.3;
2757
0b4bfd98 2758 for (Int_t i=0; i<fParticleSpecies->GetEntries(); i++)
2759 {
2760 fParticleSpecies->GetEvent(i);
eb9356d5 2761
0b4bfd98 2762
2763 Float_t* f = fParticleSpecies->GetArgs();
2764
2765 Float_t gene = 0;
2766 Float_t meas = 0;
2767
2768 for (Int_t j = 0; j < 4; j++)
2769 {
eb9356d5 2770 if (ratio[j] == 1)
2771 {
2772 gene += f[j+1];
2773 }
2774 else
2775 {
2776 for (Int_t k=0; k<f[j+1]; k++)
2777 {
2778 if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
2779 {
2780 gene += 1;
2781 }
2782 else if (ratio[j] > 1)
2783 {
2784 gene += 1;
2785 if (gRandom->Uniform() < ratio[j] - 1)
2786 gene += 1;
2787 }
2788 }
2789 }
2790
2791 if (ratio[j] == 1)
2792 {
2793 meas += f[j+1+4];
2794 }
2795 else
2796 {
2797 for (Int_t k=0; k<f[j+1+4]; k++)
2798 {
2799 if (ratio[j] < 1 && gRandom->Uniform() < ratio[j])
2800 {
2801 meas += 1;
2802 }
2803 else if (ratio[j] > 1)
2804 {
2805 meas += 1;
2806 if (gRandom->Uniform() < ratio[j] - 1)
2807 meas += 1;
2808 }
2809 }
2810 }
2811
0b4bfd98 2812 tracks += f[j+1+4];
2813 }
2814
2815 // add the ones w/o label
2816 tracks += f[9];
2817 noLabel += f[9];
2818
6d81c2de 2819 // secondaries are already part of meas!
2820 secondaries += f[10];
2821
2822 // double counted are already part of meas!
2823 doubleCount += f[11];
0b4bfd98 2824
6d81c2de 2825 // 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 2826 meas += f[9];
0b4bfd98 2827
2828 //Printf("%.f %.f %.f %.f %.f", f[5], f[6], f[7], f[8], f[9]);
2829
eb9356d5 2830 fMultiplicity->FillCorrection(f[0], gene, gene, gene, gene, meas, meas, meas);
69b09e3b 2831 // HACK all as kND = 1
eb9356d5 2832 fMultiplicity->FillGenerated(f[0], kTRUE, kTRUE, 1, gene, gene, gene, gene);
2833 fMultiplicity->FillMeasured(f[0], meas, meas, meas);
0b4bfd98 2834 }
2835
2836 //fMultiplicity->DrawHistograms();
2837
2838 TFile* file = TFile::Open(target, "UPDATE");
2839 fMultiplicity->SaveHistograms();
2840 file->Close();
2841
2842 if (num == 0)
6d81c2de 2843 {
2844 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);
2845 if ((Float_t) noLabel / tracks > 0.02)
2846 Printf("WARNING: More than 2%% of tracks without label, this might bias the study!");
2847 }
0b4bfd98 2848 }
2849}
2850
69b09e3b 2851void ParticleRatioStudy()
0b4bfd98 2852{
69b09e3b 2853 loadlibs();
0b4bfd98 2854
eb9356d5 2855 for (Int_t num = 0; num < 7; num++)
69b09e3b 2856 {
2857 TString target;
eb9356d5 2858 target.Form("chi2a_inel_species_%d.root", num);
2859
2860 correct("compositions.root", Form("Multiplicity_%d", num), "multiplicityESD.root", 1, -1, 0, -1, 2, target);
69b09e3b 2861 }
2862}
2863
eb9356d5 2864void MergeCrossSection(Int_t xsectionID, const char* output = "multiplicityMC_xsection.root")
69b09e3b 2865{
eb9356d5 2866 const char* files[] = { "multiplicitySD.root", "multiplicityDD.root", "multiplicityND.root" };
69b09e3b 2867
2868 loadlibs();
0b4bfd98 2869
2870 TFile::Open(output, "RECREATE");
2871 gFile->Close();
2872
eb9356d5 2873 AliMultiplicityCorrection* data[3];
2874 TList list;
0b4bfd98 2875
eb9356d5 2876 Float_t ref_SD = -1;
2877 Float_t ref_DD = -1;
2878 Float_t ref_ND = -1;
0b4bfd98 2879
eb9356d5 2880 Float_t error_SD = -1;
2881 Float_t error_DD = -1;
2882 Float_t error_ND = -1;
0b4bfd98 2883
eb9356d5 2884 gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawSystematics.C"));
2885 GetRelativeFractions(xsectionID, ref_SD, ref_DD, ref_ND, error_SD, error_DD, error_ND);
2886
2887 for (Int_t i=0; i<3; ++i)
2888 {
2889 TString name;
2890 name.Form("Multiplicity");
2891 if (i > 0)
2892 name.Form("Multiplicity_%d", i);
0b4bfd98 2893
eb9356d5 2894 TFile::Open(files[i]);
2895 data[i] = new AliMultiplicityCorrection(name, name);
2896 data[i]->LoadHistograms("Multiplicity");
2897 }
2898
2899 // TODO is the under/overflow bin scaled as well? --> seems like, verify anyway!
0b4bfd98 2900
eb9356d5 2901 // calculating relative
2902 Float_t sd = data[0]->GetMultiplicityINEL(0)->Integral(0, data[0]->GetMultiplicityINEL(0)->GetNbinsX()+1);
2903 Float_t dd = data[1]->GetMultiplicityINEL(0)->Integral(0, data[1]->GetMultiplicityINEL(0)->GetNbinsX()+1);
2904 Float_t nd = data[2]->GetMultiplicityINEL(0)->Integral(0, data[2]->GetMultiplicityINEL(0)->GetNbinsX()+1);
2905 Float_t total = nd + dd + sd;
2906
2907 nd /= total;
2908 sd /= total;
2909 dd /= total;
2910
2911 Printf("Ratios in the correction map are: ND=%f, DD=%f, SD=%f", nd, dd, sd);
2912
2913 Float_t ratio[3];
2914 ratio[0] = ref_SD / sd;
2915 ratio[1] = ref_DD / dd;
2916 ratio[2] = ref_ND / nd;
2917
2918 Printf("SD=%.2f, DD=%.2f, ND=%.2f",ratio[0], ratio[1], ratio[2]);
2919
2920 for (Int_t i=0; i<3; ++i)
2921 {
2922 // modify x-section
2923 for (Int_t j=0; j<AliMultiplicityCorrection::kMCHists; j++)
2924 {
2925 data[i]->GetMultiplicityVtx(j)->Scale(ratio[i]);
2926 data[i]->GetMultiplicityMB(j)->Scale(ratio[i]);
2927 data[i]->GetMultiplicityINEL(j)->Scale(ratio[i]);
2928 data[i]->GetMultiplicityNSD(j)->Scale(ratio[i]);
2929 }
0b4bfd98 2930
eb9356d5 2931 for (Int_t j=0; j<AliMultiplicityCorrection::kESDHists; j++)
2932 {
2933 data[i]->GetMultiplicityESD(j)->Scale(ratio[i]);
2934 data[i]->GetTriggeredEvents(j)->Scale(ratio[i]);
2935 data[i]->GetNoVertexEvents(j)->Scale(ratio[i]);
0b4bfd98 2936 }
eb9356d5 2937
2938 for (Int_t j=0; j<AliMultiplicityCorrection::kCorrHists; j++)
2939 data[i]->GetCorrelation(j)->Scale(ratio[i]);
2940
2941 if (i > 0)
2942 list.Add(data[i]);
2943 }
0b4bfd98 2944
eb9356d5 2945 printf("Case %d, %s: Entries in response matrix: SD: %.2f DD: %.2f ND: %.2f", xsectionID, data[0]->GetName(), data[0]->GetCorrelation(0)->Integral(), data[1]->GetCorrelation(0)->Integral(), data[2]->GetCorrelation(0)->Integral());
0b4bfd98 2946
eb9356d5 2947 data[0]->Merge(&list);
0b4bfd98 2948
eb9356d5 2949 Printf(" Total: %.2f", data[0]->GetCorrelation(0)->Integral());
0b4bfd98 2950
eb9356d5 2951 TFile::Open(output, "RECREATE");
2952 data[0]->SaveHistograms();
2953 gFile->Close();
0b4bfd98 2954
eb9356d5 2955 list.Clear();
0b4bfd98 2956
eb9356d5 2957 for (Int_t i=0; i<3; ++i)
2958 delete data[i];
0b4bfd98 2959}
44466df2 2960
2961void Rebin(const char* fileName = "multiplicityMC_3M.root", Int_t corrMatrix = 3)
2962{
2963 // rebins MC axis of correlation map, MC and histogram for corrected (for evaluation of effect of regularization)
2964 // rebin does not exist for 3D hists, so we convert to 2D and then back to 3D (loosing the vertex information)
2965
2966 Printf("WARNING: Vertex information is lost in this process. Use result only for evaluation of errors.");
2967
2968 gSystem->Load("libPWG0base");
2969
2970 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
2971
2972 TFile::Open(fileName);
2973 mult->LoadHistograms("Multiplicity");
2974
2975 // rebin correlation
2976 TH3* old = mult->GetCorrelation(corrMatrix);
2977
2978 // empty under/overflow bins in x, otherwise Project3D takes them into account
2979 for (Int_t y=1; y<=old->GetYaxis()->GetNbins(); ++y)
2980 {
2981 for (Int_t z=1; z<=old->GetZaxis()->GetNbins(); ++z)
2982 {
2983 old->SetBinContent(0, y, z, 0);
2984 old->SetBinContent(old->GetXaxis()->GetNbins()+1, y, z, 0);
2985 }
2986 }
2987
2988 TH2* response = (TH2*) old->Project3D("zy");
2989 response->RebinX(2);
2990
2991 TH3F* new = new TH3F(old->GetName(), old->GetTitle(),
2992 old->GetXaxis()->GetNbins(), old->GetXaxis()->GetBinLowEdge(1), old->GetXaxis()->GetBinUpEdge(old->GetXaxis()->GetNbins()),
2993 old->GetYaxis()->GetNbins() / 2, old->GetYaxis()->GetBinLowEdge(1), old->GetYaxis()->GetBinUpEdge(old->GetYaxis()->GetNbins()),
2994 old->GetZaxis()->GetNbins(), old->GetZaxis()->GetBinLowEdge(1), old->GetZaxis()->GetBinUpEdge(old->GetZaxis()->GetNbins()));
2995 new->Reset();
2996
2997 Int_t vtxBin = new->GetNbinsX() / 2;
2998 if (vtxBin == 0)
2999 vtxBin = 1;
3000
3001 for (Int_t i=1; i<=new->GetYaxis()->GetNbins(); i+=1)
3002 for (Int_t j=1; j<=new->GetZaxis()->GetNbins(); j+=1)
3003 new->SetBinContent(vtxBin, i, j, response->GetBinContent(i, j));
3004
3005 // rebin MC + hist for corrected
3006 for (AliMultiplicityCorrection::EventType eventType = AliMultiplicityCorrection::kTrVtx; eventType <= AliMultiplicityCorrection::kINEL; eventType++)
3007 mult->GetMultiplicityMC(corrMatrix, eventType)->RebinY(2);
3008
3009 mult->GetMultiplicityESDCorrected(corrMatrix)->Rebin(2);
3010
3011 // recreate measured from correlation matrix to get rid of vertex shift effect
3012 TH2* newMeasured = (TH2*) old->Project3D("zx");
3013 TH2* esd = mult->GetMultiplicityESD(corrMatrix);
3014 esd->Reset();
3015
3016 // transfer from TH2D to TH2F
3017 for (Int_t i=0; i<=new->GetXaxis()->GetNbins()+1; i+=1)
3018 for (Int_t j=0; j<=new->GetYaxis()->GetNbins()+1; j+=1)
3019 esd->SetBinContent(i, j, newMeasured->GetBinContent(i, j));
3020
3021 new TCanvas;
3022 new->Project3D("zy")->DrawCopy("COLZ");
3023
3024 TFile* file = TFile::Open("out.root", "RECREATE");
3025 mult->SetCorrelation(corrMatrix, new);
3026 mult->SaveHistograms();
3027 file->Close();
3028}
3029
3030void EvaluateRegularizationEffect(Int_t step, const char* fileNameRebinned = "multiplicityMC_3M_rebinned.root", const char* fileNameNormal = "multiplicityMC_3M.root", Int_t histID = 3)
3031{
3032 // due to some static members in AliMultiplicityCorrection, the session has to be restarted after changing the number of parameters, to be fixed
3033 // that is why this is done in 2 steps
3034
3035 gSystem->Load("libPWG0base");
3036
3037 Bool_t fullPhaseSpace = kFALSE;
3038
3039 if (step == 1)
3040 {
3041 // first step: unfold without regularization and rebinned histogram ("N=M")
3042 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3043 TFile::Open(fileNameRebinned);
3044 mult->LoadHistograms();
3045
3046 mult->SetRegularizationParameters(AliMultiplicityCorrection::kNone, 0, 125);
3047 mult->SetCreateBigBin(kFALSE);
3048
3049 mult->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
3050 mult->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult->GetMultiplicityVtx(histID)->ProjectionY("mymchist"));
3051
3052 TFile* file = TFile::Open("EvaluateRegularizationEffect1.root", "RECREATE");
3053 mult->SaveHistograms();
3054 file->Close();
3055 }
3056 else if (step == 2)
3057 {
3058 // second step: unfold with regularization and normal histogram
3059 AliMultiplicityCorrection* mult2 = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3060 TFile::Open(fileNameNormal);
3061 mult2->LoadHistograms();
3062
3063 mult2->SetRegularizationParameters(AliMultiplicityCorrection::kPol1, 1e4);
3064 mult2->SetCreateBigBin(kTRUE);
3065 mult2->ApplyMinuitFit(histID, fullPhaseSpace, AliMultiplicityCorrection::kTrVtx, kFALSE);
3066 mult2->DrawComparison("MinuitChi2", histID, fullPhaseSpace, kTRUE, mult2->GetMultiplicityVtx(histID)->ProjectionY("mymchist"));
3067
3068 TH1* result2 = mult2->GetMultiplicityESDCorrected(histID);
3069
3070 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3071 TFile* file = TFile::Open("EvaluateRegularizationEffect1.root");
3072 mult->LoadHistograms();
3073
3074 TH1* result1 = mult->GetMultiplicityESDCorrected(histID);
3075
3076 // compare results
3077 TCanvas* canvas = new TCanvas("EvaluateRegularizationEffect", "EvaluateRegularizationEffect", 1000, 800);
3078 canvas->Divide(2, 2);
3079
3080 canvas->cd(1);
3081 result1->SetLineColor(1);
3082 result1->DrawCopy();
3083 result2->SetLineColor(2);
3084 result2->DrawCopy("SAME");
3085 gPad->SetLogy();
3086
3087 result2->Rebin(2);
3088 result1->Scale(1.0 / result1->Integral());
3089 result2->Scale(1.0 / result2->Integral());
3090
3091 canvas->cd(2);
3092 result1->DrawCopy();
3093 result2->DrawCopy("SAME");
3094 gPad->SetLogy();
3095
3096 TH1* diff = (TH1*) result1->Clone("diff");
3097 diff->Add(result2, -1);
3098
3099 canvas->cd(3);
3100 diff->DrawCopy("HIST");
3101
3102 canvas->cd(4);
3103 diff->Divide(result1);
3104 diff->GetYaxis()->SetRangeUser(-0.3, 0.3);
3105 diff->DrawCopy("HIST");
3106
3107 Double_t chi2 = 0;
3108 for (Int_t i=1; i<=diff->GetNbinsX(); i++)
3109 chi2 += diff->GetBinContent(i) * diff->GetBinContent(i);
3110
3111 Printf("Chi2 is %e", chi2);
3112
3113 canvas->SaveAs(Form("%s.eps", canvas->GetName()));
3114 }
3115}
eb9356d5 3116
3117void MergeDistributions()
3118{
3119 loadlibs();
3120
3121 const char* dir1 = "run000104824-52_pass4";
3122 const char* dir2 = "run000104867_90_92_pass4";
3123
3124 for (Int_t evType = 0; evType < 2; evType++)
3125 {
3126 Printf("%d", evType);
3127
3128 const char* evTypeStr = ((evType == 0) ? "inel" : "nsd");
3129
3130 const char* id = "chi2a";
3131 //const char* id = "bayes";
3132
3133 TString suffix;
3134 suffix.Form("/all/spd/%s_", id);
3135 if (evType == 1)
3136 suffix.Form("/v0and/spd/%s_", id);
3137
3138 TString file1, file2;
3139 file1.Form("%s%s%%s.root", dir1, suffix.Data());
3140 file2.Form("%s%s%%s.root", dir2, suffix.Data());
3141
3142 if (1)
3143 {
3144 const char* files[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr) };
3145 Merge(2, files, Form("merged/%s_%s.root", id, evTypeStr));
3146 }
3147 else
3148 {
3149 AliMultiplicityCorrection* mult1 = AliMultiplicityCorrection::Open(Form(file1.Data(), evTypeStr));
3150 AliMultiplicityCorrection* mult2 = AliMultiplicityCorrection::Open(Form(file2.Data(), evTypeStr));
3151
3152 AliMultiplicityCorrection* target = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
3153
3154 for (Int_t etaRange = 0; etaRange < 3; etaRange++)
3155 {
3156 hist1 = mult1->GetMultiplicityESDCorrected(etaRange);
3157 hist2 = mult2->GetMultiplicityESDCorrected(etaRange);
3158 targetHist = target->GetMultiplicityESDCorrected(etaRange);
3159
3160 //hist1->Scale(1.0 / hist1->Integral());
3161 //hist2->Scale(1.0 / hist2->Integral());
3162
3163 //targetHist->SetBinContent(1, hist1->GetBinContent(1) * 0.5 + hist2->GetBinContent(1) * 0.5);
3164 targetHist->SetBinContent(1, hist1->GetBinContent(1) + hist2->GetBinContent(1));
3165 for (Int_t i=1; i<=hist1->GetNbinsX(); i++)
3166 {
3167 if (hist1->GetBinError(i) > 0 && hist2->GetBinError(i) > 0)
3168 {
3169 Float_t w1 = 1.0 / hist1->GetBinError(i) / hist1->GetBinError(i);
3170 Float_t w2 = 1.0 / hist2->GetBinError(i) / hist2->GetBinError(i);
3171
3172 Float_t average = (hist1->GetBinContent(i) * w1 + hist2->GetBinContent(i) * w2) / (w1 + w2);
3173
3174 //targetHist->SetBinContent(i, average);
3175 //targetHist->SetBinError(i, TMath::Max(mult1->GetBinError(i), mult2->GetBinError(i)));
3176
3177 targetHist->SetBinContent(i, hist1->GetBinContent(i) + hist2->GetBinContent(i));
3178 targetHist->SetBinError(i, hist1->GetBinError(i) + hist2->GetBinError(i));
3179 }
3180 }
3181 }
3182
3183 file = TFile::Open(Form("merged/%s_%s.root", id, evTypeStr), "RECREATE");
3184 target->SaveHistograms();
3185 file->Close();
3186 }
3187
3188 const char* files2[] = { Form(file1.Data(), evTypeStr), Form(file2.Data(), evTypeStr), Form("merged/%s_%s.root", id, evTypeStr) };
3189 drawAll(files2, (evType == 1), kTRUE);
3190 }
3191}
3192
3193void CompareDistributions(Int_t evType)
3194{
3195 loadlibs();
3196 gROOT->ProcessLine(gSystem->ExpandPathName(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C"));
3197
3198 const char* dir1 = "run000104824-52_pass4";
3199 const char* dir2 = "run000104867_90_92_pass4";
3200
3201 const char* evTypeStr = (evType == 0) ? "inel" : "nsd";
3202
3203 const char* suffix = "/all/spd/chi2a_";
3204 if (evType == 1)
3205 suffix = "/v0and/spd/chi2a_";
3206
3207 TString file1, file2;
3208 file1.Form("%s%s%s.root", dir1, suffix, evTypeStr);
3209 file2.Form("%s%s%s.root", dir2, suffix, evTypeStr);
3210
3211 for (Int_t hID = 0; hID < 3; hID++)
3212 CompareQualityHists(file1, file2, Form("Multiplicity/fMultiplicityESDCorrected%d", hID), 1, 1);
3213}
3214
3215void StatisticalUncertainty(Int_t methodType, Int_t etaRange, Bool_t mc = kFALSE)
3216{
3217 loadlibs();
3218
3219 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicityMC.root");
3220 AliMultiplicityCorrection* esd = AliMultiplicityCorrection::Open("multiplicityESD.root");
3221 AliMultiplicityCorrection* multTrigger = AliMultiplicityCorrection::Open("multiplicityTrigger.root");
3222
3223 TH1* mcHist = esd->GetMultiplicityMB(etaRange)->ProjectionY("mymc", 1, 1);
3224
3225 Int_t geneLimits[] = { 0, 0, 0 };
3226
3227 LoadAndInitialize(mult, esd, multTrigger, etaRange, kFALSE, geneLimits);
3228
3229 AliUnfolding::SetNbins(kBinLimits[etaRange], geneLimits[etaRange]);
3230 AliUnfolding::SetChi2Regularization(AliUnfolding::kPol0, 5);
3231 AliUnfolding::SetBayesianParameters(1, 10);
3232
3233 // background subtraction
3234 Int_t background = 0;
3235
3236 //background = 1091 + 4398; // MB1 for run 104824 - 52
3237 background = 417 + 1758; // MB1 for run 104867 - 92
3238
3239 //background = 20 + 0; // V0AND for run 104824 - 52
3240 //background = 10 + 0; // V0AND for run 104867 - 92
3241
3242 Printf("NOTE: Subtracting %d background events", background);
3243
3244 Int_t zeroBin = mult->GetTriggeredEvents(etaRange)->GetBinContent(1) - background - mult->GetMultiplicityESD(etaRange)->Integral(0, 2, 1, 1);
3245
3246 TH1* errorMeasured = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kFALSE, ((mc) ? mcHist : 0))->Clone("errorMeasured");
3247
3248 new TCanvas; errorMeasured->Draw();
3249 new TCanvas;
3250
3251 AliPWG0Helper::NormalizeToBinWidth(mult->GetMultiplicityESDCorrected(etaRange));
3252 mult->GetMultiplicityESDCorrected(etaRange)->Draw();
3253
3254 if (mc)
3255 {
3256 AliPWG0Helper::NormalizeToBinWidth(mcHist);
3257 mcHist->Scale(mult->GetMultiplicityESDCorrected(etaRange)->Integral() / mcHist->Integral());
3258 mcHist->SetLineColor(2);
3259 mcHist->Draw("SAME");
3260 }
3261 gPad->SetLogy();
3262
3263 TH1* errorResponse = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kFALSE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorResponse");
3264
3265 TH1* errorBoth = (TH1*) mult->StatisticalUncertainty((AliUnfolding::MethodType) methodType, etaRange, kFALSE, AliMultiplicityCorrection::kMB, zeroBin, kTRUE, kTRUE, ((mc) ? mcHist : 0))->Clone("errorBoth");
3266
3267 if (mc)
3268 {
3269 TH1* result = mult->GetMultiplicityESDCorrected(etaRange);
3270 DrawResultRatio(mcHist, result, "StatisticalUncertainty2.eps");
3271 }
3272
3273 TFile* file = new TFile(Form("StatisticalUncertaintySPD%s.root", (methodType == 0) ? "Chi2" : "Bayesian"), "RECREATE");
3274 errorResponse->Write();
3275 errorMeasured->Write();
3276 errorBoth->Write();
3277 file->Close();
3278}
3279
3280void ChangeResponseMatrixEfficiency(Bool_t reduceEff = kTRUE, Float_t factor = 0.01625)
3281{
3282 loadlibs();
3283
3284 AliMultiplicityCorrection* mult = AliMultiplicityCorrection::Open("multiplicity.root");
3285
3286 for (Int_t etaR = 0; etaR < 3; etaR++)
3287 {
3288 TH3* corr = mult->GetCorrelation(etaR);
3289 TH3* corrOld = (TH3*) corr->Clone("corrOld");
3290
3291 corr->Reset();
3292
3293 Int_t total = 0;
3294 Int_t change = 0;
3295 for (Int_t x=0; x<=corr->GetNbinsX()+1; x++)
3296 {
3297 for (Int_t y=0; y<=corr->GetNbinsY()+1; y++)
3298 {
3299 Printf("%d", y);
3300 for (Int_t z=0; z<=corr->GetNbinsZ()+1; z++)
3301 {
3302 Float_t tracklets = corr->GetZaxis()->GetBinCenter(z);
3303
3304 for (Int_t n=0; n<corrOld->GetBinContent(x, y, z); n++)
3305 {
3306 total += tracklets;
3307 Float_t trackletsNew = tracklets;
3308
3309 for (Int_t i=0; i<tracklets; i++)
3310 {
3311 if (gRandom->Uniform() < factor)
3312 {
3313 trackletsNew += (reduceEff) ? -1 : 1;
3314 change += (reduceEff) ? -1 : 1;
3315 }
3316 }
3317
3318 corr->Fill(corr->GetXaxis()->GetBinCenter(x), corr->GetYaxis()->GetBinCenter(y), trackletsNew);
3319 }
3320 }
3321 }
3322 }
3323
3324 Printf("%d: Changed %d out of %d total tracklets", etaR, change, total);
3325 }
3326
3327 file = TFile::Open("out.root", "RECREATE");
3328 mult->SaveHistograms();
3329 file->Close();
3330}
3331