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 <TEllipse.h>
#include <TFile.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TKey.h>
#include <TLegend.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 <vector>


struct GateSpec {
  Int_t sx; 
  Int_t sy; 
};

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 s16_s17;  
  BeamFit2D s15_s16;  
  Bool_t ok;          
  Bool_t use_s15_s16; 
};

struct TraceEvt {
  Float_t total[18];     
  Float_t total_adc[18]; 
  Float_t long_au[18];  
  Float_t short_au[18]; 
 
  UInt_t reac_mask; 
  Bool_t beam_flat; 
  Int_t both_mult;  
  ULong64_t seed_ts;
};

struct SingleRunFitResult {
  BeamEllipses pure_beam;              
  std::vector<BeamFit2D> series_gates; 
  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;
  SingleRunFillResult() : gated(0), seen(0), normed(0) {}
};

struct SimPop {
  TString file;  
  TString label; 
};

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

  void Run();

  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 JumpSigma(Int_t strip);
  static Double_t StripSigma(Int_t strip);
  static void SetJumpSigma(const Double_t *sigma);
  static void SetStripSigma(const Double_t *sigma);
  static Double_t JumpMin(Int_t reac);
 
private:
  static Double_t s_jumpSigma[18];
  static Double_t s_stripSigma[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;
  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();
  // Fill the scatters from the reservoir alone. The reservoir keeps every
  // tagged event, so a plane or axes change needs only this and not a
  // 25-minute pass over the events files.
  void ReprojectFromReservoir();
  void FillScatters(const std::vector<Int_t> &runOrder,
                    std::map<Int_t, TChain *> &chains);
 
  void PlotScatters();
 
  // Diagnostic (PLOT_PARITY_REJECTED_GRID): grid #DeltaE of events that pass
  // the cheap pre-tag cuts and are then rejected by the parity cut. Saves both
  // an a.u. view (decoded grid) and an ADC view (raw trigger channel), each on
  // a log-y axis. Reads the Grid branch directly (a self-contained pass), so it
  // works regardless of whether the scatter cache has already been filled.
  void PlotParityRejectedGrid(const std::vector<Int_t> &run_order,
                              std::map<Int_t, TChain *> &chains);
 
  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);
  // Sigma-clipped width of each strip-to-strip difference and of each strip's
  // deposit over a capped sample of `chain`, after the cheap pre-tag cuts.
  // False when too few events survive to measure them.
  static Bool_t MeasureBeamNoise(TChain *chain, Double_t *jump_sigma,
                                 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 IsPileup(const EnergyView &ev);
  static Bool_t IsNoise(const EnergyView &ev);
  static Bool_t IsOffbeam(const EnergyView &ev);
  static Bool_t IsParityAsymmetric(const EnergyView &ev);
  static Double_t SumRange(const Double_t *total, Int_t lo, Int_t hi);
  static std::vector<GateSpec> ActiveGates();
 
  static Bool_t PassesGate(const BeamFit2D &gate, const EnergyView &ev,
                           Int_t sx, Int_t sy);
 
  static BeamFit2D FindBeamGate(TChain *chain, Int_t sx, Int_t sy,
                                const std::vector<GateSpec> &prior_specs,
                                const std::vector<BeamFit2D> &prior_gates,
                                const TString &tag, const TString &subdir);
 
  static void DrawTraceSet(const std::vector<TGraph *> &traces, Int_t color);
  void DrawAltDecodeRegionTraces(Int_t reac, TCutG *cutAn, TCutG *cutAa);
  static TGraph *TraceFromTotal(const Float_t *total);
  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);
 
  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 run, TChain *chain,
              const std::vector<GateSpec> &activeGates);
  static SingleRunFillResult FillRunScatters(
      Int_t run, TChain *chain, const std::vector<GateSpec> &activeGates,
      const std::vector<BeamFit2D> &runGates, const BeamEllipses &runBeam);
  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 SmoothTrace(const Double_t *in, Double_t *out, Int_t width);
 
  // Savitzky-Golay smoothing: 3rd-degree polynomial, half-window of 2
  // (5-point convolution). Uses standard SG coefficients [-3,12,17,12,-3]/35.
  // At edges, the window shrinks and coefficients are renormalised.
  static void SavitzkyGolay(const Double_t *in, Double_t *out);
 
  // CFD-style trigger finder: locate the first strip whose beam-subtracted
  // signal (td[s]-1) exceeds both a fraction of the trace peak and a multiple
  // of the beam sigma. Returns the strip index, or -1 if no trigger fires.
  static Int_t FindTrigger(const Double_t *td, const Double_t *base,
                           Double_t beam_sigma);
 
  // 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 Float_t *total);
 
  void ClusterVarHists(Int_t reac, TCutG *cut_aa, TCutG *cut_an,
                       const TString &subdir);
 
  static TString
  SimFingerprint(const std::vector<RemixSim::SimFileSpec> &specs);
  static TString BuildFingerprint(const std::vector<Int_t> &run_order,
                                  std::map<Int_t, TChain *> &chains);
  static void YBounds(Double_t *y_lo, Double_t *y_hi);
  static TString PrettyLabel(const TString &tag);
 
public:
  static Bool_t PassesReaction(const EnergyView &ev, Int_t reac);
 
private:
  static Bool_t SimBeamGains(Double_t *gain);
  static void SimTotal(const Float_t *left, const Float_t *right,
                       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