Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/underground_physics/src/DMXEventAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/underground_physics/src/DMXEventAction.cc (Version 11.3.0) and /examples/advanced/underground_physics/src/DMXEventAction.cc (Version 9.6.p4)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 // -------------------------------------------     27 // --------------------------------------------------------------
 28 //   GEANT 4 - Underground Dark Matter Detecto     28 //   GEANT 4 - Underground Dark Matter Detector Advanced Example
 29 //                                                 29 //
 30 //      For information related to this code c     30 //      For information related to this code contact: Alex Howard
 31 //      e-mail: alexander.howard@cern.ch           31 //      e-mail: alexander.howard@cern.ch
 32 // -------------------------------------------     32 // --------------------------------------------------------------
 33 // Comments                                        33 // Comments
 34 //                                                 34 //
 35 //                  Underground Advanced           35 //                  Underground Advanced
 36 //               by A. Howard and H. Araujo        36 //               by A. Howard and H. Araujo 
 37 //                    (27th November 2001)         37 //                    (27th November 2001)
 38 //                                                 38 //
 39 // History/Additions:                              39 // History/Additions:
 40 // 16 Jan 2002  Added analysis                     40 // 16 Jan 2002  Added analysis
 41 //                                                 41 //
 42 //                                                 42 //
 43 // EventAction program                             43 // EventAction program
 44 // -------------------------------------------     44 // --------------------------------------------------------------
 45                                                    45 
 46 #include "DMXEventAction.hh"                       46 #include "DMXEventAction.hh"
 47                                                    47 
 48 // pass parameters for messengers:                 48 // pass parameters for messengers:
 49 #include "DMXRunAction.hh"                         49 #include "DMXRunAction.hh"
 50 #include "DMXPrimaryGeneratorAction.hh"            50 #include "DMXPrimaryGeneratorAction.hh"
 51                                                    51 
 52 // note DMXPmtHit.hh and DMXScintHit.hh are in     52 // note DMXPmtHit.hh and DMXScintHit.hh are included in DMXEventAction.hh
 53                                                    53 
 54 #include "DMXEventActionMessenger.hh"              54 #include "DMXEventActionMessenger.hh"
 55                                                    55 
 56 #include "G4AnalysisManager.hh"                <<  56 #ifdef G4ANALYSIS_USE
                                                   >>  57 #include "DMXAnalysisManager.hh"
                                                   >>  58 #endif
                                                   >>  59 
                                                   >>  60 #include "G4SystemOfUnits.hh"
 57 #include "G4Event.hh"                              61 #include "G4Event.hh"
 58 #include "G4EventManager.hh"                       62 #include "G4EventManager.hh"
 59 #include "G4HCofThisEvent.hh"                      63 #include "G4HCofThisEvent.hh"
 60 #include "G4VHitsCollection.hh"                    64 #include "G4VHitsCollection.hh"
 61 #include "G4TrajectoryContainer.hh"                65 #include "G4TrajectoryContainer.hh"
 62 #include "G4Trajectory.hh"                         66 #include "G4Trajectory.hh"
 63 #include "G4VVisManager.hh"                        67 #include "G4VVisManager.hh"
 64 #include "G4SDManager.hh"                          68 #include "G4SDManager.hh"
 65 #include "G4UImanager.hh"                          69 #include "G4UImanager.hh"
 66 #include "G4UnitsTable.hh"                         70 #include "G4UnitsTable.hh"
 67 #include "G4RunManager.hh"                     <<  71 #include "G4ios.hh"
 68 #include "G4Threading.hh"                      << 
 69 #include <fstream>                                 72 #include <fstream>
 70 #include <iomanip>                                 73 #include <iomanip>
 71                                                    74 
 72 #include "G4SystemOfUnits.hh"                  << 
 73                                                    75 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75 DMXEventAction::DMXEventAction()               <<  77 DMXEventAction::DMXEventAction(DMXRunAction* DMXRun,
 76   : runAct(0),genAction(0),hitsfile(0),pmtfile <<  78              DMXPrimaryGeneratorAction* DMXGenerator) 
                                                   >>  79   : runAct(DMXRun),genAction(DMXGenerator)
 77 {                                                  80 {
 78                                                    81 
 79   // create messenger                              82   // create messenger
 80   eventMessenger = new DMXEventActionMessenger     83   eventMessenger = new DMXEventActionMessenger(this);
 81                                                    84 
 82   // defaults for messenger                        85   // defaults for messenger
 83   drawColsFlag = "standard";                       86   drawColsFlag = "standard";
 84   drawTrksFlag = "all";                            87   drawTrksFlag = "all";
 85   drawHitsFlag = 1;                                88   drawHitsFlag = 1;
 86   savePmtFlag  = 0;                                89   savePmtFlag  = 0;
 87   saveHitsFlag = 1;                                90   saveHitsFlag = 1;
 88                                                    91   
 89   printModulo = 1;                                 92   printModulo = 1;
 90                                                    93 
 91   // hits collections                              94   // hits collections
 92   scintillatorCollID = -1;                         95   scintillatorCollID = -1;
 93   pmtCollID = -1;                                  96   pmtCollID = -1;
 94                                                    97 
 95   energy_pri=0;                                    98   energy_pri=0;
 96   seeds=NULL;                                      99   seeds=NULL;
 97                                                   100 
 98 }                                                 101 }
 99                                                   102 
