alpharad_new: z_analysis/readtree.c Source File

00001 
00002 /*
00003 * 
00004 * Copyright: Mathieu Trocme, mathieu.trocme@gmail.com (21 september 2009)
00005 * 
00006 * This software is a computer program whose purpose is to: simulate the AlphaRad device. 
00007 * 
00008 * This software is governed by the CeCILL license under French law and abiding by the rules of distribution of free software. 
00009 * You can use, modify and/or redistribute the software under the terms of the CeCILL license as circulated by CEA, CNRS and INRIA at the following URL "http://www.cecill.info".
00010 * As a counterpart to the access to the source code and rights to copy, modify and redistribute granted by the license, users are provided only with a limited warranty and the software's author, the holder of the economic rights, and the successive licensors have only limited liability. 
00011 * In this respect, the user's attention is drawn to the risks associated with loading, using, modifying and/or developing or reproducing the software by the user in light of its specific status of free software, that may mean that it is complicated to manipulate, and that also therefore means that it is reserved for developers and experienced professionals having in-depth computer knowledge. 
00012 * Users are therefore encouraged to load and test the software's suitability as regards their requirements in conditions enabling the security of their systems and/or data to be ensured and, more generally, to use and operate it in the same conditions as regards security. 
00013 * 
00014 * The fact that you are presently reading this means that you have had knowledge of the CeCILL license and that you accept its terms.
00015 * 
00016 */
00017 
00018 
00019 // 
00020 // Program       : alphaRad_new.cc
00021 // Task          : Analysis of simulation of the alphaRad device
00022 // Location      : $G4MYSIMU_DIR/z_analysis (See https://sbgtrac.in2p3.fr/projects/ramses/wiki/G4_MyTuto for the list of environment variables used)
00023 // Compilation   : '.L TMyRootEvent.hpp+' then '.L readtree.c+'
00024 // Execution     : 'readtree("myFileName")'
00025 // Output        : plots...
00026 // Documentation : See $G4MYSIMU_PATH/doc.html
00027 // Convention    : Root coding convention: kMyVariable for constants, gMyVariable for global variable, fMyVariable for class data member
00028 // Help          : http://www.lcsim.org/software/geant4/doxygen/html/classes.html
00029 // Author        : Mathieu Trocme (mathieu.trocme@ires.in2p3.fr)
00030 // Date          : 21 September 2009
00031 // Framework     : PhD thesis, CNRS/IPHC/DRS/RAMSES, Strasbourg, France
00032 // 
00033 
00034 
00035 #include <TROOT.h>
00036 #include <TGClient.h>
00037 #include <TSystem.h>
00038 #include <TFile.h>
00039 #include <TTree.h>
00040 #include <TParameter.h>
00041 #include <TObjString.h>
00042 #include <TH1.h>
00043 #include <TH2.h>
00044 #include <TCanvas.h>
00045 #include <TPad.h>
00046 #include <TPaveStats.h>
00047 #include <TStyle.h>
00048 #include <TLatex.h>
00049 #include <TLine.h>
00050 
00051 #include "TMyRootEvent.hpp"
00052 
00053 
00054 
00055 void Plot ( TH1D** histo, TString ytitle, TString label[], Int_t nCells, Int_t nAlphas_in, Bool_t noDummyCell, Bool_t displayDose, Double_t *integral_dose, Double_t deltaEdep ) 
00056 {
00057    Double_t ymax=-1., sizeX=0.20, sizeY=0.16 ; 
00058    TPaveStats** stats = new TPaveStats* [nCells] ; TLatex** tex_dose  = new TLatex* [nCells] ; TLine** line = new TLine* [nCells] ;
00059    /// Scaling & getting the max
00060    for (Int_t cellNumber=0; cellNumber<nCells; cellNumber++) 
00061    { 
00062       histo[cellNumber]->Scale(1./(nAlphas_in/100.)) ; 
00063       if ( histo[cellNumber]->Integral() != 0. && histo[cellNumber]->GetMaximum() > ymax ) { ymax = histo[cellNumber]->GetMaximum() ; }
00064    }
00065    ymax *= 1.1 ; 
00066    for (Int_t cellNumber=0; cellNumber<nCells; cellNumber++) 
00067    {
00068       histo[cellNumber]->SetName(label[cellNumber]) ; 
00069       histo[cellNumber]->SetTitle("") ; 
00070       histo[cellNumber]->SetXTitle(ytitle                    ) ; histo[cellNumber]->GetXaxis()->CenterTitle() ; 
00071       histo[cellNumber]->SetYTitle("Absolute efficiency [\%]") ; histo[cellNumber]->GetYaxis()->CenterTitle() ; 
00072       histo[cellNumber]->SetLineColor(cellNumber+1) ; histo[cellNumber]->SetLineStyle(cellNumber+1) ; 
00073       histo[cellNumber]->SetMaximum(ymax) ; 
00074       histo[cellNumber]->Draw((cellNumber==0)?"":"sames") ; 
00075       if ( displayDose ) {
00076          tex_dose[cellNumber] = new TLatex ( deltaEdep*0.25 , ymax*(0.90-0.06*cellNumber) , TString(Form("Total dose = %.2e Gy",integral_dose[cellNumber]))) ; tex_dose[cellNumber]->SetTextColor(cellNumber+1) ; tex_dose[cellNumber]->SetTextSize(0.035) ; tex_dose[cellNumber]->Draw() ; 
00077          line[cellNumber] = new TLine ( deltaEdep*0.15 , ymax*(0.915-0.06*cellNumber) , deltaEdep*0.23 , ymax*(0.915-0.06*cellNumber) ) ; line[cellNumber]->SetLineColor(cellNumber+1) ; line[cellNumber]->SetLineStyle(cellNumber+1) ; line[cellNumber]->SetLineWidth(2) ; line[cellNumber]->Draw() ; 
00078       }
00079       gPad->Update() ; 
00080       stats[cellNumber] = new TPaveStats() ; stats[cellNumber] = (TPaveStats*)histo[cellNumber]->GetListOfFunctions()->FindObject("stats") ;
00081       stats[cellNumber]->SetTextColor(cellNumber+1) ; stats[cellNumber]->SetLineColor(cellNumber+1) ; stats[cellNumber]->SetLineStyle(cellNumber+1) ; stats[cellNumber]->SetBorderSize(1) ; stats[cellNumber]->SetX1NDC(0.68) ; stats[cellNumber]->SetY1NDC(0.70-(noDummyCell?0.08:0.)-(sizeY+0.025)*cellNumber) ; stats[cellNumber]->SetX2NDC(stats[cellNumber]->GetX1NDC()+sizeX) ; stats[cellNumber]->SetY2NDC(stats[cellNumber]->GetY1NDC()+sizeY) ; 
00082       gPad->Modified() ;
00083    }
00084    return ;
00085 }
00086 
00087 
00088 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
00089 // MAIN
00090 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
00091 
00092 
00093 void readtree ( TString finName="data/alphaRad_1e+03_alpha_t=000mum_r=001mm_dsd=020mm.root" )
00094 {
00095 
00096    gROOT->ProcessLine(".L TMyRootEvent.hpp+") ;
00097 
00098    TH1::AddDirectory(kFALSE) ;
00099    Bool_t noDummyCell = kTRUE ; 
00100    TString finName_raw = finName ; finName_raw.ReplaceAll("data/","") ; finName_raw.ReplaceAll(".root","") ; 
00101 
00102 ////////////////////////////////////////////////////////////
00103 // Getting back the tree
00104 ////////////////////////////////////////////////////////////
00105 
00106    TFile *fin = new TFile (finName,"read") ; 
00107    TTree *tree = (TTree*)fin->Get("tree") ;
00108 
00109    TMyRootEvent* myEvent = new TMyRootEvent() ;
00110    tree->SetBranchAddress("evt",&myEvent) ; 
00111 
00112    Int_t nParticles = Int_t( tree->GetEntries() ) ;
00113    fprintf (stdout,"\nnParticles = %d",nParticles) ;  
00114    Int_t nAlphas_initial = ((TParameter<long>*)(tree->GetUserInfo()->FindObject("nInitialParticles")))->GetVal() ;  // nAlpha_initial = nEvents
00115    Int_t nAlphas_inSiO2=0, nElectrons_inSiO2=0, nAlphas_inSiepi=0, nElectrons_inSiepi=0, nAlphas_inSisub=0, nElectrons_inSisub=0 ;
00116    for (Int_t particle_idx=0; particle_idx<nParticles; particle_idx++) {
00117       tree->GetEntry(particle_idx) ; 
00118       if ( myEvent->GetParticleName() == "alpha" ) { if ( myEvent->GetEkin("alphaRad_Sensor1_SiO2_SD")>0. ) { nAlphas_inSiO2   ++ ; } if ( myEvent->GetEkin("alphaRad_Sensor1_Siepi_SD")>0. ) { nAlphas_inSiepi   ++ ; } if ( myEvent->GetEkin("alphaRad_Sensor1_Sisub_SD")>0. ) { nAlphas_inSisub   ++ ; } }
00119       if ( myEvent->GetParticleName() == "e-"    ) { if ( myEvent->GetEkin("alphaRad_Sensor1_SiO2_SD")>0. ) { nElectrons_inSiO2++ ; } if ( myEvent->GetEkin("alphaRad_Sensor1_Siepi_SD")>0. ) { nElectrons_inSiepi++ ; } if ( myEvent->GetEkin("alphaRad_Sensor1_Sisub_SD")>0. ) { nElectrons_inSisub++ ; } }
00120    }
00121    fprintf (stdout,"\nnAlpha_in = %d \nnParticles_inSiO2 = %d ?=? %d alpha + %d e- \nnParticles_inSiepi = %d alpha + %d e-\nnParticles_inSisub = %d alpha + %d e-\n",nAlphas_initial,nParticles, nAlphas_inSiO2,nElectrons_inSiO2, nAlphas_inSiepi,nElectrons_inSiepi, nAlphas_inSisub,nElectrons_inSisub) ;
00122 
00123    tree->GetEntry(0) ;
00124    Int_t nCells = myEvent->GetNCells() ;
00125    if ( noDummyCell ) { nCells -=  1 ; }
00126    fprintf (stdout,"\nnCells = %d\n",nCells) ; 
00127 
00128 
00129 ////////////////////////////////////////////////////////////
00130 // Creating histograms
00131 ////////////////////////////////////////////////////////////
00132 
00133    TLatex** tex_label = new TLatex* [nCells] ; 
00134    TString label [] = { "#alpha entering the SiO2 layer", "#alpha entering the Si-epi layer", "#alpha entering the Si-sub layer", "#alpha escaping the Si-sub layer" } ;
00135 
00136    TH1D** h_Ekin  = new TH1D* [nCells] ; Double_t Ekin_min  = 0.       , Ekin_max  = 10.      , Ekin_resolution  = 0.05            ; Int_t Ekin_nBins  = Int_t ( (Ekin_max-Ekin_min)   / Ekin_resolution  ) ;  // [MeV]
00137    TH1D** h_Theta = new TH1D* [nCells] ; Double_t Theta_min = -10.     , Theta_max = 110.     , Theta_resolution = 1.              ; Int_t Theta_nBins = Int_t ( (Theta_max-Theta_min) / Theta_resolution ) ;  // [#circ]
00138    TH1D** h_Edep  = new TH1D* [nCells] ; Double_t Edep_min  = Ekin_min , Edep_max  = Ekin_max , Edep_resolution  = Ekin_resolution ; Int_t Edep_nBins  = Int_t ( (Edep_max-Edep_min)   / Edep_resolution  ) ;  // [keV]
00139    TH1D** h_Dose  = new TH1D* [nCells] ; Double_t Dose_min  = 0.       , Dose_max  = 1e-6     , Dose_resolution  = 1e-8            ; Int_t Dose_nBins  = Int_t ( (Dose_max-Dose_min)   / Dose_resolution  ) ;  // [Gy]
00140    TH1D** h_Range = new TH1D* [nCells] ; Double_t Range_min = 0.       , Range_max = 40.      , Range_resolution = 1.              ; Int_t Range_nBins = Int_t ( (Range_max-Range_min) / Range_resolution ) ;  // [mum]
00141    TH2D** h_XY    = new TH2D* [nCells] ; Double_t XY_min    = 0.       , XY_max    = 64.      , XY_resolution    = 1.              ; Int_t XY_nBins    = Int_t ( (XY_max-XY_min)       / XY_resolution    ) ;  // [pix]
00142    TH1D** h_Z     = new TH1D* [nCells] ; Double_t Z_min     = -1.      , Z_max     = 30.      , Z_resolution     = 0.1             ; Int_t Z_nBins     = Int_t ( (Z_max-Z_min)         / Z_resolution     ) ;  // [mum]
00143    Double_t pix_x = 80., pix_y = 80. ;  // pixel size [mum]
00144 
00145    for (Int_t cellNumber=0; cellNumber<nCells; cellNumber++)
00146    {
00147       tex_label [cellNumber] = new TLatex (0.10,0.50,label[cellNumber]) ; tex_label[cellNumber]->SetTextSize(0.065) ; 
00148       h_Ekin    [cellNumber] = new TH1D ( TString(Form("h_Ekin_%02d" ,cellNumber)) , TString(Form("h_Ekin_%02d" ,cellNumber)) , Ekin_nBins , Ekin_min ,Ekin_max  ) ; 
00149       h_Theta   [cellNumber] = new TH1D ( TString(Form("h_Theta_%02d",cellNumber)) , TString(Form("h_Theta_%02d",cellNumber)) , Theta_nBins, Theta_min,Theta_max ) ; 
00150       h_Edep    [cellNumber] = new TH1D ( TString(Form("h_Edep_%02d" ,cellNumber)) , TString(Form("h_Edep_%02d" ,cellNumber)) , Edep_nBins , Edep_min ,Edep_max  ) ; 
00151       h_Dose    [cellNumber] = new TH1D ( TString(Form("h_Dose_%02d" ,cellNumber)) , TString(Form("h_Dose_%02d" ,cellNumber)) , Dose_nBins , Dose_min ,Dose_max  ) ; 
00152       h_Range   [cellNumber] = new TH1D ( TString(Form("h_Range_%02d",cellNumber)) , TString(Form("h_Range_%02d",cellNumber)) , Range_nBins, Range_min,Range_max ) ; 
00153       h_Z       [cellNumber] = new TH1D ( TString(Form("h_Z_%02d"    ,cellNumber)) , TString(Form("h_Z_%02d"    ,cellNumber)) , Z_nBins    , Z_min    ,Z_max     ) ;                                    
00154       h_XY      [cellNumber] = new TH2D ( TString(Form("h_XY_%02d"   ,cellNumber)) , TString(Form("h_XY_%02d"   ,cellNumber)) , XY_nBins   , XY_min   ,XY_max    , XY_nBins   , XY_min   ,XY_max    ) ; 
00155    }
00156 
00157 
00158 ////////////////////////////////////////////////////////////
00159 // Filling histograms
00160 ////////////////////////////////////////////////////////////
00161 
00162    Double_t x=-1., y=-1. ; 
00163    Double_t *integral_Edep = new Double_t [nCells], *integral_dose = new Double_t [nCells] ;
00164    for (Int_t cellNumber=0; cellNumber<nCells; cellNumber++) { integral_Edep[cellNumber]=0. ; integral_dose[cellNumber]=0. ; }
00165 
00166    for (Int_t particleNumber=0; particleNumber<nParticles; particleNumber++) 
00167    { 
00168       tree->GetEntry(particleNumber) ; 
00169       if ( TString(myEvent->GetParticleName()) == "alpha" ) 
00170       { 
00171          for (Int_t cellNumber=0; cellNumber<nCells; cellNumber++)
00172          {
00173             if ( myEvent->GetEkin(cellNumber) > -999. ) 
00174             { 
00175                h_Ekin  [cellNumber]->Fill(myEvent->GetEkin    (cellNumber)) ;
00176                h_Theta [cellNumber]->Fill(myEvent->GetTheta   (cellNumber)) ; 
00177                h_Edep  [cellNumber]->Fill(myEvent->GetEdep    (cellNumber)) ; integral_Edep[cellNumber] += myEvent->GetEdep (cellNumber) ; 
00178                h_Dose  [cellNumber]->Fill(myEvent->GetDose    (cellNumber)) ; integral_dose[cellNumber] += myEvent->GetDose (cellNumber) ; 
00179                h_Range [cellNumber]->Fill(myEvent->GetRangeTot(cellNumber)) ; 
00180                h_Z     [cellNumber]->Fill(myEvent->GetRz      (cellNumber)) ; 
00181                x = myEvent->GetRx(cellNumber) ; y = myEvent->GetRy(cellNumber) ; h_XY[cellNumber]->Fill(x/pix_x+XY_nBins/2,y/pix_y+XY_nBins/2) ;
00182             }
00183          }
00184       }
00185    }     
00186 
00187    fin->Close() ; 
00188 
00189 
00190 ////////////////////////////////////////////////////////////
00191 // Plotting
00192 ////////////////////////////////////////////////////////////
00193 
00194    Int_t nTallies = 3 ; 
00195    TCanvas *canvas_geom = new TCanvas ("canvas_geom","canvas_geom",10,10,750,750) ;
00196    canvas_geom->Divide(nTallies,nCells) ; canvas_geom->SetBorderMode(0) ; canvas_geom->SetFrameBorderMode(0) ; canvas_geom->SetFillColor(kWhite) ; 
00197    for (Int_t cellNumber=0; cellNumber<nCells; cellNumber++)
00198    {
00199       h_Z [cellNumber]->SetStats(kFALSE) ; h_Z [cellNumber]->SetTitle("") ;
00200          h_Z [cellNumber]->SetXTitle("Sensor Depth  [#mum]") ; h_Z [cellNumber]->GetXaxis()->CenterTitle() ;
00201          h_Z [cellNumber]->SetYTitle("Number of particles" ) ; h_Z [cellNumber]->GetYaxis()->CenterTitle() ; h_Z [cellNumber]->GetYaxis()->SetTitleOffset(1.7) ; 
00202       h_XY[cellNumber]->SetStats(kFALSE) ; h_XY[cellNumber]->SetTitle("") ;
00203          h_XY[cellNumber]->SetXTitle(TString(Form("Column number (%.0f#times%.0f #mum^{2} pixel)",pix_x,pix_y))) ; h_XY[cellNumber]->GetXaxis()->CenterTitle() ; 
00204          h_XY[cellNumber]->SetYTitle(TString(Form("Row number (%.0f#times%.0f #mum^{2} pixel)"   ,pix_x,pix_y))) ; h_XY[cellNumber]->GetYaxis()->CenterTitle() ; 
00205       canvas_geom->cd(nTallies*cellNumber+1) ; tex_label [cellNumber]->Draw()        ; gPad->SetGrid(0,0) ; gPad->SetRightMargin(0.10) ; gPad->SetLeftMargin(0.18) ;
00206       canvas_geom->cd(nTallies*cellNumber+2) ; h_Z       [cellNumber]->Draw()        ; gPad->SetGrid(0,0) ; gPad->SetRightMargin(0.10) ; gPad->SetLeftMargin(0.18) ;
00207       canvas_geom->cd(nTallies*cellNumber+3) ; h_XY      [cellNumber]->Draw("col2z") ; gPad->SetGrid(0,0) ; gPad->SetRightMargin(0.10) ; gPad->SetLeftMargin(0.18) ;
00208    }
00209    canvas_geom->SaveAs(TString(Form("results/%s__geom.eps",finName_raw.Data()))) ;
00210 
00211 /////
00212 
00213    TCanvas *canvas_kin = new TCanvas ("canvas_kin","canvas_kin",0,0,gClient->GetDisplayWidth(),gClient->GetDisplayHeight()) ;
00214    canvas_kin->Divide(2,2) ; canvas_kin->SetBorderMode(0) ; canvas_kin->SetFrameBorderMode(0) ; canvas_kin->SetFillColor(kWhite) ;  canvas_kin->SetGrid(0,0) ;    
00215       canvas_kin->cd(1) ; gPad->SetGrid(0,0) ;  
00216          Plot ( h_Ekin , "Incident Kinetic Energy  [MeV]", label, nCells, nAlphas_initial, noDummyCell, kFALSE, integral_dose, Edep_max-Edep_min ) ;
00217       canvas_kin->cd(2) ; gPad->SetGrid(0,0) ;  
00218          Plot ( h_Theta, "Incident Angle  [#circ]"       , label, nCells, nAlphas_initial, noDummyCell, kFALSE, integral_dose, Edep_max-Edep_min ) ;
00219       canvas_kin->cd(3) ; gPad->SetGrid(0,0) ;  
00220          Plot ( h_Edep , "Deposited Energy  [MeV]"       , label, nCells, nAlphas_initial, noDummyCell, kTRUE , integral_dose, Edep_max-Edep_min ) ;
00221       canvas_kin->cd(4) ; gPad->SetGrid(0,0) ;  
00222          Plot ( h_Range, "Range  [#mum]"                 , label, nCells, nAlphas_initial, noDummyCell, kFALSE, integral_dose, Edep_max-Edep_min ) ;
00223    canvas_kin->SaveAs(TString(Form("results/%s__kin.eps",finName_raw.Data()))) ;
00224 
00225 
00226 ////////////////////////////////////////////////////////////
00227 // Cleaning
00228 ////////////////////////////////////////////////////////////
00229 
00230 //    for (Int_t cellNumber=0; cellNumber<nCells; cellNumber++)
00231 //    {
00232 //       delete tex_label [cellNumber] ; tex_label [cellNumber] = NULL ; 
00233 //       delete h_Ekin    [cellNumber] ; h_Ekin    [cellNumber] = NULL ; 
00234 //       delete h_Theta   [cellNumber] ; h_Theta   [cellNumber] = NULL ; 
00235 //       delete h_Edep    [cellNumber] ; h_Edep    [cellNumber] = NULL ; 
00236 //       delete h_Dose    [cellNumber] ; h_Dose    [cellNumber] = NULL ; 
00237 //       delete h_Range   [cellNumber] ; h_Range   [cellNumber] = NULL ; 
00238 //       delete h_Z       [cellNumber] ; h_Z       [cellNumber] = NULL ; 
00239 //       delete h_XY      [cellNumber] ; h_XY      [cellNumber] = NULL ; 
00240 //    }
00241 // 
00242 //    delete [] tex_label ; tex_label = NULL ; 
00243 //    delete [] h_Ekin    ; h_Ekin    = NULL ; 
00244 //    delete [] h_Theta   ; h_Theta   = NULL ; 
00245 //    delete [] h_Edep    ; h_Edep    = NULL ; 
00246 //    delete [] h_Dose    ; h_Dose    = NULL ; 
00247 //    delete [] h_Range   ; h_Range   = NULL ; 
00248 //    delete [] h_Z       ; h_Z       = NULL ; 
00249 //    delete [] h_XY      ; h_XY      = NULL ; 
00250 // 
00251 //    delete integral_Edep ; integral_Edep = NULL ; 
00252 //    delete integral_dose ; integral_dose = NULL ; 
00253 
00254 /////
00255 
00256    return ; 
00257 }