hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.1 -r1.2
--- RunMPAlignment.java 23 May 2012 01:04:12 -0000 1.1
+++ RunMPAlignment.java 30 May 2012 18:24:44 -0000 1.2
@@ -1,14 +1,16 @@
package org.lcsim.hps.users.phansson;
+import hep.aida.*;
import org.lcsim.hps.users.phansson.MPAlignmentParameters;
-import hep.aida.IAnalysisFactory;
import java.io.IOException;
+import java.util.ArrayList;
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
+import org.lcsim.geometry.Detector;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -33,14 +35,29 @@
// note: this should be -1 for Test configurations and +1 for Full (v3.X and lower) configurations
// this is set by the _config variable (detType in HeavyPhotonDriver)
int flipSign = 1;
- boolean _DEBUG = true;
-
+ boolean _DEBUG = false;
+ double _resLayer1MinY;
+// _resLayer2MinY,_resLayer3MinY,_resLayer4MinY,_resLayer5MinY,_resLayer6MinY,_resLayer7MinY,_resLayer8MinY,_resLayer9MinY,_resLayer10MinY;
+// double _resLayer1MinX,_resLayer2MinX,_resLayer3MinX,_resLayer4MinX,_resLayer5MinX,_resLayer6MinX,_resLayer7MinX,_resLayer8MinX,_resLayer9MinX,_resLayer10MinX;
+// double _resLayer1MaxY,_resLayer2MaxY,_resLayer3MaxY,_resLayer4MaxY,_resLayer5MaxY,_resLayer6MaxY,_resLayer7MaxY,_resLayer8MaxY,_resLayer9MaxY,_resLayer10MaxY;
+// double _resLayer1MaxX,_resLayer2MaxX,_resLayer3MaxX,_resLayer4MaxX,_resLayer5MaxX,_resLayer6MaxX,_resLayer7MaxX,_resLayer8MaxX,_resLayer9MaxX,_resLayer10MaxX;
+//
+
+
+
+
+
+
+
+
+
+
public RunMPAlignment() {
nlayers[0] = 10;
_minLayers = 8;
//// if (_config.contains("Test"))
//// flipSign = -1;
- ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt");
+ ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_DEBUG);
}
@@ -51,10 +68,64 @@
_config = config;
if (_config.contains("Test"))
flipSign = -1;
- ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt");
+ ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_DEBUG);
}
+
+ public void detectorChanged(Detector detector) {
+
+
+
+
+ //ap.setResLimits();
+ for(int i=1;i<=10;++i) {
+ for(int j=0;j<3;++j) {
+ double xmin = -50;
+ double xmax = 50;
+// if(i==3 && j==0) {
+// xmin = -1;
+// xmax = 0;
+// }
+ //xmin = _resLayer1MinY;
+ //}
+ for(int side=0;side<2;++side) {
+ ap.setResLimits(side,i,j, xmin, xmax);
+
+ }
+ }
+ }
+
+ //ap.setResLimits(top/bottom,layer,direction: u/v/w, min, max);
+ ap.setResLimits(0,1,0, -0.3, 0.3);
+ ap.setResLimits(0,2,0, -0.3, 0.3);
+ ap.setResLimits(0,3,0, -0.4, 0.1);
+ ap.setResLimits(0,4,0, -0.2, 0.3);
+ ap.setResLimits(0,5,0, -0.1, 0.4);
+ ap.setResLimits(0,6,0, -0.3,0.1);
+ ap.setResLimits(0,7,0, -0.4, 0.45);
+ ap.setResLimits(0,8,0, -0.5, 0.3);
+ ap.setResLimits(0,9,0, -1, 1);
+ ap.setResLimits(0,10,0, -1, 1);
+
+ ap.setResLimits(1,1,0, -0.2,0.2);
+ ap.setResLimits(1,2,0, -0.2, 0.2);
+ ap.setResLimits(1,3,0, -0.2, 0.2);
+ ap.setResLimits(1,4,0, -0.3, 0.3);
+ ap.setResLimits(1,5,0, -0.2, 0.4);
+ ap.setResLimits(1,6,0, -0.5, 0.);
+ ap.setResLimits(1,7,0, -0.7, 0.15);
+ ap.setResLimits(1,8,0, -0.1, 0.7);
+ ap.setResLimits(1,9,0, -1, 0.7);
+ ap.setResLimits(1,10,0, -0.5, 1);
+
+
+
+
+ }
+
+
+
public void process(
EventHeader event) {
@@ -106,4 +177,9 @@
Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE, null, ex);
}
}
+
+
+ public void setResLayer1MinY(double val) {
+ this._resLayer1MinY = val;
+ }
}
hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.2 -r1.3
--- MPAlignmentParameters.java 24 May 2012 01:30:31 -0000 1.2
+++ MPAlignmentParameters.java 30 May 2012 18:24:44 -0000 1.3
@@ -1,5 +1,6 @@
package org.lcsim.hps.users.phansson;
+import hep.aida.*;
import org.lcsim.hps.users.mgraham.alignment.*;
import hep.physics.matrix.BasicMatrix;
import hep.physics.matrix.MatrixOp;
@@ -30,8 +31,10 @@
import org.lcsim.fit.helicaltrack.MultipleScatter;
import org.lcsim.fit.helicaltrack.TrackDirection;
import org.lcsim.hps.event.HPSTransformations;
+import org.lcsim.hps.monitoring.AIDAFrame;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.util.aida.AIDA;
/**
* Class to calculate and print the residuals and derivatives
@@ -55,22 +58,45 @@
private int _nlc = 5; //the five track parameters
private int _ngl = 1; //delta(u) and delta(gamma) for each plane
- private GlobalParameters _glp;
+ private GlobalParameters _globalParameterList;
+ private List<GlobalParameter> _glp = new ArrayList<GlobalParameter>();
private BasicMatrix _dfdq;
//private BasicMatrix _dfdp;
private HelicalTrackFit _trk;
private double[] _resid = new double[3];
private double[] _error = new double[3];
- private int[] _globalLabel = new int[1];
FileWriter fWriter;
PrintWriter pWriter;
Set<SiSensor> _process_sensors = new HashSet<SiSensor>();
private HPSTransformations _detToTrk;
- boolean _DEBUG = true;
+ boolean _DEBUG=false;
double smax = 1e3;
+ ResLimit _resLimits = new ResLimit();
-
- public MPAlignmentParameters(String outfile) {
+ private AIDAFrame plotterFrame;
+ private AIDA aida = AIDA.defaultInstance();
+ private IAnalysisFactory af = aida.analysisFactory();
+ //private List< List<IHistogram1D> > hres = new ArrayList<List<IHistogram1D> >();
+// private List< List<IHistogram1D> > resy_org_layallhit = new ArrayList<List<IHistogram1D> >();
+// private List< List<IHistogram1D> > resy_org_layallsinglehit = new ArrayList<List<IHistogram1D> >();
+// private List< List<IHistogram1D> > resy_org_lay1hit = new ArrayList<List<IHistogram1D> >();
+// private List< IHistogram1D > nhits_tracker = new ArrayList<IHistogram1D>();
+// private List< List<IHistogram1D> > nhits_tracker_layer = new ArrayList<List<IHistogram1D> >();
+// private List< IHistogram1D > ncl_ecal = new ArrayList<IHistogram1D>();
+// private List< IHistogram1D > selcl_ecal_e = new ArrayList<IHistogram1D>();
+// private List< IHistogram1D > cl_ecal_e = new ArrayList<IHistogram1D>();
+// private List< IHistogram2D > ncl_ecal_map = new ArrayList<IHistogram2D>();
+// private List< IHistogram2D > nselcl_ecal_map = new ArrayList<IHistogram2D>();
+// private List< IHistogram2D> resy_2d_org_layallhit = new ArrayList<IHistogram2D>();
+// private List< IHistogram2D> resy_2d_org_layallsinglehit = new ArrayList<IHistogram2D>();
+//
+
+
+
+
+
+ public MPAlignmentParameters(String outfile, boolean debug) {
+ _DEBUG = debug;
_detToTrk = new HPSTransformations();
try {
//open things up
@@ -80,11 +106,88 @@
Logger.getLogger(RunAlignment.class.getName()).log(Level.SEVERE, null, ex);
}
- _glp = new GlobalParameters();
- if(_DEBUG) _glp.print();
+ _globalParameterList = new GlobalParameters();
+ if(_DEBUG) _globalParameterList.print();
+
+
+
+
+
+
+
+ aida.tree().cd("/");
+ plotterFrame = new AIDAFrame();
+ plotterFrame.setTitle("HPS Alignment Plots");
+
+
+
+ //IHistogramFactory hf = aida.histogramFactory();
+ String[] side = {"top","bottom"};
+ String[] direction = {"u","v","w"};
+
+ for (int d=0;d<3;++d) {
+
+ for (int s=0;s<2;++s) {
+ // if(iSide==1) continue;
+ IPlotter plotter_res = af.createPlotterFactory().create();
+ plotter_res.createRegions(5,2,0);
+ plotter_res.setTitle("res_" + direction[d] + " " + side[s]);
+ IPlotter plotter_reserr = af.createPlotterFactory().create();
+ plotter_reserr.createRegions(5,2,0);
+ plotter_reserr.setTitle("reserr_" + direction[d] + " " + side[s]);
+ IPlotter plotter_respull = af.createPlotterFactory().create();
+ plotter_respull.createRegions(5,2,0);
+ plotter_respull.setTitle("respull_" + direction[d] + " " + side[s]);
+
+ for (int iLayer=1;iLayer<11;++iLayer) {
+ double xmin,xmax;
+ if(direction[d]=="v") {
+ xmin = -100;
+ xmax = 100;
+ } else {
+ xmin = -2;
+ xmax = 2;
+
+ }
+ // IHistogram h = aida.histogram1D("res_"+direction[d]+"_layer" + (iLayer) + "_" + side[iSide] , 50, -10, 10);
+ IHistogram h = aida.histogram1D("res_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, xmin, xmax);
+ IHistogram h1 = aida.histogram1D("reserr_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, 0, 1);
+ IHistogram h2 = aida.histogram1D("respull_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, -10, 10);
+
+ int region = (iLayer-1);//*2+iSide;
+
+ plotter_res.region(region).plot(h);
+ plotter_reserr.region(region).plot(h1);
+ plotter_respull.region(region).plot(h2);
+
+ }
+
+ //plotter_res.show();
+ plotterFrame.addPlotter(plotter_res);
+ plotterFrame.addPlotter(plotter_reserr);
+ plotterFrame.addPlotter(plotter_respull);
+ }
+ }
+
+ plotterFrame.pack();
+ plotterFrame.setVisible(true);
+
+
+
+
+
+ }
+
+ public void setResLimits(int l,int d, double low,double high) {
+ _resLimits.add(0,l,d, low,high); //top
+ _resLimits.add(1,l,d, low,high); //bottom
}
+ public void setResLimits(int s, int l,int d, double low,double high) {
+ _resLimits.add(s,l,d, low,high);
+ }
+
public void PrintResidualsAndDerivatives(Track track) {
SeedTrack st = (SeedTrack) track;
@@ -92,6 +195,7 @@
Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap();
_trk = seed.getHelix();
List<TrackerHit> hitsOnTrack = track.getTrackerHits();
+ pWriter.printf("TRACK\n");
for (TrackerHit hit : hitsOnTrack) {
HelicalTrackHit htc = (HelicalTrackHit) hit;
double msdrphi = msmap.get(htc).drphi();
@@ -106,7 +210,9 @@
CalculateGlobalDerivatives(cl);
CalculateResidual(cl, msdrphi, msdz);
// CalculateResidual(cl, 0,0);
- PrintStripResiduals(cl);
+ //PrintStripResiduals(cl);
+ PrintStripResidualsNew(cl);
+
}
}
//AddTarget(0.1, 0.02);
@@ -187,19 +293,23 @@
//Find out which global derivative this is
//align only layer 3 on top side!
- int side = 10000;
+ int side = 20000;
if( strip.origin().z()>0) side = 10000;
- else side = 20000;
- int l = strip.layer();
int type = 1000;
- int dir = 300;
- int label = side + type + dir + l;
//if(strip.origin().z()>0 && strip.layer()==3) {
+ //Flag to tell if this residual is affected by the given global parameter
+ boolean useGL = false;
+
+ //Clear the old parameter list
+ _glp.clear();
+
//Go through each gl parameter and see if it has non-zero contribution
- for (GlobalParameter gp : _glp.getList()) {
- //if(gp.active()==false) continue;
-
+ for (GlobalParameter gp : _globalParameterList.getList()) {
+
+
+ useGL = false;
+
dfdpTRACK[0][0] = 0;
dfdpTRACK[1][0] = 0;
dfdpTRACK[2][0] = 0;
@@ -213,6 +323,10 @@
//Calcuate dfdp derivatives depending on type of global parameter
if(gp.getType()==1000) {
+
+ //This residual is affected by a global parameter
+ useGL = true;
+
//translation
if(gp.getDirection()==100) {
//x - tracking frame --> beamline direction
@@ -232,26 +346,59 @@
dfdpTRACK[0][0] = 0; //df/dx
dfdpTRACK[1][0] = 0; //df/dy
dfdpTRACK[2][0] = 1; //df/dz
+ }
+
+ }//type
+ else if(gp.getType()==2000) {
+ //rotation
+
+ //This residual is affected by a global parameter
+ useGL = true;
+
+
+ if(gp.getDirection()==100) {
+ //x - tracking frame --> beamline direction
+ dfdpTRACK[0][0] = 0; //df_x/d(alpha)
+ dfdpTRACK[1][0] = -1; //df_y/d(alpha)
+ dfdpTRACK[2][0] = 1; //df_z/d(alpha)
+ }
+ else if(gp.getDirection()==200) {
+ //y - tracking frame --> almost unmeasured coordinate direction
+ dfdpTRACK[0][0] = 1; //df_x/d(beta)
+ dfdpTRACK[1][0] = 0; //df_y/d(beta)
+ dfdpTRACK[2][0] = -1; //df_z/d(beta)
+ }
+ else {
+ // Rotation around z
+ dfdpTRACK[0][0] = -1; //df_x/d(gamma)
+ dfdpTRACK[1][0] = 1; //df_y/d(gamma)
+ dfdpTRACK[2][0] = 0; //df_z/d(gamma)
}
}//type
- }//layer
-
+ }//layer
}//side
- //Put it into a matrix to be able to transform easily
- BasicMatrix _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
- //Get transformation matrix from tracking frame to sensor frame where the residuals are calculated
- Hep3Matrix trkToStrip = getTrackToStripRotation(strip);
- //Transform derivatives to sensor frame!
- BasicMatrix dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, _dfdpTRACK);
- //Add it to the global parameter object
- gp.setDfDp(dfdp);
- if (_DEBUG) {
- System.out.printf("dfdz = %5.5f %5.5f %5.5f GL%d name: %s\n", gp.dfdp(0), gp.dfdp(1), gp.dfdp(2), gp.getLabel(),gp.getName());
-
- }
+ if(useGL) {
+
+ //Put it into a matrix to be able to transform easily
+ BasicMatrix _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
+ //Get transformation matrix from tracking frame to sensor frame where the residuals are calculated
+ Hep3Matrix trkToStrip = getTrackToStripRotation(strip);
+ //Transform derivatives to sensor frame!
+ BasicMatrix dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, _dfdpTRACK);
+ //Add it to the global parameter object
+ gp.setDfDp(dfdp);
+ _glp.add(gp);
+ if (_DEBUG) {
+ System.out.printf("dfdp = %5.5f %5.5f %5.5f GL%d name: %s\n", gp.dfdp(0), gp.dfdp(1), gp.dfdp(2), gp.getLabel(),gp.getName());
+ }
+ } else {
+ if (_DEBUG) {
+ System.out.printf("Global parameters %s name %s didn't affect this strip on side %d and layer %d\n",gp.getLabel(),gp.getName(),side,strip.layer());
+ }
+ }
} //gps
@@ -290,7 +437,7 @@
double vmeas = 0;
double vError = (strip.vmax() - strip.vmin()) / Math.sqrt(12);
double wmeas = 0;
- double wError = 0.001;
+ double wError = 10.0/Math.sqrt(12); //0.001;
//System.out.println("strip error="+uError+"; ms error ="+msuError);
_resid[0] = umeas - umc;
_error[0] = Math.sqrt(uError * uError + msuError * msuError);
@@ -395,7 +542,7 @@
}
//String[] p = {"u-displacement"};
System.out.println("global parameter derivatives");
- for (GlobalParameter gl : _glp.getList()) {
+ for (GlobalParameter gl : _glp) {
System.out.printf("%s %5.5e %5.5e %5.5e %5d %10s\n", "", gl.dfdp(0), gl.dfdp(1), gl.dfdp(2), gl.getLabel(),gl.getName());
//System.out.printf("%s %5.5e %5.5e %5.5e %5d\n", p[j], _dfdp.e(0, j), _dfdp.e(1, j), _dfdp.e(2, j), _globalLabel[j]);
}
@@ -407,13 +554,80 @@
for (int i = 0; i < _nlc; i++) {
pWriter.printf("%5.5e %5.5e %5.5e\n", _dfdq.e(0, i), _dfdq.e(1, i), _dfdq.e(2, i));
}
- for (GlobalParameter gl: _glp.getList()) {
+ for (GlobalParameter gl: _glp) {
if(gl.active()){
pWriter.printf("%5.5e %5.5e %5.5e %5d\n", gl.dfdp(0), gl.dfdp(1), gl.dfdp(2), gl.getLabel());
}
}
}
+ private void PrintStripResidualsNew(HelicalTrackStrip strip) {
+ if (_DEBUG) {
+ int s = 1;
+ if(strip.origin().z()>0) s = 0;
+ System.out.printf("Strip Layer = %4d\n", strip.layer());
+ System.out.printf("Residuals (u,v,w) : %5.5e %5.5e %5.5e (limits: [%.1e,%.1e] [%.1e,%.1e] [%.1e,%.1e])\n"
+ , _resid[0], _resid[1], _resid[2],
+ _resLimits.getMin(s,strip.layer(),0),
+ _resLimits.getMax(s,strip.layer(),0),
+ _resLimits.getMin(s,strip.layer(),1),
+ _resLimits.getMax(s,strip.layer(),1),
+ _resLimits.getMin(s,strip.layer(),2),
+ _resLimits.getMax(s,strip.layer(),2));
+ System.out.printf("Errors (u,v,w) : %5.5e %5.5e %5.5e\n", _error[0], _error[1], _error[2]);
+ String[] q = {"d0", "z0", "slope", "phi0", "R"};
+ System.out.println("track parameter derivatives");
+ for (int i = 0; i < _nlc; i++) {
+ System.out.printf("%s %5.5e %5.5e %5.5e\n", q[i], _dfdq.e(0, i), _dfdq.e(1, i), _dfdq.e(2, i));
+ }
+ //String[] p = {"u-displacement"};
+ System.out.println("global parameter derivatives");
+ for (GlobalParameter gl : _glp) {
+ System.out.printf("%s %5.5e %5.5e %5.5e %5d %10s\n", "", gl.dfdp(0), gl.dfdp(1), gl.dfdp(2), gl.getLabel(),gl.getName());
+ //System.out.printf("%s %5.5e %5.5e %5.5e %5d\n", p[j], _dfdp.e(0, j), _dfdp.e(1, j), _dfdp.e(2, j), _globalLabel[j]);
+ }
+
+ }
+ pWriter.printf("%d\n", strip.layer());
+ String[] d = {"u","v","w"};
+ for(int j=0;j<3;++j){
+ String side = "bottom";
+ int s = 1;
+ if( strip.origin().z()>0) {
+ side = "top";
+ s = 0;
+ }
+
+ if(s!=0) continue;
+
+ if(!isAllowedResidual(s,strip.layer(),j,_resid[j])) {
+ if(_DEBUG) System.out.println("Layer " + strip.layer() + " in " + d[j] + " res " + _resid[j] + " +- " + _error[j] + " was outside allowed range");
+ continue;
+ }
+//// if(Math.abs(_resid[j]/_error[j])>3) {
+//// if(_DEBUG) System.out.println("Layer " + strip.layer() + " in " + d[j] + " res " + _resid[j] + " +- " + _error[j] + " had too high pull");
+//// continue;
+//// }
+ pWriter.printf("res_%s %5.5e %5.5e \n", d[j], _resid[j], _error[j]);
+ for (int i = 0; i < _nlc; i++) {
+ pWriter.printf("lc_%s %5.5e \n", d[j], _dfdq.e(j, i));
+ }
+ for (GlobalParameter gl: _glp) {
+ if(gl.active()){
+ pWriter.printf("gl_%s %5.5e %5d\n", d[j], gl.dfdp(j), gl.getLabel());
+ }
+ }
+
+ aida.histogram1D("res_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_resid[j]);
+ aida.histogram1D("reserr_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_error[j]);
+ aida.histogram1D("respull_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_resid[j]/_error[j]);
+
+
+ }
+ }
+
+
+
private Hep3Matrix getTrackToStripRotation(HelicalTrackStrip strip) {
//This function transforms the hit to the sensor coordinates
@@ -661,4 +875,16 @@
// System.out.println( _trk.xc()+ "; "+_trk.yc());
// System.out.println( _trk.x0()+ "; "+_trk.y0());
}
+
+
+
+ private boolean isAllowedResidual(int side,int layer,int d,double res) {
+
+ if(res <_resLimits.getMin(side,layer,d)) return false;
+ if(res >_resLimits.getMax(side,layer,d)) return false;
+ return true;
+
+ }
+
+
}