Merge pull request #675 from brownd1978/compfix

Fix MCTrajectory reference updating. Add a comparison plot for CrvSteps
Mu2e · Dec 13, 2021 · d582a8d · d582a8d
2 parents 2cdee13 + 3c163dd
commit d582a8d
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 16 deletions.
diff --git a/Compression/src/CompressDetStepMCs_module.cc b/Compression/src/CompressDetStepMCs_module.cc
@@ -157,6 +157,8 @@ class mu2e::CompressDetStepMCs : public art::EDProducer {
   std::unique_ptr<CrvStepCollection> _newCrvSteps;
   std::map<InstanceLabel, std::unique_ptr<StepPointMCCollection> > _newStepPointMCs;
   std::unique_ptr<MCTrajectoryCollection> _newMCTrajectories;
+  // temporary storage for MCTrajectories
+  std::vector<MCTrajectory> _newMCTrajs;
 
   // To create art::Ptrs to the new SimParticles and GenParticles,
   // we need their art::ProductIDs and art::EDProductGetters
@@ -443,7 +445,7 @@ void mu2e::CompressDetStepMCs::compressSimParticles(const art::Event& event) {
     if (_keepNGenerations >= 0) {
       // Go through the particles we are keeping and see if any parents are not there
       for (auto& i_keptSimPart : _simParticlesToKeep[i_product_id]) {
-	
+
 	art::Ptr<mu2e::SimParticle> i_childPtr = i_keptSimPart;
 	art::Ptr<mu2e::SimParticle> i_parentPtr = i_childPtr->parent();
 	while (i_parentPtr) {
@@ -584,15 +586,15 @@ void mu2e::CompressDetStepMCs::compressMCTrajectories(const art::Event& event) {
         const SimParticleSet& alreadyKeptSimParts = simPartsToKeep.second;
         for (const auto& alreadyKeptSimPart : alreadyKeptSimParts) {
           if (mcTrajectory.second.sim() == alreadyKeptSimPart) {
-            MCTrajectory newMCTrajectory(mcTrajectory.second);
-            _newMCTrajectories->insert(std::make_pair(mcTrajectory.second.sim(), newMCTrajectory));
+            _newMCTrajs.emplace_back(mcTrajectory.second);
+            if(_debugLevel>0 ) std::cout << "Inserting new MCTrajectory with key " << mcTrajectory.second.sim() << std::endl;
           }
         }
       }
     }
     else if (_mcTrajectoryCompressionLevel == mu2e::CompressionLevel::kNoCompression) {
-      MCTrajectory newMCTrajectory(mcTrajectory.second);
-      _newMCTrajectories->insert(std::make_pair(mcTrajectory.second.sim(), newMCTrajectory));
+      _newMCTrajs.emplace_back(mcTrajectory.second);
+      if(_debugLevel>0 ) std::cout << "Inserting new MCTrajectory with key " << mcTrajectory.second.sim() << " and recording it " << std::endl;
       recordSimParticle(mcTrajectory.second.sim());
     }
     else {
@@ -615,18 +617,17 @@ void mu2e::CompressDetStepMCs::updateStepPointMCs() {
 }
 
 void mu2e::CompressDetStepMCs::updateMCTrajectories() {
-  for (auto& i_mcTrajectory : *_newMCTrajectories) {
-    const auto& oldSimPtr = i_mcTrajectory.second.sim();
+  for (auto& i_mcTrajectory : _newMCTrajs) {
+    const auto& oldSimPtr = i_mcTrajectory.sim();
     art::Ptr<mu2e::SimParticle> newSimPtr = _simPtrRemap.at(oldSimPtr);
     if(_debugLevel>0) {
       std::cout << "Updating SimParticlePtr in MCTrajectory from " << oldSimPtr << " to " << newSimPtr << std::endl;
     }
-    i_mcTrajectory.second.sim() = newSimPtr;
-    auto mcTrajPair = _newMCTrajectories->extract(i_mcTrajectory.first);
-    mcTrajPair.key() = newSimPtr;
-    _newMCTrajectories->insert(std::move(mcTrajPair));
-    //    i_mcTrajectory.first = newSimPtr;
+    // overwrite the internal ptr and insert with the new key
+    i_mcTrajectory.sim() = newSimPtr;
+    _newMCTrajectories->insert(std::make_pair(i_mcTrajectory.sim(), i_mcTrajectory));
   }
+  _newMCTrajs.clear(); // clear the temporary storage
 }
 
 void mu2e::CompressDetStepMCs::recordSimParticle(const art::Ptr<mu2e::SimParticle>& sim_ptr) {

diff --git a/Validation/inc/ValCrvStep.hh b/Validation/inc/ValCrvStep.hh
@@ -20,11 +20,12 @@ namespace mu2e {
 
   private:
     std::string _name;
-    
+
     TH1D* _hVer;
     TH1D* _hN;
     TH1D* _hb;
     TH1D* _ht;
+    TH1D* _ht2;
     TH1D* _hlt;
     TH1D* _hE;
     TH1D* _hlE;

diff --git a/Validation/src/ValCrvStep.cc b/Validation/src/ValCrvStep.cc
@@ -7,6 +7,7 @@ int mu2e::ValCrvStep::declare(art::TFileDirectory tfs) {
   _hN = tfs.make<TH1D>( "NStep", "N Steps", 101, -9.5, 1000.5);
   _hb = tfs.make<TH1D>( "bar", "bar number", 100, -0.5, 5503.5);
   _ht = tfs.make<TH1D>( "t", "time", 100, 0.0, 2000.0);
+  _ht2 = tfs.make<TH1D>( "t2", "time", 100, -100.0, 100.0e3);
   _hlt = tfs.make<TH1D>( "lt", "log_{10} time;log_{10} time (ns)", 100, 0.0, 10.0);
   _hE = tfs.make<TH1D>( "E", "Energy",100, 0.0, 0.2);
   _hlE = tfs.make<TH1D>( "lE", "log10(Energy)",100, -9.0, 2.0);
@@ -22,14 +23,15 @@ int mu2e::ValCrvStep::declare(art::TFileDirectory tfs) {
 int mu2e::ValCrvStep::fill(const mu2e::CrvStepCollection & coll,
 				art::Event const& event) {
 
-  // increment this by 1 any time the defnitions of the histograms or the 
+  // increment this by 1 any time the defnitions of the histograms or the
   // histogram contents change, and will not match previous versions
-  _hVer->Fill(1.0);
+  _hVer->Fill(2.0);
 
-  _hN->Fill(coll.size()); 
+  _hN->Fill(coll.size());
   for(auto cs : coll) {
     _hb->Fill(cs.barIndex().asInt());
     _ht->Fill(cs.startTime());
+    _ht2->Fill(cs.startTime());
     _hlt->Fill(log10(std::max(cs.startTime(),1e-5)));
     _hE->Fill(cs.visibleEDep());
     _hlE->Fill( (cs.visibleEDep()>0.0 ? log10(cs.visibleEDep()) : -10.0)  );