#include <CLHEP/Units/SystemOfUnits.h>
using namespace CLHEP;
#include <RAT/SNOMANCerenkov.hh>

// BuildThePhysicsTable for the Cerenkov process
// ---------------------------------------------
//

void SNOMANCerenkov::BuildThePhysicsTable()
{
  if (thePhysicsTable) return;

  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
  // create new physics table
  thePhysicsTable = new G4PhysicsTable(numOfMaterials);

  // loop for materials
  for (G4int i=0 ; i < numOfMaterials; i++)
	{
      G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();

      // Retrieve vector of refraction indices for the material
      // from the material's optical properties table
      G4Material* aMaterial = (*theMaterialTable)[i];

      G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();

      if (aMaterialPropertiesTable)
        {
          G4MaterialPropertyVector* theRefractionIndexVector = aMaterialPropertiesTable->GetProperty("RINDEX");

          if (theRefractionIndexVector)
            {
              G4double meanRI = 0.0;
              if( aMaterialPropertiesTable->ConstPropertyExists( "RINDEX_MEAN" ) )
                meanRI = aMaterialPropertiesTable->GetConstProperty( "RINDEX_MEAN" );
              else // Must calculate the average
                {
                  for( size_t iR = 0; iR < theRefractionIndexVector->GetVectorLength(); iR++ )
                    meanRI += (*theRefractionIndexVector)[iR];
                  meanRI /= static_cast<G4double>( theRefractionIndexVector->GetVectorLength() );
                }

			 // Create first (photon energy, Cerenkov Integral)
			 // pair

              G4double currentPM = theRefractionIndexVector->Energy(0);
              G4double currentCAI = 0.0;

              aPhysicsOrderedFreeVector->InsertValues(currentPM , currentCAI);

              // Set previous values to current ones prior to loop

              G4double prevPM  = currentPM;
              G4double prevCAI = currentCAI;
              G4double prevRI  = meanRI;

              // loop over all (photon energy, refraction index)
              // pairs stored for this material

              for (size_t i = 1; i < theRefractionIndexVector->GetVectorLength(); i++)
                {
                  currentPM = theRefractionIndexVector->Energy(i);

                  currentCAI = 0.5*(1.0/(prevRI*prevRI) + 1.0/(meanRI*meanRI));

                  currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;

                  aPhysicsOrderedFreeVector->InsertValues(currentPM, currentCAI);

                  prevPM  = currentPM;
                  prevCAI = currentCAI;
                  prevRI  = meanRI;
                }

            }
        }

      // The Cerenkov integral for a given material
      // will be inserted in thePhysicsTable
      // according to the position of the material in
      // the material table.

      thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector);
    }
}

G4double
SNOMANCerenkov::GetAverageNumberOfPhotons(const G4double charge,
                                          const G4double beta,
                                          const G4Material* aMaterial,
                                          G4MaterialPropertyVector* Rindex) const
{
  const G4double Rfact = 369.81/(eV * cm);

  if(beta <= 0.0)return 0.0;

  G4double BetaInverse = 1./beta;

  // Vectors used in computation of Cerenkov Angle Integral:
  //  - Refraction Indices for the current material
  //  - new G4PhysicsOrderedFreeVector allocated to hold CAI's
  G4int materialIndex = aMaterial->GetIndex();

  // Retrieve the Cerenkov Angle Integrals for this material
  G4PhysicsOrderedFreeVector* CerenkovAngleIntegrals =
    (G4PhysicsOrderedFreeVector*)((*thePhysicsTable)(materialIndex));

  if(!(CerenkovAngleIntegrals->IsFilledVectorExist()))return 0.0;

  // Min and Max photon energies
  G4double Pmin = Rindex->GetMinLowEdgeEnergy();
  G4double Pmax = Rindex->GetMaxLowEdgeEnergy();

  // Max Cerenkov Angle Integral
  G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();

  G4double dp, ge;

  G4double nAverage = 0.0;
  for (size_t iR = 0; iR < Rindex->GetVectorLength();iR++)
    nAverage += (*Rindex)[iR];
  nAverage /= static_cast<G4double>( Rindex->GetVectorLength() );

  // If n_average < 1/Beta -- no photons generated
  if (nAverage < BetaInverse) {
    dp = 0;
    ge = 0;
  }
  // otherwise photons generated
  else {
    dp = Pmax - Pmin;
    ge = CAImax;
  }

  // Calculate number of photons
  G4double NumPhotons = Rfact * charge/eplus * charge/eplus *
    (dp - ge * BetaInverse*BetaInverse);
  return NumPhotons;
}