StripSumScatter.hpp

Go to the documentation of this file.

#ifndef STRIP_SUM_SCATTER_HPP
#define STRIP_SUM_SCATTER_HPP
 
#include "BeamFit2D.hpp"
#include "Constants.hpp"
#include "EventsSummary.hpp"
#include "FileSet.hpp"
#include "IOUtils.hpp"
#include "InitUtils.hpp"
#include "Normalization.hpp"
#include "PlottingUtils.hpp"
#include "RemixSim.hpp"
#include <Rtypes.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TCutG.h>
#include <TFile.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TKey.h>
#include <TLegend.h>
#include <TLine.h>
#include <TMath.h>
#include <TNamed.h>
#include <TROOT.h>
#include <TString.h>
#include <TSystem.h>
#include <TTree.h>
#include <algorithm>
#include <functional>
#include <iostream>
#include <limits>
#include <map>
#include <mutex>
#include <queue>
#include <stdexcept>
#include <string>
#include <thread>
#include <utility>
#include <vector>


enum TagCut {
  kTAG_PASS = 0,
  kCUT_UPSTREAM,   
  kCUT_JUMP,       
  kCUT_REAC_LEVEL, 
  kCUT_TAIL_RISE,  
  kCUT_RERISE,     
  kCUT_POST_ABOVE, 
  kCUT_CROSS,      
  kCUT_END_STRIP,  
  kCUT_CLIFF,      
  kCUT_OTHER_TAG,  
  kNTAGCUTS
};

extern const char *const kTagCutName[kNTAGCUTS];

struct TagThresholds {
  Double_t gate_nsigma = 0.0;   
  Double_t pileup_nsigma = 0.0; 
  Double_t noise_nsigma = 0.0;  
  Bool_t require_upstream = kTRUE;
  Double_t upstream_nsigma = 0.0;
  Double_t jump_nsigma = 0.0;
  Double_t smooth_nsigma = 0.0; 
  Int_t tail_fall_from_strip = 0;
  Double_t tail_rise_nsigma = 0.0;
  Double_t tail_return_nsigma = 0.0;
  Double_t tail_rerise_nsigma = 0.0;
  Double_t post_above_nsigma = 0.0;
  Int_t post_above_strips = 0;
  Double_t end_strip_nsigma = 0.0;
  Double_t cross_min_strip = 0.0;
  Double_t cliff_max = 0.0;
};

enum PreCut {
  kPRE_PASS = 0,
  kPRE_ALL_STRIPS, 
  kPRE_GATE,       
  kPRE_PILEUP,     
  kPRE_NOISE,      
  kPRE_SMOOTH,     
  kPRE_BOTH_MULT,  
  kNPRECUTS
};

extern const char *const kPreCutName[kNPRECUTS];

struct SelectionStep {
  enum Stage { kINPUT, kEVENT_LEVEL, kPER_STRIP, kOUTCOME };
  Stage stage;
  TString id;     
  TString name;   
  TString detail; 
  Bool_t on;      
  Int_t pre_cut;  
  Int_t tag_cut;  
};

namespace ScatterBuildRange {
const Double_t kXMin = 0.0;  
const Double_t kXMax = 40.0; 
const Double_t kYMin = 0.0;  
const Double_t kYMax =
    20.0; 
} // namespace ScatterBuildRange

struct BeamEllipses {
  BeamFit2D s0_s1;    
  BeamFit2D s1_s2;    
  BeamFit2D s1_grid;  
  BeamFit2D s16_s17;  
  BeamFit2D s15_s16;  
  Bool_t ok;          
  Bool_t use_s15_s16; 
};

struct BeamGate1D {
  Double_t mu;
  Double_t sigma;
  Bool_t ok;
  BeamGate1D() : mu(0.0), sigma(0.0), ok(kFALSE) {}
};

struct TraceEvt {
  Float_t leftdE[16];  
  Float_t rightdE[16]; 
  Float_t strip0dE;    
  Float_t strip17dE;   
  Float_t leftdE_adc[16];
  Float_t rightdE_adc[16];
  Float_t strip0_adc;
  Float_t strip17_adc;
  UInt_t reac_mask; 
  Bool_t beam_flat; 
  Int_t both_mult;  
  ULong64_t seed_ts;

  Double_t Total(Int_t strip) const {
    if (strip <= 0)
      return strip0dE;
    if (strip >= 17)
      return strip17dE;
    return Double_t(leftdE[strip - 1]) + Double_t(rightdE[strip - 1]);
  }

