MdEdepSimulator.cc

Go to the documentation of this file.

#include "MdEdepSimulator.h"
#include "PhysicsListCustomization.h"
#include "PhysicsList.h"
#include "ScintillatorConstruction.h"
#include "PrimaryGenerator.h"
#include <G4RunManager.hh>
#include <G4UImanager.hh>
#include <G4VisManager.hh>
#include <G4VisExecutive.hh>

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

#include <det/Detector.h>
#include <mdet/Counter.h>
#include <mdet/Module.h>
#include <mdet/Scintillator.h>

#include <evt/Event.h>
#include <mevt/Counter.h>
#include <mevt/Module.h>
#include <mevt/Scintillator.h>
#include <mevt/CounterSimData.h>
#include <mevt/ScintillatorSimData.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 <utl/TabularStream.h>
#include <utl/ConfigUtils.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 EdepSimulatorAG;

const mdet::Scintillator* EdepSimulator::fCurrentDetectorScintillator = 0;
mevt::Module::ScintillatorIterator EdepSimulator::fCurrentEventScintillatorIt;
mevt::ScintillatorSimData::ParticleIterator EdepSimulator::fCurrentParticleIt;

EdepSimulator::EdepSimulator() :
  fUImanager(0),
  fVisManager(0),
  fSteppingAction(0)
{
  fGeomVisOn = false;
  fTrajVisOn = false;
  fUseGlobalPhysicsList = false;
}

EdepSimulator::~EdepSimulator() { }

VModule::ResultFlag
EdepSimulator::Init()
{
  //Top branch

  Branch config = CentralConfig::GetInstance()->GetTopBranch("MdEdepSimulator");
  fLog.Configure(config);
  //Set energy units to MeV (otherwise default is eV)
  fUnits.SetEnergyDefault( utl::MeV, "MeV");
  fUnits.SetLengthDefault( utl::cm , "cm");
  fUnits.Configure(config);

  LoadConfig( config.GetChild("visualization"), "geometry", fGeomVisOn, 0 );
  LoadConfig( config.GetChild("visualization"), "trajectories", fTrajVisOn, 0 );

  //Physics list branch
  string physList;

  LoadConfig( config.GetChild("physics"), "list", physList, "QGSP_BERT");
  LoadConfig(config, "globalPhysicsList", fUseGlobalPhysicsList, 1 );

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

  LoadConfig(config, "verbosityLevel", fVerbosityLevel, 0 );
  stringstream oss (ostringstream::out);
  fUImanager  = G4UImanager::GetUIpointer();

  oss << "/run/verbose " << fVerbosityLevel;
  fUImanager->ApplyCommand( oss.str() );
  oss.str("");
  oss << "/event/verbose " << fVerbosityLevel;
  fUImanager->ApplyCommand( oss.str() );
  oss.str("");
  oss << "/tracking/verbose " << fVerbosityLevel;
  fUImanager->ApplyCommand( oss.str() );
  oss.str("");
  oss << "/stepping/verbose " << fVerbosityLevel;
  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");
    ScintillatorConstruction* scintConstruction = new ScintillatorConstruction();

    tls::Geant4Customization custom("MdEdepSimulator", scintConstruction, physicsListCustomization);

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

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

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

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

    INFO("\n++++\nConstructing G4 ground");
    fRunManager->SetUserInitialization(new ScintillatorConstruction());

    INFO("\n++++\nInitializing G4 stepping action");
    fSteppingAction = new SteppingAction();
    fRunManager->SetUserAction(fSteppingAction);

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

    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;
}

