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 11.0)


  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 #include "G4AnalysisManager.hh"
 57 #include "G4Event.hh"                              57 #include "G4Event.hh"
 58 #include "G4EventManager.hh"                       58 #include "G4EventManager.hh"
 59 #include "G4HCofThisEvent.hh"                      59 #include "G4HCofThisEvent.hh"
 60 #include "G4VHitsCollection.hh"                    60 #include "G4VHitsCollection.hh"
 61 #include "G4TrajectoryContainer.hh"                61 #include "G4TrajectoryContainer.hh"
 62 #include "G4Trajectory.hh"                         62 #include "G4Trajectory.hh"
 63 #include "G4VVisManager.hh"                        63 #include "G4VVisManager.hh"
 64 #include "G4SDManager.hh"                          64 #include "G4SDManager.hh"
 65 #include "G4UImanager.hh"                          65 #include "G4UImanager.hh"
 66 #include "G4UnitsTable.hh"                         66 #include "G4UnitsTable.hh"
 67 #include "G4RunManager.hh"                         67 #include "G4RunManager.hh"
 68 #include "G4Threading.hh"                          68 #include "G4Threading.hh"
 69 #include <fstream>                                 69 #include <fstream>
 70 #include <iomanip>                                 70 #include <iomanip>
 71                                                    71 
 72 #include "G4SystemOfUnits.hh"                      72 #include "G4SystemOfUnits.hh"
 73                                                    73 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75 DMXEventAction::DMXEventAction()                   75 DMXEventAction::DMXEventAction() 
 76   : runAct(0),genAction(0),hitsfile(0),pmtfile     76   : runAct(0),genAction(0),hitsfile(0),pmtfile(0)
 77 {                                                  77 {
 78                                                    78 
 79   // create messenger                              79   // create messenger
 80   eventMessenger = new DMXEventActionMessenger     80   eventMessenger = new DMXEventActionMessenger(this);
 81                                                    81 
 82   // defaults for messenger                        82   // defaults for messenger
 83   drawColsFlag = "standard";                       83   drawColsFlag = "standard";
 84   drawTrksFlag = "all";                            84   drawTrksFlag = "all";
 85   drawHitsFlag = 1;                                85   drawHitsFlag = 1;
 86   savePmtFlag  = 0;                                86   savePmtFlag  = 0;
 87   saveHitsFlag = 1;                                87   saveHitsFlag = 1;
 88                                                    88   
 89   printModulo = 1;                                 89   printModulo = 1;
 90                                                    90 
 91   // hits collections                              91   // hits collections
 92   scintillatorCollID = -1;                         92   scintillatorCollID = -1;
 93   pmtCollID = -1;                                  93   pmtCollID = -1;
 94                                                    94 
 95   energy_pri=0;                                    95   energy_pri=0;
 96   seeds=NULL;                                      96   seeds=NULL;
 97                                                    97 
 98 }                                                  98 }
 99                                                    99 