  Double_t TotalAdc(Int_t strip) const {
    if (strip <= 0)
      return strip0_adc;
    if (strip >= 17)
      return strip17_adc;
    return Double_t(leftdE_adc[strip - 1]) + Double_t(rightdE_adc[strip - 1]);
  }

  void Totals(Double_t *out) const {
    for (Int_t s = 0; s < 18; s++)
      out[s] = Total(s);
  }

  void TotalsAdc(Double_t *out) const {
    for (Int_t s = 0; s < 18; s++)
      out[s] = TotalAdc(s);
  }
};

struct SingleRunFitResult {
  BeamEllipses pure_beam; 
  BeamGate1D gate;        
  Bool_t ok;              
  SingleRunFitResult() : ok(kFALSE) {}
};

struct SingleRunFillResult {
  std::vector<TH2F *> scatters;
  std::vector<TraceEvt> reservoir; 
  Long64_t gated;                  
  Long64_t seen;                   
  Long64_t normed;
  std::vector<Long64_t> normed_at;
  std::vector<Long64_t> tagged;
  std::vector<std::vector<Long64_t>> tagged_var;
  std::vector<std::vector<Long64_t>> normed_var;
  std::vector<Long64_t> cut_counts;
  std::vector<Long64_t> pre_counts;
  SingleRunFillResult() : gated(0), seen(0), normed(0) {}
};

struct SimPop {
  TString file;  
  TString label; 
};

class StripSumScatter {
public:
  StripSumScatter();
  ~StripSumScatter();

  void Run();

  Bool_t Prepare();
  void SimTagReport();

  void DrawAllTaggedTraces(const TString &subdir);

  static TString CacheName();

  static void PlaneXY(const Double_t *total, Int_t reac, Double_t &x,
                      Double_t &y);
  static Int_t YLoOf(Int_t reac);
  static Int_t YHiOf(Int_t reac);

  static Double_t StripSigma(Int_t strip);
  static Double_t StripMean(Int_t strip);
  static void SetStripNoise(const Double_t *mean, const Double_t *sigma);
  static Double_t JumpMin(Int_t reac);
 
private:
  static Double_t s_stripSigma[18];
  static Double_t s_stripMean[18];
  std::map<Int_t, TH2F *> m_scatter;
  std::vector<TraceEvt> m_reservoir;
  // Normalization counts, merged over every run and persisted in the cache so
  // a cross section can be taken from it without a second pass over the data.
  Long64_t m_nSeen;
  Long64_t m_nNormed;
  std::vector<Long64_t> m_normedAt;
  std::vector<Long64_t> m_tagged;
  std::vector<std::vector<Long64_t>> m_taggedVar;
  std::vector<std::vector<Long64_t>> m_normedVar;
  std::vector<TString> m_variantNames;
  // Per-condition tag counts from the last fill, for the cut report; not
  // cached, so a cache load leaves them empty and writes no report.
  std::vector<Long64_t> m_cutCounts;
  std::vector<Long64_t> m_preCounts;
  void WriteCutReport() const;
  Double_t m_yLo[64];
  Double_t m_yHi[64];
 
  const TString kSimCacheName = "StripSumScatter_simcache.root";
 
  static Int_t ReacIndex(Int_t reac);
 
  Bool_t TryLoadCache(const TString &cacheName, const TString &fingerprint);
  void WriteCache(const TString &cacheName, const TString &fingerprint);
 
  // Fresh, empty scatters for the configured plane and build range.
  void AllocateScatters();
  void ReprojectFromReservoir();
  void FillScatters(const FileSet::GateGroups &groups);
 
  void PlotScatters();
 
  void InteractiveOverlay(Int_t reac);
 
  static void EnableEventBranches(TChain *chain);
  static Bool_t AllStripsFired(const EnergyView &ev);
  static Bool_t IsPureBeam(const EnergyView &ev, const BeamEllipses &be);
  static Bool_t MeasureBeamNoise(TChain *chain, const BeamEllipses &be,
                                 Double_t *strip_mean, Double_t *strip_sigma);
  // Strips 1..reac-1 within BEAM_UPSTREAM_NSIGMA of the beam, or the
  // requirement is off. Shared by the tag and its per-strip denominator.
  static Bool_t BeamUpstreamOf(const EnergyView &ev, Int_t reac);
  static Bool_t BeamUpstreamOf(const EnergyView &ev, Int_t reac,
                               const TagThresholds &T);
  static Double_t SumRange(const Double_t *total, Int_t lo, Int_t hi);

