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