Geant4  10.01.p01
XrayTelAnalysis.cc
Go to the documentation of this file.
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 // $Id: XrayTelAnalysis.cc 81195 2014-05-22 14:49:13Z gcosmo $
28 //
29 // Author: A. Pfeiffer (Andreas.Pfeiffer@cern.ch)
30 // (copied from his UserAnalyser class)
31 //
32 // History:
33 // -----------
34 // 19 Mar 2013 LP Migrated to G4tools
35 // 8 Nov 2002 GS Migration to AIDA 3
36 // 7 Nov 2001 MGP Implementation
37 //
38 // -------------------------------------------------------------------
39 
40 #include <fstream>
41 #include <iomanip>
42 #include "G4AutoLock.hh"
43 #include "XrayTelAnalysis.hh"
44 #include "globals.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4Track.hh"
47 #include "G4ios.hh"
48 #include "G4SteppingManager.hh"
49 #include "G4ThreeVector.hh"
50 
52 
53 namespace {
54  //Mutex to acquire access to singleton instance check/creation
55  G4Mutex instanceMutex = G4MUTEX_INITIALIZER;
56  //Mutex to acquire accss to histograms creation/access
57  //It is also used to control all operations related to histos
58  //File writing and check analysis
59  G4Mutex dataManipulationMutex = G4MUTEX_INITIALIZER;
60 }
61 
63  asciiFile(0),nEnteringTracks(0), totEnteringEnergy(0)
64 {
65  histFileName = "xraytel";
66 
67 
68  asciiFileName="xraytel.out";
69  asciiFile = new std::ofstream(asciiFileName);
70 
71  if(asciiFile->is_open())
72  (*asciiFile) << "Energy (keV) x (mm) y (mm) z (mm)" << G4endl << G4endl;
73 }
74 
76 {
77  if (asciiFile)
78  delete asciiFile;
79  if (nEnteringTracks)
80  delete nEnteringTracks;
82  delete totEnteringEnergy;
83 }
84 
85 
87 {
88  G4AutoLock l(&instanceMutex);
89  if (instance == 0)
91  return instance;
92 }
93 
94 
96 {
97  G4AutoLock l(&dataManipulationMutex);
98 
99  //reset counters: do be done only once, by the master
100  if (isMaster)
101  {
102  if (nEnteringTracks)
103  {
104  delete nEnteringTracks;
105  nEnteringTracks = 0;
106  }
107  nEnteringTracks = new std::map<G4int,G4int>;
108 
109  if (totEnteringEnergy)
110  {
111  delete totEnteringEnergy;
112  totEnteringEnergy = 0;
113  }
114  totEnteringEnergy = new std::map<G4int,G4double>;
115  }
116 
117  // Get/create analysis manager
118  G4AnalysisManager* man = G4AnalysisManager::Instance();
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 
151 {
152  G4AutoLock l(&dataManipulationMutex);
153  // Save histograms
154  G4AnalysisManager* man = G4AnalysisManager::Instance();
155  man->Write();
156  man->CloseFile();
157  // Complete clean-up
158  delete G4AnalysisManager::Instance();
159 
160  if (!isMaster)
161  return;
162 
163  //only master performs these operations
164  asciiFile->close();
165 
166  //Sequential run: output results!
167  if (nEnteringTracks->count(-2))
168  {
169  G4cout << "End of Run summary (sequential)" << G4endl << G4endl;
170  G4cout << "Total Entering Detector : " << nEnteringTracks->find(-2)->second
171  << G4endl;
172  G4cout << "Total Entering Detector Energy : "
173  << (totEnteringEnergy->find(-2)->second)/MeV
174  << " MeV"
175  << G4endl;
176  return;
177  }
178 
179 
180  //MT run: sum results
181 
182 
183  //MT build, but sequential run
184  if (nEnteringTracks->count(-1))
185  {
186  G4cout << "End of Run summary (sequential with MT build)" << G4endl << G4endl;
187  G4cout << "Total Entering Detector : " << nEnteringTracks->find(-1)->second
188  << G4endl;
189  G4cout << "Total Entering Detector Energy : "
190  << (totEnteringEnergy->find(-1)->second)/MeV
191  << " MeV"
192  << G4endl;
193  G4cout << "##########################################" << G4endl;
194  return;
195  }
196 
197  G4bool loopAgain = true;
198  G4int totEntries = 0;
199  G4double totEnergy = 0.;
200 
201  G4cout << "##########################################" << G4endl;
202  for (size_t i=0; loopAgain; i++)
203  {
204  //ok, this thread was found
205  if (nEnteringTracks->count(i))
206  {
207  G4cout << "End of Run summary (thread= " << i << ")" << G4endl;
208  G4int part = nEnteringTracks->find(i)->second;
209  G4cout << "Total Entering Detector : " << part << G4endl;
210  G4double ene = totEnteringEnergy->find(i)->second;
211  G4cout << "Total Entering Detector Energy : "
212  << ene/MeV
213  << " MeV" << G4endl;
214  totEntries += part;
215  totEnergy += ene;
216  G4cout << "##########################################" << G4endl;
217  }
218  else
219  loopAgain = false;
220  }
221  //Report global numbers
222  if (totEntries)
223  {
224  G4cout << "End of Run summary (MT) global" << G4endl << G4endl;
225  G4cout << "Total Entering Detector : " << totEntries << G4endl;
226  G4cout << "Total Entering Detector Energy : "
227  << totEnergy/MeV
228  << " MeV" << G4endl;
229  G4cout << "##########################################" << G4endl;
230  }
231 }
232 
233 void XrayTelAnalysis::analyseStepping(const G4Track& track, G4bool entering)
234 {
235  G4AutoLock l(&dataManipulationMutex);
236  eKin = track.GetKineticEnergy()/keV;
237  G4ThreeVector pos = track.GetPosition()/mm;
238  y = pos.y();
239  z = pos.z();
240  G4ThreeVector dir = track.GetMomentumDirection();
241  dirX = dir.x();
242  dirY = dir.y();
243  dirZ = dir.z();
244 
245  // Fill histograms
246  G4AnalysisManager* man = G4AnalysisManager::Instance();
247  man->FillH1(1,eKin);
248  man->FillH2(1,y,z);
249 
250  // Fill histograms and ntuple, tracks entering the detector
251  if (entering) {
252  // Fill and plot histograms
253  man->FillH1(2,eKin);
254  man->FillH2(2,y,z);
255 
256  man->FillNtupleDColumn(0,eKin);
257  man->FillNtupleDColumn(1,x);
258  man->FillNtupleDColumn(2,y);
259  man->FillNtupleDColumn(3,z);
260  man->FillNtupleDColumn(4,dirX);
261  man->FillNtupleDColumn(5,dirY);
262  man->FillNtupleDColumn(6,dirZ);
263  man->AddNtupleRow();
264  }
265 
266 
267  // Write to file
268  if (entering) {
269  if(asciiFile->is_open()) {
270  (*asciiFile) << std::setiosflags(std::ios::fixed)
271  << std::setprecision(3)
272  << std::setiosflags(std::ios::right)
273  << std::setw(10);
274  (*asciiFile) << eKin;
275  (*asciiFile) << std::setiosflags(std::ios::fixed)
276  << std::setprecision(3)
277  << std::setiosflags(std::ios::right)
278  << std::setw(10);
279  (*asciiFile) << x;
280  (*asciiFile) << std::setiosflags(std::ios::fixed)
281  << std::setprecision(3)
282  << std::setiosflags(std::ios::right)
283  << std::setw(10);
284  (*asciiFile) << y;
285  (*asciiFile) << std::setiosflags(std::ios::fixed)
286  << std::setprecision(3)
287  << std::setiosflags(std::ios::right)
288  << std::setw(10);
289  (*asciiFile) << z
290  << G4endl;
291  }
292  }
293 }
294 
296 {
297  G4AutoLock l(&dataManipulationMutex);
298  //It already exists: increase the counter
299  if (nEnteringTracks->count(threadID))
300  {
301  (nEnteringTracks->find(threadID)->second)++;
302  }
303  else //enter a new one
304  {
305  G4int tracks = 1;
306  nEnteringTracks->insert(std::make_pair(threadID,tracks));
307  }
308 
309  //It already exists: increase the counter
310  if (totEnteringEnergy->count(threadID))
311  (totEnteringEnergy->find(threadID)->second) += energy;
312  else //enter a new one
313  totEnteringEnergy->insert(std::make_pair(threadID,energy));
314 
315  return;
316 }
317 
static const double MeV
Definition: G4SIunits.hh:193
void book(G4bool isMaster)
G4bool SetFirstHistoId(G4int firstId)
void finish(G4bool isMaster)
G4int CreateH1(const G4String &name, const G4String &title, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binSchemeName="linear")
CLHEP::Hep3Vector G4ThreeVector
const G4ThreeVector & GetPosition() const
G4int CreateNtuple(const G4String &name, const G4String &title)
int G4int
Definition: G4Types.hh:78
G4bool OpenFile(const G4String &fileName="")
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:171
G4String asciiFileName
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4bool FillNtupleDColumn(G4int id, G4double value)
bool G4bool
Definition: G4Types.hh:79
void Update(G4double energy, G4int threadID)
G4int G4Mutex
Definition: G4Threading.hh:169
const G4ThreeVector & GetMomentumDirection() const
G4bool FillH2(G4int id, G4double xvalue, G4double yvalue, G4double weight=1.0)
G4double energy(const ThreeVector &p, const G4double m)
G4bool FillH1(G4int id, G4double value, G4double weight=1.0)
#define G4endl
Definition: G4ios.hh:61
std::map< G4int, G4double > * totEnteringEnergy
static const double keV
Definition: G4SIunits.hh:195
void analyseStepping(const G4Track &track, G4bool entering)
G4int CreateNtupleDColumn(const G4String &name)
double G4double
Definition: G4Types.hh:76
static const double mm
Definition: G4SIunits.hh:102
std::map< G4int, G4int > * nEnteringTracks
std::ofstream * asciiFile
static const G4double pos
static XrayTelAnalysis * getInstance()
G4int CreateH2(const G4String &name, const G4String &title, G4int nxbins, G4double xmin, G4double xmax, G4int nybins, G4double ymin, G4double ymax, const G4String &xunitName="none", const G4String &yunitName="none", const G4String &xfcnName="none", const G4String &yfcnName="none", const G4String &xbinSchemeName="linear", const G4String &ybinSchemeName="linear")
static XrayTelAnalysis * instance