FillTruth.cxx

Go to the documentation of this file.

// $Id: FillTruth.cxx,v 1.18 2008/03/10 19:58:06 rustem Exp $

// C++
#include <vector>

// ROOT
#include "TClonesArray.h"

// MINOS
#include "Conventions/BeamType.h"
#include "DataUtil/infid.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCReweight/NeugenWeightCalculator.h"
#include "MCReweight/ReweightHelpers.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "StandardNtuple/NtpStRecord.h"
#include "TruthHelperNtuple/NtpTHEvent.h"
#include "TruthHelperNtuple/NtpTHShower.h"
#include "TruthHelperNtuple/NtpTHTrack.h"

// Local
#include "PhysicsNtuple/Default.h"
#include "PhysicsNtuple/Factory.h"
#include "FillTruth.h"

REGISTER_ANP_OBJECT(AlgStore,FillTruth)

CVSID("$Id: FillTruth.cxx,v 1.18 2008/03/10 19:58:06 rustem Exp $");

using namespace std;

//-----------------------------------------------------------------------------------------
Anp::FillTruth::FillTruth()
   :fStdHep(true),
    fCheck(true),
    fHepKey(-1),
    fNTruth(0),
    fNStHep(0),
    fNRecoEvent(0),
    fNRecoTrack(0),
    fNRecoShower(0)
{    
}

//-----------------------------------------------------------------------------------------
Anp::FillTruth::~FillTruth()
{
   cout << "   Number of truths        " << fNTruth << endl
        << "   Number of stdheps       " << fNStHep << endl
        << "   Number of event truths  " << fNRecoEvent << endl
        << "   Number of track truths  " << fNRecoTrack << endl
        << "   Number of shower truths " << fNRecoShower << endl;
}

//-----------------------------------------------------------------------------------------
bool Anp::FillTruth::Run(Record &record, TObject *ptr)
{
   //
   // Get and fill Monte-Carlo truth and StdHep information
   //

   NtpStRecord *ntprec = dynamic_cast<NtpStRecord *>(ptr);
   if(!ntprec)
   {
      const MomNavigator *mom = dynamic_cast<const MomNavigator *> (ptr);
      if(mom)
      {
         ntprec = dynamic_cast<NtpStRecord *>(mom -> GetFragment("NtpStRecord")); 
      }
      else
      {
         MSG("FillAlg", Msg::kError) << "Failed to find MomNavigator pointer" << endl;
         return false;
      }      
   }

   if(!ntprec)
   {
      MSG("FillAlg", Msg::kError) << "Failed to find NtpStRecord pointer" << endl;
      return false;
   }

   //
   // Do nothing if running on data events
   //
   if(ntprec -> GetVldContext() -> GetSimFlag() != SimFlag::kMC)
   {
      return true;
   }

   const TClonesArray *stdhep_array = ntprec -> stdhep;
   if(!stdhep_array)
   {
      MSG("FillAlg", Msg::kWarning) << "TClonesArray of NtpMCStdHep doesn't exist" << endl;      
      return true;
   }

   //
   // If fStdHep option is set than fill StdHep blocks
   //
   if(fStdHep)
   {
      FillStdHep(record, *stdhep_array);
   }

   const TClonesArray *mc_array = ntprec -> mc;
   if(!mc_array)
   {
      MSG("FillAlg", Msg::kError) << "TClonesArray of NtpMCTruth doesn't exist" << endl;      
      return false;
   }

   if(!ntprec -> thevt || !ntprec -> thtrk || !ntprec -> thshw)
   {
      MSG("FillAlg", Msg::kError) << "One or more TClonesArray of NtpTH doesn't exist" << endl;      
      return false;
   }   

   const vector<const NtpTHEvent *>  thevt_vec = Anp::FillRecoTruth<NtpTHEvent> (ntprec -> thevt);
   const vector<const NtpTHTrack *>  thtrk_vec = Anp::FillRecoTruth<NtpTHTrack> (ntprec -> thtrk);
   const vector<const NtpTHShower *> thshw_vec = Anp::FillRecoTruth<NtpTHShower>(ntprec -> thshw);

   MSG("FillAlg", Msg::kDebug)
      << endl
      << "   size of NtpMCTruth array  = " << mc_array -> GetEntries() << endl
      << "   size of NtpTHEvent array  = " << thevt_vec.size() << endl
      << "   size of NtpTHTrack array  = " << thtrk_vec.size() << endl
      << "   size of NtpMCStdHep array = " << stdhep_array -> GetEntries() << endl
      << "   size of NtpTHShower array = " << thshw_vec.size() << endl;
   
   const int nentries = mc_array -> GetEntries();
   for(int index = 0; index < nentries; ++index)
   {
      const NtpMCTruth *truth_ptr  = dynamic_cast<const NtpMCTruth *>(mc_array -> At(index));
      if(!truth_ptr)
      {
         MSG("FillAlg", Msg::kWarning) << "No mc truth at " << index << endl;
         continue;
      }      
      if(index != truth_ptr -> index)
      {
         MSG("FillAlg", Msg::kWarning) << "Wrong array index" << endl;
         continue;
      }

      Truth truth = FillTruth::Fill(*truth_ptr);

      FillTruth::Fill(truth, thevt_vec);
      FillTruth::Fill(truth, thtrk_vec);
      FillTruth::Fill(truth, thshw_vec);
    
      //
      // Get parameters for neugen reweighting
      //
      MCEventInfo mcinfo;
      ReweightHelpers::MCEventInfoFilla(&mcinfo, ntprec, truth.TruthIndex());

      truth.init_state = mcinfo.initial_state;
      truth.nucleus    = mcinfo.nucleus;
      truth.had_fs     = mcinfo.had_fs;
      
      //
      // Interaction is in true fiducial volume?
      //
      truth.isfid = infid(ntprec -> GetVldContext() -> GetDetector(), SimFlag::kMC,
                          truth.VtxX(), truth.VtxY(), truth.VtxZ());
      
      record.Add(truth);
   }

   MSG("FillAlg", Msg::kDebug) 
      << endl
      << "   size of Truth vector  = " << record.GetNTruths() << endl
      << "   size of StdHep vector = " << record.GetNStdHeps() << endl;
   
   if(fCheck && !IsConsistent(record))
   {
      return false;
   }

   //
   // Add PDG code for best matching hep block to reconstucted tracks
   //
   if(fHepKey >= 0)
   {
      for(TrackIterator itrack = record.TrackBegIterator(); itrack != record.TrackEndIterator(); ++itrack)
      {
         StdHepIter ihep = record.FindStdHep(*itrack);
         if(ihep != record.StdHepEnd())
         {
            itrack -> Add(fHepKey, ihep -> GetIdHEP());
         }
         else
         {
            MSG("FillAlg", Msg::kWarning) << "Failed to find StdHep block" << endl;
         }
      }
   }
   
   return true;
}

