MdUnderGrdInjector.cc

Go to the documentation of this file.

#include "MdUnderGrdInjector.h"

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

#include <det/Detector.h>
#include <evt/ShowerSimData.h>

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

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationSimData.h>
#include <sdet/SDetector.h>

#include <utl/ConfigUtils.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/GeometryUtilities.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>

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

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


const double UnderGrdInjector::kTriggerDelay = 0.*utl::ns;


UnderGrdInjector::UnderGrdInjector() :
  fOnlyMuons(true),
  fUseT1FromSD(true),
{ }


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

  LoadConfig(topB, "OnlyMuons", fOnlyMuons, false);
  LoadConfig(topB, "UseT1FromSD", fUseT1FromSD, true);

  return eSuccess;
}


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

  if (!event.HasMEvent()) {
    ERROR("MEvent does not exist.");
    return eFailure;
  }
  // Get the MEvent
  MEvent& mEvent = event.GetMEvent();
  INFO("\n++++\nMEvent retrieved\n++++");

  utl::TimeStamp eventTime = event.GetMEvent().GetHeader().GetTime();
  double T1 = 0;

  //Retrieve muon detector
  const mdet::MDetector& mDetector = det::Detector::GetInstance().GetMDetector();

  if (fUseT1FromSD && (!event.HasSEvent())){
    ERROR("SEvent does not exist. Do not use UseT1FromSD");
    return eFailure;
  }

  ScintiBox sbox;

#ifdef AMIGA_DEBUG
  ofstream ofile ("stats_injected.ascii");
  ofile  << "CounterId  ModuleId  ScintId  Particle  X(m)  Y(m) Z(m) Time(ns)  intersected  Parent T1(ns) EventGPS(s) EventGPSNano(ns)\n";
