Some updates to PDFastSimPAR + using new OpDetBTR

larsim/MCCheater/PhotonBackTracker.cc

@@ -314,6 +314,8 @@ namespace cheat {
       std::sort(
         timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), pairSort);
     }
+    // auto const& timeSDPMap =
+    //   FindOpDetBTR(fWireReadoutGeom->OpDetFromOpChannel(opHit.OpChannel()))->timePDclockSDPsMap();
 
     //This section is a hack to make comparisons work right.
     std::vector<sim::SDP> dummyVec;
@@ -322,15 +324,18 @@ namespace cheat {
     auto start_timePair_P = &start_timePair;
     auto end_timePair_P = &end_timePair;
     //First interesting iterator.
+    // auto mapFirst = timeSDPMap.lower_bound(start_time);
     auto mapFirst = std::lower_bound(timePDclockSDPMap_SortedPointers.begin(),
                                      timePDclockSDPMap_SortedPointers.end(),
                                      start_timePair_P,
                                      pairSort);
     //Last interesting iterator.
+    // auto mapLast = timeSDPMap.upper_bound(end_time);
     auto mapLast =
       std::upper_bound(mapFirst, timePDclockSDPMap_SortedPointers.end(), end_timePair_P, pairSort);
 
     for (auto& mapitr = mapFirst; mapitr != mapLast; ++mapitr)
+      // for (auto& sdp : mapitr->second)
       for (auto& sdp : (*mapitr)->second)
         retVec.push_back(&sdp);
 

larsim/PhotonPropagation/PDFastSimPAR_module.cc

@@ -125,6 +125,7 @@ namespace phot {
       DP VUVHits{Name("VUVHits"), Comment("Configuration for UV visibility parameterization")};
       ODP VISHits{Name("VISHits"),
                   Comment("Configuration for visibile visibility parameterization")};
+      fhicl::Atom<bool> Verbose{Name("Verbose"), Comment("Print verbose information"), false};
     };
     using Parameters = art::EDProducer::Table<Config>;
 
@@ -139,6 +140,16 @@ namespace phot {
     void AddOpDetBTR(std::vector<sim::OpDetBacktrackerRecord>& opbtr,
                      std::vector<int>& ChannelMap,
                      const sim::OpDetBacktrackerRecord& btr) const;
+    void SimpleAddOpDetBTR(
+      // std::vector<sim::OpDetBacktrackerRecord>& opbtr,
+      std::map<int, sim::OBTRHelper>& opbtr,
+      std::vector<int>& ChannelMap,
+      size_t channel,
+      int trackID,
+      int time,
+      double pos[3],
+      double edeposit,
+      int num_photons = 1);
 
     bool isOpDetInSameTPC(geo::Point_t const& ScintPoint, geo::Point_t const& OpDetPoint) const;
     std::vector<geo::Point_t> opDetCenters() const;
@@ -325,6 +336,10 @@ namespace phot {
     auto opbtr = std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();
     auto phlit_ref = std::make_unique<std::vector<sim::SimPhotonsLite>>();
     auto opbtr_ref = std::make_unique<std::vector<sim::OpDetBacktrackerRecord>>();
+
+    //Helpers holding maps
+    std::map<int, sim::OBTRHelper> opbtr_helper, opbtr_helper_ref;
+
     auto& dir_phlitcol(*phlit);
     auto& ref_phlitcol(*phlit_ref);
     // SimPhotons
@@ -363,7 +378,18 @@ namespace phot {
     int num_fastdp = 0;
     int num_slowdp = 0;
 
+    mf::LogTrace("PDFastSimPAR") << "Creating SimPhotonsLite/SimPhotons from " << (*edeps).size()
+                                 << " energy deposits\n";
+
     for (auto const& edepi : *edeps) {
+
+      if (!(num_points % 1000)) {
+        mf::LogTrace("PDFastSimPAR")
+          << "SimEnergyDeposit: " << num_points << " " << edepi.TrackID() << " " << edepi.Energy()
+          << "\nStart: " << edepi.Start() << "\nEnd: " << edepi.End()
+          << "\nNF: " << edepi.NumFPhotons() << "\nNS: " << edepi.NumSPhotons()
+          << "\nSYR: " << edepi.ScintYieldRatio() << "\n";
+      }
       num_points++;
 
       int nphot_fast = edepi.NumFPhotons();
@@ -445,40 +471,79 @@ namespace phot {
 
           // SimPhotonsLite case
           if (fUseLitePhotons) {
-            sim::OpDetBacktrackerRecord tmpbtr(channel);
             if (ndetected_fast > 0 && fDoFastComponent) {
               int n = ndetected_fast;
               num_fastdp += n;
               for (int i = 0; i < n; ++i) {
                 // calculates the time at which the photon was produced
                 double dtime = edepi.StartT() + fScintTime->fastScintTime();
                 if (fIncludePropTime) dtime += transport_time[i];
+
                 int time = static_cast<int>(std::round(dtime));
-                if (Reflected)
+                if (Reflected) {
                   ++ref_phlitcol[channel].DetectedPhotons[time];
-                else
+                  SimpleAddOpDetBTR(
+                    // *opbtr_ref, PDChannelToSOCMapReflect, channel, trackID, time, pos, edeposit, 1);
+                    opbtr_helper_ref,
+                    PDChannelToSOCMapReflect,
+                    channel,
+                    trackID,
+                    time,
+                    pos,
+                    edeposit,
+                    1);
+                }
+                else {
                   ++dir_phlitcol[channel].DetectedPhotons[time];
-                tmpbtr.AddScintillationPhotons(trackID, time, 1, pos, edeposit);
+                  SimpleAddOpDetBTR(
+                    // *opbtr, PDChannelToSOCMapDirect, channel, trackID, time, pos, edeposit, 1);
+                    opbtr_helper,
+                    PDChannelToSOCMapDirect,
+                    channel,
+                    trackID,
+                    time,
+                    pos,
+                    edeposit,
+                    1);
+                }
               }
             }
             if (ndetected_slow > 0 && fDoSlowComponent) {
               int n = ndetected_slow;
               num_slowdp += n;
               for (int i = 0; i < n; ++i) {
+                // calculates the time at which the photon was produced
                 double dtime = edepi.StartT() + fScintTime->slowScintTime();
                 if (fIncludePropTime) dtime += transport_time[ndetected_fast + i];
                 int time = static_cast<int>(std::round(dtime));
-                if (Reflected)
+                if (Reflected) {
                   ++ref_phlitcol[channel].DetectedPhotons[time];
-                else
+                  SimpleAddOpDetBTR(
+                    // *opbtr_ref, PDChannelToSOCMapReflect, channel, trackID, time, pos, edeposit, 1);
+                    opbtr_helper_ref,
+                    PDChannelToSOCMapReflect,
+                    channel,
+                    trackID,
+                    time,
+                    pos,
+                    edeposit,
+                    1);
+                }
+                else {
                   ++dir_phlitcol[channel].DetectedPhotons[time];
-                tmpbtr.AddScintillationPhotons(trackID, time, 1, pos, edeposit);
+                  SimpleAddOpDetBTR(
+                    // *opbtr, PDChannelToSOCMapDirect, channel, trackID, time, pos, edeposit, 1);
+                    opbtr_helper,
+                    PDChannelToSOCMapDirect,
+                    channel,
+                    trackID,
+                    time,
+                    pos,
+                    edeposit,
+                    1);
+                }
               }
             }
-            if (Reflected)
-              AddOpDetBTR(*opbtr_ref, PDChannelToSOCMapReflect, tmpbtr);
-            else
-              AddOpDetBTR(*opbtr, PDChannelToSOCMapDirect, tmpbtr);
           }
           // SimPhotons case
           else {
@@ -531,10 +596,29 @@ namespace phot {
                                  << "\ndetected fast photons: " << num_fastdp
                                  << ", detected slow photons: " << num_slowdp;
 
+    mf::LogDebug("PDFastSimPAR") << "Number of entries in opbtrs";
+    for (auto& iopbtr : *opbtr) {
+      mf::LogDebug("PDFastSimPAR")
+        << "OpDet: " << iopbtr.OpDetNum() << " " << iopbtr.timePDclockSDPsMap().size();
+    }
+    mf::LogDebug("PDFastSimPAR") << "Number of entries in opbtrs refelected";
+    for (auto& iopbtr : *opbtr_ref) {
+      mf::LogDebug("PDFastSimPAR")
+        << "OpDet: " << iopbtr.OpDetNum() << " " << iopbtr.timePDclockSDPsMap().size();
+    }
+
     if (fUseLitePhotons) {
+
+      for (auto& iopbtr : opbtr_helper) {
+        opbtr->emplace_back(iopbtr.second);
+      }
+
       event.put(move(phlit));
       event.put(move(opbtr));
       if (fDoReflectedLight) {
+        for (auto& iopbtr : opbtr_helper_ref) {
+          opbtr_ref->emplace_back(iopbtr.second);
+        }
         event.put(move(phlit_ref), "Reflected");
         event.put(move(opbtr_ref), "Reflected");
       }
@@ -570,6 +654,26 @@ namespace phot {
     }
   }
 
+  void PDFastSimPAR::SimpleAddOpDetBTR( //std::vector<sim::OpDetBacktrackerRecord>& opbtr,
+    std::map<int, sim::OBTRHelper>& opbtr,
+    std::vector<int>& ChannelMap,
+    size_t channel,
+    int trackID,
+    int time,
+    double pos[3],
+    double edeposit,
+    int num_photons)
+  {
+    if (ChannelMap[channel] < 0) {
+      ChannelMap[channel] = opbtr.size();
+      // opbtr.push_back(sim::OpDetBacktrackerRecord(channel));
+      opbtr.emplace(channel, channel);
+    }
+    // size_t idtest = ChannelMap[channel];
+    // opbtr.at(idtest).AddScintillationPhotonsToMap(trackID, time, num_photons, pos, edeposit);
+    opbtr.at(channel).AddScintillationPhotonsToMap(trackID, time, num_photons, pos, edeposit);
+  }
+
   //......................................................................
   // calculates number of photons detected given visibility and emitted number of photons
   void PDFastSimPAR::detectedNumPhotons(std::vector<int>& DetectedNumPhotons,

larsim/PhotonPropagation/ScintTimeTools/ScintTimeLAr.cc

@@ -18,8 +18,8 @@ namespace phot {
     , SDTime{pset.get<double>("SlowDecayTime", 0.)}
     , FRTime{pset.get<double>("FastRisingTime", 0.)}
     , FDTime{pset.get<double>("FastDecayTime", 0.)}
-    , fNoFastRisingTime((FRTime >= 1.e-8) && (FRTime != -1.))
-    , fNoSlowRisingTime((SRTime >= 1.e-8) && (SRTime != -1.))
+    , fNoFastRisingTime(pset.get<bool>("NoFastRisingTime", false))
+    , fNoSlowRisingTime(pset.get<bool>("NoSlowRisingTime", false))
   {
     if (LogLevel >= 1) {
       std::cout << "ScintTimeLAr Tool configure:" << std::endl;
@@ -53,7 +53,9 @@ namespace phot {
     // If we care about simulating the rising time, it would be better
     // to simulate random numbers according a general distribution
     double d = (tau1 + tau2) / tau2;
+    int ncalls = 0;
     while (1) {
+      ++ncalls;
       double ran1 = fUniformGen->fire();
       double ran2 = fUniformGen->fire();
       double t = -tau2 * std::log(1 - ran1);