//----------------------------------------------------------------------
void Anp::FillTruth::Config(const Registry &reg)
{
   Anp::Read(reg, "FillTruthCheck",  fCheck);
   Anp::Read(reg, "FillTruthStdHep", fStdHep);

   if(reg.KeyExists("PrintConfig"))
   {
      cout << "FillTruth::Config" << endl
           << "   StdHep = " << fStdHep << endl
           << "   Check = " << fCheck << endl
           << "   HepKey = " << fHepKey << endl;
   }
}

//-----------------------------------------------------------------------------------------
const Anp::Truth Anp::FillTruth::Fill(const NtpMCTruth &ntpmc) const
{
   //
   // Fill Truth, this is just like Truth constructor, put it outside
   // of Truth class to avoid dependency on NtpMCTruth class.
   //

   Truth truth;
   
   truth.index        = ntpmc.index;
   truth.stdhep[0]    = ntpmc.stdhep[0];
   truth.stdhep[1]    = ntpmc.stdhep[1];
   truth.inu          = ntpmc.inu;
   truth.inunoosc     = ntpmc.inunoosc;
   truth.itg          = ntpmc.itg;
   truth.iboson       = ntpmc.iboson;
   truth.iresonance   = ntpmc.iresonance;
   truth.iaction      = ntpmc.iaction;
   truth.istruckq     = ntpmc.istruckq;
   truth.iflags       = ntpmc.iflags;
   truth.ndigu        = ntpmc.ndigu;
   truth.ndigv        = ntpmc.ndigv;
   truth.tphu         = ntpmc.tphu;
   truth.tphv         = ntpmc.tphv;
   truth.a            = ntpmc.a;
   truth.z            = ntpmc.z;
   truth.sigma        = ntpmc.sigma;
   truth.sigmadiff    = ntpmc.sigmadiff;
   truth.x            = ntpmc.x;
   truth.y            = ntpmc.y;
   truth.q2           = ntpmc.q2;
   truth.w2           = ntpmc.w2;
   truth.emfrac       = ntpmc.emfrac;
   truth.vtxx         = ntpmc.vtxx;
   truth.vtxy         = ntpmc.vtxy;
   truth.vtxz         = ntpmc.vtxz;   
   
   for(short i = 0; i < 4; ++i)
   {
      truth.p4neu[i]      = ntpmc.p4neu[i];
      truth.p4neunoosc[i] = ntpmc.p4neunoosc[i];
      truth.p4shw[i]      = ntpmc.p4shw[i];
      truth.p4tgt[i]      = ntpmc.p4tgt[i];
   }   
   
   truth.flux.fluxrun     = ntpmc.flux.fluxrun;
   truth.flux.fluxevtno   = ntpmc.flux.fluxevtno;
   truth.flux.ndecay      = ntpmc.flux.ndecay;
   truth.flux.ntype       = ntpmc.flux.ntype;
   truth.flux.ptype       = ntpmc.flux.ptype;
   truth.flux.tptype      = ntpmc.flux.tptype;
   truth.flux.tgen        = ntpmc.flux.tgen;

   truth.flux.ndxdznear   = ntpmc.flux.ndxdznear;
   truth.flux.ndydznear   = ntpmc.flux.ndydznear;
   truth.flux.nenergynear = ntpmc.flux.nenergynear;
   truth.flux.nwtnear     = ntpmc.flux.nwtnear;
   truth.flux.ndxdzfar    = ntpmc.flux.ndxdzfar;
   truth.flux.ndydzfar    = ntpmc.flux.ndydzfar;
   truth.flux.nenergyfar  = ntpmc.flux.nenergyfar;
   truth.flux.nwtfar      = ntpmc.flux.nwtfar;
   truth.flux.vx          = ntpmc.flux.vx;
   truth.flux.vy          = ntpmc.flux.vy;
   truth.flux.vz          = ntpmc.flux.vz;
   truth.flux.pdpx        = ntpmc.flux.pdpx;
   truth.flux.pdpy        = ntpmc.flux.pdpy;
   truth.flux.pdpz        = ntpmc.flux.pdpz;
   truth.flux.ppdxdz      = ntpmc.flux.ppdxdz;
   truth.flux.ppdydz      = ntpmc.flux.ppdydz;
   truth.flux.pppz        = ntpmc.flux.pppz;
   truth.flux.ppenergy    = ntpmc.flux.ppenergy;
   truth.flux.ppvx        = ntpmc.flux.ppvx;
   truth.flux.ppvy        = ntpmc.flux.ppvy;
   truth.flux.ppvz        = ntpmc.flux.ppvz;
   truth.flux.necm        = ntpmc.flux.necm;
   truth.flux.nimpwt      = ntpmc.flux.nimpwt;
   truth.flux.tvx         = ntpmc.flux.tvx;
   truth.flux.tvy         = ntpmc.flux.tvy;
   truth.flux.tvz         = ntpmc.flux.tvz;
   truth.flux.tpx         = ntpmc.flux.tpx;
   truth.flux.tpy         = ntpmc.flux.tpy;
   truth.flux.tpz         = ntpmc.flux.tpz;

   ++fNTruth;

   return truth;
}
//-----------------------------------------------------------------------------------------
void Anp::FillTruth::FillStdHep(Record &record, const TClonesArray &array) const
{   
   //
   // Fill all StdHep blocks from Monte-Carlo simulation
   //

   const int nstdhep = array.GetEntries();

   // initiliaze StdHep vector all at once to speed up 
   // insertion of new elements into vector
   std::vector<Anp::StdHep> &svec = record.GetStdHep();   
   svec.clear();

   // initiliaze stdhep elements
   svec.insert(svec.begin(), nstdhep, Anp::StdHep());

   for(int index = 0; index < nstdhep; ++index)
   {
      const NtpMCStdHep *stdhep_ptr  = dynamic_cast<const NtpMCStdHep *>(array.At(index));
      if(!stdhep_ptr)
      {
         MSG("FillAlg", Msg::kError) << "No stdhep at " << index << " out of " << nstdhep << endl;
         continue;
      }
      
      if(int(stdhep_ptr -> index) != index) 
      {
         MSG("FillAlg", Msg::kError) << "Wrong StdHep index" << endl;
         continue;
      }

      StdHep &stdhep = svec[index];

      stdhep.index     = stdhep_ptr -> index;
      stdhep.mc_index  = stdhep_ptr -> mc;
      stdhep.parent[0] = stdhep_ptr -> parent[0];
      stdhep.parent[1] = stdhep_ptr -> parent[1];
      stdhep.child[0]  = stdhep_ptr -> child[0];
      stdhep.child[1]  = stdhep_ptr -> child[1];
      stdhep.IstHEP    = stdhep_ptr -> IstHEP;
      stdhep.IdHEP     = stdhep_ptr -> IdHEP;
      stdhep.mass      = stdhep_ptr -> mass;

      stdhep.p4[0] = stdhep_ptr -> p4[0];
      stdhep.p4[1] = stdhep_ptr -> p4[1];
      stdhep.p4[2] = stdhep_ptr -> p4[2];
      stdhep.p4[3] = stdhep_ptr -> p4[3];

      stdhep.vtx[0] = stdhep_ptr -> vtx[0];
      stdhep.vtx[1] = stdhep_ptr -> vtx[1];
      stdhep.vtx[2] = stdhep_ptr -> vtx[2];
      stdhep.vtx[3] = stdhep_ptr -> vtx[3];

      ++fNStHep;
   }
}

