hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.2 -r1.3
--- FastTrackResidualDriver.java 4 May 2012 18:23:35 -0000 1.2
+++ FastTrackResidualDriver.java 8 May 2012 17:20:29 -0000 1.3
@@ -25,7 +25,7 @@
/**
*
- * @author phansson
+ * @author phansson+
*/
public class FastTrackResidualDriver extends Driver {
@@ -43,6 +43,8 @@
double conversionZ;
private static int crystalCols;
private static int crystalRows;
+ double ecalBeamgapCorr;
+ int ecalClusterSel;
private String outputPlotFileName = "";
@@ -52,9 +54,13 @@
private List< List<IHistogram1D> > resy_org_layallhit = 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>();
public void startOfData() {
@@ -75,26 +81,49 @@
crystalCols = 46;
crystalRows = 5;
+ ecalBeamgapCorr = 20.0;
+
// Position of the conversion
conversionZ = -1000.0;
+ ecalClusterSel = 1;
+
IHistogramFactory hf = aida.histogramFactory();
String side;
for (int iSide=0;iSide<2;++iSide) {
if(iSide==0) side="up";
else side="down";
+
+ ncl_ecal.add(hf.createHistogram1D("FT_ncl_ecal_" + side , 20, 0, 20));
+ selcl_ecal_e.add(hf.createHistogram1D("FT_selcl_ecal_e" + side , 100, 0, 20000));
+ cl_ecal_e.add(hf.createHistogram1D("FT_cl_ecal_e" + side , 100, 0, 20000));
+
+
+ nhits_tracker.add(hf.createHistogram1D("FT_nhits_tracker_" + side , 15, 0, 15));
+
+ List<IHistogram1D> nhl_list = new ArrayList<IHistogram1D>();
+ nhits_tracker_layer.add(nhl_list);
+ for (int i=0;i<5;++i) {
+ nhl_list.add(hf.createHistogram1D("FT_nhits_tracker_" + side + "_layer" + i , 15, 0, 15));
+ }
+
+
+ double res_axis[] = {-40,40};
+
+
+ resy_2d_org_layallhit.add(hf.createHistogram2D("FT_resy_2d_org_layallhit_" + side, 5, 0.5, 5.5, 60, res_axis[0], res_axis[1]));
+
+
List<IHistogram1D> list = new ArrayList<IHistogram1D>();
resy_org.add(list);
List<IHistogram1D> listLayAllHit = new ArrayList<IHistogram1D>();
resy_org_layallhit.add(listLayAllHit);
List<IHistogram1D> listLay1Hit = new ArrayList<IHistogram1D>();
resy_org_lay1hit.add(listLay1Hit);
-
- nhits_tracker.add(hf.createHistogram1D("FT_nhits_tracker_" + side , 40, 0, 40));
- ncl_ecal.add(hf.createHistogram1D("FT_ncl_ecal_" + side , 60, 0, 60));
-
-
- //Setup the ecal 2D plot
+
+
+ //Setup the ecal 2D plot
+
if (side == "up") {
ncl_ecal_map.add(hf.createHistogram2D("FT_ecal_hitmap_" + side, 92, -46.0, 46.0, 5, 0.0, 5.0));
nselcl_ecal_map.add(hf.createHistogram2D("FT_sel_ecal_hitmap_" + side, 92, -46.0, 46.0, 5, 0.0, 5.0));
@@ -105,13 +134,13 @@
for (int iLayer=1;iLayer<6;++iLayer) {
- IHistogram1D h = hf.createHistogram1D("FT_res_" + side + "_l"+iLayer, 50, -25, 25);
+ IHistogram1D h = hf.createHistogram1D("FT_res_" + side + "_l"+iLayer, 50, res_axis[0], res_axis[1]);
list.add(h);
- IHistogram1D hAll = hf.createHistogram1D("FT_res_LayAllHit_" + side + "_l"+iLayer, 50, -25, 25);
+ IHistogram1D hAll = hf.createHistogram1D("FT_res_LayAllHit_" + side + "_l"+iLayer, 50,res_axis[0], res_axis[1]);
listLayAllHit.add(hAll);
- IHistogram1D hLay1 = hf.createHistogram1D("FT_res_Lay1Hit_" + side + "_l"+iLayer, 50, -25, 25);
+ IHistogram1D hLay1 = hf.createHistogram1D("FT_res_Lay1Hit_" + side + "_l"+iLayer, 50,res_axis[0], res_axis[1]);
listLay1Hit.add(hLay1);
}
@@ -130,6 +159,11 @@
public void setConversionZ(double z) {
this.conversionZ = z;
}
+
+ public void setecalBeamgapCorr(double val) {
+ this.ecalBeamgapCorr = val;
+ System.out.println("offset set: " + ecalBeamgapCorr);
+ }
public void setDebug(boolean flag) {
this.debug = flag;
@@ -143,6 +177,11 @@
this.EcalZPosition = val;
}
+ public void setEcalClusterSel(int id) {
+ this.ecalClusterSel = id;
+ }
+
+
public void process(EventHeader event) {
++nevents;
if( debug ) {
@@ -177,17 +216,20 @@
int nhits;
int nhitsInTracker;
int nhitsInLayer1;
- Hep3Vector ecal_cl;
Hep3Vector origin_pos;
- List<Hep3Vector> ecal_cls;
- FastTrack fastTrack;
- List<SiTrackerHitStrip1D> stripList;
+ List<Integer> ecal_cls;
+ FastTrack fastTrack = new FastTrack(debug);
+ fastTrack.setEcalBeamgapCorr(ecalBeamgapCorr); //Correct for geometry
+
+ //List<SiTrackerHitStrip1D> stripList;
boolean isaxial;
String name;
SiSensor siSensor;
int layer;
String si_side;
double res;
+ int selclids[];
+ int sel_ecal_idx;
int layerIndex = -1;
for (int iSide=0;iSide<2;++iSide) {
@@ -196,31 +238,44 @@
ecal_cls = getEcalClustersForFastTracking(ecal_all_clusters, sides[iSide]);
- if (debug) System.out.println("Found " + ecal_cls.size() +" Ecal clusters on the " + sides[iSide]);
+ if (debug) {
+ System.out.println("Found " + ecal_cls.size() +" Ecal clusters on the " + sides[iSide] + ": " + ecal_cls.toString());
+ }
+
//if ( 1==1 ) return;
ncl_ecal.get(iSide).fill(ecal_cls.size());
if( ecal_cls.size() ==0 ) {
- System.out.println("No clusters on this side...");
+ if(debug) System.out.println("No clusters on this side...");
continue;
}
//Fill map of Ecal hits
- for( Hep3Vector cl : ecal_cls) {
- int clpos[] = getCrystalPair(cl);
+ for( int icl=0; icl<ecal_cls.size(); ++icl) {
+ int clid = ecal_cls.get(icl);
+ double p_cl[] = ecal_all_clusters.get(clid).getPosition();
+ int clpos[] = getCrystalPair(p_cl);
ncl_ecal_map.get(iSide).fill(clpos[0], clpos[1]);
+ cl_ecal_e.get(iSide).fill(ecal_all_clusters.get(icl).getEnergy());
}
- ecal_cl = selectCluster(ecal_cls);
-
- if (debug) System.out.println("Selected " + ecal_cl.toString() +" cluster to use as pointer of fast track ");
+ sel_ecal_idx = selectCluster(ecal_cls,ecal_all_clusters);
+
+ if (debug) System.out.println("Selected clid " + sel_ecal_idx + " is cluster " + ecal_all_clusters.get(sel_ecal_idx).toString() + " and will be used as pointer of fast track ");
+
+ if (sel_ecal_idx < 0) {
+ if (debug) System.out.println("No selected cluster!");
+ continue;
+ }
- int selclpos[] = getCrystalPair(ecal_cl);
- nselcl_ecal_map.get(iSide).fill(selclpos[0], selclpos[1]);
+ selclids = getCrystalPair(ecal_all_clusters.get(sel_ecal_idx).getPosition());
+ nselcl_ecal_map.get(iSide).fill(selclids[0], selclids[1]);
+
+ selcl_ecal_e.get(iSide).fill(ecal_all_clusters.get(sel_ecal_idx).getEnergy());
//Get "target" position i.e. the origin of the radiation
origin_pos = getFastTrackOrigin(event);
@@ -228,9 +283,9 @@
if (debug) System.out.println("Conversion started at " + origin_pos.toString());
//Create the fast track
- fastTrack = new FastTrack(origin_pos,ecal_cl,debug);
-
+ fastTrack.setTrack(origin_pos,ecal_all_clusters.get(sel_ecal_idx).getPosition());
+
if ( debug ) System.out.println(fastTrack.toString());
@@ -241,17 +296,25 @@
if( debug ) System.out.println("There are " + nhitsInTracker + " hits on this side of the tracker");
//most upstream layer nr is different for top and bottom
- int firstLayer = 1;
- if (sides[iSide] == "down") firstLayer=2;
- nhitsInLayer1 = getNAxialHitsInLayers(trackerHits,sides[iSide],firstLayer);
+// int firstLayer = 1;
+// if (sides[iSide] == "down") firstLayer=2;
+// nhitsInLayer1 = getNAxialHitsInLayers(trackerHits,sides[iSide],firstLayer);
- if( debug ) System.out.println("There are " + nhitsInLayer1 + " hits on the first layer for this side (layer=" + firstLayer +")");
+ int nhitsInLayers[] = getNAxialHitsPerLayer(trackerHits,sides[iSide]);
- nhits_tracker.get(iSide).fill(nhitsInTracker);
-
+ if( debug ) {
+ System.out.println("Number of hits per layer:");
+ for (int i=0;i<5;++i){
+ System.out.println("Layer " + (i+1) + " has " + nhitsInLayers[i] + " axial hits");
+ }
+ }
- stripList = new ArrayList<SiTrackerHitStrip1D>();
+ nhits_tracker.get(iSide).fill(nhitsInTracker);
+ for (int i=0;i<5;++i){
+ nhits_tracker_layer.get(iSide).get(i).fill(nhitsInLayers[i]);
+ }
+ //stripList = new ArrayList<SiTrackerHitStrip1D>();
for ( SiTrackerHitStrip1D stripCluster : trackerHits ) {
isaxial = isAxialHit(stripCluster);
@@ -286,6 +349,8 @@
res = fastTrack.getFastTrackResidual(stripCluster);
+ if (debug) System.out.println("-->res " + res);
+
// 1 2 3 4 5 6 7 8 9 10
// 0 1 2 3 4
// 0 1 2 3 4
@@ -298,9 +363,10 @@
if (nhitsInTracker==5) {
resy_org_layallhit.get(iSide).get(layerIndex).fill(res);
+ resy_2d_org_layallhit.get(iSide).fill((double)layerIndex+1, res, 1);
}
- if (nhitsInLayer1>0) {
+ if (nhitsInLayers[0]>0) {
resy_org_lay1hit.get(iSide).get(layerIndex).fill(res);
}
@@ -336,6 +402,9 @@
return pos;
}
+
+
+
// private List<HPSEcalCluster> getAllEcalClusters(EventHeader event) {
//
// List<HPSEcalCluster> clusters = event.get(HPSEcalCluster.class, "EcalReadoutHits");
@@ -350,72 +419,101 @@
private int[] getCrystalPair(Hep3Vector cluster) {
- double x = cluster.x();
- double y = cluster.y();
+ double pos[] = {cluster.x(),cluster.y(),0.};
+ return getCrystalPair(pos);
+ }
+
+ private int[] getCrystalPair(double[] cluster) {
+ double x = cluster[0];
+ double y = cluster[1];
int position[] = new int[2];
position[0] = (int)Math.floor(x/crystalX);
- //position[0] = crystalX * ix;
- position[1] = (int) Math.floor( (y - beamGap * Math.signum(y) ) / crystalY);
- //position[1] = crystalY * iy + beamGap * Math.signum(iy);
- //System.out.println("x " + x + " y " + y + " --> ix " + position[0] + " iy " + position[1]);
- return position;
- //IDDecoder dec = dec.getSubdetector("Ecal").getIDDecoder();
- //dec.setID(hit.getCellID());
- //int ix = dec.getValue("ix");
- //int iy = dec.getValue("iy");
+ position[1] = (int) Math.floor( (y - beamGap * Math.signum(y) ) / crystalY);
+ return position;
}
private double[] getClusterPosition(HPSEcalCluster cluster) {
CalorimeterHit hit = cluster.getSeedHit();
//return cluster.getPosition();
-
+
//IDDecoder dec = dec.getSubdetector("Ecal").getIDDecoder();
- dec.setID(hit.getCellID());
- int ix = dec.getValue("ix");
- int iy = dec.getValue("iy");
- double position[] = new double[2];
- position[0] = crystalX * ix;
- position[1] = crystalY * iy + beamGap * Math.signum(iy);
- return position;
+// dec.setID(hit.getCellID());
+// int ix = dec.getValue("ix");
+// int iy = dec.getValue("iy");
+// double position[] = new double[2];
+// position[0] = crystalX * ix;
+// position[1] = crystalY * iy + beamGap * Math.signum(iy);
+ if (debug) System.out.println("Getting cluster position");
+ double pos[] = cluster.getPosition();
+ if (debug) System.out.println("Found " + pos[0] + "," + pos[1] + "," + pos[2]);
+ return pos;
}
- private List<Hep3Vector> getEcalClustersForFastTracking(List<HPSEcalCluster> clusters, String side) {
- if(!side.equalsIgnoreCase("up") && !side.equalsIgnoreCase("down")) {
+ private List<Integer> getEcalClustersForFastTracking(List<HPSEcalCluster> clusters, String side) {
+ if(side!="up" && side!="down") {
throw new RuntimeException("This ecal side" + side + " do not exist!!");
}
- List<Hep3Vector> cls = new ArrayList<Hep3Vector>();
+ List<Integer> cls = new ArrayList<Integer>();
boolean save;
double [] pos;
- for ( HPSEcalCluster cl : clusters) {
+
+ for (int i=0;i<clusters.size();++i) {
save=false;
- pos = getClusterPosition(cl);
+ pos = getClusterPosition(clusters.get(i));
- if(pos[1]>=0 && side.equalsIgnoreCase("up")) {
+ if(pos[1]>=0 && side=="up") {
save=true;
- } else if(pos[1]<0 && side.equalsIgnoreCase("down")) {
+ }
+ if(pos[1]<0 && side=="down") {
save=true;
}
if(save) {
- cls.add(new BasicHep3Vector(pos[0],pos[1],EcalZPosition));
- //cls.add(new BasicHep3Vector(pos[0],pos[1],pos[2]));
+ //cls.add(new BasicHep3Vector(pos[0],pos[1],EcalZPosition));
+ cls.add(i);
+ //new BasicHep3Vector(pos[0],pos[1],pos[2]));
}
}
return cls;
}
- private Hep3Vector selectCluster(List<Hep3Vector> clusters) {
+ private int selectCluster(List<Integer> ids, List<HPSEcalCluster> clusters) {
//need to decide which cluster to take
- if (clusters.size()>0) {
- Hep3Vector pos = clusters.get(0);
- return pos;
- } else {
- throw new RuntimeException("No cluser positions to choose from!");
+
+ if (clusters.size()==0) {
+ throw new RuntimeException("No clusters to select from!!!");
+ }
+
+ if (ids.size()==0) {
+ throw new RuntimeException("No idx to clusters to select from!!!");
}
+ int sel_id = -1;
+
+ switch (ecalClusterSel) {
+ case 1:
+ //Require at least 3000MeV cluster
+ double E = -1.0;
+ for (int i=0;i<ids.size();++i) {
+ if (clusters.get(ids.get(i)).getEnergy()>3000.0 && clusters.get(ids.get(i)).getEnergy()>E ) {
+ sel_id = ids.get(i);
+ E = clusters.get(ids.get(i)).getEnergy();
+ }
+ }
+ break;
+ default:
+ //Use the first in list
+ if (ids.size()>0 && ids.get(0)<clusters.size()) {
+ sel_id = ids.get(0);
+ }
+ break;
+ }
+
+ return sel_id;
+
}
private List<HPSEcalCluster> getAllEcalClustersForFastTracking(EventHeader event) {
@@ -490,6 +588,48 @@
return nhits;
}
+
+ private int[] getNAxialHitsPerLayer(List<SiTrackerHitStrip1D> trackerHits, String side) {
+ int nhits=0;
+ String si_side;
+ SiSensor siSensor;
+ String name;
+ int l;
+ int i;
+ int n[] = {0,0,0,0,0};
+ for ( SiTrackerHitStrip1D stripCluster : trackerHits ) {
+
+ if(isAxialHit(stripCluster)==false) continue;
+
+ si_side = getSideFromSiCluster(stripCluster);
+
+ if( side == si_side) {
+
+ siSensor = stripCluster.getSensor();
+ name = siSensor.getName();
+ if ( name.length() < 14) {
+ System.err.println("This name is too short!!");
+ throw new RuntimeException("SiSensor name " + name + " is invalid?");
+ }
+
+ l = getLayerFromSensorName(name);
+
+ //Turn this into the layer scheme I use
+ if(side=="down") {
+ //even numbers are used
+ i = (l/2) - 1;
+ n[i] = n[i] + 1;
+ } else {
+ i = ((l+1)/2) - 1;
+ n[i] = n[i] + 1;
+ }
+
+ }
+
+ }
+
+ return n;
+ }
private String getSideFromSiCluster(SiTrackerHitStrip1D stripCluster) {
Hep3Vector posVec = stripCluster.getPositionAsVector();
@@ -531,7 +671,7 @@
IPlotter plotter_org = af.createPlotterFactory().create("HPS SVT Fast Track Residuals (All tracks)");
plotter_org.createRegions(5,2,0);
IPlotter plotter_org_LayAllHit = af.createPlotterFactory().create("HPS SVT Fast Track Residuals (All layers has hits)");
- plotter_org_LayAllHit.createRegions(5,2,0);
+ plotter_org_LayAllHit.createRegions(6,2,0);
IPlotter plotter_org_Lay1Hit = af.createPlotterFactory().create("HPS SVT Fast Track Residuals (1st layer has hits)");
plotter_org_Lay1Hit.createRegions(5,2,0);
@@ -546,30 +686,58 @@
plotter_org.region(idx).plot(resy_org.get(iSide).get(iLayer));
plotter_org_LayAllHit.region(idx).plot(resy_org_layallhit.get(iSide).get(iLayer));
plotter_org_Lay1Hit.region(idx).plot(resy_org_lay1hit.get(iSide).get(iLayer));
-
+
+
}
+ //IPlotterStyle style = plotter_org_LayAllHit.style();
+ //style.regionBoxStyle().
+
+ plotter_org_LayAllHit.region(10+iSide).style().setParameter("hist2DStyle", "colorMap");
+ plotter_org_LayAllHit.region(10+iSide).style().setParameter("colorMapScheme", "rainbow");
+ plotter_org_LayAllHit.region(10+iSide).plot(resy_2d_org_layallhit.get(iSide));
}
plotter_org.show();
plotter_org_LayAllHit.show();
plotter_org_Lay1Hit.show();
- //Hit mupltiplicity
+
+
+ //Hit multiplicity
IPlotter plotter_hitmult = af.createPlotterFactory().create();
plotter_hitmult.createRegions(1,2,0);
plotter_hitmult.region(0).plot(nhits_tracker.get(0));
plotter_hitmult.region(1).plot(nhits_tracker.get(1));
plotter_hitmult.show();
+ IPlotter plotter_hitmult_layer = af.createPlotterFactory().create();
+ plotter_hitmult_layer.createRegions(5,2,0);
+ for ( int iSide=0;iSide<2;++iSide) {
+ int nLayers = resy_org.get(iSide).size();
+ for (int iLayer=0;iLayer<nLayers;++iLayer) {
+ int idx = 2*iLayer+iSide;
+
+ //0 1 2 3 4
+ //1 3 5 7 9
+ plotter_hitmult_layer.region(idx).plot(nhits_tracker_layer.get(iSide).get(iLayer));
+
+ }
+ }
+ plotter_hitmult_layer.show();
+
//ECal plots
IPlotter plotter_ecalhitmult = af.createPlotterFactory().create();
- plotter_ecalhitmult.createRegions(3,2,0);
+ plotter_ecalhitmult.createRegions(5,2,0);
plotter_ecalhitmult.region(0).plot(ncl_ecal.get(0));
plotter_ecalhitmult.region(1).plot(ncl_ecal.get(1));
plotter_ecalhitmult.region(2).plot(ncl_ecal_map.get(0));
plotter_ecalhitmult.region(3).plot(ncl_ecal_map.get(1));
plotter_ecalhitmult.region(4).plot(nselcl_ecal_map.get(0));
plotter_ecalhitmult.region(5).plot(nselcl_ecal_map.get(1));
+ plotter_ecalhitmult.region(6).plot(cl_ecal_e.get(0));
+ plotter_ecalhitmult.region(7).plot(cl_ecal_e.get(1));
+ plotter_ecalhitmult.region(8).plot(selcl_ecal_e.get(0));
+ plotter_ecalhitmult.region(9).plot(selcl_ecal_e.get(1));
plotter_ecalhitmult.show();