/** @file GLG4Scint.cc
    For GLG4Scint class, providing advanced scintillation process.
    Distantly based on an extensively modified version of G4Scintillation.cc.

    This file is part of the GenericLAND software library.
    $Id$

    @author Glenn Horton-Smith (Tohoku) 28-Jan-1999
*/
////////////////////////////////////////////////////////////////////////
/// REVISION HISTORY:\n
///         30-Jun-2012 : Nuno Barros - corrected compilation errors in DEBUG mode. \n
///         25-Jun-2015 : Ben Land - increased max possible scintillation time. \n
///         07-Nov-2016 : Anthony LaTorre - added rise time to scintillator time
///                       distribution. \n
////////////////////////////////////////////////////////////////////////

// [see detailed class description in GLG4Scint.hh]

#include "GLG4Scint.hh"
#include "G4ios.hh"
#include "G4Timer.hh"
#include "Randomize.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIdirectory.hh"
#include "G4TrackFastVector.hh"  // for G4TrackFastVectorSize
#include <RAT/PhotonThinning.hh>
#include <RAT/fileio.hpp>
#include <RAT/Log.hh>

using CLHEP::twopi;

#include <RAT/UserTrackInformation.hh>

/////////////////////////
// Class Implementation
/////////////////////////

//////////////
// Operators
//////////////

// GLG4Scint::operator=(const GLG4Scint &right)
// {
// }

////////////////
// static data members
////////////////

G4std::vector<GLG4Scint *> GLG4Scint::fMasterVectorOfGLG4Scint;
// top level of scintillation command
G4UIdirectory* GLG4Scint::fGLG4ScintDir = 0;

// universal maximum number of secondary tracks per step for GLG4Scint
//G4int GLG4Scint::maxTracksPerStep = G4TrackFastVectorSize;
G4int GLG4Scint::fMaxTracksPerStep = 180000;

// universal mean number of true photons per secondary track in GLG4Scint
G4double GLG4Scint::fMeanPhotonsPerSecondary = 1.0;

// universal on/off flag
G4bool GLG4Scint::fDoScintillation = true;
G4bool GLG4Scint::fDoReemission = true;

// energy deposition
G4double GLG4Scint::fTotEdep = 0.0;
G4double GLG4Scint::fTotEdepQuenched = 0.0;
G4double GLG4Scint::fTotEdepTime = 0.0;
G4ThreeVector GLG4Scint::fScintCentroidSum( 0.0, 0.0, 0.0 );

//default Quenching Factor
G4double GLG4Scint::fQuenchingFactor=1.0;
//user-given (constant) quenching factor flag
G4bool GLG4Scint::fUserQF=false;
//default primary particle
G4String GLG4Scint::fPrimaryName=G4String();

G4double GLG4Scint::fPrimaryEnergy=0.0;

GLG4DummyProcess GLG4Scint::fScintProcess("Scintillation", fUserDefined);
GLG4DummyProcess GLG4Scint::fReemissionProcess("Reemission", fUserDefined);
G4std::list<GLG4DummyProcess*> GLG4Scint::fReemissionProcessVector;

unsigned int GLG4Scint::fsScintillatedCount = 0;
unsigned int GLG4Scint::fsReemittedCount = 0;

/////////////////
// Constructors
/////////////////

GLG4Scint::GLG4Scint(const G4String& tablename, G4double lowerMassLimit)
{
  fVerboseLevel= 0;
  fLowerMassLimit = lowerMassLimit;

  fMyPhysicsTable = MyPhysicsTable::FindOrBuild( tablename );
  fMyPhysicsTable -> IncUsedBy();
  if(fVerboseLevel) fMyPhysicsTable->Dump();

  // add to ordered list
  if (fMasterVectorOfGLG4Scint.size() == 0
      || lowerMassLimit >= fMasterVectorOfGLG4Scint.back()->fLowerMassLimit){

      fMasterVectorOfGLG4Scint.push_back(this);
  }
  else for (G4std::vector<GLG4Scint *>::iterator
            i = fMasterVectorOfGLG4Scint.begin(); i != fMasterVectorOfGLG4Scint.end(); i++){
    if (lowerMassLimit < (*i)->fLowerMassLimit){
      fMasterVectorOfGLG4Scint.insert(i, this);
      break;
    }
  }

  // create UI commands if necessary
  if (fGLG4ScintDir == NULL) {
     // the scintillation control commands
     new G4UIdirectory("/rat/physics/");
     fGLG4ScintDir = new G4UIdirectory("/rat/physics/scintillation/");
     fGLG4ScintDir->SetGuidance("scintillation process control.");
     G4UIcommand *cmd;
     cmd= new G4UIcommand("/rat/physics/scintillation/on", this);
     cmd->SetGuidance("Turn on scintillation");
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
     cmd= new G4UIcommand("/rat/physics/scintillation/off", this);
     cmd->SetGuidance("Turn off scintillation");
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
     cmd= new G4UIcommand("/rat/physics/scintillation/reemission", this);
     cmd->SetGuidance("Turn on/off reemission of absorbed opticalphotons");
     cmd->SetParameter(new G4UIparameter("status", 's', false));
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
     cmd= new G4UIcommand("/rat/physics/scintillation/maxTracksPerStep",this);
     cmd->SetGuidance("Set maximum number of opticalphoton tracks per step\n"
                      "(If more real photons are needed, "
                      "weight of tracked particles is increased.)\n" );
     cmd->SetParameter(new G4UIparameter("maxTracksPerStep", 'i', false));
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
     cmd= new G4UIcommand("/rat/physics/scintillation/meanPhotonsPerSecondary",this);
     cmd->SetGuidance("Set mean number of \"real\" photons per secondary\n");
     cmd->SetParameter(new G4UIparameter("meanPhotonsPerSecondary", 'd', false));
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
     cmd= new G4UIcommand("/rat/physics/scintillation/verbose",this);
     cmd->SetGuidance("Set verbose level");
     cmd->SetParameter(new G4UIparameter("level", 'i', false));
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
     cmd= new G4UIcommand("/rat/physics/scintillation/dump",this);
     cmd->SetGuidance("Dump tables");
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
     cmd= new G4UIcommand("/rat/physics/scintillation/setQF",this);
     cmd->SetGuidance("Set a constant quenching factor, default is 1");
     cmd->SetParameter(new G4UIparameter("QuenchingFactor", 'd', false));
     cmd->AvailableForStates( G4State_Idle, G4State_GeomClosed, G4State_EventProc );
  }

#ifdef RATVERBOSE
  G4cout << "GLG4Scint[" << tablename << "]"
         << " is created " << G4endl;
#endif
}

