Ignore:
Timestamp:
Jan 8, 2010, 3:02:48 PM (14 years ago)
Author:
garnier
Message:

update to geant4.9.3

Location:
trunk/examples/advanced/air_shower
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/examples/advanced/air_shower/GNUmakefile

    r807 r1230  
    2323endif
    2424
    25 ifdef G4ANALYSIS_USE
    26   CPPFLAGS += `aida-config --include`
    27   LDFLAGS += `aida-config --lib`
    28 endif
    29 
    3025
    3126.PHONY: all

  • trunk/examples/advanced/air_shower/UltraMacro.mac

    r807 r1230  
    22/run/verbose 2
    33/event/verbose    0
    4 /tracking/verbose 0
     4/tracking/verbose 1
    55/run/initialize
    6 /gps/particle opticalphoton
     6/gps/particle mu-
    77/gps/direction 0.0 0.0 -1.0
    88/gps/pos/centre  0.0 0.0 100. cm
    … …  
    1010/gps/pos/shape  Circle
    1111/gps/pos/radius 25.0 cm
    12 /gps/energy 3.0 eV
     12/gps/energy 300.0 MeV
    1313
    1414/mysetrun/SetRunID 9999
    15 /run/beamOn 100
    16 exit
     15/run/beamOn 1

  • trunk/examples/advanced/air_shower/Visualisation.mac

    r1187 r1230  
    1313#
    1414#/vis/open OGLIX
    15 /vis/open RayTracerX
    16 #/vis/open ATree
    17 #/vis/open OIX
     15#/vis/open RayTracer
    1816#/vis/open OGLSX
    19 #/vis/open DAWNFILE
     17/vis/open DAWNFILE
    2018#/vis/open HepRepFile
    2119#/vis/open VRML2FILE
    … …  
    3533/vis/scene/endOfEventAction accumulate
    3634#/vis/drawVolume
     35/gps/particle opticalphoton
     36/gps/energy 2.0 eV     
     37/mysetrun/SetRunID 9999
     38/run/beamOn 100
    3739
    38 
    39 /mysetrun/SetRunID 9999
    40 /run/beamOn 50
    41 

  • trunk/examples/advanced/air_shower/include/UltraAnalysisManager.hh

    r807 r1230  
    5858  class IAnalysisFactory;
    5959  class ITreeFactory;
    60 };
     60}
    6161#endif
    6262

  • trunk/examples/advanced/air_shower/include/UltraDetectorConstruction.hh

    r807 r1230  
    6060class UltraDetectorConstruction : public G4VUserDetectorConstruction
    6161{
     62 
    6263  public:
    6364    UltraDetectorConstruction();
    … …  
    6667  public:
    6768    G4VPhysicalVolume* Construct();
     69
     70    inline G4double GetLambdaMin() const {return lambda_min;}
     71    inline G4double GetLambdaMax() const {return lambda_max;}
    6872
    6973  private:
    … …  
    8286    UltraPMTSD*   PMTSD  ;          //pointer to the photomultiplier sensitive detector
    8387    G4SDManager*  SDmanager ;       // Sensitive Detector Manager
    84 
     88   
     89    G4double lambda_min ;
     90    G4double lambda_max ;
    8591};
    8692

  • trunk/examples/advanced/air_shower/src/UltraDetectorConstruction.cc

    r807 r1230  
    7676 // Sensitive Detector Manager
    7777 SDmanager = G4SDManager::GetSDMpointer();
     78
     79// Define wavelength limits for materials definition
     80 lambda_min = 200*nm ;
     81 lambda_max = 700*nm ;
     82
    7883}
    7984
    … …  
    225230/////////////////////////////////////////////
    226231
    227   const G4int NUMENTRIES = 32;
     232  const G4int NUMENTRIES = 2;
    228233
    229234  // Energy bins
    230   G4double X_RINDEX[NUMENTRIES] =
    231             { 2.034E-9*GeV, 2.068E-9*GeV, 2.103E-9*GeV, 2.139E-9*GeV,
    232               2.177E-9*GeV, 2.216E-9*GeV, 2.256E-9*GeV, 2.298E-9*GeV,
    233               2.341E-9*GeV, 2.386E-9*GeV, 2.433E-9*GeV, 2.481E-9*GeV,
    234               2.532E-9*GeV, 2.585E-9*GeV, 2.640E-9*GeV, 2.697E-9*GeV,
    235               2.757E-9*GeV, 2.820E-9*GeV, 2.885E-9*GeV, 2.954E-9*GeV,
    236               3.026E-9*GeV, 3.102E-9*GeV, 3.181E-9*GeV, 3.265E-9*GeV,
    237               3.353E-9*GeV, 3.446E-9*GeV, 3.545E-9*GeV, 3.649E-9*GeV,
    238               3.760E-9*GeV, 3.877E-9*GeV, 4.002E-9*GeV, 4.136E-9*GeV } ;
     235  G4double X_RINDEX[NUMENTRIES] = {h_Planck*c_light/lambda_max, h_Planck*c_light/lambda_min} ;
    239236
    240237
    241238  // Air
    242   G4double RINDEX_AIR[NUMENTRIES] =
    243             { 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
    244               1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
    245               1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
    246               1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00,
    247               1.00, 1.00, 1.00, 1.00 } ;
     239  G4double RINDEX_AIR[NUMENTRIES] = {1.00, 1.00} ;
    248240
    249241// Air refractive index at 20 oC and 1 atm (from PDG)
    … …  
    265257
    266258  const G4int N_RINDEX_QUARTZ = 2 ;
    267   G4double X_RINDEX_QUARTZ[N_RINDEX_QUARTZ] = {0.0*eV, 10.0*eV};
     259  G4double X_RINDEX_QUARTZ[N_RINDEX_QUARTZ] = {h_Planck*c_light/lambda_max, h_Planck*c_light/lambda_min} ;
    268260  G4double RINDEX_QUARTZ[N_RINDEX_QUARTZ] = {1.54, 1.54};
    269261
    … …  
    280272// Refractive index
    281273
    282   const G4int    N_RINDEX_ACRYLIC = 3 ;
    283   G4double X_RINDEX_ACRYLIC[N_RINDEX_ACRYLIC] = {320.0, 400.0, 500.0};  // Wavelength in nanometers
    284   G4double RINDEX_ACRYLIC[N_RINDEX_ACRYLIC] = {1.526, 1.507, 1.497};
    285 
    286     // Convert from nm to GeV
    287 
    288      for(G4int i=0;i<N_RINDEX_ACRYLIC; i++){
    289       X_RINDEX_ACRYLIC[i] = ((1239.84/X_RINDEX_ACRYLIC[i])*1E-9)*GeV;
    290       }
     274  const G4int NENTRIES = 11 ;
     275  G4double LAMBDA_ACRYLIC[NENTRIES] ;
     276
     277
     278  G4double RINDEX_ACRYLIC[NENTRIES] ;
     279  G4double ENERGY_ACRYLIC[NENTRIES] ;
     280
     281// Parameterization for refractive index of High Grade PMMA
     282
     283  G4double bParam[4] = {1760.7010,-1.3687,2.4388e-3,-1.5178e-6} ;
     284 
     285  for(G4int i=0;i<NENTRIES; i++){
     286 
     287    LAMBDA_ACRYLIC[i] = lambda_min + i*(lambda_max-lambda_min)/float(NENTRIES-1) ;
     288    RINDEX_ACRYLIC[i] = 0.0 ;
     289
     290    for (G4int jj=0 ; jj<4 ; jj++)
     291    {
     292      RINDEX_ACRYLIC[i] +=  (bParam[jj]/1000.0)*std::pow(LAMBDA_ACRYLIC[i]/nm,jj) ;
     293    }
     294
     295    ENERGY_ACRYLIC[i] =   h_Planck*c_light/LAMBDA_ACRYLIC[i] ;  // Convert from wavelength to energy ;
     296//  G4cout << ENERGY_ACRYLIC[i]/eV << " " << LAMBDA_ACRYLIC[i]/nm << " " << RINDEX_ACRYLIC[i] << G4endl ;
     297
     298  }
    291299
    292300  G4MaterialPropertiesTable *MPT_Acrylic = new G4MaterialPropertiesTable();
    293   MPT_Acrylic->AddProperty("RINDEX", X_RINDEX_ACRYLIC, RINDEX_ACRYLIC, N_RINDEX_ACRYLIC);
     301  MPT_Acrylic->AddProperty("RINDEX", ENERGY_ACRYLIC, RINDEX_ACRYLIC, NENTRIES);
     302
     303
     304// Absorption
     305  const G4int NENT = 25 ;
     306  G4double LAMBDAABS[NENT] =
     307  {
     308    100.0,
     309    246.528671, 260.605103, 263.853516, 266.019104, 268.726105,   
     310    271.433136, 273.598724, 276.305725, 279.554138, 300.127380,   
     311    320.159241, 340.191101, 360.764343, 381.337585, 399.745239,   
     312    421.401276, 440.891724, 460.382172, 480.414001, 500.987274,   
     313    520.477722, 540.509583, 559.458618,
     314    700.0   
     315  } ;
     316
     317  G4double ABS[NENT] =   // Transmission (in %) of  3mm thick PMMA
     318  {
     319    0.0000000,
     320    0.0000000,  5.295952,  9.657321, 19.937695, 29.283491,
     321    39.252335, 48.598133, 58.255451, 65.109039, 79.439247,
     322    85.669785, 89.719627, 91.277260, 91.588783, 91.900307,
     323    91.588783, 91.277260, 91.277260, 91.588783, 91.588783,
     324    91.900307, 91.900307, 91.588783,
     325    91.5
     326  } ;
     327
     328
     329  MPT_Acrylic->AddProperty("ABSLENGTH", new G4MaterialPropertyVector()) ;
     330  for(G4int i=0;i<NENT; i++){
     331    G4double energy    = h_Planck*c_light/(LAMBDAABS[i]*nm) ;
     332    G4double abslength ;
     333
     334    if (ABS[i] <= 0.0) {
     335      abslength = 1.0/kInfinity ;
     336    }
     337    else {
     338      abslength = -3.0*mm/(std::log(ABS[i]/100.0)) ;
     339    }
     340
     341    MPT_Acrylic->AddEntry("ABSLENGTH", energy, abslength);
     342
     343  }
     344
    294345  Acrylic->SetMaterialPropertiesTable(MPT_Acrylic);
    295 
     346 
    296347
    297348//////////////////////////////////////////////////////////////////
    … …  
    344395
    345396const G4int NUM = 2;
    346 G4double XX[NUM] = { 0.1E-9*GeV, 10.0E-9*GeV };
     397G4double XX[NUM] = {h_Planck*c_light/lambda_max, h_Planck*c_light/lambda_min} ;
    347398G4double ICEREFLECTIVITY[NUM]      = { 0.95, 0.95 };
    348399
    … …  
    403454
    404455 const G4int NUM = 2;
    405 G4double XX[NUM] = { 0.1E-9*GeV, 10.0E-9*GeV };
     456G4double XX[NUM] = {h_Planck*c_light/lambda_max, h_Planck*c_light/lambda_min} ;
    406457G4double ICEREFLECTIVITY[NUM]      = { 0.95, 0.95 };
    407458
    … …  
    566617
    567618const G4int NUM = 2;
    568 G4double XX[NUM] = { 2.030E-9*GeV, 4.144E-9*GeV };
     619G4double XX[NUM] = {h_Planck*c_light/lambda_max, h_Planck*c_light/lambda_min} ;
    569620G4double BLACKPAINTREFLECTIVITY[NUM]      = { 0.05, 0.05 };
    570621//G4double WHITEPAINTREFLECTIVITY[NUM]      = { 0.99, 0.99 };

  • trunk/examples/advanced/air_shower/src/UltraPhysicsList.cc

    r807 r1230  
    164164#include "G4PhotoElectricEffect.hh"
    165165
    166 #include "G4MultipleScattering.hh"
     166#include "G4eMultipleScattering.hh"
     167#include "G4MuMultipleScattering.hh"
     168#include "G4hMultipleScattering.hh"
    167169
    168170#include "G4eIonisation.hh"
    … …  
    196198    //electron
    197199      // Construct processes for electron
    198       pmanager->AddProcess(new G4MultipleScattering(),-1,1,1);
     200      pmanager->AddProcess(new G4eMultipleScattering(),-1,1,1);
    199201      pmanager->AddProcess(new G4eIonisation(),-1,2,2);
    200202      pmanager->AddProcess(new G4eBremsstrahlung(),-1,-1,3);
    … …  
    203205    //positron
    204206      // Construct processes for positron
    205       pmanager->AddProcess(new G4MultipleScattering(),-1,1,1);
     207      pmanager->AddProcess(new G4eMultipleScattering(),-1,1,1);
    206208      pmanager->AddProcess(new G4eIonisation(),-1,2,2);
    207209      pmanager->AddProcess(new G4eBremsstrahlung(),-1,-1,3);
    … …  
    212214    //muon
    213215     // Construct processes for muon
    214      pmanager->AddProcess(new G4MultipleScattering(),-1,1,1);
     216     pmanager->AddProcess(new G4MuMultipleScattering(),-1,1,1);
    215217     pmanager->AddProcess(new G4MuIonisation(),-1,2,2);
    216218     pmanager->AddProcess(new G4MuBremsstrahlung(),-1,-1,3);
    … …  
    221223          (particle->GetParticleName() != "chargedgeantino")) {
    222224     // all others charged particles except geantino
    223        pmanager->AddProcess(new G4MultipleScattering(),-1,1,1);
     225       pmanager->AddProcess(new G4hMultipleScattering(),-1,1,1);
    224226       pmanager->AddProcess(new G4hIonisation(),-1,2,2);
    225227     }
    … …  
    240242  // this Cerenkov Process
    241243  G4Cerenkov*   theCerenkovProcess = new G4Cerenkov("Cerenkov");
    242   // this is Scintillation process
    243   G4Scintillation* theScintillationProcess = new G4Scintillation("Scintillation");
    244244  // this absorption process inside optical media
    245245  G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
    … …  
    251251  // Chose level 0 (no verbose)
    252252  theCerenkovProcess           -> SetVerboseLevel(0);
    253   theScintillationProcess      -> SetVerboseLevel(0);
    254253  theAbsorptionProcess         -> SetVerboseLevel(0);
    255254  theRayleighScatteringProcess -> SetVerboseLevel(0);
    … …  
    262261  theCerenkovProcess->SetTrackSecondariesFirst(true);
    263262
    264   theScintillationProcess->SetTrackSecondariesFirst(true);
    265   theScintillationProcess->SetScintillationYieldFactor(1.);
    266   theScintillationProcess->SetScintillationExcitationRatio(0.0);
    267 
    268263  // Boundary model (UNIFIED OR GLISUR (OLD GEANT3)) For now only GEANT3
    269264  G4OpticalSurfaceModel themodel = unified;
    … …  
    277272
    278273    if (theCerenkovProcess->IsApplicable(*particle)) {
    279       //      pmanager->AddProcess(theCerenkovProcess);
    280       //      pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
     274      pmanager->AddProcess(theCerenkovProcess);
     275      pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
    281276    }
    282277
    283     if (theScintillationProcess->IsApplicable(*particle)) {
    284       pmanager->AddProcess(theScintillationProcess);
    285       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxAtRest);
    286       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxPostStep);
    287     }
    288278
    289279    if (particleName == "opticalphoton") {

  • trunk/examples/advanced/air_shower/src/UltraPrimaryGeneratorAction.cc

    r807 r1230  
    4343//
    4444#include "UltraPrimaryGeneratorAction.hh"
     45#include "UltraDetectorConstruction.hh"
    4546
     47#include "G4RunManager.hh"
    4648#include "G4Event.hh"
    4749#include "G4GeneralParticleSource.hh"
    … …  
    129131    G4cout << particleGun->GetCurrentSource()->GetPosDist()->GetPosDisType()    << G4endl ;
    130132
     133
     134// Check if optical photon wavelength is within limits set for material optical properties tables.
     135
     136   
     137
     138
    131139}
    132140  particleGun->GeneratePrimaryVertex(anEvent);
     141
     142    if (particleGun->GetParticleDefinition()->GetParticleName() == "opticalphoton"){
     143     
     144        const UltraDetectorConstruction * detector = 
     145        dynamic_cast<const UltraDetectorConstruction *>((G4RunManager::GetRunManager())->GetUserDetectorConstruction()) ;
     146
     147        G4double lambda_min = detector->GetLambdaMin() ;
     148        G4double lambda_max = detector->GetLambdaMax() ;
     149
     150        G4double energy = particleGun->GetParticleEnergy() ;
     151
     152        if (h_Planck*c_light/energy > lambda_max || h_Planck*c_light/energy < lambda_min){
     153               G4cerr << "Error ! Optical photon energy (" << energy/eV << " eV) out of limits set by material optical properties tables. \n"
     154              << "Please check that photon wavelength is within the following interval: ["
     155              << lambda_min/nm << ","
     156              << lambda_max/nm << "] nm"
     157              << ", i.e., ["
     158              << h_Planck*c_light/lambda_max/eV << ","
     159              << h_Planck*c_light/lambda_min/eV << "] eV"
     160              << G4endl ;
     161       
     162             G4Exception("") ;
     163        }
     164 }
     165
    133166}
    134167
Note: See TracChangeset for help on using the changeset viewer.