MdGroundPropagator.cc

Go to the documentation of this file.

#include "MdGroundPropagator.h"
#include "PhysicsListCustomization.h"
#include "PhysicsList.h"
#include "SoilConstruction.h"
#include "PrimaryGenerator.h"
#include "SteppingAction.h"

//Geant4 Managers
#include <G4RunManager.hh>
#include <G4UImanager.hh>
#include <G4VisManager.hh>
#include <G4VisExecutive.hh>

#include <G4GeometryManager.hh>
#include <G4PhysicalVolumeStore.hh>
#include <G4LogicalVolumeStore.hh>
#include <G4SolidStore.hh>

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

#include <evt/Event.h>
#include <mevt/CounterSimData.h>
#include <mdet/Counter.h>
#include <det/Detector.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/Particle.h>
#include <utl/ShowerParticleIterator.h>
#include <utl/TimeDistribution.h>
#include <utl/TimeDistributionAlgorithm.h>
#include <utl/AugerUnits.h>

#include <cstddef>
#include <iostream>
#include <fstream>
#include <sstream>

#include <CLHEP/Random/Random.h>

#include <tls/Geant4Manager.h>

using namespace utl;
using namespace fwk;
using namespace std;
using namespace evt;
using namespace mevt;
using namespace mdet;
using namespace det;
using namespace GroundPropagatorAG;

const mdet::Counter* GroundPropagator::fCurrentDetectorCounter = 0;
MEvent::CounterIterator GroundPropagator::fCurrentEventCounterIt;
CounterSimData::GrdParticleIterator GroundPropagator::fCurrentParticleIt;

GroundPropagator::GroundPropagator() :
  fEventId(""),
  fUImanager(0),
  fVisManager(0),
  fSteppingAction(0)
{
  fGeomVisOn = false;
  fTrajVisOn = false;
  fUseGlobalPhysicsList = false;

}

GroundPropagator::~GroundPropagator() { }

