Modify forward electron Vz to add target window (#464)

Author: zhiwanxu

src/main/java/org/jlab/clas/timeline/analysis/forward/forward_Tracking_EleVz.groovy

@@ -9,6 +9,7 @@ class forward_Tracking_EleVz {
 
 def data = new ConcurrentHashMap()
 def data_2nd_peak = new ConcurrentHashMap()
+def data_upstream = new ConcurrentHashMap()
 
 def is_RGD_LD2 = { run ->
   return (RunDependentCut.runIsInRange(run, 18305, 18313, true) ||
@@ -32,6 +33,11 @@ def processRun(dir, run) {
   def sigmalist_2nd_peak = []
   def chi2list_2nd_peak = []
 
+  def funclist_upstream = []
+  def meanlist_upstream = []
+  def sigmalist_upstream = []
+  def chi2list_upstream = []
+
   def histlist =   (0..<6).collect{
     def h1 = dir.getObject('/elec/H_trig_vz_mom_S'+(it+1)).projectionY()
     h1.setName("sec"+(it+1))
@@ -56,7 +62,12 @@ def processRun(dir, run) {
       }
     }
     else if (dataset == 'rgl') {
-      f1 = ForwardFitter.fitRGL(h1)
+      f1 = ForwardFitter.fitRGL(h1, 10, 30)
+      def f1_upstream = ForwardFitter.fitRGL(h1, -40, -20)
+      funclist_upstream.add(f1_upstream)
+      meanlist_upstream.add(f1_upstream.getParameter(1))
+      sigmalist_upstream.add(f1_upstream.getParameter(2).abs())
+      chi2list_upstream.add(f1_upstream.getChiSquare())
     }
     else {
       f1 = ForwardFitter.fit(h1)
@@ -80,6 +91,9 @@ def processRun(dir, run) {
   if (funclist_2nd_peak.size() > 0) {
     data_2nd_peak[run] = [run:run, hlist:histlist, flist:funclist_2nd_peak, mean:meanlist_2nd_peak, sigma:sigmalist_2nd_peak, clist:chi2list_2nd_peak]
   }
+  if (funclist_upstream.size() > 0) {
+    data_upstream[run] = [run:run, hlist:histlist, flist:funclist_upstream, mean:meanlist_upstream, sigma:sigmalist_upstream, clist:chi2list_upstream]
+  }
 }
 
 
@@ -118,6 +132,30 @@ def write() {
     if (data_2nd_peak.size() != 0) {
       out.addDataSet(grtl_2nd_peak)
     }
+
+    if (data_upstream.size() != 0) {
+      def grtl_upstream = new GraphErrors('sec'+(sec+1)+'_upstream')
+      grtl_upstream.setTitle("VZ (upstream peak value) for electrons per sector")
+      grtl_upstream.setTitleY("VZ (upstream peak value) for electrons per sector (cm)")
+      grtl_upstream.setTitleX("run number")
+
+      def grtl_window = new GraphErrors('sec'+(sec+1)+'_window')
+      grtl_window.setTitle("VZ target window length (downstream - upstream peak distance) per sector")
+      grtl_window.setTitleY("VZ target window length (cm)")
+      grtl_window.setTitleX("run number")
+
+      data_upstream.sort{it.key}.each{run,it->
+        out.cd('/'+it.run)
+        out.addDataSet(it.flist[sec])
+        grtl_upstream.addPoint(it.run, it.mean[sec], 0, 0)
+        def window_length = (data[run].mean[sec] - it.mean[sec]).abs()
+        grtl_window.addPoint(it.run, window_length, 0, 0)
+      }
+
+      out.cd('/timelines')
+      out.addDataSet(grtl_upstream)
+      out.addDataSet(grtl_window)
+    }
   }
 
   out.writeFile('forward_electron_VZ.hipo')

src/main/java/org/jlab/clas/timeline/analysis/forward/forward_Tracking_NegVz.groovy

@@ -55,7 +55,7 @@ def processRun(dir, run) {
       }
     }
     else if (dataset == 'rgl') {
-      f1 = ForwardFitter.fitRGL(h1)
+      f1 = ForwardFitter.fitRGL(h1, 10, 30)
     }
     else {
       f1 = ForwardFitter.fit(h1)

src/main/java/org/jlab/clas/timeline/analysis/forward/forward_Tracking_PosVz.groovy

@@ -55,7 +55,7 @@ def processRun(dir, run) {
       }
     }
     else if (dataset == 'rgl') {
-      f1 = ForwardFitter.fitRGL(h1)
+      f1 = ForwardFitter.fitRGL(h1, 10, 30)
     }
     else {
       f1 = ForwardFitter.fit(h1)

src/main/java/org/jlab/clas/timeline/fitter/ForwardFitter.groovy

@@ -52,30 +52,35 @@ class ForwardFitter{
   }
 
 
-  static F1D fitRGL(H1F h1) {
-    def f1 = new F1D("fit:"+h1.getName(), "[amp]*gaus(x,[mean],[sigma])", 15, 25);
-    double hAmp  = h1.getBinContent(h1.getMaximumBin());
-    double hMean = h1.getAxis().getBinCenter(h1.getMaximumBin());
-    double hRMS  = h1.getRMS(); //cm
-    f1.setRange(hMean-1.0, hMean+1.0);
+  static F1D fitRGL(H1F h1, double fit_left, double fit_right) {
+    def f1 = new F1D("fit:"+h1.getName(), "[amp]*gaus(x,[mean],[sigma])", fit_left, fit_right);
+
+    // find local max within [fit_left, fit_right] only, avoiding the other peak
+    int binLeft  = Math.max(0, h1.getAxis().getBin(fit_left))
+    int binRight = Math.min(h1.getAxis().getNBins() - 1, h1.getAxis().getBin(fit_right))
+    int maxBin = binLeft
+    for (int i = binLeft + 1; i <= binRight; i++) {
+      if (h1.getBinContent(i) > h1.getBinContent(maxBin)) maxBin = i
+    }
+    double hAmp  = h1.getBinContent(maxBin)
+    double hMean = h1.getAxis().getBinCenter(maxBin)
+
+    f1.setRange(hMean - 2.0, hMean + 2.0);
     f1.setParameter(0, hAmp);
     f1.setParameter(1, hMean);
-    //f1.setParameter(2, hRMS);
-    f1.setParameter(2, 0.3);
+    f1.setParameter(2, 1);
     MoreFitter.fit(f1,h1,"LQ");
 
     def makefit = {func->
       hMean = func.getParameter(1)
-      hRMS = func.getParameter(2).abs()
-      func.setRange(hMean-1.0,hMean+1.0)
+      func.setRange(hMean - 2.0, hMean + 2.0)
       MoreFitter.fit(func,h1,"LQ")
       return [func.getChiSquare(), (0..<func.getNPars()).collect{func.getParameter(it)}]
     }
 
     def fits1 = (0..20).collect{makefit(f1)}
     def bestfit = fits1.sort()[0]
     f1.setParameters(*bestfit[1])
-    //makefit(f1)
 
     return f1
   }

src/main/java/org/jlab/clas/timeline/histograms/GeneralMon.java

@@ -1250,7 +1250,7 @@ public GeneralMon(int reqrunNum, String reqOutputDir, float reqEB, boolean reqTi
       H_trig_theta_phi_S[s].setTitle(String.format("e sect %d",s+1));
       H_trig_theta_phi_S[s].setTitleX("#phi (^o)");
       H_trig_theta_phi_S[s].setTitleY("#theta (^o)");
-      H_trig_vz_mom_S[s] = new H2F(String.format("H_trig_vz_mom_S%d",s+1),String.format("H_trig_vz_mom_S%d",s+1),100,0,EB,100,-25,50);
+      H_trig_vz_mom_S[s] = new H2F(String.format("H_trig_vz_mom_S%d",s+1),String.format("H_trig_vz_mom_S%d",s+1),100,0,EB,200,-50,50);
       H_trig_vz_mom_S[s].setTitle(String.format("e sect %d",s+1));
       H_trig_vz_mom_S[s].setTitleX("p (GeV)");
       H_trig_vz_mom_S[s].setTitleY("vz (cm)");