#endif

  //fTrgTimes.clear();

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

    G4cerr << G4endl << "Recovering underground particles from counter " << cIt->GetId() << G4endl;

    if (!cIt->HasSimData())
      continue;

    const int cId = cIt->GetId();
    const mdet::Counter& cmDet = mDetector.GetCounter(cId);
    CounterSimData& cSimData = cIt->GetSimData();
    const unsigned long numParticles = cSimData.GetNumberOfUGrdParticles();
    G4cerr << numParticles << " particles at depth level of Counter Id :" << cId << G4endl;

    //T1 is considered as the time of the first particle arriving to the Counter. To begin with, set it at "infinitum".
    T1 = std::numeric_limits<double>::max();

    //Make all modules in the MEvent structure
    for (mdet::Counter::ModuleConstIterator mmDet = cmDet.ModulesBegin();
      mmDet != cmDet.ModulesEnd(); ++mmDet) {
      int modId = mmDet->GetId();
      if (!cIt->HasModule(modId)) {
        cIt->MakeModule(modId);
        cIt->GetModule(modId).SetCandidate();
      }

    }

    //check if sd triggered or not
    if (fUseT1FromSD) {
      int triggerFound = GetTriggerTimeFromSD(event, cmDet, T1);
      if (triggerFound) {
        cout << "SD trigger found for counter " << cId << " at T1 " << T1/utl::ns << "ns " <<endl;
        cout << "Trigger is T" << triggerFound << "\n";
        cIt->SetT1(T1);
      } else {
        cout << "XXX NO SD trigger found for counter " << cId << endl;
        cIt->SetSilent(); 
        continue;  //In real life, if no T1 no event in MD
      }
    }

    for (CounterSimData::UGrdParticleIterator pIt = cSimData.UGrdParticlesBegin();
         pIt != cSimData.UGrdParticlesEnd(); ++pIt) {

      CoordinateSystemPtr cCS = cmDet.GetLocalCoordinateSystem();
      utl::Point ppos = pIt->GetPosition();
      utl::Vector pdir = pIt->GetDirection();
      //utl::TimeInterval interval = pIt->GetTime();

      utl::Point ccenter(0,0,0, cCS);
      bool injected = false;

      for (mdet::Counter::ModuleConstIterator mmDet = cmDet.ModulesBegin();
           mmDet != cmDet.ModulesEnd(); ++mmDet) {
        int MoIdx = mmDet->GetId();
        //Make module in any case
        if (!cIt->HasModule(MoIdx))
          cIt->MakeModule(MoIdx);

        if (!numParticles)
          continue;

        //Skip non muon particles if desired
        if (fOnlyMuons && pIt->GetType() != Particle::eMuon && pIt->GetType() != Particle::eAntiMuon)
          continue;

        if (IsNeutrino(pIt))
          continue;

        if (injected)
          continue;

        //The module system of reference has its x-axis along the largest side
        CoordinateSystemPtr mCS = mmDet->GetLocalCoordinateSystem();

        double l = 0;
        double w = 0;
        double h = 0;
        mmDet->ScintillatorsBegin()->VisitShape(sbox).GetDimensions(l, w, h);

        //If particle is not over the module, then continue (WARNING: this condition depends on MModelConfig.xml)
        //NOTE: it could be that a particle is over the module BUT not over the strips because the fiber routing space in the middle
        if (fabs(ppos.GetY(mCS)) > 16.*w || fabs(ppos.GetX(mCS)) > l + 0.4*utl::m)
          continue;

        mevt::Module* mmEvt = & cIt->GetModule(MoIdx);
        /*
        cout  << "particle " << pIt->GetType() << " over "
              << "Counter " << cId << " Module " << MoIdx << " "
              << ppos.GetX(cCS)/utl::m << " m "
              << ppos.GetY(cCS)/utl::m << " m "
              << ppos.GetZ(cCS)/utl::m << " m "
              << endl;
        */
        typedef mdet::Module::ScintillatorConstIterator ScIt;

        for (ScIt smDet = mmDet->ScintillatorsBegin(), esc = mmDet->ScintillatorsEnd(); smDet != esc; ++smDet) {

          const int ScIdx = smDet->GetId();

          const CoordinateSystemPtr sCS = smDet->GetLocalCoordinateSystem();

          //If particle is over the scintillator (WARNING: this condition depends on MModelConfig.xml)
          bool intersected = fabs(ppos.GetX(sCS)) < 0.5*w && fabs(ppos.GetY(sCS)) < 0.5*l;

          if (intersected) {

            if (!mmEvt->HasScintillator(ScIdx))
              mmEvt->MakeScintillator(ScIdx);

            mevt::Scintillator* const smEvt = &mmEvt->GetScintillator(ScIdx);

            if (! smEvt->HasSimData())
              smEvt->MakeSimData();

            mevt::ScintillatorSimData* smSimData = &smEvt->GetSimData();

            Particle injectedParticle(*pIt);
            //Set new z position (impact point on the scintillator)
            const Point injectedPosition(ppos.GetX(sCS), ppos.GetY(sCS), 0.5*h, sCS);
            injectedParticle.SetPosition(injectedPosition);

            //Load particle in the proper event scintillator
            smSimData->AddParticle(injectedParticle);
            if (injectedParticle.GetType() == Particle::eMuon || injectedParticle.GetType() == Particle::eAntiMuon)
              smSimData->IncrementNumberOfInjectedMuons();

            utl::Line path(injectedPosition, pdir);

#ifdef AMIGA_DEBUG
            ofile
              << cId << " " << MoIdx << " " << ScIdx << " "
              << pIt->GetType() << " "
              << injectedPosition.GetX(cCS)/utl::m << " "
              << injectedPosition.GetY(cCS)/utl::m << " "
              << injectedPosition.GetZ(cCS)/utl::m << " "
              << interval.GetInterval()/utl::ns << " "
              << smDet->IntersectedBy(path) << " "
              << pIt->GetParent().GetType() << " "
              << cIt->GetT1()/utl::ns << " "
              << eventTime.GetGPSSecond()/utl::s << " "
              << eventTime.GetGPSNanoSecond()/utl::ns << " "
              << endl;
#endif
            injected = true;
            break;
          } // end if injected
        } // end scintillator det loop
      } // end modules det loop
    } // end particle loop over a given counter
  } // end counters loop

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

    const int cId = cIt->GetId();

    if (!cIt->HasSimData())
      continue;

    if (cIt->IsSilent())
      continue;

    T1 = cIt->GetT1();
    //utl::TimeInterval trgTime = fTrgTimes[cId];

    const mdet::Counter& cmDet = mDetector.GetCounter(cId);
    utl::TimeStamp coreTime;

    if (event.HasSimShower()) {
      const ShowerSimData& sSimShower = event.GetSimShower();
      coreTime = sSimShower.GetTimeStamp();
    }

    for (mevt::Counter::ModuleIterator mIt = cIt->ModulesBegin(), em = cIt->ModulesEnd();
         mIt != em; ++mIt) {

      const int mId = mIt->GetId();
      const mdet::Module & mmDet = cmDet.GetModule(mId);

      if(!mmDet.IsSiPM()){ //If PMT

                  //Prepares to iterate over module scintillators (mevt) and module channels (mdet)
                  //typedef mevt::Module::ScintillatorIterator ScIt;
                  typedef mdet::FrontEnd::ChannelConstIterator chIt;
                  const mdet::FrontEnd& frontEnd = mmDet.GetFrontEnd();
                  const double sRate = frontEnd.GetMeanSampleRatePeriod();
                  const double preT1BufferLength = frontEnd.GetPreT1BufferLength();
                  const double postT1BufferLength = frontEnd.GetPostT1BufferLength();

                  //shift the position of the T1 to the circular md-buffer (preT1BufferLength)
                  utl::TimeInterval timeShiftpre(sRate*preT1BufferLength);
                  utl::TimeInterval timeShiftpost(sRate*postT1BufferLength);

                  //cout << "timeShiftpre " << timeShiftpre << endl;
        /*
                  for ( ScIt sc = mIt->ScintillatorsBegin(), esc = mIt->ScintillatorsEnd();  sc != esc; ++sc) {

                        if (! sc->HasSimData() )
                          continue;//only scintillators with injected particles has SimData

                        mevt::ScintillatorSimData* ssd = & sc->GetSimData();

                  const mdet::Scintillator& scint = mmDet.GetScintillator(sc->GetId());

                        for (mevt::ScintillatorSimData::ParticleIterator pi = ssd->ParticlesBegin(), sse = ssd->ParticlesEnd(); pi != sse; ++pi)
                        {
                          utl::Particle& particle = *pi;
                          utl::TimeInterval interval = particle.GetTime();

                          //firstTime = std::min(interval.GetInterval(), T1);

                          cout
                                  << "Ctr " << cId << " (at " << counterR_i/utl::m << "m) Mod " << mId << " Sci " << sc->GetId() << " "
                                  << "PTYPE "<< pi->GetType()<< " "
                                  << "DELAY " << interval/utl::ns << " ns, "
                                  //<< "START TRACE "   << utl::TimeInterval(T1+trgTime.GetSecond()-timeShiftpre)/utl::ns << " ns "
                                  << "START TRACE " << (utl::TimeInterval(T1) + trgTime - utl::TimeInterval(timeShiftpre)) /utl::ns << " ns "
                                  << "INJECTED AT " << (utl::TimeInterval(T1)+interval) /utl::ns << " ns "
                                  << "T1 "            << utl::TimeInterval(T1)         << " ns "
                                  //<< "fTrgTimes[" << cId << "]=" << fTrgTimes[cId] << " " << " "
                                  << "trgTime "       << trgTime/utl::ns << " ns "
                                  //<< "timeShiftpre  " << timeShiftpre/utl::ns
                                  << endl;

                          //particles havewith respect to the counter T1
                          particle.SetTime( utl::TimeInterval(T1)+interval );

                         }//end particles injected in scintillator loop
                  }//end scintillator loop
        */

                   for (chIt ch = frontEnd.ChannelsBegin(), ec = frontEnd.ChannelsEnd(); ch != ec; ++ch) {

                        if (!mIt->HasChannel(ch->GetId()))
                          mIt->MakeChannel(ch->GetId());

                        // Compute the trace start wrt the station trigger T1 GPS time
                        const utl::TimeStamp traceStart = eventTime + utl::TimeInterval(T1) - timeShiftpre;
                        mIt->GetChannel(ch->GetId()).SetTraceStartTime(traceStart);
        /*
                        cout << "Setting trace start at " << mIt->GetChannel(ch->GetId()).GetTraceStartTime()    << " "
                                 << "T1 " << mIt->GetChannel(ch->GetId()).GetTraceStartTime()-eventTime+timeShiftpre << " "
                                 << "(before: " << T1 << ")"
                                 << endl;
        */

                  } // end for channel loop
      }else{//If SiPM

                  //Prepares to iterate over module channels (mdet)

                  typedef mdet::FrontEndSiPM::ChannelConstIterator chIt;
                  const mdet::FrontEndSiPM& frontEnd = mmDet.GetFrontEndSiPM();
                  const double sRate = frontEnd.GetMeanSampleRatePeriod();
                  const double preT1BufferLength = frontEnd.GetPreT1BufferLength();
                  const double postT1BufferLength = frontEnd.GetPostT1BufferLength();

                  //shift the position of the T1 to the circular md-buffer (preT1BufferLength)
                  utl::TimeInterval timeShiftpre(sRate*preT1BufferLength);
                  utl::TimeInterval timeShiftpost(sRate*postT1BufferLength);

                   for (chIt ch = frontEnd.ChannelsBegin(), ec = frontEnd.ChannelsEnd(); ch != ec; ++ch) {

                        if (!mIt->HasChannel(ch->GetId()))
                          mIt->MakeChannel(ch->GetId());

                        // Compute the trace start wrt the station trigger T1 GPS time
                        const utl::TimeStamp traceStart = eventTime + utl::TimeInterval(T1) - timeShiftpre;
                        mIt->GetChannel(ch->GetId()).SetTraceStartTime(traceStart);

                  } // end for channel loop
      }
    } // end modules loop
  } // end counter loop

  //fTrgTimes.clear();
