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 biasing/ReverseMC01/src/RMC01Analysi 27 /// \brief Implementation of the RMC01Analysis 28 // 29 // 30 ////////////////////////////////////////////// 31 // Class Name: RMC01AnalysisManage 32 // Author: L. Desorgher 33 // Organisation: SpaceIT GmbH 34 // Contract: ESA contract 21435/ 35 // Customer: ESA/ESTEC 36 ////////////////////////////////////////////// 37 38 //....oooOO0OOooo........oooOO0OOooo........oo 39 //....oooOO0OOooo........oooOO0OOooo........oo 40 41 #include "RMC01AnalysisManager.hh" 42 43 #include "RMC01AnalysisManagerMessenger.hh" 44 #include "RMC01SD.hh" 45 46 #include "G4AdjointSimManager.hh" 47 #include "G4Electron.hh" 48 #include "G4Gamma.hh" 49 #include "G4PhysicalConstants.hh" 50 #include "G4Proton.hh" 51 #include "G4RunManager.hh" 52 #include "G4SDManager.hh" 53 #include "G4SystemOfUnits.hh" 54 #include "G4THitsCollection.hh" 55 #include "G4Timer.hh" 56 57 RMC01AnalysisManager* RMC01AnalysisManager::fI 58 59 //....oooOO0OOooo........oooOO0OOooo........oo 60 61 RMC01AnalysisManager::RMC01AnalysisManager() 62 : fAccumulated_edep(0.), 63 fAccumulated_edep2(0.), 64 fMean_edep(0.), 65 fError_mean_edep(0.), 66 fRelative_error(0.), 67 fElapsed_time(0.), 68 fPrecision_to_reach(0.), 69 fStop_run_if_precision_reached(true), 70 fNb_evt_modulo_for_convergence_test(5000), 71 fEdep_rmatrix_vs_electron_prim_energy(0), 72 fElectron_current_rmatrix_vs_electron_prim 73 fGamma_current_rmatrix_vs_electron_prim_en 74 fEdep_rmatrix_vs_gamma_prim_energy(0), 75 fElectron_current_rmatrix_vs_gamma_prim_en 76 fGamma_current_rmatrix_vs_gamma_prim_energ 77 fEdep_rmatrix_vs_proton_prim_energy(0), 78 fElectron_current_rmatrix_vs_proton_prim_e 79 fProton_current_rmatrix_vs_proton_prim_ene 80 fGamma_current_rmatrix_vs_proton_prim_ener 81 fFactoryOn(false), 82 fPrimSpectrumType(EXPO), 83 fAlpha_or_E0(.5 * MeV), 84 fAmplitude_prim_spectrum(1.), 85 fEmin_prim_spectrum(1. * keV), 86 fEmax_prim_spectrum(20. * MeV), 87 fAdjoint_sim_mode(true), 88 fNb_evt_per_adj_evt(2) 89 { 90 fMsg = new RMC01AnalysisManagerMessenger(thi 91 92 //------------- 93 // Timer for convergence vector 94 //------------- 95 96 fTimer = new G4Timer(); 97 98 //--------------------------------- 99 // Primary particle ID for normalisation of 100 //--------------------------------- 101 102 fPrimPDG_ID = G4Electron::Electron()->GetPDG 103 104 fFileName[0] = "sim"; 105 } 106 107 //....oooOO0OOooo........oooOO0OOooo........oo 108 109 RMC01AnalysisManager::~RMC01AnalysisManager() 110 { 111 ; 112 } 113 114 //....oooOO0OOooo........oooOO0OOooo........oo 115 116 RMC01AnalysisManager* RMC01AnalysisManager::Ge 117 { 118 if (fInstance == 0) fInstance = new RMC01Ana 119 return fInstance; 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oo 123 124 void RMC01AnalysisManager::BeginOfRun(const G4 125 { 126 fAccumulated_edep = 0.; 127 fAccumulated_edep2 = 0.; 128 fNentry = 0.0; 129 fRelative_error = 1.; 130 fMean_edep = 0.; 131 fError_mean_edep = 0.; 132 fAdjoint_sim_mode = G4AdjointSimManager::Get 133 134 if (fAdjoint_sim_mode) { 135 fNb_evt_per_adj_evt = aRun->GetNumberOfEve 136 / G4AdjointSimManage 137 fConvergenceFileOutput.open("ConvergenceOf 138 fConvergenceFileOutput << "Normalised Edep 139 } 140 else { 141 fConvergenceFileOutput.open("ConvergenceOf 142 fConvergenceFileOutput << "Edep per event 143 } 144 fConvergenceFileOutput.setf(std::ios::scient 145 fConvergenceFileOutput.precision(6); 146 147 fTimer->Start(); 148 fElapsed_time = 0.; 149 150 Book(); 151 } 152 153 //....oooOO0OOooo........oooOO0OOooo........oo 154 155 void RMC01AnalysisManager::EndOfRun(const G4Ru 156 { 157 fTimer->Stop(); 158 G4int nb_evt = aRun->GetNumberOfEvent(); 159 G4double factor = 1. / nb_evt; 160 if (!fAdjoint_sim_mode) { 161 G4cout << "Results of forward simulation!" 162 G4cout << "edep per event [MeV] = " << fMe 163 G4cout << "error[MeV] = " << fError_mean_e 164 } 165 166 else { 167 G4cout << "Results of reverse/adjoint simu 168 G4cout << "normalised edep [MeV] = " << fM 169 G4cout << "error[MeV] = " << fError_mean_e 170 factor = 1. * G4AdjointSimManager::GetInst 171 / aRun->GetNumberOfEvent(); 172 } 173 Save(factor); 174 fConvergenceFileOutput.close(); 175 } 176 177 //....oooOO0OOooo........oooOO0OOooo........oo 178 179 void RMC01AnalysisManager::BeginOfEvent(const 180 { 181 ; 182 } 183 184 //....oooOO0OOooo........oooOO0OOooo........oo 185 186 void RMC01AnalysisManager::EndOfEvent(const G4 187 { 188 if (fAdjoint_sim_mode) 189 EndOfEventForAdjointSimulation(anEvent); 190 else 191 EndOfEventForForwardSimulation(anEvent); 192 193 // Test convergence. The error is already co 194 //-------------------------------------- 195 G4int nb_event = anEvent->GetEventID() + 1; 196 // G4double factor=1.; 197 if (fAdjoint_sim_mode) { 198 G4double n_adj_evt = nb_event / fNb_evt_pe 199 // nb_event/fNb_evt_per_adj_evt; 200 if (n_adj_evt * fNb_evt_per_adj_evt == nb_ 201 nb_event = static_cast<G4int>(n_adj_evt) 202 } 203 else 204 nb_event = 0; 205 } 206 207 if (nb_event > 100 && fStop_run_if_precision 208 G4cout << fPrecision_to_reach * 100. << "% 209 fTimer->Stop(); 210 fElapsed_time += fTimer->GetRealElapsed(); 211 fConvergenceFileOutput << fMean_edep << '\ 212 << std::endl; 213 G4RunManager::GetRunManager()->AbortRun(tr 214 } 215 216 if (nb_event > 0 && nb_event % fNb_evt_modul 217 fTimer->Stop(); 218 fElapsed_time += fTimer->GetRealElapsed(); 219 fTimer->Start(); 220 fConvergenceFileOutput << fMean_edep << '\ 221 << std::endl; 222 } 223 } 224 225 //....oooOO0OOooo........oooOO0OOooo........oo 226 227 void RMC01AnalysisManager::EndOfEventForForwar 228 { 229 G4SDManager* SDman = G4SDManager::GetSDMpoin 230 G4HCofThisEvent* HCE = anEvent->GetHCofThisE 231 RMC01DoubleWithWeightHitsCollection* edepCol 232 (RMC01DoubleWithWeightHitsCollection*)(HCE 233 234 RMC01DoubleWithWeightHitsCollection* electro 235 (RMC01DoubleWithWeightHitsCollection*)(HCE 236 237 RMC01DoubleWithWeightHitsCollection* protonC 238 (RMC01DoubleWithWeightHitsCollection*)(HCE 239 240 RMC01DoubleWithWeightHitsCollection* gammaCu 241 (RMC01DoubleWithWeightHitsCollection*)(HCE 242 243 // Total energy deposited in Event 244 //------------------------------- 245 G4double totEdep = 0; 246 std::size_t i; 247 for (i = 0; i < edepCollection->entries(); + 248 totEdep += (*edepCollection)[i]->GetValue( 249 250 if (totEdep > 0.) { 251 fAccumulated_edep += totEdep; 252 fAccumulated_edep2 += totEdep * totEdep; 253 fNentry += 1.0; 254 G4PrimaryParticle* thePrimary = anEvent->G 255 G4double E0 = thePrimary->GetG4code()->Get 256 G4double P = thePrimary->GetMomentum().mag 257 G4double prim_ekin = std::sqrt(E0 * E0 + P 258 fEdep_vs_prim_ekin->fill(prim_ekin, totEde 259 } 260 ComputeMeanEdepAndError(anEvent, fMean_edep, 261 if (fError_mean_edep > 0) fRelative_error = 262 263 // Particle current on sensitive cylinder 264 //------------------------------------- 265 266 for (i = 0; i < electronCurrentCollection->e 267 G4double ekin = (*electronCurrentCollectio 268 G4double weight = (*electronCurrentCollect 269 fElectron_current->fill(ekin, weight); 270 } 271 272 for (i = 0; i < protonCurrentCollection->ent 273 G4double ekin = (*protonCurrentCollection) 274 G4double weight = (*protonCurrentCollectio 275 fProton_current->fill(ekin, weight); 276 } 277 278 for (i = 0; i < gammaCurrentCollection->entr 279 G4double ekin = (*gammaCurrentCollection)[ 280 G4double weight = (*gammaCurrentCollection 281 fGamma_current->fill(ekin, weight); 282 } 283 } 284 285 //....oooOO0OOooo........oooOO0OOooo........oo 286 287 void RMC01AnalysisManager::EndOfEventForAdjoin 288 { 289 // Output from Sensitive volume computed dur 290 //------------------------------------------ 291 G4SDManager* SDman = G4SDManager::GetSDMpoin 292 G4HCofThisEvent* HCE = anEvent->GetHCofThisE 293 RMC01DoubleWithWeightHitsCollection* edepCol 294 (RMC01DoubleWithWeightHitsCollection*)(HCE 295 296 RMC01DoubleWithWeightHitsCollection* electro 297 (RMC01DoubleWithWeightHitsCollection*)(HCE 298 299 RMC01DoubleWithWeightHitsCollection* protonC 300 (RMC01DoubleWithWeightHitsCollection*)(HCE 301 302 RMC01DoubleWithWeightHitsCollection* gammaCu 303 (RMC01DoubleWithWeightHitsCollection*)(HCE 304 305 // Computation of total energy deposited in 306 //------------------------------- 307 G4double totEdep = 0; 308 std::size_t i; 309 for (i = 0; i < edepCollection->entries(); + 310 totEdep += (*edepCollection)[i]->GetValue( 311 312 // Output from adjoint tracking phase 313 //------------------------------------------ 314 315 G4AdjointSimManager* theAdjointSimManager = 316 317 size_t nb_adj_track = theAdjointSimManager-> 318 G4double total_normalised_weight = 0.; 319 320 // We need to loop over the adjoint tracks t 321 // surface. 322 for (std::size_t j = 0; j < nb_adj_track; ++ 323 G4int pdg_nb = theAdjointSimManager->GetFw 324 G4double prim_ekin = theAdjointSimManager- 325 G4double adj_weight = theAdjointSimManager 326 327 // Factor of normalisation to user defined 328 //---------------------------------------- 329 G4double normalised_weight = 0.; 330 if (pdg_nb == fPrimPDG_ID && prim_ekin >= 331 && prim_ekin <= fEmax_prim_spectrum) 332 normalised_weight = adj_weight * PrimDif 333 total_normalised_weight += normalised_weig 334 335 // Answer matrices 336 //------------- 337 G4H1* edep_rmatrix = 0; 338 G4H2* electron_current_rmatrix = 0; 339 G4H2* gamma_current_rmatrix = 0; 340 G4H2* proton_current_rmatrix = 0; 341 342 if (pdg_nb == G4Electron::Electron()->GetP 343 edep_rmatrix = fEdep_rmatrix_vs_electron 344 electron_current_rmatrix = fElectron_cur 345 gamma_current_rmatrix = fGamma_current_r 346 } 347 else if (pdg_nb == G4Gamma::Gamma()->GetPD 348 // gammma answer matrices 349 edep_rmatrix = fEdep_rmatrix_vs_gamma_pr 350 electron_current_rmatrix = fElectron_cur 351 gamma_current_rmatrix = fGamma_current_r 352 } 353 else if (pdg_nb == G4Proton::Proton()->Get 354 // proton answer matrices 355 edep_rmatrix = fEdep_rmatrix_vs_proton_p 356 electron_current_rmatrix = fElectron_cur 357 gamma_current_rmatrix = fGamma_current_r 358 proton_current_rmatrix = fProton_current 359 } 360 // Register histo edep vs prim ekin 361 //---------------------------------- 362 if (normalised_weight > 0) fEdep_vs_prim_e 363 // Registering answer matrix 364 //--------------------------- 365 edep_rmatrix->fill(prim_ekin, totEdep * ad 366 367 // Registering of current of particles on 368 //---------------------------------------- 369 370 for (i = 0; i < electronCurrentCollection- 371 G4double ekin = (*electronCurrentCollect 372 G4double weight = (*electronCurrentColle 373 fElectron_current->fill(ekin, weight * n 374 electron_current_rmatrix->fill(prim_ekin 375 } 376 for (i = 0; i < protonCurrentCollection->e 377 G4double ekin = (*protonCurrentCollectio 378 G4double weight = (*protonCurrentCollect 379 fProton_current->fill(ekin, weight * nor 380 proton_current_rmatrix->fill(prim_ekin, 381 } 382 for (i = 0; i < gammaCurrentCollection->en 383 G4double ekin = (*gammaCurrentCollection 384 G4double weight = (*gammaCurrentCollecti 385 fGamma_current->fill(ekin, weight * norm 386 gamma_current_rmatrix->fill(prim_ekin, e 387 } 388 } 389 390 // Registering of total energy deposited in 391 //------------------------------- 392 G4bool new_mean_computed = false; 393 if (totEdep > 0.) { 394 if (total_normalised_weight > 0.) { 395 G4double edep = totEdep * total_normalis 396 397 // Check if the edep is not wrongly too 398 //-------------------------------------- 399 G4double new_mean(0.0), new_error(0.0); 400 fAccumulated_edep += edep; 401 fAccumulated_edep2 += edep * edep; 402 fNentry += 1.0; 403 ComputeMeanEdepAndError(anEvent, new_mea 404 G4double new_relative_error = 1.; 405 if (new_error > 0) new_relative_error = 406 if (fRelative_error < 0.10 && new_relati 407 G4cout << "Potential wrong adjoint wei 408 G4cout << "The results of this event w 409 G4cout << "previous mean edep [MeV] " 410 G4cout << "previous relative error " < 411 G4cout << "new rejected mean edep [MeV 412 G4cout << "new rejected relative error 413 fAccumulated_edep -= edep; 414 fAccumulated_edep2 -= edep * edep; 415 fNentry -= 1.0; 416 return; 417 } 418 else { // accepted 419 fMean_edep = new_mean; 420 fError_mean_edep = new_error; 421 fRelative_error = new_relative_error; 422 new_mean_computed = true; 423 } 424 } 425 426 if (!new_mean_computed) { 427 ComputeMeanEdepAndError(anEvent, fMean_e 428 fRelative_error = (fMean_edep > 0.0) ? f 429 } 430 } 431 } 432 433 //....oooOO0OOooo........oooOO0OOooo........oo 434 435 G4double RMC01AnalysisManager::PrimDiffAndDirF 436 { 437 G4double flux = fAmplitude_prim_spectrum; 438 if (fPrimSpectrumType == EXPO) 439 flux *= std::exp(-prim_energy / fAlpha_or_ 440 else 441 flux *= std::pow(prim_energy, -fAlpha_or_E 442 return flux; 443 } 444 445 //....oooOO0OOooo........oooOO0OOooo........oo 446 /* 447 void RMC01AnalysisManager::WriteHisto(G4H1* a 448 G4double scaling_factor, G4String 449 { std::fstream FileOutput(fileName, std::ios:: 450 FileOutput<<header_lines; 451 FileOutput.setf(std::ios::scientific); 452 FileOutput.precision(6); 453 454 for (G4int i =0;i<G4int(anHisto->axis().bins 455 FileOutput<<anHisto->axis().bin_lower_ 456 <<'\t'<<anHisto->axis().bin_uppe 457 <<'\t'<<anHisto->bin_height(i)*s 458 <<'\t'<<anHisto->bin_error(i)*sc 459 } 460 } 461 462 //....oooOO0OOooo........oooOO0OOooo........oo 463 464 void RMC01AnalysisManager::WriteHisto(G4H2* a 465 G4double scaling_factor, G4String 466 { std::fstream FileOutput(fileName, std::ios:: 467 FileOutput<<header_lines; 468 469 FileOutput.setf(std::ios::scientific); 470 FileOutput.precision(6); 471 472 for (G4int i =0;i<G4int(anHisto->axis_x().bi 473 for (G4int j =0;j<G4int(anHisto->axis_y(). 474 FileOutput<<anHisto->axis_x().bin_lower 475 <<'\t'<<anHisto->axis_x() 476 <<'\t'<<anHisto->axis_y().b 477 <<'\t'<<anHisto->axis_y().b 478 <<'\t'<<anHisto->bin_height 479 <<'\t'<<anHisto->bin_error(i,j 480 481 } 482 } 483 } 484 */ 485 //....oooOO0OOooo........oooOO0OOooo........oo 486 487 void RMC01AnalysisManager::ComputeMeanEdepAndE 488 489 { 490 G4int nb_event = anEvent->GetEventID() + 1; 491 G4double factor = 1.; 492 if (fAdjoint_sim_mode) { 493 nb_event /= fNb_evt_per_adj_evt; 494 factor = 1. * G4AdjointSimManager::GetInst 495 } 496 497 // VI: error computation now is based on num 498 // number of events 499 // LD: This is wrong! With the use of fNentr 500 // correctly normalised. The mean and th 501 // with nb_event. The old computation ha 502 // VI: OK, but let computations be double 503 if (nb_event > 0) { 504 G4double norm = 1.0 / (G4double)nb_event; 505 mean = fAccumulated_edep * norm; 506 G4double mean_x2 = fAccumulated_edep2 * no 507 G4double zz = mean_x2 - mean * mean; 508 /* 509 G4cout << "Nevt= " << nb_event << " mean 510 << " mean_x2= " << mean_x2 << " 511 << zz << G4endl; 512 */ 513 error = factor * std::sqrt(std::max(zz, 0. 514 mean *= factor; 515 } 516 else { 517 mean = 0; 518 error = 0; 519 } 520 // G4cout << "Aend: " << mean << " " << erro 521 } 522 523 //....oooOO0OOooo........oooOO0OOooo........oo 524 525 void RMC01AnalysisManager::SetPrimaryExpSpectr 526 527 528 { 529 fPrimSpectrumType = EXPO; 530 if (particle_name == "e-") 531 fPrimPDG_ID = G4Electron::Electron()->GetP 532 else if (particle_name == "gamma") 533 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco 534 else if (particle_name == "proton") 535 fPrimPDG_ID = G4Proton::Proton()->GetPDGEn 536 else { 537 G4cout << "The particle that you did selec 538 << "list for primary [e-, gamma, pr 539 return; 540 } 541 fAlpha_or_E0 = E0; 542 fAmplitude_prim_spectrum = 543 omni_fluence / E0 / (std::exp(-Emin / E0) 544 fEmin_prim_spectrum = Emin; 545 fEmax_prim_spectrum = Emax; 546 } 547 548 //....oooOO0OOooo........oooOO0OOooo........oo 549 550 void RMC01AnalysisManager::SetPrimaryPowerLawS 551 552 553 554 { 555 fPrimSpectrumType = POWER; 556 if (particle_name == "e-") 557 fPrimPDG_ID = G4Electron::Electron()->GetP 558 else if (particle_name == "gamma") 559 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco 560 else if (particle_name == "proton") 561 fPrimPDG_ID = G4Proton::Proton()->GetPDGEn 562 else { 563 G4cout << "The particle that you did selec 564 << " list for primary [e-, gamma, p 565 return; 566 } 567 568 if (alpha == 1.) { 569 fAmplitude_prim_spectrum = omni_fluence / 570 } 571 else { 572 G4double p = 1. - alpha; 573 fAmplitude_prim_spectrum = omni_fluence / 574 } 575 576 fAlpha_or_E0 = alpha; 577 fEmin_prim_spectrum = Emin; 578 fEmax_prim_spectrum = Emax; 579 } 580 581 //....oooOO0OOooo........oooOO0OOooo........oo 582 583 void RMC01AnalysisManager::Book() 584 { 585 //---------------------- 586 // Creation of histograms 587 //---------------------- 588 589 // Energy binning of the histograms : 60 log 590 591 G4double emin = 1. * keV; 592 G4double emax = 1. * GeV; 593 594 // file_name 595 fFileName[0] = "forward_sim"; 596 if (fAdjoint_sim_mode) fFileName[0] = "adjoi 597 598 // Histo manager 599 G4AnalysisManager* theHistoManager = G4Analy 600 theHistoManager->SetDefaultFileType("root"); 601 G4String extension = theHistoManager->GetFil 602 fFileName[1] = fFileName[0] + "." + extensio 603 theHistoManager->SetFirstHistoId(1); 604 605 G4bool fileOpen = theHistoManager->OpenFile( 606 if (!fileOpen) { 607 G4cout << "\n---> RMC01AnalysisManager::Bo 608 return; 609 } 610 611 // Create directories 612 theHistoManager->SetHistoDirectoryName("hist 613 614 // Histograms for : 615 // 1)the forward simulation results 616 // 2)the Reverse MC simulation resu 617 //------------------------------------------ 618 619 G4int idHisto = 620 theHistoManager->CreateH1(G4String("Edep_v 621 60, emin, emax, 622 fEdep_vs_prim_ekin = theHistoManager->GetH1( 623 624 idHisto = theHistoManager->CreateH1(G4String 625 emax, "n 626 627 fElectron_current = theHistoManager->GetH1(i 628 629 idHisto = theHistoManager->CreateH1(G4String 630 emax, "n 631 fProton_current = theHistoManager->GetH1(idH 632 633 idHisto = theHistoManager->CreateH1(G4String 634 "none", 635 fGamma_current = theHistoManager->GetH1(idHi 636 637 // Response matrices for the adjoint simulat 638 //------------------------------------------ 639 if (fAdjoint_sim_mode) { 640 // Response matrices for external isotropi 641 //---------------------------------------- 642 643 idHisto = theHistoManager->CreateH1(G4Stri 644 G4Stri 645 emax, 646 fEdep_rmatrix_vs_electron_prim_energy = th 647 648 idHisto = theHistoManager->CreateH2( 649 G4String("Electron_current_rmatrix_vs_el 650 G4String("electron current RM vs e- pri 651 "none", "none", "none", G4String("log"), 652 653 fElectron_current_rmatrix_vs_electron_prim 654 655 idHisto = theHistoManager->CreateH2(G4Stri 656 G4Stri 657 emin, 658 G4Stri 659 660 fGamma_current_rmatrix_vs_electron_prim_en 661 662 // Response matrices for external isotropi 663 664 idHisto = theHistoManager->CreateH1(G4Stri 665 G4Stri 666 emax, 667 fEdep_rmatrix_vs_gamma_prim_energy = theHi 668 669 idHisto = theHistoManager->CreateH2(G4Stri 670 G4Stri 671 60, em 672 "none" 673 674 fElectron_current_rmatrix_vs_gamma_prim_en 675 676 idHisto = theHistoManager->CreateH2(G4Stri 677 G4Stri 678 emin, 679 G4Stri 680 681 fGamma_current_rmatrix_vs_gamma_prim_energ 682 683 // Response matrices for external isotropi 684 idHisto = theHistoManager->CreateH1(G4Stri 685 G4Stri 686 emax, 687 fEdep_rmatrix_vs_proton_prim_energy = theH 688 689 idHisto = theHistoManager->CreateH2(G4Stri 690 G4Stri 691 60, em 692 "none" 693 694 fElectron_current_rmatrix_vs_proton_prim_e 695 696 idHisto = theHistoManager->CreateH2(G4Stri 697 G4Stri 698 emin, 699 G4Stri 700 701 fGamma_current_rmatrix_vs_proton_prim_ener 702 703 idHisto = theHistoManager->CreateH2(G4Stri 704 G4Stri 705 emin, 706 G4Stri 707 708 fProton_current_rmatrix_vs_proton_prim_ene 709 } 710 fFactoryOn = true; 711 G4cout << "\n----> Histogram Tree is opened 712 } 713 714 //....oooOO0OOooo........oooOO0OOooo........oo 715 716 void RMC01AnalysisManager::Save(G4double scali 717 { 718 if (fFactoryOn) { 719 G4AnalysisManager* theHistoManager = G4Ana 720 // scaling of results 721 //----------------- 722 723 for (G4int ind = 1; ind <= theHistoManager 724 theHistoManager->SetH1Ascii(ind, true); 725 theHistoManager->ScaleH1(ind, scaling_fa 726 } 727 for (G4int ind = 1; ind <= theHistoManager 728 theHistoManager->ScaleH2(ind, scaling_fa 729 730 theHistoManager->Write(); 731 theHistoManager->CloseFile(); 732 G4cout << "\n----> Histogram Tree is saved 733 734 theHistoManager->Clear(); 735 fFactoryOn = false; 736 } 737 } 738