RPCChargeGenerator.cc

Go to the documentation of this file.

#include "RPCChargeGenerator.h"
#include <cmath>
#include <cfloat>
#include <utl/AugerUnits.h>
#include <utl/Math.h>
#include <utl/ErrorLogger.h>
#include <fwk/RandomEngineRegistry.h>
#include <CLHEP/Random/Randomize.h>
#include <boost/math/special_functions/gamma.hpp>

using namespace std;
using namespace fwk;
using namespace utl;
using namespace RPCSimulatorLX;
using CLHEP::RandGeneral;

RPCChargeGenerator::RPCChargeGenerator(double k0, double theta0, double referenceAngle) :
  fK0(k0),
  fTheta0(theta0),
  fReferenceAngle(referenceAngle),
  fMinSecTheta(1.0),
  fMaxSecTheta(11.0),
  fNumberOfChargeBins(100),
  fNumberOfChargePDF(100),
  fListChargePDF(0),
  fListMinCharge(0),
  fListMaxCharge(0),
  fIsInitialized(false),
  fRandomEngine(0)
{  
  Initialize();
}


RPCChargeGenerator::~RPCChargeGenerator()
{;}


double 
RPCChargeGenerator::GetCharge(double angle)
{

  double charge = 0.0;

  if (fIsInitialized) {

    double secTheta = 1.0/abs(cos(angle));

    if (secTheta > fMaxSecTheta)
      secTheta = fMaxSecTheta;
    else if (secTheta < fMinSecTheta)
      secTheta = fMinSecTheta;

    unsigned int iChargePDF = 
      int(rint((secTheta-fMinSecTheta)*double(fNumberOfChargePDF-1)/(fMaxSecTheta-fMinSecTheta)));

    double minCharge = fListMinCharge[iChargePDF];
    double maxCharge = fListMaxCharge[iChargePDF];

    charge = 
      minCharge + (maxCharge - minCharge)*(fListChargePDF[iChargePDF])->fire();
    
    charge *= utl::pico*coulomb;
  }
  else
    ERROR("Error in RPCChargeGenerator::GetCharge: Parameters not initialized !!!");

  return charge;
}


void RPCChargeGenerator::Initialize(void)
{
  if (fIsInitialized) {
    ERROR("RPCChargeGenerator already initialized !");
  }
  else {
    double secTheta0 = 1.0/abs(cos(fReferenceAngle));

    for (unsigned int i = 0 ; i < fNumberOfChargePDF; i++) {
      
      double secTheta = fMinSecTheta + 
        double(i)*(fMaxSecTheta-fMinSecTheta)/double(fNumberOfChargePDF-1); 
      
      //    cout << "sec theta: " << secTheta << "  theta: "  << acos(1.0/secTheta)/deg << endl;
      double k = fK0 * secTheta/secTheta0;
      double theta = fTheta0;
      
      fRandomEngine = 
        &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);
      
      // The maximum charge for the pdf is defined to contain 99.9% of the integral
      double maxCharge = InverseCumulativeGammaDistribution(0.999,k,theta);;
      double minCharge = 0.0;    
      
      //     cout << "max charge/integral " << maxCharge << "/" 
      //         << CumulativeGammaDistribution(maxCharge,k,theta)
      //         << endl;
      
      vector<double> chargeProbability;
      
      for (unsigned int jCharge = 0 ; jCharge < fNumberOfChargeBins; jCharge++) {
        double charge = 
          minCharge + jCharge*(maxCharge-minCharge)/double(fNumberOfChargeBins);
        
        chargeProbability.push_back(GammaDistribution(charge,k,theta));
        
        //      cout << charge << "   " << chargeProbability[jCharge] << endl;
      }
    
    
      const unsigned int chargeProbN = chargeProbability.size();

      RandGeneral * chargeDist = 
        new RandGeneral(&fRandomEngine->GetEngine(), &chargeProbability[0], chargeProbN);
      
      fListMinCharge.push_back(minCharge);
      fListMaxCharge.push_back(maxCharge);
      fListChargePDF.push_back(chargeDist);
    }

    fIsInitialized = true;
  }

  
}