RG-L vertex fit

dc/java/dc-alignment/pom.xml

@@ -3,7 +3,7 @@
     <modelVersion>4.0.0</modelVersion>
     <groupId>org.clas.detector</groupId>
     <artifactId>dc-alignment</artifactId>
-    <version>2.3</version>
+    <version>2.4</version>
     <packaging>jar</packaging>
 
     <properties>

dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Alignment.java

@@ -421,6 +421,12 @@ public EmbeddedCanvasTabbed analyzeFits() {
         for(EmbeddedPad pad : canvas.getCanvas("LR residuals vs. theta").getCanvasPads())
             pad.getAxisX().setRange(-2000, 2000);
 
+        canvas.addCanvas("vertex");
+        canvas.getCanvas("vertex").draw(this.getVertexGraph("fit",null));
+        canvas.getCanvas().setFont(fontName);
+        for(EmbeddedPad pad : canvas.getCanvas("vertex").getCanvasPads())
+            pad.getAxisY().setRange(-2, 2);
+
         canvas.addCanvas("residual mean and sigma");
         canvas.getCanvas("residual mean and sigma").draw(this.getSectorHistograms("fit",null, 1));
         canvas.getCanvas("residual mean and sigma").draw(this.getSectorHistograms("time",null, 4));
@@ -770,6 +776,49 @@ private DataGroup getAngularGraph(String parameter, Table alignment) {
         return residuals;        
     }
 
+    private DataGroup getVertexGraph(String parameter, Table alignment) {
+
+        Map<Integer,List<GraphErrors>> graphs = new LinkedHashMap<>();
+        for (int i = 0; i < Constants.NTARGET; i++) {
+            int il = i==0 ? i : Constants.NLAYER+i;
+            for(int it=1; it<thetaBins.length; it++) {
+                GraphErrors gr_fit = new GraphErrors("gr_fit_layer " + il + "_theta" + it);
+                for(int is=0; is<Constants.NSECTOR; is++ ) {
+                    int sector = is+1;
+                    for(int ip=1; ip<phiBins.length; ip++) {
+                        double phi = is*60 + phiBins[ip].getMean();
+                        if(phi<0) phi += 360;
+                        double shiftRes = histos.get("nominal").getParValues(parameter, sector, it, ip)[il];
+                        double errorRes = 0.0;
+                        if(!(parameter.equals("time") || parameter.equals("LR"))) {
+                            shiftRes -= this.getFittedResidual(alignment, sector, it, ip)[il];
+                            errorRes = Math.sqrt(Math.pow(histos.get("nominal").getParErrors(parameter,sector, it, ip)[il], 2)
+                                                   + 0*Math.pow(this.getFittedResidualError(alignment, sector, it, ip)[il], 2));
+                        }
+                        if(Constants.MEASWEIGHTS[is][it][ip][il]>0 || parameter.equals("time") || parameter.equals("LR"))
+                            gr_fit.addPoint(phi, shiftRes/Constants.SCALE, 0.0, errorRes/Constants.SCALE);
+                    }
+                }               
+                gr_fit.setTitle("Layer " + (il+1));
+                gr_fit.setTitleX("#phi (deg)");
+                gr_fit.setTitleY("#Deltaz (cm)");
+                gr_fit.setMarkerColor(this.markerColor[it-1]);
+                gr_fit.setMarkerSize(this.markerSize);
+                if(gr_fit.getDataSize(0)>0) {
+                    if(!graphs.containsKey(i))
+                        graphs.put(i, new ArrayList<>());
+                    graphs.get(i).add(gr_fit);
+                }                    
+            }
+        }
+        DataGroup dg = new DataGroup(1,graphs.size());
+        for(int key : graphs.keySet()) {
+            for(GraphErrors gr : graphs.get(key))
+                dg.addDataSet(gr, key);
+        }
+        return dg;        
+    }
+
     private DataGroup getSectorHistograms(String parameter, Table alignment, int icol) {
 
         DataGroup residuals = new DataGroup(6,2);

dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Histo.java

@@ -940,6 +940,9 @@ private static boolean fitVertex(int mode, H1F histo) {
             case 10:    
                 Histo.fitRGEVertex(histo);
                 break;
+            case 11:    
+                Histo.fitRGLVertex(histo);
+                break;
             default:
                 if(histo.getFunction()!=null) {
                     if(histo.getFunction().getChiSquare()>0)
@@ -1498,8 +1501,6 @@ public static void fitRGKVertex(H1F histo) {
     }
 
     /**
-=======
->>>>>>> master
      * 4-peaks vertex fitting function
      * Peaks correspond to: target windows and scattering chamber exit window
      * Initialized according to:
@@ -1624,6 +1625,58 @@ public static void fitRGEVertex(H1F histo) {
 //        if(f1_vtx.getParameter(6)<f1_vtx.getParameter(0)/4) f1_vtx.setParameter(6, 0);
     }
 
+    /**
+     * 2-peaks vertex fitting function
+     * Peaks correspond to: target windows and scattering chamber exit window
+     * Initialized according to:
+     * - chosen target length (TARGETLENGTH), 
+     * - target exit window position (TARGETPOS)     
+     * Includes a wide Gaussian to account for target gas 
+     * @param histo
+     */
+    public static void fitRGLVertex(H1F histo) {
+        int nbin = histo.getData().length;
+        double dx = histo.getDataX(1)-histo.getDataX(0);
+        //find downstream window
+        int ibin0 = Histo.getMaximumBinBetween(histo, histo.getDataX(0), (Constants.TARGETPOS-Constants.TARGETLENGTH/2));
+        //check if the found maximum is the first or second peak, ibin is tentative upstream window
+        int ibin1 = Math.max(0, ibin0 - (int)(Constants.TARGETLENGTH/dx));
+        int ibin2 = Math.min(nbin-1, ibin0 + (int)(Constants.TARGETLENGTH/dx));
+        if(histo.getBinContent(ibin1)<histo.getBinContent(ibin2)) {
+            ibin1 = ibin0;
+            ibin0 = ibin2;
+        }
+
+        double mean  = histo.getDataX(ibin0);
+        double amp   = histo.getBinContent(ibin0);
+        double sigma = 0.5;
+        double bg = histo.getBinContent((ibin1+ibin0)/2);
+        String function = "[ampU]*gaus(x,[exw]-[tl],[sigmaU])+"
+                        + "[ampD]*gaus(x,[exw],[sigmaD])+"
+                        + "[bg]*landau(x,[bgmean],[bgsigma])+"
+                        + "[p0]+[p1]*x+[p2]*x*x";
+        F1D f1_vtx   = new F1D("f"+histo.getName(), function, -10, 10);
+        f1_vtx.setLineColor(2);
+        f1_vtx.setLineWidth(2);
+        f1_vtx.setOptStat("1111111111111111");
+        f1_vtx.setParameter(0, amp);
+        f1_vtx.setParameter(1, mean);
+        f1_vtx.setParameter(2, Constants.TARGETLENGTH);
+        f1_vtx.setParLimits(2, Constants.TARGETLENGTH*0.9, Constants.TARGETLENGTH*1.1);
+        f1_vtx.setParameter(3, sigma);
+        f1_vtx.setParameter(4, amp);                
+        f1_vtx.setParameter(5, sigma);
+        f1_vtx.setParameter(6, bg);
+        f1_vtx.setParameter(7, mean-Constants.TARGETLENGTH*0.5);
+        f1_vtx.setParLimits(7, mean-Constants.TARGETLENGTH*0.9,mean-Constants.TARGETLENGTH*0.1);
+        f1_vtx.setParameter(8, Constants.TARGETLENGTH*0.25);
+        f1_vtx.setRange(Math.max(mean-Constants.TARGETLENGTH-8*sigma,histo.getDataX(0)),
+                        Math.min(mean+9*sigma,histo.getDataX(nbin-1)));
+//        histo.setFunction(f1_vtx);
+        DataFitter.fit(f1_vtx, histo, "Q"); //No options uses error for sigma
+//        if(f1_vtx.getParameter(6)<f1_vtx.getParameter(0)/4) f1_vtx.setParameter(6, 0);
+    }
+
     //This was a previous version of fitting the z vertex peaks
     public static void fitVertex(H1F histo) {
         double mean  = histo.getDataX(histo.getMaximumBin());