2 // script with functions to use AliJetSpectrumUnfolding class
7 // load all the needed libs to run wit root only
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");
26 void draw(const char* fileName = "unfolded_pwg4spec.root", const char* folder = "unfolding", Bool_t proj = kTRUE)
31 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
33 TFile::Open(fileName);
34 jetSpec->LoadHistograms();
39 canvas1 = new TCanvas("Response Map Projection", "Response Map Projection", 500, 500);
43 const Int_t NRGBs = 5;
44 const Int_t NCont = 500;
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);
56 h2 = (jetSpec->GetCorrelation())->Projection(2,3);
57 h2->SetXTitle("z^{lp}_{gen}");
58 h2->SetYTitle("z^{lp}_{rec}");
63 h3 = (jetSpec->GetCorrelation())->Projection(0,1);
64 h3->SetXTitle("E^{jet}_{gen} [GeV]");
65 h3->SetYTitle("E^{jet}_{rec} [GeV]");
68 jetSpec->DrawComparison("Draw_unfolded_pwg4spec", jetSpec->GetGenSpectrum());
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")
76 // function to load jet spectra from the output file of the task AliAnalysisTaskJetSpectrum
77 // to do the unfolding
81 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
83 TFile::Open(fileNameRec);
84 jetSpec->LoadHistograms();
86 TFile::Open(fileNameGen);
87 TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
88 jetSpec->SetGenSpectrum(hist);
90 jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
91 // last parameter = calculateErrors <- this method to calculate the errors takes a lot of time
93 TFile* file = TFile::Open(fileNameUnf,"RECREATE");
94 jetSpec->SaveHistograms();
98 //___________________________________________________________________________
99 void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile)
101 // function to test 2D Bayesian unfolding with other spectra
102 // build from a function
108 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
111 jetSpec->LoadHistograms("unfolding");
113 TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3);
114 h2->DrawCopy("colz");
116 if(getchar()=='q')return;
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;
131 //new TCanvas; func->Draw();
133 jetSpec->SetGenRecFromFunc(func);
135 TFile* file = TFile::Open(outFile,"RECREATE");
136 jetSpec->SaveHistograms();
138 h2 = (jetSpec->GetCorrelation())->Projection(0,3);
139 h2->DrawCopy("colz");
143 //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz");
144 //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz");
147 //_____________________________________________________________________________________________
148 void buildResponseMap(const char* fileName = "responseMap.root")
150 // function to build a Response Map with a gaussian distribution
153 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
155 TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)");
157 bool bPerfect = false;
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++)
164 var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx);
165 var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty);
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++)
172 var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx);
173 var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my);
177 if (var[1]==var[0] && var[3]==var[2])
178 jetSpec->GetCorrelation()->Fill(var,1);
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]));
189 TFile* file = TFile::Open(fileName,"RECREATE");
190 jetSpec->SaveHistograms();
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")
197 // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding
199 gSystem->Load("libANALYSIS");
200 gSystem->Load("libANALYSISalice");
201 gSystem->Load("libJETAN");
202 gSystem->Load("libPWG4JetTasks");
204 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
206 TFile::Open(fileNameRec);
207 jetSpec->LoadHistograms();
209 TFile::Open(fileNameGen);
210 TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
211 jetSpec->SetGenSpectrum(hist);
213 // create sigma histogram
214 TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma");
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();
222 Int_t nIterations = 1000;
224 for(Int_t i=0; i<nIterations; i++)
226 printf("iteration = %d\n", i);
228 jetSpec->GetRecSpectrum()->Reset();
229 jetSpec->GetGenSpectrum()->Reset();
230 jetSpec->GetCorrelation()->Reset();
231 jetSpec->GetUnfSpectrum()->Reset();
233 THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr");
234 TH2F* tmptrue = (TH2F*)htrue->Clone("tmptrue");
236 jetSpec->SetGenSpectrum(tmptrue);
237 jetSpec->SetCorrelation(tmpcorr);
239 // randomize reconstructed distribution
240 for (Int_t me=1; me<=hmeas->GetNbinsX(); me++)
241 for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++)
243 binContent = hmeas->GetBinContent(me,mz);
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();
251 jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent);
255 jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0);
258 for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
259 for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
261 if (htrue->GetBinContent(te,tz)!=0)
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));
272 for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
273 for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
275 binContent = sigma->GetBinContent(te,tz);
276 binContent = TMath::Sqrt(binContent/(Float_t)nIterations);
277 sigma->SetBinContent(te,tz, binContent);
280 const Int_t NRGBs = 5;
281 const Int_t NCont = 500;
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);
293 sigma->SetTitle("#sigma((U_{R} - U)/U)");
294 sigma->SetXTitle("E^{jet} [GeV]");
295 sigma->SetYTitle("z^{lp}");
298 TFile* file = TFile::Open(fileNameOut,"RECREATE");
303 //_______________________________________________________________________________________________________________
304 void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
306 // This functions avoids problems when the number of bins or the bin limits
307 // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
310 gSystem->Load("libANALYSIS");
311 gSystem->Load("libANALYSISalice");
312 gSystem->Load("libJETAN");
313 gSystem->Load("libPWG4JetTasks");
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"));
323 fhCorrelation->Add(fhCorrelation2, 1);
324 fhCorrelation->Add(fhCorrelation3, 1);
326 delete fhCorrelation2;
327 delete fhCorrelation3;
329 AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
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++)
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));
344 file = TFile::Open("gen_pwg4spec.root", "RECREATE");
345 jetSpec->SaveHistograms();
347 jetSpec->GetGenSpectrum()->Reset();
348 printf("true distribution has been set\n");
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++)
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));
365 jetSpec->GetCorrelation()->Reset();
366 jetSpec->GetCorrelation()->Sumw2();
368 for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
370 //printf("idx: %d\n",idx);
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));
388 file = TFile::Open("rec_pwg4spec.root", "RECREATE");
389 jetSpec->SaveHistograms();
392 printf("reconstructed distribution has been set\n");
393 printf("response map has been set\n");
401 // simple steering to correct a given distribution;
403 FillSpecFromFile("/home/ydelgado/pcalice014.cern.ch/macros/jets/CAFdata/histos_pwg4spec.root");
406 sprintf(name, "unfolded_pwg4spec.root");
408 unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name);
409 //draw(name, "unfolding", 1);