Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/xray_telescope/src/XrayTelAnalysis.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 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 // Author:  A. Pfeiffer (Andreas.Pfeiffer@cern.ch) 
 29 //         (copied from his UserAnalyser class)
 30 //
 31 // History:
 32 // -----------
 33 // 19 Mar 2013   LP         Migrated to G4tools
 34 //  8 Nov 2002   GS         Migration to AIDA 3
 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 instance check/creation
 54   G4Mutex instanceMutex = G4MUTEX_INITIALIZER;
 55   //Mutex to acquire accss to histograms creation/access
 56   //It is also used to control all operations related to histos 
 57   //File writing and check analysis
 58   G4Mutex dataManipulationMutex = G4MUTEX_INITIALIZER;
 59 }
 60 
 61 XrayTelAnalysis::XrayTelAnalysis() : 
 62   asciiFile(0),nEnteringTracks(0), totEnteringEnergy(0)
 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 (mm)    z (mm)" << G4endl << G4endl;  
 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 the master
 99   if (isMaster)
100     {
101       if (nEnteringTracks)    
102   {
103     delete nEnteringTracks;
104     nEnteringTracks = 0;
105   }
106       nEnteringTracks = new std::map<G4int,G4int>;
107   
108       if (totEnteringEnergy)
109   {
110     delete totEnteringEnergy;
111     totEnteringEnergy = 0;
112   }
113       totEnteringEnergy = new std::map<G4int,G4double>;
114     }
115 
116   // Get/create analysis manager
117   G4AnalysisManager* man = G4AnalysisManager::Instance();
118   man->SetDefaultFileType("root");
119 
120   // Open an output file: it is done in master and threads. The 
121   // printout is done only by the master, for tidyness
122   if (isMaster)
123     G4cout << "Opening output file " << histFileName << " ... ";
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,0.,100.);
131   man->CreateH1("h2","Energy, entering detector /keV", 500,0.,500.);
132 
133   // Book 2D histograms (notice: the numbering is independent)
134   man->CreateH2("d1","y-z, all /mm", 100,-500.,500.,100,-500.,500.); 
135   man->CreateH2("d2","y-z, entering detector /mm", 200,-50.,50.,200,-50.,50.);
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::Instance();
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 (sequential)" << G4endl << G4endl;
169       G4cout << "Total Entering Detector : " << nEnteringTracks->find(-2)->second  
170        << G4endl;
171       G4cout << "Total Entering Detector Energy : " 
172        << (totEnteringEnergy->find(-2)->second)/MeV  
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 (sequential with MT build)" << G4endl << G4endl;
186       G4cout << "Total Entering Detector : " << nEnteringTracks->find(-1)->second  
187        << G4endl;
188       G4cout << "Total Entering Detector Energy : " 
189        << (totEnteringEnergy->find(-1)->second)/MeV  
190        << " MeV"
191        << G4endl;
192       G4cout << "##########################################" << G4endl;
193       return;
194     }
195 
196   G4bool loopAgain = true;
197   G4int totEntries = 0;
198   G4double totEnergy = 0.;
199 
200   G4cout << "##########################################" << G4endl;
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= " << i << ")" << G4endl;
207     G4int part = nEnteringTracks->find(i)->second;
208     G4cout << "Total Entering Detector : " << part << G4endl;
209     G4double ene = totEnteringEnergy->find(i)->second;
210     G4cout << "Total Entering Detector Energy : " 
211      << ene/MeV  
212      << " MeV" << G4endl;
213     totEntries += part;
214     totEnergy += ene;
215     G4cout << "##########################################" << G4endl;
216   }
217       else
218   loopAgain = false;
219     }
220   //Report global numbers
221   if (totEntries)
222     {
223       G4cout << "End of Run summary (MT) global" << G4endl << G4endl;
224       G4cout << "Total Entering Detector : " << totEntries << G4endl;
225       G4cout << "Total Entering Detector Energy : " 
226        << totEnergy/MeV  
227        << " MeV" << G4endl;
228       G4cout << "##########################################" << G4endl;
229     }
230 }
231 
232 void XrayTelAnalysis::analyseStepping(const G4Track& track, G4bool entering)
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.GetMomentumDirection();
240   dirX = dir.x();
241   dirY = dir.y();
242   dirZ = dir.z();
243 
244   // Fill histograms
245   G4AnalysisManager* man = G4AnalysisManager::Instance();
246   man->FillH1(1,eKin);
247   man->FillH2(1,y,z);
248   
249   // Fill histograms and ntuple, tracks entering the detector
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::ios::fixed)
270        << std::setprecision(3)
271        << std::setiosflags(std::ios::right)
272        << std::setw(10);
273       (*asciiFile) << eKin;
274       (*asciiFile) << std::setiosflags(std::ios::fixed)
275        << std::setprecision(3)
276        << std::setiosflags(std::ios::right)
277        << std::setw(10);
278       (*asciiFile) << x;
279       (*asciiFile) << std::setiosflags(std::ios::fixed)
280        << std::setprecision(3)
281        << std::setiosflags(std::ios::right)
282        << std::setw(10);
283       (*asciiFile) << y;
284       (*asciiFile) << std::setiosflags(std::ios::fixed)
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,G4int threadID)
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(threadID,tracks));
306     }
307 
308   //It already exists: increase the counter
309   if (totEnteringEnergy->count(threadID))
310     (totEnteringEnergy->find(threadID)->second) += energy;
311   else //enter a new one    
312     totEnteringEnergy->insert(std::make_pair(threadID,energy));
313     
314   return;
315 }
316 
317