CORRFW/test/testUnfolding.C

   1
   2 void testUnfolding() {
   3   TBenchmark b;
   4
   5   b.Start("init");
   6
   7   Load();
   8
   9   AliLog::SetGlobalDebugLevel(0);
  10
  11   // get the essential
  12   AliCFDataGrid *measured    = (AliCFDataGrid*) GetMeasuredSpectrum();
  13   AliCFDataGrid *generated   = (AliCFDataGrid*) GetGeneratedSpectrum();
  14   AliCFEffGrid  *efficiency  = (AliCFEffGrid*)  GetEfficiency();
  15   THnSparse     *response    = (THnSparse*)     GetResponseMatrix();
  16
  17   // create a guessed "a priori" distribution using binning of MC
  18   THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
  19   //----
  20
  21   TF2* fit2D = new TF2("fit2D","[0]*([2]+[3]*y+[4]*y*y)*x*TMath::Exp(-x/[1])",0.1,0.8,-1.5,1.5);
  22   fit2D->SetParameter(0,1);
  23   fit2D->SetParameter(1,1);
  24   fit2D->SetParameter(2,1);
  25   fit2D->SetParameter(3,1);
  26   fit2D->SetParameter(4,1);
  27   fit2D->SetParLimits(1,0,1);
  28
  29   AliCFUnfolding unfolding("unfolding","",2,response,efficiency->GetGrid(),((AliCFGridSparse*)measured->GetData())->GetGrid()/*,guessed*/);
  30   unfolding.SetMaxNumberOfIterations(100);
  31   unfolding.SetMaxChi2PerDOF(1.e-07);
  32   unfolding.UseSmoothing(fit2D,"ren");
  33   //unfolding.UseSmoothing();
  34
  35   b.Stop("init");
  36   b.Start("unfolding");
  37   unfolding.Unfold();
  38   b.Stop("unfolding");
  39   b.Start("finish");
  40
  41   canvas->cd();
  42
  43   TH2D* h_gen = generated->Project(0,1);
  44   h_gen->SetTitle("generated");
  45   h_gen->Draw("lego2");
  46   canvas->SaveAs("/tmp/generated.eps");
  47   canvas->SaveAs("/tmp/generated.gif");
  48
  49   TH2D* h_meas = measured->Project(0,1);
  50   h_meas->SetTitle("measured");
  51   h_meas->Draw("lego2");
  52   canvas->SaveAs("/tmp/measured.eps");
  53   canvas->SaveAs("/tmp/measured.gif");
  54
  55   TH2D* h_guessed = unfolding.GetOriginalPrior()->Projection(1,0);
  56   h_guessed->SetTitle("a priori");
  57   h_guessed->Draw("lego2");
  58   canvas->SaveAs("/tmp/apriori.eps");
  59   canvas->SaveAs("/tmp/apriori.gif");
  60
  61   TH2D* h_eff = efficiency->Project(0,1);
  62   h_eff->SetTitle("efficiency");
  63   h_eff->Draw("e");
  64
  65   TH2D* h_unf = unfolding.GetUnfolded()->Projection(1,0);
  66   h_unf->SetTitle("unfolded");
  67   h_unf->Draw("lego2");
  68   fit2D->Draw("surf same");
  69   return;
  70   canvas->SaveAs("/tmp/unfolded.eps");
  71   canvas->SaveAs("/tmp/unfolded.gif");
  72
  73
  74   TH2D* ratio = (TH2D*)h_unf->Clone();
  75   ratio->SetTitle("unfolded / generated");
  76   ratio->Divide(h_unf,h_gen,1,1);
  77   ratio->GetZaxis()->SetRangeUser(0.5,1.5);
  78   //ratio->Draw("cont4z");
  79   ratio->Draw("lego2");
  80   //ratio->Draw("e");
  81 //   ratio->ProjectionY()->Draw();
  82   canvas->SaveAs("/tmp/ratio.eps");
  83   canvas->SaveAs("/tmp/ratio.gif");
  84
  85
  86   canvas->cd(7);
  87   TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
  88   orig->SetTitle("original prior");
  89   orig->Draw("lego2");
  90
  91   canvas->cd(8);
  92   TH2D* h_est = (TH2D*) unfolding.GetEstMeasured()->Projection(1,0);
  93   h_est->SetTitle("est. measured");
  94   h_est->Draw("e");
  95
  96   canvas->cd(9);
  97   unfolding.GetUnfolded()->Projection(0)->Draw();
  98
  99   canvas->cd(10);
 100   unfolding.GetUnfolded()->Projection(1)->Draw();
 101
 102   canvas->cd(11);
 103   h_gen->ProjectionX()->Draw();
 104
 105   canvas->cd(12);
 106   h_gen->ProjectionY()->Draw();
 107
 108   b.Stop("finish");
 109   Float_t x,y;
 110   b.Summary(x,y);
 111 }
 112
 113 // ====================================================================
 114
 115
 116 void Load(){
 117   gSystem->Load("libANALYSIS");
 118   gSystem->Load("libCORRFW");
 119   gStyle->SetPalette(1);
 120   gStyle->SetOptStat(0);
 121   TCanvas * canvas = new TCanvas("canvas","",600,400);
 122   canvas->SetFillColor(10);
 123   canvas->SetFrameFillColor(10);
 124 }
 125
 126 TObject* GetMeasuredSpectrum() {
 127   TFile * f = TFile::Open("test/output.root","read");
 128   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 129   AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
 130   data->SetMeasured(1);
 131   return data;
 132 }
 133 TObject* GetGeneratedSpectrum() {
 134   TFile * f = TFile::Open("test/output.root","read");
 135   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 136   AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
 137   data->SetMeasured(0);
 138   return data;
 139 }
 140 TObject* GetEfficiency() {
 141   TFile * f = TFile::Open("test/output.root","read");
 142   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 143   AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
 144   eff->CalculateEfficiency(2,0);
 145   return eff;
 146 }
 147 THnSparse* GetResponseMatrix() {
 148   TFile * f = TFile::Open("test/output.root","read");
 149   THnSparse* h = (THnSparse*)f->Get("correlation");
 150   return h;
 151 }
 152 THnSparse* CreateGuessed(const THnSparse* h) {
 153   THnSparse* guessed = (THnSparse*) h->Clone();
 154   return guessed ;
 155 }