| Geant4 Cross Reference |
1 // ******************************************* 1 // *********************************************************************
2 // To execute this macro under ROOT after your << 2 // To execute this macro under ROOT after your simulation ended,
3 // 1 - launch ROOT (usually type 'root' at y 3 // 1 - launch ROOT (usually type 'root' at your machine's prompt)
4 // 2 - type '.X plot.C' at the ROOT session 4 // 2 - type '.X plot.C' at the ROOT session prompt
5 // 5 //
6 // Author: Sebastien Incerti 6 // Author: Sebastien Incerti
7 // Date: March 2nd, 2019 7 // Date: March 2nd, 2019
8 // The Geant4-DNA collaboration 8 // The Geant4-DNA collaboration
9 // ******************************************* 9 // *********************************************************************
10 10
11 void SetLeafAddress(TNtuple* ntuple, const cha 11 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address);
12 12
13 void plot() 13 void plot()
14 { 14 {
15 gROOT->Reset(); 15 gROOT->Reset();
16 gStyle->SetPalette(1); 16 gStyle->SetPalette(1);
17 gROOT->SetStyle("Plain"); 17 gROOT->SetStyle("Plain");
18 gStyle->SetOptStat(00000); 18 gStyle->SetOptStat(00000);
19 19
20 //******************** 20 //********************
21 Int_t nbRadius = 101; 21 Int_t nbRadius = 101;
22 //******************** 22 //********************
23 23
24 TCanvas *c1 = new TCanvas ("c1","",60,60,800,8 24 TCanvas *c1 = new TCanvas ("c1","",60,60,800,800);
25 Int_t mycolor; 25 Int_t mycolor;
26 26
27 TFile f("t.root"); << 27 TFile f("t.root");
28 mycolor=4; 28 mycolor=4;
29 29
30 TNtuple* ntuple; 30 TNtuple* ntuple;
31 ntuple = (TNtuple*)f.Get("t"); << 31 ntuple = (TNtuple*)f.Get("t");
32 32
33 bool rowWise = true; 33 bool rowWise = true;
34 TBranch* eventBranch = ntuple->FindBranch("row 34 TBranch* eventBranch = ntuple->FindBranch("row_wise_branch");
35 if ( ! eventBranch ) rowWise = false; 35 if ( ! eventBranch ) rowWise = false;
>> 36 // std::cout << "rowWise: " << rowWise << std::endl
36 37
37 Double_t radius1,nofHits,nbEdep,edep,radius2,E 38 Double_t radius1,nofHits,nbEdep,edep,radius2,Einc;
38 Int_t noRadius; 39 Int_t noRadius;
39 40
40 if ( ! rowWise ) { 41 if ( ! rowWise ) {
41 ntuple->SetBranchAddress("radius1",&radius1) 42 ntuple->SetBranchAddress("radius1",&radius1);
42 ntuple->SetBranchAddress("noRadius",&noRadiu 43 ntuple->SetBranchAddress("noRadius",&noRadius);
43 ntuple->SetBranchAddress("nbHits",&nofHits); 44 ntuple->SetBranchAddress("nbHits",&nofHits);
44 ntuple->SetBranchAddress("nbScoredHits",&nbE 45 ntuple->SetBranchAddress("nbScoredHits",&nbEdep);
45 ntuple->SetBranchAddress("edep",&edep); 46 ntuple->SetBranchAddress("edep",&edep);
46 ntuple->SetBranchAddress("radius2",&radius2) 47 ntuple->SetBranchAddress("radius2",&radius2);
47 ntuple->SetBranchAddress("Einc",&Einc); 48 ntuple->SetBranchAddress("Einc",&Einc);
48 } 49 }
49 else { 50 else {
50 SetLeafAddress(ntuple, "radius1",&radius1); 51 SetLeafAddress(ntuple, "radius1",&radius1);
51 SetLeafAddress(ntuple, "noRadius",&noRadius) 52 SetLeafAddress(ntuple, "noRadius",&noRadius);
52 SetLeafAddress(ntuple, "nbHits",&nofHits); 53 SetLeafAddress(ntuple, "nbHits",&nofHits);
53 SetLeafAddress(ntuple, "nbScoredHits",&nbEde 54 SetLeafAddress(ntuple, "nbScoredHits",&nbEdep);
54 SetLeafAddress(ntuple, "edep",&edep); 55 SetLeafAddress(ntuple, "edep",&edep);
55 SetLeafAddress(ntuple, "radius2",&radius2); 56 SetLeafAddress(ntuple, "radius2",&radius2);
56 SetLeafAddress(ntuple, "Einc",&Einc); << 57 SetLeafAddress(ntuple, "Einc",&Einc);
57 } 58 }
58 59
59 Int_t nentries = (Int_t)ntuple->GetEntries(); 60 Int_t nentries = (Int_t)ntuple->GetEntries();
60 61
61 // 62 //
62 63
63 Double_t t[1000]; // 1000 is the max number of 64 Double_t t[1000]; // 1000 is the max number of radius values
64 Double_t population[1000]; 65 Double_t population[1000];
65 Double_t myRad[1000]; << 66 Double_t myRad[1000];
66 67
67 for (Int_t i=0; i<1000; i++) << 68 for (Int_t i=0; i<1000; i++)
68 { 69 {
69 t[i]=0; << 70 t[i]=0;
70 population[i]=0; 71 population[i]=0;
71 myRad[i]=0; << 72 myRad[i]=0;
72 } 73 }
73 74
74 Int_t event = 0; 75 Int_t event = 0;
75 76
76 for (Int_t i=0; i<nentries; i++) << 77 for (Int_t i=0; i<nentries; i++)
77 { 78 {
78 ntuple->GetEntry(i); 79 ntuple->GetEntry(i);
79 t[noRadius] = t[noRadius] + edep; 80 t[noRadius] = t[noRadius] + edep;
80 population[noRadius]=population[noRadius]+1; 81 population[noRadius]=population[noRadius]+1;
81 myRad[noRadius] = radius1; << 82 myRad[noRadius] = radius1;
82 } 83 }
83 84
84 // Mean 85 // Mean
85 86
86 for (Int_t j=1; j<nbRadius; j++) << 87 for (Int_t j=1; j<nbRadius; j++)
87 { 88 {
88 t[j] = t[j]/population[j]; 89 t[j] = t[j]/population[j];
89 t[j] = t[j]/(myRad[j+1]-myRad[j]); 90 t[j] = t[j]/(myRad[j+1]-myRad[j]);
90 //cout << j << " " << myRad[j] << " " << myR << 91 //cout << j << " " << myRad[j] << " " << myRad[j+1]
91 // << " " << t[j] << " " << population[j] << 92 // << " " << t[j] << " " << population[j] << endl;
92 } 93 }
93 94
94 // 95 //
95 96
96 c1->cd(1); 97 c1->cd(1);
97 98
98 TGraph* gr1 =new TGraph(nbRadius,myRad,t); 99 TGraph* gr1 =new TGraph(nbRadius,myRad,t);
99 gr1->SetMarkerColor(2); 100 gr1->SetMarkerColor(2);
100 gr1->SetMarkerStyle(20); 101 gr1->SetMarkerStyle(20);
101 gr1->SetMarkerSize(1); 102 gr1->SetMarkerSize(1);
102 gr1->SetLineColor(2); 103 gr1->SetLineColor(2);
103 gr1->SetTitle(""); 104 gr1->SetTitle("");
104 gr1->GetXaxis()->SetLimits(0.1,100); 105 gr1->GetXaxis()->SetLimits(0.1,100);
105 gr1->GetYaxis()->SetLimits(0.,22); 106 gr1->GetYaxis()->SetLimits(0.,22);
106 gr1->GetXaxis()->SetLabelSize(0.025); 107 gr1->GetXaxis()->SetLabelSize(0.025);
107 gr1->GetYaxis()->SetLabelSize(0.025); 108 gr1->GetYaxis()->SetLabelSize(0.025);
108 gr1->GetXaxis()->SetTitleSize(0.035); 109 gr1->GetXaxis()->SetTitleSize(0.035);
109 gr1->GetYaxis()->SetTitleSize(0.035); 110 gr1->GetYaxis()->SetTitleSize(0.035);
110 gr1->GetXaxis()->SetTitleOffset(1.4); 111 gr1->GetXaxis()->SetTitleOffset(1.4);
111 gr1->GetYaxis()->SetTitleOffset(1.4); 112 gr1->GetYaxis()->SetTitleOffset(1.4);
112 gr1->GetXaxis()->SetTitle("r (nm)"); 113 gr1->GetXaxis()->SetTitle("r (nm)");
113 gr1->GetYaxis()->SetTitle("t (eV/nm)"); 114 gr1->GetYaxis()->SetTitle("t (eV/nm)");
114 gr1->Draw(""); 115 gr1->Draw("");
115 gPad->SetLogx(); 116 gPad->SetLogx();
116 117
117 } 118 }
118 119
119 void SetLeafAddress(TNtuple* ntuple, const cha 120 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address) {
120 TLeaf* leaf = ntuple->FindLeaf(name); 121 TLeaf* leaf = ntuple->FindLeaf(name);
121 if ( ! leaf ) { 122 if ( ! leaf ) {
122 std::cerr << "Error in <SetLeafAddress>: u 123 std::cerr << "Error in <SetLeafAddress>: unknown leaf --> " << name << std::endl;
123 return; 124 return;
124 } 125 }
125 leaf->SetAddress(address); 126 leaf->SetAddress(address);
126 } 127 }
>> 128
127 129