Added transformation function (Marian)
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / runMultiplicitySelector.C
CommitLineData
0ab29cfa 1/* $Id$ */
2
3//
4// script to run the AliMultiplicityESDSelector
5//
6
7#include "../CreateESDChain.C"
0bd1f8a0 8#include "../PWG0Helper.C"
0ab29cfa 9
03244034 10void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "")
0ab29cfa 11{
0bd1f8a0 12 if (aProof)
cfc19dd5 13 connectProof("jgrosseo@lxb6046");
0bd1f8a0 14
15 TString libraries("libEG;libGeom;libESD;libPWG0base");
16 TString packages("PWG0base");
10ebe68d 17
0ab29cfa 18 if (aMC != kFALSE)
0bd1f8a0 19 {
20 libraries += ";libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6";
21 packages += ";PWG0dep";
22 }
23
24 if (!prepareQuery(libraries, packages, kTRUE))
25 return;
0ab29cfa 26
cfc19dd5 27 gProof->SetParallel(0);
28 gProof->SetLogLevel(4);
29 gProof->SetParallel(9999);
30
0ab29cfa 31 gROOT->ProcessLine(".L CreateCuts.C");
32 gROOT->ProcessLine(".L drawPlots.C");
33
0ab29cfa 34 // selection of esd tracks
35 AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
36 if (!esdTrackCuts)
37 {
38 printf("ERROR: esdTrackCuts could not be created\n");
39 return;
40 }
41
0bd1f8a0 42 TList inputList;
43 inputList.Add(esdTrackCuts);
44
cfc19dd5 45 TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE);
0ab29cfa 46
47 TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector");
48 AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
49
0bd1f8a0 50 selectorName += ".cxx+";
0ab29cfa 51
52 if (aDebug != kFALSE)
53 selectorName += "g";
54
03244034 55 Int_t result = executeQuery(chain, &inputList, selectorName, option);
0ab29cfa 56
57 if (result != 0)
58 {
59 printf("ERROR: Executing process failed with %d.\n", result);
60 return;
61 }
0a173978 62}
0ab29cfa 63
0a173978 64void draw(const char* fileName = "multiplicityMC.root")
65{
66 gSystem->Load("libPWG0base");
67
68 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
69
70 TFile::Open(fileName);
71 mult->LoadHistograms("Multiplicity");
72
73 mult->DrawHistograms();
0ab29cfa 74}
75
cfc19dd5 76void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
0a173978 77{
78 gSystem->Load("libPWG0base");
79
80 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
81
82 TFile::Open(fileName);
83 mult->LoadHistograms("Multiplicity");
84
cfc19dd5 85 mult->ApplyBayesianMethod(hist, kFALSE, eventType);
86 mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY());
87
88 //mult->ApplyMinuitFit(hist, kFALSE);
89 //mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
90
91 //mult->ApplyGaussianMethod(hist, kFALSE);
9ca6be09 92
93 return mult;
94}
95
cfc19dd5 96void* fitOther(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 2)
9ca6be09 97{
98 gSystem->Load("libPWG0base");
99
100 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
101
102 TFile::Open(fileNameMC);
103 mult->LoadHistograms("Multiplicity");
104
105 TFile::Open(fileNameESD);
cfc19dd5 106 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
107 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
9ca6be09 108
cfc19dd5 109 mult->SetMultiplicityESD(histID, hist);
9ca6be09 110
cfc19dd5 111 //mult->ApplyGaussianMethod(histID, kFALSE);
112 //for (Float_t f=0; f<0.1; f+=0.01)
113 //mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
114 mult->ApplyMinuitFit(histID, kFALSE);
115 mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, hist2->ProjectionY());
116
117 //mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
0a173978 118
119 return mult;
120}
cfc19dd5 121
122const char* GetRegName(Int_t type)
123{
124 switch (type)
125 {
126 case AliMultiplicityCorrection::kNone: return "None"; break;
127 case AliMultiplicityCorrection::kPol0: return "Pol0"; break;
128 case AliMultiplicityCorrection::kPol1: return "Pol1"; break;
129 case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
130 case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break;
131 }
132 return 0;
133}
134
135void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", Int_t histID = 2)
136{
137 gSystem->Load("libPWG0base");
138
139 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
140
141 TFile::Open(fileNameMC);
142 mult->LoadHistograms("Multiplicity");
143
144 TFile::Open(fileNameESD);
145 TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
146 TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
147
148 mult->SetMultiplicityESD(histID, hist);
149
150 TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600);
151 TLegend* legend = new TLegend(0.6, 0.1, 0.98, 0.4);
152 legend->SetFillColor(0);
153
154 Float_t min = 1e10;
155 Float_t max = 0;
156
157 TGraph* first = 0;
158
159 for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
160 {
161 TGraph* fitResultsMC = new TGraph;
162 fitResultsMC->SetTitle(";Weight Parameter");
163 TGraph* fitResultsRes = new TGraph;
164 fitResultsRes->SetTitle(";Weight Parameter");
165
166 fitResultsMC->SetFillColor(0);
167 fitResultsRes->SetFillColor(0);
168 fitResultsMC->SetMarkerStyle(19+type);
169 fitResultsRes->SetMarkerStyle(19+type);
170 fitResultsRes->SetMarkerColor(kRed);
171 fitResultsRes->SetLineColor(kRed);
172
173 legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetRegName(type)));
174 legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetRegName(type)));
175
176 if (first == 0)
177 first = fitResultsMC;
178
179 for (Float_t weight = 1e-2; weight < 2e4; weight *= 1e2)
180 {
181 Float_t chi2MC = 0;
182 Float_t residuals = 0;
183
184 mult->SetRegularizationParameters(type, weight);
185 mult->ApplyMinuitFit(histID, kFALSE);
186 mult->DrawComparison(Form("MinuitChi2_%d_%f", type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
187 mult->GetComparisonResults(chi2MC, residuals);
188
189 fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
190 fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
191
192 min = TMath::Min(min, TMath::Min(chi2MC, residuals));
193 max = TMath::Max(max, TMath::Max(chi2MC, residuals));
194 }
195
196 fitResultsMC->Print();
197 fitResultsRes->Print();
198
199 canvas->cd();
200 fitResultsMC->Draw(Form("%s CP", (type == AliMultiplicityCorrection::kPol0) ? "A" : "SAME"));
201 fitResultsRes->Draw("SAME CP");
202 }
203
204 gPad->SetLogx();
205 gPad->SetLogy();
206 printf("min = %f, max = %f\n", min, max);
207 if (min <= 0)
208 min = 1e-5;
209 first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
210
211 legend->Draw();
212
213 canvas->SaveAs(Form("%s.gif", canvas->GetName()));
214}
215
216void Merge(Int_t n, const char** files, const char* output)
217{
218 gSystem->Load("libPWG0base");
219
220 AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
221 TList list;
222 for (Int_t i=0; i<n; ++i)
223 {
224 TString name("Multiplicity");
225 if (i > 0)
226 name.Form("Multiplicity%d", i);
227
228 TFile::Open(files[i]);
229 data[i] = new AliMultiplicityCorrection(name, name);
230 data[i]->LoadHistograms("Multiplicity");
231 if (i > 0)
232 list.Add(data[i]);
233 }
234
235 data[0]->Merge(&list);
236
237 data[0]->DrawHistograms();
238
239 TFile::Open(output, "RECREATE");
240 data[0]->SaveHistograms();
241 gFile->Close();
242}
243
244void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
245{
246 gSystem->Load("libPWG0base");
247
248 AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
249
250 TFile::Open(fileName);
251 mult->LoadHistograms("Multiplicity");
252
253 TF1* func = 0;
254
255 if (caseNo >= 4)
256 {
257 func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
258 func->SetParNames("scaling", "averagen", "k");
259 }
260
261 switch (caseNo)
262 {
263 case 0: func = new TF1("flat", "1"); break;
264 case 1: func = new TF1("flat", "501-x"); break;
265 case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
266 case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
267 case 4: func->SetParameters(1000, 10, 2); break;
268 case 5: func->SetParameters(1000, 20, 3); break;
269 case 6: func->SetParameters(1000, 30, 4); break;
270 case 7: func->SetParameters(1000, 40, 5); break;
271
272 default: return;
273 }
274
275 mult->SetGenMeasFromFunc(func, 2);
276
277 mult->ApplyBayesianMethod(2, kFALSE);
278 mult->ApplyMinuitFit(2, kFALSE);
279 //mult->ApplyGaussianMethod(2, kFALSE);
280}