// // script with functions to use AliJetSpectrumUnfolding class // void loadlibs(){ // load all the needed libs to run wit root only gSystem->Load("libTree"); gSystem->Load("libPhysics"); gSystem->Load("libHist"); gSystem->Load("libVMC"); gSystem->Load("libSTEERBase"); gSystem->Load("libAOD"); gSystem->Load("libESD"); gSystem->Load("libANALYSIS"); gSystem->Load("libANALYSISalice"); gSystem->Load("libJETAN"); gSystem->Load("libPWG4JetTasks"); } void draw(const char* fileName = "unfolded_pwg4spec.root", const char* folder = "unfolding", Bool_t proj = kTRUE) { loadlibs(); AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder); TFile::Open(fileName); jetSpec->LoadHistograms(); if (proj) { canvas1 = new TCanvas("Response Map Projection", "Response Map Projection", 500, 500); canvas1->Divide(2); Int_t style = 1; const Int_t NRGBs = 5; const Int_t NCont = 500; Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); gStyle->SetNumberContours(NCont); gStyle->SetPalette(style); canvas1->cd(1); gPad->SetLogz(); h2 = (jetSpec->GetCorrelation())->Projection(2,3); h2->SetXTitle("z^{lp}_{gen}"); h2->SetYTitle("z^{lp}_{rec}"); h2->DrawCopy("colz"); canvas1->cd(2); gPad->SetLogz(); h3 = (jetSpec->GetCorrelation())->Projection(0,1); h3->SetXTitle("E^{jet}_{gen} [GeV]"); h3->SetYTitle("E^{jet}_{rec} [GeV]"); h3->DrawCopy("colz"); } jetSpec->DrawComparison("Draw_unfolded_pwg4spec", jetSpec->GetGenSpectrum()); return; } //________________________________________________________________________________________________________________ void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameUnf = "unfolded_pwg4spec.root") { // function to load jet spectra from the output file of the task AliAnalysisTaskJetSpectrum // to do the unfolding loadlibs(); AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder); TFile::Open(fileNameGen); jetSpec->LoadHistograms(); TFile::Open(fileNameRec); TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum"); jetSpec->SetRecSpectrum(hist); hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum"); if(hist->GetEntries()>0)jetSpec->SetGenSpectrum(hist); jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0); // last parameter = calculateErrors <- this method to calculate the errors takes a lot of time TFile* file = TFile::Open(fileNameUnf,"RECREATE"); jetSpec->SaveHistograms(); file->Close(); } //___________________________________________________________________________ void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile) { // function to test 2D Bayesian unfolding with other spectra // build from a function loadlibs(); AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding"); TFile::Open(inFile); jetSpec->LoadHistograms("unfolding"); TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3); h2->DrawCopy("colz"); c1->Update(); if(getchar()=='q')return; switch (caseNo) { case 1: func = new TF2("func", "501-x-y"); break; case 2: func = new TF2("func", "1000 * 1/(y+x+1)"); break; case 3: func = new TF2("func", "1./(y*pow(x,5.7))"); break; case 4: func = new TF2("func", "exp(-0.1*x - 0.1*y)"); break; case 5: func = new TF2("func", "x*x + (y**3)/100."); break; case 6: func = new TF2("func", "1000*y*exp(-0.1*x)"); break; case 7: func = new TF2("func", "exp(-((x-100.)/(0.3*x))**2 - ((y-0.6)/(0.8*y))**2)"); break; default: return; } //new TCanvas; func->Draw(); jetSpec->SetGenRecFromFunc(func); TFile* file = TFile::Open(outFile,"RECREATE"); jetSpec->SaveHistograms(); h2 = (jetSpec->GetCorrelation())->Projection(0,3); h2->DrawCopy("colz"); c1->Update(); file->Close(); //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz"); //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz"); } //_____________________________________________________________________________________________ void buildResponseMap(const char* fileName = "responseMap.root") { // function to build a Response Map with a gaussian distribution loadlibs(); AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding"); TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)"); bool bPerfect = false; Double_t var[4]; Float_t sigmax, sigmay; for (Int_t tx=1; tx<=jetSpec->GetCorrelation()->GetAxis(0)->GetNbins(); tx++) for (Int_t ty=1; ty<=jetSpec->GetCorrelation()->GetAxis(2)->GetNbins(); ty++) { var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx); var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty); sigmax = 0.2*var[0]; sigmay = 0.2*var[2]; func->SetParameters(var[0],sigmax,var[2],sigmay); for (Int_t mx=1; mx<=jetSpec->GetCorrelation()->GetAxis(1)->GetNbins(); mx++) for (Int_t my=1; my<=jetSpec->GetCorrelation()->GetAxis(3)->GetNbins(); my++) { var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx); var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my); if(bPerfect){ if (var[1]==var[0] && var[3]==var[2]) jetSpec->GetCorrelation()->Fill(var,1); } else { // cut at sigma if (TMath::Abs(var[1]-var[0]) < sigmax || TMath::Abs(var[3]-var[2]) < sigmay) jetSpec->GetCorrelation()->Fill(var,func->Eval(var[1],var[3])); } } } TFile* file = TFile::Open(fileName,"RECREATE"); jetSpec->SaveHistograms(); file->Close(); } //_____________________________________________________________________________ void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameOut = "Uncertainties2DBayesUnf.root") { // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding gSystem->Load("libANALYSIS"); gSystem->Load("libANALYSISalice"); gSystem->Load("libJETAN"); gSystem->Load("libPWG4JetTasks"); AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder); TFile::Open(fileNameRec); jetSpec->LoadHistograms(); TFile::Open(fileNameGen); TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum"); jetSpec->SetGenSpectrum(hist); // create sigma histogram TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma"); sigma->Reset(); THnSparseF* hcorr = (THnSparseF*)(jetSpec->GetCorrelation())->Clone(); TH2F* htrue = (TH2F*)(jetSpec->GetGenSpectrum())->Clone(); TH2F* hmeas = (TH2F*)(jetSpec->GetRecSpectrum())->Clone(); TH2F* hunfo = (TH2F*)(jetSpec->GetUnfSpectrum())->Clone(); Int_t nIterations = 1000; Float_t binContent; for(Int_t i=0; iGetRecSpectrum()->Reset(); jetSpec->GetGenSpectrum()->Reset(); jetSpec->GetCorrelation()->Reset(); jetSpec->GetUnfSpectrum()->Reset(); THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr"); TH2F* tmptrue = (TH2F*)htrue->Clone("tmptrue"); jetSpec->SetGenSpectrum(tmptrue); jetSpec->SetCorrelation(tmpcorr); // randomize reconstructed distribution for (Int_t me=1; me<=hmeas->GetNbinsX(); me++) for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++) { binContent = hmeas->GetBinContent(me,mz); if (binContent>0) { TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",binContent*0.25, binContent*1.25); poisson->SetParameters(binContent,0.); binContent = poisson->GetRandom(); delete poisson; } jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent); } // unfold jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0); // calculate sigma^2 for (Int_t te=1; te<=sigma->GetNbinsX(); te++) for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++) { if (htrue->GetBinContent(te,tz)!=0) { binContent = (jetSpec->GetUnfSpectrum()->GetBinContent(te,tz) - htrue->GetBinContent(te,tz))/htrue->GetBinContent(te,tz); binContent *= binContent; sigma->SetBinContent(te,tz, binContent + sigma->GetBinContent(te,tz)); } } } // calculate sigma for (Int_t te=1; te<=sigma->GetNbinsX(); te++) for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++) { binContent = sigma->GetBinContent(te,tz); binContent = TMath::Sqrt(binContent/(Float_t)nIterations); sigma->SetBinContent(te,tz, binContent); } const Int_t NRGBs = 5; const Int_t NCont = 500; Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); gStyle->SetNumberContours(NCont); gStyle->SetPalette(1); new TCanvas(); gPad->SetLogz(); sigma->SetTitle("#sigma((U_{R} - U)/U)"); sigma->SetXTitle("E^{jet} [GeV]"); sigma->SetYTitle("z^{lp}"); sigma->DrawCopy(); TFile* file = TFile::Open(fileNameOut,"RECREATE"); sigma->Write(); file->Close(); } //_______________________________________________________________________________________________________________ void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const char* fileNameSim = "histos_pwg4spec.root") { // This functions avoids problems when the number of bins or the bin limits // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes // are different. gSystem->Load("libANALYSIS"); gSystem->Load("libANALYSISalice"); gSystem->Load("libJETAN"); gSystem->Load("libPWG4JetTasks"); file = new TFile(fileNameSim,"read"); tlist = dynamic_cast (file->Get("pwg4spec")); THnSparseF *fhCorrelation = 0; for(int ic = 0;ic<3;++ic){ THnSparseF *hTmp = (THnSparseF*)(tlist->FindObject(Form("fhnCorrelation_%d",ic))); if(fhCorrelation==0)fhCorrelation = (THnSparseF*)hTmp->Clone("fhnCorrelation"); else fhCorrelation->Add(hTmp); } TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen")); AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding"); // generated jets (true distribution) for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++) for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++) { Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te); Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz); Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej); Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z); Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz); Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz); // merging of bins happens here! jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz)); jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz)); } Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins()); for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++) { //printf("idx: %d\n",idx); Double_t var[4]; Int_t bin[4]; Float_t BinContent = fhCorrelation->GetBinContent(idx, bin); var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]); var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]); var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]); var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]); bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]); bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]); bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]); bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]); Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin); Float_t err = jetSpec->GetCorrelation()->GetBinError(bin); // merging of bins happens here! jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx)); jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx)); } // need number of entries for correct drawing jetSpec->GetCorrelation()->SetEntries(fhCorrelation->GetEntries()); file = TFile::Open("gen_pwg4spec.root", "RECREATE"); jetSpec->SaveHistograms(); Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins()); file->Close(); jetSpec->GetGenSpectrum()->Reset(); printf("true distribution has been set\n"); // reconstructed jets (measured distribution) // Response function jetSpec->GetCorrelation()->Reset(); jetSpec->GetCorrelation()->Sumw2(); jetSpec->GetGenSpectrum()->Reset(); jetSpec->GetRecSpectrum()->Reset(); file = new TFile(fileNameReal,"read"); tlist = dynamic_cast (file->Get("pwg4spec")); TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fh2ERecZRec")); for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++) for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++) { Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me); Float_t zm = fhERecZRec->GetYaxis()->GetBinCenter(mz); Int_t bine = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec); Int_t binz = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm); Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz); Float_t err = jetSpec->GetRecSpectrum()->GetBinError(bine, binz); jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz)); jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz)); } // for control again, but now from the rec file // generated jets (true distribution) TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen")); for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++) for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++) { Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te); Float_t z = fhEGenZGen->GetYaxis()->GetBinCenter(tz); Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej); Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z); Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz); Float_t err = jetSpec->GetGenSpectrum()->GetBinError(bine,binz); jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz)); jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz)); } file = TFile::Open("rec_pwg4spec.root", "RECREATE"); jetSpec->SaveHistograms(); file->Close(); printf("reconstructed distribution has been set\n"); printf("response map has been set\n"); } void correct(){ // simple steering to correct a given distribution; loadlibs(); // rec and gen // FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root"); FillSpecFromFiles("pwg4spec_allpt_16.root","pwg4spec_allpt_16.root"); char name[100]; sprintf(name, "unfolded_pwg4spec_16.root"); unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name); //draw(name, "unfolding", 1); } void mergeJetAnaOutput(){ // This is used to merge the analysis-output from different // data samples/pt_hard bins // in case the eventweigth was set to xsection/ntrials already, this // is not needed. Both methods only work in case we do not mix different // pt_hard bins, and do not have overlapping bins const Int_t nBins = 2; // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05 // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03 // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02 // const Float_t xsection[nBins] = {2.168291,2.44068e-02}; Float_t xsection[nBins] = {0,0}; Float_t nTrials[nBins] = {0,0}; Float_t sf[nBins] = {0,0}; const char *cFile[nBins] = {"pwg4spec_0000000-1000000_LHC08v_jetjet15-50.root", "pwg4spec_0000000-1000000_LHC08r_jetjet50.root"}; TList *lIn[nBins]; TFile *fIn[nBins]; for(int ib = 0;ib < nBins;++ib){ fIn[ib] = TFile::Open(cFile[ib]); lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec"); TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials"); TProfile* hXsec = (TProfile*)lIn[ib]->FindObject("fh1Xsec"); xsection[ib] = hXsec->GetBinContent(1); nTrials[ib] = hTrials->Integral(); sf[ib] = xsection[ib]/ nTrials[ib]; } TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE"); TList *lOut = new TList(); lOut->SetName(lIn[0]->GetName()); // for the start scale all... for(int ie = 0; ie < lIn[0]->GetEntries();++ie){ TH1 *h1Add = 0; THnSparse *hnAdd = 0; Printf("%d: %s",ie, lIn[0]->At(ie)->GetName()); for(int ib = 0;ib < nBins;++ib){ // dynamic cast does not work with cint TObject *h = lIn[ib]->At(ie); if(h->InheritsFrom("TH1")){ Printf("ib %d",ib); TH1 *h1 = (TH1*)h; if(ib==0){ h1Add = (TH1*)h1->Clone(h1->GetName()); h1Add->Scale(sf[ib]); } else{ h1Add->Add(h1,sf[ib]); } } else if(h->InheritsFrom("THnSparse")){ Printf("ib %d",ib); THnSparse *hn = (THnSparse*)h; if(ib==0){ hnAdd = (THnSparse*)hn->Clone(hn->GetName()); hnAdd->Scale(sf[ib]); } else{ hnAdd->Add(hn,sf[ib]); } } }// ib if(h1Add)lOut->Add(h1Add); else if(hnAdd)lOut->Add(hnAdd); } fOut->cd(); lOut->Write(lOut->GetName(),TObject::kSingleKey); fOut->Close(); }