ROOT logo
//-----------------------------------------------------//
//                                                     //
//      Base Macro to evaluate pt-eta unfolding        //
//      for charm mesons                               //
//                                                     //
// Usage:                                              //
//                                                     //
// 1) By default need output.root file from AliCF...   //
//    task                                             //
// 2) ./output.root - used to evaluate the eff         //
// 3) To correct, the data sample should be spitted in //
//    2, in this case in GetMeasuredSpectrum() connect //
//    the second container ./output_data.root          //
// 4) Set the number of iterations for Baysian unf     //
// 5) Set min. Chi2                                    //
// 6) Select if unfold at Acceptance level or          //
//    after PPR cuts                                   //
// 7) Gaussian error propagation assumed, may be       //
//    errors overestimated - to be studied             //
//                                                     //
//-----------------------------------------------------//
//                                                     //
//  a.grelli@uu.nl- Utrecht for R.Vernaut example      //
//                                                     //
//-----------------------------------------------------//
//                                                     //
//   NOTE: by default the macro is computing the       //
//   unfolding after acceptance step. To evaluate      //
//   after PPR switch on the option in the CF task!    //
//   and change in this macro the selection steps      //
//                                                     //
//-----------------------------------------------------//

#include "TString.h"
#include <iostream>
#include <fstream>

