Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
RunMPAlignment.java | +100 | -19 | 1.2 -> 1.3 |
GlobalParameters.java | +1 | -1 | 1.2 -> 1.3 |
dataMCPlots.java | +10 | -9 | 1.1 -> 1.2 |
MPAlignmentParameters.java | +33 | -22 | 1.3 -> 1.4 |
+144 | -51 |
Only u residuals, tracks in separate halves only.
diff -u -r1.2 -r1.3 --- RunMPAlignment.java 30 May 2012 18:24:44 -0000 1.2 +++ RunMPAlignment.java 21 Aug 2012 19:34:18 -0000 1.3 @@ -1,6 +1,9 @@
package org.lcsim.hps.users.phansson; import hep.aida.*;
+import java.io.BufferedReader; +import java.io.FileNotFoundException; +import java.io.FileReader;
import org.lcsim.hps.users.phansson.MPAlignmentParameters; import java.io.IOException; import java.util.ArrayList;
@@ -10,6 +13,7 @@
import org.lcsim.event.EventHeader; import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.geometry.Detector; import org.lcsim.util.Driver; import org.lcsim.util.aida.AIDA;
@@ -29,13 +33,15 @@
String _config = ""; MPAlignmentParameters ap; int totalTracks=0;
+ int totalTracksProcessed=0; + private String _resLimitFileName="";
// flipSign is a kludge... // HelicalTrackFitter doesn't deal with B-fields in -ive Z correctly // so we set the B-field in +iveZ and flip signs of fitted tracks // 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 = false;
+ private 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;
@@ -48,7 +54,9 @@
-
+ public void setDebug(boolean v) { + this._debug = v; + }
@@ -57,7 +65,7 @@
_minLayers = 8; //// if (_config.contains("Test")) //// flipSign = -1;
- ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_DEBUG);
+ ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug);
}
@@ -68,34 +76,26 @@
_config = config; if (_config.contains("Test")) flipSign = -1;
- ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_DEBUG);
+ ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug);
} public void detectorChanged(Detector detector) {
- - - - //ap.setResLimits();
+ + //Initialize the res limits
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(side,i,j, xmin, xmax);
} } }
-
+ loadResidualLimits(); + /*
//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);
@@ -118,7 +118,7 @@
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);
-
+ */
@@ -136,7 +136,13 @@
double duRange=0.1; for (Track trk : tracklist) { totalTracks++;
- ap.PrintResidualsAndDerivatives(trk);
+ + if(hitsOnBothSides(trk)) continue; + + totalTracksProcessed++; + + + ap.PrintResidualsAndDerivatives(trk,totalTracks);
if(1==1){ aida.histogram1D("Track d0",50,-0.5,0.5).fill(trk.getTrackParameter(0));
@@ -172,14 +178,89 @@
public void endOfData() { try { System.out.println("Total Number of Tracks Found = "+totalTracks);
+ System.out.println("Total Number of Tracks Processed = "+totalTracksProcessed);
ap.closeFile(); } catch (IOException ex) { Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE, null, ex); } }
+ + public boolean hitsOnBothSides(Track track) { + List<TrackerHit> hitsOnTrack = track.getTrackerHits(); + int ihalf = hitsOnTrack.get(0).getPosition()[2]>0 ? 1 : 0; + boolean hitsOnWrongSide = false; + for (TrackerHit hit : hitsOnTrack) { + double[] pos = hit.getPosition(); + if((ihalf==0 && pos[2]>0) || (ihalf==1 && pos[2]<0)) { + hitsOnWrongSide = true; + break; + } + } + if(hitsOnWrongSide) { + System.out.println("TRACK w/ both halves hit (: chi2 "+track.getChi2()+", pX "+track.getPX()+", pY "+track.getPY()+", pZ "+track.getPZ()+")"); + System.out.printf("Hits: "); + for (TrackerHit hit : hitsOnTrack) { + double[] pos = hit.getPosition(); + System.out.printf("(%.2f,%.2f,%.2f)", pos[0],pos[1],pos[2]); + } + System.out.println(""); + } + return hitsOnWrongSide; + } + + private void loadResidualLimits() { + + FileReader fReader = null; + BufferedReader bReader = null; + try { + fReader = new FileReader(this._resLimitFileName); + bReader = new BufferedReader(fReader); + + String line; + while( (line = bReader.readLine()) != null) { + if (line.contains("#")) continue; + String[] vec = line.split("\\s+"); + if(vec.length!=5) { + System.out.println("Error: residual limits line has wrong format -> " + line); + System.exit(1); + } + try { + int side = Integer.parseInt(vec[0]); + int layer = Integer.parseInt(vec[1]); + int direction = Integer.parseInt(vec[2]); + double min = Double.parseDouble(vec[3]); + double max = Double.parseDouble(vec[4]); + ap.setResLimits(side, layer, direction, min, max); + } catch(NumberFormatException e) { + Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE,null,e); + } + } + bReader.close(); + fReader.close(); + } + catch (FileNotFoundException ex) { + Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE, null, ex); + } catch (IOException e) { + Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE,null,e); + } + + + ap.setResLimits(0,1,0, -0.3, 0.3); + + } + +
public void setResLayer1MinY(double val) { this._resLayer1MinY = val; }
+ + public void setResidualLimitFileName(String fileName) { + this._resLimitFileName = fileName; + } + + + +
}
diff -u -r1.2 -r1.3 --- GlobalParameters.java 30 May 2012 18:24:44 -0000 1.2 +++ GlobalParameters.java 21 Aug 2012 19:34:18 -0000 1.3 @@ -38,7 +38,7 @@
// 300 - z (measurement direction) // [Layer] // 1-10
-
+
boolean ON; for(int iside=0;iside<2;++iside) { for(int ilayer=1;ilayer<11;++ilayer) {
diff -u -r1.1 -r1.2 --- dataMCPlots.java 24 Jul 2012 23:27:19 -0000 1.1 +++ dataMCPlots.java 21 Aug 2012 19:34:18 -0000 1.2 @@ -109,15 +109,16 @@
if(type==1) { HashMap<Integer,DataCont> map_dc = new HashMap<Integer,DataCont>();
+ String path = "plots/20120723_ecal_dataMC/";
Integer runs[] = {1351,1354,1353}; Integer runsBkg[] = {1358,1358,1358};
- String mc[] = {"trigratefile_egs5_merged_5000mb_90na_0x016.aida", - "trigratefile_egs5_merged_20bb_90na_0x0045.aida", - "trigratefile_egs5_merged_50bb_90na_0x0018.aida"
+ String mc[] = {path+"trigratefile_egs5_merged_5000mb_90na_0x016.aida", + path+"trigratefile_egs5_merged_20bb_90na_0x0045.aida", + path+"trigratefile_egs5_merged_50bb_90na_0x0018.aida"
};
- String mcG4[] = {"trigratefile_g4_merged_5000mb_90na_0x016.aida", - "trigratefile_g4_merged_20bb_90na_0x0045.aida", - "trigratefile_g4_merged_50bb_90na_0x0018.aida"
+ String mcG4[] = {path+"trigratefile_g4_merged_5000mb_90na_0x016.aida", + path+"trigratefile_g4_merged_20bb_90na_0x0045.aida", + path+"trigratefile_g4_merged_50bb_90na_0x0018.aida"
}; double Q_normalization = 90.0; //nC i.e. 90nC is 1s of beam at 90nA
@@ -130,13 +131,13 @@
IHistogramFactory hf = aida.histogramFactory(); IDataPointSetFactory dpsf = af.createDataPointSetFactory(null); // ITree tree = aida.tree();//(ITreeFactory) af.createTreeFactory().create();
- IDataPointSet dpsObs = dpsf.create("dpsObs","Data vs MC",2);
+ IDataPointSet dpsObs = dpsf.create("dpsObs","Data",2);
IDataPointSet dpsPred = dpsf.create("dpsPred","EGS",2); IDataPointSet dpsPredG4 = dpsf.create("dpsPredG4","G4",2); //IHistogram1D hdata = aida.histogram1D("Observed",runs.length , 0, runs.length); for(int irun=0;irun<runs.length;++irun) {
- DataCont dc = getDataMC("trigratefile_run"+runs[irun]+".aida",mc[irun], "trigratefile_run"+runsBkg[irun]+".aida",outName+"_egs5_"+runs[irun],savePlots,Q_normalization);
+ DataCont dc = getDataMC(path+"trigratefile_run"+runs[irun]+".aida",mc[irun], path+"trigratefile_run"+runsBkg[irun]+".aida",outName+"_egs5_"+runs[irun],savePlots,Q_normalization);
dpsObs.addPoint(); dpsObs.point(irun).coordinate(1).setValue(dc.n_data); dpsObs.point(irun).coordinate(1).setErrorMinus(dc.en_data/2.0);
@@ -153,7 +154,7 @@
dpsPred.point(irun).coordinate(0).setErrorMinus(0.5); dpsPred.point(irun).coordinate(0).setErrorPlus(0.5);
- DataCont dcG4 = getDataMC("trigratefile_run"+runs[irun]+".aida",mcG4[irun], "trigratefile_run"+runsBkg[irun]+".aida",outName+"_g4_"+runs[irun],savePlots,Q_normalization);
+ DataCont dcG4 = getDataMC(path+"trigratefile_run"+runs[irun]+".aida",mcG4[irun], path+"trigratefile_run"+runsBkg[irun]+".aida",outName+"_g4_"+runs[irun],savePlots,Q_normalization);
dpsPredG4.addPoint(); dpsPredG4.point(irun).coordinate(1).setValue(dcG4.n_pred);
diff -u -r1.3 -r1.4 --- MPAlignmentParameters.java 30 May 2012 18:24:44 -0000 1.3 +++ MPAlignmentParameters.java 21 Aug 2012 19:34:18 -0000 1.4 @@ -188,14 +188,15 @@
_resLimits.add(s,l,d, low,high); }
- public void PrintResidualsAndDerivatives(Track track) {
+ public void PrintResidualsAndDerivatives(Track track, int itrack) {
SeedTrack st = (SeedTrack) track; SeedCandidate seed = st.getSeedCandidate(); Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap(); _trk = seed.getHelix(); List<TrackerHit> hitsOnTrack = track.getTrackerHits();
- pWriter.printf("TRACK\n");
+ String half = hitsOnTrack.get(0).getPosition()[2]>0 ? "top" : "bottom"; + pWriter.printf("TRACK %s (%d)\n",half,itrack);
for (TrackerHit hit : hitsOnTrack) { HelicalTrackHit htc = (HelicalTrackHit) hit; double msdrphi = msmap.get(htc).drphi();
@@ -283,7 +284,7 @@
//go through all the global parameters defined and calculate dfdp //The global derivatives is done in the tracking frame
- //THis means that after they are calculated I need to transform them
+ //This means that after they are calculated I need to transform them
//into the sensor frame because that is where the residuals are calculated double[][] dfdpTRACK = new double[3][1];
@@ -291,20 +292,17 @@
dfdpTRACK[1][0] = 0; //df/dy dfdpTRACK[2][0] = 0; //df/dz
- //Find out which global derivative this is - //align only layer 3 on top side! - int side = 20000; - if( strip.origin().z()>0) side = 10000; - int type = 1000; - //if(strip.origin().z()>0 && strip.layer()==3) {
+ //Which half of the detecotor? + int side = strip.origin().z()>0 ? 10000 : 20000; + int layer = strip.layer();
- //Flag to tell if this residual is affected by the given global parameter
+ //Flag to tell if this hit 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
+ //Go through each global parameter and see if it has non-zero contribution
for (GlobalParameter gp : _globalParameterList.getList()) {
@@ -318,7 +316,7 @@
// correct side
- if(gp.getLayer()==strip.layer()) {
+ if(gp.getLayer()==layer) {
//correct layer //Calcuate dfdp derivatives depending on type of global parameter
@@ -396,7 +394,7 @@
} } 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());
+ 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,layer);
} } } //gps
@@ -589,21 +587,24 @@
} pWriter.printf("%d\n", strip.layer());
- String[] d = {"u","v","w"}; - for(int j=0;j<3;++j){ - String side = "bottom"; - int s = 1;
+ // Loop over the three directions u,v,w and print residuals and derivatives + //String[] d = {"u","v","w"}; + String[] d = {"u"}; //phansson use only u direction! + String side; + int s; + for(int j=0;j<1;++j){ + side = "bottom"; + 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(!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;
@@ -615,6 +616,16 @@
for (GlobalParameter gl: _glp) { if(gl.active()){ pWriter.printf("gl_%s %5.5e %5d\n", d[j], gl.dfdp(j), gl.getLabel());
+ //x-check that side is correct + int lbl = gl.getLabel(); + Double df = Math.floor(lbl/10000); + int iside = (df.intValue() % 10) - 1; + //System.out.println("track is on " + s + " gl param is " + lbl + "(" + df + ","+iside+")" ); + if(iside!=s) { + System.out.println("ERROR track is on " + s + " while gl param is " + lbl + "(" + df + ")" ); + System.exit(1); + } +
} }
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1