4 modified files
hps-java/src/main/java/org/lcsim/hps/users/phansson
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;
+ }
+
+
+
+
}
hps-java/src/main/java/org/lcsim/hps/users/phansson
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) {
hps-java/src/main/java/org/lcsim/hps/users/phansson
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);
hps-java/src/main/java/org/lcsim/hps/users/phansson
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);
+ }
+
}
}
CVSspam 0.2.12