xray_telescope/src/XrayTelAnalysis.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // Author:  A. Pfeiffer (Andreas.Pfeiffer@cern.ch)
0029 //         (copied from his UserAnalyser class)
0030 //
0031 // History:
0032 // -----------
0033 // 19 Mar 2013   LP         Migrated to G4tools
0034 //  8 Nov 2002   GS         Migration to AIDA 3
0035 //  7 Nov 2001   MGP        Implementation
0036 //
0037 // -------------------------------------------------------------------
0038
0039 #include <fstream>
0040 #include <iomanip>
0041 #include "G4AutoLock.hh"
0042 #include "XrayTelAnalysis.hh"
0043 #include "globals.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "G4Track.hh"
0046 #include "G4ios.hh"
0047 #include "G4SteppingManager.hh"
0048 #include "G4ThreeVector.hh"
0049
0050 XrayTelAnalysis* XrayTelAnalysis::instance = 0;
0051
0052 namespace {
0053   //Mutex to acquire access to singleton instance check/creation
0054   G4Mutex instanceMutex = G4MUTEX_INITIALIZER;
0055   //Mutex to acquire accss to histograms creation/access
0056   //It is also used to control all operations related to histos
0057   //File writing and check analysis
0058   G4Mutex dataManipulationMutex = G4MUTEX_INITIALIZER;
0059 }
0060
0061 XrayTelAnalysis::XrayTelAnalysis() :
0062   asciiFile(0),nEnteringTracks(0), totEnteringEnergy(0)
0063 {
0064   histFileName = "xraytel";
0065
0066
0067   asciiFileName="xraytel.out";
0068   asciiFile = new std::ofstream(asciiFileName);
0069
0070   if(asciiFile->is_open())
0071     (*asciiFile) << "Energy (keV)  x (mm)    y (mm)    z (mm)" << G4endl << G4endl;
0072 }
0073
0074 XrayTelAnalysis::~XrayTelAnalysis()
0075 {
0076   if (asciiFile)
0077     delete asciiFile;
0078   if (nEnteringTracks)
0079     delete nEnteringTracks;
0080   if (totEnteringEnergy)
0081     delete totEnteringEnergy;
0082 }
0083
0084
0085 XrayTelAnalysis* XrayTelAnalysis::getInstance()
0086 {
0087   G4AutoLock l(&instanceMutex);
0088   if (instance == 0)
0089     instance = new XrayTelAnalysis;
0090   return instance;
0091 }
0092
0093
0094 void XrayTelAnalysis::book(G4bool isMaster)
0095 {
0096   G4AutoLock l(&dataManipulationMutex);
0097
0098   //reset counters: do be done only once, by the master
0099   if (isMaster)
0100     {
0101       if (nEnteringTracks)
0102     {
0103       delete nEnteringTracks;
0104       nEnteringTracks = 0;
0105     }
0106       nEnteringTracks = new std::map<G4int,G4int>;
0107
0108       if (totEnteringEnergy)
0109     {
0110       delete totEnteringEnergy;
0111       totEnteringEnergy = 0;
0112     }
0113       totEnteringEnergy = new std::map<G4int,G4double>;
0114     }
0115
0116   // Get/create analysis manager
0117   G4AnalysisManager* man = G4AnalysisManager::Instance();
0118   man->SetDefaultFileType("root");
0119
0120   // Open an output file: it is done in master and threads. The
0121   // printout is done only by the master, for tidyness
0122   if (isMaster)
0123     G4cout << "Opening output file " << histFileName << " ... ";
0124   man->OpenFile(histFileName);
0125   man->SetFirstHistoId(1);
0126   if (isMaster)
0127     G4cout << " done" << G4endl;
0128
0129   // Book 1D histograms
0130   man->CreateH1("h1","Energy, all /keV",  100,0.,100.);
0131   man->CreateH1("h2","Energy, entering detector /keV", 500,0.,500.);
0132
0133   // Book 2D histograms (notice: the numbering is independent)
0134   man->CreateH2("d1","y-z, all /mm", 100,-500.,500.,100,-500.,500.);
0135   man->CreateH2("d2","y-z, entering detector /mm", 200,-50.,50.,200,-50.,50.);
0136
0137   // Book ntuples
0138   man->CreateNtuple("tree", "Track ntuple");
0139   man->CreateNtupleDColumn("energy");
0140   man->CreateNtupleDColumn("x");
0141   man->CreateNtupleDColumn("y");
0142   man->CreateNtupleDColumn("z");
0143   man->CreateNtupleDColumn("dirx");
0144   man->CreateNtupleDColumn("diry");
0145   man->CreateNtupleDColumn("dirz");
0146   man->FinishNtuple();
0147 }
0148
0149
0150 void XrayTelAnalysis::finish(G4bool isMaster)
0151 {
0152   G4AutoLock l(&dataManipulationMutex);
0153   // Save histograms
0154   G4AnalysisManager* man = G4AnalysisManager::Instance();
0155   man->Write();
0156   man->CloseFile();
0157   man->Clear();
0158
0159   if (!isMaster)
0160     return;
0161
0162   //only master performs these operations
0163   asciiFile->close();
0164
0165   //Sequential run: output results!
0166   if (nEnteringTracks->count(-2))
0167     {
0168       G4cout << "End of Run summary (sequential)" << G4endl << G4endl;
0169       G4cout << "Total Entering Detector : " << nEnteringTracks->find(-2)->second
0170          << G4endl;
0171       G4cout << "Total Entering Detector Energy : "
0172          << (totEnteringEnergy->find(-2)->second)/MeV
0173          << " MeV"
0174          << G4endl;
0175       return;
0176     }
0177
0178
0179   //MT run: sum results
0180
0181
0182   //MT build, but sequential run
0183   if (nEnteringTracks->count(-1))
0184     {
0185       G4cout << "End of Run summary (sequential with MT build)" << G4endl << G4endl;
0186       G4cout << "Total Entering Detector : " << nEnteringTracks->find(-1)->second
0187          << G4endl;
0188       G4cout << "Total Entering Detector Energy : "
0189          << (totEnteringEnergy->find(-1)->second)/MeV
0190          << " MeV"
0191          << G4endl;
0192       G4cout << "##########################################" << G4endl;
0193       return;
0194     }
0195
0196   G4bool loopAgain = true;
0197   G4int totEntries = 0;
0198   G4double totEnergy = 0.;
0199
0200   G4cout << "##########################################" << G4endl;
0201   for (size_t i=0; loopAgain; i++)
0202     {
0203       //ok, this thread was found
0204       if (nEnteringTracks->count(i))
0205     {
0206       G4cout << "End of Run summary (thread= " << i << ")" << G4endl;
0207       G4int part = nEnteringTracks->find(i)->second;
0208       G4cout << "Total Entering Detector : " << part << G4endl;
0209       G4double ene = totEnteringEnergy->find(i)->second;
0210       G4cout << "Total Entering Detector Energy : "
0211          << ene/MeV
0212          << " MeV" << G4endl;
0213       totEntries += part;
0214       totEnergy += ene;
0215       G4cout << "##########################################" << G4endl;
0216     }
0217       else
0218     loopAgain = false;
0219     }
0220   //Report global numbers
0221   if (totEntries)
0222     {
0223       G4cout << "End of Run summary (MT) global" << G4endl << G4endl;
0224       G4cout << "Total Entering Detector : " << totEntries << G4endl;
0225       G4cout << "Total Entering Detector Energy : "
0226          << totEnergy/MeV
0227          << " MeV" << G4endl;
0228       G4cout << "##########################################" << G4endl;
0229     }
0230 }
0231
0232 void XrayTelAnalysis::analyseStepping(const G4Track& track, G4bool entering)
0233 {
0234   G4AutoLock l(&dataManipulationMutex);
0235   eKin = track.GetKineticEnergy()/keV;
0236   G4ThreeVector pos = track.GetPosition()/mm;
0237   y = pos.y();
0238   z = pos.z();
0239   G4ThreeVector dir = track.GetMomentumDirection();
0240   dirX = dir.x();
0241   dirY = dir.y();
0242   dirZ = dir.z();
0243
0244   // Fill histograms
0245   G4AnalysisManager* man = G4AnalysisManager::Instance();
0246   man->FillH1(1,eKin);
0247   man->FillH2(1,y,z);
0248
0249   // Fill histograms and ntuple, tracks entering the detector
0250   if (entering) {
0251     // Fill and plot histograms
0252     man->FillH1(2,eKin);
0253     man->FillH2(2,y,z);
0254
0255     man->FillNtupleDColumn(0,eKin);
0256     man->FillNtupleDColumn(1,x);
0257     man->FillNtupleDColumn(2,y);
0258     man->FillNtupleDColumn(3,z);
0259     man->FillNtupleDColumn(4,dirX);
0260     man->FillNtupleDColumn(5,dirY);
0261     man->FillNtupleDColumn(6,dirZ);
0262     man->AddNtupleRow();
0263   }
0264
0265
0266   // Write to file
0267   if (entering) {
0268     if(asciiFile->is_open()) {
0269       (*asciiFile) << std::setiosflags(std::ios::fixed)
0270            << std::setprecision(3)
0271            << std::setiosflags(std::ios::right)
0272            << std::setw(10);
0273       (*asciiFile) << eKin;
0274       (*asciiFile) << std::setiosflags(std::ios::fixed)
0275            << std::setprecision(3)
0276            << std::setiosflags(std::ios::right)
0277            << std::setw(10);
0278       (*asciiFile) << x;
0279       (*asciiFile) << std::setiosflags(std::ios::fixed)
0280            << std::setprecision(3)
0281            << std::setiosflags(std::ios::right)
0282            << std::setw(10);
0283       (*asciiFile) << y;
0284       (*asciiFile) << std::setiosflags(std::ios::fixed)
0285            << std::setprecision(3)
0286            << std::setiosflags(std::ios::right)
0287            << std::setw(10);
0288       (*asciiFile) << z
0289            << G4endl;
0290     }
0291   }
0292 }
0293
0294 void XrayTelAnalysis::Update(G4double energy,G4int threadID)
0295 {
0296   G4AutoLock l(&dataManipulationMutex);
0297   //It already exists: increase the counter
0298   if (nEnteringTracks->count(threadID))
0299     {
0300       (nEnteringTracks->find(threadID)->second)++;
0301     }
0302   else //enter a new one
0303     {
0304       G4int tracks = 1;
0305       nEnteringTracks->insert(std::make_pair(threadID,tracks));
0306     }
0307
0308   //It already exists: increase the counter
0309   if (totEnteringEnergy->count(threadID))
0310     (totEnteringEnergy->find(threadID)->second) += energy;
0311   else //enter a new one
0312     totEnteringEnergy->insert(std::make_pair(threadID,energy));
0313
0314   return;
0315 }
0316