exMPI04/src/Analysis.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 /// @file Analysis.cc
0028 /// @brief Define histograms and ntuples
0029
0030 #include "Analysis.hh"
0031
0032 #include "G4AutoDelete.hh"
0033 #include "G4SystemOfUnits.hh"
0034
0035 // Note: ntuple merging is supported only with Root format
0036 #ifndef G4MULTITHREADED
0037 #  include "G4RootMpiAnalysisManager.hh"
0038 using G4AnalysisManager = G4RootMpiAnalysisManager;
0039 #else
0040 #  include "G4RootAnalysisManager.hh"
0041 using G4AnalysisManager = G4RootAnalysisManager;
0042 #endif
0043
0044 G4ThreadLocal G4int Analysis::fincidentFlag = false;
0045 G4ThreadLocal Analysis* the_analysis = 0;
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048 Analysis* Analysis::GetAnalysis()
0049 {
0050   if (!the_analysis) {
0051     the_analysis = new Analysis();
0052     G4AutoDelete::Register(the_analysis);
0053   }
0054   return the_analysis;
0055 }
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0058 Analysis::Analysis()
0059   : fUseNtuple(true),
0060     fMergeNtuple(true),
0061     fincident_x_hist(0),
0062     fincident_map(0),
0063     fdose_hist(0),
0064     fdose_map(0),
0065     fdose_prof(0),
0066     fdose_map_prof(0),
0067     fdose_map3d(0)
0068 {}
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071 void Analysis::Book()
0072 {
0073   G4cout << "Analysis::Book start, fUseNtuple: " << fUseNtuple << G4endl;
0074
0075   G4AnalysisManager* mgr = G4AnalysisManager::Instance();
0076   mgr->SetVerboseLevel(1);
0077 #ifdef G4MULTITHREADED
0078   // MT ntuple merging
0079   mgr->SetNtupleMerging(fMergeNtuple);
0080 #endif
0081
0082   fincident_x_hist = mgr->CreateH1("incident_x", "Incident X", 100, -5 * cm, 5 * cm, "cm");
0083   fincident_map = mgr->CreateH2("incident_map", "Incident Map", 50, -5 * cm, 5 * cm, 50, -5 * cm,
0084                                 5 * cm, "cm", "cm");
0085   fdose_hist = mgr->CreateH1("dose", "Dose distribution", 500, 0, 50 * cm, "cm");
0086   fdose_map = mgr->CreateH2("dose_map", "Dose distribution", 500, 0, 50 * cm, 200, -10 * cm,
0087                             10 * cm, "cm", "cm");
0088   fdose_map3d = mgr->CreateH3("dose_map_3d", "Dose distribution", 30, 0, 50 * cm, 20, -10 * cm,
0089                               10 * cm, 20, -10 * cm, 10 * cm, "cm", "cm", "cm");
0090   fdose_prof =
0091     mgr->CreateP1("dose_prof", "Dose distribution", 300, 0, 30 * cm, 0, 100 * MeV, "cm", "MeV");
0092   fdose_map_prof = mgr->CreateP2("dose_map_prof", "Dose distribution", 300, 0, 30 * cm, 80, -4 * cm,
0093                                  4 * cm, 0, 100 * MeV, "cm", "cm", "MeV");
0094
0095   if (fUseNtuple) {
0096     mgr->CreateNtuple("Dose", "Dose distribution");  // ntuple Id = 0
0097     mgr->CreateNtupleDColumn("pz");
0098     mgr->CreateNtupleDColumn("px");
0099     mgr->CreateNtupleDColumn("dose");
0100     mgr->FinishNtuple();
0101   }
0102
0103   G4cout << "Analysis::Book finished " << G4endl;
0104 }
0105
0106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0107 Analysis::~Analysis() {}
0108
0109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0110 void Analysis::OpenFile(const G4String& fname)
0111 {
0112   G4AnalysisManager* mgr = G4AnalysisManager::Instance();
0113   mgr->OpenFile(fname.c_str());
0114 }
0115
0116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0117 void Analysis::Save()
0118 {
0119   G4AnalysisManager* mgr = G4AnalysisManager::Instance();
0120   mgr->Write();
0121 }
0122
0123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0124 void Analysis::Close(G4bool reset)
0125 {
0126   G4AnalysisManager* mgr = G4AnalysisManager::Instance();
0127   mgr->CloseFile(reset);
0128 }
0129
0130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0131 void Analysis::FillIncident(const G4ThreeVector& p)
0132 {
0133   if (!fincidentFlag) {
0134     G4AnalysisManager* mgr = G4AnalysisManager::Instance();
0135     mgr->FillH2(fincident_map, p.x(), p.y());
0136     mgr->FillH1(fincident_x_hist, p.x());
0137     fincidentFlag = true;
0138   }
0139 }
0140
0141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0142 void Analysis::FillDose(const G4ThreeVector& p, G4double dedx)
0143 {
0144   const G4double dxy = 10. * mm;
0145   if (std::abs(p.y()) < dxy) {
0146     G4AnalysisManager* mgr = G4AnalysisManager::Instance();
0147     const G4double Z0 = 25. * cm;
0148
0149     mgr->FillH2(fdose_map, p.z() + Z0, p.x(), dedx / GeV);
0150     mgr->FillP2(fdose_map_prof, p.z() + Z0, p.x(), dedx);
0151     mgr->FillH3(fdose_map3d, p.z() + Z0, p.x(), p.y(), dedx / GeV);
0152     if (std::abs(p.x()) < dxy) {
0153       mgr->FillH1(fdose_hist, p.z() + Z0, dedx / GeV);
0154       mgr->FillP1(fdose_prof, p.z() + Z0, dedx);
0155     }
0156
0157     if (fUseNtuple) {
0158       // the same fill frequency as H2 "dose_map"
0159       mgr->FillNtupleDColumn(0, p.z() + Z0);
0160       mgr->FillNtupleDColumn(1, p.x());
0161       mgr->FillNtupleDColumn(2, dedx / GeV);
0162       mgr->AddNtupleRow();
0163     }
0164   }
0165 }