In This Package:

MixInputSvc.cc

Go to the documentation of this file.

#include "MixInputSvc/MixInputSvc.h"

#include "Event/GenHeader.h"
#include "Event/SimHeader.h"
#include "Event/RegistrationSequence.h"

#include "RootIOSvc/IRootIOSvc.h"
#include "RootIOSvc/RootIOBaseCnv.h"
#include "RootIOSvc/RootIOAddress.h"
#include "RootIOSvc/RootIOTypedCnv.h"

#include "PerBaseEvent/PerRef.h"
#include "PerSimEvent/PerSimHit.h"
#include "PerGenEvent/PerGenHeader.h"
#include "PerSimEvent/PerSimHeader.h"
#include "PerSimEvent/PerSimHitHeader.h"
#include "PerSimEvent/PerSimHitCollection.h"

#include "GaudiKernel/IService.h"
#include "GaudiKernel/IRndmGen.h"
#include "GaudiKernel/IRndmEngine.h"
#include "GaudiKernel/IRndmGenSvc.h"
#include "GaudiKernel/IConversionSvc.h"
#include "GaudiKernel/IDataManagerSvc.h"

MixInputSvc::MixInputSvc(const string& name, ISvcLocator* svc)
            :RootIOCnvSvc(name, svc),/* m_ranSvc(0),*/ m_log(msgSvc(), "MixInputSvc")
{
  m_isvc    = 0;
  m_rioSvc  = 0;
  m_convSvc = 0;

  declareProperty("SampleInputs" , m_inputSamples  , "Map of sample input files");
  declareProperty("SampleTaus"   , m_timeTau       , "Map of sample Tau");
  declareProperty("SampleFileNum", m_fileEntryNum  , "Map of sample input files' numbers");
  declareProperty("BeginTime"    , m_lasttime = 0.0, "Begin time");
}

MixInputSvc::~MixInputSvc()
{
  m_log<< MSG::DEBUG<< "~MixInputSvc()"<< endreq;
}

StatusCode MixInputSvc::initialize() 
{
  m_log<< MSG::DEBUG<< "initializing @"<< (void*)this<< endreq;
  StatusCode sc = RootIOCnvSvc::initialize();
  if(sc.isFailure()) 
  {
    m_log<< MSG::ERROR<< "Parent class(RootIOCnvSvc) failed to initialize!"<< endreq;
    return sc;
  } 

  /*sc = service("RndmGenSvc", m_ranSvc, true);
  if(sc.isFailure())
  {
    m_log<< MSG::ERROR<< "Unable get RndmGenSvc service"<< m_ranSvc<< endreq;
    return sc;
    }*/

  sc = doInit();
  if(sc.isFailure())
  {
    m_log<< MSG::ERROR<< "Something wrong before mixing."<<endreq;
    return sc;
  }
  return sc;
}

StatusCode MixInputSvc::doInit()
{
  m_log<< MSG::DEBUG<< "Reset some variables at begining."<< endreq;
  if(m_inputSamples.size() != m_timeTau.size() || m_timeTau.size() == 0) 
  {
    m_log<< MSG::ERROR<< "Please check the input samples' files "
         << "and the input samples' taus!"<< endreq;
    return StatusCode::FAILURE;
  }
  else
  {
    map<string, double>::iterator itt, tend = m_timeTau.end();
    map<string, vector<string> >::iterator iter; 
    for(itt = m_timeTau.begin(); itt != tend; itt++) 
    {
      iter = m_inputSamples.find(itt->first);
      if(iter == m_inputSamples.end()) 
      {
        m_log<< MSG::ERROR<< "Configuration is not correct for mixing"<< endreq;
        return StatusCode::FAILURE;
      }
      m_fileNum[itt->first]       = 0;
      m_fileindex[itt->first]     = 0;
      m_gfileentries[itt->first]  = 0;
      m_sfileentries[itt->first]  = 0;
      m_streamentries[itt->first] = 0;
    }
    m_maintau  = 0.0;
    m_maximum  = 0.0;
    TimeStamp  lastTime(m_lasttime);
    m_lastTime = lastTime;
  }
  return StatusCode::SUCCESS;
}

