adding calculation covariance matrix to bayesian method
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / runMultiplicitySelector.C
1 /* $Id$ */
2
3 //
4 // script to run the AliMultiplicityESDSelector
5 //
6
7 #include "../CreateESDChain.C"
8 #include "../PWG0Helper.C"
9
10 void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "")
11 {
12   if (aProof)
13     connectProof("jgrosseo@lxb6046");
14
15   TString libraries("libEG;libGeom;libESD;libPWG0base");
16   TString packages("PWG0base");
17
18   if (aMC != kFALSE)
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;
26
27   gProof->SetParallel(0);
28   gProof->SetLogLevel(4);
29   gProof->SetParallel(9999);
30
31   gROOT->ProcessLine(".L CreateCuts.C");
32   gROOT->ProcessLine(".L drawPlots.C");
33
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
42   TList inputList;
43   inputList.Add(esdTrackCuts);
44
45   TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE);
46
47   TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector");
48   AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
49
50   selectorName += ".cxx+";
51
52   if (aDebug != kFALSE)
53     selectorName += "g";
54
55   Int_t result = executeQuery(chain, &inputList, selectorName, option);
56
57   if (result != 0)
58   {
59     printf("ERROR: Executing process failed with %d.\n", result);
60     return;
61   }
62 }
63
64 void 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();
74 }
75
76 void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
77 {
78   gSystem->Load("libPWG0base");
79
80   AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
81
82   TFile::Open(fileName);
83   mult->LoadHistograms("Multiplicity");
84
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);
92
93   return mult;
94 }
95
96 void* fitOther(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 2)
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);
106   TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
107   TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
108
109   mult->SetMultiplicityESD(histID, hist);
110
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);
118
119   return mult;
120 }
121
122 const 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
135 void 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
216 void 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
244 void 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 }