MdShowerRegenerator.cc

Go to the documentation of this file.

#include <utl/config.h>

#include <cmath>
#include <sstream>
#include <vector>
#include <limits>

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/RandomEngineRegistry.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <mevt/CounterSimData.h>

#include <det/Detector.h>
#include <mdet/MDetector.h>
#include <mdet/Counter.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>

#include <utl/ConfigUtils.h>
#include <utl/ErrorLogger.h>
#include <utl/GeometryUtilities.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/ParticleCases.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/RandomEngine.h>
#include <utl/TimeStamp.h>
#include <utl/Math.h>
#include <utl/Reader.h>
#include <utl/TabularStream.h>
#include <utl/TabulatedFunction.h>
//#include <utl/UTMPoint.h>
//#include <utl/ReferenceEllipsoid.h>

#include <CLHEP/Random/RandFlat.h>
#include <CLHEP/Random/RandPoisson.h>
#include <CLHEP/Random/RandGauss.h>

#include "MdShowerRegenerator.h"

using namespace MdShowerRegeneratorAG;
using namespace fwk;
using namespace evt;
using namespace mevt;
using namespace utl;
using namespace std;

using CLHEP::RandFlat;
using CLHEP::RandPoisson;

//#warning "MD_THETASIM_MAX - MD_LARGE_SIDE are brutally hard coded!!!"
#define MD_THETASIM_MAX 80*utl::degree
#define MD_SIM_AREA 30*utl::m2
#define MD_LARGE_SIDE 12*utl::m
 
namespace MdShowerRegeneratorAG {

struct InStats {
  int InRad;
  int InPhi;
  int IsIn;
};


  inline
  double
  PlaneFrontTime(const CoordinateSystemPtr showerCS, const Point& corePosition,
                 const Point& position)
  {
    return (corePosition.GetZ(showerCS) - position.GetZ(showerCS)) / kSpeedOfLight;
  }

} // namespace MdShowerRegeneratorAG


MdShowerRegenerator::MdShowerRegenerator() :
  fInnerRadiusCut(0.),
  fOuterRadiusCut(1.e6*km),
  fElectronEnergyCut(0),
  fMuonEnergyCut(0),
  fPhotonEnergyCut(0),
  fHadronEnergyCut(0),
  fMesonEnergyCut(0),
  fDeltaROverR(0),
  fDeltaPhi(0),
  fHorizontalParticleCut(0),
  fRandomEngine(0),
  fOnlyMuons( true ),
  fMuonWeightScale(1)
{
 fNMuons=0;
 fNInjected=0;
 fNHorizontalRejected=0;
}


VModule::ResultFlag
MdShowerRegenerator::Init()
{

  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("MdShowerRegenerator");

  AttributeMap useAtt;
  useAtt["use"] = string("yes");

  Branch distanceCutB = topB.GetChild("DistanceCuts");

  Branch inB = distanceCutB.GetChild("InnerRadiusCut", useAtt);
  if (inB)
    inB.GetData(fInnerRadiusCut);

  Branch outB = distanceCutB.GetChild("OuterRadiusCut", useAtt);
  if (outB)
    outB.GetData(fOuterRadiusCut);

  Branch energyCutB = topB.GetChild("EnergyCuts");
  energyCutB.GetChild("ElectronEnergyCut").GetData(fElectronEnergyCut);
  energyCutB.GetChild("MuonEnergyCut").GetData(fMuonEnergyCut);
  energyCutB.GetChild("PhotonEnergyCut").GetData(fPhotonEnergyCut);
  energyCutB.GetChild("HadronEnergyCut").GetData(fHadronEnergyCut);
  energyCutB.GetChild("MesonEnergyCut").GetData(fMesonEnergyCut);

  Branch algoB = topB.GetChild("AlgorithmParameters");
  algoB.GetChild("DeltaROverR").GetData(fDeltaROverR);
  algoB.GetChild("DeltaPhi").GetData(fDeltaPhi);
  algoB.GetChild("HorizontalParticleCut").GetData(fHorizontalParticleCut);

  if (topB.GetChild("MuonWeightScale"))
    topB.GetChild("MuonWeightScale").GetData(fMuonWeightScale);

  fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector).GetEngine();

  //This is to avoid  verification
  LoadConfig(topB, "OnlyMuons", fOnlyMuons, true );
  
  return eSuccess;
}