  static Bool_t PassesGate(const BeamGate1D &gate, const EnergyView &ev,
                           Int_t strip, Double_t nsigma);

  static Bool_t PassesEllipse(const BeamFit2D &gate, const EnergyView &ev,
                              Int_t sx, Int_t sy, Double_t nsigma);

  static BeamFit2D FindBeamGate(TChain *chain, Int_t sx, Int_t sy,
                                const TString &tag, const TString &subdir,
                                Double_t nsigma);
  static BeamGate1D FindStripGate(TChain *chain, Int_t strip,
                                  const TString &tag, const TString &subdir,
                                  Double_t nsigma, Double_t center);
 
  static void DrawTraceSet(const std::vector<TGraph *> &traces, Int_t color);
  static TGraph *TraceFromTotal(const Double_t *total);
  struct TraceClass {
    std::vector<TGraph *> traces;
    TString label;
    Int_t color;
  };
  static void DrawRegionTraces(const TString &save_name, const TString &subdir,
                               const std::vector<TGraph *> &beam,
                               const std::vector<TraceClass> &classes,
                               Double_t y_min, Double_t y_max,
                               const char *y_title,
                               Bool_t beam_sigma_measured = kTRUE);
  static void DrawRegionTraces(const TString &save_name, const TString &subdir,
                               const std::vector<TGraph *> &beam,
                               const std::vector<TGraph *> &aa,
                               const std::vector<TGraph *> &an, Double_t y_min,
                               Double_t y_max, const char *y_title,
                               Bool_t beam_sigma_measured = kTRUE,
                               const char *an_label = nullptr);
  static Int_t SimClassColor(const TString &base);
 
  static void DrawRegionMeanTraces(const TString &save_name,
                                   const TString &subdir,
                                   const std::vector<TGraph *> &beam,
                                   const std::vector<TGraph *> &aa,
                                   const std::vector<TGraph *> &an,
                                   Double_t y_min, Double_t y_max,
                                   const char *y_title);
 
  static void TraceYRange(const std::vector<TGraph *> &beam,
                          const std::vector<TGraph *> &aa,
                          const std::vector<TGraph *> &an, Double_t &y_min,
                          Double_t &y_max);
  static SingleRunFitResult FitRunGates(Int_t key, const TString &label,
                                        TChain *chain);
  static SingleRunFillResult FillRunScatters(Int_t key, const TString &label,
                                             TChain *chain,
                                             const BeamGate1D &gate,
                                             const BeamEllipses &runBeam);
  struct EventLevelVerdict {
    Bool_t gate = kFALSE;
    Bool_t pileup = kFALSE;
    Bool_t noise = kFALSE;
    Bool_t smooth_ok = kTRUE;
    Bool_t event = kFALSE;
    Int_t pre_cut = -1;
  };
  static EventLevelVerdict JudgeEventLevel(const BeamGate1D &gate,
                                           const EnergyView &ev,
                                           Int_t gate_strip, Double_t step_z,
                                           const TagThresholds &T);
  static Bool_t AnyVariantKeeps(
      const BeamGate1D &gate, const EnergyView &ev, Int_t gate_strip,
      const std::vector<std::pair<TString, TagThresholds>> &variants,
      const std::vector<Bool_t> &varies_event, Bool_t nom_event,
      std::vector<Bool_t> &var_event);
  static void TagUnderVariants(
      const EnergyView &ev, Double_t step_z, Int_t kReacMin, Int_t kReacMax,
      const std::vector<std::pair<TString, TagThresholds>> &variants,
      const std::vector<Bool_t> &var_event, std::vector<Bool_t> &pass,
      SingleRunFillResult &res);
  static Int_t TagNominal(const EnergyView &ev, const Double_t totals[18],
                          Int_t kReacMin, Int_t kReacMax, Int_t nReacStrips,
                          std::vector<Bool_t> &pass, SingleRunFillResult &res,
                          UInt_t &mask);
  static void AllocateRunScatters(Int_t run, Int_t kReacMin, Int_t kReacMax,
                                  Int_t kXBins, Int_t kYBins,
                                  std::vector<TH2F *> &scatters);
  static TCutG *PromptCut(TCanvas *c, const char *name, const char *label);
  static void SaveRegionCuts(Int_t reac, TCutG *cut_an, TCutG *cut_aa);
  static TCutG *LoadRegionCut(const char *name, Int_t reac);

