update ntuple dumper, track max radius, scint hit word, get rid of NaN's, add...

update ntuple dumper, track max radius, scint hit word, get rid of NaN's, add flag to ostoore only blinded events

PhysicsAnalysis/NtupleDumper/scripts/faser_ntuple_maker.py

@@ -44,8 +44,17 @@ parser.add_argument("--isMC", action='store_true',
 parser.add_argument("--partial", action='store_true',
                     help="Allow partial input files")
 
+parser.add_argument("--trigFilt", action='store_true',
+                    help="apply trigger event filter ")
+parser.add_argument("--scintFilt", action='store_true',
+                    help="apply scintillator event filter ")
+parser.add_argument("--trackFilt", action='store_true',
+                    help="apply track event filter ")
+
 parser.add_argument("--unblind", action='store_true',
                     help="Don't apply signal blinding (default: do)")
+parser.add_argument("--onlyblind", action='store_true',
+                    help="Only store events that were blinded (default: don't)")
  
 parser.add_argument("--fluka", action='store_true',
                     help="Add FLUKA weights to ntuple")
@@ -112,8 +121,8 @@ if filepath.is_dir():
 
         filestr = str(flist[0].resolve())
         # Use proper EOS file path here?
-        if filestr[:4] == '/eos':
-            filestr = f"root://eospublic.cern.ch/{filestr}"
+#        if filestr[:4] == '/eos':
+#            filestr = f"root://eospublic.cern.ch/{filestr}"
         filelist.append(filestr)
     # End of loop over segments
 
@@ -140,6 +149,7 @@ if filepath.is_dir():
     print(f"Last  = {lastseg}")
     print(f"Args  = {args.tag}")
     print(f"Blind = {not args.unblind}")
+    print(f"OnlyBlinded = {args.onlyblind}")
 
     # Find any tags
     tagstr = firststem.replace(f"{firstfaser}-{firstshort}-{runstr}-{firstseg}", "")
@@ -243,7 +253,7 @@ if args.isMC:
         acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile))
 
 else:
-    acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind)))
+    acc.merge(NtupleDumperAlgCfg(ConfigFlags, outfile, DoBlinding=(not args.unblind), OnlyBlinded=args.onlyblind, DoScintFilter = args.scintFilt, DoTrackFilter = args.trackFilt, DoTrigFilter = args.trigFilt))
 
 if not args.verbose:
     from AthenaConfiguration.ComponentFactory import CompFactory

PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h

@@ -52,12 +52,17 @@ private:
 
   bool waveformHitOK(const xAOD::WaveformHit* hit) const;
   void clearTree() const;
-  void setNaN() const;
+  void clearTrackTruth() const;
   void addBranch(const std::string &name,float* var);
   void addBranch(const std::string &name,unsigned int* var);
   void addWaveBranches(const std::string &name, int nchannels, int first);
   void FillWaveBranches(const xAOD::WaveformHitContainer &wave) const;
   void addCalibratedBranches(const std::string &name, int nchannels, int first);
+  double radius(const Acts::Vector3 &position) const;
+  std::optional<Acts::Vector3> positionAtMaxRadius(
+      const EventContext &ctx, const Acts::BoundTrackParameters &startParameters,
+      double targetPosition, Acts::NavigationDirection navDir = Acts::forward) const;
+
 
   ServiceHandle <ITHistSvc> m_histSvc;
 
@@ -103,7 +108,13 @@ private:
 
   BooleanProperty m_useIFT               { this, "UseIFT", false, "Use IFT tracks" };
   BooleanProperty m_doCalib              { this, "DoCalib", true, "Fill calibrated calorimeter quantities" };
