PWGHF/vertexingHF/macros/UnfoldingHF.C

   1 //-----------------------------------------------------//
   2 //                                                     //
   3 //      Base Macro to evaluate pt-eta unfolding        //
   4 //      for charm mesons                               //
   5 //                                                     //
   6 // Usage:                                              //
   7 //                                                     //
   8 // 1) By default need output.root file from AliCF...   //
   9 //    task                                             //
  10 // 2) ./output.root - used to evaluate the eff         //
  11 // 3) To correct, the data sample should be spitted in //
  12 //    2, in this case in GetMeasuredSpectrum() connect //
  13 //    the second container ./output_data.root          //
  14 // 4) Set the number of iterations for Baysian unf     //
  15 // 5) Set min. Chi2                                    //
  16 // 6) Select if unfold at Acceptance level or          //
  17 //    after PPR cuts                                   //
  18 // 7) Gaussian error propagation assumed, may be       //
  19 //    errors overestimated - to be studied             //
  20 //                                                     //
  21 //-----------------------------------------------------//
  22 //                                                     //
  23 //  a.grelli@uu.nl- Utrecht for R.Vernaut example      //
  24 //                                                     //
  25 //-----------------------------------------------------//
  26 //                                                     //
  27 //   NOTE: by default the macro is computing the       //
  28 //   unfolding after acceptance step. To evaluate      //
  29 //   after PPR switch on the option in the CF task!    //
  30 //   and change in this macro the selection steps      //
  31 //                                                     //
  32 //-----------------------------------------------------//
  33
  34 #include "TString.h"
  35 #include <iostream>
  36 #include <fstream>
  37
  38 void UnfoldingHF() {
  39
  40   // not jet impemented the "systematics" mode!
  41
  42   //TString smode(analysis_mode);
  43   //TString hist(method);
  44
  45   Load();
  46
  47   printf("==================================================================\n");
  48   printf("===========        RUNNING D MESONS UNFOLDING           ==========\n");
  49   printf("==================================================================\n");
  50
  51   // --------------------- get variables from container ----------------------//
  52
  53   AliCFDataGrid  *measured         = (AliCFDataGrid*) GetMeasuredSpectrum();
  54   AliCFDataGrid  *reconstructed    = (AliCFDataGrid*) GetReconstructedSpectrum();
  55   AliCFDataGrid  *generated        = (AliCFDataGrid*) GetGeneratedSpectrum();
  56   AliCFEffGrid   *efficiency       = (AliCFEffGrid*)  GetEfficiency();
  57
  58   // the response matrix
  59
  60   THnSparse      *response  = (THnSparse*) GetResponseMatrix();
  61
  62   // -------------------------------------------------------------------------//
  63   //                                                                          //
  64   //    The HF correction framework has 12 dimensions. This is bad for        //
  65   //    pt-eta unfolding so reduce the dimesions from 12 to 2.                //
  66   //                                                                          //
  67   //--------------------------------------------------------------------------//
  68
  69
  70   THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
  71   THnSparse* data    = CreateGuessed(((AliCFGridSparse*)measured->GetData())->GetGrid());
  72   THnSparse* Reco    = CreateGuessed(((AliCFGridSparse*)reconstructed->GetData())->GetGrid()) ;
  73
  74
  75   // best guess, traslate the 12 dimentions container over 2 (eta and pt)
  76
  77   Int_t dim[2];
  78
  79   dim[0] = 0;
  80   dim[1] = 1;
  81
  82   THnSparse* BestGuess      = guessed ->Projection(2,dim);
  83   THnSparse* DataSample     = data    ->Projection(2,dim);
  84   THnSparse* Reconstructed  = Reco    ->Projection(2,dim);
  85   THnSparse* SimpleEff      = Reco    ->Projection(2,dim);
  86
  87   SimpleEff->Divide(Reconstructed,BestGuess,1.,1.);
  88
  89
  90   // here I do the unfolding and I set the number of iterations and the chi2
  91
  92   AliCFUnfolding unfolding("unfolding","",2,response,SimpleEff,DataSample,BestGuess);
  93   unfolding.SetMaxNumberOfIterations(100);
  94   unfolding.SetMaxChi2PerDOF(0.00000000005);
  95   unfolding.Unfold();
  96
  97   THnSparse* result = unfolding.GetUnfolded();
  98
  99
 100   TH2D* h_gen     =  generated   ->Project(1,0);
 101   TH2D* h_meas    =  DataSample  ->Projection(1,0);
 102   TH2D* h_reco    =  Reconstructed ->Projection(1,0);
 103   TH2D* h_guessed =  BestGuess     ->Projection(1,0);
 104   TH2D* h_eff     =  SimpleEff     ->Projection(1,0);
 105   TH2D* h_unf     =  result        ->Projection(1,0);
 106
 107   // Apply simple efficiency correction
 108
 109
 110   TH2D* simpleC = (TH2D*)h_meas->Clone();
 111   simpleC->Divide(h_eff);
 112
 113   TCanvas * canvas3 = new TCanvas("Unfolded efficiency map","canvas 3",1000,700);
 114
 115   canvas3->cd();
 116
 117   TH2D* ratio = (TH2D*)h_unf->Clone();
 118   ratio->SetTitle("corrected using unfolding");
 119   ratio->Divide(h_unf,h_guessed,1,1);
 120
 121   ratio->Draw("cont4z");
 122
 123   TCanvas * canvas4 = new TCanvas("Simple Efficiency map","",1000,700);
 124
 125   canvas4->cd();
 126   h_meas->SetTitle("simple correction");
 127   h_meas->Divide(h_eff);
 128   h_meas->Divide(h_guessed);
 129   h_meas->Draw("cont4z");
 130
 131   TCanvas * distribution2 = new TCanvas("dist2","",1000,700);
 132
 133   distribution2->cd();
 134
 135   TH1D* gen2 = generated->Project(1); // generated pt
 136   gen2->SetLineColor(1);
 137   gen2->SetTitle("D0 #eta distribution");
 138   gen2->GetXaxis()->SetTitle("eta");
 139   gen2->SetMarkerStyle(21);
 140   gen2->SetMarkerSize(1.);
 141   gen2->SetMarkerColor(1);
 142   gen2->Draw();
 143
 144   TH1D* meas2 = measured->Project(1); // generated pt
 145   meas2->SetTitle("measurements");
 146   meas2->SetLineColor(2);
 147   meas2->SetMarkerStyle(22);
 148   meas2->SetMarkerSize(1.);
 149   meas2->SetMarkerColor(2);
 150   meas2->DrawCopy("same");
 151
 152   TH1D* unf2 = result->Projection(1); // generated pt
 153   unf2->SetTitle("measurements");
 154   unf2->SetLineColor(3);
 155   unf2->SetMarkerStyle(23);
 156   unf2->SetMarkerSize(1.);
 157   unf2->SetMarkerColor(3);
 158   unf2->DrawCopy("SAME");
 159
 160   corr2 = new TH1D();
 161   corr2->Sumw2();
 162   corr2 = simpleC->ProjectionY();
 163   corr2->SetTitle("measurements");
 164   corr2->SetLineColor(4);
 165   corr2->SetMarkerStyle(24);
 166   corr2->SetMarkerSize(1.);
 167   corr2->SetMarkerColor(4);
 168   corr2->Sumw2();
 169   corr2->DrawCopy("SAME");
 170
 171   leg = new TLegend(0.1,0.7,0.48,0.9);
 172   leg->SetHeader("Distributions");
 173   leg->AddEntry(gen2,"Generated PYTHIA","l");
 174   leg->AddEntry(meas2,"Reconstructed","l");
 175   leg->AddEntry(corr2,"Corrected for eff","l");
 176   leg->AddEntry(unf2,"Unfolded","l");
 177
 178   leg->Draw();
 179
 180   TCanvas * distribution3 = new TCanvas("dist3","",1000,700);
 181
 182   distribution3->cd();
 183
 184
 185   TH1D* gen3 = generated->Project(0); // generated pt
 186   Int_t Entries = gen3->GetEntries();
 187   gen3->SetLineColor(1);
 188   gen3->SetTitle("D0 pt distribution");
 189   gen3->GetXaxis()->SetTitle("pt (GeV/c)");
 190   gen3->SetMarkerStyle(21);
 191   gen3->SetMarkerSize(1.);
 192   gen3->SetMarkerColor(1);
 193   gen3->Draw();
 194
 195
 196   distribution3->cd();
 197   TH1D* meas3 = measured->Project(0); // generated pt
 198   meas3->SetTitle("measurements");
 199   meas3->SetLineColor(2);
 200   meas3->SetMarkerStyle(22);
 201   meas3->SetMarkerSize(1.);
 202   meas3->SetMarkerColor(2);
 203   meas3->DrawCopy("SAME");
 204
 205   TH1D* unf3 = result->Projection(0); //
 206   unf3->SetTitle("unfolded");
 207   unf3->SetLineColor(3);
 208   unf3->SetMarkerStyle(23);
 209   unf3->SetMarkerSize(1.);
 210   unf3->SetMarkerColor(3);
 211   unf3->Sumw2();
 212   unf3->DrawCopy("SAME");
 213
 214   TH1D* corr3 =  simpleC->ProjectionX();
 215   corr3->SetTitle("corrected");
 216   corr3->SetLineColor(4);
 217   corr3->SetMarkerStyle(24);
 218   corr3->SetMarkerSize(1.);
 219   corr3->SetMarkerColor(4);
 220   corr3->Sumw2();
 221   corr3->DrawCopy("SAME");
 222
 223   leg = new TLegend(0.1,0.7,0.48,0.9);
 224   leg -> SetHeader("Distributions");
 225   leg -> AddEntry(gen3,"Generated PYTHIA","l");
 226   leg -> AddEntry(meas3,"Reconstructed","l");
 227   leg -> AddEntry(corr3,"Corrected for eff","l");
 228   leg -> AddEntry(unf3,"Unfolded","l");
 229
 230   leg->Draw();
 231
 232   return;
 233 }
 234
 235 // ====================================================================
 236
 237
 238 void Load(){
 239   gSystem->Load("libANALYSIS");
 240   gSystem->Load("libCORRFW");
 241   gStyle->SetPalette(1);
 242   gStyle->SetOptStat(0);
 243 }
 244
 245 //___________________________________________________________-
 246
 247 TObject* GetMeasuredSpectrum() {
 248   TFile * f = TFile::Open("output.root","read");
 249   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 250   AliCFDataGrid* data1 = new AliCFDataGrid("data1","",*c);
 251   data1->SetMeasured(3);
 252   return data1;
 253 }
 254 TObject* GetReconstructedSpectrum() {
 255   TFile * f = TFile::Open("output.root","read");
 256   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 257   AliCFDataGrid* data2 = new AliCFDataGrid("data2","",*c);
 258   data2->SetMeasured(3);
 259
 260   return data2;
 261 }
 262 TObject* GetGeneratedSpectrum() {
 263   TFile * f = TFile::Open("output.root","read");
 264   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 265   AliCFDataGrid* data3 = new AliCFDataGrid("data3","",*c);
 266   data3->SetMeasured(1);
 267   return data3;
 268 }
 269 TObject* GetEfficiency() {
 270   TFile * f = TFile::Open("output.root","read");
 271   AliCFContainer* c = (AliCFContainer*)f->Get("container");
 272   AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
 273   eff->CalculateEfficiency(3,1);
 274   return eff;
 275 }
 276 THnSparse* GetResponseMatrix() {
 277   TFile * f = TFile::Open("output.root","read");
 278   THnSparse* h = (THnSparse*)f->Get("correlation");
 279   return h;
 280 }
 281
 282 THnSparse* GetResponseMatrixPPR() {
 283   TFile * f = TFile::Open("output.root","read");
 284   THnSparse* h = (THnSparse*)f->Get("correlationPPR");
 285   return h;
 286 }
 287
 288 THnSparse* CreateGuessed(const THnSparse* h) {
 289   THnSparse* guessed = (THnSparse*) h->Clone();
 290   return guessed ;
 291 }
 292
 293