bool 
MdShowerRegenerator::IsIn(double phi, double r, InStats & stat)     
{ 
   bool retphi = IsPhiIn( phi );
   bool retrad = IsRIn( r );

   bool isin   = retphi && retrad;

   if( retrad )
     ++(stat.InRad);

   if( retphi ) {
      ++(stat.InPhi);
/*
   std::cout << "Phi1 " << fPhi1/utl::degree  << "\n";
   std::cout << "Phi2 " << fPhi2/utl::degree  << "\n";
   std::cout << "phi " << phi/utl::degree  << "\n";
*/
   }
  
   if( isin )
     ++(stat.IsIn);

   return isin;
}

     
void 
MdShowerRegenerator::SetRadRange(double r1, double r2)
{
   fR1 = r1;
   fR2 = r2;
}
    
void 
MdShowerRegenerator::SetPhiRange( double phi, double dphi)
{

  fPhi1 = phi - dphi;
  fPhi2 = phi + dphi;

  if (fPhi1 < -kPi) {
    if (fPhi2 > kPi) {
      fPhi1 = -kPi;
      fPhi2 = kPi;
    } else
      fPhi1 += 2*kPi;
  } else if (fPhi2 > kPi)
    fPhi2 -= 2*kPi;
}

bool 
MdShowerRegenerator::IsPhiIn( double phi ) 
{
  if (fPhi1 < fPhi2)
    return fPhi1 <= phi && phi <= fPhi2;
  else
    return fPhi1 <= phi || phi <= fPhi2;
}
    
bool 
MdShowerRegenerator::IsRIn(const double r ) 
{ return fR1 <= r && r <= fR2; }

bool
MdShowerRegenerator::IsParticleEnergyLow(const int pType,
                                         const double pEnergy)
  const
{
  switch (pType) {
    case OFFLINE_PHOTON:
      return pEnergy < fPhotonEnergyCut;
    case OFFLINE_ELECTRONS:
      return pEnergy < fElectronEnergyCut;
    case OFFLINE_MUONS:
      return pEnergy < fMuonEnergyCut;
    case OFFLINE_BARYONS:
      return pEnergy < fHadronEnergyCut;
    case OFFLINE_MESONS:
      return pEnergy < fMesonEnergyCut;
    default:
      return true;
  }
}

bool
MdShowerRegenerator::IsMuon(const int pType)
{
  switch (pType) {
    case OFFLINE_MUONS:
      ++fNMuons;
      return true;
    default:
      return false;
  }
}


