Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 /// \file electromagnetic/TestEm6/src/RunActio 27 /// \brief Implementation of the RunAction cla 28 // 29 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 33 #include "RunAction.hh" 34 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" 43 #include "G4RunManager.hh" 44 #include "G4SystemOfUnits.hh" 45 #include "G4eBremsstrahlung.hh" 46 #include "G4eeToHadrons.hh" 47 #include "G4eeToHadronsModel.hh" 48 #include "Randomize.hh" 49 50 #include <sstream> 51 52 //....oooOO0OOooo........oooOO0OOooo........oo 53 54 RunAction::RunAction(DetectorConstruction* det 55 : G4UserRunAction(), fDetector(det), fProcCo 56 { 57 fMinE = 40 * GeV; 58 fMaxE = 10000 * GeV; 59 fnBin = 10000; 60 } 61 62 //....oooOO0OOooo........oooOO0OOooo........oo 63 64 RunAction::~RunAction() {} 65 66 //....oooOO0OOooo........oooOO0OOooo........oo 67 68 void RunAction::BeginOfRunAction(const G4Run* 69 { 70 G4cout << "### Run " << aRun->GetRunID() << 71 72 // get material 73 // 74 fMat = fDetector->GetMaterial(); 75 G4cout << "###RunAction::BeginOfRunAction M 76 77 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 } 130 131 //....oooOO0OOooo........oooOO0OOooo........oo 132 133 void RunAction::CountProcesses(G4String procNa 134 { 135 // does the process already encounted ? 136 // 137 size_t nbProc = fProcCounter->size(); 138 size_t i = 0; 139 while ((i < nbProc) && ((*fProcCounter)[i]-> 140 i++; 141 if (i == nbProc) fProcCounter->push_back(new 142 143 (*fProcCounter)[i]->Count(); 144 } 145 146 //....oooOO0OOooo........oooOO0OOooo........oo 147 148 void RunAction::EndOfRunAction(const G4Run*) 149 { 150 G4cout << "### RunAction::EndOfRunAction" << 151 // total number of process calls 152 // 153 G4cout << "\n Number of process calls --->"; 154 for (size_t i = 0; i < fProcCounter->size(); 155 G4String procName = (*fProcCounter)[i]->Ge 156 if (procName != "Transportation") { 157 G4int count = (*fProcCounter)[i]->GetCou 158 G4cout << "\t" << procName << " : " << c 159 } 160 } 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 } 271 272 fAnalysis->Write(); 273 fAnalysis->CloseFile(); 274 fAnalysis->Clear(); 275 276 G4cout << G4endl; 277 } 278 279 //....oooOO0OOooo........oooOO0OOooo........oo 280