100                                                   100 
101 //....oooOO0OOooo........oooOO0OOooo........oo    101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 DMXEventAction::~DMXEventAction() {               102 DMXEventAction::~DMXEventAction() {
103                                                   103 
104   if (hitsfile)                                   104   if (hitsfile)
105     {                                             105     {
106       hitsfile->close();                          106       hitsfile->close();
107       delete hitsfile;                            107       delete hitsfile;
108     }                                             108     }
109   if (pmtfile)                                    109   if (pmtfile)
110     {                                             110     {
111       pmtfile->close();                           111       pmtfile->close();
112       delete pmtfile;                             112       delete pmtfile;
113     }                                             113     }
114   delete eventMessenger;                          114   delete eventMessenger;
115 }                                                 115 }
116                                                   116 
117                                                   117 
118 //....oooOO0OOooo........oooOO0OOooo........oo    118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 void DMXEventAction::BeginOfEventAction(const     119 void DMXEventAction::BeginOfEventAction(const G4Event* evt) 
120 {                                                 120 {
121                                                   121 
122   //thread-local run action                       122   //thread-local run action
123   if (!runAct)                                    123   if (!runAct) 
124     runAct =                                      124     runAct = 
125       dynamic_cast<const DMXRunAction*>           125       dynamic_cast<const DMXRunAction*>
126       (G4RunManager::GetRunManager()->GetUserR    126       (G4RunManager::GetRunManager()->GetUserRunAction());
127                                                   127   
128   if (!genAction)                                 128   if (!genAction)
129     genAction = dynamic_cast<const DMXPrimaryG    129     genAction = dynamic_cast<const DMXPrimaryGeneratorAction*>
130       (G4RunManager::GetRunManager()->GetUserP    130       (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
131                                                   131 
132                                                   132 
133   // grab seeds                                   133   // grab seeds
134   seeds = genAction->GetEventSeeds();             134   seeds = genAction->GetEventSeeds();
135                                                   135 
136   // grab energy of primary                       136   // grab energy of primary
137   energy_pri = genAction->GetEnergyPrimary();     137   energy_pri = genAction->GetEnergyPrimary();
138                                                   138 
139   event_id = evt->GetEventID();                   139   event_id = evt->GetEventID();
140                                                   140  
141   // print this information event by event (mo    141   // print this information event by event (modulo n)   
142   if (event_id%printModulo == 0)                  142   if (event_id%printModulo == 0)
143     {                                             143     {
144       G4cout << "\n---> Begin of event: " << e    144       G4cout << "\n---> Begin of event: " << event_id << G4endl;
145       G4cout << "\n     Primary Energy: " << G    145       G4cout << "\n     Primary Energy: " << G4BestUnit(energy_pri,"Energy") 
146        << G4endl;                                 146        << G4endl;
147       //      HepRandom::showEngineStatus();      147       //      HepRandom::showEngineStatus(); 
148     }                                             148     }
149                                                   149 
150                                                   150 
151   // get ID for scintillator hits collection      151   // get ID for scintillator hits collection
152   if (scintillatorCollID==-1) {                   152   if (scintillatorCollID==-1) {
153     G4SDManager *SDman = G4SDManager::GetSDMpo    153     G4SDManager *SDman = G4SDManager::GetSDMpointer();
154     scintillatorCollID = SDman->GetCollectionI    154     scintillatorCollID = SDman->GetCollectionID("scintillatorCollection");
155   }                                               155   }
156                                                   156 
157   // get ID for pmt hits collection               157   // get ID for pmt hits collection
158   if (pmtCollID==-1) {                            158   if (pmtCollID==-1) {
159     G4SDManager *SDman = G4SDManager::GetSDMpo    159     G4SDManager *SDman = G4SDManager::GetSDMpointer();
160     pmtCollID = SDman->GetCollectionID("pmtCol    160     pmtCollID = SDman->GetCollectionID("pmtCollection");
161   }                                               161   }
162                                                   162 
163 }                                                 163 }
164                                                   164 
165                                                   165 
166 //....oooOO0OOooo........oooOO0OOooo........oo    166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167 void DMXEventAction::EndOfEventAction(const G4    167 void DMXEventAction::EndOfEventAction(const G4Event* evt) {
168                                                   168 
169   // check that both hits collections have bee    169   // check that both hits collections have been defined
170   if(scintillatorCollID<0||pmtCollID<0) return    170   if(scintillatorCollID<0||pmtCollID<0) return;
171                                                   171 
172   G4AnalysisManager* man = G4AnalysisManager::    172   G4AnalysisManager* man = G4AnalysisManager::Instance();
173                                                   173 
174   // address hits collections                     174   // address hits collections
175   DMXScintHitsCollection* SHC = NULL;             175   DMXScintHitsCollection* SHC = NULL;
176   DMXPmtHitsCollection*   PHC = NULL;             176   DMXPmtHitsCollection*   PHC = NULL;
177   G4HCofThisEvent* HCE = evt->GetHCofThisEvent    177   G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
178   if(HCE) {                                       178   if(HCE) {
179     SHC = (DMXScintHitsCollection*)(HCE->GetHC    179     SHC = (DMXScintHitsCollection*)(HCE->GetHC(scintillatorCollID));
180     PHC = (DMXPmtHitsCollection*)(HCE->GetHC(p    180     PHC = (DMXPmtHitsCollection*)(HCE->GetHC(pmtCollID));
181   }                                               181   }
182                                                   182 
183   // event summary                                183   // event summary
184   totEnergy         = 0.;                         184   totEnergy         = 0.;
185   totEnergyGammas   = 0.;                         185   totEnergyGammas   = 0.;
186   totEnergyNeutrons = 0.;                         186   totEnergyNeutrons = 0.;
187   firstParticleE    = 0.;                         187   firstParticleE    = 0.;
188   particleEnergy    = 0.;                         188   particleEnergy    = 0.;
189   firstLXeHitTime   = 0.;                         189   firstLXeHitTime   = 0.;
190   aveTimePmtHits    = 0.;                         190   aveTimePmtHits    = 0.;
191                                                   191 
192   firstParticleName = "";                         192   firstParticleName = "";
193   particleName      = "";                         193   particleName      = "";
194                                                   194 
195                                                   195 
196   // particle flags                               196   // particle flags
197   gamma_ev          = false;                      197   gamma_ev          = false;
198   neutron_ev        = false;                      198   neutron_ev        = false;
199   positron_ev       = false;                      199   positron_ev       = false;
200   electron_ev       = false;                      200   electron_ev       = false;
201   proton_ev         = false;                      201   proton_ev         = false;
202   other_ev          = false;                      202   other_ev          = false;
203   start_gamma       = false;                      203   start_gamma       = false;
204   start_neutron     = false;                      204   start_neutron     = false;
205                                                   205 
206                                                   206 
207   // scintillator hits                            207   // scintillator hits
208   if(SHC) {                                       208   if(SHC) {
209     S_hits = SHC->entries();                      209     S_hits = SHC->entries();
210                                                   210     
211     for (G4int i=0; i<S_hits; i++) {              211     for (G4int i=0; i<S_hits; i++) {
212       if(i==0) {                                  212       if(i==0) {
213   firstParticleName = (*SHC)[0]->GetParticle()    213   firstParticleName = (*SHC)[0]->GetParticle();
214   firstLXeHitTime   = (*SHC)[0]->GetTime();       214   firstLXeHitTime   = (*SHC)[0]->GetTime();
215   firstParticleE = (*SHC)[0]->GetParticleEnerg    215   firstParticleE = (*SHC)[0]->GetParticleEnergy();
216   if (event_id%printModulo == 0 && S_hits > 0)    216   if (event_id%printModulo == 0 && S_hits > 0) {
217     G4cout << "     First hit in LXe: " << fir    217     G4cout << "     First hit in LXe: " << firstParticleName << G4endl;
218     G4cout << "     Number of hits in LXe: " <    218     G4cout << "     Number of hits in LXe: " << S_hits << G4endl; 
219   }                                               219   }
220       }                                           220       }
221       hitEnergy         = (*SHC)[i]->GetEdep()    221       hitEnergy         = (*SHC)[i]->GetEdep();
222       totEnergy        += hitEnergy;              222       totEnergy        += hitEnergy;
223                                                   223       
224       particleName      = (*SHC)[i]->GetPartic    224       particleName      = (*SHC)[i]->GetParticle();
225       particleEnergy    = (*SHC)[i]->GetPartic    225       particleEnergy    = (*SHC)[i]->GetParticleEnergy();
226                                                   226 
227       if(particleName == "gamma") {               227       if(particleName == "gamma") {
228   gamma_ev = true;                                228   gamma_ev = true;
229   start_gamma = true;                             229   start_gamma = true;
230   start_neutron = false;                          230   start_neutron = false;
231       }                                           231       }
232       else if(particleName == "neutron")          232       else if(particleName == "neutron") 
233   neutron_ev = true;                              233   neutron_ev = true;
234       else if(particleName == "e+")               234       else if(particleName == "e+") 
235   positron_ev = true;                             235   positron_ev = true;
236       else if(particleName == "e-")               236       else if(particleName == "e-") 
237   electron_ev = true;                             237   electron_ev = true;
238       else if(particleName == "proton")           238       else if(particleName == "proton") 
239   proton_ev = true;                               239   proton_ev = true;
240       else {                                      240       else {
241   other_ev = true;                                241   other_ev = true;
242   start_gamma = false;                            242   start_gamma = false;
243   start_neutron = true;                           243   start_neutron = true;
244       }                                           244       }
245                                                   245 
246       if(start_gamma && !start_neutron)           246       if(start_gamma && !start_neutron) 
247   totEnergyGammas += hitEnergy;                   247   totEnergyGammas += hitEnergy;
248       if(start_neutron && !start_gamma)           248       if(start_neutron && !start_gamma) 
249   totEnergyNeutrons += hitEnergy;                 249   totEnergyNeutrons += hitEnergy;
250     }                                             250     }
251                                                   251     
252     if (event_id%printModulo == 0)                252     if (event_id%printModulo == 0)
253       G4cout << "     Total energy in LXe: "      253       G4cout << "     Total energy in LXe: " 
254        << G4BestUnit(totEnergy,"Energy") << G4    254        << G4BestUnit(totEnergy,"Energy") << G4endl;  
255                                                   255     
256   }                                               256   }
257                                                   257   
258                                                   258   
259   // PMT hits                                     259   // PMT hits
260   if(PHC) {                                       260   if(PHC) {
261     P_hits = PHC->entries();                      261     P_hits = PHC->entries();
262                                                   262     
263     // average time of PMT hits                   263     // average time of PMT hits
264     for (G4int i=0; i<P_hits; i++) {              264     for (G4int i=0; i<P_hits; i++) {
265       G4double time = ( (*PHC)[i]->GetTime() -    265       G4double time = ( (*PHC)[i]->GetTime() - firstLXeHitTime );
266       aveTimePmtHits += time / (G4double)P_hit    266       aveTimePmtHits += time / (G4double)P_hits;
267       ////      if (event_id == 0) {              267       ////      if (event_id == 0) {
268       if(P_hits >= 0) {                           268       if(P_hits >= 0) {
269   man->FillH1(7,time);                            269   man->FillH1(7,time);
270       }                                           270       }
271     }                                             271     }
272                                                   272   
273     if (event_id%printModulo == 0 && P_hits >     273     if (event_id%printModulo == 0 && P_hits > 0) {
274       G4cout << "     Average light collection    274       G4cout << "     Average light collection time: "
275        << G4BestUnit(aveTimePmtHits,"Time") <<    275        << G4BestUnit(aveTimePmtHits,"Time") << G4endl;
276       G4cout << "     Number of PMT hits (phot    276       G4cout << "     Number of PMT hits (photoelectron equivalent): " 
277        << P_hits << G4endl;                       277        << P_hits << G4endl;     
278     }                                             278     }
279     // write out (x,y,z) of PMT hits              279     // write out (x,y,z) of PMT hits
280     if (savePmtFlag)                              280     if (savePmtFlag)
281       writePmtHitsToFile(PHC);                    281       writePmtHitsToFile(PHC);
282   }                                               282   }
283                                                   283   
284                                                   284 
285   // write out event summary                      285   // write out event summary
286   if(saveHitsFlag)                                286   if(saveHitsFlag) 
287     writeScintHitsToFile();                       287     writeScintHitsToFile();
288                                                   288   
289   // draw trajectories                            289   // draw trajectories
290   if(drawColsFlag=="standard" && drawTrksFlag!    290   if(drawColsFlag=="standard" && drawTrksFlag!="none")
291     drawTracks(evt);                              291     drawTracks(evt);
292                                                   292 
293   // hits in PMT                                  293   // hits in PMT
294   if(drawHitsFlag && PHC)                         294   if(drawHitsFlag && PHC) 
295     PHC->DrawAllHits();                           295     PHC->DrawAllHits();
296                                                   296 
297   // print this event by event (modulo n)         297   // print this event by event (modulo n)   
298   if (event_id%printModulo == 0)                  298   if (event_id%printModulo == 0) 
299     G4cout << "\n---> End of event: " << event    299     G4cout << "\n---> End of event: " << event_id << G4endl << G4endl;  
300                                                   300 
301 }                                                 301 }
302                                                   302 
303                                                   303 
304                                                   304 
305 //....oooOO0OOooo........oooOO0OOooo........oo    305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306 void DMXEventAction::writeScintHitsToFile()       306 void DMXEventAction::writeScintHitsToFile() 
307 {                                                 307 {
308                                                   308 
309   G4String filename="hits.out";                   309   G4String filename="hits.out";
310   if (runAct)                                     310   if (runAct)
311     filename=runAct->GetsavehitsFile();           311     filename=runAct->GetsavehitsFile();
312                                                   312 
313                                                   313  
314                                                   314   
315   //First time it is inkoved                      315   //First time it is inkoved
316   if (!hitsfile)                                  316   if (!hitsfile)
317     {                                             317     {
318       //check that we are in a worker: returns    318       //check that we are in a worker: returns -1 in a master and -2 in sequential
319       //one file per thread is produced ending    319       //one file per thread is produced ending with ".N", with N= thread number
320       if (G4Threading::G4GetThreadId() >= 0)      320       if (G4Threading::G4GetThreadId() >= 0)
321   {                                               321   {
322     std::stringstream sss;                        322     std::stringstream sss;
323     sss << filename.c_str() << "." << G4Thread    323     sss << filename.c_str() << "." << G4Threading::G4GetThreadId();  
324     filename = sss.str();                         324     filename = sss.str();
325     //G4cout << "Filename is: " << filename <<    325     //G4cout << "Filename is: " << filename << G4endl;
326   }                                               326   }
327                                                   327       
328       hitsfile = new std::ofstream;               328       hitsfile = new std::ofstream;
329       hitsfile->open(filename);                   329       hitsfile->open(filename);
330       (*hitsfile) <<"Evt     Eprim   Etot    L    330       (*hitsfile) <<"Evt     Eprim   Etot    LXe     LXeTime PMT     PMTTime Seed1           Seed2           First   Flags" 
331          << G4endl;                               331          << G4endl;
332       (*hitsfile) <<"#       MeV     MeV     h    332       (*hitsfile) <<"#       MeV     MeV     hits    ns      hits    ns                                      hit"
333          << G4endl                                333          << G4endl
334          << G4endl;                               334          << G4endl;
335     }                                             335     }
336                                                   336 
337   if(S_hits) {                                    337   if(S_hits) {
338                                                   338 
339     if(hitsfile->is_open()) {                     339     if(hitsfile->is_open()) {
340                                                   340 
341                                                   341 
342       (*hitsfile) << std::setiosflags(std::ios    342       (*hitsfile) << std::setiosflags(std::ios::fixed)
343       << std::setprecision(4)                     343       << std::setprecision(4)
344       << std::setiosflags(std::ios::left)         344       << std::setiosflags(std::ios::left)
345       << std::setw(6)                             345       << std::setw(6)
346       << event_id << "\t"                         346       << event_id << "\t"
347       << energy_pri/MeV << "\t"                   347       << energy_pri/MeV << "\t" 
348       << totEnergy/MeV << "\t"                    348       << totEnergy/MeV << "\t"
349       << S_hits  << "\t"                          349       << S_hits  << "\t"
350       << std::setiosflags(std::ios::scientific    350       << std::setiosflags(std::ios::scientific) 
351       << std::setprecision(2)                     351       << std::setprecision(2)
352       << firstLXeHitTime/nanosecond << "\t"       352       << firstLXeHitTime/nanosecond << "\t"
353       << P_hits << "\t"                           353       << P_hits << "\t"
354       << std::setiosflags(std::ios::fixed)        354       << std::setiosflags(std::ios::fixed) 
355       << std::setprecision(4)                     355       << std::setprecision(4)
356       << aveTimePmtHits/nanosecond << "\t"        356       << aveTimePmtHits/nanosecond << "\t"
357       << *seeds     << "\t"                       357       << *seeds     << "\t"
358       << *(seeds+1) << "\t"                       358       << *(seeds+1) << "\t"
359       << firstParticleName << "\t"                359       << firstParticleName << "\t"
360       << (gamma_ev    ? "gamma " : "")            360       << (gamma_ev    ? "gamma " : "") 
361       << (neutron_ev  ? "neutron " : "")          361       << (neutron_ev  ? "neutron " : "") 
362       << (positron_ev ? "positron " : "")         362       << (positron_ev ? "positron " : "") 
363       << (electron_ev ? "electron " : "")         363       << (electron_ev ? "electron " : "") 
364       << (other_ev    ? "other " : "")            364       << (other_ev    ? "other " : "") 
365       << G4endl;                                  365       << G4endl;
366                                                   366 
367       if (event_id%printModulo == 0)              367       if (event_id%printModulo == 0)
368   G4cout << "     Event summary in file " << f    368   G4cout << "     Event summary in file " << filename << G4endl;  
369     }                                             369     }
370                                                   370 
371     G4AnalysisManager* man = G4AnalysisManager    371     G4AnalysisManager* man = G4AnalysisManager::Instance();
372     G4int firstparticleIndex = 0;                 372     G4int firstparticleIndex = 0;
373     if(firstParticleName == "gamma") firstpart    373     if(firstParticleName == "gamma") firstparticleIndex = 1;
374     else if (firstParticleName == "neutron") f    374     else if (firstParticleName == "neutron") firstparticleIndex = 2;
375     else if(firstParticleName == "electron") f    375     else if(firstParticleName == "electron") firstparticleIndex = 3;
376     else if(firstParticleName == "positron") f    376     else if(firstParticleName == "positron") firstparticleIndex = 4;
377     else{                                         377     else{
378       firstparticleIndex = 5;                     378       firstparticleIndex = 5;
379       man->FillH1(3,totEnergy);                   379       man->FillH1(3,totEnergy);
380     }                                             380     }
381                                                   381 
382     man->FillH1(4,P_hits,10); //weight            382     man->FillH1(4,P_hits,10); //weight
383     man->FillH1(5,P_hits);                        383     man->FillH1(5,P_hits);
384                                                   384 
385     man->FillH1(1,energy_pri/keV);                385     man->FillH1(1,energy_pri/keV);
386     man->FillH1(2,totEnergy/keV);                 386     man->FillH1(2,totEnergy/keV);
387     man->FillH1(6,aveTimePmtHits/ns);             387     man->FillH1(6,aveTimePmtHits/ns);
388                                                   388 
389     long seed1 = *seeds;                          389     long seed1 = *seeds;
390     long seed2 = *(seeds+1);                      390     long seed2 = *(seeds+1);    
391                                                   391     
392     //Fill ntuple #2                              392     //Fill ntuple #2
393     man->FillNtupleDColumn(2,0,event_id);         393     man->FillNtupleDColumn(2,0,event_id);
394     man->FillNtupleDColumn(2,1,energy_pri/keV)    394     man->FillNtupleDColumn(2,1,energy_pri/keV);
395     man->FillNtupleDColumn(2,2,totEnergy);        395     man->FillNtupleDColumn(2,2,totEnergy);
396     man->FillNtupleDColumn(2,3,S_hits);           396     man->FillNtupleDColumn(2,3,S_hits);
397     man->FillNtupleDColumn(2,4,firstLXeHitTime    397     man->FillNtupleDColumn(2,4,firstLXeHitTime);
398     man->FillNtupleDColumn(2,5,P_hits);           398     man->FillNtupleDColumn(2,5,P_hits);
399     man->FillNtupleDColumn(2,6,aveTimePmtHits)    399     man->FillNtupleDColumn(2,6,aveTimePmtHits);
400     man->FillNtupleDColumn(2,7,firstparticleIn    400     man->FillNtupleDColumn(2,7,firstparticleIndex);
401     man->FillNtupleDColumn(2,8,firstParticleE)    401     man->FillNtupleDColumn(2,8,firstParticleE);
402     man->FillNtupleDColumn(2,9,gamma_ev);         402     man->FillNtupleDColumn(2,9,gamma_ev);
403     man->FillNtupleDColumn(2,10,neutron_ev);      403     man->FillNtupleDColumn(2,10,neutron_ev);
404     man->FillNtupleDColumn(2,11,positron_ev);     404     man->FillNtupleDColumn(2,11,positron_ev);
405     man->FillNtupleDColumn(2,12,electron_ev);     405     man->FillNtupleDColumn(2,12,electron_ev);
406     man->FillNtupleDColumn(2,13,other_ev);        406     man->FillNtupleDColumn(2,13,other_ev);
407     man->FillNtupleDColumn(2,14,seed1);           407     man->FillNtupleDColumn(2,14,seed1);
408     man->FillNtupleDColumn(2,15,seed2);           408     man->FillNtupleDColumn(2,15,seed2);
409     man->AddNtupleRow(2);                         409     man->AddNtupleRow(2);
410                                                   410     
411   }                                               411   }
412                                                   412 
413 }                                                 413 }
414                                                   414 
415                                                   415 
416 //....oooOO0OOooo........oooOO0OOooo........oo    416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417 void DMXEventAction::writePmtHitsToFile(const     417 void DMXEventAction::writePmtHitsToFile(const DMXPmtHitsCollection* hits) 
418 {                                                 418 { 
419   G4String filename="pmt.out";                    419   G4String filename="pmt.out";
420   if (runAct)                                     420   if (runAct)
421     filename=runAct->GetsavepmtFile();            421     filename=runAct->GetsavepmtFile();
422                                                   422   
423                                                   423 
424   //first time it is invoked: create it           424   //first time it is invoked: create it
425   if (!pmtfile)                                   425   if (!pmtfile)
426     {                                             426     {
427       //check that we are in a worker: returns    427       //check that we are in a worker: returns -1 in a master and -2 in sequential
428       //one file per thread is produced ending    428       //one file per thread is produced ending with ".N", with N= thread number
429       if (G4Threading::G4GetThreadId() >= 0)      429       if (G4Threading::G4GetThreadId() >= 0)
430   {                                               430   {
431     std::stringstream sss;                        431     std::stringstream sss;
432     sss << filename.c_str() << "." << G4Thread    432     sss << filename.c_str() << "." << G4Threading::G4GetThreadId();  
433     filename = sss.str();                         433     filename = sss.str();
434     //G4cout << "Filename is: " << filename <<    434     //G4cout << "Filename is: " << filename << G4endl;
435   }                                               435   }
436       pmtfile = new std::ofstream;                436       pmtfile = new std::ofstream;
437       pmtfile->open(filename);                    437       pmtfile->open(filename);
438     }                                             438     }
439                                                   439 
440                                                   440 
441   G4double x; G4double y; G4double z;             441   G4double x; G4double y; G4double z;
442   G4AnalysisManager* man = G4AnalysisManager::    442   G4AnalysisManager* man = G4AnalysisManager::Instance();
443                                                   443 
444   if(pmtfile->is_open()) {                        444   if(pmtfile->is_open()) {
445     (*pmtfile) << "Hit#    X, mm   Y, mm   Z,     445     (*pmtfile) << "Hit#    X, mm   Y, mm   Z, mm" << G4endl;       
446     (*pmtfile) << std::setiosflags(std::ios::f    446     (*pmtfile) << std::setiosflags(std::ios::fixed)
447          << std::setprecision(3)                  447          << std::setprecision(3)
448          << std::setiosflags(std::ios::left)      448          << std::setiosflags(std::ios::left)
449          << std::setw(6);                         449          << std::setw(6);
450     for (G4int i=0; i<P_hits; i++)                450     for (G4int i=0; i<P_hits; i++)
451       {                                           451       {
452   x = ((*hits)[i]->GetPos()).x()/mm;              452   x = ((*hits)[i]->GetPos()).x()/mm;
453   y = ((*hits)[i]->GetPos()).y()/mm;              453   y = ((*hits)[i]->GetPos()).y()/mm;
454   z = ((*hits)[i]->GetPos()).z()/mm;              454   z = ((*hits)[i]->GetPos()).z()/mm;
455   (*pmtfile) << i << "\t"                         455   (*pmtfile) << i << "\t"
456        << x << "\t"                               456        << x << "\t" 
457        << y << "\t"                               457        << y << "\t"
458        << z << G4endl;                            458        << z << G4endl;
459                                                   459   
460   man->FillH2(1,x/mm, y/mm);  // fill(x,y)        460   man->FillH2(1,x/mm, y/mm);  // fill(x,y)
461   if (event_id == 0 ) {                           461   if (event_id == 0 ) {
462     man->FillH2(2,x,y);                           462     man->FillH2(2,x,y);
463   }                                               463   }
464                                                   464 
465   //Fill ntuple #3                                465   //Fill ntuple #3
466   man->FillNtupleDColumn(3,0,event_id);           466   man->FillNtupleDColumn(3,0,event_id);
467   man->FillNtupleDColumn(3,1,(G4double) i);       467   man->FillNtupleDColumn(3,1,(G4double) i);
468   man->FillNtupleDColumn(3,2,x);                  468   man->FillNtupleDColumn(3,2,x);
469   man->FillNtupleDColumn(3,3,y);                  469   man->FillNtupleDColumn(3,3,y);
470   man->FillNtupleDColumn(3,4,z);                  470   man->FillNtupleDColumn(3,4,z);
471   man->AddNtupleRow(3);                           471   man->AddNtupleRow(3);
472                                                   472 
473       }                                           473       }
474     if (event_id%printModulo == 0 && P_hits >     474     if (event_id%printModulo == 0 && P_hits > 0) 
475       G4cout << "     " << P_hits << " PMT hit    475       G4cout << "     " << P_hits << " PMT hits in " << filename << G4endl;  
476   }                                               476   }
477                                                   477   
478 }                                                 478 }
479                                                   479 
480                                                   480 
481                                                   481 
482 //....oooOO0OOooo........oooOO0OOooo........oo    482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483 void DMXEventAction::drawTracks(const G4Event*    483 void DMXEventAction::drawTracks(const G4Event* evt) {
484                                                   484 
485   if(G4VVisManager::GetConcreteInstance()) {      485   if(G4VVisManager::GetConcreteInstance()) {
486     G4UImanager::GetUIpointer()->ApplyCommand(    486     G4UImanager::GetUIpointer()->ApplyCommand("/vis/scene/notifyHandlers");    
487     G4TrajectoryContainer* trajContainer = evt    487     G4TrajectoryContainer* trajContainer = evt->GetTrajectoryContainer();
488     G4int n_trajectories = 0;                     488     G4int n_trajectories = 0;
489                                                   489 
490     if(trajContainer) n_trajectories = trajCon    490     if(trajContainer) n_trajectories = trajContainer->entries();
491     for (G4int i=0; i<n_trajectories; i++) {      491     for (G4int i=0; i<n_trajectories; i++) {
492       G4Trajectory* trj = (G4Trajectory*)(*tra    492       G4Trajectory* trj = (G4Trajectory*)(*trajContainer)[i];
493       if (drawTrksFlag == "all")                  493       if (drawTrksFlag == "all") 
494   trj->DrawTrajectory();                          494   trj->DrawTrajectory();
495       else if ((drawTrksFlag == "charged") &&     495       else if ((drawTrksFlag == "charged") && (trj->GetCharge() != 0.))
496   trj->DrawTrajectory();                          496   trj->DrawTrajectory();
497       else if ((drawTrksFlag == "noscint")        497       else if ((drawTrksFlag == "noscint") 
498          && (trj->GetParticleName() != "optica    498          && (trj->GetParticleName() != "opticalphoton"))
499   trj->DrawTrajectory();                          499   trj->DrawTrajectory();
500     }                                             500     }
501                                                   501     
502     // G4UImanager::GetUIpointer()->ApplyComma    502     // G4UImanager::GetUIpointer()->ApplyCommand("/vis/viewer/update");    
503   }                                               503   } 
504                                                   504 
505 }                                                 505 }
506                                                   506