g4simulation/g4main/PHG4ConsistencyCheck.cc

0001 #include "PHG4ConsistencyCheck.h"
0002
0003 #include "PHG4Hit.h"
0004 #include "PHG4HitContainer.h"
0005 #include "PHG4Particle.h"
0006 #include "PHG4TruthInfoContainer.h"
0007
0008 #include <phool/getClass.h>
0009
0010 #include <iostream>  // for operator<<, basic_ostream::opera...
0011 #include <map>       // for _Rb_tree_const_iterator, map<>::...
0012 #include <set>       // for set
0013 #include <utility>   // for pair
0014
0015 class PHCompositeNode;
0016
0017 using namespace std;
0018
0019 PHG4ConsistencyCheck::PHG4ConsistencyCheck(const std::string &name)
0020   : SubsysReco(name)
0021   , errorcnt(0)
0022 {
0023 }
0024
0025 int PHG4ConsistencyCheck::InitRun(PHCompositeNode * /*topNode*/)
0026 {
0027   return 0;
0028 }
0029
0030 int PHG4ConsistencyCheck::process_event(PHCompositeNode *topNode)
0031 {
0032   PHG4TruthInfoContainer *truthcont = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0033   if (!truthcont)
0034   {
0035     return 0;
0036   }
0037   PHG4TruthInfoContainer::ConstRange trange = truthcont->GetParticleRange();
0038   PHG4TruthInfoContainer::ConstIterator titer;
0039   int imax = 1000000;
0040   for (titer = trange.first; titer != trange.second; ++titer)
0041   {
0042     if (titer->first < imax)
0043     {
0044       imax = titer->first;
0045     }
0046   }
0047   cout << "min index: " << imax << endl;
0048   std::pair<std::map<int, int>::const_iterator, std::map<int, int>::const_iterator> embtrk_b_e = truthcont->GetEmbeddedTrkIds();
0049   std::map<int, int>::const_iterator embiter;
0050   for (embiter = embtrk_b_e.first; embiter != embtrk_b_e.second; ++embiter)
0051   {
0052     cout << "embedded trkid: " << embiter->first << endl;
0053   }
0054   PHG4HitContainer *ghit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_CEMC_E");
0055   if (ghit)
0056   {
0057     PHG4HitContainer::ConstIterator hit;
0058     PHG4HitContainer::ConstRange hit_begin_end = ghit->getHits();
0059     set<int> printpart;
0060     for (hit = hit_begin_end.first; hit != hit_begin_end.second; ++hit)
0061     {
0062       int trkid = hit->second->get_trkid();
0063       PHG4Particle *part = truthcont->GetParticle(trkid);
0064       if (!part)
0065       {
0066         hit->second->identify();
0067         cout << "could not locate geant particle " << trkid << " in G4HIT_CEMC_E" << endl;
0068         errorcnt++;
0069       }
0070       else
0071       {
0072         int primary_id = part->get_primary_id();
0073         if (truthcont->isEmbeded(primary_id) > 0)
0074         {
0075           if (printpart.find(primary_id) == printpart.end())
0076           {
0077             cout << "primary id " << primary_id << " is embedded" << endl;
0078             printpart.insert(primary_id);
0079             PHG4Particle *parta = truthcont->GetParticle(primary_id);
0080             parta->identify();
0081           }
0082         }
0083       }
0084     }
0085   }
0086   ghit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_SVTX");
0087   if (ghit)
0088   {
0089     PHG4HitContainer::ConstIterator hit;
0090     PHG4HitContainer::ConstRange hit_begin_end = ghit->getHits();
0091     for (hit = hit_begin_end.first; hit != hit_begin_end.second; ++hit)
0092     {
0093       int trkid = hit->second->get_trkid();
0094       PHG4Particle *part = truthcont->GetParticle(trkid);
0095       if (!part)
0096       {
0097         cout << "could not locate geant particle " << trkid << " in G4HIT_SVTX" << endl;
0098         errorcnt++;
0099       }
0100     }
0101   }
0102
0103   return 0;
0104 }