VModule::ResultFlag
GroundPropagator::Init()
{
  INFO("\n++++\nReading XML parameters");
  //Top branch
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("MdGroundPropagator");

  int verbosity = 0;
  string filename;
  string physList;
  //Verbosity branch
  Branch verB = topB.GetChild("verbosity");
  verB.GetChild("level").GetData(verbosity);
#ifdef AMIGA_DEBUG
  verB.GetChild("statsfile").GetData(filename);
  fStatFile.open( filename.c_str() );
  if( !fStatFile ) {
    G4cerr << "Error opening file " << filename << G4endl;
    return eFailure;
  }
  string strfname( filename );
  size_t found = strfname.rfind ( ".ascii" );
  if( found != string::npos ) 
    strfname = strfname.substr( 0, found) + "_undergrd" + strfname.substr( found , string::npos); 
  cout << "Ascii Underground filename " << strfname << endl;
  fStatUFile.open ( strfname.c_str() , ios_base::out );
  if( ! fStatUFile.is_open() ) {
    cout << "Error opening file: " << strfname << endl;
    return eFailure;
  }
  fStatUFile << "#Id PartName Energy(MeV) PDGCode Time(ns) TrackId X Y Z Parent\n";
#endif 
  //Soil branch
  Branch grmB = topB.GetChild("soil");
  grmB.GetChild("grammage").GetData(fGrammage);
  G4cerr << "fGrammage " << fGrammage/(utl::g/utl::cm2) << G4endl;

  //Visualization branch
  Branch visB = topB.GetChild("visualization");
  visB.GetChild("geometry").GetData(fGeomVisOn);
  visB.GetChild("trajectories").GetData(fTrajVisOn);
  //Physics list branch
  Branch phyB = topB.GetChild("physics");
  phyB.GetChild("list").GetData(physList);

  Branch globalPhysicsListB = topB.GetChild("globalPhysicsList");
  globalPhysicsListB.GetData(fUseGlobalPhysicsList);


  if ( (fGeomVisOn || fTrajVisOn) && !fVisManager){
    fVisManager = new G4VisExecutive();
  }
  if (!fRunManager) {
    if (fUseGlobalPhysicsList)
      fRunManager = tls::Geant4Manager::GetInstance().GetRunManager();
    else
      fRunManager = new G4RunManager;
  }

  fUImanager  = G4UImanager::GetUIpointer();
  stringstream oss (ostringstream::out);

  oss << "/run/verbose " << verbosity;
  fUImanager->ApplyCommand( oss.str() );

  oss.str("");
  oss << "/event/verbose " << verbosity;
  fUImanager->ApplyCommand( oss.str() );

  oss.str("");
  oss << "/tracking/verbose " << verbosity; 
  fUImanager->ApplyCommand( oss.str() );

  if(fUseGlobalPhysicsList){
    tls::Geant4Manager::GetInstance().AddVisManager(fVisManager);

    INFO("using global PhysicsList from Geant4Manager");
    PhysicsListCustomization* physicsListCustomization = new PhysicsListCustomization();

    INFO("\n++++\nConstructing G4 ground");
    SoilConstruction* soil = new SoilConstruction(fGrammage);

    tls::Geant4Customization custom("MdGroundPropagator", soil, physicsListCustomization);

    INFO("\n++++\nInitializing G4 stepping action");
    fSteppingAction = new SteppingAction(soil->GetDepth());
    custom.SetSteppingAction(fSteppingAction);

    INFO("\n++++\nInitializing G4 primary generator");
    custom.SetPrimaryGenerator(new PrimaryGenerator(filename));

    tls::Geant4Manager::GetInstance().AddCustomization(custom);

  }else{
    INFO("using PhysicsList from MdGroundPropagator");
    fRunManager->SetUserInitialization(new PhysicsList);

    INFO("\n++++\nConstructing G4 ground");
    SoilConstruction* soil = new SoilConstruction(fGrammage);
    fRunManager->SetUserInitialization(soil);
    
    INFO("\n++++\nInitializing G4 stepping action");
    fSteppingAction = new SteppingAction(soil->GetDepth());
    fRunManager->SetUserAction(fSteppingAction);

    INFO("\n++++\nInitializing G4 primary generator");
    fRunManager->SetUserAction(new PrimaryGenerator(filename));

    fUImanager->ApplyCommand("/physlist/CutsAll 1 mm");

    //configure physicsList 
    oss.str("");
    oss << "/physlist/Physics " << physList;
    fUImanager->ApplyCommand( oss.str() );
  }

  if ( fGeomVisOn || fTrajVisOn ) {
    fVisManager->Initialize();    
    fUImanager->ApplyCommand("/vis/open HepRepFile");
    fUImanager->ApplyCommand("/vis/scene/create");
    fUImanager->ApplyCommand("/vis/sceneHandler/attach");
    fUImanager->ApplyCommand("/vis/scene/add/volume");
    fUImanager->ApplyCommand("/vis/drawVolume");
    fUImanager->ApplyCommand("/vis/scene/notifyHandlers");
    fUImanager->ApplyCommand("/vis/viewer/update");
  }
  if (fTrajVisOn) {
    fUImanager->ApplyCommand("/tracking/storeTrajectory 1");
    fUImanager->ApplyCommand("/vis/scene/add/trajectories rich");    
  }

  return eSuccess;
}