void
EdepSimulator::AddUpContributions(  SteppingAction::ScintMap & scintMap )
{
  string Sname, Pname;
  int PDGCode(0), ParentID(0), TrackID(0);
  //Position
  double x(0), y(0), z(0);
  //Direction
  double px(0), py(0), pz(0);
  //Energy & Time
  double KinE(0), Time(0);

  SteppingAction::ScintMapIt sstepIt;
  DataCollection::iterator dIt;

  SteppingAction::ScintMap sstepMap = fSteppingAction->GetScintMap();
//  G4cerr << "Map size before steps integration " << sstepMap.size() << endl;
  for ( sstepIt = sstepMap.begin(); sstepIt != sstepMap.end(); ++sstepIt ) {

    DataCollection &cur = (*sstepIt).second;
//     cout << " Scintillator ID " << (*sstepIt).first << endl;

    double totalEdep = 0;
    bool first_particle = true;//to trace first particle entering in each scintillator

    for ( dIt = cur.begin(); dIt != cur.end(); ++dIt ) {

      Data data = (*dIt);
      totalEdep += data.Edep;

      if ( first_particle ){
        Sname   = data.Sname   ;
        Pname   = data.Pname   ;
        PDGCode = data.PDGCode ;
        ParentID= data.ParentID;
        TrackID = data.TrackID ;
        x       = data.x       ;
        y       = data.y       ;
        z       = data.z       ;
        px      = data.px      ;
        py      = data.py      ;
        pz      = data.pz      ;
        KinE    = data.KinE    ;
        Time    = data.Time    ;
        first_particle = false;
      }
/*
      cout << "Pname "   << data.Pname   << " - " << data.KinE/CLHEP::MeV << " (MeV)" << endl;
      cout << "PDGCode " << data.PDGCode << " - " << data.Time/CLHEP::ns  << " (ns)" << endl;
      cout << "ParentID " << data.ParentID << endl;
      cout << "TrackID " << data.TrackID << " - " << data.Sname << endl;
      cout << "Edep "    << data.Edep/CLHEP::MeV << " (MeV)" << endl;
      cout << "x " << data.x/CLHEP::m << endl;
      cout << "y " << data.y/CLHEP::m << endl;
      cout << "z " << data.z/CLHEP::m << endl;
      cout << "px " << data.px << endl;
      cout << "py " << data.py << endl;
      cout << "pz " << data.pz << endl;
*/
    }//end loop on steps

   //If scintId (key) not found, then creates the key
    if ( scintMap.find( (*sstepIt).first ) == scintMap.end()){
      DataCollection  nuovo;
//        G4cerr << "Scinti ID " << (*sstepIt).first  << " not found in scintMap... I will creat it" << G4endl;
      scintMap[ (*sstepIt).first] = nuovo;
    }
    Data partdata;
    //Fill obj with relevant step data
    partdata.Sname   = Sname;
    partdata.Pname   = Pname;
    partdata.PDGCode = PDGCode;
    partdata.ParentID= ParentID;
    partdata.TrackID = TrackID;
    partdata.x       = x;
    partdata.y       = y;
    partdata.z       = z;
    partdata.px      = px;
    partdata.py      = py;
    partdata.pz      = pz;
    partdata.KinE    = KinE;
    partdata.Edep    = totalEdep;//HERE IS ALL WE NEED ACTUALLY :-)
    partdata.Time    = Time;
    partdata.PParent = fSteppingAction->GetImpingingParticle();
    //Fill the vector with the data of a given scintillator
    scintMap[(*sstepIt).first].push_back( partdata );
  }//end loop on scintillators
  fSteppingAction->ClearMap();//Map of steps
//  G4cerr << "Cleared map (size " << fSteppingAction->GetScintMap().size() << ")" << G4endl << G4endl;
}

