[1199] | 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 | // |
---|
| 27 | // $Id: G4hTestStoppingPower.cc,v 1.18 2006/06/29 19:44:46 gunter Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-03-cand-01 $ |
---|
| 29 | // |
---|
| 30 | // ------------------------------------------------------------------- |
---|
| 31 | // GEANT 4 class file --- Copyright CERN 1998 |
---|
| 32 | // CERN Geneva Switzerland |
---|
| 33 | // |
---|
| 34 | // |
---|
| 35 | // File name: G4hTestStoppingPower.cc |
---|
| 36 | // |
---|
| 37 | // Author: Vladimir Ivanchenko |
---|
| 38 | // |
---|
| 39 | // Creation date: 28 July 2000 |
---|
| 40 | // |
---|
| 41 | // Modifications: |
---|
| 42 | // |
---|
| 43 | // ------------------------------------------------------------------- |
---|
| 44 | |
---|
| 45 | #include "globals.hh" |
---|
| 46 | #include "G4ios.hh" |
---|
| 47 | #include "G4UnitsTable.hh" |
---|
| 48 | #include <fstream> |
---|
| 49 | #include <iomanip> |
---|
| 50 | |
---|
| 51 | #include "G4Material.hh" |
---|
| 52 | #include "G4Element.hh" |
---|
| 53 | #include "G4ProcessManager.hh" |
---|
| 54 | #include "G4hLowEnergyIonisation.hh" |
---|
| 55 | #include "G4hBetheBlochModel.hh" |
---|
| 56 | #include "G4hParametrisedLossModel.hh" |
---|
| 57 | #include "G4hNuclearStoppingModel.hh" |
---|
| 58 | #include "G4QAOLowEnergyLoss.hh" |
---|
| 59 | #include "G4VParticleChange.hh" |
---|
| 60 | #include "G4ParticleChange.hh" |
---|
| 61 | #include "G4DynamicParticle.hh" |
---|
| 62 | |
---|
| 63 | #include "G4hZiegler1977p.hh" |
---|
| 64 | #include "G4hZiegler1985p.hh" |
---|
| 65 | #include "G4hZiegler1977He.hh" |
---|
| 66 | #include "G4hICRU49p.hh" |
---|
| 67 | #include "G4hICRU49He.hh" |
---|
| 68 | #include "G4hSRIM2000p.hh" |
---|
| 69 | #include "G4hSRIM2003p.hh" |
---|
| 70 | |
---|
| 71 | #include "G4hZiegler1977Nuclear.hh" |
---|
| 72 | #include "G4hZiegler1985Nuclear.hh" |
---|
| 73 | //#include "G4hMollereNuclear.hh" // exist no more |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | #include "G4Box.hh" |
---|
| 77 | #include "G4PVPlacement.hh" |
---|
| 78 | |
---|
| 79 | #include "G4Step.hh" |
---|
| 80 | #include "G4GRSVolume.hh" |
---|
| 81 | |
---|
| 82 | #include "G4ParticleTypes.hh" |
---|
| 83 | #include "G4ParticleTable.hh" |
---|
| 84 | #include "G4ParticleDefinition.hh" |
---|
| 85 | #include "G4ParticleWithCuts.hh" |
---|
| 86 | #include "G4Ions.hh" |
---|
| 87 | |
---|
| 88 | // New Histogramming (from AIDA and Anaphe): |
---|
| 89 | #include <memory> // for the auto_ptr(T> |
---|
| 90 | |
---|
| 91 | #include "AIDA/AIDA.h" |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | #include "hTest/include/G4IonC12.hh" |
---|
| 95 | #include "hTest/include/G4IonAr40.hh" |
---|
| 96 | |
---|
| 97 | #include "G4Timer.hh" |
---|
| 98 | |
---|
| 99 | int main() |
---|
| 100 | { |
---|
| 101 | // ---- HBOOK initialization |
---|
| 102 | |
---|
| 103 | G4String hFile = "htest.hbook"; |
---|
| 104 | |
---|
| 105 | //--------- Materials definition --------- |
---|
| 106 | |
---|
| 107 | G4Material* Be = new G4Material("Beryllium", 4., 9.01*g/mole, 1.848*g/cm3); |
---|
| 108 | G4Material* Graphite = new G4Material("Graphite",6.,12.0*g/mole,2.265*g/cm3); |
---|
| 109 | G4Material* Al = new G4Material("Aluminium", 13., 26.98*g/mole, 2.7 *g/cm3); |
---|
| 110 | G4Material* Si = new G4Material("Silicon", 14., 28.055*g/mole, 2.33*g/cm3); |
---|
| 111 | G4Material* LAr = new G4Material("LArgon", 18., 39.95*g/mole, 1.393*g/cm3); |
---|
| 112 | G4Material* Ti = new G4Material("Titan", 22., 47.867*g/mole, 4.54*g/cm3); |
---|
| 113 | G4Material* Fe = new G4Material("Iron", 26., 55.85*g/mole, 7.87*g/cm3); |
---|
| 114 | G4Material* Cu = new G4Material("Copper", 29., 63.55*g/mole, 8.96*g/cm3); |
---|
| 115 | G4Material* W = new G4Material("Tungsten",74., 183.85*g/mole, 19.30*g/cm3); |
---|
| 116 | G4Material* Pb = new G4Material("Lead", 82., 207.19*g/mole, 11.35*g/cm3); |
---|
| 117 | G4Material* U = new G4Material("Uranium", 92., 238.03*g/mole, 18.95*g/cm3); |
---|
| 118 | |
---|
| 119 | G4Element* H = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole); |
---|
| 120 | G4Element* O = new G4Element ("Oxygen" , "O", 8. , 16.00*g/mole); |
---|
| 121 | G4Element* C = new G4Element ("Carbon" , "C", 6. , 12.00*g/mole); |
---|
| 122 | G4Element* Cs = new G4Element ("Cesium" , "Cs", 55. , 132.905*g/mole); |
---|
| 123 | G4Element* I = new G4Element ("Iodide" , "I", 53. , 126.9044*g/mole); |
---|
| 124 | |
---|
| 125 | // G4Material* water = new G4Material ("Water" ,"H_2O", 1.*g/cm3, 2); |
---|
| 126 | G4Material* water = new G4Material ("Water" , 1.*g/cm3, 2); |
---|
| 127 | water->AddElement(H,2); |
---|
| 128 | water->AddElement(O,1); |
---|
| 129 | |
---|
| 130 | // G4Material* ethane = new G4Material ("Ethane" ,"C_2H_6", 0.4241*g/cm3, 2); |
---|
| 131 | G4Material* ethane = new G4Material ("Ethane" , 0.4241*g/cm3, 2); |
---|
| 132 | ethane->AddElement(H,6); |
---|
| 133 | ethane->AddElement(C,2); |
---|
| 134 | |
---|
| 135 | // G4Material* csi = new G4Material ("CsI" , "CsI", 4.53*g/cm3, 2); |
---|
| 136 | G4Material* csi = new G4Material ("CsI" , 4.53*g/cm3, 2); |
---|
| 137 | csi->AddElement(Cs,1); |
---|
| 138 | csi->AddElement(I,1); |
---|
| 139 | |
---|
| 140 | G4Material* el[93]; |
---|
| 141 | G4double z; |
---|
| 142 | G4int j; |
---|
| 143 | |
---|
| 144 | for ( j=1; j<93; j++) |
---|
| 145 | { |
---|
| 146 | z = G4double(j); |
---|
| 147 | el[j] = new G4Material ("atom", z , 2.0*z*g/mole, 10.0 *g/cm3); |
---|
| 148 | } |
---|
| 149 | |
---|
| 150 | const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); |
---|
| 151 | |
---|
| 152 | // create table |
---|
| 153 | G4int numOfMaterials = G4Material::GetNumberOfMaterials(); |
---|
| 154 | |
---|
| 155 | G4double dimx = 1*mm, dimy = 1*mm, dimz = 1*mm; |
---|
| 156 | G4int imat = 0; |
---|
| 157 | G4cout << "Available materials are: " << G4endl; |
---|
| 158 | for (imat = 0; imat < numOfMaterials; imat++) { |
---|
| 159 | G4cout << imat << ") " << (*theMaterialTable)[imat]->GetName() << G4endl; |
---|
| 160 | } |
---|
| 161 | |
---|
| 162 | |
---|
| 163 | // Geometry definitions |
---|
| 164 | G4Box* theFrame = new G4Box ("Frame",dimx, dimy, dimz); |
---|
| 165 | |
---|
| 166 | G4LogicalVolume* LogicalFrame = new G4LogicalVolume(theFrame, |
---|
| 167 | (*theMaterialTable)[0], |
---|
| 168 | "LFrame", 0, 0, 0); |
---|
| 169 | |
---|
| 170 | G4PVPlacement* PhysicalFrame = new G4PVPlacement(0,G4ThreeVector(), |
---|
| 171 | "PFrame",LogicalFrame,0,false,0); |
---|
| 172 | |
---|
| 173 | G4cout << PhysicalFrame << G4endl; |
---|
| 174 | |
---|
| 175 | //--------- Particle definition --------- |
---|
| 176 | |
---|
| 177 | G4ParticleDefinition* gamma = G4Gamma::Gamma(); |
---|
| 178 | G4ParticleDefinition* electron = G4Electron::Electron(); |
---|
| 179 | G4ParticleDefinition* proton = G4Proton::Proton(); |
---|
| 180 | G4ParticleDefinition* antiproton =G4AntiProton::AntiProton(); |
---|
| 181 | G4ParticleDefinition* deuteron = G4Deuteron::Deuteron(); |
---|
| 182 | G4ParticleDefinition* alpha = G4Alpha::Alpha(); |
---|
| 183 | |
---|
| 184 | G4Ions* iC12 = new G4Ions::G4Ions( |
---|
| 185 | "IonC12", 11.14945*GeV, 0.0*MeV, +6.0*eplus, |
---|
| 186 | 0, +1, 0, |
---|
| 187 | 0, 0, 0, |
---|
| 188 | "static_nucleus", 0, +12, 0, |
---|
| 189 | true, -1.0, NULL); |
---|
| 190 | |
---|
| 191 | G4Ions* iAr40 = new G4Ions::G4Ions( |
---|
| 192 | "IonAr40", 37.291*GeV, 0.0*MeV, +18.0*eplus, |
---|
| 193 | 0, +1, 0, |
---|
| 194 | 0, 0, 0, |
---|
| 195 | "static_nucleus", 0, +40, 0, |
---|
| 196 | true, -1.0, NULL); |
---|
| 197 | |
---|
| 198 | G4Ions* iFe56 = new G4Ions::G4Ions( |
---|
| 199 | "IonFe56", 52.0308*GeV, 0.0*MeV, +26.0*eplus, |
---|
| 200 | 0, +1, 0, |
---|
| 201 | 0, 0, 0, |
---|
| 202 | "static_nucleus", 0, +56, 0, |
---|
| 203 | true, -1.0, NULL); |
---|
| 204 | |
---|
| 205 | G4ParticleDefinition* ionC12 = iC12->IonsDefinition(); |
---|
| 206 | G4ParticleDefinition* ionAr40 = iAr40->IonsDefinition(); |
---|
| 207 | G4ParticleDefinition* ionFe56 = iFe56->IonsDefinition(); |
---|
| 208 | |
---|
| 209 | G4ParticleDefinition* part[7]; |
---|
| 210 | part[0] = proton; |
---|
| 211 | part[1] = antiproton; |
---|
| 212 | part[2] = deuteron; |
---|
| 213 | part[3] = alpha; |
---|
| 214 | part[4] = ionC12; |
---|
| 215 | part[5] = ionAr40; |
---|
| 216 | part[6] = ionFe56; |
---|
| 217 | |
---|
| 218 | G4double ecut = 1000.0*mm; |
---|
| 219 | G4double pcut = 0.0001*mm; |
---|
| 220 | gamma->SetCuts(ecut); |
---|
| 221 | electron->SetCuts(ecut); |
---|
| 222 | G4cout << "Cuts are following: cutElectron = " << ecut |
---|
| 223 | << " mm; cutProton = " << pcut << " mm" << G4endl; |
---|
| 224 | |
---|
| 225 | //--------- Ionisation processes definition and build physics table -- |
---|
| 226 | |
---|
| 227 | // Define models for parametrisation of electronic energy losses |
---|
| 228 | // G4VLowEnergyModel* theBetheBlochModel = |
---|
| 229 | // new G4hBetheBlochModel("Bethe-Bloch") ; |
---|
| 230 | // theProtonModel = new G4hParametrisedLossModel(theProtonTable) ; |
---|
| 231 | // theAntiProtonModel = new G4QAOLowEnergyLoss(theAntiProtonTable) ; |
---|
| 232 | |
---|
| 233 | G4hNuclearStoppingModel* theNuclearStoppingModel = |
---|
| 234 | // new G4hNuclearStoppingModel("Ziegler1977") ; |
---|
| 235 | // new G4hNuclearStoppingModel("Ziegler1985") ; |
---|
| 236 | new G4hNuclearStoppingModel("ICRU_R49") ; |
---|
| 237 | |
---|
| 238 | G4VLowEnergyModel* theIonEffChargeModel = |
---|
| 239 | new G4hIonEffChargeSquare("Ziegler1988") ; |
---|
| 240 | |
---|
| 241 | |
---|
| 242 | G4hLowEnergyIonisation* hIon[8]; |
---|
| 243 | G4ProcessManager* theProcessManager[8]; |
---|
| 244 | G4int i; |
---|
| 245 | |
---|
| 246 | G4cout << "Define processes!" << G4endl; |
---|
| 247 | |
---|
| 248 | for( i=0; i<7; i++) { |
---|
| 249 | G4cout << "Ionisation process for particle " << i << " " << part[i]->GetParticleName() << G4endl; |
---|
| 250 | theProcessManager[i] = new G4ProcessManager(part[i]); |
---|
| 251 | part[i]->SetProcessManager(theProcessManager[i]); |
---|
| 252 | hIon[i] = new G4hLowEnergyIonisation(); |
---|
| 253 | hIon[i]->SetEnlossFluc(false) ; |
---|
| 254 | |
---|
| 255 | // hIon[i]->SetBarkasOff(); |
---|
| 256 | // hIon[i]->SetNuclearStoppingOff(); |
---|
| 257 | // hIon[i]->SetStoppingPowerTableName("ICRU_R49p"); |
---|
| 258 | |
---|
| 259 | theProcessManager[i]->AddProcess(hIon[i]); |
---|
| 260 | hIon[i]->BuildPhysicsTable(*part[i]); |
---|
| 261 | } |
---|
| 262 | |
---|
| 263 | |
---|
| 264 | |
---|
| 265 | //----------- Histogram ------------------------------------- |
---|
| 266 | |
---|
| 267 | G4cout << "Fill Hbook!" << G4endl; |
---|
| 268 | |
---|
| 269 | // Creating the analysis factory |
---|
| 270 | std::auto_ptr< AIDA::IAnalysisFactory > af( AIDA_createAnalysisFactory() ); |
---|
| 271 | |
---|
| 272 | // Creating the tree factory |
---|
| 273 | std::auto_ptr< AIDA::ITreeFactory > tf( af->createTreeFactory() ); |
---|
| 274 | |
---|
| 275 | // Creating a tree mapped to a new hbook file. |
---|
| 276 | std::auto_ptr< AIDA::ITree > tree( tf->create( hFile,"hbook", false,false) ); |
---|
| 277 | std::cout << "Tree store : " << tree->storeName() << std::endl; |
---|
| 278 | |
---|
| 279 | // Creating a histogram factory, whose histograms will be handled by the tree |
---|
| 280 | std::auto_ptr< AIDA::IHistogramFactory > hf( af->createHistogramFactory( *tree ) ); |
---|
| 281 | |
---|
| 282 | // G4Material* material ; |
---|
| 283 | |
---|
| 284 | G4double minE = 1.0*eV, maxE = 10000.0*MeV, s; |
---|
| 285 | const G4int num = 200; |
---|
| 286 | G4double tkin = 0.0; |
---|
| 287 | s = (std::log10(maxE)-std::log10(minE))/num; |
---|
| 288 | |
---|
| 289 | AIDA::IHistogram1D* h[71] ; |
---|
| 290 | |
---|
| 291 | // Test on Stopping Powers for all elements |
---|
| 292 | |
---|
| 293 | h[1] = hf->createHistogram1D("1","p 40 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
| 294 | ,92,0.5,92.5) ; |
---|
| 295 | h[2] = hf->createHistogram1D("2","p 100 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
| 296 | ,92,0.5,92.5) ; |
---|
| 297 | h[3] = hf->createHistogram1D("3","p 400 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
| 298 | ,92,0.5,92.5) ; |
---|
| 299 | |
---|
| 300 | h[4] = hf->createHistogram1D("4","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
| 301 | ,92,0.5,92.5) ; |
---|
| 302 | h[5] = hf->createHistogram1D("5","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
| 303 | ,92,0.5,92.5) ; |
---|
| 304 | |
---|
| 305 | h[6] = hf->createHistogram1D("6","p 40 keV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
| 306 | ,92,0.5,92.5) ; |
---|
| 307 | h[7] = hf->createHistogram1D("7","p 100 keV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
| 308 | ,92,0.5,92.5) ; |
---|
| 309 | h[8] = hf->createHistogram1D("8","p 400 keV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
| 310 | ,92,0.5,92.5) ; |
---|
| 311 | h[9] = hf->createHistogram1D("9","p 1 MeV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
| 312 | ,92,0.5,92.5) ; |
---|
| 313 | h[10] = hf->createHistogram1D("10","p 4 MeV (keV*cm2/10^15!atoms) ICRU_49p" |
---|
| 314 | ,92,0.5,92.5) ; |
---|
| 315 | |
---|
| 316 | h[11] = hf->createHistogram1D("11","p 40 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 317 | ,92,0.5,92.5) ; |
---|
| 318 | h[12] = hf->createHistogram1D("12","p 100 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 319 | ,92,0.5,92.5) ; |
---|
| 320 | h[13] = hf->createHistogram1D("13","p 400 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 321 | ,92,0.5,92.5) ; |
---|
| 322 | h[14] = hf->createHistogram1D("14","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 323 | ,92,0.5,92.5) ; |
---|
| 324 | h[15] = hf->createHistogram1D("15","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 325 | ,92,0.5,92.5) ; |
---|
| 326 | |
---|
| 327 | h[16] = hf->createHistogram1D("16","He 10 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 328 | ,92,0.5,92.5) ; |
---|
| 329 | h[17] = hf->createHistogram1D("17","He 40 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 330 | ,92,0.5,92.5) ; |
---|
| 331 | h[18] = hf->createHistogram1D("18","He 100 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 332 | ,92,0.5,92.5) ; |
---|
| 333 | |
---|
| 334 | h[19] = hf->createHistogram1D("19","He 400 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 335 | ,92,0.5,92.5) ; |
---|
| 336 | h[20] = hf->createHistogram1D("20","He 1 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 337 | ,92,0.5,92.5) ; |
---|
| 338 | h[21] = hf->createHistogram1D("21","He 4 MeV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 339 | ,92,0.5,92.5) ; |
---|
| 340 | h[22] = hf->createHistogram1D("22","He 10 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
| 341 | ,92,0.5,92.5) ; |
---|
| 342 | h[23] = hf->createHistogram1D("23","He 40 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
| 343 | ,92,0.5,92.5) ; |
---|
| 344 | h[24] = hf->createHistogram1D("24","He 100 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
| 345 | ,92,0.5,92.5) ; |
---|
| 346 | h[25] = hf->createHistogram1D("25","He 400 keV (keV*cm2/10^15!atoms) ICRU49He" |
---|
| 347 | ,92,0.5,92.5) ; |
---|
| 348 | h[26] = hf->createHistogram1D("26","He 1 MeV (keV*cm2/10^15!atoms) ICRU49He" |
---|
| 349 | ,92,0.5,92.5) ; |
---|
| 350 | h[27] = hf->createHistogram1D("27","He 4 MeV (keV*cm2/10^15!atoms) ICRU49He" |
---|
| 351 | ,92,0.5,92.5) ; |
---|
| 352 | |
---|
| 353 | h[28]= hf->createHistogram1D("28","p in C (MeV/mm) ICRU49p" |
---|
| 354 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 355 | h[29]= hf->createHistogram1D("29","p in Al (MeV/mm) ICRU49p" |
---|
| 356 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 357 | h[30]= hf->createHistogram1D("30","p in Si (MeV/mm) ICRU49p" |
---|
| 358 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 359 | h[31]= hf->createHistogram1D("31","p in Cu (MeV/mm) ICRU49p" |
---|
| 360 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 361 | h[32]= hf->createHistogram1D("32","p in Fe (MeV/mm) ICRU49p" |
---|
| 362 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 363 | h[33]= hf->createHistogram1D("33","p in Pb (MeV/mm) ICRU49p" |
---|
| 364 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 365 | h[34]= hf->createHistogram1D("34","p in C2H6 (MeV/mm) ICRU49p" |
---|
| 366 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 367 | h[35]= hf->createHistogram1D("35","p in H2O (MeV/mm) ICRU49p" |
---|
| 368 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 369 | h[36]= hf->createHistogram1D("36","p in lAr (MeV/mm) ICRU49p" |
---|
| 370 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 371 | h[37]= hf->createHistogram1D("37","p in CsI (MeV/mm) ICRU49p" |
---|
| 372 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 373 | |
---|
| 374 | |
---|
| 375 | h[38]= hf->createHistogram1D("38","p in C (MeV/mm)Ziegler1985p" |
---|
| 376 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 377 | h[39]= hf->createHistogram1D("39","p in Al (MeV/mm)Ziegler1985p" |
---|
| 378 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 379 | h[40]= hf->createHistogram1D("40","p in Si (MeV/mm)Ziegler1985p" |
---|
| 380 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 381 | h[41]= hf->createHistogram1D("41","p in Cu (MeV/mm)Ziegler1985p" |
---|
| 382 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 383 | h[42]= hf->createHistogram1D("42","p in Fe (MeV/mm)Ziegler1985p" |
---|
| 384 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 385 | h[43]= hf->createHistogram1D("43","p in Pb (MeV/mm)Ziegler1985p" |
---|
| 386 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 387 | h[44]= hf->createHistogram1D("44","p in C2H6 (MeV/mm)Ziegler1985p" |
---|
| 388 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 389 | h[45]= hf->createHistogram1D("45","p in H2O (MeV/mm)Ziegler1985p" |
---|
| 390 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 391 | h[46]= hf->createHistogram1D("46","p in lAr (MeV/mm)Ziegler1985p" |
---|
| 392 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 393 | h[47]= hf->createHistogram1D("47","p in CsI (MeV/mm)Ziegler1985p" |
---|
| 394 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 395 | |
---|
| 396 | |
---|
| 397 | h[48]= hf->createHistogram1D("48","He effective charge for Cu" |
---|
| 398 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 399 | h[49]= hf->createHistogram1D("49","C12 effective charge in Cu" |
---|
| 400 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 401 | |
---|
| 402 | h[50]= hf->createHistogram1D("50","He in Al (MeV/(mg/cm2)) ICRU49p" |
---|
| 403 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 404 | h[51]= hf->createHistogram1D("51","C12 in Al (MeV/(mg/cm2)) ICRU49p" |
---|
| 405 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 406 | h[52]= hf->createHistogram1D("52","Ar40 in Al (MeV/(mg/cm2)) ICRU49p" |
---|
| 407 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 408 | |
---|
| 409 | // Table with the data |
---|
| 410 | h[53] = hf->createHistogram1D("53","Data p 40 keV (keV*cm2/10^15!atoms) Ziegler1977p" |
---|
| 411 | ,92,0.5,92.5) ; |
---|
| 412 | h[54] = hf->createHistogram1D("54","Data He 40 keV (keV*cm2/10^15!atoms) Ziegler1977He" |
---|
| 413 | ,92,0.5,92.5) ; |
---|
| 414 | |
---|
| 415 | h[55] = hf->createHistogram1D("55","p 40 keV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
| 416 | ,92,0.5,92.5) ; |
---|
| 417 | h[56] = hf->createHistogram1D("56","p 100 keV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
| 418 | ,92,0.5,92.5) ; |
---|
| 419 | h[57] = hf->createHistogram1D("57","p 400 keV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
| 420 | ,92,0.5,92.5) ; |
---|
| 421 | h[58] = hf->createHistogram1D("58","p 1 MeV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
| 422 | ,92,0.5,92.5) ; |
---|
| 423 | h[59] = hf->createHistogram1D("59","p 4 MeV (keV*cm2/10^15!atoms) Ziegler1985p" |
---|
| 424 | ,92,0.5,92.5) ; |
---|
| 425 | // Histo for Antiproton |
---|
| 426 | |
---|
| 427 | h[60]= hf->createHistogram1D("60","pbar in C (MeV/mm) QAOLoss" |
---|
| 428 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 429 | h[61]= hf->createHistogram1D("61","pbar in Al (MeV/mm) QAOLoss" |
---|
| 430 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 431 | h[62]= hf->createHistogram1D("62","pbar in Si (MeV/mm) QAOLoss" |
---|
| 432 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 433 | h[63]= hf->createHistogram1D("63","pbar in Cu (MeV/mm) QAOLoss" |
---|
| 434 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 435 | h[64]= hf->createHistogram1D("64","pbar in Fe (MeV/mm) QAOLoss" |
---|
| 436 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 437 | h[65]= hf->createHistogram1D("65","pbar in Pb (MeV/mm) QAOLoss" |
---|
| 438 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 439 | h[66]= hf->createHistogram1D("66","pbar in C2H6 (MeV/mm) QAOLoss" |
---|
| 440 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 441 | h[67]= hf->createHistogram1D("67","pbar in H2O (MeV/mm) QAOLoss" |
---|
| 442 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 443 | h[68]= hf->createHistogram1D("68","pbar in lAr (MeV/mm) QAOLoss" |
---|
| 444 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 445 | h[69]= hf->createHistogram1D("69","pbar in CsI (MeV/mm) QAOLoss" |
---|
| 446 | ,num,std::log10(minE),std::log10(maxE)) ; |
---|
| 447 | |
---|
| 448 | h[70] = hf->createHistogram1D("70","p 6.5 MeV (keV*cm2/10^15!atoms) Ziegler77p" |
---|
| 449 | ,92,0.5,92.5) ; |
---|
| 450 | |
---|
| 451 | |
---|
| 452 | // G4VhElectronicStoppingPower* Z77p = new G4hZiegler1977p() ; |
---|
| 453 | // G4VhElectronicStoppingPower* Z85p = new G4hZiegler1985p() ; |
---|
| 454 | G4VhElectronicStoppingPower* Z77p = new G4hSRIM2000p() ; |
---|
| 455 | G4VhElectronicStoppingPower* Z85p = new G4hSRIM2003p() ; |
---|
| 456 | G4VhElectronicStoppingPower* Z77He = new G4hZiegler1977He() ; |
---|
| 457 | G4VhElectronicStoppingPower* I49p = new G4hICRU49p() ; |
---|
| 458 | G4VhElectronicStoppingPower* I49He = new G4hICRU49He() ; |
---|
| 459 | |
---|
| 460 | G4double de; |
---|
| 461 | |
---|
| 462 | // Ziegler1977p |
---|
| 463 | G4double tau = 40.0*keV ; |
---|
| 464 | G4double q2; |
---|
| 465 | G4double rateMass=4.0026/1.007276; |
---|
| 466 | |
---|
| 467 | for ( j=1; j<93; j++) |
---|
| 468 | { |
---|
| 469 | z = G4double(j); |
---|
| 470 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
| 471 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
| 472 | h[1]->fill(z,de) ; |
---|
| 473 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
| 474 | h[55]->fill(z,de) ; |
---|
| 475 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
| 476 | h[6]->fill(z,de) ; |
---|
| 477 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
| 478 | h[11]->fill(z,de) ; |
---|
| 479 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 480 | h[17]->fill(z,de) ; |
---|
| 481 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 482 | h[23]->fill(z,de) ; |
---|
| 483 | } |
---|
| 484 | tau = 100.0*keV ; |
---|
| 485 | for ( j=1; j<93; j++) |
---|
| 486 | { |
---|
| 487 | z = G4double(j); |
---|
| 488 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
| 489 | h[2]->fill(z,de) ; |
---|
| 490 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
| 491 | h[56]->fill(z,de) ; |
---|
| 492 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
| 493 | h[7]->fill(z,de) ; |
---|
| 494 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
| 495 | h[12]->fill(z,de) ; |
---|
| 496 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
| 497 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 498 | h[18]->fill(z,de) ; |
---|
| 499 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 500 | h[24]->fill(z,de) ; |
---|
| 501 | } |
---|
| 502 | |
---|
| 503 | tau = 400.0*keV ; |
---|
| 504 | for ( j=1; j<93; j++) |
---|
| 505 | { |
---|
| 506 | z = G4double(j); |
---|
| 507 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
| 508 | h[3]->fill(z,de) ; |
---|
| 509 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
| 510 | h[57]->fill(z,de) ; |
---|
| 511 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
| 512 | h[8]->fill(z,de) ; |
---|
| 513 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
| 514 | h[13]->fill(z,de) ; |
---|
| 515 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
| 516 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 517 | h[19]->fill(z,de) ; |
---|
| 518 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 519 | h[25]->fill(z,de) ; |
---|
| 520 | } |
---|
| 521 | tau = 1000.0*keV ; |
---|
| 522 | for ( j=1; j<93; j++) |
---|
| 523 | { |
---|
| 524 | z = G4double(j); |
---|
| 525 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
| 526 | h[4]->fill(z,de) ; |
---|
| 527 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
| 528 | h[58]->fill(z,de) ; |
---|
| 529 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
| 530 | h[9]->fill(z,de) ; |
---|
| 531 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
| 532 | h[14]->fill(z,de) ; |
---|
| 533 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
| 534 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 535 | h[20]->fill(z,de) ; |
---|
| 536 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 537 | h[26]->fill(z,de) ; |
---|
| 538 | } |
---|
| 539 | tau = 4000*keV ; |
---|
| 540 | for ( j=1; j<93; j++) |
---|
| 541 | { |
---|
| 542 | z = G4double(j); |
---|
| 543 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
| 544 | h[5]->fill(z,de) ; |
---|
| 545 | de = Z85p->ElectronicStoppingPower(z, tau) ; |
---|
| 546 | h[59]->fill(z,de) ; |
---|
| 547 | de = I49p->ElectronicStoppingPower(z, tau) ; |
---|
| 548 | h[10]->fill(z,de) ; |
---|
| 549 | de = Z77He->ElectronicStoppingPower(z, tau) ; |
---|
| 550 | h[15]->fill(z,de) ; |
---|
| 551 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
| 552 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 553 | h[21]->fill(z,de) ; |
---|
| 554 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 555 | h[27]->fill(z,de) ; |
---|
| 556 | } |
---|
| 557 | tau = 6.5*MeV ; |
---|
| 558 | for ( j=1; j<93; j++) |
---|
| 559 | { |
---|
| 560 | z = G4double(j); |
---|
| 561 | de = Z77p->ElectronicStoppingPower(z, tau) ; |
---|
| 562 | h[70]->fill(z,de) ; |
---|
| 563 | } |
---|
| 564 | tau = 10.0*keV ; |
---|
| 565 | for ( j=1; j<93; j++) |
---|
| 566 | { |
---|
| 567 | z = G4double(j); |
---|
| 568 | q2 = theIonEffChargeModel->TheValue(part[3],el[j],tkin); |
---|
| 569 | de = (Z77He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 570 | h[16]->fill(z,de) ; |
---|
| 571 | de = (I49He->ElectronicStoppingPower(z, tau/rateMass))*q2 ; |
---|
| 572 | h[22]->fill(z,de) ; |
---|
| 573 | } |
---|
| 574 | |
---|
| 575 | for ( j=1; j<93; j++) |
---|
| 576 | { |
---|
| 577 | static G4double p40[92] = { |
---|
| 578 | 6.22, 6.55, 7.61, 10.4, 12.9, 13.9, 15.9, 14.6, 11.7, 11.1, |
---|
| 579 | 14.2, 20.6, 20.7, 22.3, 18.1, 19.3, 27.5, 30.5, 27.9, 29.8, |
---|
| 580 | 28.3, 26.7, 25.1, 22.5, 19.8, 20.1, 18.0, 20.1, 20.4, 23.4, |
---|
| 581 | 27.5, 29.0, 29.3, 31.0, 31.1, 34.8, 31.5, 35.0, 35.5, 37.3, |
---|
| 582 | 38.3, 35.9, 38.0, 34.4, 33.4, 29.8, 30.9, 32.7, 34.9, 36.0, |
---|
| 583 | 40.9, 39.1, 43.1, 45.8, 40.8, 44.1, 44.9, 42.0, 41.0, 40.0, |
---|
| 584 | 39.0, 38.0, 37.1, 38.2, 35.3, 31.5, 29.9, 29.1, 28.3, 27.5, |
---|
| 585 | 28.1, 28.9, 27.3, 26.3, 29.9, 29.1, 28.4, 25.8, 28.0, 24.9, |
---|
| 586 | 27.3, 30.7, 34.3, 35.4, 35.6, 35.5, 39.8, 42.9, 43.7, 44.1, |
---|
| 587 | 42.4, 41.8} ; |
---|
| 588 | h[53]->fill(double(j),p40[j-1]) ; |
---|
| 589 | |
---|
| 590 | static G4double he40[92] = { |
---|
| 591 | 11.8, 16.7, 21.9, 24.1, 34.7, 36.1, 45.5, 46.1, 45.8, 44.9, |
---|
| 592 | 46.9, 51.7, 54.9, 63.2, 63.8, 67.4, 79.8, 80.8, 78.0, 83.2, |
---|
| 593 | 85.4, 84.2, 90.6, 85.7, 84.1, 86.4, 82.0, 77.7, 74.4, 75.1, |
---|
| 594 | 75.6, 80.7, 84.9, 87.0, 88.6, 103.0, 103.0, 113.0, 116.0, 123.0, |
---|
| 595 | 120.0, 115.0, 118.0, 115.0, 113.0, 112.0, 111.0, 110.0, 116.0, 117.0, |
---|
| 596 | 119.0, 123.0, 122.0, 139.0, 138.0, 143.0, 152.0, 141.0, 139.0, 137.0, |
---|
| 597 | 135.0, 134.0, 130.0, 137.0, 129.0, 128.0, 124.0, 125.0, 120.0, 118.0, |
---|
| 598 | 119.0, 120.0, 121.0, 122.0, 125.0, 127.0, 128.0, 120.0, 125.0, 128.0, |
---|
| 599 | 136.0, 138.0, 144.0, 148.0, 151.0, 162.0, 166.0, 177.0, 179.0, 183.0, |
---|
| 600 | 177.0, 176.0 } ; |
---|
| 601 | h[54]->fill(double(j),he40[j-1]) ; |
---|
| 602 | } |
---|
| 603 | |
---|
| 604 | G4cout << "Stopping Power histograms are filled!" << G4endl; |
---|
| 605 | |
---|
| 606 | // dedx |
---|
| 607 | for (j = 0 ; j < num-1 ; j++) { |
---|
| 608 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
| 609 | de = hIon[0]->ComputeDEDX(part[0],Graphite,tkin) ; |
---|
| 610 | h[28]->fill(std::log10(tkin),de) ; |
---|
| 611 | de = hIon[0]->ComputeDEDX(part[0],Al,tkin) ; |
---|
| 612 | h[29]->fill(std::log10(tkin),de) ; |
---|
| 613 | de = hIon[0]->ComputeDEDX(part[0],Si,tkin) ; |
---|
| 614 | h[30]->fill(std::log10(tkin),de) ; |
---|
| 615 | de = hIon[0]->ComputeDEDX(part[0],Cu,tkin) ; |
---|
| 616 | h[31]->fill(std::log10(tkin),de) ; |
---|
| 617 | de = hIon[0]->ComputeDEDX(part[0],Fe,tkin) ; |
---|
| 618 | h[32]->fill(std::log10(tkin),de) ; |
---|
| 619 | de = hIon[0]->ComputeDEDX(part[0],Pb,tkin) ; |
---|
| 620 | h[33]->fill(std::log10(tkin),de) ; |
---|
| 621 | de = hIon[0]->ComputeDEDX(part[0],ethane,tkin) ; |
---|
| 622 | // G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl; |
---|
| 623 | h[34]->fill(std::log10(tkin),de) ; |
---|
| 624 | de = hIon[0]->ComputeDEDX(part[0],water,tkin) ; |
---|
| 625 | // de += theNuclearStoppingModel->TheValue(part[1],water,tkin) ; |
---|
| 626 | // G4cout << "water : E = " << tkin << "; dedx = " << de << G4endl; |
---|
| 627 | h[35]->fill(std::log10(tkin),de) ; |
---|
| 628 | de = hIon[0]->ComputeDEDX(part[0],LAr,tkin) ; |
---|
| 629 | h[36]->fill(std::log10(tkin),de) ; |
---|
| 630 | de = hIon[0]->ComputeDEDX(part[0],csi,tkin) ; |
---|
| 631 | h[37]->fill(std::log10(tkin),de) ; |
---|
| 632 | } |
---|
| 633 | |
---|
| 634 | G4cout << "Proton's dEdx histograms are filled!" << G4endl; |
---|
| 635 | |
---|
| 636 | for (j = 0 ; j < num-1 ; j++) { |
---|
| 637 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
| 638 | de = hIon[0]->ComputeDEDX(part[1],Graphite,tkin) ; |
---|
| 639 | h[60]->fill(std::log10(tkin),de) ; |
---|
| 640 | de = hIon[0]->ComputeDEDX(part[1],Al,tkin) ; |
---|
| 641 | h[61]->fill(std::log10(tkin),de) ; |
---|
| 642 | de = hIon[0]->ComputeDEDX(part[1],Si,tkin) ; |
---|
| 643 | h[62]->fill(std::log10(tkin),de) ; |
---|
| 644 | de = hIon[0]->ComputeDEDX(part[1],Cu,tkin) ; |
---|
| 645 | h[63]->fill(std::log10(tkin),de) ; |
---|
| 646 | de = hIon[0]->ComputeDEDX(part[1],Fe,tkin) ; |
---|
| 647 | h[64]->fill(std::log10(tkin),de) ; |
---|
| 648 | de = hIon[0]->ComputeDEDX(part[1],Pb,tkin) ; |
---|
| 649 | h[65]->fill(std::log10(tkin),de) ; |
---|
| 650 | de = hIon[0]->ComputeDEDX(part[1],ethane,tkin) ; |
---|
| 651 | // G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl; |
---|
| 652 | h[66]->fill(std::log10(tkin),de) ; |
---|
| 653 | de = hIon[0]->ComputeDEDX(part[1],water,tkin) ; |
---|
| 654 | // G4cout << "water : E = " << tkin << "; dedx = " << de << G4endl; |
---|
| 655 | h[67]->fill(std::log10(tkin),de) ; |
---|
| 656 | de = hIon[0]->ComputeDEDX(part[1],LAr,tkin) ; |
---|
| 657 | h[68]->fill(std::log10(tkin),de) ; |
---|
| 658 | de = hIon[0]->ComputeDEDX(part[1],csi,tkin) ; |
---|
| 659 | h[69]->fill(std::log10(tkin),de) ; |
---|
| 660 | } |
---|
| 661 | |
---|
| 662 | G4cout << "AntiProton's dEdx histograms are filled!" << G4endl; |
---|
| 663 | |
---|
| 664 | G4double mProt = part[0]->GetPDGMass()*1.007276; |
---|
| 665 | G4double fact = cm/(2700.0*MeV) ; // to MeV/mg/cm^2 |
---|
| 666 | |
---|
| 667 | for (j = 0 ; j < num-1 ; j++) { |
---|
| 668 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
| 669 | de = theIonEffChargeModel->TheValue(part[3],Cu,tkin) ; |
---|
| 670 | // G4cout << "E = " << tkin << "; dedx = " << de << G4endl; |
---|
| 671 | h[48]->fill(std::log10(tkin),de) ; |
---|
| 672 | de = theIonEffChargeModel->TheValue(part[4],Cu,tkin) ; |
---|
| 673 | h[49]->fill(std::log10(tkin),de) ; |
---|
| 674 | G4double tRed = tkin * (part[3]->GetPDGMass())/mProt ; |
---|
| 675 | de = hIon[3]->ComputeDEDX(part[3],Al,tRed) ; |
---|
| 676 | de += theNuclearStoppingModel->TheValue(part[3],Al,tRed) ; |
---|
| 677 | h[50]->fill(std::log10(tkin),de*fact) ; |
---|
| 678 | tRed = tkin * (part[4]->GetPDGMass())/mProt ; |
---|
| 679 | de = hIon[4]->ComputeDEDX(part[4],Al,tRed) ; |
---|
| 680 | //de += theNuclearStoppingModel->TheValue(part[4],Al,tRed) ; |
---|
| 681 | h[51]->fill(std::log10(tkin),de*fact) ; |
---|
| 682 | tRed = tkin * (part[5]->GetPDGMass())/mProt ; |
---|
| 683 | de = hIon[5]->ComputeDEDX(part[5],Al,tRed) ; |
---|
| 684 | //de += theNuclearStoppingModel->TheValue(part[5],Al,tRed) ; |
---|
| 685 | h[52]->fill(std::log10(tkin),de*fact) ; |
---|
| 686 | } |
---|
| 687 | |
---|
| 688 | G4cout << "Ions's dEdx histograms are filled!" << G4endl; |
---|
| 689 | |
---|
| 690 | theProcessManager[7] = new G4ProcessManager(part[0]); |
---|
| 691 | part[0]->SetProcessManager(theProcessManager[7]); |
---|
| 692 | hIon[7] = new G4hLowEnergyIonisation(); |
---|
| 693 | hIon[7]->SetElectronicStoppingPowerModel(part[0],"Ziegler1985p"); |
---|
| 694 | hIon[7]->SetEnlossFluc(false) ; |
---|
| 695 | //hIon[7]->SetNuclearStoppingOff(); |
---|
| 696 | // hIon[7]->SetBarkasOff(); |
---|
| 697 | theProcessManager[7]->AddProcess(hIon[7]); |
---|
| 698 | hIon[7]->BuildPhysicsTable(*part[0]); |
---|
| 699 | |
---|
| 700 | G4cout << "Ziegler's dEdx histograms will be filled" << G4endl; |
---|
| 701 | |
---|
| 702 | for (j = 0 ; j < num-1 ; j++) { |
---|
| 703 | tkin = std::pow(10.0,(std::log10(minE) + (G4double(j)+0.5)*s)); |
---|
| 704 | de = hIon[7]->ComputeDEDX(part[0],Graphite,tkin) ; |
---|
| 705 | h[38]->fill(std::log10(tkin),de) ; |
---|
| 706 | de = hIon[7]->ComputeDEDX(part[0],Al,tkin) ; |
---|
| 707 | h[39]->fill(std::log10(tkin),de) ; |
---|
| 708 | de = hIon[7]->ComputeDEDX(part[0],Si,tkin) ; |
---|
| 709 | h[40]->fill(std::log10(tkin),de) ; |
---|
| 710 | de = hIon[7]->ComputeDEDX(part[0],Cu,tkin) ; |
---|
| 711 | h[41]->fill(std::log10(tkin),de) ; |
---|
| 712 | de = hIon[7]->ComputeDEDX(part[0],Fe,tkin) ; |
---|
| 713 | h[42]->fill(std::log10(tkin),de) ; |
---|
| 714 | de = hIon[7]->ComputeDEDX(part[0],Pb,tkin) ; |
---|
| 715 | h[43]->fill(std::log10(tkin),de) ; |
---|
| 716 | de = hIon[7]->ComputeDEDX(part[0],ethane,tkin) ; |
---|
| 717 | // G4cout << "ethane: E = " << tkin << "; dedx = " << de << G4endl; |
---|
| 718 | h[44]->fill(std::log10(tkin),de) ; |
---|
| 719 | de = hIon[7]->ComputeDEDX(part[0],water,tkin) ; |
---|
| 720 | // G4cout << "water: E = " << tkin << "; dedx = " << de << G4endl; |
---|
| 721 | h[45]->fill(std::log10(tkin),de) ; |
---|
| 722 | de = hIon[7]->ComputeDEDX(part[0],LAr,tkin) ; |
---|
| 723 | h[46]->fill(std::log10(tkin),de) ; |
---|
| 724 | de = hIon[7]->ComputeDEDX(part[0],csi,tkin) ; |
---|
| 725 | h[47]->fill(std::log10(tkin),de) ; |
---|
| 726 | } |
---|
| 727 | |
---|
| 728 | //---------------------- Fill Ntuple ------------------------ |
---|
| 729 | /* |
---|
| 730 | G4cout << "Fill Ntuple!" << G4endl; |
---|
| 731 | |
---|
| 732 | // ---- primary ntuple ------ |
---|
| 733 | HepTuple* ntuple = hbookManager->ntuple("dEdx ntuple"); |
---|
| 734 | assert (ntuple != 0); |
---|
| 735 | |
---|
| 736 | for( i = 0; i < num; i++) { |
---|
| 737 | tkin = std::pow(10,(std::log10(minE) + i*s)); |
---|
| 738 | |
---|
| 739 | for ( G4int j = 0 ; j < numOfMaterials; j++ ) { |
---|
| 740 | |
---|
| 741 | material = (*theMaterialTable)[j] ; |
---|
| 742 | // get elements in the actual material, |
---|
| 743 | const G4ElementVector* theElementVector = material->GetElementVector() ; |
---|
| 744 | const G4double* theAtomicNumDensityVector = |
---|
| 745 | material->GetAtomicNumDensityVector() ; |
---|
| 746 | const G4int NumberOfElements = material->GetNumberOfElements() ; |
---|
| 747 | |
---|
| 748 | // loop for the elements in the material |
---|
| 749 | // to find out average values Z, vF, lF |
---|
| 750 | G4double z = 0.0, norm = 0.0 ; |
---|
| 751 | |
---|
| 752 | if( 1 == NumberOfElements ) { |
---|
| 753 | z = material->GetZ() ; |
---|
| 754 | |
---|
| 755 | } else { |
---|
| 756 | for (G4int iel=0; iel<NumberOfElements; iel++) { |
---|
| 757 | const G4Element* element = (*theElementVector)(iel) ; |
---|
| 758 | G4double z2 = element->GetZ() ; |
---|
| 759 | const G4double weight = theAtomicNumDensityVector[iel] ; |
---|
| 760 | norm += weight ; |
---|
| 761 | z += z2 * weight ; |
---|
| 762 | } |
---|
| 763 | z /= norm ; |
---|
| 764 | } |
---|
| 765 | |
---|
| 766 | for (G4int k=0 ; k<7; k++) { |
---|
| 767 | G4double dedx = hIon[k]->ComputeDEDX(part[k],material,tkin) ; |
---|
| 768 | G4double q = theIonEffChargeModel->TheValue(part[k],material,tkin); |
---|
| 769 | ntuple->column("part",k); |
---|
| 770 | ntuple->column("mat",j); |
---|
| 771 | ntuple->column("ie",i); |
---|
| 772 | ntuple->column("zeff",z); |
---|
| 773 | ntuple->column("tkin",tkin/MeV); |
---|
| 774 | ntuple->column("mass",(part[k]->GetPDGMass())/MeV); |
---|
| 775 | ntuple->column("char",(part[k]->GetPDGCharge())/eplus); |
---|
| 776 | ntuple->column("cha2",q); |
---|
| 777 | ntuple->column("dedx",dedx*mm/MeV); |
---|
| 778 | ntuple->dumpData(); |
---|
| 779 | } |
---|
| 780 | } |
---|
| 781 | } |
---|
| 782 | */ |
---|
| 783 | //----------- End of work ------------------------------------- |
---|
| 784 | |
---|
| 785 | std::cout << "Committing..." << std::endl; |
---|
| 786 | tree->commit(); |
---|
| 787 | std::cout << "Closing the tree..." << std::endl; |
---|
| 788 | tree->close(); |
---|
| 789 | |
---|
| 790 | G4cout << "Ntuple and Hbook are saved" << G4endl; |
---|
| 791 | |
---|
| 792 | // delete materials and elements |
---|
| 793 | delete Be; |
---|
| 794 | delete Graphite; |
---|
| 795 | delete Al; |
---|
| 796 | delete Si; |
---|
| 797 | delete LAr; |
---|
| 798 | delete Fe; |
---|
| 799 | delete Cu; |
---|
| 800 | delete W; |
---|
| 801 | delete Pb; |
---|
| 802 | delete U; |
---|
| 803 | delete H; |
---|
| 804 | delete C; |
---|
| 805 | delete Cs; |
---|
| 806 | delete I; |
---|
| 807 | delete Ti; |
---|
| 808 | delete O; |
---|
| 809 | delete water; |
---|
| 810 | delete ethane; |
---|
| 811 | delete csi; |
---|
| 812 | G4cout << "Materials are deleted" << G4endl; |
---|
| 813 | delete theIonEffChargeModel; |
---|
| 814 | delete Z77p; |
---|
| 815 | delete Z77He; |
---|
| 816 | delete I49p; |
---|
| 817 | delete I49He; |
---|
| 818 | |
---|
| 819 | cout<<"END OF THE MAIN PROGRAM"<<G4endl; |
---|
| 820 | } |
---|