CORRFW/test/testUnfolding.C

   1
   2 void testUnfolding() {
   3   Load();
   4
   5   // get the essential
   6   AliCFDataGrid *measured    = (AliCFDataGrid*) GetMeasuredSpectrum();
   7   AliCFDataGrid *generated   = (AliCFDataGrid*) GetGeneratedSpectrum();
   8   AliCFEffGrid  *efficiency  = (AliCFEffGrid*)  GetEfficiency();
   9   THnSparse     *response    = (THnSparse*)     GetResponseMatrix();
  10
  11   // create a guessed "a priori" distribution using binning of MC
  12   THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
  13   //----
  14   AliCFUnfolding unfolding("unfolding","",2,response,efficiency->GetGrid(),((AliCFGridSparse*)measured->GetData())->GetGrid(),guessed);
  15   unfolding.SetMaxNumberOfIterations(100);
  16   unfolding.SetMaxChi2PerDOF(1.e-07);
  17   //unfolding.UseSmoothing();
  18   unfolding.Unfold();
  19
  20   THnSparse* result = unfolding.GetUnfolded();
  21   //----
  22
  23   TCanvas * canvas = new TCanvas("canvas","",1000,700);
  24   canvas->Divide(3,3);
  25
  26   canvas->cd(1);
  27   TH2D* h_gen = generated->Project(0,1);
  28   h_gen->SetTitle("generated");
  29   h_gen->Draw("lego2");
  30
  31   canvas->cd(2);
  32   TH2D* h_meas = measured->Project(0,1);
  33   h_meas->SetTitle("measured");
  34   h_meas->Draw("lego2");
  35
  36   canvas->cd(3);
  37   TH2D* h_guessed = guessed->Projection(1,0);
  38   h_guessed->SetTitle("a priori");
  39   h_guessed->Draw("lego2");
  40
  41   canvas->cd(4);
  42   TH2D* h_eff = efficiency->Project(0,1);
  43   h_eff->SetTitle("efficiency");
  44   h_eff->Draw("lego2");
  45
  46   canvas->cd(5);
  47   TH2D* h_unf = result->Projection(1,0);
  48   h_unf->SetTitle("unfolded");
  49   h_unf->Draw("lego2");
  50
  51   canvas->cd(6);
  52   TH2D* ratio = (TH2D*)h_unf->Clone();
  53   ratio->SetTitle("unfolded / generated");
  54   ratio->Divide(h_unf,h_gen,1,1);
  55 //   ratio->Draw("cont4z");
  56 //   ratio->Draw("surf2");
  57   ratio->Draw("lego2");
  58
  59   canvas->cd(7);
  60   TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
  61   orig->SetTitle("original prior");
  62   orig->Draw("lego2");
  63
  64   canvas->cd(8);
  65   AliCFDataGrid* corrected = (AliCFDataGrid*)measured->Clone();
  66   corrected->ApplyEffCorrection(*efficiency);
  67   TH2D* corr = corrected->Project(0,1);
  68   corr->SetTitle("simple correction");
  69   corr->Draw("lego2");
  70
  71   canvas->cd(9);
  72   TH2D* ratio2 = (TH2D*) corr->Clone();
  73   ratio2->Divide(corr,h_gen,1,1);
  74   ratio2->SetTitle("simple correction / generated");
  75   ratio2->Draw("cont4z");
  76
  77   return;
  78 }
  79
  80 // ====================================================================
  81
  82
  83 void Load(){
  84   gSystem->Load("libANALYSIS");
  85   gSystem->Load("libCORRFW");
  86   gStyle->SetPalette(1);
  87   gStyle->SetOptStat(0);
  88 }
  89
  90 TObject* GetMeasuredSpectrum() {
  91   TFile * f = TFile::Open("test/output.root","read");
  92   AliCFContainer* c = (AliCFContainer*)f->Get("container");
  93   AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
  94   data->SetMeasured(1);
  95   return data;
  96 }
  97 TObject* GetGeneratedSpectrum() {
  98   TFile * f = TFile::Open("test/output.root","read");
  99   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 100   AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
 101   data->SetMeasured(0);
 102   return data;
 103 }
 104 TObject* GetEfficiency() {
 105   TFile * f = TFile::Open("test/output.root","read");
 106   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 107   AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
 108   eff->CalculateEfficiency(2,0);
 109   return eff;
 110 }
 111 THnSparse* GetResponseMatrix() {
 112   TFile * f = TFile::Open("test/output.root","read");
 113   THnSparse* h = (THnSparse*)f->Get("correlation");
 114   return h;
 115 }
 116 THnSparse* CreateGuessed(const THnSparse* h) {
 117   THnSparse* guessed = (THnSparse*) h->Clone();
 118   return guessed ;
 119 }