Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Mat_map.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Mat_map.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <TROOT.h>
10 
11 #include <string>
12 #include <fstream>
13 #include <iostream>
14 #include <sstream>
15 
16 using namespace std ;
17 
19 
20 void Draw_ratio(TCanvas* c, TProfile* h1, TProfile* h2, TLegend* leg, std::string name){
21 
22  // Upper plot will be in pad1
23  TPad *pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
24 
25  pad1->SetBottomMargin(0); // Upper and lower plot are joined
26  pad1->SetGridx(); // Vertical grid
27  pad1->Draw(); // Draw the upper pad: pad1
28  pad1->cd(); // pad1 becomes the current pad
29  h1->Draw("E");
30  h2->Draw("E SAME");
31  leg->Draw("SAME");
32 
33  double max_hist [4];
34  max_hist[0] = h1->GetMaximum()+h1->GetMaximum()*0.05;
35  max_hist[1] = h2->GetMaximum()+h2->GetMaximum()*0.05;
36 
37  h1->GetYaxis()->SetTitleOffset(1.5);
38  h1->GetYaxis()->SetRangeUser(0, *max_element( begin(max_hist),end(max_hist) ) );
39 
40  // lower plot will be in pad
41  c->cd(); // Go back to the main canvas before defining pad2
42  TPad *pad2 = new TPad("pad2", "pad2", 0, 0.05, 1, 0.3);
43  pad2->SetTopMargin(0);
44  pad2->SetBottomMargin(0.2);
45  pad2->SetGridx(); // vertical grid
46  pad2->Draw();
47  pad2->cd(); // pad2 becomes the current pad
48 
49  TLine line = TLine(h1->GetXaxis()->GetXmin(),1,h1->GetXaxis()->GetXmax(),1);
50  line.SetLineColor(kRed);
51  line.SetLineWidth(1);
52 
53  // Define the ratio plot
54  TProfile *h5 = (TProfile*)h2->Clone( (name+"_ratio").c_str() );
55  h5->SetLineColor(kBlack);
56 
57  h5->SetStats(0); // No statistics on lower plot
58  h5->Divide(h1);
59 
60  double maxi = min( max( fabs(h5->GetMinimum()-0.1*h5->GetMinimum()),h5->GetMaximum()+0.1*h5->GetMaximum() ), 10. );
61 
62  h5->SetMinimum( 0.5 ); // Define Y ..
63  h5->SetMaximum( 1.1 ); // .. range
64 
65  h5->SetMarkerStyle(7);
66  h5->Draw("hist P"); // Draw the ratio plot
67  line.Draw("SAME");
68 
69  // Y axis ratio plot settings
70  h5->GetYaxis()->SetTitle("Ratio Validation/Geantino ");
71  h5->GetYaxis()->SetNdivisions(505);
72  h5->GetYaxis()->SetTitleSize(15);
73  h5->GetYaxis()->SetTitleFont(43);
74  h5->GetYaxis()->SetTitleOffset(1.55);
75  h5->GetYaxis()->SetLabelFont(43);
76  h5->GetYaxis()->SetLabelSize(15);
77 
78  // X axis ratio plot settings
79  h5->GetXaxis()->SetTitleSize(20);
80  h5->GetXaxis()->SetTitleFont(43);
81  h5->GetXaxis()->SetTitleOffset(3);
82  h5->GetXaxis()->SetLabelFont(43);
83  h5->GetXaxis()->SetLabelSize(15);
84 
85  std::string name_axis = h1->GetXaxis()->GetTitle();
86 
87  c->Print( (name+"/Ratio_Val_geant_mat_X0_"+name_axis+".pdf").c_str());
88 
89  delete h5;
90 
91 }
92 
96 
97 void Mat_map(std::string Val = "", std::string geantino = "", std::string name = ""){
98 
99  gStyle->SetOptStat(0);
100  gStyle->SetOptTitle(0);
101 
102  TProfile * Val_X0_Eta = new TProfile("Val_X0_Eta","Val_X0_Eta",160,-4,4);
103  TProfile * Val_X0_Phi = new TProfile("Val_X0_Phi","Val_X0_Phi",160,-4,4);
104  TH2F * Val_X0_Eta_spread = new TH2F("Val_X0_Eta_spread","Val_X0_Eta_spread",160,-4,4,160,0,4);
105  TH2F * Val_X0_Phi_spread = new TH2F("Val_X0_Phi_spread","Val_X0_Phi_spread",160,-4,4,160,0,4);
106  TH2F * Val_X0 = new TH2F("Val_X0","Val_X0",160,-4,4,160,-4,4);
107 
108  TProfile * geantino_X0_Eta = new TProfile("geantino_X0_Eta","geantino_X0_Eta",160,-4,4);
109  TProfile * geantino_X0_Phi = new TProfile("geantino_X0_Phi","geantino_X0_Phi",160,-4,4);
110  TH2F * geantino_X0_Eta_spread = new TH2F("geantino_X0_Eta_spread","geantino_X0_Eta_spread",160,-4,4,160,0,4);
111  TH2F * geantino_X0_Phi_spread = new TH2F("geantino_X0_Phi_spread","geantino_X0_Phi_spread",160,-4,4,160,0,4);
112  TH2F * geantino_X0 = new TH2F("geantino_X0","geantinol_X0",160,-4,4,160,-4,4);
113 
114  TH1F * comp_X0_Eta = new TH1F("comp_X0_Eta","comp_X0_Eta",160,-4,4);
115  TH1F * comp_X0_Phi = new TH1F("comp_X0_Phi","comp_X0_Phi",160,-4,4);
116  TH2F * comp_X0 = new TH2F("comp_X0","comp_X0",160,-4,4,160,-4,4);
117 
118  TChain *Val_file = new TChain("material-tracks");
119  TChain *geantino_file = new TChain("material-tracks");
120 
121  // Define line corresponding to the different eta value
122  TLine *eta_0 = new TLine(0,-1200,0,1200);
123  eta_0->SetLineColor(kRed);
124 
125  TLine *eta_1p = new TLine(-1250,-1064,1250,1064);
126  eta_1p->SetLineColor(kRed);
127  TLine *eta_2p = new TLine(-3000,-827,3000,827);
128  eta_2p->SetLineColor(kRed);
129  TLine *eta_3p = new TLine(-3000,-300,3000,300);
130  eta_3p->SetLineColor(kRed);
131  TLine *eta_4p = new TLine(-3000,-110,3000,110);
132  eta_4p->SetLineColor(kRed);
133 
134  TLine *eta_1n = new TLine(-1250,1064,1250,-1064);
135  eta_1n->SetLineColor(kRed);
136  TLine *eta_2n = new TLine(-3000,827,3000,-827);
137  eta_2n->SetLineColor(kRed);
138  TLine *eta_3n = new TLine(-3000,300,3000,-300);
139  eta_3n->SetLineColor(kRed);
140  TLine *eta_4n = new TLine(-3000,110,3000,-110);
141  eta_4n->SetLineColor(kRed);
142 
143  if(Val != ""){
144  Val_file->Add(Val.c_str());
145 
146  // 2D map for Validation input
147  TCanvas *VM = new TCanvas("VM","Validation Map") ;
148  Val_file->Draw("mat_y:mat_z","fabs(mat_x)<1");
149 
150  eta_0->Draw("Same");
151  eta_1p->Draw("Same");
152  eta_1n->Draw("Same");
153  eta_2p->Draw("Same");
154  eta_2n->Draw("Same");
155  eta_3p->Draw("Same");
156  eta_3n->Draw("Same");
157  eta_4p->Draw("Same");
158  eta_4n->Draw("Same");
159 
160  VM->Print( (name+"/Val_mat_map.png").c_str());
161  //VM->Print( (name+"/Val_mat_map.pdf").c_str());
162 
163  // X0 as function of Eta for Validation input
164  TCanvas *VM_X0_Eta = new TCanvas("VM_X0_Eta","Validation X0 Eta") ;
165  Val_file->Draw("t_X0:v_eta>>Val_X0_Eta","","profile");
166  Val_X0_Eta->SetMarkerStyle(7);
167  Val_X0_Eta->Draw("HIST PC");
168  Val_X0_Eta->GetXaxis()->SetTitle("Eta");
169  Val_X0_Eta->GetYaxis()->SetTitle("X0");
170  VM_X0_Eta->Print( (name+"/Val_mat_Eta_X0.pdf").c_str());
171 
172  // X0 as function of Eta for Validation input
173  TCanvas *VM_X0_Eta_spread = new TCanvas("VM_X0_Eta_spread","Validation X0 Eta") ;
174  Val_file->Draw("t_X0:v_eta>>Val_X0_Eta_spread","","");
175  Val_X0_Eta_spread->GetXaxis()->SetTitle("Eta");
176  Val_X0_Eta_spread->GetYaxis()->SetTitle("X0");
177  VM_X0_Eta_spread->Print( (name+"/Val_X0_Eta_spread.pdf").c_str());
178 
179  // X0 as function of Phi for Validation input
180  TCanvas *VM_X0_Phi = new TCanvas("VM_X0_Phi","Validation X0 Phi") ;
181  Val_file->Draw("t_X0:v_phi>>Val_X0_Phi","","profile");
182  Val_X0_Phi->SetMarkerStyle(7);
183  Val_X0_Phi->Draw("HIST PC");
184  Val_X0_Phi->GetXaxis()->SetTitle("Phi");
185  Val_X0_Phi->GetYaxis()->SetTitle("X0");
186  VM_X0_Phi->Print( (name+"/Val_mat_Phi_X0.pdf").c_str());
187 
188  // X0 as function of Phi for Validation input
189  TCanvas *VM_X0_Phi_spread = new TCanvas("VM_X0_Phi_spread","Validation X0 Phi") ;
190  Val_file->Draw("t_X0:v_phi>>Val_X0_Phi_spread","","");
191  Val_X0_Phi_spread->GetXaxis()->SetTitle("Phi");
192  Val_X0_Phi_spread->GetYaxis()->SetTitle("X0");
193  VM_X0_Phi_spread->Print( (name+"/Val_mat_Phi_X0_spread.pdf").c_str());
194 
195  // 2D map of X0 for Validation input
196  TCanvas *VM_2D = new TCanvas("VM_2D","Validation X0 2D") ;
197  Val_file->Draw("v_phi:v_eta:t_X0>>Val_X0","","COLZ");
198  Val_X0->GetXaxis()->SetTitle("Eta");
199  Val_X0->GetYaxis()->SetTitle("Phi");
200  Val_X0->GetZaxis()->SetTitle("X0");
201  VM_2D->Print( (name+"/Val_mat_X0.png").c_str());
202  }
203 
204  if(geantino != ""){
205  geantino_file->Add(geantino.c_str());
206 
207  // 2D map for Geantino input
208  TCanvas *GM = new TCanvas("GM","Geantino Map") ;
209  geantino_file->Draw("mat_y:mat_z","fabs(mat_x)<1");
210 
211  eta_0->Draw("Same");
212  eta_1p->Draw("Same");
213  eta_1n->Draw("Same");
214  eta_2p->Draw("Same");
215  eta_2n->Draw("Same");
216  eta_3p->Draw("Same");
217  eta_3n->Draw("Same");
218  eta_4p->Draw("Same");
219  eta_4n->Draw("Same");
220 
221  GM->Print( (name+"/geant_mat_map.png").c_str());
222  //GM->Print( (name+"/geant_mat_map.pdf").c_str());
223 
224  // X0 as function of Eta for Geantino input
225  TCanvas *GM_X0_Eta = new TCanvas("GM_X0_Eta","Geantino X0 Eta") ;
226  geantino_file->Draw("t_X0:v_eta>>geantino_X0_Eta","","profile");
227  geantino_X0_Eta->SetMarkerStyle(7);
228  geantino_X0_Eta->Draw("HIST PC");
229  geantino_X0_Eta->GetXaxis()->SetTitle("Eta");
230  geantino_X0_Eta->GetYaxis()->SetTitle("X0");
231  GM_X0_Eta->Print( (name+"/geant_mat_Eta_X0.pdf").c_str());
232 
233  // X0 as function of Eta for Geantino input
234  TCanvas *GM_X0_Eta_spread = new TCanvas("GM_X0_Eta_spread","Geantino X0 Eta") ;
235  geantino_file->Draw("t_X0:v_eta>>geantino_X0_Eta_spread","","");
236  geantino_X0_Eta_spread->GetXaxis()->SetTitle("Eta");
237  geantino_X0_Eta_spread->GetYaxis()->SetTitle("X0");
238  GM_X0_Eta_spread->Print( (name+"/geant_mat_Eta_X0_spread.pdf").c_str());
239 
240  // X0 as function of Phi for Geantino input
241  TCanvas *GM_X0_Phi = new TCanvas("GM_X0_Phi","Geantino X0 Phi") ;
242  geantino_file->Draw("t_X0:v_phi>>geantino_X0_Phi","","profile");
243  geantino_X0_Phi->SetMarkerStyle(7);
244  geantino_X0_Phi->Draw("HIST PC");
245  geantino_X0_Phi->GetXaxis()->SetTitle("Phi");
246  geantino_X0_Phi->GetYaxis()->SetTitle("X0");
247  GM_X0_Phi->Print( (name+"/geant_mat_Phi_X0.pdf").c_str());
248 
249  // X0 as function of Phi for Geantino input
250  TCanvas *GM_X0_Phi_spread = new TCanvas("GM_X0_Phi_spread","Geantino X0 Phi") ;
251  geantino_file->Draw("t_X0:v_phi>>geantino_X0_Phi_spread","","");
252  geantino_X0_Phi_spread->GetXaxis()->SetTitle("Phi");
253  geantino_X0_Phi_spread->GetYaxis()->SetTitle("X0");
254  GM_X0_Phi_spread->Print( (name+"/geant_mat_Phi_X0_spread.pdf").c_str());
255 
256  // 2D map of X0 for Geantino input
257  TCanvas *GM_2D = new TCanvas("GM_2D","Geantino X0 2D") ;
258  geantino_file->Draw("v_phi:v_eta:t_X0>>geantino_X0","","COLZ");
259  geantino_X0->GetXaxis()->SetTitle("Eta");
260  geantino_X0->GetYaxis()->SetTitle("Phi");
261  geantino_X0->GetZaxis()->SetTitle("X0");
262  GM_2D->Print( (name+"/geant_mat_X0.png").c_str());
263  }
264 
265  if(Val != "" && geantino != ""){
266 
267  Val_X0_Eta->SetMarkerColor(kBlack);
268  Val_X0_Eta->SetLineColor(kBlack);
269  geantino_X0_Eta->SetMarkerColor(kRed);
270  geantino_X0_Eta->SetLineColor(kRed);
271 
272  Val_X0_Phi->SetMarkerColor(kBlack);
273  Val_X0_Phi->SetLineColor(kBlack);
274  geantino_X0_Phi->SetMarkerColor(kRed);
275  geantino_X0_Phi->SetLineColor(kRed);
276 
277  TLegend* leg = new TLegend(0.1,0.15,0.25,0.30);
278  leg->AddEntry(Val_X0_Eta,"Validation X0");
279  leg->AddEntry(geantino_X0_Eta,"Geantino X0");
280 
281  TLegend* leg2 = new TLegend(0.1,0.15,0.25,0.30);
282  leg2->AddEntry(Val_X0_Phi,"Validation X0");
283  leg2->AddEntry(geantino_X0_Phi,"Geantino X0");
284 
285  // X0 differences as function of eta of the Validation and Geantino input
286  comp_X0_Eta->Add(Val_X0_Eta);
287  comp_X0_Eta->Add(geantino_X0_Eta,-1);
288  comp_X0_Eta->Divide(Val_X0_Eta);
289  comp_X0_Eta->Scale(100);
290  comp_X0_Eta->GetXaxis()->SetTitle("Eta");
291  comp_X0_Eta->GetYaxis()->SetTitle("(Validation X0 - geantino X0) / Validation X0 [%]");
292  TCanvas *Diff_Eta = new TCanvas("Diff_Eta","Comparison geantino_Validation Eta") ;
293  comp_X0_Eta->Draw();
294  Diff_Eta->Print( (name+"/Comp_Val_geant_mat_X0_Eta.pdf").c_str());
295 
296  // X0 comparison as function of eta of the Validation and Geantino input
297  TCanvas *Comp_Eta_spread = new TCanvas("Comp_Eta_spread","Comparison geantino_Validation Eta") ;
298  Val_X0_Eta_spread->Draw();
299  geantino_X0_Eta_spread->SetMarkerColor(kRed);
300  geantino_X0_Eta_spread->Draw("SAME");
301  leg->Draw("SAME");
302  Comp_Eta_spread->Print( (name+"/Comp_Val_geant_mat_X0_Eta_spread.pdf").c_str());
303 
304  // X0 differences as function of phi of the Validation and Geantino input
305  comp_X0_Phi->Add(Val_X0_Phi);
306  comp_X0_Phi->Add(geantino_X0_Phi,-1);
307  comp_X0_Phi->Divide(Val_X0_Phi);
308  comp_X0_Phi->Scale(100);
309  comp_X0_Phi->GetXaxis()->SetTitle("Phi");
310  comp_X0_Phi->GetYaxis()->SetTitle("(Validation X0 - geantino X0) / Validation X0 [%]");
311  TCanvas *Diff_Phi = new TCanvas("Diff_Phi","Comparison geantino_Validation") ;
312  comp_X0_Phi->Draw();
313  Diff_Phi->Print( (name+"/Comp_Val_geant_mat_X0_Phi.pdf").c_str());
314 
315  // X0 comparison as function of eta of the Validation and Geantino input
316  TCanvas *Comp_Phi_spread = new TCanvas("Comp_Phi_spread","Comparison geantino_Validation Phi") ;
317  Val_X0_Phi_spread->Draw();
318  geantino_X0_Phi_spread->SetMarkerColor(kRed);
319  geantino_X0_Phi_spread->Draw("SAME");
320  leg2->Draw("SAME");
321  Comp_Phi_spread->Print( (name+"/Comp_Val_geant_mat_X0_Phi_spread.pdf").c_str());
322 
323  Float_t score = 0;
324  for(int i=0; i<Val_X0->GetXaxis()->GetNbins(); i++){
325  score += (Val_X0->GetBinContent(i+1)-geantino_X0->GetBinContent(i+1))*
326  (Val_X0->GetBinContent(i+1)-geantino_X0->GetBinContent(i+1))/
327  geantino_X0->GetBinError(i+1)*geantino_X0->GetBinError(i+1);
328  }
329  std::cout << "Score : " << score << std::endl;
330 
331  TText *t = new TText(0,10, ("Chi^2 score : " + to_string(score)).c_str());
332  t->Draw();
333 
334  // X0 ratio plot as function of eta of the Validation and Geantino input
335  TCanvas *Comp_Eta = new TCanvas("Comp_Eta","Ratio geantino_Validation Eta") ;
336  Val_X0_Eta->SetMarkerStyle(7);
337  geantino_X0_Eta->SetMarkerStyle(7);
338  Val_X0_Eta->SetMarkerColor(kBlack);
339  geantino_X0_Eta->SetMarkerColor(kRed);
340  Val_X0_Eta->SetLineColor(kBlack);
341  geantino_X0_Eta->SetLineColor(kRed);
342 
343  Draw_ratio(Comp_Eta, geantino_X0_Eta, Val_X0_Eta, leg, name);
344 
345  // X0 ratio plot as function of phi of the Validation and Geantino input
346  TCanvas *Comp_Phi = new TCanvas("Comp_Phi","Ratio geantino_Validation Phi") ;
347  Val_X0_Phi->SetMarkerStyle(7);
348  geantino_X0_Phi->SetMarkerStyle(7);
349  Val_X0_Phi->SetMarkerColor(kBlack);
350  geantino_X0_Phi->SetMarkerColor(kRed);
351  Val_X0_Phi->SetLineColor(kBlack);
352  geantino_X0_Phi->SetLineColor(kRed);
353 
354  Draw_ratio(Comp_Phi, geantino_X0_Phi, Val_X0_Phi, leg, name);
355  }
356 
357  return;
358 
359 }