[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / UnfoldingHF.C

//-----------------------------------------------------//
//                                                     //
//      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 ;
}
Commit	Line	Data
a5d92cb4	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