// GLG4Scint::GLG4Scint(const GLG4Scint &right)
// {
// }

////////////////
// Destructors
////////////////

GLG4Scint::~GLG4Scint()
{
  fMyPhysicsTable -> DecUsedBy();
  for (G4std::vector<GLG4Scint *>::iterator i=fMasterVectorOfGLG4Scint.begin();
           i != fMasterVectorOfGLG4Scint.end(); i++)
    if (*i == this) {
      fMasterVectorOfGLG4Scint.erase(i);
      break;
    }
}

////////////
// Methods
////////////

//Sets the quenching factor
void GLG4Scint::SetQuenchingFactor(G4double qf=1.0){ fQuenchingFactor = qf; }

// PostStepDoIt
// -------------
//

#ifdef RATDEBUG
G4double GLG4Scint_tottime = 0.0;
G4int GLG4Scint_num_calls= 0;
G4int GLG4Scint_num_phots= 0;
#endif

G4VParticleChange*
GLG4Scint::PostPostStepDoIt(const G4Track& aTrack, const G4Step& aStep)

  // This routine is called for each step of any particle
  // in a scintillator.  For accurate energy deposition, must be called
  // from user-supplied UserSteppingAction, which also must stack
  // any particles created.  A pseudo-Poisson-distributed number of
  // photons is generated according to the scintillation yield formula,
  // distributed evenly along the track segment and uniformly into 4pi.

