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(fileNameGen);
84 jetSpec->LoadHistograms();
86 TFile::Open(fileNameRec);
87 TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum");
88 jetSpec->SetRecSpectrum(hist);
89 hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
90 if(hist->GetEntries()>0)jetSpec->SetGenSpectrum(hist);
92 jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
93 // last parameter = calculateErrors <- this method to calculate the errors takes a lot of time
95 TFile* file = TFile::Open(fileNameUnf,"RECREATE");
96 jetSpec->SaveHistograms();
100 //___________________________________________________________________________
101 void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile)
103 // function to test 2D Bayesian unfolding with other spectra
104 // build from a function
110 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
113 jetSpec->LoadHistograms("unfolding");
115 TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3);
116 h2->DrawCopy("colz");
118 if(getchar()=='q')return;
123 case 1: func = new TF2("func", "501-x-y"); break;
124 case 2: func = new TF2("func", "1000 * 1/(y+x+1)"); break;
125 case 3: func = new TF2("func", "1./(y*pow(x,5.7))"); break;
126 case 4: func = new TF2("func", "exp(-0.1*x - 0.1*y)"); break;
127 case 5: func = new TF2("func", "x*x + (y**3)/100."); break;
128 case 6: func = new TF2("func", "1000*y*exp(-0.1*x)"); break;
129 case 7: func = new TF2("func", "exp(-((x-100.)/(0.3*x))**2 - ((y-0.6)/(0.8*y))**2)"); break;
133 //new TCanvas; func->Draw();
135 jetSpec->SetGenRecFromFunc(func);
137 TFile* file = TFile::Open(outFile,"RECREATE");
138 jetSpec->SaveHistograms();
140 h2 = (jetSpec->GetCorrelation())->Projection(0,3);
141 h2->DrawCopy("colz");
145 //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz");
146 //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz");
149 //_____________________________________________________________________________________________
150 void buildResponseMap(const char* fileName = "responseMap.root")
152 // function to build a Response Map with a gaussian distribution
155 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
157 TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)");
159 bool bPerfect = false;
162 Float_t sigmax, sigmay;
163 for (Int_t tx=1; tx<=jetSpec->GetCorrelation()->GetAxis(0)->GetNbins(); tx++)
164 for (Int_t ty=1; ty<=jetSpec->GetCorrelation()->GetAxis(2)->GetNbins(); ty++)
166 var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx);
167 var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty);
170 func->SetParameters(var[0],sigmax,var[2],sigmay);
171 for (Int_t mx=1; mx<=jetSpec->GetCorrelation()->GetAxis(1)->GetNbins(); mx++)
172 for (Int_t my=1; my<=jetSpec->GetCorrelation()->GetAxis(3)->GetNbins(); my++)
174 var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx);
175 var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my);
179 if (var[1]==var[0] && var[3]==var[2])
180 jetSpec->GetCorrelation()->Fill(var,1);
184 if (TMath::Abs(var[1]-var[0]) < sigmax || TMath::Abs(var[3]-var[2]) < sigmay)
185 jetSpec->GetCorrelation()->Fill(var,func->Eval(var[1],var[3]));
191 TFile* file = TFile::Open(fileName,"RECREATE");
192 jetSpec->SaveHistograms();
196 //_____________________________________________________________________________
197 void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameOut = "Uncertainties2DBayesUnf.root")
199 // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding
201 gSystem->Load("libANALYSIS");
202 gSystem->Load("libANALYSISalice");
203 gSystem->Load("libJETAN");
204 gSystem->Load("libPWG4JetTasks");
206 AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
208 TFile::Open(fileNameRec);
209 jetSpec->LoadHistograms();
211 TFile::Open(fileNameGen);
212 TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
213 jetSpec->SetGenSpectrum(hist);
215 // create sigma histogram
216 TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma");
219 THnSparseF* hcorr = (THnSparseF*)(jetSpec->GetCorrelation())->Clone();
220 TH2F* htrue = (TH2F*)(jetSpec->GetGenSpectrum())->Clone();
221 TH2F* hmeas = (TH2F*)(jetSpec->GetRecSpectrum())->Clone();
222 TH2F* hunfo = (TH2F*)(jetSpec->GetUnfSpectrum())->Clone();
224 Int_t nIterations = 1000;
226 for(Int_t i=0; i<nIterations; i++)
228 printf("iteration = %d\n", i);
230 jetSpec->GetRecSpectrum()->Reset();
231 jetSpec->GetGenSpectrum()->Reset();
232 jetSpec->GetCorrelation()->Reset();
233 jetSpec->GetUnfSpectrum()->Reset();
235 THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr");
236 TH2F* tmptrue = (TH2F*)htrue->Clone("tmptrue");
238 jetSpec->SetGenSpectrum(tmptrue);
239 jetSpec->SetCorrelation(tmpcorr);
241 // randomize reconstructed distribution
242 for (Int_t me=1; me<=hmeas->GetNbinsX(); me++)
243 for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++)
245 binContent = hmeas->GetBinContent(me,mz);
248 TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",binContent*0.25, binContent*1.25);
249 poisson->SetParameters(binContent,0.);
250 binContent = poisson->GetRandom();
253 jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent);
257 jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0);
260 for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
261 for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
263 if (htrue->GetBinContent(te,tz)!=0)
265 binContent = (jetSpec->GetUnfSpectrum()->GetBinContent(te,tz) -
266 htrue->GetBinContent(te,tz))/htrue->GetBinContent(te,tz);
267 binContent *= binContent;
268 sigma->SetBinContent(te,tz, binContent + sigma->GetBinContent(te,tz));
274 for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
275 for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
277 binContent = sigma->GetBinContent(te,tz);
278 binContent = TMath::Sqrt(binContent/(Float_t)nIterations);
279 sigma->SetBinContent(te,tz, binContent);
282 const Int_t NRGBs = 5;
283 const Int_t NCont = 500;
285 Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
286 Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
287 Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
288 Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
289 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
290 gStyle->SetNumberContours(NCont);
291 gStyle->SetPalette(1);
295 sigma->SetTitle("#sigma((U_{R} - U)/U)");
296 sigma->SetXTitle("E^{jet} [GeV]");
297 sigma->SetYTitle("z^{lp}");
300 TFile* file = TFile::Open(fileNameOut,"RECREATE");
305 //_______________________________________________________________________________________________________________
306 void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const char* fileNameSim = "histos_pwg4spec.root")
308 // This functions avoids problems when the number of bins or the bin limits
309 // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
312 gSystem->Load("libANALYSIS");
313 gSystem->Load("libANALYSISalice");
314 gSystem->Load("libJETAN");
315 gSystem->Load("libPWG4JetTasks");
317 file = new TFile(fileNameSim,"read");
318 tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
320 THnSparseF *fhCorrelation = 0;
321 for(int ic = 0;ic<3;++ic){
322 THnSparseF *hTmp = (THnSparseF*)(tlist->FindObject(Form("fhnCorrelation_%d",ic)));
323 if(fhCorrelation==0)fhCorrelation = (THnSparseF*)hTmp->Clone("fhnCorrelation");
324 else fhCorrelation->Add(hTmp);
326 TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
328 AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
330 // generated jets (true distribution)
331 for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
332 for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
334 Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
335 Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
336 Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
337 Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
338 Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
339 Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
340 // merging of bins happens here!
341 jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
342 jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
346 Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
347 for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
349 //printf("idx: %d\n",idx);
352 Float_t BinContent = fhCorrelation->GetBinContent(idx, bin);
353 var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]);
354 var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]);
355 var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]);
356 var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]);
357 bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]);
358 bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]);
359 bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]);
360 bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]);
361 Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin);
362 Float_t err = jetSpec->GetCorrelation()->GetBinError(bin);
363 // merging of bins happens here!
364 jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
365 jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
367 // need number of entries for correct drawing
368 jetSpec->GetCorrelation()->SetEntries(fhCorrelation->GetEntries());
371 file = TFile::Open("gen_pwg4spec.root", "RECREATE");
372 jetSpec->SaveHistograms();
373 Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
375 jetSpec->GetGenSpectrum()->Reset();
376 printf("true distribution has been set\n");
378 // reconstructed jets (measured distribution)
382 jetSpec->GetCorrelation()->Reset();
383 jetSpec->GetCorrelation()->Sumw2();
384 jetSpec->GetGenSpectrum()->Reset();
385 jetSpec->GetRecSpectrum()->Reset();
388 file = new TFile(fileNameReal,"read");
389 tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
392 TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fh2ERecZRec"));
394 for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
395 for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
397 Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
398 Float_t zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
399 Int_t bine = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
400 Int_t binz = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
401 Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
402 Float_t err = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
403 jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
404 jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
408 // for control again, but now from the rec file
409 // generated jets (true distribution)
410 TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
411 for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
412 for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
414 Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
415 Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
416 Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
417 Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
418 Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
419 Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
420 jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
421 jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
425 file = TFile::Open("rec_pwg4spec.root", "RECREATE");
426 jetSpec->SaveHistograms();
429 printf("reconstructed distribution has been set\n");
430 printf("response map has been set\n");
435 // simple steering to correct a given distribution;
438 // FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root");
439 FillSpecFromFiles("pwg4spec_allpt_16.root","pwg4spec_allpt_16.root");
442 sprintf(name, "unfolded_pwg4spec_16.root");
444 unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name);
445 //draw(name, "unfolding", 1);
449 void mergeJetAnaOutput(){
450 // This is used to merge the analysis-output from different
451 // data samples/pt_hard bins
452 // in case the eventweigth was set to xsection/ntrials already, this
453 // is not needed. Both methods only work in case we do not mix different
454 // pt_hard bins, and do not have overlapping bins
456 const Int_t nBins = 2;
457 // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
458 // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
459 // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
460 // const Float_t xsection[nBins] = {2.168291,2.44068e-02};
461 Float_t xsection[nBins] = {0,0};
462 Float_t nTrials[nBins] = {0,0};
463 Float_t sf[nBins] = {0,0};
465 const char *cFile[nBins] = {"pwg4spec_0000000-1000000_LHC08v_jetjet15-50.root",
466 "pwg4spec_0000000-1000000_LHC08r_jetjet50.root"};
471 for(int ib = 0;ib < nBins;++ib){
472 fIn[ib] = TFile::Open(cFile[ib]);
473 lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec");
474 TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials");
475 TProfile* hXsec = (TProfile*)lIn[ib]->FindObject("fh1Xsec");
476 xsection[ib] = hXsec->GetBinContent(1);
477 nTrials[ib] = hTrials->Integral();
478 sf[ib] = xsection[ib]/ nTrials[ib];
481 TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE");
482 TList *lOut = new TList();
483 lOut->SetName(lIn[0]->GetName());
484 // for the start scale all...
485 for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
487 THnSparse *hnAdd = 0;
488 Printf("%d: %s",ie, lIn[0]->At(ie)->GetName());
489 for(int ib = 0;ib < nBins;++ib){
490 // dynamic cast does not work with cint
491 TObject *h = lIn[ib]->At(ie);
492 if(h->InheritsFrom("TH1")){
496 h1Add = (TH1*)h1->Clone(h1->GetName());
497 h1Add->Scale(sf[ib]);
500 h1Add->Add(h1,sf[ib]);
503 else if(h->InheritsFrom("THnSparse")){
505 THnSparse *hn = (THnSparse*)h;
507 hnAdd = (THnSparse*)hn->Clone(hn->GetName());
508 hnAdd->Scale(sf[ib]);
511 hnAdd->Add(hn,sf[ib]);
517 if(h1Add)lOut->Add(h1Add);
518 else if(hnAdd)lOut->Add(hnAdd);
521 lOut->Write(lOut->GetName(),TObject::kSingleKey);