#ifdef AMIGA_DEBUG
  ofile.close();
#endif

  return eSuccess;
}


bool
UnderGrdInjector::IsNeutrino(mevt::ScintillatorSimData::ParticleIterator pIt)
{
  return
    pIt->GetType() == Particle::eNuElectron     ||
    pIt->GetType() == Particle::eAntiNuElectron ||
    pIt->GetType() == Particle::eNuMuon         ||
    pIt->GetType() == Particle::eAntiNuMuon     ||
    pIt->GetType() == Particle::eNuTau          ||
    pIt->GetType() == Particle::eAntiNuTau;
}


bool
UnderGrdInjector::IsChargedParticle(mevt::ScintillatorSimData::ParticleIterator pIt)
{
  return
    pIt->GetType() == Particle::ePhoton         ||
    /*neutrinos are skipped at the particleList level
    pIt->GetType() == Particle::eNuElectron     ||
    pIt->GetType() == Particle::eAntiNuElectron ||
    pIt->GetType() == Particle::eNuMuon         ||
    pIt->GetType() == Particle::eAntiNuMuon     ||
    pIt->GetType() == Particle::eNuTau          ||
    pIt->GetType() == Particle::eAntiNuTau      ||*/
    pIt->GetType() == Particle::ePiZero         ||
    pIt->GetType() == Particle::eKaon0L         ||
    pIt->GetType() == Particle::eKaon0S         ||
    pIt->GetType() == Particle::eNeutron        ||
    pIt->GetType() == Particle::eAntiNeutron;
}


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


