Ignore:
Timestamp:
Nov 5, 2010, 4:08:39 PM (14 years ago)
Author:
garnier
Message:

update ti head

Location:
trunk/examples/advanced/microbeam
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/examples/advanced/microbeam/History

    r1337 r1342  
    11-------------------------------------------------------------------
    2 $Id: History,v 1.28 2010/06/10 09:54:05 vnivanch Exp $
     2$Id: History,v 1.29 2010/10/07 14:03:11 sincerti Exp $
    33-------------------------------------------------------------------
    44
    … …  
    1010                       --------------------
    1111
     1207 October 2010 - S. Incerti - tag microbeam-V09-03-07
     13- Updated for use of AIDA fof histograms.
     14
    121510 June 2010 - V.Ivanchenko - tag microbeam-V09-03-06
    1316- Added option "ionGasModel" to the Physics List
    … …  
    1720
    182109 June 2010 - S. Incerti - tag microbeam-V09-03-04
    19 - Switched to physics_list library.
    20 microbeam.mac adapted.
     22- Switched to physics_list library. microbeam.mac adapted.
    2123
    222406 June 2010 - S. Incerti - tag microbeam-V09-03-03

  • trunk/examples/advanced/microbeam/Microbeam.cc

    r1337 r1342  
    2525//
    2626// -------------------------------------------------------------------
    27 // $Id: Microbeam.cc,v 1.10 2010/06/09 18:13:46 sincerti Exp $
     27// $Id: Microbeam.cc,v 1.11 2010/10/07 14:03:11 sincerti Exp $
    2828// -------------------------------------------------------------------
    2929//  GEANT4 - Microbeam example
    … …  
    5353#include "MicrobeamSteppingAction.hh"
    5454#include "MicrobeamSteppingVerbose.hh"
     55#include "MicrobeamHistoManager.hh"
    5556
    5657int main(int argc,char** argv) {
    … …  
    8283  runManager->SetUserAction(KinAct);
    8384 
    84   MicrobeamRunAction* RunAct = new MicrobeamRunAction(detector);
     85  MicrobeamHistoManager*  histo = new MicrobeamHistoManager();
     86
     87  MicrobeamRunAction* RunAct = new MicrobeamRunAction(detector,histo);
    8588 
    8689  runManager->SetUserAction(RunAct);
    87   runManager->SetUserAction(new MicrobeamEventAction(RunAct));
     90  runManager->SetUserAction(new MicrobeamEventAction(RunAct,histo));
    8891  runManager->SetUserAction(new MicrobeamTrackingAction(RunAct));
    89   runManager->SetUserAction(new MicrobeamSteppingAction(RunAct,detector));
     92  runManager->SetUserAction(new MicrobeamSteppingAction(RunAct,detector,histo));
    9093 
    9194  // initialize G4 kernel
    … …  
    9699 
    97100  // local user files created by the simulation
    98   system ("rm -rf dose.txt");
    99   system ("rm -rf 3DDose.txt");
    100   system ("rm -rf stoppingPower.txt");
    101   system ("rm -rf range.txt");
    102   system ("rm -rf beamPosition.txt");
     101  system ("rm -rf microbeam.root");
    103102       
    104103  if (argc==1)   // define UI session for interactive mode.

  • trunk/examples/advanced/microbeam/README

    r1337 r1342  
    11-------------------------------------------------------------------
    2 $Id: README,v 1.11 2010/06/09 18:13:46 sincerti Exp $
     2$Id: README,v 1.12 2010/10/07 14:03:11 sincerti Exp $
    33-------------------------------------------------------------------
    44
    … …  
    1717* e-mail:incerti@cenbg.in2p3.fr
    1818
    19 Last modified by S. Incerti, 09/06/2010
     19Last modified by S. Incerti, 07/10/2010
    2020
    2121---->0. INTRODUCTION.                                                   
    … …  
    8080- MONTE CARLO MICRODOSIMETRY FOR TARGETED IRRADIATION OF INDIVIDUAL CELLS USING
    8181A MICROBEAM FACILITY
    82 By S. Incerti, T. Pouthier, H. Seznec, Ph. Moretto, O. Boissonnade,
    83 T. M. H. Ha, F. Andersson, Ph. Barberet, C. Habchi and D. T. Nguyen
    84 In preparation (2007)
     82By S. Incerti, H. Seznec, M. Simon, Ph. Barberet, C. Habchi, Ph. Moretto
     83Published in Rad. Prot. Dos. 133, 1 (2009) 2-11
    8584
    8685- MONTE CARLO SIMULATION OF THE CENBG MICROBEAM AND NANOBEAM LINES WITH THE
    … …  
    8887By S. Incerti, Q. Zhang, F. Andersson, Ph. Moretto, G.W. Grime,
    8988M.J. Merchant, D.T. Nguyen, C. Habchi, T. Pouthier and H. Seznec
    90 In press in Nucl.Instrum.Meth.B, 2007
     89Published in Nucl. Instrum. and Meth. B 260 (2007) 20-27
    9190
    9291- A COMPARISON OF CELLULAR IRRADIATION TECHNIQUES WITH ALPHA PARTICLES USING
    … …  
    9493By S. Incerti, N. Gault, C. Habchi, J.L.. Lefaix, Ph. Moretto, J.L.. Poncy,
    9594T. Pouthier, H. Seznec. Dec 2006. 3pp.
    96 Published in Rad.Prot.Dos.,1-3,2006 (Micros 2005 special issue).
     95Published in Rad. Prot. Dos. 122, 1-4, (2006) 327-329
    9796
    9897- GEANT4 SIMULATION OF THE NEW CENBG MICRO AND NANO PROBES FACILITY
    … …  
    124123
    125124
    126 ---->3. SET-UP
    127                                                                        
    128 - a standard Geant4 example GNUmakefile is provided                     
     125------->3 VISUALIZATION
    129126
    130 The following section gives the necessary environment variables.                     
    131 
    132 ------->>3.1  ENVIRONMENT VARIABLES
    133 
    134 All variables are defined with their default value.
    135 
    136  - G4SYSTEM = Linux-g++
    137 
    138  - G4INSTALL              points to the installation directory of GEANT4;
    139 
    140  - G4LIB                  point to the compiled libraries of GEANT4;
    141 
    142  - G4WORKDIR              points to the work directory;
    143 
    144  - CLHEP_BASE_DIR         points to the installation directory of CHLEP;
    145 
    146  - G4LEDATA               points to the low energy electromagnetic libraries;
    147 
    148  - LD_LIBRARY_PATH = $CLHEP_BASE_DIR/lib
    149 
    150  - G4LEVELGAMMADATA       points to the photoevaporation library;
    151 
    152  - NeutronHPCrossSections points to the neutron data files;
    153 
    154  - G4RADIOACTIVEDATA      points to the libraries for radio-active decay
    155                           hadronic processes;
    156  
    157 However, the $G4LEVELGAMMADATA, $NeutronHPCrossSections and $G4RADIOACTIVEDATA
    158 variables do not need to be defined for this example.
    159 
    160 Once these variables have been set, simply type gmake to compile the Microbeam
    161 example.
    162 
    163 ------->>3.2  VISUALIZATION
    164 
    165 The user can visualize the targeted cell with OpenGL, DAWN and vrml,
    166 as chosen in the microbeam.mac file. OpenGL is the default viewer. The
    167 cytoplasm in shown in red and the nucleus in green.
    168 
     127The user can visualize the targeted cell by uncommenting the following line in
     128microbeam.mac:
     129#/control/execute vis.mac
    169130
    170131---->4. HOW TO RUN THE EXAMPLE                                         
    171132
    172 In interactive mode, run:
     133The variable G4ANALYSIS_USE must be set to 1.
    173134
    174 > $G4WORDIR/bin/Linux-g++/Microbeam
     135In order to generate histograms, at least one of the AIDA implementations should be
     136 available.
     137 
     138The code should be compiled with gmake and run with :
    175139
    176 The macro microbeam.mac is executed by default.
    177 All visualisation commands are in the vis.mac macro file.
     140> $G4WORDIR/bin/$G4SYSTEM/Microbeam
     141
     142The macro file microbeam.mac is read by default.
    178143
    179144
    … …  
    185150---->6. SIMULATION OUTPUT AND RESULT ANALYZIS                                   
    186151
    187 This example does not need any external analysis package.
    188 The output results consists in several .txt files:
     152The output results consist in several a microbeam.root file, containing several
     153ntuples:
    189154
    190 * dose.txt : gives the total deposited dose in the cell nucleus and in the cell
     155* total deposited dose in the cell nucleus and in the cell
    191156cytoplasm by each incident alpha particle;
    192157
    193 * 3DDose.txt : gives the average on the whole run of the dose deposited per
     158* average on the whole run of the dose deposited per
    194159Voxel per incident alpha particle;
    195160
    196 * range.txt : indicates the final stopping (x,y,z) position of the incident
     161* final stopping (x,y,z) position of the incident
    197162alpha particle within the irradiated medium (cell or culture medium)
    198163
    199 * stoppingPower.txt : gives the actual stopping power dE/dx of the incident
     164* stopping power dE/dx of the incident
    200165alpha particle just before penetrating into the targeted cell;
    201166
    202 * beamPosition.txt : gives the beam transverse position distribution(X and Y)
     167* beam transverse position distribution(X and Y)
    203168just before penetrating into the targeted cell;
    204169
    205 These files can be easily analyzed using for example the provided ROOT macro
     170These results can be easily analyzed using for example the provided ROOT macro
    206171file plot.C; to do so :
    207172* be sure to have ROOT installed on your machine
    … …  
    210175* under your ROOT session, type in : .X plot.C to execute the macro file
    211176
    212 A graphical output obtained with this macro for 40000 incident alpha particles
    213 is shown in the file microbeam.gif
    214 
    215 The simulation predicts that 95% of the incident alpha particles detected by the
    216 gas detector are located within a circle of 10 um in diameter on the target, in
    217 nice agreement with experimental measurements performed on the CENBG setup.
    218177
    219178---------------------------------------------------------------------------

  • trunk/examples/advanced/microbeam/include/MicrobeamEventAction.hh

    r807 r1342  
    2525//
    2626// -------------------------------------------------------------------
    27 // $Id: MicrobeamEventAction.hh,v 1.5 2006/06/29 16:05:05 gunter Exp $
     27// $Id: MicrobeamEventAction.hh,v 1.6 2010/10/07 14:03:11 sincerti Exp $
    2828// -------------------------------------------------------------------
    2929
    … …  
    3535
    3636class MicrobeamRunAction;
     37class MicrobeamHistoManager;
    3738class G4Event;
    3839
    … …  
    4344  public:
    4445 
    45     MicrobeamEventAction(MicrobeamRunAction*);
     46    MicrobeamEventAction(MicrobeamRunAction*, MicrobeamHistoManager *);
    4647   ~MicrobeamEventAction();
    4748
    … …  
    5556 
    5657    MicrobeamRunAction*      Run;
     58    MicrobeamHistoManager*   Histo;     
    5759    G4String                 drawFlag;
    5860    G4int                    printModulo;         

  • trunk/examples/advanced/microbeam/include/MicrobeamRunAction.hh

    r807 r1342  
    2525//
    2626// -------------------------------------------------------------------
    27 // $Id: MicrobeamRunAction.hh,v 1.5 2006/06/29 16:05:13 gunter Exp $
     27// $Id: MicrobeamRunAction.hh,v 1.6 2010/10/07 14:03:11 sincerti Exp $
    2828// -------------------------------------------------------------------
    2929
    … …  
    3535
    3636#include "MicrobeamDetectorConstruction.hh"
     37#include "MicrobeamHistoManager.hh"
    3738
    3839//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    4445public:
    4546 
    46   MicrobeamRunAction(MicrobeamDetectorConstruction*);
     47  MicrobeamRunAction(MicrobeamDetectorConstruction*, MicrobeamHistoManager*);
    4748  ~MicrobeamRunAction();
    4849
    … …  
    8182
    8283  MicrobeamDetectorConstruction* Detector;   
     84  MicrobeamHistoManager* Histo;
    8385  MicrobeamPhantomConfiguration myMicrobeamPhantomConfiguration; 
    8486

  • trunk/examples/advanced/microbeam/include/MicrobeamSteppingAction.hh

    r807 r1342  
    2525//
    2626// -------------------------------------------------------------------
    27 // $Id: MicrobeamSteppingAction.hh,v 1.5 2006/06/29 16:05:15 gunter Exp $
     27// $Id: MicrobeamSteppingAction.hh,v 1.6 2010/10/07 14:03:11 sincerti Exp $
    2828// -------------------------------------------------------------------
    2929
    … …  
    3535#include "MicrobeamRunAction.hh"
    3636#include "MicrobeamDetectorConstruction.hh"
     37#include "MicrobeamHistoManager.hh"
    3738
    3839//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    4142{
    4243public:
    43   MicrobeamSteppingAction(MicrobeamRunAction* ,MicrobeamDetectorConstruction*);
     44  MicrobeamSteppingAction(MicrobeamRunAction* ,MicrobeamDetectorConstruction*,
     45  MicrobeamHistoManager *);
    4446  ~MicrobeamSteppingAction();
    4547 
    … …  
    4951  MicrobeamRunAction*            Run;
    5052  MicrobeamDetectorConstruction* Detector;
     53  MicrobeamHistoManager*         Histo;
    5154  G4float massPhantom;
    5255

  • trunk/examples/advanced/microbeam/microbeam.out

    r1337 r1342  
    55
    66*************************************************************
    7  Geant4 version Name: geant4-09-03-ref-06    (25-June-2010)
     7 Geant4 version Name: geant4-09-03-ref-08    (25-June-2010)
    88                      Copyright : Geant4 Collaboration
    99                      Reference : NIM A 506 (2003), 250-303
    … …  
    137137### G4GoudsmitSaundersonTable loading PDF data
    138138### G4GoudsmitSaundersonMscModel loading ELSEPA data
    139 G4NeutronInelasticXS::G4NeutronInelasticXS: Initialise
    140 G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise
    141139PhysicsList::SetCuts:CutLength : 10 nm
    142140G4ProcessTable::Insert : arguments are 0 pointer
    … …  
    265263          BetheBloch :     Emin=          2 MeV        Emax=   10 TeV
    266264
     265msc:   for kaon+    SubType= 10
     266      Lambda tables from 100 eV  to 10 TeV in 220 bins, spline: 1
     267      RangeFactor= 0.2, stepLimitType: 0, latDisplacement: 1
     268      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
     269          UrbanMsc90 :     Emin=          0 eV         Emax=   10 TeV
     270
    267271hIoni:   for  kaon+    SubType= 2
    268272      dE/dx and range tables from 100 eV  to 10 TeV in 220 bins
    … …  
    364368      ===== EM models for the G4Region  DefaultRegionForTheWorld ======
    365369  eCoulombScattering :     Emin=          0 eV         Emax=   10 TeV
    366 G4NeutronInelasticXS::BuildPhysicsTable:
    367 neutron
    368 G4NEUTRONXSDATA environment variable not set
    369 G4NeutronCaptureXS::BuildPhysicsTable:
    370 neutron
    371 G4NEUTRONXSDATA environment variable not set
    372370
    373371hIoni:   for  pi+    SubType= 2
    … …  
    563561
    564562-> Event # 3 generated
    565    ===> Sorry, the incident alpha particle has missed the targeted cell !
     563   ===> The incident alpha particle has reached the targeted cell :
     564   -----> total absorbed dose within Nucleus   is (Gy) = 0.6101508140564
     565   -----> total absorbed dose within Cytoplasm is (Gy) = 0.035633563995361
    566566
    567567-> Event # 4 generated
    568    ===> Sorry, the incident alpha particle has missed the targeted cell !
     568   ===> The incident alpha particle has reached the targeted cell :
     569   -----> total absorbed dose within Nucleus   is (Gy) = 0.32617983222008
     570   -----> total absorbed dose within Cytoplasm is (Gy) = 0.037377636879683
    569571
    570572-> Event # 5 generated
    571    ===> Sorry, the incident alpha particle has missed the targeted cell !
     573   ===> The incident alpha particle has reached the targeted cell :
     574   -----> total absorbed dose within Nucleus   is (Gy) = 0.35554075241089
     575   -----> total absorbed dose within Cytoplasm is (Gy) = 0.024478180333972
    572576
    573577-> Event # 6 generated
    574    ===> The incident alpha particle has reached the targeted cell :
    575    -----> total absorbed dose within Nucleus   is (Gy) = 0.32255533337593
    576    -----> total absorbed dose within Cytoplasm is (Gy) = 0.022024262696505
     578   ===> Sorry, the incident alpha particle has missed the targeted cell !
    577579
    578580-> Event # 7 generated
    … …  
    589591
    590592ERROR: G4VisCommandsViewerUpdate::SetNewValue: no current viewer.
    591 -> Total number of particles detected by the gas detector : 1
     593-> Total number of particles detected by the gas detector : 3
    592594
    593595Idle> Graphics systems deleted.

  • trunk/examples/advanced/microbeam/plot.C

    r1230 r1342  
    11// -------------------------------------------------------------------
    2 // $Id: plot.C,v 1.5 2009/04/30 10:23:57 sincerti Exp $
     2// $Id: plot.C,v 1.6 2010/10/07 14:03:11 sincerti Exp $
    33// -------------------------------------------------------------------
    44//
    … …  
    2424gROOT->SetStyle("Plain");
    2525Double_t scale;
    26        
     26
     27
    2728c1 = new TCanvas ("c1","",20,20,1200,900);
    2829c1.Divide(4,3);
    29 
    30 FILE * fp = fopen("dose.txt","r");
    31 Float_t nD,cD;
    32 Int_t ncols=0;
    33 Int_t nlines = 0;
    34 
    35 TH1F *h1  = new TH1F("Absorbed dose distribution in Nucleus","Dose distribution in Nucleus",100,0.001,0.5);
    36 TH1F *h10 = new TH1F("Absorbed dose distribution in Cytoplasm","Dose distribution in Cytoplasm",100,0.001,0.2);
    37 
    38 while (1)
    39    {
    40       ncols = fscanf(fp,"%f %f",&nD,&cD);
    41       if (ncols < 0) break;
    42       h1->Fill(nD);
    43       h10->Fill(cD);
    44       nlines++;
    45    }
    46 fclose(fp);
    47 
    48 c1.cd(2);
    49         scale = 1/h1->Integral();
    50         h1->Scale(scale);
    51         h1->Draw();
    52         h1->GetXaxis()->SetLabelSize(0.025);
    53         h1->GetYaxis()->SetLabelSize(0.025);
    54         h1->GetXaxis()->SetTitleSize(0.035);
    55         h1->GetYaxis()->SetTitleSize(0.035);
    56         h1->GetXaxis()->SetTitleOffset(1.4);
    57         h1->GetYaxis()->SetTitleOffset(1.4);
    58         h1->GetXaxis()->SetTitle("Absorbed dose (Gy)");
    59         h1->GetYaxis()->SetTitle("Fraction of events");
    60         h1->SetLineColor(3);
    61         h1->SetFillColor(3);
    62 
    63 //*****************
    64 // DOSE IN CYTOPLASM
    65 //*****************
    66 
    67 c1.cd(6);
    68         scale = 1/h10->Integral();
    69         h10->Scale(scale);
    70         h10->Draw();
    71         h10->GetXaxis()->SetLabelSize(0.025);
    72         h10->GetYaxis()->SetLabelSize(0.025);
    73         h10->GetXaxis()->SetTitleSize(0.035);
    74         h10->GetYaxis()->SetTitleSize(0.035);
    75         h10->GetXaxis()->SetTitleOffset(1.4);
    76         h10->GetYaxis()->SetTitleOffset(1.4);
    77         h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
    78         h10->GetYaxis()->SetTitle("Fraction of events");
    79         h10->SetLineColor(2);
    80         h10->SetFillColor(2);
    81 
    82 //********************************
    83 // STOPPING POWER AT CELL ENTRANCE
    84 //********************************
    85  
    86 gStyle->SetOptStat(0000);
    87 gStyle->SetOptFit();
    88 gStyle->SetPalette(1);
    89 gROOT->SetStyle("Plain");
    90 
    91 Float_t d;
    92 FILE * fp = fopen("stoppingPower.txt","r");
    93 
    94 TH1F *h2 = new TH1F("Beam stopping Power at cell entrance","h1",200,0,300);
    95 while (1)
    96    {
    97       ncols = fscanf(fp,"%f",&d);
    98       if (ncols < 0) break;
    99       h2->Fill(d);
    100       nlines++;
    101    }
    102 fclose(fp);
    103    
    104 c1.cd(9);
    105         scale = 1/h2->Integral();
    106         h2->Scale(scale);
    107         h2->Draw();
    108         h2->GetXaxis()->SetLabelSize(0.025);
    109         h2->GetYaxis()->SetLabelSize(0.025);
    110         h2->GetXaxis()->SetTitleSize(0.035);
    111         h2->GetYaxis()->SetTitleSize(0.035);
    112         h2->GetXaxis()->SetTitleOffset(1.4);
    113         h2->GetYaxis()->SetTitleOffset(1.4);
    114         h2->GetXaxis()->SetTitle("dE/dx (keV/µm)");
    115         h2->GetYaxis()->SetTitle("Fraction of events");
    116         h2->SetTitle("dE/dx at cell entrance");
    117         h2->SetFillColor(4);
    118         h2->SetLineColor(4);
    119         h2->Fit("gaus");
    120         gaus->SetLineColor(6);
    121         h2->Fit("gaus");
    122 
    123 //**************
    124 // RANGE IN CELL
    125 //**************
    126 
    127 Float_t Xc,Zc,X1,Y1,Z1,x,z,X2,Y2,Z2;
    128 Float_t d;
    129 
    130 // X position of target in World
    131 Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
    132 
    133 // Z position of target in World
    134 Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
    135 
    136 // Line alignment (cf MicrobeamEMField.cc)
    137 Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
    138 Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
    139 
    140 FILE * fp = fopen("range.txt","r");
    141 
    142 TNtuple *ntuple = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
    143 
    144 while (1)
    145    {
    146       ncols = fscanf(fp,"%f %f %f",&X1,&Y1,&Z1);
    147       if (ncols < 0) break;
    148       x = X1-Xc;
    149       z = Z1-Zc;
    150       Z2 = z*cos(10*TMath::Pi()/180)-x*sin(10*TMath::Pi()/180);
    151       X2 = z*sin(10*TMath::Pi()/180)+x*cos(10*TMath::Pi()/180);
    152       Y2 = Y1;
    153      
    154       ntuple->Fill(Z2,Y2,X2);
    155       nlines++;
    156    }
    157 fclose(fp);
    158      
    159 c1.cd(10);
    160   ntuple->Draw("X2:Z2","abs(X2)<50","surf3");
    161   gPad->SetLogz();
    162   htemp->GetXaxis()->SetLabelSize(0.025);
    163   htemp->GetYaxis()->SetLabelSize(0.025);
    164   htemp->GetZaxis()->SetLabelSize(0.025);
    165   htemp->GetXaxis()->SetTitleSize(0.035);
    166   htemp->GetYaxis()->SetTitleSize(0.035);
    167   htemp->GetXaxis()->SetTitleOffset(1.4);
    168   htemp->GetYaxis()->SetTitleOffset(1.4);
    169   htemp->GetXaxis()->SetTitle("Z (µm)");
    170   htemp->GetYaxis()->SetTitle("X (µm)");
    171   htemp->SetTitle("Range in cell");
    17230
    17331//*********************
    17432// INTENSITY HISTOGRAMS
    17533//*********************
     34jump:
    17635
    17736FILE * fp = fopen("phantom.dat","r");
    … …  
    18948
    19049TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox");
    191 TNtuple *ntupleZX = new TNtuple("CYTOPASM","ntuple","Z:X:vox");
    192 TNtuple *ntupleYX = new TNtuple("CYTOPASM","ntuple","Y:X:vox");
     50TNtuple *ntupleZX = new TNtuple("CYTOPLASM","ntuple","Z:X:vox");
     51TNtuple *ntupleYX = new TNtuple("CYTOPLASM","ntuple","Y:X:vox");
    19352
    19453nlines=0;
    … …  
    308167  hist->SetTitle("Nucleus intensity on transverse section");
    309168
    310 //****************
    311 // ENERGY DEPOSITS
    312 //****************
    313 
     169//
     170
     171TFile f("microbeam.root");
     172
     173TNtuple* ntuple0;
     174TNtuple* ntuple1;
     175TNtuple* ntuple2;
     176TNtuple* ntuple3;
     177TNtuple* ntuple4;
     178
     179ntuple0 = (TNtuple*)f->Get("ntuple0");
     180ntuple1 = (TNtuple*)f->Get("ntuple1");
     181ntuple2 = (TNtuple*)f->Get("ntuple2");
     182ntuple3 = (TNtuple*)f->Get("ntuple3");
     183ntuple4 = (TNtuple*)f->Get("ntuple4");
     184
     185TH1F *h1  = new TH1F("h1","Dose distribution in Nucleus",100,0.001,1.);
     186TH1F *h10 = new TH1F("h10","Dose distribution in Cytoplasm",100,0.001,.2);
     187
     188c1.cd(2);
     189
     190        ntuple3->Project("h1","doseN");
     191        scale = 1/h1->Integral();
     192        h1->Scale(scale);
     193        h1->Draw();
     194        h1->GetXaxis()->SetLabelSize(0.025);
     195        h1->GetYaxis()->SetLabelSize(0.025);
     196        h1->GetXaxis()->SetTitleSize(0.035);
     197        h1->GetYaxis()->SetTitleSize(0.035);
     198        h1->GetXaxis()->SetTitleOffset(1.4);
     199        h1->GetYaxis()->SetTitleOffset(1.4);
     200        h1->GetXaxis()->SetTitle("Absorbed dose (Gy)");
     201        h1->GetYaxis()->SetTitle("Fraction of events");
     202        h1->SetLineColor(3);
     203        h1->SetFillColor(3);
     204
     205//*****************
     206// DOSE IN CYTOPLASM
     207//*****************
     208
     209c1.cd(6);
     210        ntuple3->Project("h10","doseC");
     211        scale = 1/h10->Integral();
     212        h10->Scale(scale);
     213        h10->Draw();
     214        h10->GetXaxis()->SetLabelSize(0.025);
     215        h10->GetYaxis()->SetLabelSize(0.025);
     216        h10->GetXaxis()->SetTitleSize(0.035);
     217        h10->GetYaxis()->SetTitleSize(0.035);
     218        h10->GetXaxis()->SetTitleOffset(1.4);
     219        h10->GetYaxis()->SetTitleOffset(1.4);
     220        h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
     221        h10->GetYaxis()->SetTitle("Fraction of events");
     222        h10->SetLineColor(2);
     223        h10->SetFillColor(2);
     224
     225//********************************
     226// STOPPING POWER AT CELL ENTRANCE
     227//********************************
     228 
    314229gStyle->SetOptStat(0000);
    315230gStyle->SetOptFit();
    … …  
    317232gROOT->SetStyle("Plain");
    318233
    319 FILE * fp = fopen("3DDose.txt","r");
    320 
    321 TNtuple *ntuple2 = new TNtuple("CELL","ntuple","yVox:xVox:dose");
    322 TNtuple *ntuple3 = new TNtuple("CELL","ntuple","xVox:zVox:dose");
    323 TNtuple *ntuplezyx = new TNtuple("DOSE","ntuple","zVox:yVox:xVox:dose");
    324 
    325 while (1)
    326    {
    327       ncols = fscanf(fp,"%f %f %f %f",&xVox, &yVox, &zVox, &dose);
    328       if (ncols < 0) break;
    329      
    330       xVox= xVox*voxelSizeX;
    331       yVox= yVox*voxelSizeY;
    332       zVox= zVox*voxelSizeZ;
    333      
    334       ntuple2->Fill(yVox,xVox,dose);
    335       ntuple3->Fill(xVox,zVox,dose);
    336       ntuplezyx->Fill(zVox,yVox,xVox,dose);
    337    }
    338 fclose(fp);
     234Float_t d;
     235
     236TH1F *h2 = new TH1F("h2","Beam stopping power at cell entrance",200,0,300);
     237   
     238c1.cd(9);
     239        ntuple0->Project("h2","sp");
     240        scale = 1/h2->Integral();
     241        h2->Scale(scale);
     242        h2->Draw();
     243        h2->GetXaxis()->SetLabelSize(0.025);
     244        h2->GetYaxis()->SetLabelSize(0.025);
     245        h2->GetXaxis()->SetTitleSize(0.035);
     246        h2->GetYaxis()->SetTitleSize(0.035);
     247        h2->GetXaxis()->SetTitleOffset(1.4);
     248        h2->GetYaxis()->SetTitleOffset(1.4);
     249        h2->GetXaxis()->SetTitle("dE/dx (keV/µm)");
     250        h2->GetYaxis()->SetTitle("Fraction of events");
     251        h2->SetTitle("dE/dx at cell entrance");
     252        h2->SetFillColor(4);
     253        h2->SetLineColor(4);
     254        h2->Fit("gaus");
     255        gaus->SetLineColor(6);
     256        h2->Fit("gaus");
     257
     258
     259//**************
     260// RANGE IN CELL
     261//**************
     262
     263Double_t Xc,Zc,X1,Y1,Z1,X2,Y2,Z2;
     264
     265// X position of target in World
     266Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
     267
     268// Z position of target in World
     269Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
     270
     271// Line alignment (cf MicrobeamEMField.cc)
     272Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
     273Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
     274
     275TNtuple *ntupleR = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
     276Double_t x,y,z,xx,zz;
     277ntuple2->SetBranchAddress("x",&x);
     278ntuple2->SetBranchAddress("y",&y);
     279ntuple2->SetBranchAddress("z",&z);
     280Int_t nentries = (Int_t)ntuple2->GetEntries();
     281for (Int_t i=0;i<nentries;i++)
     282{
     283      ntuple2->GetEntry(i);
     284      X1=x;
     285      Y1=y;
     286      Z1=z;
     287      xx = X1-Xc;
     288      zz = Z1-Zc;
     289      Z2 = zz*cos(10*TMath::Pi()/180)-xx*sin(10*TMath::Pi()/180);
     290      X2 = zz*sin(10*TMath::Pi()/180)+xx*cos(10*TMath::Pi()/180);
     291      Y2 = Y1;   
     292      ntupleR->Fill(Z2,Y2,X2);
     293}
     294     
     295c1.cd(10);
     296  ntupleR->Draw("X2:Z2","abs(X2)<50","surf3");
     297  gPad->SetLogz();
     298  /*
     299  htemp->GetXaxis()->SetLabelSize(0.025);
     300  htemp->GetYaxis()->SetLabelSize(0.025);
     301  htemp->GetZaxis()->SetLabelSize(0.025);
     302  htemp->GetXaxis()->SetTitleSize(0.035);
     303  htemp->GetYaxis()->SetTitleSize(0.035);
     304  htemp->GetXaxis()->SetTitleOffset(1.4);
     305  htemp->GetYaxis()->SetTitleOffset(1.4);
     306  htemp->GetXaxis()->SetTitle("Z (µm)");
     307  htemp->GetYaxis()->SetTitle("X (µm)");
     308  htemp->SetTitle("Range in cell");
     309  */
     310
     311//****************
     312// ENERGY DEPOSITS
     313//****************
     314
     315
     316gStyle->SetOptStat(0000);
     317gStyle->SetOptFit();
     318gStyle->SetPalette(1);
     319gROOT->SetStyle("Plain");
    339320
    340321c1.cd(11);
    341322  TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
    342   ntuple2->Draw("yVox:xVox>>hist","dose","contz");
     323  ntuple4->Draw("y*0.359060:x*0.359060>>hist","doseV","contz");
    343324  gPad->SetLogz();
    344325  hist->Draw("contz");
    … …  
    353334  hist->GetYaxis()->SetTitle("X (µm)");
    354335  hist->SetTitle("Mean energy deposit -transverse- (z axis in eV)");
     336
    355337 
    356338c1.cd(12);
    357339  TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
    358   ntuple3->Draw("xVox:zVox>>hist","dose","contz");
     340  ntuple4->Draw("x*0.359060:z*0.162810>>hist","doseV","contz");
    359341  gPad->SetLogz();
    360342  hist->Draw("contz");
    … …  
    379361gROOT->SetStyle("Plain");
    380362
    381 Float_t bx, by;
    382 FILE * fp = fopen("beamPosition.txt","r");
    383 
    384 TH1F *h77 = new TH1F("Beam X transverse position at cell entrance","h1",200,-10,10);
    385 TH1F *h88 = new TH1F("Beam Y transverse position at cell entrance","h1",200,-10,10);
    386 while (1)
    387    {
    388       ncols = fscanf(fp,"%f %f",&bx, &by);
    389       if (ncols < 0) break;
    390       h77->Fill(bx);
    391       h88->Fill(by);
    392       nlines++;
    393     }
    394 fclose(fp);
    395    
     363TH1F *h77 = new TH1F("hx","h1",200,-10,10);
     364TH1F *h88 = new TH1F("hy","h1",200,-10,10);
     365   
    396366c1.cd(4);
    397         scale = 1/h77->Integral();
     367        ntuple1->Project("hx","x");
     368        scale = 1/h77->Integral();
    398369        h77->Scale(scale);
    399370        h77->Draw();
    … …  
    413384
    414385c1.cd(8);
    415         scale = 1/h88->Integral();
     386        ntuple1->Project("hy","y");
     387        scale = 1/h88->Integral();
    416388        h88->Scale(scale);
    417389        h88->Draw();

  • trunk/examples/advanced/microbeam/src/MicrobeamEventAction.cc

    r1337 r1342  
    2525//
    2626// -------------------------------------------------------------------
    27 // $Id: MicrobeamEventAction.cc,v 1.8 2010/06/07 03:18:01 sincerti Exp $
     27// $Id: MicrobeamEventAction.cc,v 1.9 2010/10/07 14:03:11 sincerti Exp $
    2828// -------------------------------------------------------------------
    2929
    … …  
    3333#include "MicrobeamEventAction.hh"
    3434#include "MicrobeamRunAction.hh"
     35#include "MicrobeamHistoManager.hh"
    3536
    3637//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    3738
    38 MicrobeamEventAction::MicrobeamEventAction(MicrobeamRunAction* run)
    39   :Run(run),drawFlag("all"),printModulo(10000)
     39MicrobeamEventAction::MicrobeamEventAction(MicrobeamRunAction* run,
     40MicrobeamHistoManager * his)
     41:Run(run),Histo(his),drawFlag("all"),printModulo(10000)
    4042{}
    4143
    … …  
    6466if (Run->GetDoseN()>0 || Run->GetDoseC()>0)
    6567{
    66         FILE* myFile;
    67         myFile=fopen("dose.txt","a");
    68         fprintf(myFile,"%e %e\n",Run->GetDoseN(),Run->GetDoseC());
    69         fclose (myFile);
     68        Histo->FillNtuple(3,0,Run->GetDoseN());
     69        Histo->FillNtuple(3,1,Run->GetDoseC());
     70        Histo->AddRowNtuple(3);                 
     71
    7072        G4cout << "   ===> The incident alpha particle has reached the targeted cell :" << G4endl;
    7173        G4cout << "   -----> total absorbed dose within Nucleus   is (Gy) = " << Run->GetDoseN() << G4endl;

  • trunk/examples/advanced/microbeam/src/MicrobeamRunAction.cc

    r807 r1342  
    2525//
    2626// -------------------------------------------------------------------
    27 // $Id: MicrobeamRunAction.cc,v 1.6 2006/06/29 16:05:37 gunter Exp $
     27// $Id: MicrobeamRunAction.cc,v 1.7 2010/10/07 14:03:11 sincerti Exp $
    2828// -------------------------------------------------------------------
    2929
    … …  
    3636//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    3737
    38 MicrobeamRunAction::MicrobeamRunAction(MicrobeamDetectorConstruction* det)
    39 :Detector(det)
     38MicrobeamRunAction::MicrobeamRunAction(MicrobeamDetectorConstruction* det,
     39MicrobeamHistoManager * his)
     40:Detector(det),Histo(his)
    4041{   
    4142  saveRndm = 0; 
    … …  
    5556{ 
    5657 
     58  // Histograms
     59  Histo->book();
     60
    5761  // save Rndm status
    5862  if (saveRndm > 0)
    … …  
    99103  }   
    100104 
    101   FILE *myFile;
    102   myFile=fopen("3DDose.txt","a");
    103105  for (G4int i=0; i<nbOfPixels; i++)
    104106  { 
    … …  
    106108    v = mapVoxels[i];
    107109    if ( (GetNumEvent()+1) !=0)
    108       fprintf (myFile,"%f %f %f %f \n", v.x(), v.y(), v.z(), dose3DDose[i]/(GetNumEvent()+1));
     110      {
     111          Histo->FillNtuple(4,0,v.x());
     112          Histo->FillNtuple(4,1,v.y());
     113          Histo->FillNtuple(4,2,v.z());
     114          Histo->FillNtuple(4,3,dose3DDose[i]/(GetNumEvent()+1));
     115          Histo->AddRowNtuple(4);                           
     116      }
    109117  }
    110   fclose (myFile);                             
    111    
     118   
    112119  G4cout << "-> Total number of particles detected by the gas detector : " << GetNbOfHitsGas() << G4endl; 
    113120  G4cout << G4endl;   
     121 
     122  //save histograms     
     123  Histo->save();
     124
    114125}

  • trunk/examples/advanced/microbeam/src/MicrobeamSteppingAction.cc

    r1230 r1342  
    2525//
    2626// -------------------------------------------------------------------
    27 // $Id: MicrobeamSteppingAction.cc,v 1.9 2008/06/16 07:46:11 sincerti Exp $
     27// $Id: MicrobeamSteppingAction.cc,v 1.10 2010/10/07 14:03:11 sincerti Exp $
    2828// -------------------------------------------------------------------
    2929
    … …  
    3636//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    3737
    38 MicrobeamSteppingAction::MicrobeamSteppingAction(MicrobeamRunAction* run,MicrobeamDetectorConstruction* det)
    39 :Run(run),Detector(det)
     38MicrobeamSteppingAction::MicrobeamSteppingAction(MicrobeamRunAction* run,MicrobeamDetectorConstruction* det,
     39MicrobeamHistoManager* his)
     40:Run(run),Detector(det),Histo(his)
    4041{ }
    4142
    … …  
    8687        {
    8788       
    88         if(aStep->GetTotalEnergyDeposit()>0)
     89         if( (aStep->GetPreStepPoint()->GetKineticEnergy() - aStep->GetPostStepPoint()->GetKineticEnergy() ) >0)
    8990         {
    90           FILE *myFile;
    91           myFile=fopen("stoppingPower.txt","a");       
    92           fprintf(myFile,"%e \n",
    93           (aStep->GetTotalEnergyDeposit()/keV)/(aStep->GetStepLength()/micrometer));   
    94           fclose (myFile);
     91          Histo->FillNtuple(0,0,aStep->GetPreStepPoint()->GetKineticEnergy()/keV);
     92          Histo->FillNtuple(0,1,
     93            (aStep->GetPreStepPoint()->GetKineticEnergy() -
     94             aStep->GetPostStepPoint()->GetKineticEnergy())/ keV/(aStep->GetStepLength()/micrometer));
     95          Histo->AddRowNtuple(0);
    9596         }
    9697
    97          // Average dE over step syggested by Michel Maire
     98         // Average dE over step suggested by Michel Maire
    9899
    99100         G4StepPoint* p1 = aStep->GetPreStepPoint();
    … …  
    111112         // end
    112113
    113          FILE *myFile;
    114          myFile=fopen("beamPosition.txt","a");
    115          fprintf
    116          (
    117           myFile,"%e %e \n",
    118           localPosition.x()/micrometer,
    119           localPosition.y()/micrometer
    120          );
    121          fclose (myFile);                               
     114         Histo->FillNtuple(1,0,localPosition.x()/micrometer);
     115         Histo->FillNtuple(1,1,localPosition.y()/micrometer);
     116         Histo->AddRowNtuple(1);                       
    122117
    123118        }
    … …  
    142137       
    143138        {               
    144          FILE *myFile;
    145          myFile=fopen("range.txt","a");
    146          fprintf
    147          ( myFile,"%e %e %e\n",
    148           (aStep->GetTrack()->GetPosition().x())/micrometer,
    149           (aStep->GetTrack()->GetPosition().y())/micrometer,
    150           (aStep->GetTrack()->GetPosition().z())/micrometer
    151          );
    152          fclose (myFile);                               
     139         Histo->FillNtuple(2,0,aStep->GetPostStepPoint()->GetPosition().x()/micrometer);
     140         Histo->FillNtuple(2,1,aStep->GetPostStepPoint()->GetPosition().y()/micrometer);
     141         Histo->FillNtuple(2,2,aStep->GetPostStepPoint()->GetPosition().z()/micrometer);
     142         Histo->AddRowNtuple(2);                       
    153143        }
    154144
    155145// TOTAL DOSE DEPOSIT AND DOSE DEPOSIT WITHIN A PHANTOM VOXEL
     146// FOR ALL PARTICLES
    156147
    157148if (aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()  == "physicalNucleus")
Note: See TracChangeset for help on using the changeset viewer.