StatusCode MixInputSvc::queryInterface(const InterfaceID& riid, void** ppint)
{
  if(IID_IMixInputSvc.versionMatch(riid)) 
  {
    m_log<< MSG::DEBUG<< "queryInterface("<< riid
         << ")-->(IMixInputSvc*)"<< (void*)this<< endreq;
    *ppint = (IMixInputSvc*)this;
  }
  else 
  {
    return RootIOCnvSvc::queryInterface(riid, ppint);
  }
  addRef();
  return StatusCode::SUCCESS;
}

StatusCode MixInputSvc::prepareStream()
{
  m_log<< MSG::INFO<< "PrepareStream............"<< endreq;
  //Need to prepare or not
  if(0 == m_sampleInputStreamMap.size()) 
  {
    m_log<< MSG::DEBUG<< "Need to prepare streams from mixing sources."<< endreq;
    map<string, vector<string> >::iterator miter, mbegin, mend;
    mbegin = m_inputSamples.begin();
    mend   = m_inputSamples.end();
    m_log<< MSG::INFO<< "Initialize the mixing samples here."<< endreq;
    
    //Initialize the files of each sample
    for(miter = mbegin; miter != mend; miter++) 
    {
      m_inputStreamMap.clear();

      vector<string> files = miter->second;
      setInputFiles(files);
      doInitialize();

      m_sampleInputStreamMap[miter->first] = inputStreams();
      m_simMap[miter->first] = m_inputStreamMap[DayaBay::SimHeaderLocation::Default];
      m_genMap[miter->first] = m_inputStreamMap[DayaBay::GenHeaderLocation::Default];
    }

    //Initialize the file numbers of each sample
    for(miter = mbegin; miter != mend; miter++)
    {
      int filenum = 0;
      vector<string> files = miter->second;
      for(vector<string>::iterator vit = files.begin(); vit != files.end(); vit++)
      {
        filenum++;
      }
      m_fileNum[miter->first] = filenum;
    }
    m_log<< MSG::INFO<< "Initialize the mixing samples successfully."<< endreq;

    //Initialize the file index, file entry and stream entry for each sample
    m_log<< MSG::INFO<< "Initialize the mixing samples' file index, "
         << "file entry and stream entry here."<< endreq;    

    for(map<string, RootInputStream*>::iterator it = m_simMap.begin(); it != m_simMap.end(); it++)
    {
      RootInputStream* ris      = it->second;
      //m_fileindex[it->first]     = flatdis((double)ris->getFiles().size());
      m_fileindex[it->first]    = flatdis((double)m_fileNum[it->first]);
      m_sfileentries[it->first] = flatdis((double)m_fileEntryNum[it->first]);
      /*bool ifjump                = ris->getFiles().jump(m_fileindex[it->first]);
      if(ifjump == false)
      {
        m_log<< MSG::ERROR<< "Can't jump to the specified file!"<< endreq;
        return StatusCode::FAILURE;
        }*/
      m_log<< MSG::INFO<< "Now needing to jump to "<< it->first<< ": "
           << m_inputSamples[it->first][m_fileindex[it->first]]<< " file"<< endreq;
      m_log<< MSG::INFO<< "Now the file entry is: "<< m_sfileentries[it->first]<< endreq;

      ris->setFileEntry(m_fileindex[it->first], m_sfileentries[it->first], false);
      m_streamentries[it->first] = ris->getEntry();
    }

    m_log<< MSG::INFO<< "SteamMap Size = "<< m_sampleInputStreamMap.size()
         << "/"<< m_simMap.size()<< endreq;
    
    //Initialize the maintau, the time nod by the samples
    m_log<< MSG::INFO<< "Initialize the mixing samples' maintau, time nods here."<< endreq;    
    for(map<string, RootInputStream*>::iterator mit = m_simMap.begin(); mit != m_simMap.end(); mit++)
    { 
      m_maintau += 1/m_timeTau[mit->first];
      m_nod[mit->first] = m_maintau;
      m_log<< MSG::INFO<< "m_nod["<< mit->first<< "] is: "<< m_nod[mit->first]<< endreq;
      m_maximum  = m_maintau;
    }
    m_maintau = 1/m_maintau;
    m_log<< MSG::INFO<< "The main tau is: "<< m_maintau<< endreq;
  }
  else
  {
    m_log<< MSG::INFO<< "Streams prepared! Nothing to do!"<< endreq;
  }

  //if (m_uniSample.initialize(m_ranSvc, Rndm::Flat(0.0, m_maximum)).isFailure()) 
  //{
  //  m_log<< MSG::FATAL<< "Failed to initialize uniform random numbers" << endreq;
  //  return StatusCode::FAILURE;
  //}
  //if (m_exp.initialize(m_ranSvc, Rndm::Exponential(m_maintau)).isFailure()) 
  //{
  //  m_log<< MSG::FATAL<< "Failed to initialize Gaussian random numbers" << endreq;
  //  return StatusCode::FAILURE;
  //}
  return StatusCode::SUCCESS;
}

