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 }