{
#ifdef RATDEBUG
  G4Timer timer;
  timer.Start();
  GLG4Scint_num_calls ++;
#endif
  G4bool flagReemission = false;
  {
  // prepare to generate an event, organizing to
  // check for things that cause an early exit.
  fAParticleChange.Initialize(aTrack);
  const G4Material* aMaterial = aTrack.GetMaterial();
  const MyPhysicsTable::Entry* physicsEntry = fMyPhysicsTable->GetEntry(aMaterial->GetIndex() );

  if (aTrack.GetDefinition() == G4OpticalPhoton::OpticalPhoton()) {
     flagReemission= fDoReemission
     && aTrack.GetTrackStatus() == fStopAndKill
     && aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary;
     if (!flagReemission)
        goto PostStepDoIt_DONE;
  }

  G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();

  if (TotalEnergyDeposit <= 0.0 && !flagReemission)
     goto PostStepDoIt_DONE;

  // get pointer to the physics entry
  if (!physicsEntry)
     goto PostStepDoIt_DONE;

  //finds E-dependent QF, unless the user provided an E-independent one
  if (!fUserQF && *(fMyPhysicsTable->fName)==fPrimaryName)
    {
      if (physicsEntry->fQuenchingArray)
        {
          for( unsigned int iEntry = 0; iEntry < physicsEntry->fQuenchingArray->GetVectorLength(); iEntry++ )
            {
              //preparations
              G4double CurrentEnergy=physicsEntry->fQuenchingArray->Energy( iEntry );
              G4double PreviousEnergy=CurrentEnergy;
              G4double PrimEn=GetPrimaryEnergy();
              G4double slope;
              //if the primary is below the energy range with a QF, then use QF of lowest energy
              if (PrimEn < CurrentEnergy)
                SetQuenchingFactor(physicsEntry->fQuenchingArray->Value(CurrentEnergy));
              else
                { //find 1st energy above primary, if available
                  while ((PrimEn > CurrentEnergy) && iEntry++ )
                    {
                      PreviousEnergy=CurrentEnergy;
                      CurrentEnergy=physicsEntry->fQuenchingArray->Energy( iEntry );
                    }
                  //if primary energy above range or in quenching array, use QF of last energy
                  if (PrimEn >= CurrentEnergy)
                    SetQuenchingFactor(physicsEntry->fQuenchingArray->Value(CurrentEnergy));
                  else //otherwise interpolates QF
                    {
                      slope=(physicsEntry->fQuenchingArray->Value(CurrentEnergy)-physicsEntry->fQuenchingArray->Value(PreviousEnergy))/(CurrentEnergy-PreviousEnergy);
                      SetQuenchingFactor(slope*(PrimEn-PreviousEnergy)+physicsEntry->fQuenchingArray->Value(PreviousEnergy));
                    }
                }
            }
        }
      else
        SetQuenchingFactor(1.0); // Default to 1 if no table and no user quenching factor
    }

  // Retrieve the Light Yield or Scintillation Integral for this material

  G4double ScintillationYield=physicsEntry->fLightYield;
  G4PhysicsOrderedFreeVector* ScintillationIntegral = physicsEntry->fSpectrumIntegral;
  G4PhysicsOrderedFreeVector* ReemissionIntegral = physicsEntry->fReemissionIntegral;

  if (!ScintillationIntegral)
     goto PostStepDoIt_DONE;

  // If no LY defined Max Scintillation Integral == ScintillationYield
  if (!ScintillationYield)
     ScintillationYield= ScintillationIntegral->GetMaxValue();

  // set positions, directions, etc.
  G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();

  G4ThreeVector x0 = pPreStepPoint->GetPosition();
  G4ThreeVector p0 = pPreStepPoint->GetMomentumDirection();
  G4double      t0 = pPreStepPoint->GetGlobalTime();

  // Finally ready to start generating the event
  // figure out how many photons we want to make

  G4int numSecondaries;
  G4double weight, p_reemission=0, num_comp=0;

  // Marker to record which component absorbed the photon
  int absorbed=0;

  if (flagReemission) {
     G4MaterialPropertiesTable* mpt_scint= aMaterial->GetMaterialPropertiesTable();
     if (mpt_scint->ConstPropertyExists("NUM_COMP"))
       num_comp= mpt_scint->GetConstProperty("NUM_COMP");

     // Reemission from a single component material is treated the same way as before
     if(!num_comp){
       G4MaterialPropertyVector* mpv_scint_reemission= mpt_scint->GetProperty("REEMISSION_PROB");
       if(mpv_scint_reemission == 0)
         goto PostStepDoIt_DONE;
       p_reemission= mpv_scint_reemission->Value(aTrack.GetKineticEnergy());
     }
     // Otherwise we have multi-component material and need a different approach
     else{
       const G4double tot = mpt_scint->GetProperty("ABSLENGTH")->Value( aTrack.GetKineticEnergy() );
       const G4double rand = G4UniformRand();
       G4double prob=0;
       for(int j=0; j<num_comp; j++){
         // calculate the probability that each component absorbs
         char temp[32];
         sprintf(temp, "ABSLENGTH%d", j);

         G4MaterialPropertyVector* mpv_absl= mpt_scint->GetProperty(temp);
         G4double absl=DBL_MAX;
         if( mpv_absl ){
           absl= mpv_absl->Value(aTrack.GetKineticEnergy());
           prob += tot/absl;

           if( rand <= prob ){
             // if the component absorbed the photon, get the reemission probability
             char temp2[32];
             sprintf(temp2, "REEMISSION_PROB%d", j);
             G4MaterialPropertyVector *mpv_reem = mpt_scint->GetProperty(temp2);
             p_reemission = mpv_reem->Value(aTrack.GetKineticEnergy());
             absorbed = j;
             goto ENDLOOP;
           }
         }
       }
     }
  ENDLOOP:
     /*
      This approach was adopted by DEAP due to the oddities of their TPB wavelengthshifter
      but we never expect >1 photon from reemission of LAB ppo.
      numSecondaries= (G4int)(CLHEP::RandPoisson::shoot(p_reemission));
      if (numSecondaries==0)
      goto PostStepDoIt_DONE;
     */
     // use original method for SNO+
     if (G4UniformRand() >= p_reemission) goto PostStepDoIt_DONE;
          numSecondaries = 1;
     weight= aTrack.GetWeight();
  }
  else { // apply Birk's law
     G4double birksConstant = physicsEntry->fBirksConstant;
     G4double QuenchedTotalEnergyDeposit= TotalEnergyDeposit;
     if ( birksConstant != 0.0 ) {
             G4double dE_dx = TotalEnergyDeposit /  aStep.GetStepLength();
             QuenchedTotalEnergyDeposit/= (1.0 + birksConstant * dE_dx);
     }

     // track total edep, quenched edep
     fTotEdep += TotalEnergyDeposit;
     fTotEdepQuenched += QuenchedTotalEnergyDeposit;
     fTotEdepTime = t0;
     fScintCentroidSum += QuenchedTotalEnergyDeposit *( x0 + p0*(0.5*aStep.GetStepLength()) );

     // now we are done if we are not actually making photons here
     if ( !fDoScintillation )
             goto PostStepDoIt_DONE;

     // calculate MeanNumPhotons
     G4double MeanNumPhotons = ScintillationYield*GetQuenchingFactor()
          * QuenchedTotalEnergyDeposit
          * (1.0 + birksConstant * (physicsEntry->fRefdEdx) );

     if (MeanNumPhotons <= 0.0)
             goto PostStepDoIt_DONE;

     // randomize number of TRACKS (not photons)
     // this gets statistics right for number of PE after applying
     // boolean random choice to final absorbed track (change from
     // old method of applying binomial random choice to final absorbed
     // track, which did want poissonian number of photons divided
     // as evenly as possible into tracks)
     // Note for weight=1, there's no difference between tracks and photons.
     G4double MeanNumTracks= MeanNumPhotons/fMeanPhotonsPerSecondary
          / RAT::PhotonThinning::GetFactor();

     G4double resolutionScale= physicsEntry->fResolutionScale;
     if (MeanNumTracks > 12.0)
             numSecondaries= (G4int)(CLHEP::RandGauss::shoot(MeanNumTracks,
                                 resolutionScale * sqrt(MeanNumTracks)));
     else {
             if (resolutionScale > 1.0)
                MeanNumTracks = CLHEP::RandGauss::shoot (MeanNumTracks, sqrt(resolutionScale*resolutionScale-1.0 )*MeanNumTracks);
             numSecondaries =(G4int)( CLHEP::RandPoisson::shoot(MeanNumTracks) );
     }

     weight= fMeanPhotonsPerSecondary;
     if (numSecondaries > fMaxTracksPerStep) {
             // it's probably better to just set meanPhotonsPerSecondary to
             // a big number if you want a small number of secondaries, but
             // this feature is retained for backwards compatibility.
             weight= weight * numSecondaries/fMaxTracksPerStep;
             numSecondaries= fMaxTracksPerStep;
     }
  }

  // if there are no photons, then we're all done now
  if (numSecondaries <= 0) {
          // return unchanged particle and no secondaries
          fAParticleChange.SetNumberOfSecondaries(0);
          goto PostStepDoIt_DONE;
  }

  // Okay, we will make at least one secondary.
  // Notify the proper authorities.
  fAParticleChange.SetNumberOfSecondaries(numSecondaries);
  if (!flagReemission)
     if (aTrack.GetTrackStatus() == fAlive)
             fAParticleChange.ProposeTrackStatus(fSuspend);

  // now look up waveform information we need to add the secondaries
  G4PhysicsOrderedFreeVector* WaveformIntegral = physicsEntry->fTimeIntegral;
  G4PhysicsOrderedFreeVector* ReemitWaveformIntegral = physicsEntry->fReemissionTimeIntegral;
  G4std::vector<G4PhysicsOrderedFreeVector*> ReemitVector = physicsEntry->fReemissionTimeVector;
  G4std::vector<G4PhysicsOrderedFreeVector*> ReemitSpectrumVector = physicsEntry->fReemissionSpectrumVector;

  for (G4int iSecondary = 0; iSecondary < numSecondaries; iSecondary++) {
      // Determine photon momentum
            G4double sampledMomentum;
           if ( !flagReemission ){
              // normal scintillation
              G4double CIIvalue = G4UniformRand()* ScintillationIntegral->GetMaxValue();
              sampledMomentum = ScintillationIntegral->GetEnergy(CIIvalue);
#ifdef RATVERBOSE
              RAT::debug << "GLG4Scint: sampledMomentum = " << sampledMomentum << newline;
              RAT::debug << "GLG4Scint: CIIvalue =        " << CIIvalue << newline;
#endif
           }
           else { // reemission
              G4bool this_is_REALLY_STUPID;
              G4double CIIvalue;
              if(num_comp){
                // Get reemission spectrum for each component
                CIIvalue = G4UniformRand() * ReemitSpectrumVector[absorbed]->
                GetValue(aTrack.GetKineticEnergy(), this_is_REALLY_STUPID );
                if (CIIvalue == 0.0) {
                  // return unchanged particle and no secondaries
                  fAParticleChange.SetNumberOfSecondaries(0);
                  goto PostStepDoIt_DONE;
                }

                sampledMomentum= ReemitSpectrumVector[absorbed]->GetEnergy(CIIvalue);
              }
              else{
                // Get reemission spectrum for a single component
                CIIvalue = G4UniformRand() * ReemissionIntegral->
                GetValue(aTrack.GetKineticEnergy(), this_is_REALLY_STUPID );
                if (CIIvalue == 0.0) {
                  // return unchanged particle and no secondaries
                  fAParticleChange.SetNumberOfSecondaries(0);
                  goto PostStepDoIt_DONE;
                }

                sampledMomentum= ReemissionIntegral->GetEnergy(CIIvalue);
              }
              fAParticleChange.ProposeLocalEnergyDeposit
              ( aTrack.GetKineticEnergy() - sampledMomentum );
              if (sampledMomentum > aTrack.GetKineticEnergy()) {
                goto PostStepDoIt_DONE;
        }
#ifdef RATDEBUG
              if (sampledMomentum > aTrack.GetKineticEnergy()) {
                RAT::warn << "GLG4Scint: Error in GLG4Scint: sampled reemitted photon momentum "
                    << sampledMomentum
                    << " is greater than track energy "
                    << aTrack.GetKineticEnergy() << newline;
              }
              if (fVerboseLevel>1) {
                 RAT::debug << "GLG4Scint: oldMomentum = " <<aTrack.GetKineticEnergy() << newline
                            << "GLG4Scint: reemittedSampledMomentum = " << sampledMomentum << newline
                            << "GLG4Scint: reemittedCIIvalue =        " << CIIvalue << newline;
              }
#endif
           }
           // Generate random photon direction
           G4double cost = 1. - 2.*G4UniformRand();
           G4double sint = sqrt(1.-cost*cost);  // FIXED BUG from G4Scint

           G4double phi = twopi*G4UniformRand();
           G4double sinp = sin(phi);
           G4double cosp = cos(phi);

           G4double px = sint*cosp;
           G4double py = sint*sinp;
           G4double pz = cost;

           // Create photon momentum direction vector
           G4ParticleMomentum photonMomentum(px, py, pz);

           // Determine polarization of new photon
           G4double sx = cost*cosp;
           G4double sy = cost*sinp;
           G4double sz = -sint;

           G4ThreeVector photonPolarization(sx, sy, sz);

           G4ThreeVector perp = photonMomentum.cross(photonPolarization);

           phi = twopi*G4UniformRand();
           sinp = sin(phi);
           cosp = cos(phi);

           photonPolarization = cosp * photonPolarization + sinp * perp;

           photonPolarization = photonPolarization.unit();

           // Generate a new photon:
            G4DynamicParticle* aScintillationPhoton = new G4DynamicParticle  (G4OpticalPhoton::OpticalPhoton(),        photonMomentum);
           aScintillationPhoton->SetPolarization (photonPolarization.x(),
                                                  photonPolarization.y(),
                                                  photonPolarization.z());

           aScintillationPhoton->SetKineticEnergy(sampledMomentum);

           // Generate new G4Track object:
           G4ThreeVector aSecondaryPosition;
           G4double deltaTime;
           if (flagReemission) {
              deltaTime= pPostStepPoint->GetGlobalTime() - t0;
              aSecondaryPosition= pPostStepPoint->GetPosition();
           }
           else {
              G4double delta= G4UniformRand() * aStep.GetStepLength();
              aSecondaryPosition = x0 + delta * p0;

              // start deltaTime based on where on the track it happened
              deltaTime = delta /
              ((pPreStepPoint->GetVelocity()+ pPostStepPoint->GetVelocity())/2.);
           }

           // delay for scintillation time

           if ( flagReemission ) {
              if(num_comp){
                // Get the re-emission timing of the absorbing component
                G4double WFvalue = G4UniformRand() * ReemitVector[absorbed]->GetMaxValue();
                G4double sampledDelayTime = ReemitVector[absorbed]->GetEnergy(WFvalue);
                deltaTime += sampledDelayTime;
              }
              else{
                // Re-emission timing for a single components
                G4double WFvalue = G4UniformRand() * ReemitWaveformIntegral->GetMaxValue();
                G4double sampledDelayTime = ReemitWaveformIntegral->GetEnergy(WFvalue);
                deltaTime += sampledDelayTime;
              }
           }
           else if ( WaveformIntegral ) {
              // Get the scintillation timing
              G4double WFvalue = G4UniformRand() * WaveformIntegral->GetMaxValue();
              G4double sampledDelayTime = WaveformIntegral->GetEnergy(WFvalue);
              deltaTime += sampledDelayTime;
           }

           // set secondary time
           G4double aSecondaryTime = t0 + deltaTime;

           // create secondary track
           G4Track* aSecondaryTrack = new G4Track(aScintillationPhoton,
                                                 aSecondaryTime,
                                                 aSecondaryPosition);
           // Add the information to the track history (if it exists)
           RAT::UserTrackInformation *trackHistory = new RAT::UserTrackInformation();
           // Copy the history of the parent into this one
           if (aTrack.GetParticleDefinition() == G4OpticalPhoton::OpticalPhoton()) {
             // This new track is produced from an existing photon.
             // Copy the parent history into it
             trackHistory->ImportHistory(dynamic_cast<RAT::UserTrackInformation *>(aTrack.GetUserInformation()));
           }
           aSecondaryTrack->SetWeight( weight );
           aSecondaryTrack->SetParentID(aTrack.GetTrackID());
           if (flagReemission){
              if(num_comp){
                // If multi-component, get absorbing component and assign
                // the creator process "Reemission_from_comp.."
                std::list<GLG4DummyProcess*>::iterator it;
                it = fReemissionProcessVector.begin();
                advance(it,absorbed);
                aSecondaryTrack->SetCreatorProcess( (*it ) );
              }
              else
                aSecondaryTrack->SetCreatorProcess( &fReemissionProcess );
           }
           else
              aSecondaryTrack->SetCreatorProcess( &fScintProcess );

           // Add the information to the track history (if it exists)
           trackHistory->AddProcessToHistory(aSecondaryTrack);

           aSecondaryTrack->SetUserInformation(trackHistory);
           // add the secondary to the ParticleChange object
           fAParticleChange.SetSecondaryWeightByProcess( true ); // recommended
           fAParticleChange.AddSecondary(aSecondaryTrack);

           aSecondaryTrack->SetWeight( weight );
           // The above line is necessary because AddSecondary() overrides
           // our setting of the secondary track weight, in Geant4.3.1 & earlier.
           // (and also later, at least until Geant4.7 (and beyond?)
           //  -- maybe not required if SetWeightByProcess(true) called,
           //  but we do both, just to be sure)
  }
    // done iSecondary loop
}
  PostStepDoIt_DONE:
  if( !flagReemission )
    fsScintillatedCount += fAParticleChange.GetNumberOfSecondaries();
  else
    fsReemittedCount += fAParticleChange.GetNumberOfSecondaries();
#ifdef RATDEBUG
  timer.Stop();
  GLG4Scint_tottime += timer.GetUserElapsed();
  GLG4Scint_num_phots += fAParticleChange.GetNumberOfSecondaries();
#endif
#ifdef RATVERBOSE
  if (fVerboseLevel>1) {
     G4cout << "\n Exiting from GLG4Scint::DoIt -- NumberOfSecondaries = "
             << fAParticleChange.GetNumberOfSecondaries() << G4endl;
  }
#endif

  return &fAParticleChange;
}


