ROOT logo

void testUnfolding() {
  TBenchmark b;

  b.Start("init");

  Load();

  AliLog::SetGlobalDebugLevel(0);

  // get the essential
  AliCFDataGrid *measured    = (AliCFDataGrid*) GetMeasuredSpectrum();
  AliCFDataGrid *generated   = (AliCFDataGrid*) GetGeneratedSpectrum();
  AliCFEffGrid  *efficiency  = (AliCFEffGrid*)  GetEfficiency();
  THnSparse     *response    = (THnSparse*)     GetResponseMatrix();
  
  // create a guessed "a priori" distribution using binning of MC
  THnSparse* guessed = CreateGuessed(((AliCFGridSparse*)generated->GetData())->GetGrid()) ;
  //----

  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);
  fit2D->SetParameter(0,1);
  fit2D->SetParameter(1,1);
  fit2D->SetParameter(2,1);
  fit2D->SetParameter(3,1);
  fit2D->SetParameter(4,1);
  fit2D->SetParLimits(1,0,1);

  AliCFUnfolding unfolding("unfolding","",2,response,efficiency->GetGrid(),((AliCFGridSparse*)measured->GetData())->GetGrid()/*,guessed*/);
  unfolding.SetMaxNumberOfIterations(100);
  unfolding.SetMaxChi2PerDOF(1.e-07);
  unfolding.UseSmoothing(fit2D,"ren");
  //unfolding.UseSmoothing();
  
  b.Stop("init");
  b.Start("unfolding");
  unfolding.Unfold();
  b.Stop("unfolding");
  b.Start("finish");

  canvas->cd();

  TH2D* h_gen = generated->Project(0,1);
  h_gen->SetTitle("generated");
  h_gen->Draw("lego2");
  canvas->SaveAs("/tmp/generated.eps");
  canvas->SaveAs("/tmp/generated.gif");

  TH2D* h_meas = measured->Project(0,1);
  h_meas->SetTitle("measured");
  h_meas->Draw("lego2");
  canvas->SaveAs("/tmp/measured.eps");
  canvas->SaveAs("/tmp/measured.gif");

  TH2D* h_guessed = unfolding.GetOriginalPrior()->Projection(1,0);
  h_guessed->SetTitle("a priori");
  h_guessed->Draw("lego2");
  canvas->SaveAs("/tmp/apriori.eps");
  canvas->SaveAs("/tmp/apriori.gif");

  TH2D* h_eff = efficiency->Project(0,1);
  h_eff->SetTitle("efficiency");
  h_eff->Draw("e");

  TH2D* h_unf = unfolding.GetUnfolded()->Projection(1,0);
  h_unf->SetTitle("unfolded");
  h_unf->Draw("lego2");
  fit2D->Draw("surf same");
  return;
  canvas->SaveAs("/tmp/unfolded.eps");
  canvas->SaveAs("/tmp/unfolded.gif");


  TH2D* ratio = (TH2D*)h_unf->Clone();
  ratio->SetTitle("unfolded / generated");
  ratio->Divide(h_unf,h_gen,1,1);
  ratio->GetZaxis()->SetRangeUser(0.5,1.5);
  //ratio->Draw("cont4z");
  ratio->Draw("lego2");
  //ratio->Draw("e");
//   ratio->ProjectionY()->Draw();
  canvas->SaveAs("/tmp/ratio.eps");
  canvas->SaveAs("/tmp/ratio.gif");
   

  canvas->cd(7);
  TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
  orig->SetTitle("original prior");
  orig->Draw("lego2");

  canvas->cd(8);
  TH2D* h_est = (TH2D*) unfolding.GetEstMeasured()->Projection(1,0);
  h_est->SetTitle("est. measured");
  h_est->Draw("e");

  canvas->cd(9);
  unfolding.GetUnfolded()->Projection(0)->Draw();

  canvas->cd(10);
  unfolding.GetUnfolded()->Projection(1)->Draw();

  canvas->cd(11);
  h_gen->ProjectionX()->Draw();

  canvas->cd(12);
  h_gen->ProjectionY()->Draw();

  b.Stop("finish");
  Float_t x,y;
  b.Summary(x,y);
}

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


void Load(){
  gSystem->Load("libANALYSIS");
  gSystem->Load("libCORRFW");
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  TCanvas * canvas = new TCanvas("canvas","",600,400);
  canvas->SetFillColor(10);
  canvas->SetFrameFillColor(10);
}