VModule::ResultFlag
EdepSimulator::Run(Event& event)
{
  if (!event.HasMEvent()) {
    ERROR("\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("MdEdepSimulator");
    }
  }

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

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

    const int counterId = cIt->GetId();
    const mdet::Counter& counter = mDetector.GetCounter(counterId);
    for (mevt::Counter::ModuleIterator im = cIt->ModulesBegin(), em = cIt->ModulesEnd(); im != em; ++im) {

      const int modId = im->GetId();
      const mdet::Module& module = counter.GetModule(modId);

      SteppingAction::ScintMap scintMap;

      for (mevt::Module::ScintillatorIterator is = im->ScintillatorsBegin(), es = im->ScintillatorsEnd();
            is != es; ++is) {
        if (is->HasSimData()) {
          //const mdet::Scintillator& scint = module.GetScintillator(is->GetId());
          mevt::ScintillatorSimData& sSimData = is->GetSimData();
          fCurrentEventScintillatorIt = is;
          const int scintId = is->GetId();
          fCurrentDetectorScintillator = &module.GetScintillator(scintId);

          if( sSimData.GetNumberOfParticles() != 0 ) {
            fLog("Counter ", 1, false)( counterId )(" ");
            fLog("Module ", 1, false)( modId )(" ");
            fLog("Scint. ", 1)( scintId )(" has : ")( sSimData.GetNumberOfParticles())(" particles on top" );
          }

          fSteppingAction->SetScintillatorId( scintId );//Store ID of current scintillator (to distinguish it from neighours)

          for (ScintillatorSimData::ParticleIterator pIt = sSimData.ParticlesBegin();
                pIt != sSimData.ParticlesEnd(); ++pIt) {

            fCurrentParticleIt = pIt;
            if ( fCurrentParticleIt->HasParent() ) {
              utl::Particle * parent = new Particle(fCurrentParticleIt->GetParent() );
              fSteppingAction->SetImpingingParticle(  parent  );
            }

            //Run G4 simulation
            fRunManager->BeamOn(1);

            //Integrate all particle steps plus secondary particles data
            //in central and (if needed) in lateral scintillators
            //where particles might have deposited energy
            AddUpContributions( scintMap );

          }//end loop on particles in a given scintillator

          //Remove already processes particles in this scintillator
          sSimData.ClearParticleList();

        }//end if is->HasSimData()
      }//scintillator loop

      SteppingAction::ScintMapIt scintIt;

       // Table for the inserted particles.
       const unsigned int particleTableLogLevel = 2;
       std::unique_ptr<utl::TabularStream> table;
       if (fLog.GetLevel() >= particleTableLogLevel) {
         //-----------------------------------1-2-3-4-5-6-7-8-9-10-11-12
         table.reset(new utl::TabularStream("|r|r|r|r|r|r|r|r|r|r|r|r|"));
         fLog.ApplyConfigurationOn(*table);
         *table << utl::hline
                << '#'           << utl::endc                                         //1
                << "name"        << utl::endc                                         //2
                << "parentID"    << utl::endc                                         //3
                << "trackID"     << utl::endc                                         //4
                << "kin. energy" << "(" << fUnits.GetEnergyName() << ")" << utl::endc //5
                << "x"           << "(" << fUnits.GetLengthName() << ")" << utl::endc //6
                << "y"           << "(" << fUnits.GetLengthName() << ")" << utl::endc //7
                << "z"           << "(" << fUnits.GetLengthName() << ")" << utl::endc //8
                << "energy dep." << "(" << fUnits.GetEnergyName() << ")" << utl::endc //9
                << "px  "        << utl::endc                                         //10
                << "py  "        << utl::endc                                         //11
                << "pz  "        << utl::endr                                         //12
               << utl::hline;
       }

      for ( scintIt = scintMap.begin(); scintIt != scintMap.end(); ++scintIt ) {

        DataCollection::iterator dIt;
        DataCollection &cur = (*scintIt).second;
        int ScIdx = (*scintIt).first;
        /*
        cout << "\n*************\n";
        cout << "Found Scintillator ID " << ScIdx << endl;
        */
        if (!im->HasScintillator(ScIdx)) {
          im->MakeScintillator(ScIdx);
          im->GetScintillator(ScIdx).MakeSimData();
        }

        for (mevt::Module::ScintillatorIterator is = im->ScintillatorsBegin(), es = im->ScintillatorsEnd();
              is != es; ++is) {

          const int scintId = is->GetId();

          if( (*scintIt).first != scintId )
            continue;

          const CoordinateSystemPtr ScintillatorCS = module.GetScintillator(scintId).GetLocalCoordinateSystem();
          mevt::ScintillatorSimData& sSimData = is->GetSimData();

          fLog("Particles that deposited energy in scintillator ", 5 )(scintId)(" : ")( cur.size() );

          int n = 0;
          for ( dIt = cur.begin(); dIt != cur.end(); ++dIt, ++n ) {

            Data data = (*dIt);

            // only accumulate in the table if it will be used at last:
            if (fLog.GetLevel() >= particleTableLogLevel) {
              *table << (n+1)                                                            << utl::endc
                     << fLog << data.Pname << " (PDG code " << data.PDGCode << ")"       << utl::endc
                     << fLog << data.ParentID                                            << utl::endc
                     << fLog << data.TrackID                                             << utl::endc
                     << fLog << data.KinE/CLHEP::MeV * utl::MeV / fUnits.GetEnergyUnit() << utl::endc
                     << fLog << data.x   /CLHEP::m   * utl::m   / fUnits.GetLengthUnit() << utl::endc
                     << fLog << data.y   /CLHEP::m   * utl::m   / fUnits.GetLengthUnit() << utl::endc
                     << fLog << data.z   /CLHEP::m   * utl::m   / fUnits.GetLengthUnit() << utl::endc
                     << fLog << data.Edep/CLHEP::MeV * utl::MeV / fUnits.GetEnergyUnit() << utl::endc
                     << fLog << data.px                                                  << utl::endc
                     << fLog << data.py                                                  << utl::endc
                     << fLog << data.pz                                                  << utl::endr;
            }

            Point  pPos( data.x/CLHEP::m * utl::m, data.y/CLHEP::m * utl::m, data.z/CLHEP::m * utl::m, ScintillatorCS);
            Vector pDir( data.px, data.py, data.pz, ScintillatorCS );
            //Creates new particle (uses kinetic energy to pass the energy deposited by the particle)
            //#warning "ENERGY DEPOSITED BY THE PARTICLE IS NOW PASSED in the WEIGHT field (a la MARTA :-))"
            Particle part( data.PDGCode, utl::Particle::eShower, pPos, pDir, data.Time/CLHEP::ns*utl::ns, data.Edep/CLHEP::MeV*utl::MeV, data.KinE/CLHEP::MeV * utl::MeV );


            if ( data.PParent )
        part.SetParent( *data.PParent );

            sSimData.AddParticle( part );

          }//end loop on steps

        }//end of mevt::scintillator loop
      }//end of scintiMap loop
      scintMap.clear();
      scintMap.erase( scintMap.begin(), scintMap.end() );

      if (fLog.GetLevel() >= particleTableLogLevel) {
        *table << utl::hline;
        fLog("Particle table after G4 Scintillator simulation:", particleTableLogLevel);
        fLog(table->Str(), particleTableLogLevel);
      }
    }//end of module loop
  }//end loop on counters with particles on ground
  cout << endl;

  return eSuccess;
}

VModule::ResultFlag
EdepSimulator::Finish()
{
  if (fUseGlobalPhysicsList) {
    tls::Geant4Manager::GetInstance().NotifyDelete();
  } else {
    if (fRunManager) {
      delete fRunManager;
      fRunManager = 0;
    }
  }

  return eSuccess;
}