VModule::ResultFlag
GroundPropagator::Run(Event& event)
{
  INFO("\n++++\nStarting MdGroundPropagator run");

  if (!event.HasMEvent()) {
    ERROR("\n++++\nMEvent does not exist.");
    return eFailure;
  }

  // Get the MEvent
  MEvent& mEvent = event.GetMEvent();

  if (mEvent.GetNumberOfCounters() > 0) {
    INFO("\n++++\nActivating G4 for this run");
    if(!fUseGlobalPhysicsList){
      fRunManager->Initialize();
    } else {
      tls::Geant4Manager::GetInstance().Customize("MdGroundPropagator");
    }
  }

  // Get the MDetector
  const MDetector& mDetector = Detector::GetInstance().GetMDetector();

  for (MEvent::CounterIterator cIt = mEvent.CountersBegin();
       cIt != mEvent.CountersEnd(); ++cIt) {

    if (cIt->HasSimData()) {

      CounterSimData& cSimData = cIt->GetSimData();

      const int counterId = cIt->GetId();
      const unsigned long numParticles = cSimData.GetNumberOfGrdParticles();

      G4cerr << numParticles << " particles in Counter Id :" << counterId << G4endl;

      if (numParticles)
        G4cerr << "Id "<< counterId << "->inj. part " << numParticles << " : " << flush;
      else
        continue;

      fCurrentEventCounterIt = cIt;
      fCurrentDetectorCounter = &mDetector.GetCounter(counterId);

      for (CounterSimData::GrdParticleIterator pIt = cSimData.GrdParticlesBegin();
           pIt != cSimData.GrdParticlesEnd(); ++pIt) {

        fCurrentParticleIt = pIt;
        //utl::Particle * parent = new Particle(*fCurrentParticleIt);
        smartPointer parent( new Particle(*fCurrentParticleIt) );
        fSteppingAction->SetImpingingParticle(  parent  ); 

        //Run G4 simulation
        fRunManager->BeamOn(1);
      }//end loop on particles

      const CoordinateSystemPtr CounterCS = fCurrentDetectorCounter->GetLocalCoordinateSystem();
      //Remove particles on-ground (NO MORE NECESSARY WE HAVE TWO LISTS NOW)
      //cSimData.ClearGrdParticleList();
      
      //Loop on particles that reached the MD underground level
      //and store them into the CounterSimData class
      cout << "Counter " << counterId << " has " << fSteppingAction->fPartColl.size()  << " particles on top\n";
      SteppingAction::PartCollectionIt pCollIt;
      for(pCollIt = fSteppingAction->fPartColl.begin(); pCollIt != fSteppingAction->fPartColl.end(); ++pCollIt)
      {
        PartData pData = (*pCollIt);
        /*
        cout << "pData.Pname "   << pData.Pname   << " - " << pData.KinE/CLHEP::MeV << " (MeV)" << endl;
        cout << "pData.PDGCode " << pData.PDGCode << " - " << pData.Time/CLHEP::ns  << " (ns)" << endl;
        cout << "pData.PParent->GetType() " << pData.PParent->GetType() << " - " << (pData.PParent->GetType()==utl::Particle::eMuon||pData.PParent->GetType()==utl::Particle::eAntiMuon?"SHOWER MUON!!!!":"OTHER PARTICLE")<< endl;
        cout << "pData.TrackID " << pData.TrackID << endl;
      
        cout << "pData.x " << pData.x/CLHEP::m << endl;
        cout << "pData.y " << pData.y/CLHEP::m << endl;
        cout << "pData.z " << pData.z/CLHEP::m << endl;

        cout << "pData.px " << pData.px << endl;
        cout << "pData.py " << pData.py << endl;
        cout << "pData.pz " << pData.pz << endl;
        */
#ifdef AMIGA_DEBUG
        fStatUFile << counterId              << " "
                   << pData.Pname            << " " 
                   << pData.KinE/CLHEP::MeV  << " "
                   << pData.PDGCode          << " " 
                   << pData.Time/CLHEP::ns   << " "
                   << pData.TrackID          << " "
                   << pData.x/CLHEP::m       << " " 
                   << pData.y/CLHEP::m       << " " 
                   << pData.z/CLHEP::m       << " "
                   << pData.PParent->GetType()  << " "
                   << endl; 
#endif

        Point  pPos( pData.x/CLHEP::m * utl::m, pData.y/CLHEP::m * utl::m, pData.z/CLHEP::m * utl::m, CounterCS);
        Vector pDir( pData.px, pData.py, pData.pz, CounterCS  );
         //Creates new particle
        Particle part( pData.PDGCode, utl::Particle::eShower, pPos, pDir, pData.Time/CLHEP::ns*utl::ns, 1, pData.KinE/CLHEP::MeV * utl::MeV );
        part.SetParent( *pData.PParent ); 
        cSimData.AddUGrdParticle( part );
      }
     
      fSteppingAction->fPartColl.clear(); 
    }//end if cIt->HasSimData()
  }//end loop on counters with particles on ground

  return eSuccess;
}

VModule::ResultFlag
GroundPropagator::Finish()
{
#ifdef AMIGA_DEBUG
  fStatFile.close();
  fStatUFile.close();
#endif
  if (fUseGlobalPhysicsList) {
    tls::Geant4Manager::GetInstance().NotifyDelete();
  } else {
    if (fRunManager) {
      delete fRunManager;
      fRunManager = 0;
    }
  }

  return eSuccess;
}