////////////////////////////////////////////////////////////////
// the generic (static) PostPostStepDoIt
G4VParticleChange*
GLG4Scint::
GenericPostPostStepDoIt(const G4Step *pStep)
{
  G4Track *track= pStep->GetTrack();
  G4double mass= track->GetDynamicParticle()->GetMass();
  G4std::vector<GLG4Scint *>::iterator it = fMasterVectorOfGLG4Scint.begin();
  for (int i= fMasterVectorOfGLG4Scint.size(); (i--) > 1; ) {
      it++;
      if ( mass < (*it)->fLowerMassLimit ) {
              return (*(--it))->PostPostStepDoIt(*track,*pStep);
      }
  }
  return (*it)->PostPostStepDoIt(*track, *pStep);
}



////////////////////////////////////////////////////////////////
// build physics tables for the scintillation process
// --------------------------------------------------
//

static G4PhysicsOrderedFreeVector*
Integrate_MPV_to_POFV( G4MaterialPropertyVector* inputVector )
{
    G4PhysicsOrderedFreeVector *aPhysicsOrderedFreeVector
      = new G4PhysicsOrderedFreeVector();

    // Retrieve the first intensity point in vector
    // of (photon momentum, intensity) pairs

    unsigned int i = 0;
    G4double currentIN = (*inputVector)[i];

    if (currentIN >= 0.0)
      {

        // Create first (photon momentum, Scintillation
        // Integral pair

        G4double currentPM = inputVector->Energy(i);

        G4double currentCII = 0.0;

        aPhysicsOrderedFreeVector->
          InsertValues(currentPM , currentCII);

        // Set previous values to current ones prior to loop
        G4double prevPM  = currentPM;
        G4double prevCII = currentCII;
        G4double prevIN  = currentIN;

        // loop over all (photon momentum, intensity)
        // pairs stored for this material
        while(i < inputVector->GetVectorLength()-1)
          {
              i++;
              currentPM = inputVector->Energy(i);
              currentIN=(*inputVector)[i];

              currentCII = 0.5 * (prevIN + currentIN);

              currentCII = prevCII +
                (currentPM - prevPM) * currentCII;

              aPhysicsOrderedFreeVector->
                InsertValues(currentPM, currentCII);

              prevPM  = currentPM;
              prevCII = currentCII;
              prevIN  = currentIN;
          }

      }

    return aPhysicsOrderedFreeVector;
}

