LArSoft  v08_35_01
Liquid Argon Software toolkit - http://larsoft.org/
plot.C
Go to the documentation of this file.
1 // *********************************************************************
2 // To execute this macro under ROOT after your simulation ended,
3 // 1 - launch ROOT (usually type 'root' at your machine's prompt)
4 // 2 - type '.X plot.C' at the ROOT session prompt
5 //
6 // Author: Sebastien Incerti, CNRS, France
7 // Date: 25 Feb. 2015
8 // The Geant4-DNA collaboration
9 // *********************************************************************
10 
11 {
12 
13 gROOT->Reset();
14 gStyle->SetPalette(1);
15 gROOT->SetStyle("Plain");
16 gStyle->SetOptStat(00000);
17 
18 //***************************************
19 //***************************************
20 // MAKE YOUR SELECTION OF B, min and max
21 // for log histograms
22 //***************************************
23 //***************************************
24 
25 Int_t B=1000; // bins per decade
26 Int_t min=-2; // minimum x-axis value as 10^min
27 Int_t max=2; // maximum x-axis value as 10^max
28 
29 //
30 FILE *fp = fopen("yz.root","r");
31 if( fp ) {
32  // exists
33  cout << "*** Notice: the output file yz.root exists ***"<< endl;
34  fclose(fp);
35 } else {
36  cout << "*** Notice: the output file yz.root does not exist ***"<< endl;
37  cout << "*** it will be created from merged ROOT files ***"<< endl;
38  system ("rm -rf yz.root");
39  system ("hadd yz.root yz_*.root");}
40 //
41 
42 c1 = new TCanvas ("c1","",60,60,800,600);
43 c1.Divide(3,2);
44 
45 /*
46 // for testing only
47 FILE * fp = fopen("yz.txt","r");
48 Float_t radius,y,z;
49 Int_t ncols = 0;
50 Int_t nlines = 0;
51 
52 TNtuple *ntuple = new TNtuple("ntuple","micro","radius:y:z");
53 while (1)
54  {
55  ncols = fscanf(fp,"%f %f %f",&radius,&y,&z);
56  if (ncols < 0) break;
57  ntuple->Fill(radius,y,z);
58  nlines++;
59  }
60 fclose(fp);
61 */
62 
63 TFile f("yz.root");
64 
65 TNtuple* ntuple;
66 ntuple = (TNtuple*)f->Get("yz");
67 
68 // --->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->
69 
70 //plot f(y)
71 c1.cd(1);
72 
73 ntuple->Draw("y>>hfy","","");
74 
75 hfy->Scale (1./(hfy->GetEntries()*hfy->GetBinWidth(1))); // DIVIDE BY BIN WIDTH !!!
76 hfy->SetTitle("f(y) (um/keV)");
77 hfy->GetXaxis()->SetTitle("y (keV/um)");
78 hfy->SetFillColor(1);
79 hfy->Draw("");
80 
81 //check normalization
82 
83 Double_t norm=0;
84 
85 for (Int_t j=0;j<hfy->GetNbinsX(); j++)
86  norm=norm+hfy->GetBinContent(j)*hfy->GetBinWidth(1); // MULTIPLY BY BIN WIDTH !!!
87 
88 cout << endl;
89 cout << "**** I - Results from lin-lin histograms ****" << endl;
90 cout << endl;
91 cout << "---> sum of f(y)dy =" << norm << endl;
92 
93 //plot y*f(y)
94 c1.cd(2);
95 
96 TH1F *hyfy = (TH1F*)hfy->Clone("hyfy");
97 
98 for (Int_t i=0;i<hyfy->GetNbinsX(); i++)
99  hyfy->SetBinContent(i,hyfy->GetBinCenter(i)*hyfy->GetBinContent(i));
100 
101 hyfy->SetLineColor(2);
102 hyfy->SetFillColor(2);
103 hyfy->SetTitle("y*f(y)");
104 hyfy->GetXaxis()->SetTitle("y (keV/um)");
105 hyfy->Draw("");
106 
107 // --->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->
108 
109 //calculate yF
110 
111 Double_t yF=0;
112 
113 for (Int_t j=0;j<hyfy->GetNbinsX(); j++)
114  yF=yF+hyfy->GetBinContent(j)*hyfy->GetBinWidth(1); // MULTIPLY BY BIN WIDTH !!!
115 
116 cout << "---> yF=" << yF << " keV/um" << endl;
117 
118 //plot y*f(y)/yF = d(y) (cf. Burigo et al., NIMB 320 (2014) 89-99)
119 c1.cd(4);
120 
121 TH1F *hdy = (TH1F*)hyfy->Clone("hdy");
122 
123 hdy->Scale (1./yF);
124 
125 hdy->SetLineColor(4);
126 hdy->SetFillColor(4);
127 hdy->SetTitle("d(y) (um/keV)");
128 hdy->GetXaxis()->SetTitle("y (keV/um)");
129 hdy->Draw("");
130 
131 //check normalization
132 
133 Double_t normD=0;
134 
135 for (Int_t j=0;j<hdy->GetNbinsX(); j++)
136  normD=normD+hdy->GetBinContent(j)*hdy->GetBinWidth(1); // MULTIPLY BY BIN WIDTH !!!
137 
138 cout << "---> sum of d(y)dy =" << normD << endl;
139 
140 //plot y*d(y)
141 c1.cd(5);
142 
143 TH1F *hydy = (TH1F*)hdy->Clone("hydy");
144 
145 for (Int_t k=0;k<hydy->GetNbinsX(); k++)
146  hydy->SetBinContent(k,hydy->GetBinCenter(k)*hdy->GetBinContent(k));
147 
148 hydy->SetLineColor(3);
149 hydy->SetFillColor(3);
150 hydy->SetTitle("y*d(y)");
151 hydy->GetXaxis()->SetTitle("y (keV/um)");
152 hydy->Draw("");
153 
154 //calculate yD
155 
156 Double_t yD=0;
157 
158 for (Int_t l=0;l<hydy->GetNbinsX(); l++)
159  yD=yD+hydy->GetBinContent(l)*hydy->GetBinWidth(1); // MULTIPLY BY BIN WIDTH !!!
160 
161 cout << "---> yD=" << yD << " keV/um" << endl;
162 
163 // --->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->
164 // --->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->
165 // --->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->
166 
167 //log plot of y*f(y)
168 
169 //goto end;
170 
171 c1.cd(3);
172 
173 cout << endl;
174 cout << "**** II - Results from log-lin histograms ****" << endl;
175 cout << endl;
176 cout << " You have selected "<< B << " bins per decade, from 10^" << min << " (a.u.) to 10^"<< max << " (a.u.) "<< endl;
177 cout << endl;
178 
179 
180 Int_t bins=B*(max-min);
181 
182 TH1F *hlogyfy = new TH1F("1","1",bins,min,max);
183 TAxis *axis = hlogyfy->GetXaxis();
184 
185 TH1F *hlogy2fy = new TH1F("2","2",bins,min,max);
186 TAxis *axis2 = hlogy2fy->GetXaxis();
187 
188 TH1F *hlogydy = new TH1F("3","3",bins,min,max);
189 TAxis *axis3 = hlogydy->GetXaxis();
190 
191 TH1F *hlogy3fy = new TH1F("4","4",bins,min,max);
192 TAxis *axis4 = hlogy3fy->GetXaxis();
193 
194 
195 Axis_t from = axis->GetXmin();
196 Axis_t to = axis->GetXmax();
197 Axis_t width = (to - from) / bins;
198 //cout << "*** width=" << width << endl;
199 //cout << "*** bins=" << bins << endl;
200 Axis_t *new_bins = new Axis_t[bins + 1];
201 for (int i = 0; i <= bins; i++)
202 {
203  new_bins[i] = TMath::Power(10, from + i * width);
204 }
205 
206 axis->Set(bins, new_bins);
207 axis2->Set(bins, new_bins);
208 axis3->Set(bins, new_bins);
209 axis4->Set(bins, new_bins);
210 
211 delete new_bins;
212 
213 //
214 /*
215 //for testing only
216 FILE * fp2 = fopen("yz.txt","r");
217 while (1)
218  {
219  ncols = fscanf(fp2,"%f %f %f",&radius,&y,&z);
220  if (ncols < 0) break;
221  hlogyfy->Fill(y);
222  hlogy2fy->Fill(y);
223  hlogydy->Fill(y);
224  hlogy3fy->Fill(y);
225  hlogzfz->Fill(z);
226  hlogz2fz->Fill(z);
227  hlogzdz->Fill(z);
228  hlogz3fz->Fill(z);
229  nlines++;
230  }
231 fclose(fp2);
232 */
233 //
234 
235 Double_t radius,y,z;
236 
237 ntuple->SetBranchAddress("radius",&radius);
238 ntuple->SetBranchAddress("y",&y);
239 ntuple->SetBranchAddress("z",&z);
240 Int_t nentries = (Int_t)ntuple->GetEntries();
241 for (Int_t i=0; i<nentries; i++)
242 {
243  ntuple->GetEntry(i);
244  hlogyfy->Fill(y);
245  hlogy2fy->Fill(y);
246  hlogydy->Fill(y);
247  hlogy3fy->Fill(y);
248 
249 }
250 
251 //
252 
253 hlogyfy->Scale(1./( hlogyfy->GetEntries() ));
254 hlogy2fy->Scale(1./( hlogy2fy->GetEntries() ));
255 hlogydy->Scale(1./( hlogydy->GetEntries() ));
256 hlogy3fy->Scale(1./( hlogy3fy->GetEntries() ));
257 
258 
259 //plot y*f(y)
260 
261 gPad->SetLogx();
262 
263 for (Int_t i=0;i<hlogyfy->GetNbinsX(); i++)
264 {
265  hlogyfy->SetBinContent(i,
266 //(log(10)/B)*
267 hlogyfy->GetBinCenter(i)
268 *hlogyfy->GetBinContent(i)
269 /(TMath::Power(10, from + (i+1) * width)-TMath::Power(10, from + i * width))
270 );
271 }
272 
273 hlogyfy->SetLineColor(2);
274 hlogyfy->SetFillColor(2);
275 hlogyfy->SetTitle("y*f(y)");
276 hlogyfy->GetXaxis()->SetTitle("y (keV/um)");
277 hlogyfy->Draw("");
278 
279 //check normalization
280 
281 Double_t normLogfy=0;
282 
283 for (Int_t j=0;j<hlogyfy->GetNbinsX(); j++)
284  normLogfy=normLogfy+(log(10)/B)*hlogyfy->GetBinContent(j);
285 
286 cout << "---> sum of Log f(y)dy =" << normLogfy << endl;
287 
288 //calculate yF
289 //requires plot of log y2fy and integration of this plot
290 
291 for (Int_t i=0;i<hlogy2fy->GetNbinsX(); i++)
292 {
293  hlogy2fy->SetBinContent(i,
294 //(log(10)/B)*
295 hlogy2fy->GetBinCenter(i)
296 *hlogy2fy->GetBinCenter(i)
297 *hlogy2fy->GetBinContent(i)
298 /(TMath::Power(10, from + (i+1) * width)-TMath::Power(10, from + i * width))
299 );
300 }
301 
302 Double_t logyF=0;
303 
304 for (Int_t j=0;j<hlogy2fy->GetNbinsX(); j++)
305  logyF=logyF+(log(10)/B)*hlogy2fy->GetBinContent(j);
306 
307 cout << "---> yF =" << logyF << " keV/um" << endl;
308 
309 // --->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->
310 
311 //plot y*d(y)
312 
313 c1.cd(6);
314 
315 for (Int_t i=0;i<hlogydy->GetNbinsX(); i++)
316 {
317  hlogydy->SetBinContent(i,
318 //(log(10)/B)*
319 hlogydy->GetBinCenter(i)
320 *hlogydy->GetBinCenter(i)
321 *hlogydy->GetBinContent(i)
322 /(TMath::Power(10, from + (i+1) * width)-TMath::Power(10, from + i * width))
323 /logyF
324 );
325 }
326 
327 hlogydy->SetLineColor(3);
328 hlogydy->SetFillColor(3);
329 hlogydy->SetTitle("y*d(y)");
330 hlogydy->GetXaxis()->SetTitle("y (keV/um)");
331 gPad->SetLogx();
332 hlogydy->Draw("");
333 
334 //check normalization
335 
336 Double_t normLogdy=0;
337 
338 for (Int_t j=0;j<hlogydy->GetNbinsX(); j++)
339  normLogdy=normLogdy+(log(10)/B)*hlogydy->GetBinContent(j);
340 
341 cout << "---> sum of Log d(y)dy =" << normLogdy << endl;
342 
343 //calculate yD
344 //requires plot of log y3fy and integration of this plot
345 
346 for (Int_t i=0;i<hlogy3fy->GetNbinsX(); i++)
347 {
348  hlogy3fy->SetBinContent(i,
349 //(log(10)/B)*
350 hlogy3fy->GetBinCenter(i)
351 *hlogy3fy->GetBinCenter(i)
352 *hlogy3fy->GetBinCenter(i)
353 *hlogy3fy->GetBinContent(i)
354 /(TMath::Power(10, from + (i+1) * width)-TMath::Power(10, from + i * width))
355 /logyF
356 );
357 }
358 
359 Double_t logyD=0;
360 
361 for (Int_t j=0;j<hlogy3fy->GetNbinsX(); j++)
362  logyD=logyD+(log(10)/B)*hlogy3fy->GetBinContent(j);
363 
364 cout << "---> yD =" << logyD << " keV/um" << endl;
365 
366 // --->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->--->
367 
368 //
369 end:
370 cout << endl;
371 }
372 
Double_t y
Definition: plot.C:279
system("rm -rf microbeam.root")
TFile f("microbeam.root")
Double_t norm
Definition: plot.C:83
Int_t B
Definition: plot.C:25
Double_t z
Definition: plot.C:279
TNtuple * ntuple
Definition: plot.C:20
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:72
Int_t max
Definition: plot.C:27
Double_t normLogfy
Definition: plot.C:281
Int_t nentries
Definition: plot.C:283
FILE * fp
Definition: plot.C:36
Double_t radius
Definition: plot.C:235
fclose(fp)
Int_t min
Definition: plot.C:26
c1
Definition: plot.C:28
delete new_bins
Definition: plot.C:211
Double_t logyF
Definition: plot.C:302