Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // 28 // Author: A. Pfeiffer (Andreas.Pfeiffer@cern 29 // (copied from his UserAnalyser class 30 // 31 // History: 32 // ----------- 33 // 19 Mar 2013 LP Migrated to G4tool 34 // 8 Nov 2002 GS Migration to AIDA 35 // 7 Nov 2001 MGP Implementation 36 // 37 // ------------------------------------------- 38 39 #include <fstream> 40 #include <iomanip> 41 #include "G4AutoLock.hh" 42 #include "XrayTelAnalysis.hh" 43 #include "globals.hh" 44 #include "G4SystemOfUnits.hh" 45 #include "G4Track.hh" 46 #include "G4ios.hh" 47 #include "G4SteppingManager.hh" 48 #include "G4ThreeVector.hh" 49 50 XrayTelAnalysis* XrayTelAnalysis::instance = 0 51 52 namespace { 53 //Mutex to acquire access to singleton insta 54 G4Mutex instanceMutex = G4MUTEX_INITIALIZER; 55 //Mutex to acquire accss to histograms creat 56 //It is also used to control all operations 57 //File writing and check analysis 58 G4Mutex dataManipulationMutex = G4MUTEX_INIT 59 } 60 61 XrayTelAnalysis::XrayTelAnalysis() : 62 asciiFile(0),nEnteringTracks(0), totEntering 63 { 64 histFileName = "xraytel"; 65 66 67 asciiFileName="xraytel.out"; 68 asciiFile = new std::ofstream(asciiFileName) 69 70 if(asciiFile->is_open()) 71 (*asciiFile) << "Energy (keV) x (mm) y 72 } 73 74 XrayTelAnalysis::~XrayTelAnalysis() 75 { 76 if (asciiFile) 77 delete asciiFile; 78 if (nEnteringTracks) 79 delete nEnteringTracks; 80 if (totEnteringEnergy) 81 delete totEnteringEnergy; 82 } 83 84 85 XrayTelAnalysis* XrayTelAnalysis::getInstance( 86 { 87 G4AutoLock l(&instanceMutex); 88 if (instance == 0) 89 instance = new XrayTelAnalysis; 90 return instance; 91 } 92 93 94 void XrayTelAnalysis::book(G4bool isMaster) 95 { 96 G4AutoLock l(&dataManipulationMutex); 97 98 //reset counters: do be done only once, by t 99 if (isMaster) 100 { 101 if (nEnteringTracks) 102 { 103 delete nEnteringTracks; 104 nEnteringTracks = 0; 105 } 106 nEnteringTracks = new std::map<G4int,G4i 107 108 if (totEnteringEnergy) 109 { 110 delete totEnteringEnergy; 111 totEnteringEnergy = 0; 112 } 113 totEnteringEnergy = new std::map<G4int,G 114 } 115 116 // Get/create analysis manager 117 G4AnalysisManager* man = G4AnalysisManager:: 118 man->SetDefaultFileType("root"); 119 120 // Open an output file: it is done in master 121 // printout is done only by the master, for 122 if (isMaster) 123 G4cout << "Opening output file " << histFi 124 man->OpenFile(histFileName); 125 man->SetFirstHistoId(1); 126 if (isMaster) 127 G4cout << " done" << G4endl; 128 129 // Book 1D histograms 130 man->CreateH1("h1","Energy, all /keV", 100, 131 man->CreateH1("h2","Energy, entering detecto 132 133 // Book 2D histograms (notice: the numbering 134 man->CreateH2("d1","y-z, all /mm", 100,-500. 135 man->CreateH2("d2","y-z, entering detector / 136 137 // Book ntuples 138 man->CreateNtuple("tree", "Track ntuple"); 139 man->CreateNtupleDColumn("energy"); 140 man->CreateNtupleDColumn("x"); 141 man->CreateNtupleDColumn("y"); 142 man->CreateNtupleDColumn("z"); 143 man->CreateNtupleDColumn("dirx"); 144 man->CreateNtupleDColumn("diry"); 145 man->CreateNtupleDColumn("dirz"); 146 man->FinishNtuple(); 147 } 148 149 150 void XrayTelAnalysis::finish(G4bool isMaster) 151 { 152 G4AutoLock l(&dataManipulationMutex); 153 // Save histograms 154 G4AnalysisManager* man = G4AnalysisManager:: 155 man->Write(); 156 man->CloseFile(); 157 man->Clear(); 158 159 if (!isMaster) 160 return; 161 162 //only master performs these operations 163 asciiFile->close(); 164 165 //Sequential run: output results! 166 if (nEnteringTracks->count(-2)) 167 { 168 G4cout << "End of Run summary (sequentia 169 G4cout << "Total Entering Detector : " < 170 << G4endl; 171 G4cout << "Total Entering Detector Energ 172 << (totEnteringEnergy->find(-2)->second 173 << " MeV" 174 << G4endl; 175 return; 176 } 177 178 179 //MT run: sum results 180 181 182 //MT build, but sequential run 183 if (nEnteringTracks->count(-1)) 184 { 185 G4cout << "End of Run summary (sequentia 186 G4cout << "Total Entering Detector : " < 187 << G4endl; 188 G4cout << "Total Entering Detector Energ 189 << (totEnteringEnergy->find(-1)->second 190 << " MeV" 191 << G4endl; 192 G4cout << "############################# 193 return; 194 } 195 196 G4bool loopAgain = true; 197 G4int totEntries = 0; 198 G4double totEnergy = 0.; 199 200 G4cout << "################################# 201 for (size_t i=0; loopAgain; i++) 202 { 203 //ok, this thread was found 204 if (nEnteringTracks->count(i)) 205 { 206 G4cout << "End of Run summary (thread= " < 207 G4int part = nEnteringTracks->find(i)->sec 208 G4cout << "Total Entering Detector : " << 209 G4double ene = totEnteringEnergy->find(i)- 210 G4cout << "Total Entering Detector Energy 211 << ene/MeV 212 << " MeV" << G4endl; 213 totEntries += part; 214 totEnergy += ene; 215 G4cout << "############################### 216 } 217 else 218 loopAgain = false; 219 } 220 //Report global numbers 221 if (totEntries) 222 { 223 G4cout << "End of Run summary (MT) globa 224 G4cout << "Total Entering Detector : " < 225 G4cout << "Total Entering Detector Energ 226 << totEnergy/MeV 227 << " MeV" << G4endl; 228 G4cout << "############################# 229 } 230 } 231 232 void XrayTelAnalysis::analyseStepping(const G4 233 { 234 G4AutoLock l(&dataManipulationMutex); 235 eKin = track.GetKineticEnergy()/keV; 236 G4ThreeVector pos = track.GetPosition()/mm; 237 y = pos.y(); 238 z = pos.z(); 239 G4ThreeVector dir = track.GetMomentumDirecti 240 dirX = dir.x(); 241 dirY = dir.y(); 242 dirZ = dir.z(); 243 244 // Fill histograms 245 G4AnalysisManager* man = G4AnalysisManager:: 246 man->FillH1(1,eKin); 247 man->FillH2(1,y,z); 248 249 // Fill histograms and ntuple, tracks enteri 250 if (entering) { 251 // Fill and plot histograms 252 man->FillH1(2,eKin); 253 man->FillH2(2,y,z); 254 255 man->FillNtupleDColumn(0,eKin); 256 man->FillNtupleDColumn(1,x); 257 man->FillNtupleDColumn(2,y); 258 man->FillNtupleDColumn(3,z); 259 man->FillNtupleDColumn(4,dirX); 260 man->FillNtupleDColumn(5,dirY); 261 man->FillNtupleDColumn(6,dirZ); 262 man->AddNtupleRow(); 263 } 264 265 266 // Write to file 267 if (entering) { 268 if(asciiFile->is_open()) { 269 (*asciiFile) << std::setiosflags(std::io 270 << std::setprecision(3) 271 << std::setiosflags(std::ios::right) 272 << std::setw(10); 273 (*asciiFile) << eKin; 274 (*asciiFile) << std::setiosflags(std::io 275 << std::setprecision(3) 276 << std::setiosflags(std::ios::right) 277 << std::setw(10); 278 (*asciiFile) << x; 279 (*asciiFile) << std::setiosflags(std::io 280 << std::setprecision(3) 281 << std::setiosflags(std::ios::right) 282 << std::setw(10); 283 (*asciiFile) << y; 284 (*asciiFile) << std::setiosflags(std::io 285 << std::setprecision(3) 286 << std::setiosflags(std::ios::right) 287 << std::setw(10); 288 (*asciiFile) << z 289 << G4endl; 290 } 291 } 292 } 293 294 void XrayTelAnalysis::Update(G4double energy,G 295 { 296 G4AutoLock l(&dataManipulationMutex); 297 //It already exists: increase the counter 298 if (nEnteringTracks->count(threadID)) 299 { 300 (nEnteringTracks->find(threadID)->second 301 } 302 else //enter a new one 303 { 304 G4int tracks = 1; 305 nEnteringTracks->insert(std::make_pair(t 306 } 307 308 //It already exists: increase the counter 309 if (totEnteringEnergy->count(threadID)) 310 (totEnteringEnergy->find(threadID)->second 311 else //enter a new one 312 totEnteringEnergy->insert(std::make_pair(t 313 314 return; 315 } 316 317