| Geant4 Cross Reference |
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 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address);
12
13 void plot()
14 {
15 gROOT->Reset();
16 gStyle->SetPalette(1);
17 gROOT->SetStyle("Plain");
18 gStyle->SetOptStat(00000);
19
20 //***************************************
21 //***************************************
22 // MAKE YOUR SELECTIONS
23 // for histograms
24 //***************************************
25 //***************************************
26
27 Int_t linB=100; // linear histo: nb of bins in x - 1000 is best for integration
28 Double_t ymin=1; // minimum x-axis value
29 Double_t ymax=300; // maximum x-axis value
30
31 //***************************************
32 //***************************************
33
34 //Uncomment for manual merging of histograms
35 //system ("rm -rf yz.root");
36 //system ("hadd yz.root yz_*.root");
37
38 TCanvas* c1 = new TCanvas ("c1","",60,60,800,800);
39 Int_t mycolor;
40
41 TFile* f = new TFile("yz.root");
42 mycolor=4;
43
44 TNtuple* ntuple;
45 ntuple = (TNtuple*)f->Get("yz");
46 bool rowWise = true;
47 TBranch* eventBranch = ntuple->FindBranch("row_wise_branch");
48 if ( ! eventBranch ) rowWise = false;
49 // std::cout << "rowWise: " << rowWise << std::endl;
50
51 Double_t radius,eventID,nofHits,nbEdep,y,z,Einc;
52 if ( ! rowWise ) {
53 ntuple->SetBranchAddress("radius",&radius);
54 ntuple->SetBranchAddress("eventID",&eventID);
55 ntuple->SetBranchAddress("nbHits",&nofHits);
56 ntuple->SetBranchAddress("nbScoredHits",&nbEdep);
57 ntuple->SetBranchAddress("y",&y);
58 ntuple->SetBranchAddress("z",&z);
59 ntuple->SetBranchAddress("Einc",&Einc);
60 }
61 else {
62 SetLeafAddress(ntuple, "radius",&radius);
63 SetLeafAddress(ntuple, "eventID",&eventID);
64 SetLeafAddress(ntuple, "nbHits",&nofHits);
65 SetLeafAddress(ntuple, "nbScoredHits",&nbEdep);
66 SetLeafAddress(ntuple, "y",&y);
67 SetLeafAddress(ntuple, "z",&z);
68 SetLeafAddress(ntuple, "Einc",&Einc);
69 }
70
71 //plot f(y)
72
73 c1->cd(1);
74
75 TH1F *hfyw = new TH1F ("hfyw","hfyw",linB,0,ymax);
76
77 Int_t nentries = (Int_t)ntuple->GetEntries();
78 Double_t population=0;
79 Double_t yLocalMin=1e100;
80 Double_t yLocalMax=0;
81
82 Double_t yF_anal=0;
83 Double_t yD_anal=0;
84
85 for (Int_t i=0; i<nentries; i++)
86 {
87 ntuple->GetEntry(i);
88
89 hfyw->Fill(y,nofHits/nbEdep);
90 if (yLocalMin>y) yLocalMin=y;
91 if (yLocalMax<y) yLocalMax=y;
92 population=population+nofHits/nbEdep;
93 yF_anal = yF_anal + (nofHits/nbEdep)*y;
94 yD_anal = yD_anal + (nofHits/nbEdep)*y*y;
95 }
96
97 cout << "**** Results ****" << endl;
98 cout << endl;
99 cout << "---> yF =" << yF_anal/population << " keV/um" << endl;
100 cout << "---> yD =" << (yD_anal/population)/(yF_anal/population) << " keV/um" << endl;
101 cout << endl;
102 cout << "---> Limits: " << endl;
103 cout << " * min value of y = " << yLocalMin << " keV/um" << endl;
104 cout << " * max value of y = " << yLocalMax << " keV/um" << endl;
105
106 if ( (yLocalMax>ymax) || (yLocalMin<ymin) )
107 {
108 cout << "WARNING: please check your histogram limits ! " << endl;
109 }
110
111 gPad->SetLogy();
112 hfyw->Scale (1./(population*hfyw->GetBinWidth(1)));
113 hfyw->SetTitle("f(y) (um/keV)");
114 hfyw->GetXaxis()->SetTitle("y (keV/um)");
115 hfyw->SetFillColor(2);
116 hfyw->SetLineColor(2);
117 hfyw->Draw("HIST");
118 }
119
120 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address) {
121 TLeaf* leaf = ntuple->FindLeaf(name);
122 if ( ! leaf ) {
123 std::cerr << "Error in <SetLeafAddress>: unknown leaf --> " << name << std::endl;
124 return;
125 }
126 leaf->SetAddress(address);
127 }
128