void UnfoldingHF() {

  // not jet impemented the "systematics" mode!

  //TString smode(analysis_mode);
  //TString hist(method);
  
  Load();

  printf("==================================================================\n");
  printf("===========        RUNNING D MESONS UNFOLDING           ==========\n");
  printf("==================================================================\n");

  // --------------------- get variables from container ----------------------//

  AliCFDataGrid  *measured         = (AliCFDataGrid*) GetMeasuredSpectrum();
  AliCFDataGrid  *reconstructed    = (AliCFDataGrid*) GetReconstructedSpectrum();
  AliCFDataGrid  *generated        = (AliCFDataGrid*) GetGeneratedSpectrum();
  AliCFEffGrid   *efficiency       = (AliCFEffGrid*)  GetEfficiency();
  
  // the response matrix

  THnSparse      *response  = (THnSparse*) GetResponseMatrix();
    
  // -------------------------------------------------------------------------//
  //                                                                          //
  //    The HF correction framework has 12 dimensions. This is bad for        //
  //    pt-eta unfolding so reduce the dimesions from 12 to 2.                //
  //                                                                          //
  //--------------------------------------------------------------------------//
  

  THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
  THnSparse* data    = CreateGuessed(((AliCFGridSparse*)measured->GetData())->GetGrid());
  THnSparse* Reco    = CreateGuessed(((AliCFGridSparse*)reconstructed->GetData())->GetGrid()) ;
  
 
  // best guess, traslate the 12 dimentions container over 2 (eta and pt)
  
  Int_t dim[2];

  dim[0] = 0;
  dim[1] = 1;
 
  THnSparse* BestGuess      = guessed ->Projection(2,dim);
  THnSparse* DataSample     = data    ->Projection(2,dim);
  THnSparse* Reconstructed  = Reco    ->Projection(2,dim);
  THnSparse* SimpleEff      = Reco    ->Projection(2,dim);

  SimpleEff->Divide(Reconstructed,BestGuess,1.,1.);
 

  // here I do the unfolding and I set the number of iterations and the chi2

  AliCFUnfolding unfolding("unfolding","",2,response,SimpleEff,DataSample,BestGuess);
  unfolding.SetMaxNumberOfIterations(100);
  unfolding.SetMaxChi2PerDOF(0.00000000005);
  unfolding.Unfold();

  THnSparse* result = unfolding.GetUnfolded();
  

  TH2D* h_gen     =  generated   ->Project(1,0);
  TH2D* h_meas    =  DataSample  ->Projection(1,0);
  TH2D* h_reco    =  Reconstructed ->Projection(1,0);
  TH2D* h_guessed =  BestGuess     ->Projection(1,0);
  TH2D* h_eff     =  SimpleEff     ->Projection(1,0);
  TH2D* h_unf     =  result        ->Projection(1,0);

  // Apply simple efficiency correction 


  TH2D* simpleC = (TH2D*)h_meas->Clone();
  simpleC->Divide(h_eff);
  
  TCanvas * canvas3 = new TCanvas("Unfolded efficiency map","canvas 3",1000,700);

  canvas3->cd();

  TH2D* ratio = (TH2D*)h_unf->Clone();
  ratio->SetTitle("corrected using unfolding");
  ratio->Divide(h_unf,h_guessed,1,1);
 
  ratio->Draw("cont4z");

  TCanvas * canvas4 = new TCanvas("Simple Efficiency map","",1000,700);

  canvas4->cd();
  h_meas->SetTitle("simple correction");
  h_meas->Divide(h_eff);
  h_meas->Divide(h_guessed);
  h_meas->Draw("cont4z");  

  TCanvas * distribution2 = new TCanvas("dist2","",1000,700);

  distribution2->cd();

  TH1D* gen2 = generated->Project(1); // generated pt
  gen2->SetLineColor(1);
  gen2->SetTitle("D0 #eta distribution");
  gen2->GetXaxis()->SetTitle("eta");
  gen2->SetMarkerStyle(21);
  gen2->SetMarkerSize(1.);
  gen2->SetMarkerColor(1);
  gen2->Draw();

  TH1D* meas2 = measured->Project(1); // generated pt
  meas2->SetTitle("measurements");
  meas2->SetLineColor(2);
  meas2->SetMarkerStyle(22);
  meas2->SetMarkerSize(1.);
  meas2->SetMarkerColor(2);
  meas2->DrawCopy("same");

  TH1D* unf2 = result->Projection(1); // generated pt
  unf2->SetTitle("measurements");
  unf2->SetLineColor(3);
  unf2->SetMarkerStyle(23);
  unf2->SetMarkerSize(1.);
  unf2->SetMarkerColor(3);
  unf2->DrawCopy("SAME");

  corr2 = new TH1D();
  corr2->Sumw2();
  corr2 = simpleC->ProjectionY();
  corr2->SetTitle("measurements");
  corr2->SetLineColor(4);
  corr2->SetMarkerStyle(24);
  corr2->SetMarkerSize(1.);
  corr2->SetMarkerColor(4);
  corr2->Sumw2();
  corr2->DrawCopy("SAME");

  leg = new TLegend(0.1,0.7,0.48,0.9);
  leg->SetHeader("Distributions");
  leg->AddEntry(gen2,"Generated PYTHIA","l");
  leg->AddEntry(meas2,"Reconstructed","l");
  leg->AddEntry(corr2,"Corrected for eff","l");
  leg->AddEntry(unf2,"Unfolded","l");

  leg->Draw();

  TCanvas * distribution3 = new TCanvas("dist3","",1000,700);

  distribution3->cd();


  TH1D* gen3 = generated->Project(0); // generated pt
  Int_t Entries = gen3->GetEntries();
  gen3->SetLineColor(1);
  gen3->SetTitle("D0 pt distribution");
  gen3->GetXaxis()->SetTitle("pt (GeV/c)");
  gen3->SetMarkerStyle(21);
  gen3->SetMarkerSize(1.);
  gen3->SetMarkerColor(1);
  gen3->Draw();


  distribution3->cd();
  TH1D* meas3 = measured->Project(0); // generated pt
  meas3->SetTitle("measurements");
  meas3->SetLineColor(2);
  meas3->SetMarkerStyle(22);
  meas3->SetMarkerSize(1.);
  meas3->SetMarkerColor(2);
  meas3->DrawCopy("SAME");

  TH1D* unf3 = result->Projection(0); // 
  unf3->SetTitle("unfolded");
  unf3->SetLineColor(3);
  unf3->SetMarkerStyle(23);
  unf3->SetMarkerSize(1.);
  unf3->SetMarkerColor(3);
  unf3->Sumw2();
  unf3->DrawCopy("SAME");

  TH1D* corr3 =  simpleC->ProjectionX();
  corr3->SetTitle("corrected");
  corr3->SetLineColor(4);
  corr3->SetMarkerStyle(24);
  corr3->SetMarkerSize(1.);
  corr3->SetMarkerColor(4);
  corr3->Sumw2();
  corr3->DrawCopy("SAME");

  leg = new TLegend(0.1,0.7,0.48,0.9);
  leg -> SetHeader("Distributions");
  leg -> AddEntry(gen3,"Generated PYTHIA","l");
  leg -> AddEntry(meas3,"Reconstructed","l");
  leg -> AddEntry(corr3,"Corrected for eff","l");
  leg -> AddEntry(unf3,"Unfolded","l");

  leg->Draw();

  return;
}

