diff --git a/dc/java/dc-alignment/pom.xml b/dc/java/dc-alignment/pom.xml index 8c9ebaa..6835c00 100644 --- a/dc/java/dc-alignment/pom.xml +++ b/dc/java/dc-alignment/pom.xml @@ -3,7 +3,7 @@ 4.0.0 org.clas.detector dc-alignment - 2.4 + 2.5 jar @@ -57,7 +57,7 @@ jhep-maven - https://clasweb.jlab.org/jhep/maven + https://clasweb.jlab.org/.jhep/maven diff --git a/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Alignment.java b/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Alignment.java index 16c3b78..eeb7e5a 100644 --- a/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Alignment.java +++ b/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Alignment.java @@ -201,6 +201,14 @@ else if(Math.abs(vertex)==6) Constants.initTargetPars(Constants.RGCSUMMER2022); else if(Math.abs(vertex)==7) Constants.initTargetPars(Constants.RGDFALL2023); + else if(Math.abs(vertex)==8) + Constants.initTargetPars(Constants.RGASPRING2018); + else if(Math.abs(vertex)==9) + Constants.initTargetPars(Constants.RGKSPRING2024); + else if(Math.abs(vertex)==10) + Constants.initTargetPars(Constants.RGESPRING2024); + else if(Math.abs(vertex)==11) + Constants.initTargetPars(Constants.RGLSPRING2025); } if(vertex<0) vertex=0; } @@ -231,13 +239,13 @@ public void addHistoSet(String name, Histo histo) { this.histos.put(name, histo); } - public void analyzeHistos(int resFit, int vertexFit, String vertexPar, boolean test) { + public void analyzeHistos(int resFit, int vertexFit, String vertexPar, String test) { this.printConfig(resFit, "resFit", ""); this.printConfig(vertexFit, "vertexFit", ""); this.initMeasurementWeights(1.0); this.initVertexPar(vertexFit, vertexPar); for(String key : histos.keySet()) { - if(test && !key.equals("nominal")) continue; + if(histos.containsKey(test) && !key.equals(test)) continue; LOGGER.log(LEVEL,"\nAnalyzing histos for variation " + key); histos.get(key).analyzeHisto(resFit, vertexFit); for(int i=0; i0) { if(!graphs.containsKey(i)) @@ -865,7 +886,7 @@ private DataGroup getRegionHistograms(String parameter, Table alignment, int ico residuals.addDataSet(hi_res, ir*Constants.NSECTOR + is); for(int it=1; it0 && pars.length<=DEFAULT.length) { TARGETPOS = pars[0]; TARGETCENTER = pars[0]; + System.out.println("[CONFIG] target parameters set to:"); + System.out.println(" - TARGETPOS = " + TARGETPOS); if(pars.length>1) { TARGETLENGTH = pars[1]; TARGETCENTER = pars[0]-pars[1]/2; + System.out.println(" - TARGETLENGTH = " + TARGETLENGTH); } - if(pars.length>2) + if(pars.length>2) { WINDOWDIST = pars[2]; - if(pars.length>3) + System.out.println(" - WINDOWDIST = " + WINDOWDIST); + } + if(pars.length>3) { SCEXIT = pars[3]; + System.out.println(" - SCEXIT = " + SCEXIT); + } } else { System.out.println("[WARNING] wrong number of target parameters. Number is " + pars.length + " instead of [1:4]"); } - System.out.println("[CONFIG] target parameters set to:"); - System.out.println(" - TARGETPOS = " + TARGETPOS); - System.out.println(" - TARGETLENGTH = " + TARGETLENGTH); - System.out.println(" - WINDOWDIST = " + WINDOWDIST); - System.out.println(" - SCEXIT = " + SCEXIT); } } diff --git a/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Histo.java b/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Histo.java index 7516f33..ccf4425 100644 --- a/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Histo.java +++ b/dc/java/dc-alignment/src/main/java/org/clas/dc/alignment/Histo.java @@ -42,21 +42,23 @@ public class Histo { private DataGroup offset = null; private DataGroup calib = null; private DataGroup alpha = null; + private DataGroup track = null; private DataGroup doca = null; - private DataGroup[] wires = null; + private DataGroup[][] wires = null; private DataGroup[][][] residuals = null; // indices are theta bin, phi bin and sector, datagroup is 6x6 and contains layers private DataGroup[][][] time = null; // indices are theta bin, phi bin and sector, datagroup is 6x6 and contains layers private DataGroup[][][] leftright = null; // indices are theta bin, phi bin and sector, datagroup is 6x6 and contains layers private DataGroup[][] vertex = null; // indices are theta bin, phi bin and sector, datagroup is 6x6 and contains sectors - private double[][][][] zeroes = null; - private double[][][][] timeValues = null; - private double[][][][] timeSigmas = null; - private double[][][][] lrValues = null; - private double[][][][] parValues = null; - private double[][][][] parErrors = null; - private double[][][][] parSigmas = null; - private double[][] beamOffset= {{0, 0}, {0, 0}}; + private double[][][][] zeroes = null; + private double[][][][] timeValues = null; + private double[][][][] timeSigmas = null; + private double[][][][] lrValues = null; + private double[][][][] parValues = null; + private double[][][][] parErrors = null; + private double[][][][] parSigmas = null; + private boolean[][][][] parStatus = null; + private double[][] beamOffset= {{0, 0}, {0, 0}}; private Bin[] thetaBins = null; private Bin[] phiBins = null; @@ -68,7 +70,7 @@ public class Histo { private List shiftedFiles = null; private SchemaFactory schema = null; - private boolean tres = false; + private boolean extra = false; private boolean shift = false; private static final Logger LOGGER = Logger.getLogger(Constants.LOGGERNAME); @@ -77,17 +79,19 @@ public Histo(String name, List files, Bin[] thetabins, Bin[] phibins, do this.name = name; this.thetaBins = thetabins; this.phiBins = phibins; + this.minVtx = vertexrange[0]; + this.maxVtx = vertexrange[1]; this.nominalFiles = files; this.createHistos(optstats); } - public Histo(String name, List files, Bin[] thetabins, Bin[] phibins, boolean time, double[] vertexrange, String optstats) { + public Histo(String name, List files, Bin[] thetabins, Bin[] phibins, boolean extraplots, double[] vertexrange, String optstats) { this.name = name; this.thetaBins = thetabins; this.phiBins = phibins; this.minVtx = vertexrange[0]; this.maxVtx = vertexrange[1]; - this.tres = time; + this.extra = extraplots; this.nominalFiles = files; this.createHistos(optstats); } @@ -110,12 +114,12 @@ public Histo(String name, boolean shift, Bin[] thetabins, Bin[] phibins, double[ this.createHistos(optstats); } - public Histo(String name, boolean shift, Bin[] thetabins, Bin[] phibins, boolean time, double[] vertexrange, String optstats) { + public Histo(String name, boolean shift, Bin[] thetabins, Bin[] phibins, boolean extraplots, double[] vertexrange, String optstats) { this.name = name; this.thetaBins = thetabins; this.phiBins = phibins; this.shift = shift; - this.tres = time; + this.extra = extraplots; this.createHistos(optstats); } @@ -125,8 +129,7 @@ private void createHistos(String optStats) { + (phiBins.length-1) + " phi bins " + "for variation " + name); this.residuals = new DataGroup[nSector][thetaBins.length][phiBins.length]; - if(tres) { - this.wires = new DataGroup[nSector]; + if(extra) { this.time = new DataGroup[nSector][thetaBins.length][phiBins.length]; this.leftright = new DataGroup[nSector][thetaBins.length][phiBins.length]; } @@ -135,6 +138,7 @@ private void createHistos(String optStats) { this.parValues = new double[nSector][thetaBins.length][phiBins.length][nLayer+nTarget]; this.parErrors = new double[nSector][thetaBins.length][phiBins.length][nLayer+nTarget]; this.parSigmas = new double[nSector][thetaBins.length][phiBins.length][nLayer+nTarget]; + this.parStatus = new boolean[nSector][thetaBins.length][phiBins.length][nLayer+nTarget]; this.timeValues = new double[nSector][thetaBins.length][phiBins.length][nLayer+nTarget]; this.timeSigmas = new double[nSector][thetaBins.length][phiBins.length][nLayer+nTarget]; this.lrValues = new double[nSector][thetaBins.length][phiBins.length][nLayer+nTarget]; @@ -152,15 +156,54 @@ private void createHistos(String optStats) { maxVtx = Constants.VDFMAX; } - this.calib = new DataGroup(nSector, nSLayer); - this.alpha = new DataGroup(nSector, nSLayer); - this.doca = new DataGroup(nSector, nSLayer); + if(extra) { + this.calib = new DataGroup(nSector, nSLayer); + this.alpha = new DataGroup(nSector, nSLayer); + this.track = new DataGroup(nSector, nSLayer); + this.doca = new DataGroup(nSector, nSLayer); + this.wires = new DataGroup[2][nSector]; + for(int ires=0; ires<2; ires++) { + for(int is=0; is0) this.leftright[sector-1][it][ip].getH1F("hi-rL" + hit.layer).fill(hit.time); @@ -611,10 +629,12 @@ public void analyzeHisto(int fit, int vertexFit) { Histo.fitVertex(vertexFit,electron.getH1F("hi-vtx")); for(int is=0; is24) this.parErrors[is][it][ip][l] *= 2; @@ -656,10 +678,12 @@ public void analyzeHisto(int fit, int vertexFit) { } H1F hvtx = vertex[it][ip].getH1F("hi-S"+s); double dx = hvtx.getDataX(1)-hvtx.getDataX(0); + System.out.print(String.format("\tsector=%1d theta bin=%1d phi bin=%1d layer=%2d",s,it,ip,0)); if(Histo.fitVertex(vertexFit, hvtx)) { this.parValues[is][it][ip][0] = hvtx.getFunction().getParameter(1)*Constants.SCALE; this.parErrors[is][it][ip][0] = hvtx.getFunction().parameter(1).error()*Constants.SCALE; this.parSigmas[is][it][ip][0] = hvtx.getFunction().getParameter(2)*Constants.SCALE; + this.parStatus[is][it][ip][0] = true; if(!shift) { this.parValues[is][it][ip][0] -= Constants.TARGETPOS*Constants.SCALE; this.parErrors[is][it][ip][0] = Math.max(this.parErrors[is][it][ip][0], Constants.SCALE*dx/2); @@ -675,6 +699,7 @@ public void analyzeHisto(int fit, int vertexFit) { this.parValues[is][it][ip][nLayer+nTarget-1] = (hvtx.getFunction().getParameter(iscw)-Constants.SCEXIT)*Constants.SCALE; this.parErrors[is][it][ip][nLayer+nTarget-1] = Math.max(hvtx.getFunction().parameter(iscw).error()*Constants.SCALE, Constants.SCALE*dx); this.parSigmas[is][it][ip][nLayer+nTarget-1] = hvtx.getFunction().getParameter(2)*Constants.SCALE; + this.parStatus[is][it][ip][nLayer+nTarget-1] = true; } else { Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-1]=0; @@ -683,6 +708,7 @@ public void analyzeHisto(int fit, int vertexFit) { this.parValues[is][it][ip][nLayer+nTarget-2] = (hvtx.getFunction().getParameter(itl)-Constants.TARGETLENGTH)*Constants.SCALE; this.parErrors[is][it][ip][nLayer+nTarget-2] = Math.max(hvtx.getFunction().parameter(itl).error()*Constants.SCALE, Constants.SCALE*dx); this.parSigmas[is][it][ip][nLayer+nTarget-2] = hvtx.getFunction().getParameter(2)*Constants.SCALE; + this.parStatus[is][it][ip][nLayer+nTarget-2] = true; } else { Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-2]=0; @@ -694,6 +720,7 @@ public void analyzeHisto(int fit, int vertexFit) { Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-2]=0; Constants.MEASWEIGHTS[is][it][ip][nLayer+nTarget-1]=0; } + System.out.print("\r"); } } } @@ -747,6 +774,16 @@ public double[][] getBeamOffset() { return this.beamOffset; } + public boolean[] getParStatus(int sector, int itheta, int iphi) { + if(sector<1 || sector>6) + throw new IllegalArgumentException("Error: invalid sector="+sector); + if(itheta<0 || itheta>=thetaBins.length) + throw new IllegalArgumentException("Error: invalid theta bin="+itheta); + if(iphi<0 || iphi>phiBins.length) + throw new IllegalArgumentException("Error: invalid phi bin="+iphi); + return this.parStatus[sector-1][itheta][iphi]; + } + public double[] getParValues(int sector, int itheta, int iphi) { return this.getParValues("", sector, itheta, iphi); } @@ -835,18 +872,30 @@ public EmbeddedCanvasTabbed getElectronPlots() { } public EmbeddedCanvasTabbed plotHistos() { - EmbeddedCanvasTabbed canvas = new EmbeddedCanvasTabbed("Calibration", "Alpha", "Doca"); - canvas.getCanvas("Calibration").draw(calib); - canvas.getCanvas("Alpha").draw(alpha); - for(EmbeddedPad pad : canvas.getCanvas("Alpha").getCanvasPads()) - pad.getAxisZ().setLog(true); - canvas.getCanvas("Doca").draw(doca); - if(tres) { + EmbeddedCanvasTabbed canvas = null; + if(this.extra) { + canvas = new EmbeddedCanvasTabbed("Calibration", "Alpha", "Track", "Doca"); + canvas.getCanvas("Calibration").draw(calib); + canvas.getCanvas("Alpha").draw(alpha); + for(EmbeddedPad pad : canvas.getCanvas("Alpha").getCanvasPads()) + pad.getAxisZ().setLog(true); + canvas.getCanvas("Track").draw(track); + canvas.getCanvas("Doca").draw(doca); + for(int is=0; is0 ? histo.getBinContent((ibin1+ibin2)/2) : histo.getBinContent(0); + String function = "[ampD]*gaus(x,[exw],[sigmaD])+" + + "[ampDL]*landau(x,[exw],[sigmaD]*2.5)+" + + "[f0]+([p0]+[p1]*x+[p2]*x*x)"; F1D f1_vtx = new F1D("f"+histo.getName(), function, -10, 10); + f1_vtx.setOptStat("1111111111"); + int np = f1_vtx.getNPars(); + if(ibin1>0) { + f1_vtx = new F1D("f"+histo.getName(), function + + "*erf(x,[exw]-[tl]-0.9,-[sigmaUL])" + + "+[ampU]*gaus(x,[exw]-[tl],[sigmaU])" + + "+[ampU]*gaus(x,[exw]-[tl]-0.9,[sigmaU])", -10, 10); + f1_vtx.setOptStat("11111111111111"); + f1_vtx.setParameter(np+0, meanD-meanU-Constants.WINDOWDIST/2);//Constants.TARGETLENGTH); + f1_vtx.setParLimits(np+0, Constants.TARGETLENGTH*0.9, Constants.TARGETLENGTH*1.1); +// f1_vtx.setParameter(np+1, Constants.WINDOWDIST); +// f1_vtx.setParLimits(np+1, Constants.WINDOWDIST*0.99, Constants.WINDOWDIST*1.01); + f1_vtx.setParameter(np+1, sigma*5); + f1_vtx.setParStep(np+1, 0); + f1_vtx.setParameter(np+2, ampU/2); + f1_vtx.setParameter(np+3, sigma*2); +// f1_vtx.setParameter(np+4, ampU*0.1); + } 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); + f1_vtx.setParameter(0, ampD); + f1_vtx.setParameter(1, meanD); + f1_vtx.setParameter(2, sigma); + f1_vtx.setParameter(3, ampD*0.1); + f1_vtx.setParameter(4, 0); + f1_vtx.setParameter(5, bg/4); + f1_vtx.setRange(Math.max(meanU-12*sigma,histo.getDataX(0)), + Math.min(meanD+8*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)Constants.CHI2MAX) && ibin1>0) { + meanU = f1_vtx.getParameter(1)-f1_vtx.getParameter(np+1); + meanD = f1_vtx.getParameter(1); + sigma = f1_vtx.getParameter(np+2); + f1_vtx.setRange(Math.max(meanU-6*sigma,histo.getDataX(0)), + Math.min(meanD+8*sigma,histo.getDataX(nbin-1))); + DataFitter.fit(f1_vtx, histo, "Q"); //No options uses error for sigma + } } //This was a previous version of fitting the z vertex peaks @@ -1779,7 +1849,7 @@ public static int getMaximumBinBetween(H1F histo, double min, double max) { double x_val; double y_max_temp; double y_max = 0; - int max_bin_num = histo.getMaximumBin(); + int max_bin_num = 0; for (int i = 0; i < nbin; i++) { x_val_temp = histo.getAxis().getBinCenter(i); if (x_val_temp >= min && x_val_temp <= max) { @@ -1879,14 +1949,32 @@ private static double getRMSIDataSet(IDataSet data, double min, double max) { } return Math.sqrt(sum / (double) nEntries); } + + public static boolean rebin(H1F histo, int factor) { + + int nbin = histo.getData().length; + if(nbin histos = new LinkedHashMap<>(); + + + public Kinematics(double energy, String vertex) { + this.initGraphics(); + this.initVertexPar(vertex); + this.setEbeam(energy); + } + + private void initVertexPar(String pars) { + if(!pars.isEmpty()) { + double[] parValues = new double[pars.split(":").length]; + for(int i=0; i0) targetPos = parValues[0]; + if(parValues.length>1) targetLength = parValues[1]; + if(parValues.length>2) scWindow = parValues[2]+targetPos; + } + } + + public double getEbeam() { + return ebeam; + } + + public void setEbeam(double ebeam) { + this.ebeam = ebeam; + } + + public double getMaxQ2() { + double maxQ2 = 3; + if(ebeam>6) maxQ2 = 5; + return maxQ2; + } + + public double getMaxW() { + double maxW = ebeam*0.4+0.6; + if(ebeam>6.5) maxW = 4.5; + return maxW; + } + + private void initGraphics() { + GStyle.getAxisAttributesX().setTitleFontSize(18); + GStyle.getAxisAttributesX().setLabelFontSize(14); + GStyle.getAxisAttributesY().setTitleFontSize(18); + GStyle.getAxisAttributesY().setLabelFontSize(14); + GStyle.getAxisAttributesZ().setLabelFontSize(14); + GStyle.getAxisAttributesX().setLabelFontName(this.fontName); + GStyle.getAxisAttributesY().setLabelFontName(this.fontName); + GStyle.getAxisAttributesZ().setLabelFontName(this.fontName); + GStyle.getAxisAttributesX().setTitleFontName(this.fontName); + GStyle.getAxisAttributesY().setTitleFontName(this.fontName); + GStyle.getAxisAttributesZ().setTitleFontName(this.fontName); + GStyle.setGraphicsFrameLineWidth(1); + GStyle.getH1FAttributes().setLineWidth(2); + } + + public void createHistos() { + // General + H2F hi_q2w = new H2F("hi_q2w","hi_q2w",100, 0.6, this.getMaxW(), 100, 0.0, this.getMaxQ2()); + hi_q2w.setTitleX("W (GeV)"); + hi_q2w.setTitleY("Q2 (GeV2)"); + H1F hi_w = new H1F("hi_w","hi_w",500, 0.6, this.getMaxW()); + hi_w.setTitleX("W (GeV)"); + hi_w.setTitleY("Counts"); + H2F hi_w_phi = new H2F("hi_w_phi","hi_w_phi",100, -180.0, 180.0, 250, 0.6, this.getMaxW()); + hi_w_phi.setTitleX("#phi (deg)"); + hi_w_phi.setTitleY("W (GeV)"); + H2F hi_z_phi = new H2F("hi_z_phi","hi_z_phi",100, -180.0, 180.0, 100, 20, 30); + hi_z_phi.setTitleX("#phi (deg)"); + hi_z_phi.setTitleY("z (cm)"); + H2F hi_w_theta = new H2F("hi_w_theta","hi_w_theta",100, 5.0, 15.0, 250, 0.6, this.getMaxW()); + hi_w_theta.setTitleX("#theta (deg)"); + hi_w_theta.setTitleY("W (GeV)"); + H2F hi_el = new H2F("hi_el","hi_el",100, 0.5, ebeam+0.5, 100, 0.0, 35.0); + hi_el.setTitleX("p (GeV)"); + hi_el.setTitleY("#theta (deg)"); + hi_el.setTitle("Electron"); + F1D f1_el = new F1D("f1_el", "2*(180/3.14)*atan(sqrt(0.93832*([e0]-x)/2/[e0]/x))", ebeam*0.75, ebeam*0.99); + f1_el.setParameter(0, ebeam); + DataGroup dg_general = new DataGroup(3,2); + dg_general.addDataSet(hi_q2w, 0); + dg_general.addDataSet(hi_w_theta, 1); + dg_general.addDataSet(hi_el, 2); + dg_general.addDataSet(f1_el, 2); + dg_general.addDataSet(hi_w, 3); + dg_general.addDataSet(hi_w_phi, 4); + dg_general.addDataSet(hi_z_phi, 5); + histos.put("General", dg_general); + // Elastic + H1F hi_ela = new H1F("hi_w","hi_w",500, 0.6, 1.3); + hi_ela.setTitleX("W (GeV)"); + hi_ela.setTitleY("Counts"); + H2F hi_ela_phi = new H2F("hi_w_phi","hi_w_phi",100, -180.0, 180.0, 250, 0.6, 1.3); + hi_ela_phi.setTitleX("#phi (deg)"); + hi_ela_phi.setTitleY("W (GeV)"); + H2F hi_ela_theta = new H2F("hi_w_theta","hi_w_theta",100, 5.0, 15.0, 250, 0.6, 1.3); + hi_ela_theta.setTitleX("#theta (deg)"); + hi_ela_theta.setTitleY("W (GeV)"); + H2F hi_ela_p = new H2F("hi_w_p","hi_w_p",100, ebeam*0.8, ebeam, 100, 0.6, 1.3); + hi_ela_p.setTitleX("p (GeV)"); + hi_ela_p.setTitleY("W (GeV)"); + hi_ela_p.setTitle("Electron"); + DataGroup dg_elastic = new DataGroup(2,2); + dg_elastic.addDataSet(hi_ela, 0); + dg_elastic.addDataSet(hi_ela_p, 1); + dg_elastic.addDataSet(hi_ela_theta, 2); + dg_elastic.addDataSet(hi_ela_phi, 3); + histos.put("Elastic", dg_elastic); + // vertex + String[] names= ["pi+", "pi-", "el"]; + DataGroup dg_vertex = new DataGroup(3, 4); + for(int i=0; i1) rmin = -0.5; + H1F hi_mmass = new H1F("mxt" + i, "", 200, rmin, rmax); + hi_mmass.setTitleX(a2pimx[i-1]); + hi_mmass.setTitleY("Counts"); + // hi_mmass.setLineColor(col); + H1F hi_imass = new H1F("mit" + i, "", 200, 0.0, 3.0); + hi_imass.setTitleX("M(#pi#pi) (GeV)"); + hi_imass.setTitleY("Counts"); + // hi_imass.setLineColor(col); + dg_2pi.addDataSet(hi_mmass, i-1); + dg_2pi.addDataSet(hi_imass, i+2); + } + histos.put("2pi", dg_2pi); + // 1 pi + DataGroup dg_1pi = new DataGroup(2,2); + String[] atitle = ["M_X(ep#rarrow e'^#pi+X) (GeV)","M_^X2(ep#rarrow e'pX) (Ge^V2)"]; + for(int i=1; i<=2; i++) { + double rmin = 0.5; + double rmax = 2; + if(i>1) rmin = -0.5; + H2F hi_mw = new H2F("W" + i, "", 100, -180, 180, 100, rmin, rmax); + hi_mw.setTitleX("W (GeV)"); + hi_mw.setTitleY(atitle[i-1]); + // hi_w.setLineColor(col); + H1F hi_mmass = new H1F("mxt" + i, "", 400, rmin, rmax); + hi_mmass.setTitleX(atitle[i-1]); + hi_mmass.setTitleY("Counts"); + // hi_mmass.setLineColor(col); + dg_1pi.addDataSet(hi_mw, i-1); + dg_1pi.addDataSet(hi_mmass, i+1); + } + histos.put("1pi", dg_1pi); + // Proton + DataGroup dg_proton = new DataGroup(4,2); + H2F hi_pr_p_dp = new H2F("hi_p_dp", "", 100, 1, 4, 100, ebeam*0.75, ebeam*1.2); + hi_pr_p_dp.setTitleX("p (GeV)"); hi_pr_p_dp.setTitleY("Ebeam (GeV)"); + H2F hi_pr = new H2F("hi_pr","hi_pr",100, 0, ebeam/2-0.1, 100, 20.0, 80.0); + hi_pr.setTitleX("p (GeV)"); + hi_pr.setTitleY("#theta (deg)"); + hi_pr.setTitle("Proton"); + F1D f1_pr = new F1D("f1_pr", "(180/3.14)*acos(([e0]*[e0]+x*x-pow(([e0]+0.93832-sqrt(x*x+0.9382*0.9382)),2))/2/[e0]/x)", ebeam*0.08, ebeam*0.4); + f1_pr.setParameter(0, ebeam); + H2F hi_phi = new H2F("hi_phi", "hi_phi", 200, -180, 180, 200, -180, 180); + hi_phi.setTitleX("El #phi (deg)"); + hi_phi.setTitleY("Pr #phi (deg)"); + // hi_dphi.setOptStat("1110"); + H2F hi_dphi = new H2F("hi_dphi", "hi_dphi", 200, -180, 180, 200, -10, 10); + hi_dphi.setTitleX("Pr #phi (deg)"); + hi_dphi.setTitleY("#Delta#phi (deg)"); + H2F hi_dpr = new H2F("hi_dpr", "hi_dpr", 200, -180, 180, 200, -10, 10); + hi_dpr.setTitleX("Pr #phi (deg)"); + hi_dpr.setTitleY("#Delta p (GeV)"); + H2F hi_pr_phi_dphi = new H2F("hi_phi_dphi", "", 100, -180.0, 180.0, 100, 170.0, 190.0); + hi_pr_phi_dphi.setTitleX("#phi (deg)"); hi_pr_phi_dphi.setTitleY("#Delta#phi (deg)"); + H2F hi_pr_phi_dz = new H2F("hi_phi_dz", "", 100, -180.0, 180.0, 100, -5.0, 5.0); + hi_pr_phi_dz.setTitleX("#phi (deg)"); hi_pr_phi_dz.setTitleY("#Deltaz (cm)"); + H1F hi_pr_dp = new H1F("hi_dp", "Ebeam (GeV)", "Counts", 100, ebeam*0.75, ebeam*1.2); + H1F hi_pr_dtheta = new H1F("hi_dtheta", "#Delta#theta (deg)", "Counts", 100, -15.0, 15.0); + H1F hi_pr_dphi = new H1F("hi_dphi", "#Delta#phi (deg)", "Counts", 100, 170.0, 190.0); + H1F hi_pr_dz = new H1F("hi_dz", "#Deltaz (cm)", "Counts", 100, -5.0, 5.0); + dg_proton.addDataSet(hi_pr_p_dp, 0); + dg_proton.addDataSet(hi_pr, 1); + dg_proton.addDataSet(f1_pr, 1); + dg_proton.addDataSet(hi_pr_phi_dphi, 2); + dg_proton.addDataSet(hi_pr_phi_dz, 3); + dg_proton.addDataSet(hi_pr_dp, 4); + dg_proton.addDataSet(hi_pr_dtheta, 5); + dg_proton.addDataSet(hi_pr_dphi, 6); + dg_proton.addDataSet(hi_pr_dz, 7); + histos.put("Proton", dg_proton); + // Phi + DataGroup dg_phi = new DataGroup(2,3); + for(int sector=1; sector <= 6; sector++) { + H1F hi_dphi_sec = new H1F("hi_dphi_" + sector, "hi_dphi_" + sector, 100, 160.0, 200.0); + hi_dphi_sec.setTitleX("#Delta#phi (deg)"); + hi_dphi_sec.setTitleY("Counts"); + hi_dphi_sec.setTitle("Sector " + sector); + dg_phi.addDataSet(hi_dphi_sec, sector-1); + } + histos.put("Phi", dg_phi); + // Beam + DataGroup dg_beam = new DataGroup(2,3); + for(int sector=1; sector <= 6; sector++) { + H1F hi_beam_sec = new H1F("hi_beam_" + sector, "hi_beam_" + sector, 100, ebeam*0.75, ebeam*1.2); + hi_beam_sec.setTitleX("Beam Energy (GeV)"); + hi_beam_sec.setTitleY("Counts"); + hi_beam_sec.setTitle("Sector " + sector); + dg_beam.addDataSet(hi_beam_sec, sector-1); + } + histos.put("Beam", dg_beam); + } + + public void processEvent(DataEvent event) { + + DataBank recEvent = null; + DataBank recPart = null; + DataBank recScint = null; + DataBank recCal = null; + DataBank recTracks = null; + DataBank recTraj = null; + + if(event.hasBank("REC::Event")) recEvent = event.getBank("REC::Event"); + if(event.hasBank("REC::Particle")) recPart = event.getBank("REC::Particle"); + if(event.hasBank("REC::Scintillator")) recScint = event.getBank("REC::Scintillator"); + if(event.hasBank("REC::Calorimeter")) recCal = event.getBank("REC::Calorimeter"); + if(event.hasBank("REC::Track")) recTracks = event.getBank("REC::Track"); + if(event.hasBank("REC::Traj")) recTraj = event.getBank("REC::Traj"); + + Particle recEl = null; + Particle recPr = null; + Particle recPip = null; + Particle recPim = null; + Particle recPro = null; + LorentzVector virtualPhoton = null; + LorentzVector hadronSystem = null; + LorentzVector virtualPhotonP = null; + LorentzVector hadronSystemP = null; + Particle beam = Particle.createWithPid(11, 0,0,ebeam, 0,0,0); + Particle target = Particle.createWithPid(2212, 0,0,0, 0,0,0); + if(event.hasBank("REC::Particle")==true && event.hasBank("REC::Track")){ + DataBank bank = event.getBank("REC::Particle"); + DataBank track = event.getBank("REC::Track"); + int rows = bank.rows(); + for(int loop = 0; loop < rows; loop++){ + int status = (int) Math.floor(Math.abs(bank.getShort("status", loop))/1000); /// + if(loop==0 && bank.getInt("pid", loop)==11 && status==2) { + recEl = new Particle( + bank.getInt("pid", loop), + bank.getFloat("px", loop), + bank.getFloat("py", loop), + bank.getFloat("pz", loop), + bank.getFloat("vx", loop), + bank.getFloat("vy", loop), + bank.getFloat("vz", loop)); + for(int j=0; j0){ + histos.get("General").getH2F("hi_q2w").fill(hadronSystem.mass(),-virtualPhoton.mass2()); + histos.get("General").getH1F("hi_w").fill(hadronSystem.mass()); + histos.get("General").getH2F("hi_z_phi").fill(Math.toDegrees(recEl.phi()), recEl.vz()); + histos.get("General").getH2F("hi_w_phi").fill(Math.toDegrees(recEl.phi()), hadronSystem.mass()); + histos.get("General").getH2F("hi_w_theta").fill(Math.toDegrees(recEl.theta()), hadronSystem.mass()); + histos.get("General").getH2F("hi_el").fill(recEl.p(),Math.toDegrees(recEl.theta())); + histos.get("Elastic").getH1F("hi_w").fill(hadronSystem.mass()); + histos.get("Elastic").getH2F("hi_w_phi").fill(Math.toDegrees(recEl.phi()), hadronSystem.mass()); + histos.get("Elastic").getH2F("hi_w_theta").fill(Math.toDegrees(recEl.theta()), hadronSystem.mass()); + histos.get("Elastic").getH2F("hi_w_p").fill(recEl.p(),hadronSystem.mass()); + histos.get("W").getH1F("hi_w_" + secEl).fill(hadronSystem.mass()); + histos.get("Vertex").getH1F("h1_vz_el").fill(recEl.vz()); + histos.get("Vertex").getH2F("h2_theta_el").fill(recEl.vz(), Math.toDegrees(recEl.theta())); + histos.get("Vertex").getH2F("h2_phi_el").fill(recEl.vz()), Math.toDegrees(recEl.phi()); + histos.get("Vertex").getH2F("h2_p_el").fill(recEl.vz(), recEl.p()); + histos.get("Vertex-Sector").getH1F("h1_vz_el_"+secEl).fill(recEl.vz()); + } + if(hadronSystem.mass()<1.1) { + if(recPr != null) { + double phPr = Math.toDegrees(recPr.phi()); + if(phPr < phEl) phPr +=360; + histos.get("Proton").getH2F("hi_pr").fill(recPr.p(),Math.toDegrees(recPr.theta())); + // dg_proton.getH2F("hi_phi").fill(Math.toDegrees(recEl.phi()),Math.toDegrees(recPr.phi())); + // dg_proton.getH2F("hi_dphi").fill(Math.toDegrees(recPr.phi()),phPr-phEl-180); + histos.get("Phi").getH1F("hi_dphi_" + secEl).fill(phPr-phEl); + histos.get("Proton").getH1F("hi_dphi").fill(phPr-phEl); + histos.get("Proton").getH2F("hi_phi_dphi").fill(Math.toDegrees(recPr.phi()),phPr-phEl); + histos.get("Proton").getH2F("hi_phi_dz").fill(Math.toDegrees(recPr.phi()),recPr.vz()-recEl.vz()); + histos.get("Proton").getH1F("hi_dz").fill(recPr.vz()-recEl.vz()); + if(Math.abs(phPr-phEl-180)<10 && recPr.getProperty("NDF")>=0) { + histos.get("Proton").getH2F("hi_p_dp").fill(recPr.p(),-0.93832+recPr.e()+recEl.p()); + histos.get("Proton").getH1F("hi_dp").fill(-0.93832+recPr.e()+recEl.p()); + histos.get("Proton").getH1F("hi_dtheta").fill(Math.toDegrees(recPr.theta())-histos.get("Proton").getF1D("f1_pr").evaluate(recPr.p())); + histos.get("Beam").getH1F("hi_beam_" + secEl).fill((-0.93832+recPr.e()+recEl.p())); + if(lundfile!=null) { + PhysicsEvent ev = writeToLund(ebeam, recEl); + if(ev!=null) lundfile << ev.toLundString(); + } + } + } + } + if(recPip!=null) { + histos.get("Vertex").getH1F("h1_vz_pi+").fill(recPip.vz()); + histos.get("Vertex").getH2F("h2_theta_pi+").fill(recPip.vz(), Math.toDegrees(recPip.theta())); + histos.get("Vertex").getH2F("h2_phi_pi+").fill(recPip.vz(), Math.toDegrees(recPip.phi())); + histos.get("Vertex").getH2F("h2_p_pi+").fill(recPip.vz(), recPip.p()); + histos.get("Vertex-Sector").getH1F("h1_vz_pi+_"+((int) recPip.getProperty("sector"))).fill(recPip.vz()); + } + if(recPim!=null) { + histos.get("Vertex").getH1F("h1_vz_pi-").fill(recPim.vz()); + histos.get("Vertex").getH2F("h2_theta_pi-").fill(recPim.vz(), Math.toDegrees(recPim.theta())); + histos.get("Vertex").getH2F("h2_phi_pi-").fill(recPim.vz(), Math.toDegrees(recPim.phi())); + histos.get("Vertex").getH2F("h2_p_pi-").fill(recPim.vz(), recPim.p()); + histos.get("Vertex-Sector").getH1F("h1_vz_pi-_"+((int) recPim.getProperty("sector"))).fill(recPim.vz()); + } + } + if(recEl!=null && recPip!=null && recPim!=null) { + recPro = new Particle(); + recPro.copy(target); + recPro.combine(beam, +1); + recPro.combine(recEl, -1); + recPro.combine(recPip, -1); + recPro.combine(recPim, -1); + Particle rho = new Particle(); + rho.copy(recPip); + rho.combine(recPim, +1); + histos.get("2pi").getH1F("mxt1").fill(recPro.mass()); + histos.get("2pi").getH1F("mit1").fill(rho.mass()); + } + else if(recEl!=null && recPip!=null && recPro!=null && recPim==null) { + recPim = new Particle(); + recPim.copy(target); + recPim.combine(beam, +1); + recPim.combine(recEl, -1); + recPim.combine(recPip, -1); + recPim.combine(recPro, -1); + Particle rho = new Particle(); + rho.copy(recPip); + rho.combine(recPim, +1); + histos.get("2pi").getH1F("mxt2").fill(recPim.mass2()); + histos.get("2pi").getH1F("mit2").fill(rho.mass()); + } + else if(recEl!=null && recPim!=null && recPro!=null && recPip==null) { + recPip = new Particle(); + recPip.copy(target); + recPip.combine(beam, +1); + recPip.combine(recEl, -1); + recPip.combine(recPim, -1); + recPip.combine(recPro, -1); + Particle rho = new Particle(); + rho.copy(recPip); + rho.combine(recPim, +1); + histos.get("2pi").getH1F("mxt3").fill(recPip.mass2()); + histos.get("2pi").getH1F("mit3").fill(rho.mass()); + } + else if(recEl!=null && recPip!=null && recPim==null && recPro==null) { + Particle neutron = new Particle(); + neutron.copy(target); + neutron.combine(beam, +1); + neutron.combine(recEl, -1); + neutron.combine(recPip, -1); + Particle W = new Particle(); + W.copy(target); + W.combine(beam, +1); + W.combine(recEl, -1); + histos.get("1pi").getH2F("W1").fill(Math.toDegrees(recEl.phi()), neutron.mass()); + histos.get("1pi").getH1F("mxt1").fill(neutron.mass()); + } + else if(recEl!=null && recPro!=null && recPim==null && recPip==null) { + Particle pizero = new Particle(); + pizero.copy(target); + pizero.combine(beam, +1); + pizero.combine(recEl, -1); + pizero.combine(recPro, -1); + Particle W = new Particle(); + W.copy(target); + W.combine(beam, +1); + W.combine(recEl, -1); + histos.get("1pi").getH2F("W2").fill(Math.toDegrees(recEl.phi()), pizero.mass2()); + histos.get("1pi").getH1F("mxt2").fill(pizero.mass2()); + } + } + + } + + public void analyzeHistos() { + Logger.getLogger("org.freehep.math.minuit").setLevel(Level.WARNING); + + fitW(histos.get("General").getH1F("hi_w"), 0.8, 1.1, 0.05); + fitW(histos.get("Elastic").getH1F("hi_w"), 0.8, 1.1, 0.05); + fitW(histos.get("2pi").getH1F("mxt1"), 0.8, 1.1, 0.05); + fitW(histos.get("2pi").getH1F("mxt2"), -0.1, 0.2, 0.05); + fitW(histos.get("2pi").getH1F("mxt3"), -0.1, 0.2, 0.05); + fitW(histos.get("1pi").getH1F("mxt1"), 0.8, 1.1, 0.05); + fitW(histos.get("1pi").getH1F("mxt2"), -0.1, 0.2, 0.05); + for(int sector=1; sector <= 6; sector++) { + fitW(histos.get("W").getH1F("hi_w_" + sector), 0.8, 1.1, 0.05); + fitGauss(histos.get("Phi").getH1F("hi_dphi_" + sector)); + fitGauss(histos.get("Beam").getH1F("hi_beam_" + sector)); + } + fitGauss(histos.get("Proton").getH1F("hi_dp")); + fitGauss(histos.get("Proton").getH1F("hi_dtheta")); + fitGauss(histos.get("Proton").getH1F("hi_dphi")); + fitGauss(histos.get("Proton").getH1F("hi_dz")); +// fitW(histos.get("Vertex").getH1F("h1_vz_el"), scWindow+targetPos-2,scWindow+targetPos+2, 1); +// fitW(histos.get("Vertex").getH1F("h1_vz_pi+"), scWindow+targetPos-2,scWindow+targetPos+2, 1); +// fitW(histos.get("Vertex").getH1F("h1_vz_pi-"), scWindow+targetPos-2,scWindow+targetPos+2, 1); + } + + public EmbeddedCanvasTabbed drawHistos(String optStats) { + + EmbeddedCanvasTabbed canvas = null; + + String[] titles = (String[]) histos.keySet().toArray(); + for(int i=0; i dsList = histos.get(key).getData(i); + for(IDataSet ds : dsList){ + // System.out.println("\t --> " + ds.getName()); + dir.addDataSet(ds); + } + } + } + dir.writeFile(filename); + } + + public void readHistos(String filename) { + this.createHistos(); + TDirectory dir = new TDirectory(); + dir.readFile(filename); + System.out.println(dir.getDirectoryList()); + for(String key : histos.keySet()) { + String folder = key + "/"; + System.out.println("Reading from: " + folder); + DataGroup group = this.histos.get(key); + int nrows = group.getRows(); + int ncols = group.getColumns(); + int nds = nrows*ncols; + DataGroup newGroup = new DataGroup(ncols,nrows); + for(int i = 0; i < nds; i++){ + List dsList = group.getData(i); + for(IDataSet ds : dsList){ + System.out.println("\t --> " + ds.getName()); + newGroup.addDataSet(dir.getObject(folder, ds.getName()),i); + } + } + this.histos.replace(key, newGroup); + } + } + + private void fitW(H1F hiw, double min, double max, double sigma) { + + F1D f1w = new F1D("f1_w", "[amp]*gaus(x,[mean],[sigma])", min, max); + // get histogram maximum in the rane 0.7-1.1 + int i1=hiw.getXaxis().getBin(min); + int i2=hiw.getXaxis().getBin(max); + double hiMax=0; + int imax=i1; + for(int i=i1; i<=i2; i++) { + if(hiMax inputList = null; + + Kinematics analysis = null; + + if(parser.getCommand().equals("-process")) { + int maxEvents = parser.getOptionParser("-process").getOption("-nevent").intValue(); + double beamEnergy = parser.getOptionParser("-process").getOption("-beam").doubleValue(); + String vertexPar = parser.getOptionParser("-process").getOption("-vertpar").stringValue(); + String namePrefix = parser.getOptionParser("-process").getOption("-o").stringValue(); + String histoName = "histo.hipo"; + if(!namePrefix.isEmpty()) { + histoName = namePrefix + "_" + histoName; + } + optStats = parser.getOptionParser("-process").getOption("-stats").stringValue(); + openWindow = parser.getOptionParser("-process").getOption("-display").intValue()!=0; + if(!openWindow) System.setProperty("java.awt.headless", "true"); + + inputList = parser.getOptionParser("-process").getInputList(); + if(inputList.isEmpty()==true){ + parser.printUsage(); + System.out.println("\n >>>> error: no input file is specified....\n"); + System.exit(0); + } + + analysis = new Kinematics(beamEnergy, vertexPar); + + ProgressPrintout progress = new ProgressPrintout(); + + int counter = -1; + analysis.createHistos(); + for(String inputFile : inputList){ + HipoDataSource reader = new HipoDataSource(); + reader.open(inputFile); + + + while (reader.hasEvent()) { + + counter++; + + DataEvent event = reader.getNextEvent(); + analysis.processEvent(event); + + progress.updateStatus(); + if(maxEvents>0){ + if(counter>=maxEvents) break; + } + } + progress.showStatus(); + reader.close(); + } + analysis.analyzeHistos(); + analysis.saveHistos(histoName); + } + + if(parser.getCommand().equals("-plot")) { + double beamEnergy = parser.getOptionParser("-plot").getOption("-beam").doubleValue(); + String vertexPar = parser.getOptionParser("-plot").getOption("-vertpar").stringValue(); + optStats = parser.getOptionParser("-plot").getOption("-stats").stringValue(); + openWindow = parser.getOptionParser("-plot").getOption("-display").intValue()!=0; + if(!openWindow) System.setProperty("java.awt.headless", "true"); + + inputList = parser.getOptionParser("-plot").getInputList(); + if(inputList.isEmpty()==true){ + parser.printUsage(); + System.out.println("\n >>>> error: no input file is specified....\n"); + System.exit(0); + } + + analysis = new Kinematics(beamEnergy, vertexPar); + analysis.readHistos(inputList.get(0)); + analysis.analyzeHistos(); + } + + + if(openWindow) { + JFrame frame = new JFrame(inputList.get(0)); + frame.setSize(1200, 800); + frame.add(analysis.drawHistos(optStats)); + frame.setLocationRelativeTo(null); + frame.setVisible(true); + } + } +}