API tools faq
paste
Advertisement
luizrosalba

GEANT4 - Run.cpp

luizrosalba
Nov 17th, 2015
109
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 20.74 KB | None | 0 0
raw download clone embed print report
  1. //
  2. // ********************************************************************
  3. // * License and Disclaimer *
  4. // * *
  5. // * The Geant4 software is copyright of the Copyright Holders of *
  6. // * the Geant4 Collaboration. It is provided under the terms and *
  7. // * conditions of the Geant4 Software License, included in the file *
  8. // * LICENSE and available at http://cern.ch/geant4/license . These *
  9. // * include a list of copyright holders. *
  10. // * *
  11. // * Neither the authors of this software system, nor their employing *
  12. // * institutes,nor the agencies providing financial support for this *
  13. // * work make any representation or warranty, express or implied, *
  14. // * regarding this software system or assume any liability for its *
  15. // * use. Please see the license in the file LICENSE and URL above *
  16. // * for the full disclaimer and the limitation of liability. *
  17. // * *
  18. // * This code implementation is the result of the scientific and *
  19. // * technical work of the GEANT4 collaboration. *
  20. // * By using, copying, modifying or distributing the software (or *
  21. // * any work based on the software) you agree to acknowledge its *
  22. // * use in resulting scientific publications, and indicate your *
  23. // * acceptance of all terms of the Geant4 Software license. *
  24. // ********************************************************************
  25. //
  26. /// \file electromagnetic/TestEm11/src/Run.cc
  27. /// \brief Implementation of the Run class
  28. //
  29. // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $
  30. //
  31. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  32. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  33.  
  34. #include "Run.hh"
  35. #include "DetectorConstruction.hh"
  36. #include "PrimaryGeneratorAction.hh"
  37.  
  38. #include "G4EmCalculator.hh"
  39. #include "G4SystemOfUnits.hh"
  40. #include "G4UnitsTable.hh"
  41. #include "G4HCofThisEvent.hh"
  42. #include "G4Step.hh"
  43. #include "G4ThreeVector.hh"
  44. #include "G4ios.hh"
  45. #include "G4SDManager.hh"
  46. #include "G4VSensitiveDetector.hh"
  47. #include "G4THitsMap.hh"
  48. #include <iomanip>
  49.  
  50. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  51.  
  52. Run::Run(DetectorConstruction* det)
  53. : G4Run(),
  54. fDetector(det),
  55. fParticle(0), fEkin(0.)
  56. {
  57. fEnergyDeposit = fEnergyDeposit2 = 0.;
  58. fTrakLenCharged = fTrakLenCharged2 = 0.;
  59. fTrakLenNeutral = fTrakLenNeutral2 = 0.;
  60. fNbStepsCharged = fNbStepsCharged2 = 0.;
  61. fNbStepsNeutral = fNbStepsNeutral2 = 0.;
  62. fMscProjecTheta = fMscProjecTheta2 = 0.;
  63. fMscThetaCentral = 0.;
  64.  
  65. fNbGamma = fNbElect = fNbPosit = 0;
  66.  
  67. fTransmit[0] = fTransmit[1] = fReflect[0] = fReflect[1] = 0;
  68.  
  69. fMscEntryCentral = 0;
  70.  
  71. fEnergyLeak[0] = fEnergyLeak[1] = fEnergyLeak2[0] = fEnergyLeak2[1] = 0.;
  72.  
  73.  
  74. hitsCollID = -1;
  75. /// gerando o vetor de energias discretas
  76. energia_maxima=30;
  77.  
  78. tamanho = 10000; /// 100 = duas casas deciamis
  79. tamanho_vetor = int(tamanho); /// calcula o tamanho do vetor
  80. /// pegando tres n decimais do valor da energia onde tamanho = 10^n
  81. energias_disc.assign (tamanho_vetor,0.000); /// cria o vetor e preenche com zeros
  82. energias_continuo.assign (tamanho_vetor,0.0000);
  83.  
  84. //Calibra_vetores_energia( energias_disc , energias_continuo , energia_maxima);
  85.  
  86.  
  87. //// chamando a colecao com o nome correto
  88. G4SDManager * SDman = G4SDManager::GetSDMpointer();
  89. if(hitsCollID<0){
  90. G4String colNam;
  91. hitsCollID = SDman->GetCollectionID(colNam="SD");
  92. }
  93. }
  94.  
  95. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  96.  
  97. Run::~Run()
  98. { }
  99.  
  100. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  101.  
  102. void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
  103. {
  104. fParticle = particle;
  105. fEkin = energy;
  106.  
  107. //compute theta0
  108. fMscThetaCentral = 3*ComputeMscHighland();
  109. }
  110.  
  111. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  112.  
  113. void Run::Merge(const G4Run* run)
  114. {
  115. const Run* localRun = static_cast<const Run*>(run);
  116.  
  117. // pass information about primary particle
  118. fParticle = localRun->fParticle;
  119. fEkin = localRun->fEkin;
  120.  
  121. fMscThetaCentral = localRun->fMscThetaCentral;
  122.  
  123. // accumulate sums
  124. //
  125. fEnergyDeposit += localRun->fEnergyDeposit;
  126. fEnergyDeposit2 += localRun->fEnergyDeposit2;
  127. fTrakLenCharged += localRun->fTrakLenCharged;
  128. fTrakLenCharged2 += localRun->fTrakLenCharged2;
  129. fTrakLenNeutral += localRun->fTrakLenNeutral;
  130. fTrakLenNeutral2 += localRun->fTrakLenNeutral2;
  131. fNbStepsCharged += localRun->fNbStepsCharged;
  132. fNbStepsCharged2 += localRun->fNbStepsCharged2;
  133. fNbStepsNeutral += localRun->fNbStepsNeutral;
  134. fNbStepsNeutral2 += localRun->fNbStepsNeutral2;
  135.  
  136. fNbGamma += localRun->fNbGamma;
  137. fNbElect += localRun->fNbElect;
  138. fNbPosit += localRun->fNbPosit;
  139.  
  140. fTransmit[0] += localRun->fTransmit[0];
  141. fTransmit[1] += localRun->fTransmit[1];
  142. fReflect[0] += localRun->fReflect[0];
  143. fReflect[1] += localRun->fReflect[1];
  144.  
  145. fMscEntryCentral += localRun->fMscEntryCentral;
  146.  
  147. fEnergyLeak[0] += localRun->fEnergyLeak[0];
  148. fEnergyLeak[1] += localRun->fEnergyLeak[1];
  149. fEnergyLeak2[0] += localRun->fEnergyLeak2[0];
  150. fEnergyLeak2[1] += localRun->fEnergyLeak2[1];
  151.  
  152.  
  153. G4cout<< " Rodando o Merge " << G4endl;
  154.  
  155.  
  156. /// unindo os resultados de cada thread no master energias_disc pertence ao master
  157. G4cout << " Unindo os resultados das threads no master " << G4endl;
  158. for (unsigned int i = 0 ; i < localRun->energias_disc.size() ; i++ )
  159. {
  160. energias_disc[i] += localRun->energias_disc.at(i);
  161. }
  162. G4cout << " Criando o vetor continuo no master " << G4endl;
  163.  
  164. RealizaConvolucao(energias_disc);
  165.  
  166.  
  167. G4Run::Merge(run);
  168. }
  169.  
  170.  
  171. std::vector<double> Run::detector_efficiency(int Z,double thickness, int Z_window,double window_thick, double Emin, double &Emax, int &spectrum_size)
  172. {
  173.  
  174. int i;
  175. double density,density1,density_air=0.00120479,attenuation, attenuation_window,attenuation_air,energy, air_thickness;
  176.  
  177. XRayInit();
  178.  
  179. std::vector <double> spectrum_saida ;
  180.  
  181. air_thickness=1.0; /// em cm
  182.  
  183. Emax = energia_maxima ; /// le Emax do arquivo (primeira linha do arquivo)
  184.  
  185. if (Emax == 0 ) G4cout << " Erro na energia maxima " << G4endl;
  186.  
  187. spectrum_size = energias_continuo.size(); /// define o numero de linhas como o tamanho do espectro
  188.  
  189. // G4cout << spectrum_size << G4endl;
  190.  
  191. density=2.33;
  192. density1=1.848;
  193. double calculos(0.0);
  194.  
  195. /// problemas aqui
  196. for (i=0;i<spectrum_size;i++)
  197. {
  198. energy=Emin +(Emax-Emin)*i/spectrum_size;
  199. if (energy>0.0) {
  200. attenuation=CS_Total(Z, energy); /// atenuacao do
  201. attenuation_window=CS_Total(Z_window, energy);
  202. attenuation_air=0.755267*CS_Total(7, energy)+0.231781*CS_Total(8, energy)+ 0.012827*CS_Total(18, energy)+0.000124*CS_Total(6, energy);
  203. /// calcula um novo espectro usando os dados da densidade do ar, atenua��o do ar
  204. /// e largura da janela
  205. calculos = energias_disc.at(i) *(1.0-exp(-attenuation*thickness*density)) *exp(-attenuation_window*window_thick*density1 -attenuation_air*air_thickness*density_air);
  206. spectrum_saida.push_back(calculos);
  207. }
  208. else
  209. {
  210. spectrum_saida.push_back(0.0);
  211. }
  212. }
  213.  
  214. /// imprime o espectro da eficiencia
  215. std::ofstream imprime_arq;
  216. imprime_arq.open("eficiencia");
  217. for (i=0;i<spectrum_saida.size();i++)
  218. {
  219. imprime_arq << (Emin +(Emax-Emin)*i/spectrum_size) << " " << spectrum_saida.at(i) << std::endl;
  220. }
  221. imprime_arq.close();
  222. /// retronarna o espectro da eficiencia
  223. return spectrum_saida;
  224. }
  225.  
  226.  
  227. int Run::detector_convolution_SDD(std::vector<double> spectrum_eff, char *spectrum_out, double Emin, double Emax, int spectrum_size)
  228. {
  229. std::ofstream saida;
  230. int i,j;
  231. double sigma_var;
  232.  
  233.  
  234. double factor;
  235. double E;
  236. int conv_lenght=60; /// ?
  237.  
  238. std::vector <double> conv_nucleus , spectrum_conv ;
  239. //conv_nucleus=(double *)calloc(conv_lenght,sizeof(double));
  240. //spectrum_conv=(double *)calloc(spectrum_size,sizeof(double));
  241.  
  242. conv_nucleus.assign(conv_lenght,0.0);/// zerando
  243. spectrum_conv.assign(spectrum_size,0.0);/// zerando
  244.  
  245. for (i=0;i<energias_continuo.size();i++) spectrum_conv[i]=0.0;
  246.  
  247. for (i=0;i<spectrum_eff.size();i++)
  248. {
  249. sigma_var=sqrt(125.0*125.0-120*120.0+2440.0*i*0.03)/2.3548/30;
  250. factor=1.0/sigma_var/sqrt(2.0*M_PI);
  251. for (j=0;j<conv_lenght;j++) conv_nucleus[j]=factor*exp(-1.0*j*j/(2.0*sigma_var*sigma_var));
  252. for (j=-conv_lenght;j<conv_lenght;j++)
  253. if ( ((i+j)>=0) && ((i+j)<spectrum_size))
  254. spectrum_conv[i+j]+=spectrum_eff.at(i)*conv_nucleus[abs(j)];
  255. }
  256.  
  257. /// escreve o espectro continuo
  258. saida.open(spectrum_out, std::ofstream::out);
  259.  
  260. for (i=0;i<spectrum_size;i++)
  261. {
  262. E = Emin +(Emax-Emin)*i/spectrum_size;
  263. ///cout << E << " " << spectrum_conv[i] << endl;
  264. saida << E << " " << spectrum_conv[i] << std::endl;
  265.  
  266. }
  267. saida.close();
  268.  
  269.  
  270.  
  271.  
  272. return 0;
  273. }
  274.  
  275. void Run::RealizaConvolucao(std::vector<int> energias_disc)
  276. {
  277. ////////////////////////
  278. /// inicio da convolucao
  279. /////////////////////////
  280. char *nome_espectro_saida ;
  281. nome_espectro_saida = new char[8];
  282. nome_espectro_saida="en_cont.out";
  283. int spectrum_size= energias_disc.size();
  284. double Emin=0, Emax=energia_maxima;
  285. std::vector<double> spectrum_eff;
  286. int tamvet = energias_disc.size();
  287.  
  288. ///Material do Detector 14 = Silicio , Espessura do detector , Material da Janela 4= Berilio,
  289. // Espessura da janela = 0 pois o programa já considera a janela de berilio na simulaçao
  290. spectrum_eff= detector_efficiency(14, 0.05, 4, 0.0, Emin, Emax, spectrum_size);
  291. detector_convolution_SDD(spectrum_eff, nome_espectro_saida, Emin, Emax, spectrum_size );
  292.  
  293. /////////
  294. // fim da convolucao arquivo nome_espectro_saida escrito
  295. ////////
  296. }
  297.  
  298.  
  299. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  300.  
  301. void Run::PrintSummary()
  302. {
  303. // compute mean and rms
  304. //
  305. G4int TotNbofEvents = numberOfEvent;
  306. if (TotNbofEvents == 0) return;
  307.  
  308. G4double EnergyBalance = fEnergyDeposit + fEnergyLeak[0] + fEnergyLeak[1];
  309. EnergyBalance /= TotNbofEvents;
  310.  
  311. fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
  312. G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
  313. if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
  314. else rmsEdep = 0.;
  315.  
  316. fTrakLenCharged /= TotNbofEvents; fTrakLenCharged2 /= TotNbofEvents;
  317. G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged*fTrakLenCharged;
  318. if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
  319. else rmsTLCh = 0.;
  320.  
  321. fTrakLenNeutral /= TotNbofEvents; fTrakLenNeutral2 /= TotNbofEvents;
  322. G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral*fTrakLenNeutral;
  323. if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
  324. else rmsTLNe = 0.;
  325.  
  326. fNbStepsCharged /= TotNbofEvents; fNbStepsCharged2 /= TotNbofEvents;
  327. G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged*fNbStepsCharged;
  328. if (rmsStCh>0.) rmsStCh = std::sqrt(rmsTLCh/TotNbofEvents);
  329. else rmsStCh = 0.;
  330.  
  331. fNbStepsNeutral /= TotNbofEvents; fNbStepsNeutral2 /= TotNbofEvents;
  332. G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral*fNbStepsNeutral;
  333. if (rmsStNe>0.) rmsStNe = std::sqrt(rmsTLCh/TotNbofEvents);
  334. else rmsStNe = 0.;
  335.  
  336. G4double Gamma = (G4double)fNbGamma/TotNbofEvents;
  337. G4double Elect = (G4double)fNbElect/TotNbofEvents;
  338. G4double Posit = (G4double)fNbPosit/TotNbofEvents;
  339.  
  340. G4double transmit[2];
  341. transmit[0] = 100.*fTransmit[0]/TotNbofEvents;
  342. transmit[1] = 100.*fTransmit[1]/TotNbofEvents;
  343.  
  344. G4double reflect[2];
  345. reflect[0] = 100.*fReflect[0]/TotNbofEvents;
  346. reflect[1] = 100.*fReflect[1]/TotNbofEvents;
  347.  
  348. G4double rmsMsc = 0., tailMsc = 0.;
  349. if (fMscEntryCentral > 0) {
  350. fMscProjecTheta /= fMscEntryCentral; fMscProjecTheta2 /= fMscEntryCentral;
  351. rmsMsc = fMscProjecTheta2 - fMscProjecTheta*fMscProjecTheta;
  352. if (rmsMsc > 0.) { rmsMsc = std::sqrt(rmsMsc); }
  353. if(fTransmit[1] > 0.0) {
  354. tailMsc = 100.- (100.*fMscEntryCentral)/(2*fTransmit[1]);
  355. }
  356. }
  357.  
  358. fEnergyLeak[0] /= TotNbofEvents; fEnergyLeak2[0] /= TotNbofEvents;
  359. G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0]*fEnergyLeak[0];
  360. if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
  361. else rmsEl0 = 0.;
  362.  
  363. fEnergyLeak[1] /= TotNbofEvents; fEnergyLeak2[1] /= TotNbofEvents;
  364. G4double rmsEl1 = fEnergyLeak2[1] - fEnergyLeak[1]*fEnergyLeak[1];
  365. if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
  366. else rmsEl1 = 0.;
  367.  
  368.  
  369. //Stopping Power from input Table.
  370. //
  371. G4Material* material = fDetector->GetSubstratoMaterial();
  372. G4double length = fDetector->GetSubstratoThickness();
  373. G4double density = material->GetDensity();
  374. G4String partName = fParticle->GetParticleName();
  375.  
  376. G4EmCalculator emCalculator;
  377. G4double dEdxTable = 0., dEdxFull = 0.;
  378. if (fParticle->GetPDGCharge()!= 0.) {
  379. dEdxTable = emCalculator.GetDEDX(fEkin,fParticle,material);
  380. dEdxFull = emCalculator.ComputeTotalDEDX(fEkin,fParticle,material);
  381. }
  382. G4double stopTable = dEdxTable/density;
  383. G4double stopFull = dEdxFull /density;
  384.  
  385. //Stopping Power from simulation.
  386. //
  387. G4double meandEdx = fEnergyDeposit/length;
  388. G4double stopPower = meandEdx/density;
  389.  
  390. G4cout << "\n ======================== run summary ======================\n";
  391.  
  392. G4int prec = G4cout.precision(3);
  393.  
  394. G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
  395. << G4BestUnit(fEkin,"Energy") << " through "
  396. << G4BestUnit(length,"Length") << " of "
  397. << material->GetName() << " (density: "
  398. << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
  399.  
  400. G4cout.precision(4);
  401.  
  402. G4cout << "\n Total energy deposit in absorber per event = "
  403. << G4BestUnit(fEnergyDeposit,"Energy") << " +- "
  404. << G4BestUnit(rmsEdep, "Energy")
  405. << G4endl;
  406.  
  407. G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
  408. << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
  409. << G4endl;
  410.  
  411. G4cout << "\n From formulas :" << G4endl;
  412. G4cout << " restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
  413. << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
  414. << G4endl;
  415.  
  416. G4cout << " full dEdx = " << dEdxFull/(MeV/cm) << " MeV/cm"
  417. << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
  418. << G4endl;
  419.  
  420. G4cout << "\n Leakage : primary = "
  421. << G4BestUnit(fEnergyLeak[0],"Energy") << " +- "
  422. << G4BestUnit(rmsEl0, "Energy")
  423. << " secondaries = "
  424. << G4BestUnit(fEnergyLeak[1],"Energy") << " +- "
  425. << G4BestUnit(rmsEl1, "Energy")
  426. << G4endl;
  427.  
  428. G4cout << " Energy balance : edep + eleak = "
  429. << G4BestUnit(EnergyBalance,"Energy")
  430. << G4endl;
  431.  
  432. G4cout << "\n Total track length (charged) in absorber per event = "
  433. << G4BestUnit(fTrakLenCharged,"Length") << " +- "
  434. << G4BestUnit(rmsTLCh, "Length") << G4endl;
  435.  
  436. G4cout << " Total track length (neutral) in absorber per event = "
  437. << G4BestUnit(fTrakLenNeutral,"Length") << " +- "
  438. << G4BestUnit(rmsTLNe, "Length") << G4endl;
  439.  
  440. G4cout << "\n Number of steps (charged) in absorber per event = "
  441. << fNbStepsCharged << " +- " << rmsStCh << G4endl;
  442.  
  443. G4cout << " Number of steps (neutral) in absorber per event = "
  444. << fNbStepsNeutral << " +- " << rmsStNe << G4endl;
  445.  
  446. G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
  447. << "; electrons = " << Elect
  448. << "; positrons = " << Posit << G4endl;
  449.  
  450. G4cout << "\n Number of events with the primary particle transmitted = "
  451. << transmit[1] << " %" << G4endl;
  452.  
  453. G4cout << " Number of events with at least 1 particle transmitted "
  454. << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
  455.  
  456. G4cout << "\n Number of events with the primary particle reflected = "
  457. << reflect[1] << " %" << G4endl;
  458.  
  459. G4cout << " Number of events with at least 1 particle reflected "
  460. << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
  461.  
  462. // compute width of the Gaussian central part of the MultipleScattering
  463. //
  464. G4cout << "\n MultipleScattering:"
  465. << "\n rms proj angle of transmit primary particle = "
  466. << rmsMsc/mrad << " mrad (central part only)" << G4endl;
  467.  
  468. G4cout << " computed theta0 (Highland formula) = "
  469. << ComputeMscHighland()/mrad << " mrad" << G4endl;
  470.  
  471. G4cout << " central part defined as +- "
  472. << fMscThetaCentral/mrad << " mrad; "
  473. << " Tail ratio = " << tailMsc << " %" << G4endl;
  474.  
  475. // reset default precision
  476. G4cout.precision(prec);
  477. }
  478.  
  479. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  480.  
  481. G4double Run::ComputeMscHighland()
  482. {
  483. //compute the width of the Gaussian central part of the MultipleScattering
  484. //projected angular distribution.
  485. //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
  486.  
  487. G4double t = (fDetector->GetSubstratoThickness())
  488. /(fDetector->GetSubstratoMaterial()->GetRadlen());
  489. if (t < DBL_MIN) return 0.;
  490.  
  491. G4double T = fEkin;
  492. G4double M = fParticle->GetPDGMass();
  493. G4double z = std::abs(fParticle->GetPDGCharge()/eplus);
  494.  
  495. G4double bpc = T*(T+2*M)/(T+M);
  496. G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
  497. return teta0;
  498. }
  499.  
  500. double Run::Calibracao(int _pos , int Emin , int Emax , int spectrum_size )const
  501. {
  502. return (Emin + ((Emax-Emin)*_pos)/ (spectrum_size*1.0)) ;
  503. //return (_pos*6.40/630);
  504. }
  505.  
  506.  
  507. //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  508. //Method to be overwritten by the user for recording events in this (thread-local) run.
  509. // este metodo é chamado em todos os eventos mesmo quando nao há colisao
  510. void Run::RecordEvent(const G4Event* evt)
  511. {
  512.  
  513.  
  514. /// Mostra os resultados
  515. G4HCofThisEvent * HCE = evt->GetHCofThisEvent();
  516. DetectorHitsCollection* HC = 0;
  517. if(HCE)
  518. {
  519. HC = (DetectorHitsCollection*)(HCE->GetHC(hitsCollID));
  520. }
  521. if ( HC )
  522. {
  523. int n_hit = HC->entries();
  524. for ( int i = 0 ; i < n_hit; i++)
  525. {
  526.  
  527. //G4ThreeVector position = (*HC)[i]->GetPos();
  528. G4double energy = (*HC)[i]->GetEdep();
  529. //G4cout << "---- Hit # " << i << G4endl;
  530. //if (energy<10)
  531. {
  532.  
  533. // G4cout<< energy<<G4endl;
  534.  
  535. }
  536. G4double energyant , energydep;
  537. //G4cout<< "Valor real " << energy ;
  538. for (unsigned int j=1;j<(energias_disc.size()-1); j ++)
  539. {
  540. energyant = Calibracao(j-1,0,energia_maxima,energias_disc.size());
  541. energydep = Calibracao(j+1,0,energia_maxima,energias_disc.size());
  542.  
  543. if (energyant<=energy && energydep>=energy)
  544. {
  545. energias_disc[j]+=1;
  546. //G4cout<< " Incluido no canal " << j << " de energia " << Calibracao(j,0,30,energias_disc.size()) << G4endl ;
  547. }
  548. }
  549. // posicao_no_vetor = G4int (ceil(energy/CLHEP::MeV*tamanho)); /// prob
  550. // energias_disc[posicao_no_vetor]+=1;
  551.  
  552. //// MELHORAR !!! aqui o a energia da particula deve
  553. /// acrescentar uma unidade no vetor en_continua
  554. /// deve-se olhar para o vetor en_continua
  555. /// buscar a posicao mais proxima a energia de colisao no detector
  556. /// e adicionar uma unidade a energia de colisao noa eh discreta !!
  557.  
  558. //posicao_no_vetor = G4int (ceil(energy/CLHEP::MeV*tamanho)); /// problemas !
  559. //posicao_no_vetor = G4int (ceil(energy/CLHEP::MeV*tamanho)); /// problemas !
  560. //G4cout << " Position " << position/CLHEP::mm << " [cm] " << " Energy " << energy/CLHEP::MeV << " [keV] " <<G4endl;
  561.  
  562. //energias_disc[posicao_no_vetor]+=1;
  563. }
  564. }
  565.  
  566. G4Run::RecordEvent(evt);
  567.  
  568. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!