StatusCode MixInputSvc::getSimHeader(GenHeader*& gh, SimHeader*& sh)
{
  m_log<< MSG::INFO<< "Now get one SimHeader(GenHeader) from the samples!"<< endreq;
  
  //Initialize Service
  StatusCode sc = this->service("RootIOCnvSvc", m_isvc, true);
  if(sc.isFailure())
  {
    m_log<< MSG::ERROR<< "Failed to get RootIOCnvSvc as IService"<< endreq;
    return sc;
  }

  //Initialize RootIOSvc
  sc = m_isvc->queryInterface(IRootIOSvc::interfaceID(), (void**)&m_rioSvc);
  if(sc.isFailure()) 
  {
    m_log<< MSG::ERROR<< "Conversion service RootIOCnvSvc"
         << " does not implement IRootIOCnvSvc"<< endreq;
    return sc;
  }

  //Initialize ConversionSvc
  sc = m_isvc->queryInterface(IConversionSvc::interfaceID(), (void**)&m_convSvc);
  if(sc.isFailure())
  {
    m_log<< MSG::ERROR<< "Conversion service RootIOCnvSvc"
         << " does not implement IConversionSvc"<< endreq;
    return sc;
  }   

  double randomnum = flatdis(m_maximum);
  //double randomnum = m_uniSample();
  map<string, RootInputStream*>::iterator it;

  //Loop on the samples
  for(it = m_simMap.begin(); it != m_simMap.end(); ) 
  {
    //using randomnum to decide the input sample
    if(randomnum < m_nod[it->first])
    {
      //SetEntry to find that if reaching the end of a file
      RootInputStream* ris = it->second;
      RootInputFile*   rif = ris->getFiles().current();

      RootInputStream* gris = m_genMap[it->first];
      RootInputFile*   grif = gris->getFiles().current();
      
      if(m_sfileentries[it->first] >= rif->entries())
      {
        m_sfileentries[it->first] = flatdis((double)m_fileEntryNum[it->first]);
        m_fileindex[it->first]    = flatdis((double)m_fileNum[it->first]);
        /*bool ifjump              = ris->getFiles().jump(m_fileindex[it->first]);
        if(ifjump == false)
        {
          m_log<< MSG::ERROR<< "Can't jump to specified file!"<< endreq;
          return StatusCode::FAILURE;
          }*/
        m_log<< MSG::INFO<< "Now needing to jump to "<< it->first<< ": "
             << m_inputSamples[it->first][m_fileindex[it->first]]<< " file."<< endreq;
        m_log<< MSG::INFO<< "Now the file entry is: "<< m_sfileentries[it->first]<< endreq;

        ris->setFileEntry(m_fileindex[it->first], m_sfileentries[it->first], true);
        m_streamentries[it->first] = ris->getEntry();
        rif = ris->getFiles().current();
      }
      else
      {
        //SetEntry 
        m_log << MSG::INFO<< "SetEntry manually: "<< " The sample: "<< it->first
              << "  File entry: "     << m_sfileentries[it->first]
              << "  Stream entry: "   << m_streamentries[it->first]
              << "  RootInputStream: "<< ris
              << "  Old entry: "      << ris->getEntry()<< endreq;
        ris->setFileEntry(m_fileindex[it->first], m_sfileentries[it->first], true);
        rif = ris->getFiles().current();
        //bool ifend = ris->setEntry(m_streamentries[it->first] ,true);
        //if(ifend == false) return StatusCode::FAILURE;
      }

      m_log<< MSG::INFO<< "Now the file is: "      << rif->filename()          << endreq;
      m_log<< MSG::INFO<< "Now the file index is: "<< m_fileindex[it->first]   << endreq;
      m_log<< MSG::INFO<< "Now the file entry is: "<< m_sfileentries[it->first]<< endreq;

      //Read in one PerSimHeader
      m_log<< MSG::INFO<< "Casting TObject to PerSimHeader."<< endreq;
      PerSimHeader* psh = (PerSimHeader*)((void*)(ris->obj()));
      m_log<< MSG::INFO<< "======= PerSimHeader: "<< psh<< endreq;
      if(!psh)
      { 
        m_log<< MSG::ERROR<< "Failed to retrive pointer to PerSimHeader : Terminated!"<< endreq;
        return StatusCode::FAILURE;
      }

      //Important parameters
      TimeStamp newtime(nextTime(m_maintau));
      //TimeStamp newtime(m_exp());
      TimeStamp pretime = m_lastTime;

      //Do a simple test of PerSimHeader
      vector<PerSimHit*>::size_type vpsh = 0;
      PerSimHitHeader* pshh = psh->hits;
      map<short int, PerSimHitCollection*> pshc = pshh->hitCollections;
      for(map<short int, PerSimHitCollection*>::iterator mpit = pshc.begin(); mpit != pshc.end(); mpit++)
      {
        vpsh += (mpit->second)->hits.size();
        m_log<< MSG::INFO<< "Detector : "<< mpit->first
             << "-- PerSimHit numbers : "<< (mpit->second)->hits.size()<< endreq;
      }
      
      if(0 == vpsh) //There is an empty event
      {
        m_log<< MSG::DEBUG<< "There is no SimHit in SimHeader!"<< endreq;
        m_sfileentries[it->first]++;
        m_streamentries[it->first]++;
        m_log<< MSG::DEBUG<< "Set new file entry as: "<< m_sfileentries[it->first]<< endreq;
        continue;
      }
      else
      {
        m_log<< MSG::INFO<< "PerSimHeader from "<< it->first<< " with "<< vpsh<< " PerSimHits!"<< endreq;
        m_log<< MSG::INFO<< "Generate new time["<< newtime.GetSeconds()
             << "sec] for sample ["<< it->first<< "]"<< endreq;

        //Change the TimeStamp of PerSimHeader(TimeStamp in Context, earliest, latest)
        pretime.Add(newtime);  
        TimeStamp temptime = psh->context.GetTimeStamp();
        psh->context.SetTimeStamp(pretime);
        TimeStamp eT  = temptime - psh->earliest;
        psh->earliest = pretime  - eT;
        TimeStamp Tl  = temptime - psh->latest;
        psh->latest   = pretime  - Tl;

        //Get the PerSimHeader's PerGenHeaders
        vector<PerRef> iheaders = psh->inputHeaders;
        for(vector<PerRef>::iterator prit = iheaders.begin(); prit != iheaders.end(); prit++)
        {
          m_gfileentries[it->first] = (*prit).entry();
          gris->setFileEntry(m_fileindex[it->first], m_gfileentries[it->first], true);
          grif = gris->getFiles().current();
          m_log<< MSG::INFO<< "Casting TObject to PerGenHeader."<< endreq;
          PerGenHeader* pgh = (PerGenHeader*)((void*)(gris->obj()));
          m_log<< MSG::INFO<< "======= PerGenHeader: "<< pgh<< endreq;
          if(!pgh)
          { 
            m_log<< MSG::ERROR<< "Failed to retrive pointer to PerGenHeader : Terminated!"<< endreq;
            return StatusCode::FAILURE;
          }
          
          //Convert PerGenHeader to GenHeader
          IConverter* gcnv = m_convSvc->converter(gh->classID());
          if(!gcnv) 
          {
            m_log<< MSG::ERROR<< "Failed to get GenHeader converter"<< endreq;
            return StatusCode::FAILURE;
          }
          RootIOTypedCnv<PerGenHeader, GenHeader>* ghCnv = 0;
          ghCnv = dynamic_cast< RootIOTypedCnv<PerGenHeader, GenHeader>* >(gcnv);
          if(!ghCnv)
          {
            m_log<< MSG::ERROR<< "Failed to cast GenHeader converter"<< endreq;
            return StatusCode::FAILURE;
          }

          sc = ghCnv->PerToTran(*pgh, *gh);
          if(sc.isFailure())
          {
            m_log<< MSG::ERROR<< "Failed to translate PerGenHeader to GenHeader!" 
                 << endreq;
            return sc;
          }
        }
        
        m_log<< MSG::INFO<< "GenHeaders: "<< gh<< endreq;

        //Convert PerSimHeader to SimHeader
        IConverter* cnv = m_convSvc->converter(sh->classID());
        if(!cnv) 
        {
          m_log<< MSG::ERROR<< "Failed to get SimHeader converter"<< endreq;
          return StatusCode::FAILURE;
        }
        RootIOTypedCnv<PerSimHeader, SimHeader>* shCnv = 0;
        shCnv = dynamic_cast< RootIOTypedCnv<PerSimHeader, SimHeader>* >(cnv);
        if(!shCnv)
        {
          m_log<< MSG::ERROR<< "Failed to cast SimHeader converter"<< endreq;
          return StatusCode::FAILURE;
        }
        sc = shCnv->PerToTran(*psh, *sh);
        if(sc.isFailure())
        {
          m_log<< MSG::ERROR<< "Failed to translate PerSimHeader to SimHeader!" 
               << endreq;
          return sc;
        }

        //Need to do a test of SimHeader(GenHeader)
        //Check the TimeStamps in SimHeader
        m_log<< MSG::INFO<< "SimHeader's TimeStamp: "<< sh->timeStamp().GetSec()<< "--"
             << sh->timeStamp().GetNanoSec()<< endreq;
        m_log<< MSG::INFO<< "SimHeader's Earliest: "<< sh->earliest().GetSec()<< "--"
             << sh->earliest().GetNanoSec()<< endreq;
        m_log<< MSG::INFO<< "SimHeader's Latest: "<< sh->latest().GetSec()<< "--"
             << sh->latest().GetNanoSec()<< endreq;

        //Check the SimHits in SimHeader
        const SimHitHeader* shh = sh->hits();
        m_log<< MSG::INFO<< "SimHitHeader's address: "<< shh<< endreq;
        map<short int, SimHitCollection*> shc = shh->hitCollection();
        for(map<short int, SimHitCollection*>::iterator mit = shc.begin(); mit != shc.end(); mit++)
        {
          vector<SimHit*> vsh = mit->second->collection();
          m_log<< MSG::INFO<< "In Det"<< mit->first<< " : SimHit's num = "<< vsh.size()<< endreq;
          /*double earliesttime = 0.0;
          double latesttime   = 0.0;
          for(vector<SimHit*>::iterator vit = vsh.begin(); vit != vsh.end(); vit++)
          {
            if(vit == vsh.begin() || earliesttime > (*vit)->hitTime())
            {
              earliesttime = (*vit)->hitTime();
            }
            if(vit == vsh.begin() || latesttime < (*vit)->hitTime()) 
            {
              latesttime = (*vit)->hitTime();
            }
            //m_log<< MSG::INFO<< "SimHit's time: "<< (*vit)->hitTime()<< endreq;           
          }
          
          TimeStamp detale(earliesttime/1.0e9);
          detale.Add(sh->timeStamp());
          TimeStamp detall(latesttime/1.0e9);
          detall.Add(sh->timeStamp());

          cout<< "Timeinterval between earliest: "<< (detale.GetSec() - sh->earliest().GetSec())
              << (detale.GetNanoSec() - sh->earliest().GetNanoSec())<< endl;
          cout<< "Timeinterval between latest: "<< (detall.GetSec() - sh->latest().GetSec())
          << (detall.GetNanoSec() - sh->latest().GetNanoSec())<< endl;*/
        }
        m_sfileentries[it->first]++;
        m_streamentries[it->first]++;
        m_lastTime = pretime;
    
        break;
      }
    }
    it++;
  }
  return StatusCode::SUCCESS;
}

double MixInputSvc::flatdis(double maximum)
{
  //Rndm::Numbers flat(m_ranSvc, Rndm::Flat(0.0, maximum));
  //if(flat) return flat();
  //else
  //{
  //m_log<< MSG::ERROR<< "Random Engine was not correct."<< endreq;
  //return 0.0;
  //}
  double tempnum = m_randomNum.Uniform(0.0, maximum);
  return tempnum;
}

double MixInputSvc::nextTime(double tau)
{
  //Rndm::Numbers exponential(m_ranSvc, Rndm::Exponential(tau));
  //if(exponential)  return  exponential();
  //else 
  //{ 
  //m_log<< MSG::ERROR<< "Random Engine was not correct."<< endreq;
  //return 0.0; 
  //}
  double tempnum = m_randomNum.Exp(tau);
  return tempnum;
}

---