//-----------------------------------------------------------------------------------------
void Anp::FillTruth::Fill(Truth &truth, const vector<const NtpTHEvent *> &vec)
{
   //
   // Fill all "reconstructed event to truth" records for this Truth object
   //

   unsigned int count = 0;

   for(vector<const NtpTHEvent *>::const_iterator it = vec.begin(); it != vec.end(); ++it)
   {
      const NtpTHEvent *th_event = *it;

      if(th_event -> neumc != truth.TruthIndex() && !truth.IsChild(th_event -> neustdhep))
      {
         continue;
      }       

      const TruthReco truth_reco(TruthReco::kEvent,
                                 th_event -> index,
                                 th_event -> neumc,
                                 th_event -> neustdhep,
                                 th_event -> purity, 
                                 th_event -> completeall, 
                                 th_event -> completeslc);
      
      //
      // Fill event truth record in this Truth
      //
      truth.fReco.push_back(truth_reco);

      ++count;
   }

   fNRecoEvent += count;

   MSG("FillAlg", Msg::kDebug) << "MC truth record with index " << truth.TruthIndex()
                               << " matches "<< count << " reconstructed events" << endl;
}

//-----------------------------------------------------------------------------------------
void Anp::FillTruth::Fill(Truth &truth, const vector<const NtpTHTrack *> &vec)
{
   //
   // Fill all "reconstructed track to truth" records for this Truth object
   //

   unsigned int count = 0;

   for(vector<const NtpTHTrack *>::const_iterator it = vec.begin(); it != vec.end(); ++it)
   {
      const NtpTHTrack *th_track = *it;

      if(th_track -> neumc != truth.TruthIndex() && 
         !truth.IsChild(th_track -> neustdhep) &&
         !truth.IsChild(th_track -> trkstdhep))
      {
         continue;
      }
      
      const TruthReco truth_reco(TruthReco::kTrack,
                                 th_track -> index,
                                 th_track -> neumc,
                                 th_track -> trkstdhep,
                                 th_track -> purity, 
                                 th_track -> completeall, 
                                 th_track -> completeslc);

      //
      // Fill event truth record in this Truth
      //
      truth.fReco.push_back(truth_reco);

      ++count;
   } 

   fNRecoTrack += count;

   MSG("FillAlg", Msg::kDebug) << "MC truth record with index " << truth.TruthIndex()
                               << " matches "<< count << " reconstructed tracks" << endl;  
}