-  BooleanProperty m_doBlinding        { this, "DoBlinding", true, "Blinding will not output events with Calo signal > 10 GeV e-" };
+  BooleanProperty m_doBlinding           { this, "DoBlinding", true, "Blinding will not output events with no veto signal adn Calo signal > 20 GeV e-" };
+  BooleanProperty m_onlyBlinded          { this, "OnlyBlinded", false, "Only events that would be blinded are saved" };
+  BooleanProperty m_doTrigFilter         { this, "DoTrigFilter", false, "Only events that pass trigger cuts are passed" };
+  BooleanProperty m_doScintFilter        { this, "DoScintFilter", false, "Only events that pass scintillator coincidence cuts are passed" };
+  BooleanProperty m_doTrackFilter        { this, "DoTrackFilter", false, "Only events that have >= 1 long track are passed" };
+  BooleanProperty m_doScintOrTrackFilter { this, "DoScintOrTrackFilter", false, "Only events that pass scintillator coincidence cuts or have >= 1 long track are passed" };
+
   BooleanProperty m_useFlukaWeights      { this, "UseFlukaWeights", false, "Flag to weight events according to value stored in HepMC::GenEvent" };
   BooleanProperty m_useGenieWeights      { this, "UseGenieWeights", false, "Flag to weight events according to Genie luminosity" };
   IntegerProperty m_flukaCollisions      { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." };
@@ -112,6 +123,7 @@ private:
   DoubleProperty  m_minMomentum          { this, "MinMomentum", 50000.0, "Write out all truth particles with a momentum larger MinMomentum"};
 
   DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 25000.0, "Blind events with Ecal energy above threshold (in MeV)"}; 
+  DoubleProperty m_caloMC_FudgeFactor    {this, "CaloFudgeFactorMC", 1.088, "Correct the MC energy calibration by this fudge factor"};
 
   double m_baseEventCrossSection {1.0};
   const double kfemtoBarnsPerMilliBarn {1.0e12};
@@ -151,6 +163,8 @@ private:
   mutable float m_wave_baseline_rms[15];
   mutable unsigned int m_wave_status[15];
  
+  mutable unsigned int m_scintHit;
+ 
   mutable float m_calibrated_nMIP[15];
   mutable float m_calibrated_E_dep[15];
   mutable float m_calibrated_E_EM[15];
@@ -159,6 +173,8 @@ private:
   mutable float m_calo_total_E_dep;
   mutable float m_calo_total_E_EM;
 
+  mutable float m_calo_total_E_EM_fudged;
+
   mutable float m_preshower_total_nMIP;
   mutable float m_preshower_total_E_dep;
 
@@ -170,9 +186,6 @@ private:
   mutable float m_Preshower12_Edep;
   mutable float m_Preshower13_Edep;
 
-  mutable float m_MIP_sim_Edep_calo;
-  mutable float m_MIP_sim_Edep_preshower;
-
   mutable float m_clock_phase;
 
   mutable unsigned int m_station0Clusters;
@@ -247,6 +260,10 @@ private:
   mutable std::vector<double> m_yCalo;
   mutable std::vector<double> m_thetaxCalo;
   mutable std::vector<double> m_thetayCalo;
+  mutable std::vector<double> m_x_atMaxRadius;
+  mutable std::vector<double> m_y_atMaxRadius;
+  mutable std::vector<double> m_z_atMaxRadius;
+  mutable std::vector<double> m_r_atMaxRadius;
 
   mutable std::vector<int> m_t_pdg; // pdg code of the truth matched particle
   mutable std::vector<int> m_t_barcode; // barcode of the truth matched particle
@@ -324,9 +341,6 @@ private:
 
   mutable std::vector<int> m_truth_pdg; // pdg of first 10 truth particles 
 
-
-
-
   mutable double m_truthLeptonMomentum;
   mutable int    m_truthBarcode;
   mutable int    m_truthPdg;
@@ -334,6 +348,8 @@ private:
 
   mutable int    m_eventsPassed = 0;
 
+  mutable std::array<double, 4> m_tracking_station_positions {}; // get average position of each tracking station for track extrapolation
+
 };
 
 inline const ServiceHandle <ITHistSvc> &NtupleDumperAlg::histSvc() const {