Geant4 Cross Reference |
1 // ******************************************* 1 2 // To execute this macro under ROOT after your 3 // 1 - launch ROOT (usually type 'root' at y 4 // 2 - type '.X plot.C' at the ROOT session 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 cha 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 bi 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 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("r 48 if ( ! eventBranch ) rowWise = false; 49 // std::cout << "rowWise: " << rowWise << s 50 51 Double_t radius,eventID,nofHits,nbEdep,y,z,E 52 if ( ! rowWise ) { 53 ntuple->SetBranchAddress("radius",&radius) 54 ntuple->SetBranchAddress("eventID",&eventI 55 ntuple->SetBranchAddress("nbHits",&nofHits 56 ntuple->SetBranchAddress("nbScoredHits",&n 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",&nbE 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, 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 << 100 cout << "---> yD =" << (yD_anal/population)/ 101 cout << endl; 102 cout << "---> Limits: " << endl; 103 cout << " * min value of y = " << yLocal 104 cout << " * max value of y = " << yLocal 105 106 if ( (yLocalMax>ymax) || (yLocalMin<ymin) ) 107 { 108 cout << "WARNING: please check your histog 109 } 110 111 gPad->SetLogy(); 112 hfyw->Scale (1./(population*hfyw->GetBinWidt 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 cha 121 TLeaf* leaf = ntuple->FindLeaf(name); 122 if ( ! leaf ) { 123 std::cerr << "Error in <SetLeafAddress>: u 124 return; 125 } 126 leaf->SetAddress(address); 127 } 128