  static void SavitzkyGolay(const Double_t *in, Double_t *out);
 
  // Build a TGraph from Savitzky-Golay-smoothed per-strip totals. Input is
  // the raw normed array; smoothing is applied internally before graph build.
  static TGraph *SmoothedTraceFromTotal(const Double_t *total);
 
  static TString
  SimFingerprint(const std::vector<RemixSim::SimFileSpec> &specs);
  static TString BuildFingerprint(const FileSet::GateGroups &groups);
  static void YBounds(Double_t *y_lo, Double_t *y_hi);
  static TString PrettyLabel(const TString &tag);
 
public:
  static Int_t ResolveTag(const EnergyView &ev,
                          const std::vector<Bool_t> &pass);
  static Bool_t Condition(TagCut c, const EnergyView &ev, Int_t reac,
                          const TagThresholds &T);
  static Bool_t ConditionActive(TagCut c, const TagThresholds &T);
  static TagCut RejectReason(const EnergyView &ev, Int_t reac);
  static TagCut RejectReason(const EnergyView &ev, Int_t reac,
                             const TagThresholds &T);
  static TagCut TailReason(const EnergyView &ev, Int_t reac,
                           const TagThresholds &T);
  static Bool_t IsPileup(const EnergyView &ev);
  static Bool_t IsNoise(const EnergyView &ev);
  static Bool_t IsPileup(const EnergyView &ev, Double_t nsigma);
  static Bool_t IsNoise(const EnergyView &ev, Double_t nsigma);
  static Double_t MaxStepNSigma(const EnergyView &ev);
  static Bool_t IsSmooth(const EnergyView &ev, Double_t nsigma);
  static TagThresholds NominalThresholds();
  static std::vector<std::pair<TString, TagThresholds>> ThresholdVariants();

  static std::vector<SelectionStep> DescribeSelection();
 
private:
  static Bool_t SimBeamGains(Double_t *gain);
  // Run n indexed tasks on a pool of workers, pulling from a shared queue.
  static void RunIndexedParallel(Int_t n, Int_t workers,
                                 const std::function<void(Int_t)> &task);
  // One fill task: a processed events file and the index of the gate group
  // (run) it belongs to.
  struct FillTask {
    Int_t run_idx;
    TString path;
  };
  // The fill phase's state shared by the workers: the per-file results, the
  // merge order and the totals merged so far.
  struct FillPhase {
    std::vector<SingleRunFillResult> fills;
    std::vector<Bool_t> filled;
    Int_t next_merge = 0;
    std::mutex merge_mutex;
    Int_t nBeamKept = 0;
    Long64_t totalGated = 0;
    Long64_t totalSeen = 0;
  };
  // Phase 1: one run's beam ellipses and strip gate; the first groups of the
  // epoch draw their gate figures.
  void FitGateRun(Int_t i, std::vector<SingleRunFitResult> &fits,
                  const std::vector<Int_t> &runOrder,
                  const std::vector<TString> &labelVec,
                  const std::vector<TChain *> &chainVec);
  // Merge one file's fill into the totals: the scatters, the reservoir
  // (pure-beam events are capped globally, not per task) and the counts.
  void MergeRunFill(Int_t t, std::vector<SingleRunFillResult> &fills,
                    Int_t &nBeamKept, Long64_t &totalGated,
                    Long64_t &totalSeen);
  // Phase 2: one file of the fill, then the finished prefixes merged into
  // the totals in task order as each prefix finishes, freeing on the spot.
  void FillScatterFile(Int_t t, const std::vector<FillTask> &tasks,
                       const std::vector<Int_t> &runOrder,
                       const std::vector<TString> &labelVec,
                       const std::vector<SingleRunFitResult> &fits,
                       FillPhase &phase);
  static void SimTotal(const RemixSim::Event &e, const Double_t *gain,
                       Double_t *total);
  static TGraph *SimPopScatter(const TString &file, Int_t reac,
                               const Double_t *gain, Long64_t max_points);
  static std::vector<TGraph *>
  SimPopTraces(const TString &file, const Double_t *gain, Long64_t max_traces);
  void SimTraceOverlay();
  static Bool_t LoadSimCache(const TString &fp,
                             std::map<Int_t, std::vector<TGraph *>> &by_strip);
  static void
  WriteSimCache(const TString &fp,
                const std::map<Int_t, std::vector<TGraph *>> &by_strip);
  void SimOverlay();
};
 
#endif