//-----------------------------------------------------------------------------------------
void Anp::FillTruth::Fill(Truth &truth, const vector<const NtpTHShower *> &vec)
{
   //
   // Fill all "reconstructed shower to truth" records for this Truth object
   //

   unsigned int count = 0;

   for(vector<const NtpTHShower *>::const_iterator it = vec.begin(); it != vec.end(); ++it)
   {
      const NtpTHShower *th_shower = *it;

      if(th_shower -> neumc != truth.TruthIndex() && !truth.IsChild(th_shower -> neustdhep))
      {
         continue;
      }
      
      const TruthReco truth_reco(TruthReco::kShower,
                                 th_shower -> index,
                                 th_shower -> neumc,
                                 th_shower -> neustdhep,
                                 th_shower -> purity, 
                                 th_shower -> completeall, 
                                 th_shower -> completeslc);

      //
      // Fill event truth record in this Truth
      //
      truth.fReco.push_back(truth_reco);

      ++count;
   } 

   fNRecoShower += count;

   MSG("FillAlg", Msg::kDebug) << "MC truth record with index " << truth.TruthIndex()
                               << " matches "<< count << " reconstructed showers" << endl;  
}

//-----------------------------------------------------------------------------------------
bool Anp::FillTruth::IsConsistent(const Record &record) const
{
   //
   // Run a loop over all truth records in a snarl. For each event and track truth record
   // check that there exists a corresponding reco event or track. Small fraction of 
   // tracks are not included with any reco event. Such tracks for now are ignored.
   //

   for(TruthIter truth = record.TruthBeg(); truth != record.TruthEnd(); ++truth)
   {
      for(TruthRecoIter reco = truth -> RecoBeg(); reco != truth -> RecoEnd(); ++reco)
      {
         short found = 0;
         
         if(record.StdHepBeg() != record.StdHepEnd())
         {
            for(StdHepIter hep = record.StdHepBeg(); hep != record.StdHepEnd(); ++hep)
            {
               if(reco -> StdHepIndex() == hep -> StdHepIndex())
               {
                  ++found;
               }
            }

            if(found != 1)
            {
               MSG("FillAlg", Msg::kError) << "Failed to match StdHep to TruthReco" << endl;
               reco -> Print();
               return false;
            }
         }

         found = std::count(record.EventBeg(), record.EventEnd(), *reco);
         if(reco -> IsEvent() && record.EventBeg() != record.EventEnd() && found != 1)
         {
            MSG("FillAlg", Msg::kError) << "No match for event to TruthReco: found = " << found << endl;
            reco -> Print();
            return false;
         }
         
         found = std::count(record.TrackBeg(), record.TrackEnd(), *reco);
         if(reco -> IsTrack() && record.TrackBeg() != record.TrackEnd() && found != 1)
         {
            MSG("FillAlg", Msg::kError) << "Failed to match track to TruthReco" << endl;
            reco -> Print();
            return false;
         }
      }
   }

   //
   // The loop below runs over all reco events and tracks. For each reco event or track 
   // run a loop over all Truth records. Check that for each reco event or track there 
   // is only one matching TruthEvent or TruthTrack object.
   //

   for(EventIter event = record.EventBeg(); event != record.EventEnd(); ++event)
   {
      short found = 0;
      for(TruthIter truth = record.TruthBeg(); truth != record.TruthEnd(); ++truth)
      { 
         //
         // check that some Truth record matches event event
         //
         found += std::count(truth -> RecoBeg(), truth -> RecoEnd(), *event);
      }

      //
      // every event has its truth info filled?
      //
      if(found != 1)
      {
         MSG("FillAlg", Msg::kError) << "Failed to find truth info for a reco event" << endl;
         return false;
      }
   }
   
   for(TrackIter track = record.TrackBeg(); track != record.TrackEnd(); ++track)
   {
      short found = 0;
      for(TruthIter truth = record.TruthBeg(); truth != record.TruthEnd(); ++truth)
      { 
         //
         // check that some Truth record matches track event
         //
         found += std::count(truth -> RecoBeg(), truth -> RecoEnd(), *track);
      }
      
      //
      // every track has its truth info filled?
      //
      if(found != 1)
      {
         MSG("FillAlg", Msg::kError) << "Found " << found <<" truth records for reco track" << endl;
         return false;
      }
   }

   for(ShowerIter shower = record.ShowerBeg(); shower != record.ShowerEnd(); ++shower)
   {
      short found = 0;
      for(TruthIter truth = record.TruthBeg(); truth != record.TruthEnd(); ++truth)
      { 
         //
         // check that some Truth record matches shower event
         //
         found += std::count(truth -> RecoBeg(), truth -> RecoEnd(), *shower);
      }
      
      //
      // every shower has its truth info filled?
      //
      if(found != 1)
      {
         MSG("FillAlg", Msg::kError) << "Found " << found <<" truth records for reco shower" << endl;
         return false;
      }
   }     
   
   return true;
}

//-----------------------------------------------------------------------------------------
template<class T> 
   vector<const T *> Anp::FillRecoTruth(const TClonesArray *array)
{
   //
   // Collect and return records that match true and reconstructed objects
   //

   vector<const T *> tvec;

   if(!array)
   {
      cerr << "FillRecoTruth<T> - TClonesArray of reco-to-truth records does not exist" << endl;
      return tvec;
   }

   const int nentries = array -> GetEntries();
   for(int index = 0; index < nentries; ++index)
   {
      const T *treco = dynamic_cast<const T *>(array -> At(index));
      if(!treco)
      {
         cerr << "FillRecoTruth<T> - dynamic_cast failed at index: "<< index << endl;
         continue;
      }      

      if(treco -> index != index)
      {
         cerr << "FillRecoTruth<T> - wrong index at index: "<< index << endl;
         continue;
      } 
      
      tvec.push_back(treco);
   }
   
   return tvec;
}

---