Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 /// \file electromagnetic/TestEm6/src/RunActio 26 /// \file electromagnetic/TestEm6/src/RunAction.cc 27 /// \brief Implementation of the RunAction cla 27 /// \brief Implementation of the RunAction class 28 // 28 // 29 // << 29 // $Id$ >> 30 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 33 33 #include "RunAction.hh" 34 #include "RunAction.hh" 34 35 35 #include "G4AnnihiToMuPair.hh" << 36 #include "G4EmCalculator.hh" << 37 #include "G4MuonMinus.hh" << 38 #include "G4ParticleDefinition.hh" << 39 #include "G4ParticleTable.hh" << 40 #include "G4PhysicalConstants.hh" << 41 #include "G4Positron.hh" << 42 #include "G4Run.hh" 36 #include "G4Run.hh" 43 #include "G4RunManager.hh" 37 #include "G4RunManager.hh" 44 #include "G4SystemOfUnits.hh" << 45 #include "G4eBremsstrahlung.hh" << 46 #include "G4eeToHadrons.hh" << 47 #include "G4eeToHadronsModel.hh" << 48 #include "Randomize.hh" << 49 38 50 #include <sstream> << 39 #include "Randomize.hh" 51 40 52 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 53 42 54 RunAction::RunAction(DetectorConstruction* det << 43 RunAction::RunAction() 55 : G4UserRunAction(), fDetector(det), fProcCo << 56 { 44 { 57 fMinE = 40 * GeV; << 45 // Create analysis manager 58 fMaxE = 10000 * GeV; << 46 // The choice of analysis technology is done via selection of a namespace 59 fnBin = 10000; << 47 // >> 48 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); >> 49 >> 50 // Open an output file >> 51 // >> 52 G4String fileName = "testem6"; >> 53 analysisManager->OpenFile(fileName); >> 54 analysisManager->SetVerboseLevel(1); >> 55 G4String extension = analysisManager->GetFileType(); >> 56 fileName = fileName + "." + extension; >> 57 >> 58 // Creating histograms >> 59 // >> 60 analysisManager->SetFirstHistoId(1); >> 61 analysisManager->CreateH1("1","1/(1+(theta+[g]+)**2)",100, 0 ,1.); >> 62 analysisManager->CreateH1("2","log10(theta+ [g]+)", 100,-3.,1.); >> 63 analysisManager->CreateH1("3","log10(theta- [g]-)", 100,-3.,1.); >> 64 analysisManager->CreateH1("4","log10(theta+ [g]+ -theta- [g]-)", 100,-3.,1.); >> 65 analysisManager->CreateH1("5","xPlus" ,100,0.,1.); >> 66 analysisManager->CreateH1("6","xMinus",100,0.,1.); >> 67 >> 68 G4cout << "\n----> Histogram file is opened in " << fileName << G4endl; 60 } 69 } 61 70 62 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 63 72 64 RunAction::~RunAction() {} << 73 RunAction::~RunAction() >> 74 { >> 75 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); >> 76 >> 77 analysisManager->Write(); >> 78 analysisManager->CloseFile(); >> 79 >> 80 delete G4AnalysisManager::Instance(); >> 81 } 65 82 66 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 67 84 68 void RunAction::BeginOfRunAction(const G4Run* 85 void RunAction::BeginOfRunAction(const G4Run* aRun) 69 { << 86 { 70 G4cout << "### Run " << aRun->GetRunID() << 87 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; 71 << 88 72 // get material << 89 // save Rndm status 73 // << 90 G4RunManager::GetRunManager()->SetRandomNumberStore(true); 74 fMat = fDetector->GetMaterial(); << 91 CLHEP::HepRandom::showEngineStatus(); 75 G4cout << "###RunAction::BeginOfRunAction M << 76 << 77 fProcCounter = new ProcessesCount; 92 fProcCounter = new ProcessesCount; 78 << 79 fAnalysis = G4AnalysisManager::Instance(); << 80 fAnalysis->SetDefaultFileType("root"); << 81 << 82 // Open an output file << 83 // << 84 std::stringstream tmp; << 85 tmp << "testem6_" << aRun->GetRunID(); << 86 G4String fileName = tmp.str(); << 87 fAnalysis->OpenFile(fileName); << 88 fAnalysis->SetVerboseLevel(2); << 89 fAnalysis->SetActivation(true); << 90 << 91 // Creating histograms 1,2,3,4,5,6 << 92 // << 93 fAnalysis->SetFirstHistoId(1); << 94 fAnalysis->CreateH1("h1", "1/(1+(theta+[g]+) << 95 fAnalysis->CreateH1("h2", "log10(theta+ [g]+ << 96 fAnalysis->CreateH1("h3", "log10(theta- [g]- << 97 fAnalysis->CreateH1("h4", "log10(theta+ [g]+ << 98 fAnalysis->CreateH1("h5", "xPlus", 100, 0., << 99 fAnalysis->CreateH1("h6", "xMinus", 100, 0., << 100 << 101 // creating histogram 7,8,9,10,11 (CrossSect << 102 // << 103 G4double minBin = std::log10(fMinE / GeV); << 104 G4double maxBin = std::log10(fMaxE / GeV); << 105 fAnalysis->CreateH1("h7", "CrossSectionPerAt << 106 maxBin); << 107 fAnalysis->CreateH1("h8", "CrossSectionPerAt << 108 maxBin); << 109 fAnalysis->CreateH1("h9", "CrossSectionPerAt << 110 maxBin); << 111 fAnalysis->CreateH1("h10", "Theoretical Cros << 112 fnBin, minBin, maxBin); << 113 fAnalysis->CreateH1("h11", "Theoretical Cros << 114 minBin, maxBin); << 115 << 116 // creating histogram 12,13,14,15,16(CrossSe << 117 // << 118 fAnalysis->CreateH1("h12", "CrossSectionPerV << 119 fAnalysis->CreateH1("h13", "CrossSectionPerV << 120 fAnalysis->CreateH1("h14", "CrossSectionPerV << 121 fAnalysis->CreateH1("h15", "CrossSectionPerV << 122 maxBin); << 123 fAnalysis->CreateH1("h16", "CrossSectionPerV << 124 << 125 // creating histogram 17 (R ratio) << 126 fAnalysis->CreateH1("h17", "R : eeToHadr/eeT << 127 << 128 G4cout << "\n----> Histogram file is opened << 129 } 93 } 130 94 131 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 132 96 133 void RunAction::CountProcesses(G4String procNa 97 void RunAction::CountProcesses(G4String procName) 134 { 98 { 135 // does the process already encounted ? << 99 //does the process already encounted ? 136 // << 137 size_t nbProc = fProcCounter->size(); 100 size_t nbProc = fProcCounter->size(); 138 size_t i = 0; 101 size_t i = 0; 139 while ((i < nbProc) && ((*fProcCounter)[i]-> << 102 while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++; 140 i++; << 103 if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName)); 141 if (i == nbProc) fProcCounter->push_back(new << 104 142 << 143 (*fProcCounter)[i]->Count(); 105 (*fProcCounter)[i]->Count(); 144 } 106 } 145 107 146 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 147 109 148 void RunAction::EndOfRunAction(const G4Run*) 110 void RunAction::EndOfRunAction(const G4Run*) 149 { 111 { 150 G4cout << "### RunAction::EndOfRunAction" << << 112 // show Rndm status 151 // total number of process calls << 113 CLHEP::HepRandom::showEngineStatus(); 152 // << 114 //total number of process calls >> 115 G4double countTot = 0.; 153 G4cout << "\n Number of process calls --->"; 116 G4cout << "\n Number of process calls --->"; 154 for (size_t i = 0; i < fProcCounter->size(); << 117 for (size_t i=0; i< fProcCounter->size();i++) { 155 G4String procName = (*fProcCounter)[i]->Ge << 118 G4String procName = (*fProcCounter)[i]->GetName(); 156 if (procName != "Transportation") { << 119 if (procName != "Transportation") { 157 G4int count = (*fProcCounter)[i]->GetCou << 120 G4int count = (*fProcCounter)[i]->GetCounter(); 158 G4cout << "\t" << procName << " : " << c << 121 G4cout << "\t" << procName << " : " << count; 159 } << 122 countTot += count; 160 } << 123 } 161 G4cout << G4endl; << 162 << 163 // instance EmCalculator << 164 // << 165 G4EmCalculator emCal; << 166 emCal.SetVerbose(0); << 167 << 168 // get positron << 169 // << 170 G4String positronName = "e+"; << 171 G4ParticleDefinition* positron = G4ParticleT << 172 << 173 // process name << 174 // << 175 G4String annihilName = "annihil"; << 176 G4String annihiToMuName = "AnnihiToMuPair"; << 177 G4String annihiToHadrName = "ee2hadr"; << 178 G4String BremName = "eBrem"; << 179 G4String IoniName = "eIoni"; << 180 << 181 // get AnnihiToMuPair << 182 // << 183 G4AnnihiToMuPair* annihiToMu = << 184 reinterpret_cast<G4AnnihiToMuPair*>(emCal. << 185 << 186 // parameters for ComputeCrossSection << 187 // << 188 G4double atomicZ = 1.; << 189 G4double atomicA = 2.; << 190 << 191 // set parameters for theory << 192 // << 193 const G4ParticleDefinition* muon = G4MuonMin << 194 G4double Mu = muon->GetPDGMass(); << 195 G4double Me = electron_mass_c2; << 196 G4double Re = classic_electr_radius; << 197 G4double Ru = Re * Me / Mu; << 198 G4double Eth = 2 * Mu * Mu / Me - Me; << 199 G4PhysicsLogVector v(fMinE, fMaxE, fnBin, fa << 200 << 201 // Compute CrossSections and Fill histgrams << 202 // << 203 for (G4int i = 0; i <= fnBin; ++i) { << 204 G4double energy = v.Energy(i); << 205 G4double x = std::log10(energy / GeV); << 206 << 207 // CrossSectionPerAtom << 208 // << 209 G4double crs_annihiToMu = annihiToMu->Comp << 210 // G4cout << "crs_annihiToMu(mkb)=" << crs << 211 fAnalysis->FillH1(7, x, crs_annihiToMu / m << 212 << 213 G4double crs_annihil = << 214 emCal.ComputeCrossSectionPerAtom(energy, << 215 fAnalysis->FillH1(8, x, crs_annihil / micr << 216 << 217 G4double crs_annihiToHadr = << 218 emCal.ComputeCrossSectionPerAtom(energy, << 219 fAnalysis->FillH1(9, x, crs_annihiToHadr / << 220 << 221 // CrossSectionPerVolume << 222 // << 223 G4double crsVol_Brem = << 224 emCal.ComputeCrossSectionPerVolume(energ << 225 fAnalysis->FillH1(12, x, crsVol_Brem * mm) << 226 << 227 G4double crsVol_Ioni = << 228 emCal.ComputeCrossSectionPerVolume(energ << 229 fAnalysis->FillH1(13, x, crsVol_Ioni * mm) << 230 << 231 G4double crsVol_annihiToMu = annihiToMu->C << 232 fAnalysis->FillH1(14, x, crsVol_annihiToMu << 233 << 234 G4double crsVol_annihil = << 235 emCal.ComputeCrossSectionPerVolume(energ << 236 fAnalysis->FillH1(15, x, crsVol_annihil * << 237 << 238 G4double crsVol_annihiToHadr = << 239 emCal.ComputeCrossSectionPerVolume(energ << 240 fAnalysis->FillH1(16, x, crsVol_annihiToHa << 241 << 242 // R ratio << 243 // << 244 G4double RR = 0.0; << 245 if (crsVol_annihiToMu > 0.) RR = crsVol_an << 246 fAnalysis->FillH1(17, x, RR); << 247 << 248 // Theoretical calculation << 249 // << 250 G4double X1 = energy / Me; << 251 if (X1 > 1 && i % 1000 == 0) { << 252 G4double crs_annihil_theory = << 253 atomicZ * pi * Re * Re << 254 * ((X1 * X1 + 4 * X1 + 1) * G4Log(X1 + << 255 - (X1 + 3) / std::sqrt(X1 * X1 - 1) << 256 / (X1 + 1); << 257 fAnalysis->FillH1(10, x, crs_annihil_the << 258 } << 259 << 260 G4double X2 = Eth / energy; << 261 if (X2 < 1. && i % 1000 == 0) { << 262 G4double crs_annihiToMu_theory = << 263 atomicZ * pi * Ru * Ru / 3 * X2 * (1 + << 264 fAnalysis->FillH1(11, x, crs_annihiToMu_ << 265 } << 266 << 267 // if(i%1000==0)G4cout <<"###energy:" << e << 268 // << crs_annihiToMu << "/crs_ToTw << 269 // <<"/crs_ToToHadr:"<<crs_annihiT << 270 } 124 } 271 << 272 fAnalysis->Write(); << 273 fAnalysis->CloseFile(); << 274 fAnalysis->Clear(); << 275 << 276 G4cout << G4endl; 125 G4cout << G4endl; >> 126 277 } 127 } 278 128 279 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 280 130