Merge pull request #57 from mccumbermike/calo_eval_looping

blackcathj · blackcathj · commit 5c3eee008033 · 2015-09-19T00:05:53.000-04:00
Calorimeter evaluator modifications from Mike
diff --git a/simulation/g4simulation/g4eval/CaloEvaluator.C b/simulation/g4simulation/g4eval/CaloEvaluator.C
@@ -21,7 +21,7 @@
 #include <TFile.h>
 
 #include <iostream>
-//#include <vector>
+#include <set>
 //#include <algorithm>
 #include <cmath>
 
@@ -31,6 +31,8 @@ CaloEvaluator::CaloEvaluator(const string &name, const string &caloname, const s
   : SubsysReco(name),
     _caloname(caloname),
     _ievent(0),
+    _truth_trace_embed_flags(),
+    _reco_e_threshold(0.0), // 0 GeV before reco is traced
     _do_gpoint_eval(true),
     _do_gshower_eval(true),
     _do_tower_eval(true),
@@ -90,7 +92,7 @@ int CaloEvaluator::process_event(PHCompositeNode *topNode) {
   //---------------------------
 
   fillOutputNtuples(topNode);
-
+  
   //--------------------------------------------------
   // Print out the ancestry information for this event
   //--------------------------------------------------
@@ -368,6 +370,11 @@ void CaloEvaluator::fillOutputNtuples(PHCompositeNode *topNode) {
 	 iter != map.end(); 
 	 ++iter) {
       PHG4Particle* primary = iter->second;
+
+      if (!_truth_trace_embed_flags.empty()) {
+	if (_truth_trace_embed_flags.find(trutheval->get_embed(primary)) ==
+	    _truth_trace_embed_flags.end()) continue;
+      }
       
       float gparticleID = primary->get_track_id();
       float gflavor     = primary->get_pid();
@@ -466,6 +473,8 @@ void CaloEvaluator::fillOutputNtuples(PHCompositeNode *topNode) {
     for (rtiter = begin_end.first; rtiter !=  begin_end.second; ++rtiter) {
       RawTower *tower = rtiter->second;
 
+      if (tower->get_energy() < _reco_e_threshold) continue;
+      
       float towerid = tower->get_id();
       float ieta    = tower->get_bineta();
       float iphi    = tower->get_binphi();
@@ -573,6 +582,8 @@ void CaloEvaluator::fillOutputNtuples(PHCompositeNode *topNode) {
     for (unsigned int icluster = 0; icluster < clusters->size(); icluster++) {
       RawCluster *cluster = clusters->getCluster(icluster);
 
+      if (cluster->get_energy() < _reco_e_threshold) continue;
+      
       float clusterID = cluster->get_id();
       float ntowers   = cluster->getNTowers();
       float eta       = cluster->get_eta();
diff --git a/simulation/g4simulation/g4eval/CaloEvaluator.h b/simulation/g4simulation/g4eval/CaloEvaluator.h
@@ -13,6 +13,7 @@
 #include <TNtuple.h>
 #include <TFile.h>
 
+#include <set>
 #include <string>
 
 /// \class CaloEvaluator
@@ -36,12 +37,30 @@ class CaloEvaluator : public SubsysReco {
   int process_event(PHCompositeNode *topNode);
   int End(PHCompositeNode *topNode);
 
+  // forward trace only off of primaries with embed flags set
+  void add_truth_tracing_embed_flag(int flag) {
+    _truth_trace_embed_flags.insert(flag);
+  }
+
+  // backward trace only if reco meets an energy threshold requirement
+  void set_reco_tracing_energy_threshold(float thresh) {
+    _reco_e_threshold = thresh;
+  }
+
+  void set_do_gpoint_eval(bool b) {_do_gpoint_eval = b;}
+  void set_do_gshower_eval(bool b) {_do_gshower_eval = b;}
+  void set_do_tower_eval(bool b) {_do_tower_eval = b;}
+  void set_do_cluster_eval(bool b) {_do_cluster_eval = b;}
+  
  private:
 
   std::string _caloname;
   
   unsigned long _ievent;
 
+  std::set<int> _truth_trace_embed_flags;
+  float _reco_e_threshold;
+
   //----------------------------------
   // evaluator output ntuples
 
diff --git a/simulation/g4simulation/g4eval/CaloRawClusterEval.C b/simulation/g4simulation/g4eval/CaloRawClusterEval.C
@@ -64,22 +64,21 @@ std::set<PHG4Hit*> CaloRawClusterEval::all_truth_hits(RawCluster* cluster) {
   }
   
   std::set<PHG4Hit*> truth_hits;
-
+  
   // loop over all the clustered towers
   RawCluster::TowerConstRange begin_end = cluster->get_towers();
   RawCluster::TowerConstIterator iter;
   for (iter = begin_end.first; iter != begin_end.second; ++iter)
     { 
     RawTower* tower = _towers->getTower(iter->first);
     
-    std::set<PHG4Hit*> new_hits = _towereval.all_truth_hits(tower);
-    std::set<PHG4Hit*> union_hits;
+    std::set<PHG4Hit*> new_hits = get_rawtower_eval()->all_truth_hits(tower);
 
-    std::set_union(truth_hits.begin(),truth_hits.end(),
-		   new_hits.begin(),new_hits.end(),
-		   std::inserter(union_hits,union_hits.begin()));
-
-    std::swap(truth_hits,union_hits); // swap union into truth_particles    
+    for (std::set<PHG4Hit*>::iterator iter = new_hits.begin();
+	 iter != new_hits.end();
+	 ++iter) {
+      truth_hits.insert(*iter);
+    }
   }
 
   if (_do_cache) _cache_all_truth_hits.insert(make_pair(cluster,truth_hits));
@@ -106,13 +105,12 @@ std::set<PHG4Particle*> CaloRawClusterEval::all_truth_primaries(RawCluster* clus
     { 
     RawTower* tower = _towers->getTower(iter->first);
     std::set<PHG4Particle*> new_primaries = _towereval.all_truth_primaries(tower);
-    std::set<PHG4Particle*> union_primaries;
 
-    std::set_union(truth_primaries.begin(),truth_primaries.end(),
-		   new_primaries.begin(),new_primaries.end(),
-		   std::inserter(union_primaries,union_primaries.begin()));
-
-    std::swap(truth_primaries,union_primaries); // swap union into truth_particles    
+    for (std::set<PHG4Particle*>::iterator iter = new_primaries.begin();
+	 iter != new_primaries.end();
+	 ++iter) {
+      truth_primaries.insert(*iter);
+    }
   }
 
   if (_do_cache) _cache_all_truth_primaries.insert(make_pair(cluster,truth_primaries));
@@ -141,6 +139,7 @@ PHG4Particle* CaloRawClusterEval::max_truth_primary_by_energy(RawCluster* cluste
 
     PHG4Particle* primary = *iter;
     float e = get_energy_contribution(cluster,primary);
+    if (isnan(e)) continue;
     if (e > max_e) {
       max_e = e;
       max_primary = primary;      
@@ -154,15 +153,18 @@ PHG4Particle* CaloRawClusterEval::max_truth_primary_by_energy(RawCluster* cluste
 
 std::set<RawCluster*> CaloRawClusterEval::all_clusters_from(PHG4Particle* primary) { 
 
-  CaloTruthEval* trutheval = _towereval.get_truth_eval();
-  if (!trutheval->is_primary(primary)) return std::set<RawCluster*>();
-  
-  if ((_do_cache) &&
-      (_cache_all_clusters_from_primary.find(primary) != _cache_all_clusters_from_primary.end())) {
-    return _cache_all_clusters_from_primary[primary];
+  if (!get_truth_eval()->is_primary(primary)) return std::set<RawCluster*>();
+
+  primary = get_truth_eval()->get_primary_particle(primary);
+
+  if (_do_cache) {
+    std::map<PHG4Particle*,std::set<RawCluster*> >::iterator iter =
+      _cache_all_clusters_from_primary.find(primary);
+    if (iter != _cache_all_clusters_from_primary.end()) {
+      return iter->second;
+    }
   }
   
-
   std::set<RawCluster*> clusters;
   
   // loop over all the clusters
@@ -191,8 +193,9 @@ std::set<RawCluster*> CaloRawClusterEval::all_clusters_from(PHG4Particle* primar
 
 RawCluster* CaloRawClusterEval::best_cluster_from(PHG4Particle* primary) {
 
-  CaloTruthEval* trutheval = _towereval.get_truth_eval();
-  if (!trutheval->is_primary(primary)) return NULL;
+  if (!get_truth_eval()->is_primary(primary)) return NULL;
+
+  primary = get_truth_eval()->get_primary_particle(primary);
       
   if (_do_cache) {
     std::map<PHG4Particle*,RawCluster*>::iterator iter =
@@ -203,13 +206,14 @@ RawCluster* CaloRawClusterEval::best_cluster_from(PHG4Particle* primary) {
   }
   
   RawCluster* best_cluster = NULL;
-  float best_energy = 0.0;  
+  float best_energy = FLT_MAX*-1.0;  
   std::set<RawCluster*> clusters = all_clusters_from(primary);
   for (std::set<RawCluster*>::iterator iter = clusters.begin();
        iter != clusters.end();
        ++iter) {
     RawCluster* cluster = *iter;
     float energy = get_energy_contribution(cluster,primary);
+    if(isnan(energy)) continue;
     if (energy > best_energy) {
       best_cluster = cluster;
       best_energy = energy;
@@ -224,8 +228,10 @@ RawCluster* CaloRawClusterEval::best_cluster_from(PHG4Particle* primary) {
 // overlap calculations
 float CaloRawClusterEval::get_energy_contribution(RawCluster* cluster, PHG4Particle* primary) {
 
-  CaloTruthEval* trutheval = _towereval.get_truth_eval();
-  if (!trutheval->is_primary(primary)) return NAN;
+  if (!get_truth_eval()->is_primary(primary)) return NAN;
+
+  // reduce cache misses by using only pointer from PrimaryMap
+  primary = get_truth_eval()->get_primary_particle(primary);
   
   if (_do_cache) {
     std::map<std::pair<RawCluster*,PHG4Particle*>,float>::iterator iter =
@@ -241,7 +247,8 @@ float CaloRawClusterEval::get_energy_contribution(RawCluster* cluster, PHG4Parti
        iter != g4hits.end();
        ++iter) {
     PHG4Hit* g4hit = *iter;
-    if (g4hit->get_trkid() == primary->get_track_id()) {
+    PHG4Particle* candidate = get_truth_eval()->get_primary_particle(g4hit);
+    if (candidate->get_track_id() == primary->get_track_id()) {
       energy += g4hit->get_edep();
     }
   }
diff --git a/simulation/g4simulation/g4eval/CaloRawTowerEval.C b/simulation/g4simulation/g4eval/CaloRawTowerEval.C
@@ -105,9 +105,7 @@ std::set<PHG4Particle*> CaloRawTowerEval::all_truth_primaries(RawTower* tower) {
        iter != g4hits.end();
        ++iter) {
     PHG4Hit* g4hit = *iter;
-    PHG4Particle* particle = _truthinfo->GetHit( g4hit->get_trkid() );
-    PHG4Particle* primary = _trutheval.get_primary_particle( particle );
-    if (!primary) continue;
+    PHG4Particle* primary = get_truth_eval()->get_primary_particle( g4hit );    
     truth_primaries.insert(primary);
   }
 
@@ -131,12 +129,14 @@ PHG4Particle* CaloRawTowerEval::max_truth_primary_by_energy(RawTower* tower) {
   PHG4Particle* max_primary = NULL;
   float max_e = FLT_MAX*-1.0;
   std::set<PHG4Particle*> primaries = all_truth_primaries(tower);
+
   for (std::set<PHG4Particle*>::iterator iter = primaries.begin();
        iter != primaries.end();
        ++iter) {
 
     PHG4Particle* primary = *iter;
     float e = get_energy_contribution(tower,primary);
+    if (isnan(e)) continue;
     if (e > max_e) {
       max_e = e;
       max_primary = primary;      
@@ -151,6 +151,9 @@ PHG4Particle* CaloRawTowerEval::max_truth_primary_by_energy(RawTower* tower) {
 std::set<RawTower*> CaloRawTowerEval::all_towers_from(PHG4Particle* primary) { 
 
   if (!_trutheval.is_primary(primary)) return std::set<RawTower*>();
+
+  // use primary map pointer
+  primary = get_truth_eval()->get_primary_particle(primary);
   
   if (_do_cache) {
     std::map<PHG4Particle*,std::set<RawTower*> >::iterator iter =
@@ -189,7 +192,9 @@ std::set<RawTower*> CaloRawTowerEval::all_towers_from(PHG4Particle* primary) {
 RawTower* CaloRawTowerEval::best_tower_from(PHG4Particle* primary) {
 
   if (!_trutheval.is_primary(primary)) return NULL;
-      
+
+  primary = get_truth_eval()->get_primary_particle(primary);
+  
   if (_do_cache) {
     std::map<PHG4Particle*,RawTower*>::iterator iter =
       _cache_best_tower_from_primary.find(primary);
@@ -199,13 +204,14 @@ RawTower* CaloRawTowerEval::best_tower_from(PHG4Particle* primary) {
   }
   
   RawTower* best_tower = NULL;
-  float best_energy = 0.0;  
+  float best_energy = FLT_MAX*-1.0;  
   std::set<RawTower*> towers = all_towers_from(primary);
   for (std::set<RawTower*>::iterator iter = towers.begin();
        iter != towers.end();
        ++iter) {
     RawTower* tower = *iter;
     float energy = get_energy_contribution(tower,primary);
+    if (isnan(energy)) continue;
     if (energy > best_energy) {
       best_tower = tower;
       best_energy = energy;
@@ -221,6 +227,9 @@ RawTower* CaloRawTowerEval::best_tower_from(PHG4Particle* primary) {
 float CaloRawTowerEval::get_energy_contribution(RawTower* tower, PHG4Particle* primary) {
 
   if (!_trutheval.is_primary(primary)) return NAN;
+
+  // reduce cache misses by using only pointer from PrimaryMap
+  primary = get_truth_eval()->get_primary_particle(primary);
   
   if (_do_cache) {
     std::map<std::pair<RawTower*,PHG4Particle*>, float>::iterator iter =
@@ -236,7 +245,8 @@ float CaloRawTowerEval::get_energy_contribution(RawTower* tower, PHG4Particle* p
        iter != g4hits.end();
        ++iter) {
     PHG4Hit* g4hit = *iter;
-    if (g4hit->get_trkid() == primary->get_track_id()) {
+    PHG4Particle* candidate = get_truth_eval()->get_primary_particle(g4hit);
+    if (candidate->get_track_id() == primary->get_track_id()) {
       energy += g4hit->get_edep();
     }
   }
diff --git a/simulation/g4simulation/g4eval/CaloTruthEval.C b/simulation/g4simulation/g4eval/CaloTruthEval.C
@@ -23,7 +23,6 @@ CaloTruthEval::CaloTruthEval(PHCompositeNode* topNode,std::string caloname)
     _cache_all_truth_hits_g4particle(),
     _cache_get_primary_particle_g4particle(),
     _cache_get_primary_particle_g4hit(),
-    _cache_is_primary(),
     _cache_get_shower_from_primary(),
     _cache_get_shower_moliere_radius(),
     _cache_get_shower_energy_deposit() {
@@ -35,7 +34,6 @@ void CaloTruthEval::next_event(PHCompositeNode* topNode) {
   _cache_all_truth_hits_g4particle.clear();
   _cache_get_primary_particle_g4particle.clear();
   _cache_get_primary_particle_g4hit.clear();
-  _cache_is_primary.clear();
   _cache_get_shower_from_primary.clear();
   _cache_get_shower_moliere_radius.clear();
   _cache_get_shower_energy_deposit.clear();
@@ -72,7 +70,6 @@ std::set<PHG4Hit*> CaloTruthEval::all_truth_hits(PHG4Particle* particle) {
 
 PHG4Particle* CaloTruthEval::get_parent_particle(PHG4Hit* g4hit) {
 
-  // expensive call, but only because it gets called a ton (i think)
   PHG4Particle* particle = _truthinfo->GetHit( g4hit->get_trkid() );
   return particle;
 }
@@ -86,11 +83,14 @@ PHG4Particle* CaloTruthEval::get_primary_particle(PHG4Particle* particle) {
       return iter->second;
     }
   }
-  
-  PHG4Particle* returnval = particle;
-  if ((returnval->get_primary_id() != returnval->get_track_id()) ||
-      (returnval->get_primary_id() != -1)) {
-    returnval = _truthinfo->GetHit( particle->get_primary_id() );
+
+  // always report the primary from the Primary Map regardless if a
+  // primary from the full Map was the argument
+  PHG4Particle* returnval = NULL;
+  if (particle->get_primary_id() != -1) {
+    returnval = _truthinfo->GetPrimaryHit( particle->get_primary_id() );
+  } else {
+    returnval = _truthinfo->GetPrimaryHit( particle->get_track_id() );
   }
 
   if (_do_cache) _cache_get_primary_particle_g4particle.insert(make_pair(particle,returnval));
@@ -133,23 +133,13 @@ PHG4VtxPoint* CaloTruthEval::get_vertex(PHG4Particle* particle) {
 
 bool CaloTruthEval::is_primary(PHG4Particle* particle) {
 
-  if (_do_cache) {
-    std::map<PHG4Particle*,bool>::iterator iter =
-      _cache_is_primary.find(particle);
-    if (iter != _cache_is_primary.end()) {
-      return iter->second;
-    }
-  }
-
-  bool is_primary = false;
-  if (particle->get_primary_id() == particle->get_track_id()) {
-    is_primary = true;
-  } else if (particle->get_primary_id() == -1) {
-    is_primary = true;
+  bool is_primary = true;
+  if (!_truthinfo->GetPrimaryHit(particle->get_track_id())) {
+    // does the particle id appear in the primary map?
+    // this way either copy (in Map or Primary Map) will report correctly
+    is_primary = false;
   }
   
-  if (_do_cache) _cache_is_primary.insert(make_pair(particle,is_primary));
-  
   return is_primary;
 }
 
diff --git a/simulation/g4simulation/g4eval/CaloTruthEval.h b/simulation/g4simulation/g4eval/CaloTruthEval.h
@@ -48,7 +48,6 @@ class CaloTruthEval {
   std::map<PHG4Particle*,std::set<PHG4Hit*> > _cache_all_truth_hits_g4particle;
   std::map<PHG4Particle*,PHG4Particle*>       _cache_get_primary_particle_g4particle;
   std::map<PHG4Hit*,PHG4Particle*>            _cache_get_primary_particle_g4hit;
-  std::map<PHG4Particle*,bool>                _cache_is_primary;
   std::map<PHG4Particle*,std::set<PHG4Hit*> > _cache_get_shower_from_primary;
   std::map<PHG4Particle*,float>               _cache_get_shower_moliere_radius;
   std::map<PHG4Particle*,float>               _cache_get_shower_energy_deposit;
diff --git a/simulation/g4simulation/g4main/PHG4TruthInfoContainer.cc b/simulation/g4simulation/g4main/PHG4TruthInfoContainer.cc
diff --git a/simulation/g4simulation/g4main/PHG4TruthInfoContainer.h b/simulation/g4simulation/g4main/PHG4TruthInfoContainer.h