int
UnderGrdInjector::GetTriggerTimeFromSD(evt::Event& event, const mdet::Counter& cmDet, double& T1)
{
  if (!event.HasSEvent()) {
    ERROR("SEvent does not exist.");
    return 0;
  }
  sevt::SEvent& sEvent = event.GetSEvent();

  int partnerId = cmDet.GetAssociatedTankId();

  if (!sEvent.HasStation(partnerId)) {
    ostringstream message;
    message << "partner station with id " << partnerId << "was not found in SEvent\n";
    WARNING(message);
    return 0;
  }

  sevt::Station& sStation = sEvent.GetStation(partnerId);
  const sdet::Station& detStation = Detector::GetInstance().GetSDetector().GetStation(sStation);

  if (!sStation.HasSimData())
    return 0;

  int ret = 0;
  sevt::StationSimData& sSimData = sStation.GetSimData();
  utl::TimeStamp firstT1Sd;

  for (sevt::StationSimData::TriggerTimeIterator it = sSimData.TriggerTimesBegin(); it != sSimData.TriggerTimesEnd(); ++it) {

    if (it == sSimData.TriggerTimesBegin())
      firstT1Sd = *it;

    utl::TimeStamp trigTime = *it;
    const sevt::StationTriggerData& trig = sSimData.GetTriggerData(trigTime);
    cout << "IsT1 " << trig.IsT1() << " - " << "IsT2 " << trig.IsT2() << endl;

    if (!(trig.IsT2() || trig.IsT1()))
      continue;

    ret = 1;
    if (trig.IsT2())
       ret = 2;

    if (trigTime < firstT1Sd) {
      cout << "There is a previous trigger\n";
      firstT1Sd = trigTime;
    }
    cout << "Trigger time: " << trigTime << endl;
    cout << "Algorithm: "     << trig.GetAlgorithmName() << endl;

  } // end loop on TriggerTimeIterator

  if (!ret)
    return 0;

  utl::TimeStamp eventTime = event.GetMEvent().GetHeader().GetTime();
  //see sdet::StationTriggerAlgorithm.h and TankTriggerSimulator.cc
  utl::TimeInterval shift = utl::TimeInterval((detStation.GetFADCTraceLength() - detStation.GetLatchBin())*detStation.GetFADCBinSize());
  utl::TimeStamp TstartTrace = firstT1Sd - shift;
  //Time interval from Eventime to TstartTrace
  T1 = TstartTrace - eventTime;

  cout << "eventTime "   << eventTime   << endl;
  cout << "firstT1Sd "   << firstT1Sd   << endl;
  cout << "TstartTrace " << TstartTrace << endl;
  cout << "shift "       << shift       << endl;
  cout << "T1 "          << T1          << endl;

/*
  //T1 is absolute GPS T1 with respect to event time
  utl::TimeInterval shiftToEvent = firstT1Sd - eventTime;

  //Keep track of how much the local T1 GPS is delayed with
  //respect to core time (event time). This values
  //is dependent of radial distance (through particle delay times)
  fTrgTimes[cmDet.GetId()] = shiftToEvent-shift;
  T1       = shiftToEvent.GetInterval();

  cout << "eventTime "    << eventTime << endl;
  cout << "firstT1Sd " << firstT1Sd << " " << endl;
  cout << "shift "        << shift << " " << endl;
  cout << "shiftToEvent " << shiftToEvent << " " << endl;
  cout << "shiftToEvent - shift " << shiftToEvent-shift << " " << endl;
  cout << "fTrgTimes[" << cmDet.GetId() << "]=" << fTrgTimes[cmDet.GetId()] << " " << endl;
*/

  return ret;
}