// ====================================================================


void Load(){
  gSystem->Load("libANALYSIS");
  gSystem->Load("libCORRFW");
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
}

//___________________________________________________________-

TObject* GetMeasuredSpectrum() {
  TFile * f = TFile::Open("output.root","read");
  AliCFContainer* c = (AliCFContainer*)f->Get("container");
  AliCFDataGrid* data1 = new AliCFDataGrid("data1","",*c);
  data1->SetMeasured(3);
  return data1;
}
TObject* GetReconstructedSpectrum() {
  TFile * f = TFile::Open("output.root","read");
  AliCFContainer* c = (AliCFContainer*)f->Get("container");
  AliCFDataGrid* data2 = new AliCFDataGrid("data2","",*c);
  data2->SetMeasured(3);
 
  return data2;
}
TObject* GetGeneratedSpectrum() {
  TFile * f = TFile::Open("output.root","read");
  AliCFContainer* c = (AliCFContainer*)f->Get("container");
  AliCFDataGrid* data3 = new AliCFDataGrid("data3","",*c);
  data3->SetMeasured(1);
  return data3;
}
TObject* GetEfficiency() {
  TFile * f = TFile::Open("output.root","read");
  AliCFContainer* c = (AliCFContainer*)f->Get("container");
  AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
  eff->CalculateEfficiency(3,1);
  return eff;
}
THnSparse* GetResponseMatrix() {
  TFile * f = TFile::Open("output.root","read");
  THnSparse* h = (THnSparse*)f->Get("correlation");
  return h;
}

THnSparse* GetResponseMatrixPPR() {
  TFile * f = TFile::Open("output.root","read");
  THnSparse* h = (THnSparse*)f->Get("correlationPPR");
  return h;
}