TObject* GetMeasuredSpectrum() {
  TFile * f = TFile::Open("test/output.root","read");
  AliCFContainer* c = (AliCFContainer*)f->Get("container");
  AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
  data->SetMeasured(1);
  return data;
}
TObject* GetGeneratedSpectrum() {
  TFile * f = TFile::Open("test/output.root","read");
  AliCFContainer* c = (AliCFContainer*)f->Get("container");
  AliCFDataGrid* data = new AliCFDataGrid("data","",*c);
  data->SetMeasured(0);
  return data;
}
TObject* GetEfficiency() {
  TFile * f = TFile::Open("test/output.root","read");
  AliCFContainer* c = (AliCFContainer*)f->Get("container");
  AliCFEffGrid* eff = new AliCFEffGrid("eff","",*c);
  eff->CalculateEfficiency(2,0);
  return eff;
}
THnSparse* GetResponseMatrix() {
  TFile * f = TFile::Open("test/output.root","read");
  THnSparse* h = (THnSparse*)f->Get("correlation");
  return h;
}
THnSparse* CreateGuessed(const THnSparse* h) {
  THnSparse* guessed = (THnSparse*) h->Clone();
  return guessed ;
}
 testUnfolding.C:1
 testUnfolding.C:2
 testUnfolding.C:3
 testUnfolding.C:4
 testUnfolding.C:5
 testUnfolding.C:6
 testUnfolding.C:7
 testUnfolding.C:8
 testUnfolding.C:9
 testUnfolding.C:10
 testUnfolding.C:11
 testUnfolding.C:12
 testUnfolding.C:13
 testUnfolding.C:14
 testUnfolding.C:15
 testUnfolding.C:16
 testUnfolding.C:17
 testUnfolding.C:18
 testUnfolding.C:19
 testUnfolding.C:20
 testUnfolding.C:21
 testUnfolding.C:22
 testUnfolding.C:23
 testUnfolding.C:24
 testUnfolding.C:25
 testUnfolding.C:26
 testUnfolding.C:27
 testUnfolding.C:28
 testUnfolding.C:29
 testUnfolding.C:30
 testUnfolding.C:31
 testUnfolding.C:32
 testUnfolding.C:33
 testUnfolding.C:34
 testUnfolding.C:35
 testUnfolding.C:36
 testUnfolding.C:37
 testUnfolding.C:38
 testUnfolding.C:39
 testUnfolding.C:40
 testUnfolding.C:41
 testUnfolding.C:42
 testUnfolding.C:43
 testUnfolding.C:44
 testUnfolding.C:45
 testUnfolding.C:46
 testUnfolding.C:47
 testUnfolding.C:48
 testUnfolding.C:49
 testUnfolding.C:50
 testUnfolding.C:51
 testUnfolding.C:52
 testUnfolding.C:53
 testUnfolding.C:54
 testUnfolding.C:55
 testUnfolding.C:56
 testUnfolding.C:57
 testUnfolding.C:58
 testUnfolding.C:59
 testUnfolding.C:60
 testUnfolding.C:61
 testUnfolding.C:62
 testUnfolding.C:63
 testUnfolding.C:64
 testUnfolding.C:65
 testUnfolding.C:66
 testUnfolding.C:67
 testUnfolding.C:68
 testUnfolding.C:69
 testUnfolding.C:70
 testUnfolding.C:71
 testUnfolding.C:72
 testUnfolding.C:73
 testUnfolding.C:74
 testUnfolding.C:75
 testUnfolding.C:76
 testUnfolding.C:77
 testUnfolding.C:78
 testUnfolding.C:79
 testUnfolding.C:80
 testUnfolding.C:81
 testUnfolding.C:82
 testUnfolding.C:83
 testUnfolding.C:84
 testUnfolding.C:85
 testUnfolding.C:86
 testUnfolding.C:87
 testUnfolding.C:88
 testUnfolding.C:89
 testUnfolding.C:90
 testUnfolding.C:91
 testUnfolding.C:92
 testUnfolding.C:93
 testUnfolding.C:94
 testUnfolding.C:95
 testUnfolding.C:96
 testUnfolding.C:97
 testUnfolding.C:98
 testUnfolding.C:99
 testUnfolding.C:100
 testUnfolding.C:101
 testUnfolding.C:102
 testUnfolding.C:103
 testUnfolding.C:104
 testUnfolding.C:105
 testUnfolding.C:106
 testUnfolding.C:107
 testUnfolding.C:108
 testUnfolding.C:109
 testUnfolding.C:110
 testUnfolding.C:111
 testUnfolding.C:112
 testUnfolding.C:113
 testUnfolding.C:114
 testUnfolding.C:115
 testUnfolding.C:116
 testUnfolding.C:117
 testUnfolding.C:118
 testUnfolding.C:119
 testUnfolding.C:120
 testUnfolding.C:121
 testUnfolding.C:122
 testUnfolding.C:123
 testUnfolding.C:124
 testUnfolding.C:125
 testUnfolding.C:126
 testUnfolding.C:127
 testUnfolding.C:128
 testUnfolding.C:129
 testUnfolding.C:130
 testUnfolding.C:131
 testUnfolding.C:132
 testUnfolding.C:133
 testUnfolding.C:134
 testUnfolding.C:135
 testUnfolding.C:136
 testUnfolding.C:137
 testUnfolding.C:138
 testUnfolding.C:139
 testUnfolding.C:140
 testUnfolding.C:141
 testUnfolding.C:142
 testUnfolding.C:143
 testUnfolding.C:144
 testUnfolding.C:145
 testUnfolding.C:146
 testUnfolding.C:147
 testUnfolding.C:148
 testUnfolding.C:149
 testUnfolding.C:150
 testUnfolding.C:151
 testUnfolding.C:152
 testUnfolding.C:153
 testUnfolding.C:154
 testUnfolding.C:155
 testUnfolding.C:156