Adding the macros and the class for 2-dim unfolding (Ydalia)
[u/mrichter/AliRoot.git] / PWG4 / macros / runJetSpectrumUnfolding.C
1 //
2 // script with functions to use AliJetSpectrumUnfolding class
3 //
4
5
6 void loadlibs(){
7   // load all the needed libs to run wit root only
8
9   gSystem->Load("libTree");
10   gSystem->Load("libPhysics");
11   gSystem->Load("libHist");
12   gSystem->Load("libVMC");
13   gSystem->Load("libSTEERBase");
14   gSystem->Load("libAOD");
15   gSystem->Load("libESD");
16   gSystem->Load("libANALYSIS");
17   gSystem->Load("libANALYSISalice");
18   gSystem->Load("libJETAN");
19   gSystem->Load("libPWG4JetTasks");
20
21
22 }
23
24
25
26 void draw(const char* fileName = "unfolded_pwg4spec.root", const char* folder = "unfolding", Bool_t proj = kTRUE)
27 {
28
29   loadlibs();
30
31   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
32
33   TFile::Open(fileName);
34   jetSpec->LoadHistograms();
35   
36   
37   if (proj)
38   {
39     canvas1 = new TCanvas("Response Map Projection", "Response Map Projection", 500, 500);
40     canvas1->Divide(2);
41   
42     Int_t style = 1;
43     const Int_t NRGBs = 5;
44     const Int_t NCont = 500;
45
46     Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
47     Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
48     Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
49     Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
50     TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
51     gStyle->SetNumberContours(NCont);
52     gStyle->SetPalette(style);
53
54     canvas1->cd(1);
55     gPad->SetLogz();
56     h2 = (jetSpec->GetCorrelation())->Projection(2,3);
57     h2->SetXTitle("z^{lp}_{gen}");
58     h2->SetYTitle("z^{lp}_{rec}");    
59     h2->DrawCopy("colz");
60   
61     canvas1->cd(2);
62     gPad->SetLogz();  
63     h3 = (jetSpec->GetCorrelation())->Projection(0,1);
64     h3->SetXTitle("E^{jet}_{gen} [GeV]");
65     h3->SetYTitle("E^{jet}_{rec} [GeV]");    
66     h3->DrawCopy("colz");
67   }
68   jetSpec->DrawComparison("Draw_unfolded_pwg4spec", jetSpec->GetGenSpectrum());
69
70   return;
71 }
72
73 //________________________________________________________________________________________________________________
74 void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameUnf = "unfolded_pwg4spec.root")
75 {
76   // function to load jet spectra from the output file of the task AliAnalysisTaskJetSpectrum
77   // to do the unfolding
78
79   loadlibs();
80
81   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
82
83   TFile::Open(fileNameRec);
84   jetSpec->LoadHistograms();
85
86   TFile::Open(fileNameGen);
87   TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
88   jetSpec->SetGenSpectrum(hist);
89
90   jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
91   // last parameter = calculateErrors  <- this method to calculate the errors takes a lot of time
92    
93   TFile* file = TFile::Open(fileNameUnf,"RECREATE");
94   jetSpec->SaveHistograms();
95   file->Close();
96 }
97
98 //___________________________________________________________________________
99 void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile)
100 {
101   // function to test 2D Bayesian unfolding with other spectra
102   // build from a function
103
104   loadlibs();
105
106
107
108   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
109
110   TFile::Open(inFile);
111   jetSpec->LoadHistograms("unfolding");
112   
113   TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3);
114   h2->DrawCopy("colz");
115   c1->Update();
116   if(getchar()=='q')return;
117   
118
119   switch (caseNo)
120   {
121     case 1: func = new TF2("func", "501-x-y"); break;
122     case 2: func = new TF2("func", "1000 * 1/(y+x+1)"); break;
123     case 3: func = new TF2("func", "1./(y*pow(x,5.7))"); break;
124     case 4: func = new TF2("func", "exp(-0.1*x - 0.1*y)"); break;
125     case 5: func = new TF2("func", "x*x + (y**3)/100."); break;
126     case 6: func = new TF2("func", "1000*y*exp(-0.1*x)"); break;
127     case 7: func = new TF2("func", "exp(-((x-100.)/(0.3*x))**2 - ((y-0.6)/(0.8*y))**2)"); break;
128     default: return;
129   }
130
131   //new TCanvas; func->Draw();
132
133   jetSpec->SetGenRecFromFunc(func);
134               
135   TFile* file = TFile::Open(outFile,"RECREATE");
136   jetSpec->SaveHistograms();
137   
138   h2 = (jetSpec->GetCorrelation())->Projection(0,3);
139   h2->DrawCopy("colz");
140   c1->Update();
141   file->Close();
142
143   //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz");
144   //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz");
145 }
146
147 //_____________________________________________________________________________________________
148 void buildResponseMap(const char* fileName = "responseMap.root")
149 {
150   // function to build a Response Map with a gaussian distribution
151   loadlibs();
152
153   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
154
155   TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)");
156   
157   bool bPerfect = false;
158   
159   Double_t var[4];
160   Float_t sigmax, sigmay;
161   for (Int_t tx=1; tx<=jetSpec->GetCorrelation()->GetAxis(0)->GetNbins(); tx++)
162     for (Int_t ty=1; ty<=jetSpec->GetCorrelation()->GetAxis(2)->GetNbins(); ty++)
163     {
164       var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx);
165       var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty);
166       sigmax = 0.2*var[0];
167       sigmay = 0.2*var[2];
168       func->SetParameters(var[0],sigmax,var[2],sigmay);
169       for (Int_t mx=1; mx<=jetSpec->GetCorrelation()->GetAxis(1)->GetNbins(); mx++)
170         for (Int_t my=1; my<=jetSpec->GetCorrelation()->GetAxis(3)->GetNbins(); my++)
171         {
172           var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx);
173           var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my);
174
175
176           if(bPerfect){
177             if (var[1]==var[0] && var[3]==var[2])
178               jetSpec->GetCorrelation()->Fill(var,1);
179           }
180           else {
181             // cut at  sigma
182             if (TMath::Abs(var[1]-var[0]) < sigmax || TMath::Abs(var[3]-var[2]) < sigmay)
183               jetSpec->GetCorrelation()->Fill(var,func->Eval(var[1],var[3]));
184           }
185         }
186     }
187
188
189   TFile* file = TFile::Open(fileName,"RECREATE");
190   jetSpec->SaveHistograms();
191   file->Close();
192 }
193
194 //_____________________________________________________________________________
195 void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameOut = "Uncertainties2DBayesUnf.root")
196 {
197   // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding
198
199   gSystem->Load("libANALYSIS");
200   gSystem->Load("libANALYSISalice");
201   gSystem->Load("libJETAN");
202   gSystem->Load("libPWG4JetTasks");
203
204   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
205
206   TFile::Open(fileNameRec);
207   jetSpec->LoadHistograms();
208
209   TFile::Open(fileNameGen);
210   TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
211   jetSpec->SetGenSpectrum(hist);
212     
213   // create sigma histogram
214   TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma");
215   sigma->Reset();
216     
217   THnSparseF* hcorr = (THnSparseF*)(jetSpec->GetCorrelation())->Clone();
218   TH2F*       htrue = (TH2F*)(jetSpec->GetGenSpectrum())->Clone();  
219   TH2F*       hmeas = (TH2F*)(jetSpec->GetRecSpectrum())->Clone();    
220   TH2F*       hunfo = (TH2F*)(jetSpec->GetUnfSpectrum())->Clone();      
221   
222   Int_t nIterations = 1000;
223   Float_t binContent;
224   for(Int_t i=0; i<nIterations; i++)
225   {
226     printf("iteration = %d\n", i);
227     // reset histograms
228     jetSpec->GetRecSpectrum()->Reset();
229     jetSpec->GetGenSpectrum()->Reset();
230     jetSpec->GetCorrelation()->Reset();
231     jetSpec->GetUnfSpectrum()->Reset();
232   
233     THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr"); 
234     TH2F*       tmptrue = (TH2F*)htrue->Clone("tmptrue");  
235     
236     jetSpec->SetGenSpectrum(tmptrue);
237     jetSpec->SetCorrelation(tmpcorr);
238   
239     // randomize reconstructed distribution
240     for (Int_t me=1; me<=hmeas->GetNbinsX(); me++)
241       for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++)
242       {
243         binContent = hmeas->GetBinContent(me,mz);
244         if (binContent>0)
245         {
246           TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",binContent*0.25, binContent*1.25);
247           poisson->SetParameters(binContent,0.);
248           binContent = poisson->GetRandom();
249           delete poisson;
250         }  
251         jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent);
252       } 
253         
254     // unfold
255     jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0);
256     
257     // calculate sigma^2
258     for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
259       for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
260       {
261         if (htrue->GetBinContent(te,tz)!=0)
262         {
263           binContent = (jetSpec->GetUnfSpectrum()->GetBinContent(te,tz) -
264                         htrue->GetBinContent(te,tz))/htrue->GetBinContent(te,tz);
265           binContent *= binContent;
266           sigma->SetBinContent(te,tz, binContent + sigma->GetBinContent(te,tz));
267         }  
268       } 
269   }
270  
271   // calculate sigma   
272   for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
273     for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
274     {
275       binContent = sigma->GetBinContent(te,tz);
276       binContent = TMath::Sqrt(binContent/(Float_t)nIterations);
277       sigma->SetBinContent(te,tz, binContent);
278     } 
279           
280   const Int_t NRGBs = 5;
281   const Int_t NCont = 500;
282
283   Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
284   Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
285   Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
286   Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
287   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
288   gStyle->SetNumberContours(NCont);   
289   gStyle->SetPalette(1);  
290
291   new TCanvas();
292   gPad->SetLogz();
293   sigma->SetTitle("#sigma((U_{R} - U)/U)");
294   sigma->SetXTitle("E^{jet} [GeV]");
295   sigma->SetYTitle("z^{lp}");  
296   sigma->DrawCopy();
297   
298   TFile* file = TFile::Open(fileNameOut,"RECREATE");
299   sigma->Write();
300   file->Close();   
301 }
302
303 //_______________________________________________________________________________________________________________
304 void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
305 {
306   // This functions avoids problems when the number of bins or the bin limits
307   // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
308   // are different.
309
310   gSystem->Load("libANALYSIS");
311   gSystem->Load("libANALYSISalice");
312   gSystem->Load("libJETAN");
313   gSystem->Load("libPWG4JetTasks");
314
315   file = new TFile(fileNameSpec,"read");
316   tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
317   THnSparseF *fhCorrelation  = (THnSparseF*)(tlist->FindObject("fhCorrelation_less5tracks"));
318   THnSparseF *fhCorrelation2 = (THnSparseF*)(tlist->FindObject("fhCorrelation_5to10tracks"));
319   THnSparseF *fhCorrelation3 = (THnSparseF*)(tlist->FindObject("fhCorrelation_more10tracks"));
320   TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fhEGenZGen"));
321   TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fhERecZRec"));  
322
323   fhCorrelation->Add(fhCorrelation2, 1);
324   fhCorrelation->Add(fhCorrelation3, 1);
325
326   delete fhCorrelation2;
327   delete fhCorrelation3;
328
329   AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
330
331   // generated jets (true distribution)
332   for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
333     for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
334     {
335        Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
336        Float_t  z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
337        Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
338        Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
339        Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
340        Float_t err  = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
341        jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
342        jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
343     }
344   file = TFile::Open("gen_pwg4spec.root", "RECREATE");
345   jetSpec->SaveHistograms();
346   file->Close();
347   jetSpec->GetGenSpectrum()->Reset();
348   printf("true distribution has been set\n");
349
350   // reconstructed jets (measured distribution)
351   for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
352     for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
353     {
354        Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
355        Float_t   zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
356        Int_t bine   = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
357        Int_t binz   = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
358        Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
359        Float_t err  = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
360        jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
361        jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
362     }
363
364   // Response function
365   jetSpec->GetCorrelation()->Reset();
366   jetSpec->GetCorrelation()->Sumw2();
367         
368   for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
369   {
370     //printf("idx: %d\n",idx);
371     Double_t var[4];
372     Int_t bin[4];
373     Float_t BinContent = fhCorrelation->GetBinContent(idx, bin);
374     var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]);
375     var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]);
376     var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]);
377     var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]);
378     bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]);
379     bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]);    
380     bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]);
381     bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]);        
382     Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin);
383     Float_t err  = jetSpec->GetCorrelation()->GetBinError(bin);
384     jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
385     jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
386   }
387
388   file = TFile::Open("rec_pwg4spec.root", "RECREATE");
389   jetSpec->SaveHistograms();
390   file->Close();
391
392   printf("reconstructed distribution has been set\n");    
393   printf("response map has been set\n");
394   
395 }
396
397
398
399
400 void correct(){
401   // simple steering to correct a given distribution;
402   loadlibs();
403   FillSpecFromFile("/home/ydelgado/pcalice014.cern.ch/macros/jets/CAFdata/histos_pwg4spec.root");
404
405   char name[100];
406   sprintf(name, "unfolded_pwg4spec.root");
407
408   unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name);
409   //draw(name, "unfolding", 1); 
410
411 }