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