PWG3/hfe/testUnfolding.C

   1 void Load();
   2 TList *GetResults(const Char_t *testfile);
   3 TObject* GetMeasuredSpectrum(TList *fContainer);
   4 TObject* GetGeneratedSpectrum(TList *fContainer);
   5 TObject* GetEfficiency(TList *fContainer);
   6 THnSparse* GetResponseMatrix(TList *fContainer);
   7 THnSparse* CreateGuessed(const THnSparse* h);
   8
   9 void testUnfolding(const char *testfile = "HFEtask.root") {
  10   Load();
  11   TList *results = GetResults(testfile);
  12   if(!results){
  13     Error("No output objects: Calculation will terminate here");
  14     return;
  15   }
  16   // get the essential
  17   AliCFDataGrid *measured    = (AliCFDataGrid*) GetMeasuredSpectrum(testfile);
  18   AliCFDataGrid *generated   = (AliCFDataGrid*) GetGeneratedSpectrum(testfile);
  19   AliCFEffGrid  *efficiency  = (AliCFEffGrid*)  GetEfficiency(testfile);
  20   THnSparse     *response    = (THnSparse*)     GetResponseMatrix(testfile);
  21
  22   // create a guessed "a priori" distribution using binning of MC
  23   THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
  24   //----
  25   AliCFUnfolding unfolding("unfolding","",3,response,efficiency->GetGrid(),((AliCFGridSparse*)measured->GetData())->GetGrid(),guessed);
  26   unfolding.SetMaxNumberOfIterations(100);
  27   unfolding.SetMaxChi2PerDOF(1.e-07);
  28   //unfolding.UseSmoothing();
  29   unfolding.Unfold();
  30
  31   THnSparse* result = unfolding.GetUnfolded();
  32   //----
  33
  34   TCanvas * canvas = new TCanvas("canvas","",1000,700);
  35   canvas->Divide(3,3);
  36
  37   canvas->cd(1);
  38   TH2D* h_gen = generated->Project(0,1);
  39   h_gen->SetTitle("generated");
  40   //printf("c1\n");
  41   h_gen->Draw("lego2");
  42
  43   canvas->cd(2);
  44   TH2D* h_meas = measured->Project(0,1);
  45   h_meas->SetTitle("measured");
  46   //printf("c2\n");
  47   h_meas->Draw("lego2");
  48
  49   canvas->cd(3);
  50   TH2D* h_guessed = guessed->Projection(1,0);
  51   h_guessed->SetTitle("a priori");
  52   //printf("c3\n");
  53   h_guessed->Draw("lego2");
  54
  55   canvas->cd(4);
  56   TH2D* h_eff = efficiency->Project(0,1);
  57   h_eff->SetTitle("efficiency");
  58   //printf("c4\n");
  59   h_eff->Draw("lego2");
  60
  61   canvas->cd(5);
  62   TH2D* h_unf = result->Projection(1,0);
  63   h_unf->SetTitle("unfolded");
  64   //printf("c5\n");
  65   h_unf->Draw("lego2");
  66
  67   canvas->cd(6);
  68   TH2D* ratio = (TH2D*)h_unf->Clone();
  69   ratio->SetTitle("unfolded / generated");
  70   ratio->Divide(h_unf,h_gen,1,1);
  71 //   ratio->Draw("cont4z");
  72 //   ratio->Draw("surf2");
  73   //printf("c6\n");
  74   ratio->Draw("lego2");
  75
  76   canvas->cd(7);
  77   TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
  78   orig->SetTitle("original prior");
  79   //printf("c7\n");
  80   orig->Draw("lego2");
  81
  82   canvas->cd(8);
  83   AliCFDataGrid* corrected = (AliCFDataGrid*)measured->Clone();
  84   //corrected->ApplyEffCorrection(*efficiency);
  85   TH2D* corr = corrected->Project(0,1);
  86   corr->SetTitle("simple correction");
  87   //printf("c8\n");
  88   corr->Draw("lego2");
  89
  90   canvas->cd(9);
  91   TH2D* ratio2 = (TH2D*) corr->Clone();
  92   ratio2->Divide(corr,h_gen,1,1);
  93   ratio2->SetTitle("simple correction / generated");
  94   //printf("c9\n");
  95   ratio2->Draw("lego2");
  96
  97   return;
  98 }
  99
 100 // ====================================================================
 101
 102
 103 void Load(){
 104   gSystem->Load("libANALYSIS");
 105   gSystem->Load("libCORRFW");
 106   gStyle->SetPalette(1);
 107   gStyle->SetOptStat(0);
 108 }
 109
 110 TList *GetResults(const Char_t *testfile){
 111   //
 112   // read output
 113   //
 114   TFile *f = TFile::Open(testfile);
 115   if(!f || f->IsZombie()){
 116     Error("File not readable");
 117     return 0x0;
 118   }
 119   TList *l = dynamic_cast<TList *>(f->Get("Results"));
 120   if(l){
 121     Error("Output container not found");
 122     f->Close(); delete f;
 123     return 0x0;
 124   }
 125   TList *returnlist = dynamic_cast<TList *>(l->Clone());
 126   f->Close(); delete f;
 127   return returnlist;
 128 }
 129
 130 TObject* GetMeasuredSpectrum(TList *fContainer) {
 131   AliCFContainer* c = (AliCFContainer*)fContainer->FindObject("container");
 132   AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
 133   //data->SetMeasured(AliHFEcuts::kStepHFEcuts + 1);
 134   data->SetMeasured(5);
 135   return data;
 136 }
 137
 138 TObject* GetGeneratedSpectrum(TList *fContainer) {
 139   AliCFContainer* c = (AliCFContainer*)fContainer->FindObject("container");
 140   AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
 141   //data->SetMeasured(AliHFEcuts::kStepMCGenerated);
 142   data->SetMeasured(0);
 143   return data;
 144 }
 145
 146 TObject* GetEfficiency(TList *fContainer) {
 147   AliCFContainer* c = (AliCFContainer*)fContainer->FindObject("container");
 148   AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
 149   //eff->CalculateEfficiency(AliHFEcuts::kStepHFEcuts + 2,AliHFEcuts::kStepMCGenerated);
 150   eff->CalculateEfficiency(6,0);
 151   return eff;
 152 }
 153
 154 THnSparse* GetResponseMatrix(TList *fContainer) {
 155   THnSparse* h = (THnSparse*)f->Get("correlation");
 156   return h;
 157 }
 158
 159 THnSparse* CreateGuessed(const THnSparse* h) {
 160   THnSparse* guessed = (THnSparse*) h->Clone();
 161   return guessed ;
 162 }