ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
XrayTelAnalysis.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file XrayTelAnalysis.cc
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 
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 
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 
75 {
76  if (asciiFile)
77  delete asciiFile;
78  if (nEnteringTracks)
79  delete nEnteringTracks;
81  delete totEnteringEnergy;
82 }
83 
84 
86 {
87  G4AutoLock l(&instanceMutex);
88  if (instance == 0)
90  return instance;
91 }
92 
93 
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
118 
119  // Open an output file: it is done in master and threads. The
120  // printout is done only by the master, for tidyness
121  if (isMaster)
122  G4cout << "Opening output file " << histFileName << " ... ";
123  man->OpenFile(histFileName);
124  man->SetFirstHistoId(1);
125  if (isMaster)
126  G4cout << " done" << G4endl;
127 
128  // Book 1D histograms
129  man->CreateH1("h1","Energy, all /keV", 100,0.,100.);
130  man->CreateH1("h2","Energy, entering detector /keV", 500,0.,500.);
131 
132  // Book 2D histograms (notice: the numbering is independent)
133  man->CreateH2("d1","y-z, all /mm", 100,-500.,500.,100,-500.,500.);
134  man->CreateH2("d2","y-z, entering detector /mm", 200,-50.,50.,200,-50.,50.);
135 
136  // Book ntuples
137  man->CreateNtuple("tree", "Track ntuple");
138  man->CreateNtupleDColumn("energy");
139  man->CreateNtupleDColumn("x");
140  man->CreateNtupleDColumn("y");
141  man->CreateNtupleDColumn("z");
142  man->CreateNtupleDColumn("dirx");
143  man->CreateNtupleDColumn("diry");
144  man->CreateNtupleDColumn("dirz");
145  man->FinishNtuple();
146 }
147 
148 
150 {
151  G4AutoLock l(&dataManipulationMutex);
152  // Save histograms
154  man->Write();
155  man->CloseFile();
156  // Complete clean-up
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 
233 {
234  G4AutoLock l(&dataManipulationMutex);
235  eKin = track.GetKineticEnergy()/keV;
236  G4ThreeVector pos = track.GetPosition()/mm;
237  y = pos.y();
238  z = pos.z();
240  dirX = dir.x();
241  dirY = dir.y();
242  dirZ = dir.z();
243 
244  // Fill histograms
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 
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