eic-opticks/u4/Local_G4Cerenkov_modified.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // $Id: G4Cerenkov.cc 108508 2018-02-15 15:54:35Z gcosmo $
0028 //
0029 ////////////////////////////////////////////////////////////////////////
0030 // Cerenkov Radiation Class Implementation
0031 ////////////////////////////////////////////////////////////////////////
0032 //
0033 // File:        G4Cerenkov.cc
0034 // Description: Discrete Process -- Generation of Cerenkov Photons
0035 // Version:     2.1
0036 // Created:     1996-02-21
0037 // Author:      Juliet Armstrong
0038 // Updated:     2007-09-30 by Peter Gumplinger
0039 //              > change inheritance to G4VDiscreteProcess
0040 //              GetContinuousStepLimit -> GetMeanFreePath (StronglyForced)
0041 //              AlongStepDoIt -> PostStepDoIt
0042 //              2005-08-17 by Peter Gumplinger
0043 //              > change variable name MeanNumPhotons -> MeanNumberOfPhotons
0044 //              2005-07-28 by Peter Gumplinger
0045 //              > add G4ProcessType to constructor
0046 //              2001-09-17, migration of Materials to pure STL (mma)
0047 //              2000-11-12 by Peter Gumplinger
0048 //              > add check on CerenkovAngleIntegrals->IsFilledVectorExist()
0049 //              in method GetAverageNumberOfPhotons
0050 //              > and a test for MeanNumberOfPhotons <= 0.0 in DoIt
0051 //              2000-09-18 by Peter Gumplinger
0052 //              > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
0053 //                        aSecondaryTrack->SetTouchable(0);
0054 //              1999-10-29 by Peter Gumplinger
0055 //              > change: == into <= in GetContinuousStepLimit
0056 //              1997-08-08 by Peter Gumplinger
0057 //              > add protection against /0
0058 //              > G4MaterialPropertiesTable; new physics/tracking scheme
0059 //
0060 // mail:        gum@triumf.ca
0061 //
0062 ////////////////////////////////////////////////////////////////////////
0063
0064 #include "G4ios.hh"
0065 #include "G4PhysicalConstants.hh"
0066 #include "G4SystemOfUnits.hh"
0067 #include "G4Poisson.hh"
0068 #include "G4EmProcessSubType.hh"
0069
0070 #include "G4LossTableManager.hh"
0071 #include "G4MaterialCutsCouple.hh"
0072 #include "G4ParticleDefinition.hh"
0073
0074 #include "G4Version.hh"
0075
0076 #include "Local_G4Cerenkov_modified.hh"
0077
0078 #ifdef INSTRUMENTED
0079 #include "OpticksDebug.hh"
0080 #include "OpticksRandom.hh"
0081 #endif
0082
0083 #ifdef STANDALONE
0084 #include "SLOG.hh"
0085 #include "U4.hh"
0086 #endif
0087
0088
0089 /////////////////////////
0090 // Class Implementation
0091 /////////////////////////
0092
0093 //G4bool G4Cerenkov::fTrackSecondariesFirst = false;
0094 //G4double G4Cerenkov::fMaxBetaChange = 0.;
0095 //G4int G4Cerenkov::fMaxPhotons = 0;
0096
0097   //////////////
0098   // Operators
0099   //////////////
0100
0101 // G4Cerenkov::operator=(const G4Cerenkov &right)
0102 // {
0103 // }
0104
0105   /////////////////
0106   // Constructors
0107   /////////////////
0108
0109 Local_G4Cerenkov_modified::Local_G4Cerenkov_modified(const G4String& processName, G4ProcessType type)
0110            : G4VProcess(processName, type),
0111              fTrackSecondariesFirst(false),
0112              fMaxBetaChange(0.0),
0113              fMaxPhotons(0),
0114              fStackingFlag(true),
0115 #ifdef INSTRUMENTED
0116              override_fNumPhotons(0),
0117 #endif
0118              fNumPhotons(0)
0119 {
0120   SetProcessSubType(fCerenkov);
0121
0122   thePhysicsTable = nullptr;
0123
0124   if (verboseLevel>0) {
0125      G4cout << GetProcessName() << " is created " << G4endl;
0126   }
0127 }
0128
0129 // G4Cerenkov::G4Cerenkov(const G4Cerenkov &right)
0130 // {
0131 // }
0132
0133   ////////////////
0134   // Destructors
0135   ////////////////
0136
0137 Local_G4Cerenkov_modified::~Local_G4Cerenkov_modified()
0138 {
0139   if (thePhysicsTable != nullptr) {
0140      thePhysicsTable->clearAndDestroy();
0141      delete thePhysicsTable;
0142   }
0143 }
0144
0145   ////////////
0146   // Methods
0147   ////////////
0148
0149 G4bool Local_G4Cerenkov_modified::IsApplicable(const G4ParticleDefinition& aParticleType)
0150 {
0151   G4bool result = false;
0152
0153   if (aParticleType.GetPDGCharge() != 0.0 &&
0154       aParticleType.GetPDGMass() != 0.0 &&
0155       aParticleType.GetParticleName() != "chargedgeantino" &&
0156       !aParticleType.IsShortLived() ) { result = true; }
0157
0158   return result;
0159 }
0160
0161 void Local_G4Cerenkov_modified::SetTrackSecondariesFirst(const G4bool state)
0162 {
0163   fTrackSecondariesFirst = state;
0164 }
0165
0166 void Local_G4Cerenkov_modified::SetMaxBetaChangePerStep(const G4double value)
0167 {
0168   fMaxBetaChange = value*CLHEP::perCent;
0169 }
0170
0171 void Local_G4Cerenkov_modified::SetMaxNumPhotonsPerStep(const G4int NumPhotons)
0172 {
0173   fMaxPhotons = NumPhotons;
0174 }
0175
0176 void Local_G4Cerenkov_modified::BuildPhysicsTable(const G4ParticleDefinition&)
0177 {
0178   if (!thePhysicsTable) BuildThePhysicsTable();
0179 }
0180
0181 // PostStepDoIt
0182 // -------------
0183 //
0184 G4VParticleChange* Local_G4Cerenkov_modified::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
0185
0186 // This routine is called for each tracking Step of a charged particle
0187 // in a radiator. A Poisson-distributed number of photons is generated
0188 // according to the Cerenkov formula, distributed evenly along the track
0189 // segment and uniformly azimuth w.r.t. the particle direction. The
0190 // parameters are then transformed into the Master Reference System, and
0191 // they are added to the particle change.
0192
0193 {
0194   ////////////////////////////////////////////////////
0195   // Should we ensure that the material is dispersive?
0196   ////////////////////////////////////////////////////
0197
0198
0199   aParticleChange.Initialize(aTrack);
0200
0201   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0202   const G4Material* aMaterial = aTrack.GetMaterial();
0203
0204   G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
0205   G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
0206
0207   G4ThreeVector x0 = pPreStepPoint->GetPosition();
0208   G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
0209   G4double t0 = pPreStepPoint->GetGlobalTime();
0210
0211   G4MaterialPropertiesTable* aMaterialPropertiesTable =
0212                                aMaterial->GetMaterialPropertiesTable();
0213   if (!aMaterialPropertiesTable) return pParticleChange;
0214
0215   G4MaterialPropertyVector* Rindex =
0216                 aMaterialPropertiesTable->GetProperty(kRINDEX);
0217   if (!Rindex) return pParticleChange;
0218
0219   // particle charge
0220   G4double charge = aParticle->GetDefinition()->GetPDGCharge();
0221
0222   // particle beta
0223   G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta())*0.5;
0224
0225   fNumPhotons = 0;
0226
0227   G4double MeanNumberOfPhotons =
0228                      GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
0229
0230   if (MeanNumberOfPhotons <= 0.0) {
0231
0232      // return unchanged particle and no secondaries
0233
0234      aParticleChange.SetNumberOfSecondaries(0);
0235
0236      return pParticleChange;
0237
0238   }
0239
0240   G4double step_length = aStep.GetStepLength();
0241
0242   MeanNumberOfPhotons = MeanNumberOfPhotons * step_length;
0243
0244   fNumPhotons = (G4int) G4Poisson(MeanNumberOfPhotons);
0245
0246 #ifdef INSTRUMENTED
0247    if( override_fNumPhotons > 0 )
0248    {
0249        fNumPhotons = override_fNumPhotons ;
0250    }
0251 #endif
0252
0253
0254   // calculate the fNumPhotons1 and fNumPhotons2 {
0255
0256   // }
0257
0258
0259   // if ( fNumPhotons <= 0 || !fStackingFlag ) {
0260   if ( fNumPhotons <= 0 ) {
0261
0262      // return unchanged particle and no secondaries
0263
0264      aParticleChange.SetNumberOfSecondaries(0);
0265
0266      fNumPhotons1 = 0;
0267      fNumPhotons2 = 0;
0268
0269      return pParticleChange;
0270
0271   }
0272
0273   ////////////////////////////////////////////////////////////////
0274 #if G4VERSION_NUMBER < 1100
0275   G4double Pmin = Rindex->GetMinLowEdgeEnergy();
0276   G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
0277 #else
0278   G4double Pmin = Rindex->GetMinEnergy();
0279   G4double Pmax = Rindex->GetMaxEnergy();
0280 #endif
0281   G4double dp = Pmax - Pmin;
0282
0283
0284 #ifdef FLOAT_TEST
0285   float nMax = Rindex->GetMaxValue();
0286 #else
0287   G4double nMax = Rindex->GetMaxValue();
0288 #endif
0289   if (Rindex) {
0290       // nMax = Rindex->GetMaxValue();
0291       size_t ri_sz = Rindex->GetVectorLength();
0292
0293       for (size_t i = 0; i < ri_sz; ++i) {
0294           if ((*Rindex)[i]>nMax) nMax = (*Rindex)[i];
0295       }
0296   }
0297
0298   G4double BetaInverse = 1./beta;
0299
0300 #ifdef FLOAT_TEST
0301   float maxCos = BetaInverse / nMax;
0302   float maxSin2 = (1.f - maxCos) * (1.f + maxCos);
0303 #else
0304   G4double maxCos = BetaInverse / nMax;
0305   G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
0306 #endif
0307
0308   G4double beta1 = pPreStepPoint ->GetBeta();
0309   G4double beta2 = pPostStepPoint->GetBeta();
0310
0311   G4double MeanNumberOfPhotons1 =
0312                      GetAverageNumberOfPhotons(charge,beta1,aMaterial,Rindex);
0313   G4double MeanNumberOfPhotons2 =
0314                      GetAverageNumberOfPhotons(charge,beta2,aMaterial,Rindex);
0315
0316   fNumPhotons1 = MeanNumberOfPhotons1;
0317   fNumPhotons2 = MeanNumberOfPhotons2;
0318
0319
0320 #ifdef INSTRUMENTED
0321   par->append( BetaInverse, "BetaInverse" );
0322   par->append( beta       , "beta" );
0323   par->append( Pmin       , "Pmin" );
0324   par->append( Pmax       , "Pmax" );
0325
0326   par->append( nMax       , "nMax" );
0327   par->append( maxCos     , "maxCos" );
0328   par->append( maxSin2    , "maxSin2" );
0329   par->append( fNumPhotons, "fNumPhotons" );
0330 #endif
0331
0332
0333
0334   if (MeanNumberOfPhotons1 <= 0.0 or MeanNumberOfPhotons2<=0.0) {
0335
0336      // return unchanged particle and no secondaries
0337
0338      aParticleChange.SetNumberOfSecondaries(0);
0339
0340      // Force no secondaries
0341      fNumPhotons = 0;
0342      fNumPhotons1 = 0;
0343      fNumPhotons2 = 0;
0344
0345      return pParticleChange;
0346
0347   }
0348
0349   // if stacking is false, then stop the generation of Cerenkov photons
0350   if (!fStackingFlag) {
0351      aParticleChange.SetNumberOfSecondaries(0);
0352
0353      return pParticleChange;
0354   }
0355
0356   ////////////////////////////////////////////////////////////////
0357
0358 #ifdef STANDALONE
0359   fNumPhotons = std::min( fNumPhotons, 5 );   // artifical reduction for debugging convenience
0360 #endif
0361
0362   aParticleChange.SetNumberOfSecondaries(fNumPhotons);
0363
0364   if (fTrackSecondariesFirst) {
0365      if (aTrack.GetTrackStatus() == fAlive )
0366          aParticleChange.ProposeTrackStatus(fSuspend);
0367   }
0368
0369
0370 #ifdef STANDALONE
0371   U4::CollectGenstep_G4Cerenkov_modified(
0372       &aTrack,
0373       &aStep,
0374       fNumPhotons,
0375       BetaInverse,
0376       Pmin,
0377       Pmax,
0378       maxCos,
0379       maxSin2,
0380       MeanNumberOfPhotons1,
0381       MeanNumberOfPhotons2
0382   );
0383 #endif
0384
0385 #ifdef STANDALONE
0386     U4::GenPhotonAncestor(&aTrack);
0387 #endif
0388
0389   for (G4int i = 0; i < fNumPhotons; i++) {
0390       // Determine photon energy
0391
0392 #ifdef STANDALONE
0393       U4::GenPhotonBegin(i);
0394 #endif
0395
0396 #ifdef FLOAT_TEST
0397       float rand(0.f);
0398       float rand0(0.f);
0399       float rand1(0.f) ;
0400       float sampledEnergy(0.f);
0401       float sampledRI(0.f);
0402       float cosTheta(0.f) ;
0403       float sin2Theta(0.f) ;
0404 #else
0405       G4double rand(0.);
0406       G4double rand0(0.);
0407       G4double rand1(0.) ;
0408       G4double sampledEnergy(0.);
0409       G4double sampledRI(0.);
0410       G4double cosTheta(0.);
0411       G4double sin2Theta(0.);
0412 #endif
0413
0414 #ifdef INSTRUMENTED
0415       unsigned head_count = 0 ;
0416       unsigned tail_count = 0 ;
0417       unsigned continue_count = 0 ;
0418       unsigned condition_count = 0 ;
0419       int seqidx = -1 ;
0420       if(rnd)
0421       {
0422           rnd->setSequenceIndex(i);
0423           seqidx = rnd->getSequenceIndex();
0424
0425           if(i < 10) std::cout
0426               << " i " << std::setw(6) << i
0427               << " seqidx " << std::setw(7) << seqidx
0428               << " Pmin/eV " << std::fixed << std::setw(10) << std::setprecision(5) << Pmin/eV
0429               << " Pmax/eV " << std::fixed << std::setw(10) << std::setprecision(5) << Pmax/eV
0430               << " dp/eV " << std::fixed << std::setw(10) << std::setprecision(5) << dp/eV
0431               << " maxSin2 "  << std::fixed << std::setw(10) << std::setprecision(5) << maxSin2
0432               << std::endl
0433               ;
0434
0435       }
0436 #endif
0437
0438       // sample an energy
0439
0440       do {
0441 #ifdef INSTRUMENTED
0442          head_count += 1 ;
0443 #endif
0444          rand0 = G4UniformRand();
0445          sampledEnergy = Pmin + rand0 * dp;
0446          sampledRI = Rindex->Value(sampledEnergy);
0447
0448
0449 #ifdef SKIP_CONTINUE
0450 #else
0451          // check if n(E) > 1/beta
0452          if (sampledRI < BetaInverse) {
0453 #ifdef INSTRUMENTED
0454              continue_count += 1 ;
0455 #endif
0456              continue;
0457          }
0458
0459 #endif
0460
0461 #ifdef INSTRUMENTED
0462          tail_count += 1 ;
0463 #endif
0464
0465
0466          cosTheta = BetaInverse / sampledRI;
0467
0468          // G4cout << "TAO --> L" << __LINE__ << ": "
0469          //        << " BetaInverse : " << BetaInverse
0470          //        << " sampledRI : " << sampledRI
0471          //        << " cosTheta : " << cosTheta
0472          //        << G4endl;
0473
0474 #ifdef FLOAT_TEST
0475          sin2Theta = (1.f - cosTheta)*(1.f + cosTheta);
0476 #else
0477          sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
0478 #endif
0479
0480          rand1 = G4UniformRand();
0481 #ifdef ONE_RAND
0482          rand1 = 1.0 ;
0483 #endif
0484
0485         // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
0486 #ifdef INSTRUMENTED
0487
0488          if( i < 10 ) std::cout
0489              << " tc " << std::setw(6) << tail_count
0490              << " u0 " << std::fixed << std::setw(10) << std::setprecision(5) << rand0
0491              << " eV " << std::fixed << std::setw(10) << std::setprecision(5) << sampledEnergy/eV
0492              << " ri " << std::fixed << std::setw(10) << std::setprecision(5) << sampledRI
0493              << " ct " << std::fixed << std::setw(10) << std::setprecision(5) << cosTheta
0494              << " s2 " << std::fixed << std::setw(10) << std::setprecision(5) << sin2Theta
0495              << " rand1*maxSin2 " << std::fixed << std::setw(10) << std::setprecision(5) << rand1*maxSin2
0496              << " rand1*maxSin2 - sin2Theta " <<  std::fixed << std::setw(10) << std::setprecision(5) << rand1*maxSin2 - sin2Theta
0497              << " loop " << ( rand1*maxSin2 > sin2Theta ? "Y" : "N" )
0498              << std::endl
0499              ;
0500
0501
0502       } while ( looping_condition(condition_count) && rand1*maxSin2 > sin2Theta  );
0503 #else
0504       } while (rand1*maxSin2 > sin2Theta);
0505 #endif
0506
0507 #ifdef INSTRUMENTED
0508         G4double sampledEnergy_eV = sampledEnergy/eV ;
0509         G4double sampledWavelength_nm = h_Planck*c_light/sampledEnergy/nm ;
0510
0511         gen->append( sampledEnergy_eV ,       "sampledEnergy" );
0512         gen->append( sampledWavelength_nm ,    "sampledWavelength" );
0513         gen->append( sampledRI ,               "sampledRI" );
0514         gen->append( cosTheta ,                "cosTheta" );
0515
0516         gen->append( sin2Theta ,               "sin2Theta" );
0517         gen->append( head_count ,     tail_count,       "head_tail" );
0518         gen->append( continue_count , condition_count,  "continue_condition" );
0519         gen->append( BetaInverse , "BetaInverse" );
0520
0521         if(rnd)
0522         {
0523            rnd->setSequenceIndex(-1);
0524         }
0525 #endif
0526
0527
0528
0529
0530       // Generate random position of photon on cone surface
0531       // defined by Theta
0532
0533       rand = G4UniformRand();
0534
0535       G4double phi = twopi*rand;
0536       G4double sinPhi = std::sin(phi);
0537       G4double cosPhi = std::cos(phi);
0538
0539       // calculate x,y, and z components of photon energy
0540       // (in coord system with primary particle direction
0541       //  aligned with the z axis)
0542
0543       G4double sinTheta = std::sqrt(sin2Theta);
0544       G4double px = sinTheta*cosPhi;
0545       G4double py = sinTheta*sinPhi;
0546       G4double pz = cosTheta;
0547
0548       // Create photon momentum direction vector
0549       // The momentum direction is still with respect
0550       // to the coordinate system where the primary
0551       // particle direction is aligned with the z axis
0552
0553       G4ParticleMomentum photonMomentum(px, py, pz);
0554
0555       // Rotate momentum direction back to global reference
0556       // system
0557
0558       photonMomentum.rotateUz(p0);
0559
0560       // Determine polarization of new photon
0561
0562       G4double sx = cosTheta*cosPhi;
0563       G4double sy = cosTheta*sinPhi;
0564       G4double sz = -sinTheta;
0565
0566       G4ThreeVector photonPolarization(sx, sy, sz);
0567
0568       // Rotate back to original coord system
0569
0570       photonPolarization.rotateUz(p0);
0571
0572       // Generate a new photon:
0573
0574       G4DynamicParticle* aCerenkovPhoton =
0575         new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),photonMomentum);
0576
0577       aCerenkovPhoton->SetPolarization(photonPolarization.x(),
0578                                        photonPolarization.y(),
0579                                        photonPolarization.z());
0580
0581       aCerenkovPhoton->SetKineticEnergy(sampledEnergy);
0582
0583       // Generate new G4Track object:
0584
0585       G4double NumberOfPhotons, N;
0586
0587       do {
0588          rand = G4UniformRand();
0589          NumberOfPhotons = MeanNumberOfPhotons1 - rand *
0590                                 (MeanNumberOfPhotons1-MeanNumberOfPhotons2);
0591          N = G4UniformRand() *
0592                         std::max(MeanNumberOfPhotons1,MeanNumberOfPhotons2);
0593         // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
0594       } while (N > NumberOfPhotons);
0595
0596       G4double delta = rand * aStep.GetStepLength();
0597
0598       G4double deltaTime = delta / (pPreStepPoint->GetVelocity()+
0599                                       rand*(pPostStepPoint->GetVelocity()-
0600                                             pPreStepPoint->GetVelocity())*0.5);
0601
0602       G4double aSecondaryTime = t0 + deltaTime;
0603
0604       G4ThreeVector aSecondaryPosition = x0 + rand * aStep.GetDeltaPosition();
0605
0606       G4Track* aSecondaryTrack =
0607                new G4Track(aCerenkovPhoton,aSecondaryTime,aSecondaryPosition);
0608
0609       aSecondaryTrack->SetTouchableHandle(
0610                                aStep.GetPreStepPoint()->GetTouchableHandle());
0611
0612       aSecondaryTrack->SetParentID(aTrack.GetTrackID());
0613
0614       aParticleChange.AddSecondary(aSecondaryTrack);
0615
0616 #ifdef STANDALONE
0617       U4::GenPhotonEnd(i, aSecondaryTrack);
0618 #endif
0619
0620
0621
0622   }
0623
0624   if (verboseLevel>0) {
0625      G4cout <<"\n Exiting from Local_G4Cerenkov_modified::DoIt -- NumberOfSecondaries = "
0626       << aParticleChange.GetNumberOfSecondaries() << G4endl;
0627   }
0628
0629     return pParticleChange;
0630 }
0631
0632
0633 #ifdef INSTRUMENTED
0634 bool Local_G4Cerenkov_modified::looping_condition(unsigned& count)
0635 {
0636     count += 1 ;
0637     return true ;
0638 }
0639
0640
0641 #endif
0642
0643
0644 // BuildThePhysicsTable for the Cerenkov process
0645 // ---------------------------------------------
0646 //
0647
0648 void Local_G4Cerenkov_modified::BuildThePhysicsTable()
0649 {
0650   if (thePhysicsTable) return;
0651
0652   const G4MaterialTable* theMaterialTable=
0653   G4Material::GetMaterialTable();
0654   G4int numOfMaterials = G4Material::GetNumberOfMaterials();
0655
0656   // create new physics table
0657
0658   thePhysicsTable = new G4PhysicsTable(numOfMaterials);
0659
0660   // loop for materials
0661
0662   for (G4int i=0 ; i < numOfMaterials; i++) {
0663
0664       G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = 0;
0665
0666       // Retrieve vector of refraction indices for the material
0667       // from the material's optical properties table
0668
0669       G4Material* aMaterial = (*theMaterialTable)[i];
0670
0671       G4MaterialPropertiesTable* aMaterialPropertiesTable =
0672                                       aMaterial->GetMaterialPropertiesTable();
0673
0674       if (aMaterialPropertiesTable) {
0675          aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
0676          G4MaterialPropertyVector* theRefractionIndexVector =
0677                               aMaterialPropertiesTable->GetProperty(kRINDEX);
0678
0679          if (theRefractionIndexVector) {
0680
0681             // Retrieve the first refraction index in vector
0682             // of (photon energy, refraction index) pairs
0683
0684             G4double currentRI = (*theRefractionIndexVector)[0];
0685
0686             if (currentRI > 1.0) {
0687
0688                // Create first (photon energy, Cerenkov Integral)
0689                // pair
0690
0691                G4double currentPM = theRefractionIndexVector->Energy(0);
0692                G4double currentCAI = 0.0;
0693
0694                aPhysicsOrderedFreeVector->InsertValues(currentPM , currentCAI);
0695
0696                // Set previous values to current ones prior to loop
0697
0698                G4double prevPM  = currentPM;
0699                G4double prevCAI = currentCAI;
0700                G4double prevRI  = currentRI;
0701
0702                // loop over all (photon energy, refraction index)
0703                // pairs stored for this material
0704
0705                for (size_t ii = 1;
0706                            ii < theRefractionIndexVector->GetVectorLength();
0707                            ++ii) {
0708                    currentRI = (*theRefractionIndexVector)[ii];
0709                    currentPM = theRefractionIndexVector->Energy(ii);
0710
0711                    currentCAI = 0.5*(1.0/(prevRI*prevRI) +
0712                                      1.0/(currentRI*currentRI));
0713
0714                    currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
0715
0716                    aPhysicsOrderedFreeVector->
0717                                          InsertValues(currentPM, currentCAI);
0718
0719                    prevPM  = currentPM;
0720                    prevCAI = currentCAI;
0721                    prevRI  = currentRI;
0722                }
0723
0724             }
0725          }
0726       }
0727
0728       // The Cerenkov integral for a given material
0729       // will be inserted in thePhysicsTable
0730       // according to the position of the material in
0731       // the material table.
0732
0733       thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector);
0734
0735   }
0736 }
0737
0738 // GetMeanFreePath
0739 // ---------------
0740 //
0741
0742 G4double Local_G4Cerenkov_modified::GetMeanFreePath(const G4Track&,
0743                                            G4double,
0744                                            G4ForceCondition*)
0745 {
0746   return 1.;
0747 }
0748
0749 G4double Local_G4Cerenkov_modified::PostStepGetPhysicalInteractionLength(
0750                                            const G4Track& aTrack,
0751                                            G4double,
0752                                            G4ForceCondition* condition)
0753 {
0754   *condition = NotForced;
0755   G4double StepLimit = DBL_MAX;
0756
0757   const G4Material* aMaterial = aTrack.GetMaterial();
0758   G4int materialIndex = aMaterial->GetIndex();
0759
0760   // If Physics Vector is not defined no Cerenkov photons
0761   //    this check avoid string comparison below
0762
0763   if(!(*thePhysicsTable)[materialIndex]) { return StepLimit; }
0764
0765   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0766   const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
0767
0768   G4double kineticEnergy = aParticle->GetKineticEnergy();
0769   const G4ParticleDefinition* particleType = aParticle->GetDefinition();
0770   G4double mass = particleType->GetPDGMass();
0771
0772   // particle beta
0773   G4double beta = aParticle->GetTotalMomentum() /
0774                   aParticle->GetTotalEnergy();
0775   // particle gamma
0776   G4double gamma = aParticle->GetTotalEnergy()/mass;
0777
0778
0779   G4MaterialPropertiesTable* aMaterialPropertiesTable =
0780                                       aMaterial->GetMaterialPropertiesTable();
0781
0782   G4MaterialPropertyVector* Rindex = NULL;
0783
0784   if (aMaterialPropertiesTable)
0785                      Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
0786
0787   G4double nMax;
0788   if (Rindex) {
0789       // nMax = Rindex->GetMaxValue();
0790       size_t ri_sz = Rindex->GetVectorLength();
0791       nMax = (*Rindex)[0];
0792       for (size_t i = 1; i < ri_sz; ++i) {
0793           if ((*Rindex)[i]>nMax) nMax = (*Rindex)[i];
0794       }
0795   } else {
0796      return StepLimit;
0797   }
0798
0799   G4double BetaMin = 1./nMax;
0800
0801   if ( BetaMin >= 1. ) return StepLimit;
0802
0803   G4double GammaMin = 1./std::sqrt(1.-BetaMin*BetaMin);
0804
0805   if (gamma < GammaMin ) return StepLimit;
0806
0807   G4double kinEmin = mass*(GammaMin-1.);
0808
0809   G4double RangeMin = G4LossTableManager::Instance()->GetRange(particleType,
0810                                                                kinEmin,
0811                                                                couple);
0812   G4double Range    = G4LossTableManager::Instance()->GetRange(particleType,
0813                                                                kineticEnergy,
0814                                                                couple);
0815
0816   G4double Step = Range - RangeMin;
0817 //  if (Step < 1.*um ) return StepLimit;
0818
0819
0820   if (Step > 0. && Step < StepLimit) StepLimit = Step;
0821
0822   // If user has defined an average maximum number of photons to
0823   // be generated in a Step, then calculate the Step length for
0824   // that number of photons.
0825
0826   if (fMaxPhotons > 0) {
0827      // particle charge
0828      const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
0829
0830      G4double MeanNumberOfPhotons =
0831                       GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
0832
0833      Step = 0.;
0834      if (MeanNumberOfPhotons > 0.0) Step = fMaxPhotons / MeanNumberOfPhotons;
0835
0836      if (Step > 0. && Step < StepLimit) StepLimit = Step;
0837   }
0838
0839   // If user has defined an maximum allowed change in beta per step
0840   if (fMaxBetaChange > 0.) {
0841
0842      G4double dedx = G4LossTableManager::Instance()->GetDEDX(particleType,
0843                                                              kineticEnergy,
0844                                                              couple);
0845
0846      G4double deltaGamma = gamma - 1./std::sqrt(1.-beta*beta*
0847                                                 (1.-fMaxBetaChange)*
0848                                                 (1.-fMaxBetaChange));
0849
0850      Step = mass * deltaGamma / dedx;
0851
0852      if (Step > 0. && Step < StepLimit) StepLimit = Step;
0853
0854   }
0855
0856   *condition = StronglyForced;
0857   return StepLimit;
0858 }
0859
0860 // GetAverageNumberOfPhotons
0861 // -------------------------
0862 // This routine computes the number of Cerenkov photons produced per
0863 // GEANT-unit (millimeter) in the current medium.
0864 //             ^^^^^^^^^^
0865
0866 G4double
0867   Local_G4Cerenkov_modified::GetAverageNumberOfPhotons(const G4double charge,
0868                                         const G4double beta,
0869                       const G4Material* aMaterial,
0870                       G4MaterialPropertyVector* Rindex) //const
0871 {
0872
0873   const G4double Rfact = 369.81/(eV * cm);
0874
0875 #ifdef X_INSTRUMENTED
0876   std::cout
0877        << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons"
0878        << " Rfact " << std::fixed << std::setw(10) << std::setprecision(5) << Rfact
0879        << " eV " << std::fixed << std::setw(10) << std::setprecision(7) << eV
0880        << " cm " << std::fixed << std::setw(10) << std::setprecision(5) << cm
0881        << " charge " << std::fixed << std::setw(10) << std::setprecision(5) << charge
0882        << std::endl
0883        ;
0884 #endif
0885
0886   if(beta <= 0.0)return 0.0;
0887
0888   G4double BetaInverse = 1./beta;
0889
0890   // Vectors used in computation of Cerenkov Angle Integral:
0891   //  - Refraction Indices for the current material
0892   //  - new G4PhysicsOrderedFreeVector allocated to hold CAI's
0893
0894   G4int materialIndex = aMaterial->GetIndex();
0895
0896   // Retrieve the Cerenkov Angle Integrals for this material
0897
0898   G4PhysicsOrderedFreeVector* CerenkovAngleIntegrals =
0899              (G4PhysicsOrderedFreeVector*)((*thePhysicsTable)(materialIndex));
0900 #if G4VERSION_NUMBER < 1100
0901   if(!(CerenkovAngleIntegrals->IsFilledVectorExist()))return 0.0;
0902 #else
0903   G4int length = CerenkovAngleIntegrals->GetVectorLength();
0904   if(0 == length)     return 0.0;
0905 #endif
0906   /*
0907   // Min and Max photon energies
0908   G4double Pmin = Rindex->GetMinLowEdgeEnergy();
0909   G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
0910
0911   // Min and Max Refraction Indices
0912   G4double nMin = Rindex->GetMinValue();
0913   G4double nMax = Rindex->GetMaxValue();
0914
0915   // Max Cerenkov Angle Integral
0916   G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();
0917
0918   G4double dp, ge;
0919
0920   // If n(Pmax) < 1/Beta -- no photons generated
0921
0922   if (nMax < BetaInverse) {
0923      dp = 0.0;
0924      ge = 0.0;
0925   }
0926
0927   // otherwise if n(Pmin) >= 1/Beta -- photons generated
0928
0929   else if (nMin > BetaInverse) {
0930      dp = Pmax - Pmin;
0931      ge = CAImax;
0932   }
0933
0934   // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
0935   // we need to find a P such that the value of n(P) == 1/Beta.
0936   // Interpolation is performed by the GetEnergy() and
0937   // Value() methods of the G4MaterialPropertiesTable and
0938   // the GetValue() method of G4PhysicsVector.
0939
0940   else {
0941      Pmin = Rindex->GetEnergy(BetaInverse);
0942      dp = Pmax - Pmin;
0943
0944      // need boolean for current implementation of G4PhysicsVector
0945      // ==> being phased out
0946      G4bool isOutRange;
0947      G4double CAImin = CerenkovAngleIntegrals->GetValue(Pmin, isOutRange);
0948      ge = CAImax - CAImin;
0949
0950      if (verboseLevel>0) {
0951         G4cout << "CAImin = " << CAImin << G4endl;
0952         G4cout << "ge = " << ge << G4endl;
0953       }
0954       //////// old version ///////////
0955
0956   }
0957
0958         */
0959      ///////// new version ///////////
0960
0961     G4int    cross_num;
0962     // G4double cross_up[10];   // max crossing point number : 10
0963     // G4double cross_down[10];   // max crossing point number : 10
0964     std::vector<double> the_energies_threshold; // there are several pairs (ranges) at the threshold of 1/beta
0965     // [ E_l_0, E_r_0, E_l_1, E_r_1, ... ]
0966     // the energies between [E_l_i, E_r_i] is above the threshold.
0967
0968      cross_num = 0;
0969      size_t vec_length = Rindex->GetVectorLength();
0970
0971      G4double maxRI=(*Rindex)[0]; G4double minRI=(*Rindex)[0];
0972      for (size_t ii = 1;
0973              ii < vec_length;
0974              ++ii) {
0975          G4double currentRI = (*Rindex)[ii];
0976          if (currentRI > maxRI) maxRI = currentRI;
0977          if (currentRI < minRI) minRI = currentRI;
0978      }
0979
0980
0981
0982      if (BetaInverse <= minRI) { // All range is OK!
0983
0984          // cross_up[0] = Rindex->Energy(0);
0985          // cross_down[0] = Rindex->Energy(vec_length-1);
0986          cross_num = 1;
0987
0988          the_energies_threshold.push_back(Rindex->Energy(0));
0989          the_energies_threshold.push_back(Rindex->Energy(vec_length-1));
0990
0991          //G4cout << "Range [ " << cross_up[0] << ", " << cross_down[0] << "]" << G4endl;
0992
0993      } else if (BetaInverse >= maxRI) { // Out of Range
0994          cross_num = 0;
0995
0996      } else {   // between min and max
0997
0998          // below is impl by Tao Lin
0999          double currentRI = (*Rindex)[0];
1000          double currentPM = Rindex->Energy(0);
1001
1002          // first point
1003          if (currentRI >= BetaInverse) {
1004              the_energies_threshold.push_back(currentPM);
1005          }
1006
1007          // middle points
1008          if (vec_length>2) {
1009              for (size_t ii = 1; ii < vec_length-1; ++ii) {
1010                  double prevRI = (*Rindex)[ii-1];
1011                  double prevPM = Rindex->Energy(ii-1);
1012                  double currentRI = (*Rindex)[ii];
1013                  double currentPM = Rindex->Energy(ii);
1014
1015                  // two case here
1016                  if ( (prevRI >= BetaInverse and currentRI < BetaInverse)
1017                       or (prevRI < BetaInverse and currentRI >= BetaInverse) ) {
1018                      double energy_threshold = (BetaInverse-prevRI)/(currentRI-prevRI)*(currentPM-prevPM) + prevPM;
1019                      the_energies_threshold.push_back(energy_threshold);
1020                  }
1021
1022              }
1023          }
1024
1025          // last point
1026          currentRI = (*Rindex)[vec_length-1];
1027          currentPM = Rindex->Energy(vec_length-1);
1028          if (currentRI >= BetaInverse) {
1029              the_energies_threshold.push_back(currentPM);
1030          }
1031
1032          if ((the_energies_threshold.size()%2) != 0) {
1033              G4cerr << "ERROR: missing endpoint for the_energies_threshold? "
1034                     << "The size of the_energies_threshold is "
1035                     << the_energies_threshold.size()
1036                     << G4endl;
1037          }
1038
1039          cross_num = the_energies_threshold.size() / 2;
1040          // for (int i = 0; i < cross_num; ++i) {
1041          //     cross_up[i] = the_energies_threshold[i*2];
1042          //     cross_down[i] = the_energies_threshold[i*2+1];
1043          // }
1044
1045          // END
1046
1047          // // below is impl by Miao Yu
1048
1049          // G4double up_leftx[10], up_lefty[10];
1050          // G4double up_rightx[10], up_righty[10];
1051          // G4double down_leftx[10], down_lefty[10];
1052          // G4double down_rightx[10], down_righty[10];
1053          // G4double extremex[10], extremey[10];
1054          // G4int num0 = 0;
1055          // G4int num1 = 0;
1056          // G4int num2 = 0;
1057
1058          // double currentRI = (*Rindex)[0];
1059          // double currentPM = Rindex->Energy(0);
1060          // if (currentRI == BetaInverse) {
1061          //     extremex[num2] =  currentPM;
1062          //     extremey[num2] = currentRI;
1063          //     num2++;
1064          // }
1065
1066          // for (size_t ii = 1;
1067          //         ii < vec_length;
1068          //         ++ii) {
1069          //     double prevRI = (*Rindex)[ii-1];
1070          //     double prevPM = Rindex->Energy(ii-1);
1071          //     double currentRI = (*Rindex)[ii];
1072          //     double currentPM = Rindex->Energy(ii);
1073
1074          //     if(prevRI < BetaInverse and currentRI > BetaInverse) { // up
1075          //         up_leftx[num0] = prevPM;
1076          //         up_rightx[num0] = currentPM;
1077          //         up_lefty[num0] = prevRI;
1078          //         up_righty[num0] = currentRI;
1079          //         num0++;
1080          //     } else if(prevRI > BetaInverse and currentRI < BetaInverse) {
1081          //         down_leftx[num1] = prevPM;
1082          //         down_rightx[num1] = currentPM;
1083          //         down_lefty[num1] = prevRI;
1084          //         down_righty[num1] = currentRI;
1085          //         num1++;
1086          //     } else if (currentRI == BetaInverse) {
1087          //         extremex[num2] = currentPM;
1088          //         extremey[num2] = currentRI;
1089          //         num2++;
1090          //     }
1091          // }
1092
1093          // if (num0 - num1 == 1) // up > down
1094          // {
1095          //     down_leftx[num1] = Rindex->Energy(vec_length-1);
1096          //     down_rightx[num1] = down_leftx[num1];
1097          //     down_lefty[num1] = (*Rindex)[vec_length-1];
1098          //     down_righty[num1] = down_lefty[num1];
1099          //     num1++;
1100          // } else if(num1 - num0 == 1) {
1101          //     up_leftx[num0] = Rindex->Energy(0);
1102          //     up_rightx[num0] = up_leftx[num0];
1103          //     up_lefty[num0] = (*Rindex)[0];
1104          //     up_righty[num0] = up_lefty[num0];
1105          //     num0++;
1106          // }
1107
1108          // if(num0 != num1) G4cout << "Error: Vector Length Mismatching ! " << G4endl;
1109
1110          // // linear-interpolation for crossing points:
1111          // for (int i=0; i<num0; i++) {
1112
1113          //     if (up_leftx[i] == up_rightx[i]) cross_up[i] = up_leftx[i];
1114          //     else cross_up[i] = (BetaInverse-up_lefty[i])/(up_righty[i] - up_lefty[i])*(up_rightx[i]-up_leftx[i]) + up_leftx[i];
1115          //     if (down_leftx[i] == down_rightx[i]) cross_down[i] = down_leftx[i];
1116          //     else cross_down[i] = (BetaInverse-down_lefty[i])/(down_righty[i] - down_lefty[i])*(down_rightx[i]-down_leftx[i]) + down_leftx[i];
1117          //     cross_num++;
1118
1119          //     //G4cout << "Range [ " << cross_up[i] << ", " << cross_down[i] << "]" << G4endl;
1120          // }
1121      }
1122      G4double dp1 = 0; G4double ge1 = 0;
1123      for (int i=0; i<cross_num; i++) {
1124         dp1 += the_energies_threshold[2*i+1] - the_energies_threshold[2*i];
1125         G4bool isOutRange;
1126         ge1 += CerenkovAngleIntegrals->GetValue(the_energies_threshold[2*i+1], isOutRange)
1127                - CerenkovAngleIntegrals->GetValue(the_energies_threshold[2*i], isOutRange);
1128      }
1129
1130   // Calculate number of photons
1131   //G4double NumPhotons = Rfact * charge/eplus * charge/eplus *
1132   //                               (dp - ge * BetaInverse*BetaInverse);
1133   G4double NumPhotons = Rfact * charge/eplus * charge/eplus *
1134          (dp1 - ge1 * BetaInverse*BetaInverse);
1135
1136
1137 #ifdef X_INSTRUMENTED
1138   std::cout
1139        << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons"
1140        << " BetaInverse " << std::fixed << std::setw(10) << std::setprecision(5) << BetaInverse
1141        << " maxRI " << std::fixed << std::setw(10) << std::setprecision(5) << maxRI
1142        << " minRI " << std::fixed << std::setw(10) << std::setprecision(5) << minRI
1143        << " cross_num " << cross_num
1144        << " dp1 " << std::fixed << std::setw(10) << std::setprecision(5) << dp1
1145        << " dp1/eV " << std::fixed << std::setw(10) << std::setprecision(5) << dp1/eV
1146        << " ge1 " << std::fixed << std::setw(10) << std::setprecision(5) << ge1
1147        << " NumPhotons " << std::fixed << std::setw(10) << std::setprecision(5) << NumPhotons
1148        << std::endl
1149        ;
1150
1151   for(int i=0 ; i < cross_num ; i++)
1152   {
1153
1154       G4bool isOutRange;
1155       G4double cai0 = CerenkovAngleIntegrals->GetValue(the_energies_threshold[2*i+0], isOutRange);
1156       G4double cai1 = CerenkovAngleIntegrals->GetValue(the_energies_threshold[2*i+1], isOutRange);
1157
1158       std::cout
1159            << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons"
1160            << " the_energies_threshold[2*i+0]/eV " << std::fixed << std::setw(10) << std::setprecision(5) << the_energies_threshold[2*i+0]/eV
1161            << " the_energies_threshold[2*i+1]/eV " << std::fixed << std::setw(10) << std::setprecision(5) << the_energies_threshold[2*i+1]/eV
1162            << " cai0 " << std::fixed << std::setw(20) << std::setprecision(10) << cai0
1163            << " cai1 " << std::fixed << std::setw(20) << std::setprecision(10) << cai1
1164            << std::endl
1165            ;
1166   }
1167 #endif
1168
1169
1170   return NumPhotons;
1171 }
1172
1173 /**
1174 Local_G4Cerenkov_modified::s2Integrate
1175 ----------------------------------
1176
1177 See:
1178
1179 * ~/opticks/ana/ckn.py
1180 * ~/np/NP.hh NP::trapz
1181
1182 Comparing the speed of _s2 with _asis shows that this takes
1183 an almost constant time no matter what the BetaInverse
1184 unlike _asis which varies greatly.  The only situation where
1185 _s2 is slower is when the number of photons is zero.
1186
1187 Hence could short circuit that situation to make this always faster.
1188 But its better to do that in the caller scope as max Rindex is
1189 already available there.
1190
1191 s2 numerical integral using trapezoidal approximation,
1192
1193 * x,y: s2 zero crossings corresponding to RINDEX(E) == BetaInverse
1194 * B,C,D : trapezoids
1195 * A,E   : edge triangles
1196
1197                                  *
1198                                 /|\
1199                                / | \
1200                               /  |  \
1201                              /   |   \
1202                             /    |    \
1203                            /     |     \
1204                           /      |      \
1205                  *-------*       |       *
1206                 /|       |       |       |\
1207                / |       |       |       | \
1208               / A|  B    |   C   |   D   |E \
1209      -----0--x---1-------2-------3-------4---y----5----------
1210             /                                 \
1211
1212 * determine x,y energy crossings from BetaInverse-RINDEX(E) crossings, as that is linear
1213   so the linear interpolation to give crossings will be precise.
1214   Using s2 zero crossing not as precise, as non-linear.
1215
1216
1217 Formerly did::
1218
1219    G4double en_cross =  cross ? (s2_1*en_0 - s2_0*en_1)/(s2_1 - s2_0) : -1. ;
1220
1221 But it is more precise to find en_cross from BetaInverse-ri crossings directly
1222 rather than the derived and non-linear s2
1223
1224 **/
1225
1226 G4double Local_G4Cerenkov_modified::GetAverageNumberOfPhotons_s2(const G4double charge, const G4double beta, const G4Material*, G4MaterialPropertyVector* Rindex) const
1227 {
1228     if(beta <= 0.0)return 0.0;
1229     G4double BetaInverse = 1./beta;
1230     G4double half(0.5) ;
1231     G4double s2integral(0.) ;
1232
1233     for(unsigned i=0 ; i < Rindex->GetVectorLength()-1 ; i++)
1234     {
1235         G4double en_0 = Rindex->Energy(i) ;
1236         G4double en_1 = Rindex->Energy(i+1) ;
1237
1238         G4double ri_0 = (*Rindex)[i] ;
1239         G4double ri_1 = (*Rindex)[i+1] ;
1240
1241         G4double ct_0 = BetaInverse/ri_0 ;
1242         G4double ct_1 = BetaInverse/ri_1 ;
1243
1244         G4double s2_0 = (1.-ct_0)*(1.+ct_0) ;
1245         G4double s2_1 = (1.-ct_1)*(1.+ct_1) ;
1246
1247         G4bool cross = s2_0*s2_1 < 0. ;
1248         G4double en_cross =  cross ? en_0 + (BetaInverse - ri_0)*(en_1 - en_0)/(ri_1 - ri_0) : -1. ;
1249
1250         if( s2_0 <= 0. and s2_1 <= 0. )  // no CK
1251         {
1252             // noop
1253         }
1254         else if( s2_0 < 0. and s2_1 > 0. )  // s2 becomes +ve within the bin
1255         {
1256             s2integral +=  (en_1 - en_cross)*s2_1*half ;
1257         }
1258         else if( s2_0 >= 0. and s2_1 >= 0. )   // s2 +ve across full bin
1259         {
1260             s2integral += (en_1 - en_0)*(s2_0 + s2_1)*half ;
1261         }
1262         else if( s2_0 > 0. and s2_1 < 0. )  // s2 becomes -ve within the bin
1263         {
1264             s2integral +=  (en_cross - en_0)*s2_0*half ;
1265         }
1266         else
1267         {
1268             std::cout
1269                 << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons_s2"
1270                 << " FATAL "
1271                 << " s2_0 " << std::fixed << std::setw(10) << std::setprecision(5) << s2_0
1272                 << " s2_1 " << std::fixed << std::setw(10) << std::setprecision(5) << s2_1
1273                 << " en_0/eV " << std::fixed << std::setw(10) << std::setprecision(5) << en_0/eV
1274                 << " en_1/eV " << std::fixed << std::setw(10) << std::setprecision(5) << en_1/eV
1275                 << " ri_0 " << std::fixed << std::setw(10) << std::setprecision(5) << ri_0
1276                 << " ri_1 " << std::fixed << std::setw(10) << std::setprecision(5) << ri_1
1277                 << std::endl
1278                 ;
1279             assert(0);
1280         }
1281     }
1282     const G4double Rfact = 369.81/(eV * cm);
1283     G4double NumPhotons = Rfact * charge/eplus * charge/eplus * s2integral ;
1284
1285 #ifdef X_INSTRUMENTED
1286   std::cout
1287        << "Local_G4Cerenkov_modified::GetAverageNumberOfPhotons_s2"
1288        << " Rfact " << std::fixed << std::setw(10) << std::setprecision(5) << Rfact
1289        << " eV " << std::fixed << std::setw(10) << std::setprecision(7) << eV
1290        << " cm " << std::fixed << std::setw(10) << std::setprecision(5) << cm
1291        << " mm " << std::fixed << std::setw(10) << std::setprecision(5) << mm
1292        << " charge " << std::fixed << std::setw(10) << std::setprecision(5) << charge
1293        << " s2integral " << std::fixed << std::setw(10) << std::setprecision(5) << s2integral
1294        << " NumPhotons " << std::fixed << std::setw(10) << std::setprecision(5) << NumPhotons
1295        << std::endl
1296        ;
1297 #endif
1298     return NumPhotons ;
1299 }
1300
1301
1302 void Local_G4Cerenkov_modified::DumpPhysicsTable() const
1303 {
1304   G4int PhysicsTableSize = thePhysicsTable->entries();
1305   G4PhysicsOrderedFreeVector *v;
1306
1307   for (G4int i = 0 ; i < PhysicsTableSize ; i++ ) {
1308       v = (G4PhysicsOrderedFreeVector*)(*thePhysicsTable)[i];
1309       v->DumpValues();
1310   }
1311 }
1312
1313