////////////////////////////////////////////////////////////////
// MyPhysicsTable (nested class) definitions
////////////////////////////////////////////////////////////////

////////////////
// "static" members of the class
// [N.B. don't use "static" keyword here, because it means something
//  entirely different in this context.]
////////////////

GLG4Scint::MyPhysicsTable*
GLG4Scint::MyPhysicsTable::fHead = NULL;


////////////////
// constructor
////////////////

GLG4Scint::MyPhysicsTable::MyPhysicsTable()
{
  fName=0; fNext=0; fUsedByCount=0; fData=0; fLength=0;
}

////////////////
// destructor
////////////////

GLG4Scint::MyPhysicsTable::~MyPhysicsTable()
{
  if (fUsedByCount != 0) {
     RAT::warn << "GLG4Scint: Error, GLG4Scint::MyPhysicsTable is being deleted with "
          "used_by_count=" << fUsedByCount << newline;
     return;
  }
  if (fName) delete fName;
  if (fData) delete[] fData;
}


////////////////
// member functions
////////////////

void GLG4Scint::MyPhysicsTable::Dump(void) const
{
  G4cout << " GLG4Scint::MyPhysicsTable {\n"
            "  fName=" << (*fName) << G4endl
         << "  fLength=" << fLength << G4endl
         << "  fUsedByCount=" << fUsedByCount << G4endl;
  for (G4int i=0; i<fLength; i++) {
      G4cout << "  data[" << i << "]= { // "
             << (*G4Material::GetMaterialTable())[i]->GetName() << G4endl;
      G4cout << "   spectrumIntegral=";
      if (fData[i].fSpectrumIntegral)
          (fData[i].fSpectrumIntegral)->DumpValues();
      else
          G4cout << "NULL" << G4endl;

      G4cout << "   reemissionIntegral=";
      if (fData[i].fReemissionIntegral)
          (fData[i].fReemissionIntegral)->DumpValues();
      else
          G4cout << "NULL" << G4endl;

      G4cout << "   timeIntegral=";
      if (fData[i].fTimeIntegral)
          (fData[i].fTimeIntegral)->DumpValues();
      else
          G4cout << "NULL" << G4endl;

      G4cout << "   resolutionScale=" << fData[i].fResolutionScale
             << "   birksConstant=" << fData[i].fBirksConstant
             << "   ref_dE_dx=" << fData[i].fRefdEdx << G4endl
             << "   light yield=" << fData[i].fLightYield << G4endl;

      G4cout << "Quenching = ";
      if (fData[i].fQuenchingArray!=NULL)
          fData[i].fQuenchingArray->DumpValues();
      else
          G4cout << "NULL" << G4endl << "  }\n";
  }
  G4cout << " }\n";

}

GLG4Scint::MyPhysicsTable *
GLG4Scint::MyPhysicsTable::FindOrBuild(const G4String& name)
{
  // head should always exist and should always be the default (name=="")
  if (fHead == NULL) {
     fHead= new MyPhysicsTable;
     fHead->Build("");
  }

  MyPhysicsTable *rover= fHead;
  while (rover) {
         if ( name == *(rover->fName) )
            return rover;
        rover= rover->fNext;
  }

  rover= new MyPhysicsTable;
  rover->Build(name);
  rover->fNext= fHead->fNext;  // always keep head pointing to default
  fHead->fNext= rover;

  return rover;
}

void GLG4Scint::MyPhysicsTable::Build(const G4String& newname)
{
  if (fName) delete fName;
  if (fData) delete[] fData;

  fName= new G4String(newname);

  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
  fLength = G4Material::GetNumberOfMaterials();

  fData = new Entry [fLength];

  // create new physics tables

  for (G4int i=0 ; i < fLength; i++) {

      // look for material properties table entry.
      const G4Material* aMaterial = (*theMaterialTable)[i];

      // ask data[i] to Build itself
      fData[i].Build(*fName, i, aMaterial->GetMaterialPropertiesTable() );

  }
}

////////////////
// constructor for Entry
////////////////

