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