dna/microyz/plot.C

0001 // *********************************************************************
0002 // To execute this macro under ROOT after your simulation ended,
0003 //   1 - launch ROOT (usually type 'root' at your machine's prompt)
0004 //   2 - type '.X plot.C' at the ROOT session prompt
0005 //
0006 // Author: Sebastien Incerti, CNRS, France
0007 // Date: 25 Feb. 2015
0008 // The Geant4-DNA collaboration
0009 // *********************************************************************
0010
0011 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address);
0012
0013 void plot()
0014 {
0015   gROOT->Reset();
0016   gStyle->SetPalette(1);
0017   gROOT->SetStyle("Plain");
0018   gStyle->SetOptStat(00000);
0019
0020   //***************************************
0021   //***************************************
0022   // MAKE YOUR SELECTIONS
0023   // for histograms
0024   //***************************************
0025   //***************************************
0026
0027   Int_t linB=100;    // linear histo: nb of bins in x - 1000 is best for integration
0028   Double_t ymin=1;   // minimum x-axis value
0029   Double_t ymax=300; // maximum x-axis value
0030
0031   //***************************************
0032   //***************************************
0033
0034   //Uncomment for manual merging of histograms
0035   //system ("rm -rf yz.root");
0036   //system ("hadd yz.root yz_*.root");
0037
0038   TCanvas* c1 = new TCanvas ("c1","",60,60,800,800);
0039   Int_t mycolor;
0040
0041   TFile* f = new TFile("yz.root");
0042   mycolor=4;
0043
0044   TNtuple* ntuple;
0045   ntuple = (TNtuple*)f->Get("yz");
0046   bool rowWise = true;
0047   TBranch* eventBranch = ntuple->FindBranch("row_wise_branch");
0048   if ( ! eventBranch ) rowWise = false;
0049   // std::cout <<  "rowWise: " << rowWise << std::endl;
0050
0051   Double_t radius,eventID,nofHits,nbEdep,y,z,Einc;
0052   if ( ! rowWise ) {
0053     ntuple->SetBranchAddress("radius",&radius);
0054     ntuple->SetBranchAddress("eventID",&eventID);
0055     ntuple->SetBranchAddress("nbHits",&nofHits);
0056     ntuple->SetBranchAddress("nbScoredHits",&nbEdep);
0057     ntuple->SetBranchAddress("y",&y);
0058     ntuple->SetBranchAddress("z",&z);
0059     ntuple->SetBranchAddress("Einc",&Einc);
0060   }
0061   else {
0062     SetLeafAddress(ntuple, "radius",&radius);
0063     SetLeafAddress(ntuple, "eventID",&eventID);
0064     SetLeafAddress(ntuple, "nbHits",&nofHits);
0065     SetLeafAddress(ntuple, "nbScoredHits",&nbEdep);
0066     SetLeafAddress(ntuple, "y",&y);
0067     SetLeafAddress(ntuple, "z",&z);
0068     SetLeafAddress(ntuple, "Einc",&Einc);
0069   }
0070
0071   //plot f(y)
0072
0073   c1->cd(1);
0074
0075   TH1F *hfyw = new TH1F ("hfyw","hfyw",linB,0,ymax);
0076
0077   Int_t nentries = (Int_t)ntuple->GetEntries();
0078   Double_t population=0;
0079   Double_t yLocalMin=1e100;
0080   Double_t yLocalMax=0;
0081
0082   Double_t yF_anal=0;
0083   Double_t yD_anal=0;
0084
0085   for (Int_t i=0; i<nentries; i++)
0086   {
0087    ntuple->GetEntry(i);
0088
0089    hfyw->Fill(y,nofHits/nbEdep);
0090    if (yLocalMin>y) yLocalMin=y;
0091    if (yLocalMax<y) yLocalMax=y;
0092    population=population+nofHits/nbEdep;
0093    yF_anal = yF_anal + (nofHits/nbEdep)*y;
0094    yD_anal = yD_anal + (nofHits/nbEdep)*y*y;
0095   }
0096
0097   cout << "**** Results ****" << endl;
0098   cout << endl;
0099   cout << "---> yF =" << yF_anal/population << " keV/um" << endl;
0100   cout << "---> yD =" << (yD_anal/population)/(yF_anal/population) << " keV/um" << endl;
0101   cout << endl;
0102   cout << "---> Limits: " << endl;
0103   cout << "     * min value of y = " << yLocalMin << " keV/um" << endl;
0104   cout << "     * max value of y = " << yLocalMax << " keV/um" << endl;
0105
0106   if ( (yLocalMax>ymax) || (yLocalMin<ymin) )
0107   {
0108     cout << "WARNING: please check your histogram limits ! " << endl;
0109   }
0110
0111   gPad->SetLogy();
0112   hfyw->Scale (1./(population*hfyw->GetBinWidth(1)));
0113   hfyw->SetTitle("f(y) (um/keV)");
0114   hfyw->GetXaxis()->SetTitle("y (keV/um)");
0115   hfyw->SetFillColor(2);
0116   hfyw->SetLineColor(2);
0117   hfyw->Draw("HIST");
0118 }
0119
0120 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address) {
0121   TLeaf* leaf = ntuple->FindLeaf(name);
0122   if ( ! leaf ) {
0123     std::cerr << "Error in <SetLeafAddress>: unknown leaf --> " << name << std::endl;
0124     return;
0125   }
0126   leaf->SetAddress(address);
0127 }