GLG4Scint::MyPhysicsTable::Entry::Entry()
{
  fSpectrumIntegral= fReemissionIntegral = fTimeIntegral=
                                       fReemissionTimeIntegral= NULL;
  fOwnSpectrumIntegral= fOwnTimeIntegral= fOwnReemissionTimeIntegral= 0;
  fResolutionScale= 1.0;
  fLightYield=0.0;
  fBirksConstant= fRefdEdx= 0.0;
}

////////////////
// destructor for Entry
////////////////

GLG4Scint::MyPhysicsTable::Entry::~Entry()
{
  if (fSpectrumIntegral && fOwnSpectrumIntegral) {
     delete fSpectrumIntegral;
     delete fReemissionIntegral;
  }
  if (fTimeIntegral && fOwnTimeIntegral) delete fTimeIntegral;
  if (fReemissionTimeIntegral && fOwnReemissionTimeIntegral) delete fReemissionTimeIntegral;
}

////////////////
// Build for Entry
////////////////

void GLG4Scint::MyPhysicsTable::Entry::Build(const G4String& name,
                                             int material_index,
                                             G4MaterialPropertiesTable
                                             *aMaterialPropertiesTable)
{
  // delete old data, if any
  if (fSpectrumIntegral && fOwnSpectrumIntegral) {
     delete fSpectrumIntegral;
     delete fReemissionIntegral;
  }
  if (fTimeIntegral && fOwnTimeIntegral) delete fTimeIntegral;
  if (fReemissionTimeIntegral && fOwnReemissionTimeIntegral) delete fReemissionTimeIntegral;

  // set defaults
  fSpectrumIntegral= fReemissionIntegral= fTimeIntegral= fReemissionTimeIntegral= NULL;
  fResolutionScale= 1.0;
  fBirksConstant= fRefdEdx= 0.0;
  fLightYield=0.0;
  fQuenchingArray=NULL;
  // exit, leaving default values, if no material properties

  if (!aMaterialPropertiesTable)
     return;

  // Retrieve vector of scintillation wavelength intensity
  // for the material from the material's optical
  // properties table ("SCINTILLATION")
  std::stringstream property_string;

  // strncpy(property_string, ("SCINTILLATION"+name).c_str(),
  //         sizeof(property_string));
  property_string.str("");
  property_string << "SCINTILLATION" << name;
  G4MaterialPropertyVector* theScintillationLightVector=
      aMaterialPropertiesTable->GetProperty((property_string.str()).c_str());
  property_string.str("");
  property_string << "SCINTILLATION_WLS" << name;
  G4MaterialPropertyVector* theReemissionLightVector=
      aMaterialPropertiesTable->GetProperty((property_string.str()).c_str());
  if (theScintillationLightVector && !theReemissionLightVector) {
    RAT::warn << "GLG4Scint: No reemission spectrum for " << (name==""?"material":name) << ", "
              << "setting equal to primary scintillation spectrum" << newline;
     theReemissionLightVector=theScintillationLightVector;
  }

  double num_components=0;
  if(aMaterialPropertiesTable->ConstPropertyExists("NUM_COMP"))
    num_components = aMaterialPropertiesTable->GetConstProperty("NUM_COMP");

  if(num_components)
    for(int cnt=0; cnt<num_components; cnt++){
      property_string.str("");
      property_string << "SCINTILLATION_WLS" << cnt;
      G4MaterialPropertyVector* reemissionVector= aMaterialPropertiesTable->GetProperty((property_string.str()).c_str());
      if(reemissionVector){
        G4PhysicsOrderedFreeVector* tempVector= Integrate_MPV_to_POFV( reemissionVector );
        fReemissionSpectrumVector.push_back(tempVector);
      }
    }

  if (theScintillationLightVector) {
     if (aMaterialPropertiesTable->ConstPropertyExists("LIGHT_YIELD"))
        fLightYield=aMaterialPropertiesTable->GetConstProperty("LIGHT_YIELD");
     else {
        RAT::warn << "GLG4Scint: No light yield parameter for " << (name==""?"material":name) << ", "
                  << "assuming it is implicit in the scintillation integral" << newline;
     }
     // find the integral

     fSpectrumIntegral = Integrate_MPV_to_POFV( theScintillationLightVector );
     fReemissionIntegral = Integrate_MPV_to_POFV( theReemissionLightVector );
     fOwnSpectrumIntegral= 1;

  }
  else {
     // use default integral (possibly null)
     if (aMaterialPropertiesTable->ConstPropertyExists("LIGHT_YIELD"))
       fLightYield=aMaterialPropertiesTable->GetConstProperty("LIGHT_YIELD");
     fSpectrumIntegral=MyPhysicsTable::GetDefault()->GetEntry(material_index)->fSpectrumIntegral;
     fReemissionIntegral = fSpectrumIntegral;
     fOwnSpectrumIntegral = 0;
  }

  // Retrieve vector of scintillation time profile
  // for the material from the material's optical
  // properties table ("SCINTWAVEFORM")

  //strncpy(property_string, ("SCINTWAVEFORM"+name).c_str(),
  //        sizeof(property_string));
  property_string.str("");
  property_string << "SCINTWAVEFORM" << name;
  G4MaterialPropertyVector* theWaveForm =
      aMaterialPropertiesTable->GetProperty((property_string.str()).c_str());

  double rise_time = 0.0;

  if (aMaterialPropertiesTable->ConstPropertyExists("SCINT_RISE_TIME")) {
    rise_time = aMaterialPropertiesTable->GetConstProperty("SCINT_RISE_TIME");
  }

  RAT::Log::Assert(rise_time >= 0.0, "GLG4Scint::MyPhysicsTable::Entry::Build(): "
                                     "rise time must be greater than or equal to 0.");

  if (theWaveForm) {

     // do we have time-series or decay-time data?
     if (theWaveForm->GetMinLowEdgeEnergy() >= 0.0) {
             // we have digitized waveform (time-series) data
             // find the integral

             fTimeIntegral = Integrate_MPV_to_POFV( theWaveForm );
             fOwnTimeIntegral= 1;
     }
     else {
       // we have decay-time data.
       // sanity-check user's values:
       // issue a warning if they are nonsense, but continue
       if ( theWaveForm->Energy(theWaveForm->GetVectorLength() - 1) > 0.0 ) {
                G4cerr << "GLG4Scint::MyPhysicsTable::Entry::Build():  "
                       << "SCINTWAVEFORM" << name
                       << " has both positive and negative X values.  "
                          " Undefined results will ensue!\n";
             }

             /* Set the bin width to 100 times smaller than the smallest
              * decay constant. */
             G4double mintime = -1.0*(theWaveForm->GetMaxLowEdgeEnergy());
             G4double bin_width = mintime/100;

             /* Set the maximum time for the PDF to 30 times the longest
              * decay constant. */
             G4double maxtime = -30.0*(theWaveForm->GetMinLowEdgeEnergy());
             int nbins = ((int) (maxtime/bin_width)) + 1;

             G4double *tval= new G4double[nbins];
             G4double *ival= new G4double[nbins];

             /* Set the time array, and zero out the CDF. */
             for (int i=0; i < nbins; i++) {
                 tval[i]= i*maxtime/nbins;
                 ival[i]= 0.0;
             }

             G4double ampl, decy;
             for (size_t j=0; j < theWaveForm->GetVectorLength(); j++) {
               ampl = (*theWaveForm)[j];
               decy = -theWaveForm->Energy(j);
               for (int i=0; i < nbins; i++) {
                 if (rise_time != 0.0) {
                   ival[i] += ampl*(decy*(1.0-exp(-tval[i]/decy))+rise_time*(exp(-tval[i]/rise_time)-1))/(decy-rise_time);
                 } else {
                   ival[i] += ampl*(1.0-exp(-tval[i]/decy));
                 }
               }
             }

             /* Divide the CDF by the value at the end of the array to make sure
              * it is normalized to 1. */
             for (int i=0; i < nbins; i++) {
                 ival[i] /= ival[nbins-1];
             }

             fTimeIntegral= new G4PhysicsOrderedFreeVector(tval, ival, nbins);
             fOwnTimeIntegral= 1;

             // in Geant4.0.0, G4PhysicsOrderedFreeVector makes its own copy
             // of any array passed to its constructor, so ...
             delete [] tval;
             delete [] ival;
     }
  }
  else {
     // use default integral (possibly null)
     fTimeIntegral = MyPhysicsTable::GetDefault()->GetEntry(material_index)->fTimeIntegral;
     fOwnTimeIntegral= 0;

  }

  // Next get the re-emission time of each component, if it exists
  double num_comp=0;
  if(aMaterialPropertiesTable->ConstPropertyExists("NUM_COMP"))
    num_comp = aMaterialPropertiesTable->GetConstProperty("NUM_COMP");

  if(num_comp)
    for(int cnt=0; cnt<num_comp; cnt++){
      property_string.str("");
      property_string << "REEMITWAVEFORM" << cnt;
      G4MaterialPropertyVector* theReemitWaveForm =
        aMaterialPropertiesTable->GetProperty((property_string.str()).c_str());

      if (theReemitWaveForm) {
        // do we have time-series or decay-time data?
        if (theReemitWaveForm->GetMinLowEdgeEnergy() >= 0.0) {
          // we have digitized waveform (time-series) data
          // find the integral
          fReemissionTimeVector.push_back( Integrate_MPV_to_POFV( theReemitWaveForm ));
          fOwnReemissionTimeVector= 1;
        }
        else {
          // we have decay-time data.
          // sanity-check user's values:
          // issue a warning if they are nonsense, but continue
          if (theReemitWaveForm->GetMaxLowEdgeEnergy() > 0.0) {
            G4cerr << "GLG4Scint::MyPhysicsTable::Entry::Build():  "
                   << "REEMITWAVEFORM" << name
                   << " has both positive and negative X values.  "
                      " Undefined results will ensue!\n";
          }

          G4double maxtime= -3.0*(theReemitWaveForm->GetMinLowEdgeEnergy());
          G4double mintime= -1.0*(theReemitWaveForm->GetMaxLowEdgeEnergy());
          G4double bin_width = mintime/100;
          int nbins= ((int) (maxtime/bin_width)) + 1;
          G4double *tval= new G4double[nbins];
          G4double *ival= new G4double[nbins];
          for (int ii=0; ii<nbins; ii++) {
            tval[ii]= ii*maxtime/nbins;
            ival[ii]= 0.0;
          }

          for (unsigned int j=0; j < theReemitWaveForm->GetVectorLength(); j++) {
            G4double ampl = (*theReemitWaveForm)[j];
            G4double decy = theReemitWaveForm->Energy(j);
            for (int ii=0; ii<nbins; ii++) {
              ival[ii] += ampl*(1.0-exp(tval[ii]/decy));
            }
          }

          for (int ii=0; ii<nbins; ii++) {
            ival[ii]/= ival[nbins-1];
          }

          fReemissionTimeIntegral= new G4PhysicsOrderedFreeVector(tval, ival, nbins);
          fReemissionTimeVector.push_back( fReemissionTimeIntegral );
          fOwnReemissionTimeVector= 1;

          // in Geant4.0.0, G4PhysicsOrderedFreeVector makes its own copy
          // of any array passed to its constructor, so ...
          delete [] tval;
          delete [] ival;
        }
      }
      else {
        // if reemission time is not specified, set it to be the same as scintillation time
        fReemissionTimeVector.push_back(fTimeIntegral);
        fOwnReemissionTimeVector= 0;
      }
    }
  else{
    // In case of single component re-emission
    G4MaterialPropertyVector* theReemitWaveForm =
      aMaterialPropertiesTable->GetProperty("REEMITWAVEFORM");

    if (theReemitWaveForm) {
      // do we have time-series or decay-time data?
      if (theReemitWaveForm->GetMinLowEdgeEnergy() >= 0.0) {
        // we have digitized waveform (time-series) data
        // find the integral
        fReemissionTimeIntegral = Integrate_MPV_to_POFV( theReemitWaveForm );
        fOwnReemissionTimeIntegral= 1;
      }
      else {
        // we have decay-time data.
        // sanity-check user's values:
        // issue a warning if they are nonsense, but continue
        if (theReemitWaveForm->GetMaxLowEdgeEnergy() > 0.0) {
          G4cerr << "GLG4Scint::MyPhysicsTable::Entry::Build():  "
                 << "REEMITWAVEFORM" << name
                 << " has both positive and negative X values.  "
                    " Undefined results will ensue!\n";
        }

        G4double maxtime= -3.0*(theReemitWaveForm->GetMinLowEdgeEnergy());
        G4double mintime= -1.0*(theReemitWaveForm->GetMaxLowEdgeEnergy());
        G4double bin_width = mintime/100;
        int nbins= ((int) (maxtime/bin_width)) + 1;
        G4double *tval= new G4double[nbins];
        G4double *ival= new G4double[nbins];
        for (int ii=0; ii<nbins; ii++) {
          tval[ii]= ii*maxtime/nbins;
          ival[ii]= 0.0;
        }

        for (unsigned int j=0; j < theReemitWaveForm->GetVectorLength(); j++) {
          G4double ampl = (*theReemitWaveForm)[j];
          G4double decy = theReemitWaveForm->Energy(j);
          for (int ii=0; ii<nbins; ii++) {
            ival[ii] += ampl*(1.0-exp(tval[ii]/decy));
          }
        }

        for (int ii=0; ii<nbins; ii++) {
          ival[ii]/= ival[nbins-1];
        }

        fReemissionTimeIntegral= new G4PhysicsOrderedFreeVector(tval, ival, nbins);
        fOwnReemissionTimeIntegral= 1;

        // in Geant4.0.0, G4PhysicsOrderedFreeVector makes its own copy
        // of any array passed to its constructor, so ...
        delete [] tval;
        delete [] ival;
      }
    }
    else {
      // if reemission time is not specified, set it to be the same as scintillation time
      fReemissionTimeIntegral = fTimeIntegral;
      fOwnReemissionTimeIntegral= 0;
    }
  }

  // Retrieve vector of scintillation "modifications"
  // for the material from the material's optical
  // properties table ("SCINTMOD")
  //strncpy(property_string, ("SCINTMOD"+name).c_str(),
  //        sizeof(property_string));
  property_string.str("");
  property_string << "SCINTMOD" << name;
  G4MaterialPropertyVector* theScintModVector =
     aMaterialPropertiesTable->GetProperty((property_string.str()).c_str());
  if (theScintModVector == NULL) {
     // use default if not particle-specific value given
     theScintModVector = aMaterialPropertiesTable->GetProperty("SCINTMOD");
  }

  if (theScintModVector) {

     // parse the entries in ScintMod
     //  ResolutionScale= ScintMod(0);
     //  BirksConstant= ScintMod(1);
     //  Ref_dE_dx= ScintMod(2);

    for (unsigned int jscint=0; jscint < theScintModVector->GetVectorLength(); jscint++) {
      G4double key = theScintModVector->Energy(jscint);
      G4double value = (*theScintModVector)[jscint];

          if (key == 0.0) {
             fResolutionScale= value;
          }
          else if (key == 1.0) {
             fBirksConstant= value;
          }
          else if (key == 2.0) {
             fRefdEdx= value;
          }
          else {
               G4cerr<< "GLG4Scint::MyPhysicsTable::Entry::Build"
                     ":  Warning, unknown key " << key
                  << "in SCINTMOD" << name << G4endl;
          }
     }
  }

  property_string.str("");
  property_string << "QF" << name;
  fQuenchingArray = aMaterialPropertiesTable->GetProperty((property_string.str()).c_str());
  //if (fQuenchingArray!=NULL)
  //   fQuenchingArray->DumpVector();

  // Create a dummy process describing reemission for each component
  // up to the max number of components in any material in RAT
  if(fReemissionProcessVector.size() < num_comp)
    for(int i=fReemissionProcessVector.size(); i<num_comp; i++){
      std::stringstream process_name;
      process_name << "Reemission_from_comp" << i;
      GLG4DummyProcess* ReemissionCompProcess = new GLG4DummyProcess((process_name.str()).c_str(), fUserDefined);
      fReemissionProcessVector.push_back(ReemissionCompProcess);
    }
}