100                                                   103 
101 //....oooOO0OOooo........oooOO0OOooo........oo    104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 DMXEventAction::~DMXEventAction() {               105 DMXEventAction::~DMXEventAction() {
103                                                   106 
104   if (hitsfile)                                << 
105     {                                          << 
106       hitsfile->close();                       << 
107       delete hitsfile;                         << 
108     }                                          << 
109   if (pmtfile)                                 << 
110     {                                          << 
111       pmtfile->close();                        << 
112       delete pmtfile;                          << 
113     }                                          << 
114   delete eventMessenger;                          107   delete eventMessenger;
115 }                                                 108 }
116                                                   109 
117                                                   110 
118 //....oooOO0OOooo........oooOO0OOooo........oo    111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 void DMXEventAction::BeginOfEventAction(const  << 112 void DMXEventAction::BeginOfEventAction(const G4Event* evt) {
120 {                                              << 
121                                                << 
122   //thread-local run action                    << 
123   if (!runAct)                                 << 
124     runAct =                                   << 
125       dynamic_cast<const DMXRunAction*>        << 
126       (G4RunManager::GetRunManager()->GetUserR << 
127                                                << 
128   if (!genAction)                              << 
129     genAction = dynamic_cast<const DMXPrimaryG << 
130       (G4RunManager::GetRunManager()->GetUserP << 
131                                                   113 
132                                                   114 
133   // grab seeds                                   115   // grab seeds
134   seeds = genAction->GetEventSeeds();             116   seeds = genAction->GetEventSeeds();
135                                                   117 
136   // grab energy of primary                       118   // grab energy of primary
137   energy_pri = genAction->GetEnergyPrimary();     119   energy_pri = genAction->GetEnergyPrimary();
138                                                   120 
139   event_id = evt->GetEventID();                   121   event_id = evt->GetEventID();
140                                                   122  
141   // print this information event by event (mo    123   // print this information event by event (modulo n)   
142   if (event_id%printModulo == 0)                  124   if (event_id%printModulo == 0)
143     {                                             125     {
144       G4cout << "\n---> Begin of event: " << e    126       G4cout << "\n---> Begin of event: " << event_id << G4endl;
145       G4cout << "\n     Primary Energy: " << G    127       G4cout << "\n     Primary Energy: " << G4BestUnit(energy_pri,"Energy") 
146        << G4endl;                                 128        << G4endl;
147       //      HepRandom::showEngineStatus();      129       //      HepRandom::showEngineStatus(); 
148     }                                             130     }
149                                                   131 
150                                                   132 
151   // get ID for scintillator hits collection      133   // get ID for scintillator hits collection
152   if (scintillatorCollID==-1) {                   134   if (scintillatorCollID==-1) {
153     G4SDManager *SDman = G4SDManager::GetSDMpo    135     G4SDManager *SDman = G4SDManager::GetSDMpointer();
154     scintillatorCollID = SDman->GetCollectionI    136     scintillatorCollID = SDman->GetCollectionID("scintillatorCollection");
155   }                                               137   }
156                                                   138 
157   // get ID for pmt hits collection               139   // get ID for pmt hits collection
158   if (pmtCollID==-1) {                            140   if (pmtCollID==-1) {
159     G4SDManager *SDman = G4SDManager::GetSDMpo    141     G4SDManager *SDman = G4SDManager::GetSDMpointer();
160     pmtCollID = SDman->GetCollectionID("pmtCol    142     pmtCollID = SDman->GetCollectionID("pmtCollection");
161   }                                               143   }
162                                                   144 
163 }                                                 145 }
164                                                   146 
165                                                   147 
166 //....oooOO0OOooo........oooOO0OOooo........oo    148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167 void DMXEventAction::EndOfEventAction(const G4    149 void DMXEventAction::EndOfEventAction(const G4Event* evt) {
168                                                   150 
169   // check that both hits collections have bee    151   // check that both hits collections have been defined
170   if(scintillatorCollID<0||pmtCollID<0) return    152   if(scintillatorCollID<0||pmtCollID<0) return;
171                                                   153 
172   G4AnalysisManager* man = G4AnalysisManager:: << 
173                                                << 
174   // address hits collections                     154   // address hits collections
175   DMXScintHitsCollection* SHC = NULL;             155   DMXScintHitsCollection* SHC = NULL;
176   DMXPmtHitsCollection*   PHC = NULL;             156   DMXPmtHitsCollection*   PHC = NULL;
177   G4HCofThisEvent* HCE = evt->GetHCofThisEvent    157   G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
178   if(HCE) {                                       158   if(HCE) {
179     SHC = (DMXScintHitsCollection*)(HCE->GetHC    159     SHC = (DMXScintHitsCollection*)(HCE->GetHC(scintillatorCollID));
180     PHC = (DMXPmtHitsCollection*)(HCE->GetHC(p    160     PHC = (DMXPmtHitsCollection*)(HCE->GetHC(pmtCollID));
181   }                                               161   }
182                                                   162 
183   // event summary                                163   // event summary
184   totEnergy         = 0.;                         164   totEnergy         = 0.;
185   totEnergyGammas   = 0.;                         165   totEnergyGammas   = 0.;
186   totEnergyNeutrons = 0.;                         166   totEnergyNeutrons = 0.;
187   firstParticleE    = 0.;                         167   firstParticleE    = 0.;
188   particleEnergy    = 0.;                         168   particleEnergy    = 0.;
189   firstLXeHitTime   = 0.;                         169   firstLXeHitTime   = 0.;
190   aveTimePmtHits    = 0.;                         170   aveTimePmtHits    = 0.;
191                                                   171 
192   firstParticleName = "";                         172   firstParticleName = "";
193   particleName      = "";                         173   particleName      = "";
194                                                   174 
195                                                   175 
196   // particle flags                               176   // particle flags
197   gamma_ev          = false;                      177   gamma_ev          = false;
198   neutron_ev        = false;                      178   neutron_ev        = false;
199   positron_ev       = false;                      179   positron_ev       = false;
200   electron_ev       = false;                      180   electron_ev       = false;
201   proton_ev         = false;                      181   proton_ev         = false;
202   other_ev          = false;                      182   other_ev          = false;
203   start_gamma       = false;                      183   start_gamma       = false;
204   start_neutron     = false;                      184   start_neutron     = false;
205                                                   185 
206                                                   186 
207   // scintillator hits                            187   // scintillator hits
208   if(SHC) {                                       188   if(SHC) {
209     S_hits = SHC->entries();                      189     S_hits = SHC->entries();
210                                                   190     
211     for (G4int i=0; i<S_hits; i++) {              191     for (G4int i=0; i<S_hits; i++) {
212       if(i==0) {                                  192       if(i==0) {
213   firstParticleName = (*SHC)[0]->GetParticle()    193   firstParticleName = (*SHC)[0]->GetParticle();
214   firstLXeHitTime   = (*SHC)[0]->GetTime();       194   firstLXeHitTime   = (*SHC)[0]->GetTime();
215   firstParticleE = (*SHC)[0]->GetParticleEnerg    195   firstParticleE = (*SHC)[0]->GetParticleEnergy();
216   if (event_id%printModulo == 0 && S_hits > 0)    196   if (event_id%printModulo == 0 && S_hits > 0) {
217     G4cout << "     First hit in LXe: " << fir    197     G4cout << "     First hit in LXe: " << firstParticleName << G4endl;
218     G4cout << "     Number of hits in LXe: " <    198     G4cout << "     Number of hits in LXe: " << S_hits << G4endl; 
219   }                                               199   }
220       }                                           200       }
221       hitEnergy         = (*SHC)[i]->GetEdep()    201       hitEnergy         = (*SHC)[i]->GetEdep();
222       totEnergy        += hitEnergy;              202       totEnergy        += hitEnergy;
223                                                   203       
224       particleName      = (*SHC)[i]->GetPartic    204       particleName      = (*SHC)[i]->GetParticle();
225       particleEnergy    = (*SHC)[i]->GetPartic    205       particleEnergy    = (*SHC)[i]->GetParticleEnergy();
226                                                   206 
227       if(particleName == "gamma") {               207       if(particleName == "gamma") {
228   gamma_ev = true;                                208   gamma_ev = true;
229   start_gamma = true;                             209   start_gamma = true;
230   start_neutron = false;                          210   start_neutron = false;
231       }                                           211       }
232       else if(particleName == "neutron")          212       else if(particleName == "neutron") 
233   neutron_ev = true;                              213   neutron_ev = true;
234       else if(particleName == "e+")               214       else if(particleName == "e+") 
235   positron_ev = true;                             215   positron_ev = true;
236       else if(particleName == "e-")               216       else if(particleName == "e-") 
237   electron_ev = true;                             217   electron_ev = true;
238       else if(particleName == "proton")           218       else if(particleName == "proton") 
239   proton_ev = true;                               219   proton_ev = true;
240       else {                                      220       else {
241   other_ev = true;                                221   other_ev = true;
242   start_gamma = false;                            222   start_gamma = false;
243   start_neutron = true;                           223   start_neutron = true;
244       }                                           224       }
245                                                   225 
246       if(start_gamma && !start_neutron)           226       if(start_gamma && !start_neutron) 
247   totEnergyGammas += hitEnergy;                   227   totEnergyGammas += hitEnergy;
248       if(start_neutron && !start_gamma)           228       if(start_neutron && !start_gamma) 
249   totEnergyNeutrons += hitEnergy;                 229   totEnergyNeutrons += hitEnergy;
250     }                                             230     }
251                                                   231     
252     if (event_id%printModulo == 0)                232     if (event_id%printModulo == 0)
253       G4cout << "     Total energy in LXe: "      233       G4cout << "     Total energy in LXe: " 
254        << G4BestUnit(totEnergy,"Energy") << G4    234        << G4BestUnit(totEnergy,"Energy") << G4endl;  
255                                                   235     
256   }                                               236   }
257                                                   237   
258                                                   238   
259   // PMT hits                                     239   // PMT hits
260   if(PHC) {                                       240   if(PHC) {
261     P_hits = PHC->entries();                      241     P_hits = PHC->entries();
262                                                   242     
263     // average time of PMT hits                   243     // average time of PMT hits
264     for (G4int i=0; i<P_hits; i++) {              244     for (G4int i=0; i<P_hits; i++) {
265       G4double time = ( (*PHC)[i]->GetTime() -    245       G4double time = ( (*PHC)[i]->GetTime() - firstLXeHitTime );
266       aveTimePmtHits += time / (G4double)P_hit    246       aveTimePmtHits += time / (G4double)P_hits;
267       ////      if (event_id == 0) {              247       ////      if (event_id == 0) {
268       if(P_hits >= 0) {                           248       if(P_hits >= 0) {
269   man->FillH1(7,time);                         << 249 #ifdef G4ANALYSIS_USE
270       }                                        << 250   // pass first event for histogram record:
                                                   >> 251   DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
                                                   >> 252   analysis->HistTime(time);
                                                   >> 253 #endif
271     }                                             254     }
                                                   >> 255   }
272                                                   256   
273     if (event_id%printModulo == 0 && P_hits >     257     if (event_id%printModulo == 0 && P_hits > 0) {
274       G4cout << "     Average light collection    258       G4cout << "     Average light collection time: "
275        << G4BestUnit(aveTimePmtHits,"Time") <<    259        << G4BestUnit(aveTimePmtHits,"Time") << G4endl;
276       G4cout << "     Number of PMT hits (phot    260       G4cout << "     Number of PMT hits (photoelectron equivalent): " 
277        << P_hits << G4endl;                    << 261        << P_hits << G4endl;
                                                   >> 262       //#ifdef G4ANALYSIS_USE
                                                   >> 263       // plot histograms, interactively:
                                                   >> 264       //G4bool plotevent=runAct->Getplotevent();      
                                                   >> 265       //if(plotevent) {
                                                   >> 266       //  DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
                                                   >> 267       //  analysis->PlotHistosInter(P_hits);
                                                   >> 268       //}
                                                   >> 269       //#endif
278     }                                             270     }
279     // write out (x,y,z) of PMT hits              271     // write out (x,y,z) of PMT hits
280     if (savePmtFlag)                              272     if (savePmtFlag)
281       writePmtHitsToFile(PHC);                    273       writePmtHitsToFile(PHC);
282   }                                               274   }
283                                                   275   
284                                                   276 
285   // write out event summary                      277   // write out event summary
286   if(saveHitsFlag)                                278   if(saveHitsFlag) 
287     writeScintHitsToFile();                       279     writeScintHitsToFile();
288                                                   280   
289   // draw trajectories                            281   // draw trajectories
290   if(drawColsFlag=="standard" && drawTrksFlag!    282   if(drawColsFlag=="standard" && drawTrksFlag!="none")
291     drawTracks(evt);                              283     drawTracks(evt);
292                                                   284 
293   // hits in PMT                                  285   // hits in PMT
294   if(drawHitsFlag && PHC)                         286   if(drawHitsFlag && PHC) 
295     PHC->DrawAllHits();                           287     PHC->DrawAllHits();
296                                                   288 
297   // print this event by event (modulo n)         289   // print this event by event (modulo n)   
298   if (event_id%printModulo == 0)                  290   if (event_id%printModulo == 0) 
299     G4cout << "\n---> End of event: " << event    291     G4cout << "\n---> End of event: " << event_id << G4endl << G4endl;  
300                                                   292 
301 }                                                 293 }
302                                                   294 
303                                                   295 
304                                                   296 
305 //....oooOO0OOooo........oooOO0OOooo........oo    297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306 void DMXEventAction::writeScintHitsToFile()    << 298 void DMXEventAction::writeScintHitsToFile(void) {
307 {                                              << 
308                                                   299 
309   G4String filename="hits.out";                << 300   //  G4String filename="hits.out";
310   if (runAct)                                  << 301   G4String filename=runAct->GetsavehitsFile();
311     filename=runAct->GetsavehitsFile();        << 302   std::ofstream hitsfile(filename, std::ios::app);
312                                                << 303   if(!event_id) {
313                                                << 304     std::ofstream hitsfile(filename);
314                                                << 305     hitsfile <<"Evt     Eprim   Etot    LXe     LXeTime PMT     PMTTime Seed1           Seed2           First   Flags" 
315   //First time it is inkoved                   << 306        << G4endl;
316   if (!hitsfile)                               << 307     hitsfile <<"#       MeV     MeV     hits    ns      hits    ns                                      hit"
317     {                                          << 308        << G4endl
318       //check that we are in a worker: returns << 309        << G4endl;
319       //one file per thread is produced ending << 310   }
320       if (G4Threading::G4GetThreadId() >= 0)   << 
321   {                                            << 
322     std::stringstream sss;                     << 
323     sss << filename.c_str() << "." << G4Thread << 
324     filename = sss.str();                      << 
325     //G4cout << "Filename is: " << filename << << 
326   }                                            << 
327                                                << 
328       hitsfile = new std::ofstream;            << 
329       hitsfile->open(filename);                << 
330       (*hitsfile) <<"Evt     Eprim   Etot    L << 
331          << G4endl;                            << 
332       (*hitsfile) <<"#       MeV     MeV     h << 
333          << G4endl                             << 
334          << G4endl;                            << 
335     }                                          << 
336                                                   311 
337   if(S_hits) {                                    312   if(S_hits) {
338                                                   313 
339     if(hitsfile->is_open()) {                  << 314     if(hitsfile.is_open()) {
340                                                   315 
341                                                   316 
342       (*hitsfile) << std::setiosflags(std::ios << 317       hitsfile << std::setiosflags(std::ios::fixed)
343       << std::setprecision(4)                  << 318          << std::setprecision(4)
344       << std::setiosflags(std::ios::left)      << 319          << std::setiosflags(std::ios::left)
345       << std::setw(6)                          << 320          << std::setw(6)
346       << event_id << "\t"                      << 321          << event_id << "\t"
347       << energy_pri/MeV << "\t"                << 322          << energy_pri/MeV << "\t" 
348       << totEnergy/MeV << "\t"                 << 323          << totEnergy/MeV << "\t"
349       << S_hits  << "\t"                       << 324          << S_hits  << "\t"
350       << std::setiosflags(std::ios::scientific << 325          << std::setiosflags(std::ios::scientific) 
351       << std::setprecision(2)                  << 326          << std::setprecision(2)
352       << firstLXeHitTime/nanosecond << "\t"    << 327          << firstLXeHitTime/nanosecond << "\t"
353       << P_hits << "\t"                        << 328          << P_hits << "\t"
354       << std::setiosflags(std::ios::fixed)     << 329          << std::setiosflags(std::ios::fixed) 
355       << std::setprecision(4)                  << 330          << std::setprecision(4)
356       << aveTimePmtHits/nanosecond << "\t"     << 331          << aveTimePmtHits/nanosecond << "\t"
357       << *seeds     << "\t"                    << 332          << *seeds     << "\t"
358       << *(seeds+1) << "\t"                    << 333          << *(seeds+1) << "\t"
359       << firstParticleName << "\t"             << 334          << firstParticleName << "\t"
360       << (gamma_ev    ? "gamma " : "")         << 335          << (gamma_ev    ? "gamma " : "") 
361       << (neutron_ev  ? "neutron " : "")       << 336          << (neutron_ev  ? "neutron " : "") 
362       << (positron_ev ? "positron " : "")      << 337          << (positron_ev ? "positron " : "") 
363       << (electron_ev ? "electron " : "")      << 338          << (electron_ev ? "electron " : "") 
364       << (other_ev    ? "other " : "")         << 339          << (other_ev    ? "other " : "") 
365       << G4endl;                               << 340          << G4endl;
366                                                   341 
367       if (event_id%printModulo == 0)              342       if (event_id%printModulo == 0)
368   G4cout << "     Event summary in file " << f    343   G4cout << "     Event summary in file " << filename << G4endl;  
                                                   >> 344       hitsfile.close();
369     }                                             345     }
370                                                   346 
371     G4AnalysisManager* man = G4AnalysisManager << 347 #ifdef G4ANALYSIS_USE
372     G4int firstparticleIndex = 0;              << 
373     if(firstParticleName == "gamma") firstpart << 
374     else if (firstParticleName == "neutron") f << 
375     else if(firstParticleName == "electron") f << 
376     else if(firstParticleName == "positron") f << 
377     else{                                      << 
378       firstparticleIndex = 5;                  << 
379       man->FillH1(3,totEnergy);                << 
380     }                                          << 
381                                                << 
382     man->FillH1(4,P_hits,10); //weight         << 
383     man->FillH1(5,P_hits);                     << 
384                                                << 
385     man->FillH1(1,energy_pri/keV);             << 
386     man->FillH1(2,totEnergy/keV);              << 
387     man->FillH1(6,aveTimePmtHits/ns);          << 
388                                                << 
389     long seed1 = *seeds;                          348     long seed1 = *seeds;
390     long seed2 = *(seeds+1);                      349     long seed2 = *(seeds+1);    
391                                                << 350     // pass event summary to analysis manager for booking into histos and ntple
392     //Fill ntuple #2                           << 351     DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
393     man->FillNtupleDColumn(2,0,event_id);      << 352     analysis->analyseScintHits(event_id,energy_pri,totEnergy,S_hits,firstLXeHitTime,P_hits,aveTimePmtHits,firstParticleName,firstParticleE,gamma_ev,neutron_ev,positron_ev,electron_ev,other_ev,seed1,seed2);
394     man->FillNtupleDColumn(2,1,energy_pri/keV) << 353 #endif
395     man->FillNtupleDColumn(2,2,totEnergy);     << 
396     man->FillNtupleDColumn(2,3,S_hits);        << 
397     man->FillNtupleDColumn(2,4,firstLXeHitTime << 
398     man->FillNtupleDColumn(2,5,P_hits);        << 
399     man->FillNtupleDColumn(2,6,aveTimePmtHits) << 
400     man->FillNtupleDColumn(2,7,firstparticleIn << 
401     man->FillNtupleDColumn(2,8,firstParticleE) << 
402     man->FillNtupleDColumn(2,9,gamma_ev);      << 
403     man->FillNtupleDColumn(2,10,neutron_ev);   << 
404     man->FillNtupleDColumn(2,11,positron_ev);  << 
405     man->FillNtupleDColumn(2,12,electron_ev);  << 
406     man->FillNtupleDColumn(2,13,other_ev);     << 
407     man->FillNtupleDColumn(2,14,seed1);        << 
408     man->FillNtupleDColumn(2,15,seed2);        << 
409     man->AddNtupleRow(2);                      << 
410                                                << 
411   }                                               354   }
412                                                   355 
413 }                                                 356 }
414                                                   357 
415                                                   358 
416 //....oooOO0OOooo........oooOO0OOooo........oo    359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417 void DMXEventAction::writePmtHitsToFile(const  << 360 void DMXEventAction::writePmtHitsToFile(const DMXPmtHitsCollection* hits) {
418 {                                              << 
419   G4String filename="pmt.out";                 << 
420   if (runAct)                                  << 
421     filename=runAct->GetsavepmtFile();         << 
422                                                << 
423                                                << 
424   //first time it is invoked: create it        << 
425   if (!pmtfile)                                << 
426     {                                          << 
427       //check that we are in a worker: returns << 
428       //one file per thread is produced ending << 
429       if (G4Threading::G4GetThreadId() >= 0)   << 
430   {                                            << 
431     std::stringstream sss;                     << 
432     sss << filename.c_str() << "." << G4Thread << 
433     filename = sss.str();                      << 
434     //G4cout << "Filename is: " << filename << << 
435   }                                            << 
436       pmtfile = new std::ofstream;             << 
437       pmtfile->open(filename);                 << 
438     }                                          << 
439                                                << 
440                                                   361 
                                                   >> 362   //  G4String filename="pmt.out";
                                                   >> 363   G4String filename=runAct->GetsavepmtFile();
                                                   >> 364   std::ofstream pmtfile(filename, std::ios::app);
441   G4double x; G4double y; G4double z;             365   G4double x; G4double y; G4double z;
442   G4AnalysisManager* man = G4AnalysisManager:: << 
443                                                   366 
444   if(pmtfile->is_open()) {                     << 367   if(pmtfile.is_open()) {
445     (*pmtfile) << "Hit#    X, mm   Y, mm   Z,  << 368     pmtfile << "Hit#    X, mm   Y, mm   Z, mm" << G4endl;       
446     (*pmtfile) << std::setiosflags(std::ios::f << 369     pmtfile << std::setiosflags(std::ios::fixed)
447          << std::setprecision(3)               << 370       << std::setprecision(3)
448          << std::setiosflags(std::ios::left)   << 371       << std::setiosflags(std::ios::left)
449          << std::setw(6);                      << 372       << std::setw(6);
450     for (G4int i=0; i<P_hits; i++)                373     for (G4int i=0; i<P_hits; i++)
451       {                                           374       {
452   x = ((*hits)[i]->GetPos()).x()/mm;              375   x = ((*hits)[i]->GetPos()).x()/mm;
453   y = ((*hits)[i]->GetPos()).y()/mm;              376   y = ((*hits)[i]->GetPos()).y()/mm;
454   z = ((*hits)[i]->GetPos()).z()/mm;              377   z = ((*hits)[i]->GetPos()).z()/mm;
455   (*pmtfile) << i << "\t"                      << 378   pmtfile << i << "\t"
456        << x << "\t"                            << 379     << x << "\t" 
457        << y << "\t"                            << 380     << y << "\t"
458        << z << G4endl;                         << 381     << z << G4endl;
459                                                << 382 #ifdef G4ANALYSIS_USE
460   man->FillH2(1,x/mm, y/mm);  // fill(x,y)     << 383     // pass pmt hit summary to analysis manager for booking in ntples
461   if (event_id == 0 ) {                        << 384   DMXAnalysisManager* analysis = DMXAnalysisManager::getInstance();
462     man->FillH2(2,x,y);                        << 385   analysis->analysePMTHits(event_id,i,x,y,z);
463   }                                            << 386 #endif
464                                                << 
465   //Fill ntuple #3                             << 
466   man->FillNtupleDColumn(3,0,event_id);        << 
467   man->FillNtupleDColumn(3,1,(G4double) i);    << 
468   man->FillNtupleDColumn(3,2,x);               << 
469   man->FillNtupleDColumn(3,3,y);               << 
470   man->FillNtupleDColumn(3,4,z);               << 
471   man->AddNtupleRow(3);                        << 
472                                                << 
473       }                                           387       }
474     if (event_id%printModulo == 0 && P_hits >     388     if (event_id%printModulo == 0 && P_hits > 0) 
475       G4cout << "     " << P_hits << " PMT hit    389       G4cout << "     " << P_hits << " PMT hits in " << filename << G4endl;  
                                                   >> 390     pmtfile.close();
476   }                                               391   }
477                                                   392   
478 }                                                 393 }
479                                                   394 
480                                                   395 
481                                                   396 
482 //....oooOO0OOooo........oooOO0OOooo........oo    397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483 void DMXEventAction::drawTracks(const G4Event*    398 void DMXEventAction::drawTracks(const G4Event* evt) {
484                                                   399 
485   if(G4VVisManager::GetConcreteInstance()) {      400   if(G4VVisManager::GetConcreteInstance()) {
486     G4UImanager::GetUIpointer()->ApplyCommand(    401     G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers");    
487     G4TrajectoryContainer* trajContainer = evt    402     G4TrajectoryContainer* trajContainer = evt->GetTrajectoryContainer();
488     G4int n_trajectories = 0;                     403     G4int n_trajectories = 0;
489                                                   404 
490     if(trajContainer) n_trajectories = trajCon    405     if(trajContainer) n_trajectories = trajContainer->entries();
491     for (G4int i=0; i<n_trajectories; i++) {      406     for (G4int i=0; i<n_trajectories; i++) {
492       G4Trajectory* trj = (G4Trajectory*)(*tra    407       G4Trajectory* trj = (G4Trajectory*)(*trajContainer)[i];
493       if (drawTrksFlag == "all")                  408       if (drawTrksFlag == "all") 
494   trj->DrawTrajectory();                          409   trj->DrawTrajectory();
495       else if ((drawTrksFlag == "charged") &&     410       else if ((drawTrksFlag == "charged") && (trj->GetCharge() != 0.))
496   trj->DrawTrajectory();                          411   trj->DrawTrajectory();
497       else if ((drawTrksFlag == "noscint")        412       else if ((drawTrksFlag == "noscint") 
498          && (trj->GetParticleName() != "optica    413          && (trj->GetParticleName() != "opticalphoton"))
499   trj->DrawTrajectory();                          414   trj->DrawTrajectory();
500     }                                             415     }
501                                                   416     
502     // G4UImanager::GetUIpointer()->ApplyComma    417     // G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");    
503   }                                               418   } 
504                                                   419 
505 }                                                 420 }
506                                                   421