Added possibility to scale track momenta by fixed amount (i.e. correc…

…ting bField)
JeffersonLab · Apr 25, 2024 · ca4d7f7 · ca4d7f7
1 parent cb4c776
commit ca4d7f7
Show file tree

Hide file tree

Showing 4 changed files with 67 additions and 5 deletions.
diff --git a/processors/include/NewVertexAnaProcessor.h b/processors/include/NewVertexAnaProcessor.h
@@ -152,6 +152,8 @@ class NewVertexAnaProcessor : public Processor {
         json v0proj_fits_;//!< json object v0proj
         double eleTrackTimeBias_ = 0.0;
         double posTrackTimeBias_ = 0.0;
+
+        double bFieldScaleFactor_ = -1;
         int current_run_number_{-999}; //!< track current run number
 };
 

diff --git a/processors/src/NewVertexAnaProcessor.cxx b/processors/src/NewVertexAnaProcessor.cxx
@@ -53,6 +53,9 @@ void NewVertexAnaProcessor::configure(const ParameterSet& parameters) {
         eleTrackTimeBias_ = parameters.getDouble("eleTrackTimeBias",eleTrackTimeBias_);
         posTrackTimeBias_ = parameters.getDouble("posTrackTimeBias",posTrackTimeBias_);
 
+	//bField scale factor (ONLY to correct for mistakes in ntuplizing). If < 0  is not used
+	bFieldScaleFactor_ = parameters.getDouble("bFieldScaleFactor",-1);
+
     }
     catch (std::runtime_error& error)
     {
@@ -355,21 +358,42 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
         if (debug_) std::cout << "got parts" << std::endl;
         Track ele_trk = ele->getTrack();
         Track pos_trk = pos->getTrack();
+
+	if (debug_) {
+	  std::cout<<"Check Ele/Pos Track momenta"<<std::endl;
+	  std::cout<<ele_trk.getP()<<" "<<pos_trk.getP()<<std::endl;
+	  std::cout<<ele_trk.getOmega()<<" "<<pos_trk.getOmega()<<std::endl;
+	}
 
-        //Beam Position Corrections
+	// Beam Position Corrections
         ele_trk.applyCorrection("z0", beamPosCorrections_.at(1));
         pos_trk.applyCorrection("z0", beamPosCorrections_.at(1));
-        //Track Time Corrections
+        // Track Time Corrections
         ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
         pos_trk.applyCorrection("track_time", posTrackTimeBias_);
 
-	//Correct for the momentum bias
+	// Correct for the momentum bias
 
 	if (biasingTool_) {
+
+	  // Correct for wrong track momentum - Bug Fix
+	  // In case there was mis-configuration during reco/hpstr-ntuple step, correct
+	  // the momentum magnitude here using the right bField for the data taking year
+
+	  if (bFieldScaleFactor_ > 0) {
+	    biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_);
+	    biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_);
+	  }
+
 	  biasingTool_->updateWithBiasP(ele_trk);
           biasingTool_->updateWithBiasP(pos_trk);
 	}
-
+
+	if (debug_) {
+	  std::cout<<"Corrected Ele/Pos Track momenta"<<std::endl;
+	  std::cout<<ele_trk.getP()<<" "<<pos_trk.getP()<<std::endl;
+	}
+
         double invm_smear = 1.;
         if (smearingTool_) {
           double unsmeared_prod = ele_trk.getP()*pos_trk.getP();
@@ -620,6 +644,7 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
     //TODO Clean this up => Cuts should be implemented in each region?
     //TODO Bring the preselection out of this stupid loop
 
+
 
     if (debug_) std::cout << "start regions" << std::endl;
     //TODO add yields. => Quite terrible way to loop.
@@ -671,8 +696,15 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
             //Track Time Corrections
             ele_trk.applyCorrection("track_time",eleTrackTimeBias_);
             pos_trk.applyCorrection("track_time", posTrackTimeBias_);
-
+	    
 	    if (biasingTool_) {
+
+	      //Correct the wrong Bfield first
+	      if (bFieldScaleFactor_ > 0) {
+		biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_);
+		biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_);
+	      }
+
 	      biasingTool_->updateWithBiasP(ele_trk);
 	      biasingTool_->updateWithBiasP(pos_trk);
 	    }
@@ -1076,6 +1108,17 @@ bool NewVertexAnaProcessor::process(IEvent* ievent) {
 	    // Track Momentum bias
 
 	    if (biasingTool_) {
+
+	      // Correct for wrong track momentum - Bug Fix
+	      // In case there was mis-configuration during reco/hpstr-ntuple step, correct
+	      // the momentum magnitude here using the right bField for the data taking year
+
+	      if (bFieldScaleFactor_ > 0) {
+		biasingTool_->updateWithBiasP(ele_trk,bFieldScaleFactor_);
+		biasingTool_->updateWithBiasP(pos_trk,bFieldScaleFactor_);
+	      }
+
+
 	      biasingTool_->updateWithBiasP(ele_trk);
 	      biasingTool_->updateWithBiasP(pos_trk);
 	    }

diff --git a/utils/include/TrackBiasingTool.h b/utils/include/TrackBiasingTool.h
@@ -24,6 +24,11 @@ class TrackBiasingTool {
 		   const std::string& tracks = "KalmanFullTracks");
 
   double biasTrackP(const Track& track);
+
+  //Update the track P with a specific scale Factor
+  void updateWithBiasP(Track& trk, double scaleFactor);
+
+  //Update the track P with scale Factors according to the internal calibration plots
   void updateWithBiasP(Track& trk);
 
 private:

diff --git a/utils/src/TrackBiasingTool.cxx b/utils/src/TrackBiasingTool.cxx
@@ -57,3 +57,15 @@ void TrackBiasingTool::updateWithBiasP(Track& trk) {
   trk.setMomentum(momentum);
 
 }
+
+// This will recompute the track momentum from curvature and store it in the track
+void TrackBiasingTool::updateWithBiasP(Track& trk, double scaleFactor) {
+
+  std::vector<double> momentum = trk.getMomentum();
+
+  for (double& coordinate : momentum)
+    coordinate *= scaleFactor ;
+
+  trk.setMomentum(momentum);    
+
+}