VModule::ResultFlag
MdShowerRegenerator::Run(Event& event)
{
  INFO("Performing Unthinning...");

  if (!event.HasSimShower()) {
    ERROR("Current event does not have a simulated shower.");
    return eFailure;
  }
  ShowerSimData& simShower = event.GetSimShower();

  if( !simShower.HasGroundParticles()){
    ERROR("Current event does not have a simulated shower with particles.");
    return eFailure;
  } else {
    INFO("Retrieving ground particles from shower");
  }

  if (simShower.GetMinRadiusCut() < fInnerRadiusCut) {
    ostringstream msg;
    msg << "Setting new value for fMinRadiusCut: " << fInnerRadiusCut;
    WARNING(msg);
    simShower.SetMinRadiusCut(fInnerRadiusCut);
  }
  
  //Retrieve muon detector
  const mdet::MDetector& mDetector = det::Detector::GetInstance().GetMDetector();

  if (!event.HasMEvent())
    event.MakeMEvent();

  MEvent& mEvent = event.GetMEvent();

  if(event.HasSEvent()){
    WARNING("SEvent found. InnerRadiusCut from MdShowerRegenerator.xml meaningless!");
  }
  /* Get the shower coordinate system. 
   * The origin of this coordinate system is the core position of the shower. 
   * The z axis points opposite to the direction of motion of the primary 
   * particle (i.e., its 3-momentum is anti-parallel to the unit vector in z direction). 
   * The y axis is in the horizontal plane as defined by the local coordinate system. 
   * The x axis completes the system to a right-handed coordinate system.
   */
  const CoordinateSystemPtr showerCS     = simShower.GetShowerCoordinateSystem();

   /* Get the auger coordinate system associated to the shower core position.
    * Get the local coordinate system corresponding to the shower core position. 
    * This is an Auger Coordinate System, i.e., 
    * the x-y plane is horizontal and the z axis points up.
    */
  const CoordinateSystemPtr grParticleCS = simShower.GetLocalCoordinateSystem();
 
  /* Site coordinate system */
  const CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();

  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const double showerZenith = (-simShower.GetDirection()).GetTheta(localCS); 
  const double samplingAreaFactor = 4. * fDeltaPhi * fDeltaROverR / cos(showerZenith);

  //Loop over all counters
  for (mdet::MDetector::CounterConstIterator cIt = mDetector.CountersBegin();
       cIt != mDetector.CountersEnd(); ++cIt) {

    const int mId = cIt->GetId();
    
    if (!mEvent.HasCounter(mId))
      mEvent.MakeCounter(mId);
    
    mevt::Counter& eCounter = mEvent.GetCounter(mId);

    if (!eCounter.HasSimData())
      eCounter.MakeSimData();
    
    CounterSimData& sSim = eCounter.GetSimData();
    
    //check sd if station is inside min radius
    //using the distance of Sd instead of Md to always have pairs of simulators when simulating in XD (Md + Sd), the difference is of the order of few meters and is zero for vertical showers (for now, see MModelConfig.xml)

    bool sd_exists = false;
    if(event.HasSEvent()){
      sevt::SEvent& sEvent = event.GetSEvent();
      int partnerId = cIt->GetAssociatedTankId();
      if(sEvent.HasStation(partnerId)){
        sevt::Station& sStation = sEvent.GetStation(partnerId);
        if(sStation.HasSimData()){
          sevt::StationSimData& sSimData = sStation.GetSimData();
          sd_exists = true;
          if(sSimData.IsInsideMinRadius()){
            INFO("Associated Sd-Tank was flagged InsideMinRadiusCut. Forcing the same for Md-Counter.");
            sSim.SetIsInsideMinRadius();
            continue;
          }
        }    
      }
    }
    if(!sd_exists){ //if simulating in Md-only mode, use fInnerRadiusCut from xml
      const Point& countPos = cIt->GetPosition();
      double cR = countPos.GetRho(showerCS);
      if (cR < fInnerRadiusCut) {
        sSim.SetIsInsideMinRadius();
        continue;
      }
    }
  }
  InStats stat;
  stat.InRad = 0;
  stat.InPhi = 0;
  stat.IsIn  = 0;

  ShowerParticleIterator particleIt = simShower.GroundParticlesBegin();
  for (; particleIt != simShower.GroundParticlesEnd(); ++particleIt) {

    const Point& partPos = particleIt->GetPosition();

    if (IsParticleEnergyLow(particleIt->GetType(), particleIt->GetKineticEnergy()))
      continue;
    
    if( fOnlyMuons && !IsMuon(particleIt->GetType()) )
      continue;
    
    //Loop over all counters
    for (mdet::MDetector::CounterConstIterator cIt = mDetector.CountersBegin();
         cIt != mDetector.CountersEnd(); ++cIt) {
      
      const int mId = cIt->GetId();
      
      mevt::Counter& eCounter = mEvent.GetCounter(mId);
      
      CounterSimData& sSim = eCounter.GetSimData();
      
      if(sSim.IsInsideMinRadius())
        continue;

      //GetPosition actually is from Station obj which is a friend of Counter obj
      const Point& countPos    = cIt->GetPosition();
      //      const UTMPoint utm( countPos, ReferenceEllipsoid::GetWGS84());
      
      double cR = countPos.GetRho(showerCS);
      
      if (cR > fOuterRadiusCut)
        continue;
      
      double r   = partPos.GetRho(showerCS);
      double phi = partPos.GetPhi(showerCS);

      SetPhiRange( countPos.GetPhi(showerCS), fDeltaPhi );
      SetRadRange( cR*(1.-fDeltaROverR),cR*(1.+fDeltaROverR));
      
      if( !IsIn(phi, r, stat) )
        continue;
      
      //scale the weight of muons accordingly the fMuonWeightScale parameter
      const double pWeight =
        (particleIt->GetType() == utl::Particle::eMuon ||
         particleIt->GetType() == utl::Particle::eAntiMuon?
         fMuonWeightScale*particleIt->GetWeight():
         particleIt->GetWeight());
      
      /*
      std::cout << "Counter Id: " << mId << "\n";
      std::cout << "particleIt->GetType() "  << particleIt->GetType() << "\n";
      std::cout << "Particle original weitgh " << pWeight << "\n";
      std::cout << "fPhi1 " << fPhi1/utl::degree  << "\n";
      std::cout << "fPhi2 " << fPhi2/utl::degree  << "\n";
      std::cout << "fDeltaPhi " << fDeltaPhi/utl::degree  << "\n";
      std::cout << "phi " << phi/utl::degree  << "\n";
      */
      
      Particle newParticle(*particleIt);
      //Set Parent (same particle at this level)
      newParticle.SetParent( newParticle );

      // Uses the depth of the first module (arbitratry)
      const double cDepth = cIt->ModulesBegin()->GetDepth();
      //std::cout << "cDepth " << cDepth/utl::m << " (m)" << std::endl;
      
      const Vector& pDirection  = newParticle.GetDirection();
      const double pCosine      = -1.*pDirection.GetZ(grParticleCS);
      //const double pSine        = pDirection.GetRho(grParticleCS);
      
      //Avoid upgoing or too inclined particles
      if (pCosine < 0 || pCosine < cos(MD_THETASIM_MAX)){
        /*
        std::cout << "Horizontal cut at cos(" << MD_THETASIM_MAX/utl::degree << ")=" << cos(MD_THETASIM_MAX) << std::endl;
        std::cout << "pCosine " << pCosine << std::endl;  
        */
        fNHorizontalRejected++;
        continue;
      }
      
      // particle time relative to the core time
      const double pTime = particleIt->GetTime().GetInterval();
      
      // Nominal area of counter detectors
      //const double RadiusCounter = sqrt( MD_SIM_AREA / kPi );
      const double RadiusCounter = MD_LARGE_SIDE;
      //const double ExtraRadius   = cDepth * pSine/pCosine;
      const double ExtraRadius   = cDepth * tan(MD_THETASIM_MAX);
      /*
      std::cout << "RadiusCounter  " << RadiusCounter/utl::m << " (m)" << std::endl;
      std::cout << "ExtraRadius    " << ExtraRadius/utl::m   << " (m)" << std::endl;
      std::cout << "Particle theta " << acos(pCosine)/utl::degree   << " (deg)" << std::endl;
      */
      const double samplingArea   = samplingAreaFactor * Sqr(cR);
      const double pWeightDensity = pWeight / samplingArea;
      const double effArea        = kPi * Sqr( RadiusCounter/utl::m+ExtraRadius/utl::m );
      
      const double avgN = pWeightDensity * effArea;
      
      unsigned int n(0);
      double weight(0);
      n = RandPoisson::shoot(fRandomEngine, avgN);
      /*
      std::cout << "Reweight: " << pWeightDensity << std::endl;
      std::cout << "effective area: " << effArea << std::endl;
      std::cout << "Average number of particles: " <<  avgN  << std::endl;
      std::cout << "Poisson number to propagate: " <<  n     << std::endl;
      */
      weight = 1;
      newParticle.SetWeight(weight);
      
      const CoordinateSystemPtr localCS = cIt->GetLocalCoordinateSystem();
      
      for (unsigned int i = 0; i < n; ++i) {


        //Uniform distribution on a circle 
        const double r   = (RadiusCounter+ExtraRadius) * sqrt( RandFlat::shoot(fRandomEngine, 0, 1) );
        const double phi = 2.*kPi * RandFlat::shoot(fRandomEngine, 0, 1);
        
        //Set new position (impact point)
        const Point pNewPosition(r, phi, 0, localCS, Point::kCylindrical);
        newParticle.SetPosition(pNewPosition);
        
        //Set new direction (impinging vector) 
        const Vector pNewDirection( pDirection.GetR(grParticleCS), pDirection.GetTheta(grParticleCS), pDirection.GetPhi(grParticleCS), grParticleCS, Vector::kSpherical );      
        newParticle.SetDirection(pNewDirection);
        
        //Set new time
        double pShiftedTime = pTime + PlaneFrontTime(showerCS, partPos, pNewPosition);
        newParticle.SetTime(pShiftedTime);
        /*
        std::cout << "injected radius: " << r/utl::m        << " m"   << std::endl;
        std::cout << "injected phi   : " << phi/utl::degree << " deg" << std::endl;
        std::cout << "pTime : " << pTime/utl::ns << " ns " << std::endl;
        std::cout << "pShiftedTime : " << pShiftedTime/utl::ns << " ns " << std::endl;
        */
        sSim.AddGrdParticle(newParticle);
        ++fNInjected;
        
      } // for n

    } // loop on mdet::Counters
  }//loop on shower particles
  
  std::cout << "Total number of muons in the shower                " << fNMuons << std::endl;
  std::cout << "Total number of muons injected for ground prop     " << fNInjected << std::endl;
  std::cout << "IsIn       " << stat.IsIn << std::endl;
  std::cout << "InDeltaR   " << stat.InRad   << std::endl;
  std::cout << "InDeltaPhi " << stat.InPhi << std::endl;
  std::cout << "muons rejected because of cut in zenith at " << MD_THETASIM_MAX/utl::degree<< "deg "<< fNHorizontalRejected << std::endl;

  return eSuccess;
}


VModule::ResultFlag
MdShowerRegenerator::Finish()
{
  return eSuccess;
}