void
GLG4Scint::SetNewValue(G4UIcommand * command, G4String newValues)
{
  G4String commandName= command -> GetCommandName();
  if (commandName == "on") {
     fDoScintillation= true;
  }
  else if (commandName == "off") {
     fDoScintillation= false;
  }
  else if (commandName == "reemission") {
     char *endptr;
     G4int i= strtol((const char *)newValues, &endptr, 0);
     if (*endptr != '\0')  { // non-numerical argument
       if ( !(i = strcmp((const char *)newValues, "on"))) {
          fDoReemission = true;
       }
       else if ( !(i = strcmp((const char *)newValues, "off"))) {
          fDoReemission = false;
       }
       else {
          G4cerr << "Command /glg4scint/reemission given unknown parameter "
                 << '\"' << newValues << '\"' << G4endl
                 << "  old value unchanged: "
                 << ( fDoReemission ? "on" : "off" ) << G4endl;
       }
     }
     else {
        fDoReemission= (i != 0);
     }
  }
  else if (commandName == "maxTracksPerStep") {
     G4int i= strtol((const char *)newValues, NULL, 0);
     if (i > 0) {
        fMaxTracksPerStep= i;
     }
     else {
        G4cerr << "Value must be greater than 0, old value unchanged" << G4endl;
     }
  }
  else if (commandName == "meanPhotonsPerSecondary") {
     G4double d= strtod((const char *)newValues, NULL);
     if (d >= 1.0) {
        fMeanPhotonsPerSecondary= d;
     }
     else {
        G4cerr << "Value must be >= 1.0, old value unchanged" << G4endl;
     }
  }
  else if (commandName == "verbose") {
     fVerboseLevel= strtol((const char *)newValues, NULL, 0);
  }
  else if (commandName == "dump") {
     G4std::vector<GLG4Scint *>::iterator it =fMasterVectorOfGLG4Scint.begin();
     for (; it != fMasterVectorOfGLG4Scint.end(); it++) {
         (*it)->DumpInfo();
     }
  }
  else if (commandName=="setQF"){
    G4double d= strtod((const char *)newValues, NULL);
    if (d<=1.0) {
       SetQuenchingFactor(d);
       fUserQF=true;
    }
    else {
       G4cerr << "The quenching factor is <= 1.0, old value unchanged" << G4endl;
    }
  }
  else {
     G4cerr << "No GLG4Scint command named " << commandName << G4endl;
  }
  return;
}

G4String
GLG4Scint::GetCurrentValue(G4UIcommand * command)
{
  G4String commandName= command -> GetCommandName();

  if (commandName == "on" || commandName == "off") {
     return fDoScintillation ? "on" : "off";
  }
  else if (commandName == "reemission") {
     return fDoReemission ? "1" : "0";
  }
  else if (commandName == "maxTracksPerStep") {
     char outbuff[64];
     sprintf(outbuff, "%d", fMaxTracksPerStep);
     return G4String(outbuff);
  }
  else if (commandName == "meanPhotonsPerSecondary") {
     char outbuff[64];
     sprintf(outbuff, "%g", fMeanPhotonsPerSecondary);
     return G4String(outbuff);
  }
  else if (commandName == "verbose") {
     char outbuff[64];
     sprintf(outbuff, "%d", fVerboseLevel);
     return G4String(outbuff);
  }
  else if (commandName == "dump") {
     return "?/glg4scint/dump not supported";
  }
  else if(commandName=="setQF"){
     char outbuff[64];
     sprintf(outbuff, "%g", GetQuenchingFactor());
     return G4String(outbuff);
  }
  else {
     return (commandName+" is not a valid GLG4Scint command");
  }
}