THnSparse* CreateGuessed(const THnSparse* h) {
  THnSparse* guessed = (THnSparse*) h->Clone();
  return guessed ;
}


 UnfoldingHF.C:1
 UnfoldingHF.C:2
 UnfoldingHF.C:3
 UnfoldingHF.C:4
 UnfoldingHF.C:5
 UnfoldingHF.C:6
 UnfoldingHF.C:7
 UnfoldingHF.C:8
 UnfoldingHF.C:9
 UnfoldingHF.C:10
 UnfoldingHF.C:11
 UnfoldingHF.C:12
 UnfoldingHF.C:13
 UnfoldingHF.C:14
 UnfoldingHF.C:15
 UnfoldingHF.C:16
 UnfoldingHF.C:17
 UnfoldingHF.C:18
 UnfoldingHF.C:19
 UnfoldingHF.C:20
 UnfoldingHF.C:21
 UnfoldingHF.C:22
 UnfoldingHF.C:23
 UnfoldingHF.C:24
 UnfoldingHF.C:25
 UnfoldingHF.C:26
 UnfoldingHF.C:27
 UnfoldingHF.C:28
 UnfoldingHF.C:29
 UnfoldingHF.C:30
 UnfoldingHF.C:31
 UnfoldingHF.C:32
 UnfoldingHF.C:33
 UnfoldingHF.C:34
 UnfoldingHF.C:35
 UnfoldingHF.C:36
 UnfoldingHF.C:37
 UnfoldingHF.C:38
 UnfoldingHF.C:39
 UnfoldingHF.C:40
 UnfoldingHF.C:41
 UnfoldingHF.C:42
 UnfoldingHF.C:43
 UnfoldingHF.C:44
 UnfoldingHF.C:45
 UnfoldingHF.C:46
 UnfoldingHF.C:47
 UnfoldingHF.C:48
 UnfoldingHF.C:49
 UnfoldingHF.C:50
 UnfoldingHF.C:51
 UnfoldingHF.C:52
 UnfoldingHF.C:53
 UnfoldingHF.C:54
 UnfoldingHF.C:55
 UnfoldingHF.C:56
 UnfoldingHF.C:57
 UnfoldingHF.C:58
 UnfoldingHF.C:59
 UnfoldingHF.C:60
 UnfoldingHF.C:61
 UnfoldingHF.C:62
 UnfoldingHF.C:63
 UnfoldingHF.C:64
 UnfoldingHF.C:65
 UnfoldingHF.C:66
 UnfoldingHF.C:67
 UnfoldingHF.C:68
 UnfoldingHF.C:69
 UnfoldingHF.C:70
 UnfoldingHF.C:71
 UnfoldingHF.C:72
 UnfoldingHF.C:73
 UnfoldingHF.C:74
 UnfoldingHF.C:75
 UnfoldingHF.C:76
 UnfoldingHF.C:77
 UnfoldingHF.C:78
 UnfoldingHF.C:79
 UnfoldingHF.C:80
 UnfoldingHF.C:81
 UnfoldingHF.C:82
 UnfoldingHF.C:83
 UnfoldingHF.C:84
 UnfoldingHF.C:85
 UnfoldingHF.C:86
 UnfoldingHF.C:87
 UnfoldingHF.C:88
 UnfoldingHF.C:89
 UnfoldingHF.C:90
 UnfoldingHF.C:91
 UnfoldingHF.C:92
 UnfoldingHF.C:93
 UnfoldingHF.C:94
 UnfoldingHF.C:95
 UnfoldingHF.C:96
 UnfoldingHF.C:97
 UnfoldingHF.C:98
 UnfoldingHF.C:99
 UnfoldingHF.C:100
 UnfoldingHF.C:101
 UnfoldingHF.C:102
 UnfoldingHF.C:103
 UnfoldingHF.C:104
 UnfoldingHF.C:105
 UnfoldingHF.C:106
 UnfoldingHF.C:107
 UnfoldingHF.C:108
 UnfoldingHF.C:109
 UnfoldingHF.C:110
 UnfoldingHF.C:111
 UnfoldingHF.C:112
 UnfoldingHF.C:113
 UnfoldingHF.C:114
 UnfoldingHF.C:115
 UnfoldingHF.C:116
 UnfoldingHF.C:117
 UnfoldingHF.C:118
 UnfoldingHF.C:119
 UnfoldingHF.C:120
 UnfoldingHF.C:121
 UnfoldingHF.C:122
 UnfoldingHF.C:123
 UnfoldingHF.C:124
 UnfoldingHF.C:125
 UnfoldingHF.C:126
 UnfoldingHF.C:127
 UnfoldingHF.C:128
 UnfoldingHF.C:129
 UnfoldingHF.C:130
 UnfoldingHF.C:131
 UnfoldingHF.C:132
 UnfoldingHF.C:133
 UnfoldingHF.C:134
 UnfoldingHF.C:135
 UnfoldingHF.C:136
 UnfoldingHF.C:137
 UnfoldingHF.C:138
 UnfoldingHF.C:139
 UnfoldingHF.C:140
 UnfoldingHF.C:141
 UnfoldingHF.C:142
 UnfoldingHF.C:143
 UnfoldingHF.C:144
 UnfoldingHF.C:145
 UnfoldingHF.C:146
 UnfoldingHF.C:147
 UnfoldingHF.C:148
 UnfoldingHF.C:149
 UnfoldingHF.C:150
 UnfoldingHF.C:151
 UnfoldingHF.C:152
 UnfoldingHF.C:153
 UnfoldingHF.C:154
 UnfoldingHF.C:155
 UnfoldingHF.C:156
 UnfoldingHF.C:157
 UnfoldingHF.C:158
 UnfoldingHF.C:159
 UnfoldingHF.C:160
 UnfoldingHF.C:161
 UnfoldingHF.C:162
 UnfoldingHF.C:163
 UnfoldingHF.C:164
 UnfoldingHF.C:165
 UnfoldingHF.C:166
 UnfoldingHF.C:167
 UnfoldingHF.C:168
 UnfoldingHF.C:169
 UnfoldingHF.C:170
 UnfoldingHF.C:171
 UnfoldingHF.C:172
 UnfoldingHF.C:173
 UnfoldingHF.C:174
 UnfoldingHF.C:175
 UnfoldingHF.C:176
 UnfoldingHF.C:177
 UnfoldingHF.C:178
 UnfoldingHF.C:179
 UnfoldingHF.C:180
 UnfoldingHF.C:181
 UnfoldingHF.C:182
 UnfoldingHF.C:183
 UnfoldingHF.C:184
 UnfoldingHF.C:185
 UnfoldingHF.C:186
 UnfoldingHF.C:187
 UnfoldingHF.C:188
 UnfoldingHF.C:189
 UnfoldingHF.C:190
 UnfoldingHF.C:191
 UnfoldingHF.C:192
 UnfoldingHF.C:193
 UnfoldingHF.C:194
 UnfoldingHF.C:195
 UnfoldingHF.C:196
 UnfoldingHF.C:197
 UnfoldingHF.C:198
 UnfoldingHF.C:199
 UnfoldingHF.C:200
 UnfoldingHF.C:201
 UnfoldingHF.C:202
 UnfoldingHF.C:203
 UnfoldingHF.C:204
 UnfoldingHF.C:205
 UnfoldingHF.C:206
 UnfoldingHF.C:207
 UnfoldingHF.C:208
 UnfoldingHF.C:209
 UnfoldingHF.C:210
 UnfoldingHF.C:211
 UnfoldingHF.C:212
 UnfoldingHF.C:213
 UnfoldingHF.C:214
 UnfoldingHF.C:215
 UnfoldingHF.C:216
 UnfoldingHF.C:217
 UnfoldingHF.C:218
 UnfoldingHF.C:219
 UnfoldingHF.C:220
 UnfoldingHF.C:221
 UnfoldingHF.C:222
 UnfoldingHF.C:223
 UnfoldingHF.C:224
 UnfoldingHF.C:225
 UnfoldingHF.C:226
 UnfoldingHF.C:227
 UnfoldingHF.C:228
 UnfoldingHF.C:229
 UnfoldingHF.C:230
 UnfoldingHF.C:231
 UnfoldingHF.C:232
 UnfoldingHF.C:233
 UnfoldingHF.C:234
 UnfoldingHF.C:235
 UnfoldingHF.C:236
 UnfoldingHF.C:237
 UnfoldingHF.C:238
 UnfoldingHF.C:239
 UnfoldingHF.C:240
 UnfoldingHF.C:241
 UnfoldingHF.C:242
 UnfoldingHF.C:243
 UnfoldingHF.C:244
 UnfoldingHF.C:245
 UnfoldingHF.C:246
 UnfoldingHF.C:247
 UnfoldingHF.C:248
 UnfoldingHF.C:249
 UnfoldingHF.C:250
 UnfoldingHF.C:251
 UnfoldingHF.C:252
 UnfoldingHF.C:253
 UnfoldingHF.C:254
 UnfoldingHF.C:255
 UnfoldingHF.C:256
 UnfoldingHF.C:257
 UnfoldingHF.C:258
 UnfoldingHF.C:259
 UnfoldingHF.C:260
 UnfoldingHF.C:261
 UnfoldingHF.C:262
 UnfoldingHF.C:263
 UnfoldingHF.C:264
 UnfoldingHF.C:265
 UnfoldingHF.C:266
 UnfoldingHF.C:267
 UnfoldingHF.C:268
 UnfoldingHF.C:269
 UnfoldingHF.C:270
 UnfoldingHF.C:271
 UnfoldingHF.C:272
 UnfoldingHF.C:273
 UnfoldingHF.C:274
 UnfoldingHF.C:275
 UnfoldingHF.C:276
 UnfoldingHF.C:277
 UnfoldingHF.C:278
 UnfoldingHF.C:279
 UnfoldingHF.C:280
 UnfoldingHF.C:281
 UnfoldingHF.C:282
 UnfoldingHF.C:283
 UnfoldingHF.C:284
 UnfoldingHF.C:285
 UnfoldingHF.C:286
 UnfoldingHF.C:287
 UnfoldingHF.C:288
 UnfoldingHF.C:289
 UnfoldingHF.C:290
 UnfoldingHF.C:291
 UnfoldingHF.C:292
 UnfoldingHF.C:293
 UnfoldingHF.C:294