| [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 | }
|
|---|
