Commit in lcsim/src/org/lcsim/rpc on MAIN | |||
PadDigiDriver.java | +568 | added 1.1 | |
RPCDigiDriver.java | -777 | 1.1 removed | |
+568 | -777 |
rename
diff -N PadDigiDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ PadDigiDriver.java 4 Feb 2013 22:14:39 -0000 1.1 @@ -0,0 +1,568 @@
+package org.lcsim.rpc; + +import hep.physics.vec.BasicHep3Vector; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; + +import java.util.ArrayList; +import java.util.List; + +import org.lcsim.event.CalorimeterHit; +import org.lcsim.event.EventHeader; +import org.lcsim.event.EventHeader.LCMetaData; +import org.lcsim.event.SimCalorimeterHit; +import org.lcsim.geometry.IDDecoder; +import org.lcsim.util.Driver; +import org.lcsim.util.lcio.LCIOConstants; + +// FIXME: The code to for charge deposition on barrel and endcap is essentially copied and pasted. +// This should be one method called with parameters. +// TODO: Represent output data as RawCalorimeterHit possible. Charge should be converted to int from double. +// TODO: Move histograms to separate Driver. +// TODO: Add settable input parameters: +// -q0 (what is it???) +public class PadDigiDriver extends Driver +{ + // Charge distribution data. + ChargeDistributionTables chargeTables; + + // charge cut + // units??? + double minimumCharge = 0.400; + + // distance cut in mm + double distanceCut = 0.90; + + // Collection names. + private String barrelInputCollection = "HcalBarrelHits"; + private String endcapInputCollection = "HcalEndcapHits"; + private static String defaultOutputTag = "_pad_digi"; + private String barrelOutputCollection = barrelInputCollection + defaultOutputTag; + private String endcapOutputCollection = endcapInputCollection + defaultOutputTag; + + public void setMinimumCharge(double minimumCharge) + { + this.minimumCharge = minimumCharge; + } + + public void setDistanceCut(double distanceCut) + { + this.distanceCut = distanceCut; + } + + public void setBarrelInputCollection(String barrelInputCollection) + { + this.barrelInputCollection = barrelInputCollection; + } + + public void setEndcapInputCollection(String endcapInputCollection) + { + this.endcapInputCollection = endcapInputCollection; + } + + public void setBarrelOutputCollection(String barrelOutputCollection) + { + this.barrelOutputCollection = barrelOutputCollection; + } + + public void setEndcapOutputCollection(String endcapOutputCollection) + { + this.endcapOutputCollection = endcapOutputCollection; + } + + public void setup() throws Exception + { + // Initialize charge distrubtion tables. + chargeTables = new ChargeDistributionTables(); + } + + public void startOfData() + { + try { + setup(); + } catch (Exception e) { + throw new RuntimeException(e); + } + } + + public void process(EventHeader event) + { + // parameter for slope = 0.084 + // FIXME: What are these parameters? Should be settable on Driver? + double q0 = -0.062; // ????? + + List<SimCalorimeterHit> barrelHits = event.get(SimCalorimeterHit.class, barrelInputCollection); + System.out.println("number of hits in barrel = " + barrelHits.size()); + + List<SimCalorimeterHit> endcapHits = event.get(SimCalorimeterHit.class, endcapInputCollection); + System.out.println("number of hits in endcap = " + endcapHits.size()); + + IDDecoder decoder = event.getMetaData(barrelHits).getIDDecoder(); + + /* + * // id = 0 doesn't exist --> id=0 can be used to indicate no cell found decoder.setID(0); + * Hep3Vector Pos = decoder.getPositionVector(); + * System.out.println("position for id(0) = "+Pos); + */ + long[] ordered_ids = new long[25]; + + ArrayList<Hep3Vector> used_points = new ArrayList<Hep3Vector>(); + int count_points = 0; + + ArrayList<RPCHit> br_rpchits = new ArrayList<RPCHit>(); // hit storage for barrel hits after rpc + // simulation + ArrayList<RPCHit> ec_rpchits = new ArrayList<RPCHit>(); // hit storage for endcap hits after rpc + // simulation + + int count = 0; + for (SimCalorimeterHit hit : barrelHits) { + decoder.setID(hit.getCellID()); + Hep3Vector cellPos = decoder.getPositionVector(); + int ncontribs = hit.getMCParticleCount(); + + long[] nb_ids = decoder.getNeighbourIDs(0, 2, 2); + + for (int i = 0; i < 25; i++) { + ordered_ids[i] = 0; // make sure to have to way to tell if a neighboring cell is not + // available (id = 0) + } + + // these are two independent, orthogonal vector bases in the readout pads plane + Hep3Vector vx = new BasicHep3Vector(); + Hep3Vector vy = new BasicHep3Vector(); + + // here are the two compoments of the local point position within a cell, according to the + // directions defined by the above bases + double pdx = 0.; + double pdy = 0.; + + long[] nb_final_ids = new long[24]; + int cnt = 0; + if (nb_ids.length > 0) // for all neighbors, need to remove cells that are not on the same + // detector sector (and thus, not on the same RPC) + { + ArrayList<Hep3Vector> bases = new ArrayList<Hep3Vector>(); + + cnt = 0; + for (int i = 0; i < nb_ids.length; i++) { + decoder.setID(nb_ids[i]); + Hep3Vector cellPos_i = decoder.getPositionVector(); + Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); + double dist = Math.sqrt((cellPos.x() - cellPos_i.x()) * (cellPos.x() - cellPos_i.x()) + (cellPos.y() - cellPos_i.y()) * (cellPos.y() - cellPos_i.y()) + (cellPos.z() - cellPos_i.z()) * (cellPos.z() - cellPos_i.z())); + + if ((Math.abs(dist - 10.) < 0.0001) || (Math.abs(dist - 14.142136) < 0.0001) || (Math.abs(dist - 20.) < 0.0001) || (Math.abs(dist - 22.36068) < 0.0001) || (Math.abs(dist - 28.284271) < 0.0001)) { + nb_final_ids[cnt] = nb_ids[i]; + cnt++; + } else { + // System.out.println(" dist = "+dist+" cnt = "+cnt); + } + + if (Math.abs(dist - 10.) < 0.0001) + bases.add(vi0); + + } + + // System.out.println(" bases.size() = "+bases.size()); + if (bases.size() < 3) + System.out.println(" problem finding enough vector bases, bases.size() = " + bases.size()); + vx = bases.get(0); + vy = bases.get(1); + if ((Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) - 10.) < 0.0001) + vy = bases.get(2); + + if (Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) > 0.0001) + System.out.println(" vx*vy (barrel) = " + ((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))); + + for (int i = 0; i < cnt; i++) { + decoder.setID(nb_final_ids[i]); + Hep3Vector cellPos_i = decoder.getPositionVector(); + Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); + + double dx = ((vx.x()) * (vi0.x()) + (vx.y()) * (vi0.y()) + (vx.z()) * (vi0.z())) / 10.; + double dy = ((vy.x()) * (vi0.x()) + (vy.y()) * (vi0.y()) + (vy.z()) * (vi0.z())) / 10.; + int idx = (int) Math.floor(dx / 10. + 0.00001); + int idy = (int) Math.floor(dy / 10. + 0.00001); + int ind = idx + 2 + (idy + 2) * 5; + + ordered_ids[ind] = nb_final_ids[i]; + } + ordered_ids[12] = hit.getCellID(); // beside neighboring cells, also need to set id for + // the hit cell itself. + + } + + count++; + // if (count < 10000) { + if (count > -1) { + // if (event.getEventNumber() == 1) System.out.println("barrel hit number: "+count); + + for (int i = 0; i < ncontribs; i++) { + float stepPos[] = hit.getStepPosition(i); + Hep3Vector stepVec = new BasicHep3Vector(stepPos[0], stepPos[1], stepPos[2]); + + int ifound = 0; + count_points++; + if (used_points.size() > 0) { + for (Hep3Vector ivec : used_points) { + Hep3Vector tmpv = VecOp.sub(ivec, stepVec); + if (Math.sqrt(VecOp.dot(tmpv, tmpv)) < distanceCut) + ifound = 1; + } + + if (ifound == 0) + used_points.add(stepVec); + } else { + used_points.add(stepVec); + } + + if (ifound == 0) // only deal with points that survived dcut + { + Hep3Vector pdr = VecOp.sub(stepVec, cellPos); + pdx = VecOp.dot(pdr, vx) / 10. + 5.; + pdy = VecOp.dot(pdr, vy) / 10. + 5.; + + if ((pdx < 0.) || (pdx > 10.) || (pdy < 0.) || (pdy > 10.)) + System.out.println("out of bound point position!!!!!!!!!!!!!! barrel pdx = " + pdx + " pdy = " + pdy); + + double rnd = Math.random(); // generate total charge here + if (rnd == 1.) + rnd = rnd - 0.0000001; + int i1 = (int) Math.floor(rnd * 1000.); + int i2 = i1 + 1; + + double icharge = 0.; + if (Math.abs(i2 * 1.0 - rnd * 1000.) < 0.00001) { + //icharge = tabCharge[i2]; + icharge = chargeTables.getCharge(i2); + } else { + double k = (rnd * 1000. - i1 * 1.) / (i2 * 1. - rnd * 1000.); + icharge = (chargeTables.getCharge(i1) + chargeTables.getCharge(i2) * k) / (1 + k); + } + + icharge = icharge / 1100. - q0; // done generating total charge + + pdx = pdx / 10.; // change pdx into cm unit + pdy = pdy / 10.; // change pdx into cm unit + + double[] charge = new double[25]; // start charge distribution calulation + int ix1 = (int) Math.floor(pdx / 0.5 * 100); + int iy1 = (int) Math.floor(pdy / 0.5 * 100); + if (ix1 == 200) + ix1 = 199; + if (iy1 == 200) + iy1 = 199; + int ix2 = ix1 + 1; + int iy2 = iy1 + 1; + + double mid1 = 0.; + double mid2 = 0.; + for (int k = 0; k < 25; k++) { + if (Math.abs(ix2 * 1.0 - pdx * 200.) < 0.00001) { + mid1 = chargeTables.getPadCharge(ix2, iy1, k); + mid2 = chargeTables.getPadCharge(ix2, iy2, k); + } else { + double k1 = (pdx * 200. - ix1 * 1.) / (ix2 * 1. - pdx * 200.); + mid1 = (chargeTables.getPadCharge(ix1, iy1, k) + chargeTables.getPadCharge(ix2, iy1, k) * k1) / (1 + k1); + mid2 = (chargeTables.getPadCharge(ix1, iy2, k) + chargeTables.getPadCharge(ix2, iy2, k) * k1) / (1 + k1); + } + + if (Math.abs(iy2 * 1.0 - pdy * 200.) < 0.00001) { + charge[k] = icharge * mid2; + } else { + double k2 = (pdy * 200. - iy1 * 1.) / (iy2 * 1. - pdy * 200.); + charge[k] = icharge * (mid1 + mid2 * k2) / (1 + k2); + } + + // if (event.getEventNumber() < 1) + // System.out.println("barrel charge["+k+"] = "+charge[k]); + + if ((Math.abs(ordered_ids[k]) > 0) && (charge[k] > 0)) // pad exist on the + // neighber or hit pad + // list, and there's + // charge to deposit + { + if (br_rpchits.size() == 0) // first hit, just add + { + RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); + br_rpchits.add(tmphit); + } else // need to see if the pad is already there + { + int hfound = 0; + int j = 0; + while ((hfound == 0) && (j < br_rpchits.size())) { + if (br_rpchits.get(j).getCid() == ordered_ids[k]) { + hfound = 1; + br_rpchits.get(j).addCharge(charge[k]); + } + + j++; + } + + if (hfound == 0) { + RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); + br_rpchits.add(tmphit); + } + } + } + } // done with charge distribution calulation + + } + } + } + } + + used_points.clear(); + count_points = 0; + + decoder = event.getMetaData(endcapHits).getIDDecoder(); + + count = 0; + for (SimCalorimeterHit hit : endcapHits) { + decoder.setID(hit.getCellID()); + Hep3Vector cellPos = decoder.getPositionVector(); + int ncontribs = hit.getMCParticleCount(); + + long[] nb_ids = decoder.getNeighbourIDs(0, 2, 2); + + for (int i = 0; i < 25; i++) { + ordered_ids[i] = 0; // make sure to have to way to tell if a neighboring cell is not + // available (id = 0) + } + + // these are two independent, orthogonal vector bases in the readout pads plane + Hep3Vector vx = new BasicHep3Vector(); + Hep3Vector vy = new BasicHep3Vector(); + + // here are the two compoments of the local point position within a cell, according to the + // directions defined by the above bases + double pdx = 0.; + double pdy = 0.; + + long[] nb_final_ids = new long[24]; + int cnt = 0; + if (nb_ids.length > 0) // for all neighbors, need to remove cells that are not on the same + // detector sector (and thus, not on the same RPC) + { + ArrayList<Hep3Vector> bases = new ArrayList<Hep3Vector>(); + + cnt = 0; + for (int i = 0; i < nb_ids.length; i++) { + decoder.setID(nb_ids[i]); + Hep3Vector cellPos_i = decoder.getPositionVector(); + Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); + double dist = Math.sqrt((cellPos.x() - cellPos_i.x()) * (cellPos.x() - cellPos_i.x()) + (cellPos.y() - cellPos_i.y()) * (cellPos.y() - cellPos_i.y()) + (cellPos.z() - cellPos_i.z()) * (cellPos.z() - cellPos_i.z())); + + if ((Math.abs(dist - 10.) < 0.0001) || (Math.abs(dist - 14.142136) < 0.0001) || (Math.abs(dist - 20.) < 0.0001) || (Math.abs(dist - 22.36068) < 0.0001) || (Math.abs(dist - 28.284271) < 0.0001)) { + nb_final_ids[cnt] = nb_ids[i]; + cnt++; + } else { + // System.out.println(" dist = "+dist+" cnt = "+cnt); + } + + if (Math.abs(dist - 10.) < 0.0001) + bases.add(vi0); + + } + + if (bases.size() < 3) + System.out.println(" problem finding enough vector bases, bases.size() = " + bases.size()); + vx = bases.get(0); + vy = bases.get(1); + if ((Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) - 10.) < 0.0001) + vy = bases.get(2); + + if (Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) > 0.0001) + System.out.println(" vx*vy (endcap) = " + ((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))); + + for (int i = 0; i < cnt; i++) { + decoder.setID(nb_final_ids[i]); + Hep3Vector cellPos_i = decoder.getPositionVector(); + Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); + + double dx = ((vx.x()) * (vi0.x()) + (vx.y()) * (vi0.y()) + (vx.z()) * (vi0.z())) / 10.; + double dy = ((vy.x()) * (vi0.x()) + (vy.y()) * (vi0.y()) + (vy.z()) * (vi0.z())) / 10.; + int idx = (int) Math.floor(dx / 10. + 0.00001); + int idy = (int) Math.floor(dy / 10. + 0.00001); + int ind = idx + 2 + (idy + 2) * 5; + + ordered_ids[ind] = nb_final_ids[i]; + } + ordered_ids[12] = hit.getCellID(); // beside neighboring cells, also need to set id for + // the hit cell itself. + + } + + count++; + // if (count < 10000) { + if (count > -1) { + // System.out.println("endcap hit number: "+count); + + for (int i = 0; i < ncontribs; i++) { + float stepPos[] = hit.getStepPosition(i); + Hep3Vector stepVec = new BasicHep3Vector(stepPos[0], stepPos[1], stepPos[2]); + + int ifound = 0; + count_points++; + if (used_points.size() > 0) { + for (Hep3Vector ivec : used_points) { + // if + // (Math.sqrt((ivec.x()-stepVec.x())*(ivec.x()-stepVec.x())+(ivec.y()-stepVec.y())*(ivec.y()-stepVec.y())+(ivec.z()-stepVec.z())*(ivec.z()-stepVec.z())) + // < dcut) ifound = 1; + Hep3Vector tmpv = VecOp.sub(ivec, stepVec); + if (Math.sqrt(VecOp.dot(tmpv, tmpv)) < distanceCut) + ifound = 1; + } + + if (ifound == 0) + used_points.add(stepVec); + } else { + used_points.add(stepVec); + } + + if (ifound == 0) // only deal with points that survived dcut + { + Hep3Vector pdr = VecOp.sub(stepVec, cellPos); + pdx = VecOp.dot(pdr, vx) / 10. + 5.; + pdy = VecOp.dot(pdr, vy) / 10. + 5.; + + if ((pdx < 0.) || (pdx > 10.) || (pdy < 0.) || (pdy > 10.)) + System.out.println("out of bound point position; endcap pdx = " + pdx + " pdy = " + pdy); + + double rnd = Math.random(); // generate total charge here + if (rnd == 1.) + rnd = rnd - 0.0000001; + int i1 = (int) Math.floor(rnd * 1000.); + int i2 = i1 + 1; + + double icharge = 0.; + if (Math.abs(i2 * 1.0 - rnd * 1000.) < 0.00001) { + icharge = chargeTables.getCharge(i2); + } else { + double k = (rnd * 1000. - i1 * 1.) / (i2 * 1. - rnd * 1000.); + icharge = (chargeTables.getCharge(i1) + chargeTables.getCharge(i2) * k) / (1 + k); + } + + icharge = icharge / 1100. - q0; // done generating total charge + + pdx = pdx / 10.; // change pdx into cm unit + pdy = pdy / 10.; // change pdx into cm unit + + double[] charge = new double[25]; // start charge distribution calulation + int ix1 = (int) Math.floor(pdx / 0.5 * 100); + int iy1 = (int) Math.floor(pdy / 0.5 * 100); + if (ix1 == 200) + ix1 = 199; + if (iy1 == 200) + iy1 = 199; + int ix2 = ix1 + 1; + int iy2 = iy1 + 1; + + double mid1 = 0.; + double mid2 = 0.; + for (int k = 0; k < 25; k++) { + if (Math.abs(ix2 * 1.0 - pdx * 200.) < 0.00001) { + mid1 = chargeTables.getPadCharge(ix2, iy1, k); + mid2 = chargeTables.getPadCharge(ix2, iy2, k); + } else { + double k1 = (pdx * 200. - ix1 * 1.) / (ix2 * 1. - pdx * 200.); + mid1 = (chargeTables.getPadCharge(ix1, iy1, k) + chargeTables.getPadCharge(ix2, iy1, k) * k1) / (1 + k1); + mid2 = (chargeTables.getPadCharge(ix1, iy2, k) + chargeTables.getPadCharge(ix2, iy2, k) * k1) / (1 + k1); + } + + if (Math.abs(iy2 * 1.0 - pdy * 200.) < 0.00001) { + charge[k] = icharge * mid2; + } else { + double k2 = (pdy * 200. - iy1 * 1.) / (iy2 * 1. - pdy * 200.); + charge[k] = icharge * (mid1 + mid2 * k2) / (1 + k2); + } + + // if (event.getEventNumber() < 3) + // System.out.println("endcap charge["+k+"] = "+charge[k]); + + if ((Math.abs(ordered_ids[k]) > 0) && (charge[k] > 0)) // pad exist on the + // neighber or hit pad + // list, and there's + // charge to deposit + { + if (ec_rpchits.size() == 0) // first hit, just add + { + RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); + ec_rpchits.add(tmphit); + } else // need to see if the pad is already there + { + int hfound = 0; + int j = 0; + while ((hfound == 0) && (j < ec_rpchits.size())) { + if (ec_rpchits.get(j).getCid() == ordered_ids[k]) { + hfound = 1; + ec_rpchits.get(j).addCharge(charge[k]); + } + + j++; + } + + if (hfound == 0) { + RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); + ec_rpchits.add(tmphit); + } + } + } + + } // done with charge distribution calulation + + } + } + } + } + + used_points.clear(); + count_points = 0; + + //System.out.println("before thr cut, barrel cells with charge = " + br_rpchits.size()); + //System.out.println("before thr cut, endcap cells with charge = " + ec_rpchits.size()); + + ArrayList<RPCHit> br_final = new ArrayList<RPCHit>(); // barrel hits after rpc simulation + ArrayList<RPCHit> ec_final = new ArrayList<RPCHit>(); // endcap hits after rpc simulation + + decoder = event.getMetaData(barrelHits).getIDDecoder(); + for (int i = 0; i < br_rpchits.size(); i++) { + if (br_rpchits.get(i).getCharge() > minimumCharge) { + br_final.add(br_rpchits.get(i)); + } + } + + decoder = event.getMetaData(endcapHits).getIDDecoder(); + for (int i = 0; i < ec_rpchits.size(); i++) { + if (ec_rpchits.get(i).getCharge() > minimumCharge) { + ec_final.add(ec_rpchits.get(i)); + } + } + + //br_rpchits.clear(); + //br_final.clear(); + //ec_rpchits.clear(); + //ec_final.clear(); + + int flag = 0; + flag = flag | (1 << LCIOConstants.CHBIT_LONG); + flag = flag | (1 << LCIOConstants.RCHBIT_ID1); + + // Store final collection of barrel hits in event. + LCMetaData barrelMetaData = event.getMetaData(barrelHits); + event.put( + barrelOutputCollection, + br_final, + CalorimeterHit.class, + flag, + barrelMetaData.getName()); + + // Store final collection of endcap hits in event. + LCMetaData endcapMetaData = event.getMetaData(endcapHits); + event.put( + endcapOutputCollection, + ec_final, + CalorimeterHit.class, + flag, + endcapMetaData.getName()); + } +}
\ No newline at end of file
diff -N RPCDigiDriver.java --- RPCDigiDriver.java 1 Feb 2013 01:54:56 -0000 1.1 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,777 +0,0 @@
-package org.lcsim.rpc; - -import hep.aida.IHistogram1D; -import hep.aida.IHistogram2D; -import hep.physics.vec.BasicHep3Vector; -import hep.physics.vec.Hep3Vector; -import hep.physics.vec.VecOp; - -import java.io.IOException; -import java.util.ArrayList; -import java.util.List; - -import org.lcsim.event.CalorimeterHit; -import org.lcsim.event.EventHeader; -import org.lcsim.event.EventHeader.LCMetaData; -import org.lcsim.event.SimCalorimeterHit; -import org.lcsim.geometry.IDDecoder; -import org.lcsim.util.Driver; -import org.lcsim.util.aida.AIDA; -import org.lcsim.util.lcio.LCIOConstants; - -// FIXME: The code to for charge deposition on barrel and endcap is essentially copied and pasted. -// This should be one method called with parameters. -// TODO: Represent output data as RawCalorimeterHit possible. Charge should be converted to int from double. -// TODO: Move histograms to separate class or Driver. -// TODO: Add settable input parameters: -// -input barrel collection -// -input endcap collection -// -output barrel collection -// -output endcap collection -// -distance cut -// -charge cut -// -q0 (what is it???) -public class RPCDigiDriver extends Driver -{ - // Charge distribution data. - ChargeDistributionTables chargeTables; - - // ??? - int count_b = 0; - - // Collection names. - private String hcalBarrelCollectionName = "HcalBarrelHits"; - private String hcalEndcapCollectionName = "HcalEndcapHits"; - private static String collectionOutputTag = "_pad_digi"; - - // True to fill histograms; false to skip. - private static final boolean fillHistograms = true; - - // Declare histograms: - - IHistogram1D c1_cell_x; - IHistogram1D c1_cell_y; - IHistogram1D c1_cell_z; - IHistogram1D c1_step_x; - IHistogram1D c1_step_y; - IHistogram1D c1_step_z; - - IHistogram1D c1_nb_ct_br; - IHistogram1D c1_nb_ct_ec; - - IHistogram1D c1_nb_ct_br_final; - IHistogram1D c1_nb_ct_ec_final; - - IHistogram2D c2_cell_xy_br; - IHistogram2D c2_cell_xy_ec; - IHistogram2D c2_cell_xy_ec_nnb; - - IHistogram2D[] c2_xz_cell; - IHistogram2D[] c2_yz_cell; - - IHistogram2D[] c2_xz_step; - IHistogram2D[] c2_yz_step; - - /* - public void setHcalBarrelCollectionName(String hcalBarrelCollectionName) - { - this.hcalBarrelCollectionName = hcalBarrelCollectionName; - } - - public void setHcalEndcapCollectionName(String hcalEndcapCollectionName) - { - this.hcalEndcapCollectionName = hcalEndcapCollectionName; - } - */ - - public void setup() throws Exception - { - // Initialize charge distrubtion tables. - chargeTables = new ChargeDistributionTables(); - - // Define histograms. - defineHistograms(); - } - - public void defineHistograms() throws IOException - { - AIDA aida = AIDA.defaultInstance(); - - c1_cell_x = aida.histogram1D("cell pos x", 200, -5000., 5000.); - c1_cell_y = aida.histogram1D("cell pos y", 200, -5000., 5000.); - c1_cell_z = aida.histogram1D("cell pos z", 200, -5000., 5000.); - c1_step_x = aida.histogram1D("step pos x", 200, -8., 8.); - c1_step_y = aida.histogram1D("step pos y", 200, -8., 8.); - c1_step_z = aida.histogram1D("step pos z", 200, -8., 8.); - - c1_nb_ct_br = aida.histogram1D("number of neighbour count: barrel", 50, -0.5, 49.5); - c1_nb_ct_ec = aida.histogram1D("number of neighbour count: endcap", 50, -0.5, 49.5); - - c1_nb_ct_br_final = aida.histogram1D("number of neighbour count (final): barrel", 50, -0.5, 49.5); - c1_nb_ct_ec_final = aida.histogram1D("number of neighbour count (final): endcap", 50, -0.5, 49.5); - - c2_cell_xy_br = aida.histogram2D("br hits with too many nb: xy distribution", 200, -3000., 3000., 200, -3000., 3000.); - c2_cell_xy_ec = aida.histogram2D("ec hits with too many nb: xy distribution", 200, -3000., 3000., 200, -3000., 3000.); - c2_cell_xy_ec_nnb = aida.histogram2D("ec hits with no nb: xy distribution", 200, -3000., 3000., 200, -3000., 3000.); - - c2_xz_cell = new IHistogram2D[10]; - c2_yz_cell = new IHistogram2D[10]; - - c2_xz_step = new IHistogram2D[10]; - c2_yz_step = new IHistogram2D[10]; - - c2_xz_cell[0] = aida.histogram2D("event 0: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[1] = aida.histogram2D("event 1: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[2] = aida.histogram2D("event 2: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[3] = aida.histogram2D("event 3: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[4] = aida.histogram2D("event 4: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[5] = aida.histogram2D("event 5: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[6] = aida.histogram2D("event 6: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[7] = aida.histogram2D("event 7: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[8] = aida.histogram2D("event 8: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_cell[9] = aida.histogram2D("event 9: xz_cell", 200, -3500., 3500., 200, -3500., 3500.); - - c2_yz_cell[0] = aida.histogram2D("event 0: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[1] = aida.histogram2D("event 1: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[2] = aida.histogram2D("event 2: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[3] = aida.histogram2D("event 3: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[4] = aida.histogram2D("event 4: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[5] = aida.histogram2D("event 5: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[6] = aida.histogram2D("event 6: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[7] = aida.histogram2D("event 7: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[8] = aida.histogram2D("event 8: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_cell[9] = aida.histogram2D("event 9: yz_cell", 200, -3500., 3500., 200, -3500., 3500.); - - c2_xz_step[0] = aida.histogram2D("event 0: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[1] = aida.histogram2D("event 1: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[2] = aida.histogram2D("event 2: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[3] = aida.histogram2D("event 3: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[4] = aida.histogram2D("event 4: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[5] = aida.histogram2D("event 5: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[6] = aida.histogram2D("event 6: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[7] = aida.histogram2D("event 7: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[8] = aida.histogram2D("event 8: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_xz_step[9] = aida.histogram2D("event 9: xz_step", 200, -3500., 3500., 200, -3500., 3500.); - - c2_yz_step[0] = aida.histogram2D("event 0: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[1] = aida.histogram2D("event 1: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[2] = aida.histogram2D("event 2: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[3] = aida.histogram2D("event 3: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[4] = aida.histogram2D("event 4: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[5] = aida.histogram2D("event 5: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[6] = aida.histogram2D("event 6: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[7] = aida.histogram2D("event 7: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[8] = aida.histogram2D("event 8: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - c2_yz_step[9] = aida.histogram2D("event 9: yz_step", 200, -3500., 3500., 200, -3500., 3500.); - } - - public void startOfData() - { - try { - setup(); - } catch (Exception e) { - throw new RuntimeException(e); - } - } - - public void process(EventHeader event) - { - // parameter for slope = 0.084 - // FIXME: What are these parameters? Should be settable on Driver? - double q0 = -0.062; // ????? - double thr = 0.400; // charge cut - double dcut = 0.90; // in mm unit - - List<SimCalorimeterHit> barrelHits = event.get(SimCalorimeterHit.class, hcalBarrelCollectionName); - System.out.println("number of hits in barrel = " + barrelHits.size()); - - List<SimCalorimeterHit> endcapHits = event.get(SimCalorimeterHit.class, hcalEndcapCollectionName); - System.out.println("number of hits in endcap = " + endcapHits.size()); - - IDDecoder decoder = event.getMetaData(barrelHits).getIDDecoder(); - - /* - * // id = 0 doesn't exist --> id=0 can be used to indicate no cell found decoder.setID(0); Hep3Vector - * Pos = decoder.getPositionVector(); System.out.println("position for id(0) = "+Pos); - */ - long[] ordered_ids = new long[25]; - - ArrayList<Hep3Vector> used_points = new ArrayList<Hep3Vector>(); - int count_points = 0; - - ArrayList<RPCHit> br_rpchits = new ArrayList<RPCHit>(); // hit storage for barrel hits after rpc - // simulation - ArrayList<RPCHit> ec_rpchits = new ArrayList<RPCHit>(); // hit storage for endcap hits after rpc - // simulation - - int count = 0; - for (SimCalorimeterHit hit : barrelHits) { - decoder.setID(hit.getCellID()); - Hep3Vector cellPos = decoder.getPositionVector(); - int ncontribs = hit.getMCParticleCount(); - - long[] nb_ids = decoder.getNeighbourIDs(0, 2, 2); - if (fillHistograms) - c1_nb_ct_br.fill(nb_ids.length); - - for (int i = 0; i < 25; i++) { - ordered_ids[i] = 0; // make sure to have to way to tell if a neighboring cell is not - // available (id = 0) - } - - // these are two independent, orthogonal vector bases in the readout pads plane - Hep3Vector vx = new BasicHep3Vector(); - Hep3Vector vy = new BasicHep3Vector(); - - // here are the two compoments of the local point position within a cell, according to the - // directions defined by the above bases - double pdx = 0.; - double pdy = 0.; - - long[] nb_final_ids = new long[24]; - int cnt = 0; - if (nb_ids.length > 0) // for all neighbors, need to remove cells that are not on the same - // detector sector (and thus, not on the same RPC) - { - ArrayList<Hep3Vector> bases = new ArrayList<Hep3Vector>(); - - cnt = 0; - for (int i = 0; i < nb_ids.length; i++) { - decoder.setID(nb_ids[i]); - Hep3Vector cellPos_i = decoder.getPositionVector(); - Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); - double dist = Math.sqrt((cellPos.x() - cellPos_i.x()) * (cellPos.x() - cellPos_i.x()) + (cellPos.y() - cellPos_i.y()) * (cellPos.y() - cellPos_i.y()) + (cellPos.z() - cellPos_i.z()) * (cellPos.z() - cellPos_i.z())); - - if ((Math.abs(dist - 10.) < 0.0001) || (Math.abs(dist - 14.142136) < 0.0001) || (Math.abs(dist - 20.) < 0.0001) || (Math.abs(dist - 22.36068) < 0.0001) || (Math.abs(dist - 28.284271) < 0.0001)) { - nb_final_ids[cnt] = nb_ids[i]; - cnt++; - } else { - // System.out.println(" dist = "+dist+" cnt = "+cnt); - } - - if (Math.abs(dist - 10.) < 0.0001) - bases.add(vi0); - - } - - // System.out.println(" bases.size() = "+bases.size()); - if (bases.size() < 3) - System.out.println(" problem finding enough vector bases, bases.size() = " + bases.size()); - vx = bases.get(0); - vy = bases.get(1); - if ((Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) - 10.) < 0.0001) - vy = bases.get(2); - /* - * for (int i = 0; i < bases.size(); i++) { - * System.out.println(" vect "+i+" = "+bases.get(i)); } System.out.println(" vx = "+vx); - * System.out.println(" vy = "+vy); - */ - if (Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) > 0.0001) - System.out.println(" vx*vy (barrel) = " + ((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))); - - for (int i = 0; i < cnt; i++) { - decoder.setID(nb_final_ids[i]); - Hep3Vector cellPos_i = decoder.getPositionVector(); - Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); - - double dx = ((vx.x()) * (vi0.x()) + (vx.y()) * (vi0.y()) + (vx.z()) * (vi0.z())) / 10.; - double dy = ((vy.x()) * (vi0.x()) + (vy.y()) * (vi0.y()) + (vy.z()) * (vi0.z())) / 10.; - int idx = (int) Math.floor(dx / 10. + 0.00001); - int idy = (int) Math.floor(dy / 10. + 0.00001); - int ind = idx + 2 + (idy + 2) * 5; - - // if (event.getEventNumber() == 1) - // System.out.println("barrel i = "+i+" dx = "+dx+" dy = "+dy+" idx = "+idx+" idy = "+idy+" ind = "+ind); - - ordered_ids[ind] = nb_final_ids[i]; - } - ordered_ids[12] = hit.getCellID(); // beside neighboring cells, also need to set id for - // the hit cell itself. - - } - if (fillHistograms) - c1_nb_ct_br_final.fill(cnt); - - count++; - // if (count < 10000) { - if (count > -1) { - // if (event.getEventNumber() == 1) System.out.println("barrel hit number: "+count); - - for (int i = 0; i < ncontribs; i++) { - float stepPos[] = hit.getStepPosition(i); - Hep3Vector stepVec = new BasicHep3Vector(stepPos[0], stepPos[1], stepPos[2]); - - int ifound = 0; - count_points++; - if (used_points.size() > 0) { - for (Hep3Vector ivec : used_points) { - // if - // (Math.sqrt((ivec.x()-stepVec.x())*(ivec.x()-stepVec.x())+(ivec.y()-stepVec.y())*(ivec.y()-stepVec.y())+(ivec.z()-stepVec.z())*(ivec.z()-stepVec.z())) - // < dcut) ifound = 1; - Hep3Vector tmpv = VecOp.sub(ivec, stepVec); - if (Math.sqrt(VecOp.dot(tmpv, tmpv)) < dcut) - ifound = 1; - // if (event.getEventNumber() == 1) - // System.out.println(" sqrt(VecOp.dot(tmpv, tmpv)) = "+Math.sqrt(VecOp.dot(tmpv, - // tmpv))+" ivec = "+ivec+" stepVec = "+stepVec); - } - - if (ifound == 0) - used_points.add(stepVec); - } else { - used_points.add(stepVec); - } - - if (ifound == 0) // only deal with points that survived dcut - { - Hep3Vector pdr = VecOp.sub(stepVec, cellPos); - pdx = VecOp.dot(pdr, vx) / 10. + 5.; - pdy = VecOp.dot(pdr, vy) / 10. + 5.; - // if (event.getEventNumber() < 1) - // System.out.println("barrel pdx = "+pdx+" pdy = "+pdy); - if ((pdx < 0.) || (pdx > 10.) || (pdy < 0.) || (pdy > 10.)) - System.out.println("out of bound point position!!!!!!!!!!!!!! barrel pdx = " + pdx + " pdy = " + pdy); - - double rnd = Math.random(); // generate total charge here - if (rnd == 1.) - rnd = rnd - 0.0000001; - int i1 = (int) Math.floor(rnd * 1000.); - int i2 = i1 + 1; - - double icharge = 0.; - if (Math.abs(i2 * 1.0 - rnd * 1000.) < 0.00001) { - //icharge = tabCharge[i2]; - icharge = chargeTables.getCharge(i2); - } else { - double k = (rnd * 1000. - i1 * 1.) / (i2 * 1. - rnd * 1000.); - icharge = (chargeTables.getCharge(i1) + chargeTables.getCharge(i2) * k) / (1 + k); - } - - icharge = icharge / 1100. - q0; // done generating total charge - // if (event.getEventNumber() == 1) - // System.out.println("barrel icharge = "+icharge); - - pdx = pdx / 10.; // change pdx into cm unit - pdy = pdy / 10.; // change pdx into cm unit - // if (event.getEventNumber() < 1) - // System.out.println("barrel pdx = "+pdx+" pdy = "+pdy); - - double[] charge = new double[25]; // start charge distribution calulation - int ix1 = (int) Math.floor(pdx / 0.5 * 100); - int iy1 = (int) Math.floor(pdy / 0.5 * 100); - if (ix1 == 200) - ix1 = 199; - if (iy1 == 200) - iy1 = 199; - int ix2 = ix1 + 1; - int iy2 = iy1 + 1; - - double mid1 = 0.; - double mid2 = 0.; - for (int k = 0; k < 25; k++) { - if (Math.abs(ix2 * 1.0 - pdx * 200.) < 0.00001) { - mid1 = chargeTables.getPadCharge(ix2, iy1, k); - mid2 = chargeTables.getPadCharge(ix2, iy2, k); - } else { - double k1 = (pdx * 200. - ix1 * 1.) / (ix2 * 1. - pdx * 200.); - mid1 = (chargeTables.getPadCharge(ix1, iy1, k) + chargeTables.getPadCharge(ix2, iy1, k) * k1) / (1 + k1); - mid2 = (chargeTables.getPadCharge(ix1, iy2, k) + chargeTables.getPadCharge(ix2, iy2, k) * k1) / (1 + k1); - } - - if (Math.abs(iy2 * 1.0 - pdy * 200.) < 0.00001) { - charge[k] = icharge * mid2; - } else { - double k2 = (pdy * 200. - iy1 * 1.) / (iy2 * 1. - pdy * 200.); - charge[k] = icharge * (mid1 + mid2 * k2) / (1 + k2); - } - - // if (event.getEventNumber() < 1) - // System.out.println("barrel charge["+k+"] = "+charge[k]); - - if ((Math.abs(ordered_ids[k]) > 0) && (charge[k] > 0)) // pad exist on the - // neighber or hit pad - // list, and there's - // charge to deposit - { - if (br_rpchits.size() == 0) // first hit, just add - { - RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); - br_rpchits.add(tmphit); - } else // need to see if the pad is already there - { - int hfound = 0; - int j = 0; - while ((hfound == 0) && (j < br_rpchits.size())) { - if (br_rpchits.get(j).getCid() == ordered_ids[k]) { - hfound = 1; - br_rpchits.get(j).addCharge(charge[k]); - } - - j++; - } - - if (hfound == 0) { - RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); - br_rpchits.add(tmphit); - } - } - } - // if (event.getEventNumber() == 1) - // System.out.println("ordered_ids[k] = "+ordered_ids[k]); - // if (event.getEventNumber() == 1) - // System.out.println("barrel cells with charge = "+br_rpchits.size()); - - } // done with charge distribution calulation - - } - - //if (event.getEventNumber() < 10) { - // if ((hits.size() > 0)&&(hits_ec.size() > 0)&&(count_b < 10)) { - // System.out.println(" barrel: cell pos = " + cellPos); // Print cell position. - // System.out.println(" barrel: absolute step pos = " + stepVec); // Print - // absolute step position. - // System.out.println(" barrel: relative step pos = " + VecOp.sub(stepVec, - // cellPos)); // Print relative step position. - - if (fillHistograms) { - c2_xz_cell[event.getEventNumber()].fill(cellPos.z(), cellPos.x()); - c2_xz_step[count_b].fill(stepVec.z(), stepVec.x()); - - c2_yz_cell[event.getEventNumber()].fill(cellPos.z(), cellPos.y()); - c2_yz_step[count_b].fill(stepVec.z(), stepVec.y()); - //} - - c1_cell_x.fill(cellPos.x()); - c1_cell_y.fill(cellPos.y()); - c1_cell_z.fill(cellPos.z()); - c1_step_x.fill(VecOp.sub(stepVec, cellPos).x()); - c1_step_y.fill(VecOp.sub(stepVec, cellPos).y()); - c1_step_z.fill(VecOp.sub(stepVec, cellPos).z()); - } - } - } - } - - // System.out.println("before d cut, barrel points = "+count_points); - // System.out.println("after d cut, barrel points = "+used_points.size()); - used_points.clear(); - count_points = 0; - // System.out.println("before dealing with endcap hits, before thr cut, barrel cells with charge = "+br_rpchits.size()); - - decoder = event.getMetaData(endcapHits).getIDDecoder(); - - count = 0; - for (SimCalorimeterHit hit : endcapHits) { - decoder.setID(hit.getCellID()); - Hep3Vector cellPos = decoder.getPositionVector(); - int ncontribs = hit.getMCParticleCount(); - - long[] nb_ids = decoder.getNeighbourIDs(0, 2, 2); - if (fillHistograms) - c1_nb_ct_ec.fill(nb_ids.length); - - for (int i = 0; i < 25; i++) { - ordered_ids[i] = 0; // make sure to have to way to tell if a neighboring cell is not - // available (id = 0) - } - - // these are two independent, orthogonal vector bases in the readout pads plane - Hep3Vector vx = new BasicHep3Vector(); - Hep3Vector vy = new BasicHep3Vector(); - - // here are the two compoments of the local point position within a cell, according to the - // directions defined by the above bases - double pdx = 0.; - double pdy = 0.; - - long[] nb_final_ids = new long[24]; - int cnt = 0; - if (nb_ids.length > 0) // for all neighbors, need to remove cells that are not on the same - // detector sector (and thus, not on the same RPC) - { - ArrayList<Hep3Vector> bases = new ArrayList<Hep3Vector>(); - - cnt = 0; - for (int i = 0; i < nb_ids.length; i++) { - decoder.setID(nb_ids[i]); - Hep3Vector cellPos_i = decoder.getPositionVector(); - Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); - double dist = Math.sqrt((cellPos.x() - cellPos_i.x()) * (cellPos.x() - cellPos_i.x()) + (cellPos.y() - cellPos_i.y()) * (cellPos.y() - cellPos_i.y()) + (cellPos.z() - cellPos_i.z()) * (cellPos.z() - cellPos_i.z())); - - if ((Math.abs(dist - 10.) < 0.0001) || (Math.abs(dist - 14.142136) < 0.0001) || (Math.abs(dist - 20.) < 0.0001) || (Math.abs(dist - 22.36068) < 0.0001) || (Math.abs(dist - 28.284271) < 0.0001)) { - nb_final_ids[cnt] = nb_ids[i]; - cnt++; - } else { - // System.out.println(" dist = "+dist+" cnt = "+cnt); - } - - if (Math.abs(dist - 10.) < 0.0001) - bases.add(vi0); - - } - - // System.out.println(" bases.size() = "+bases.size()); - if (bases.size() < 3) - System.out.println(" problem finding enough vector bases, bases.size() = " + bases.size()); - vx = bases.get(0); - vy = bases.get(1); - if ((Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) - 10.) < 0.0001) - vy = bases.get(2); - /* - * for (int i = 0; i < bases.size(); i++) { - * System.out.println(" vect "+i+" = "+bases.get(i)); } System.out.println(" vx = "+vx); - * System.out.println(" vy = "+vy); - */ - if (Math.sqrt(Math.abs((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))) > 0.0001) - System.out.println(" vx*vy (endcap) = " + ((vx.x()) * (vy.x()) + (vx.y()) * (vy.y()) + (vx.z()) * (vy.z()))); - - for (int i = 0; i < cnt; i++) { - decoder.setID(nb_final_ids[i]); - Hep3Vector cellPos_i = decoder.getPositionVector(); - Hep3Vector vi0 = VecOp.sub(cellPos_i, cellPos); - - double dx = ((vx.x()) * (vi0.x()) + (vx.y()) * (vi0.y()) + (vx.z()) * (vi0.z())) / 10.; - double dy = ((vy.x()) * (vi0.x()) + (vy.y()) * (vi0.y()) + (vy.z()) * (vi0.z())) / 10.; - int idx = (int) Math.floor(dx / 10. + 0.00001); - int idy = (int) Math.floor(dy / 10. + 0.00001); - int ind = idx + 2 + (idy + 2) * 5; - - // if (event.getEventNumber() < 50) - // System.out.println("endcap i = "+i+" dx = "+dx+" dy = "+dy+" idx = "+idx+" idy = "+idy+" ind = "+ind); - - ordered_ids[ind] = nb_final_ids[i]; - } - ordered_ids[12] = hit.getCellID(); // beside neighboring cells, also need to set id for - // the hit cell itself. - - } - if (fillHistograms) - c1_nb_ct_ec_final.fill(cnt); - - count++; - // if (count < 10000) { - if (count > -1) { - // System.out.println("endcap hit number: "+count); - - for (int i = 0; i < ncontribs; i++) { - float stepPos[] = hit.getStepPosition(i); - Hep3Vector stepVec = new BasicHep3Vector(stepPos[0], stepPos[1], stepPos[2]); - - int ifound = 0; - count_points++; - if (used_points.size() > 0) { - for (Hep3Vector ivec : used_points) { - // if - // (Math.sqrt((ivec.x()-stepVec.x())*(ivec.x()-stepVec.x())+(ivec.y()-stepVec.y())*(ivec.y()-stepVec.y())+(ivec.z()-stepVec.z())*(ivec.z()-stepVec.z())) - // < dcut) ifound = 1; - Hep3Vector tmpv = VecOp.sub(ivec, stepVec); - if (Math.sqrt(VecOp.dot(tmpv, tmpv)) < dcut) - ifound = 1; - } - - if (ifound == 0) - used_points.add(stepVec); - } else { - used_points.add(stepVec); - } - - if (ifound == 0) // only deal with points that survived dcut - { - Hep3Vector pdr = VecOp.sub(stepVec, cellPos); - pdx = VecOp.dot(pdr, vx) / 10. + 5.; - pdy = VecOp.dot(pdr, vy) / 10. + 5.; - // if (event.getEventNumber() < 3) - // System.out.println("endcap pdx = "+pdx+" pdy = "+pdy); - if ((pdx < 0.) || (pdx > 10.) || (pdy < 0.) || (pdy > 10.)) - System.out.println("out of bound point position!!!!!!!!!!!!!! endcap pdx = " + pdx + " pdy = " + pdy); - - double rnd = Math.random(); // generate total charge here - if (rnd == 1.) - rnd = rnd - 0.0000001; - int i1 = (int) Math.floor(rnd * 1000.); - int i2 = i1 + 1; - - double icharge = 0.; - if (Math.abs(i2 * 1.0 - rnd * 1000.) < 0.00001) { - icharge = chargeTables.getCharge(i2); - } else { - double k = (rnd * 1000. - i1 * 1.) / (i2 * 1. - rnd * 1000.); - icharge = (chargeTables.getCharge(i1) + chargeTables.getCharge(i2) * k) / (1 + k); - } - - icharge = icharge / 1100. - q0; // done generating total charge - // if (event.getEventNumber() < 3) - // System.out.println("endcap icharge = "+icharge); - - pdx = pdx / 10.; // change pdx into cm unit - pdy = pdy / 10.; // change pdx into cm unit - // if (event.getEventNumber() < 3) - // System.out.println("endcap pdx = "+pdx+" pdy = "+pdy); - - double[] charge = new double[25]; // start charge distribution calulation - int ix1 = (int) Math.floor(pdx / 0.5 * 100); - int iy1 = (int) Math.floor(pdy / 0.5 * 100); - if (ix1 == 200) - ix1 = 199; - if (iy1 == 200) - iy1 = 199; - int ix2 = ix1 + 1; - int iy2 = iy1 + 1; - - double mid1 = 0.; - double mid2 = 0.; - for (int k = 0; k < 25; k++) { - if (Math.abs(ix2 * 1.0 - pdx * 200.) < 0.00001) { - mid1 = chargeTables.getPadCharge(ix2, iy1, k); - mid2 = chargeTables.getPadCharge(ix2, iy2, k); - } else { - double k1 = (pdx * 200. - ix1 * 1.) / (ix2 * 1. - pdx * 200.); - mid1 = (chargeTables.getPadCharge(ix1, iy1, k) + chargeTables.getPadCharge(ix2, iy1, k) * k1) / (1 + k1); - mid2 = (chargeTables.getPadCharge(ix1, iy2, k) + chargeTables.getPadCharge(ix2, iy2, k) * k1) / (1 + k1); - } - - if (Math.abs(iy2 * 1.0 - pdy * 200.) < 0.00001) { - charge[k] = icharge * mid2; - } else { - double k2 = (pdy * 200. - iy1 * 1.) / (iy2 * 1. - pdy * 200.); - charge[k] = icharge * (mid1 + mid2 * k2) / (1 + k2); - } - - // if (event.getEventNumber() < 3) - // System.out.println("endcap charge["+k+"] = "+charge[k]); - - if ((Math.abs(ordered_ids[k]) > 0) && (charge[k] > 0)) // pad exist on the - // neighber or hit pad - // list, and there's - // charge to deposit - { - if (ec_rpchits.size() == 0) // first hit, just add - { - RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); - ec_rpchits.add(tmphit); - } else // need to see if the pad is already there - { - int hfound = 0; - int j = 0; - while ((hfound == 0) && (j < ec_rpchits.size())) { - if (ec_rpchits.get(j).getCid() == ordered_ids[k]) { - hfound = 1; - ec_rpchits.get(j).addCharge(charge[k]); - } - - j++; - } - - if (hfound == 0) { - RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]); - ec_rpchits.add(tmphit); - } - } - } - - } // done with charge distribution calulation - - } - - //if (event.getEventNumber() < 10) { - // if ((hits.size() > 0)&&(hits_ec.size() > 0)&&(count_b < 10)) { - // System.out.println(" endcap: cell pos = " + cellPos); // Print cell position. - // System.out.println(" endcap: absolute step pos = " + stepVec); // Print - // absolute step position. - // System.out.println(" endcap: relative step pos = " + VecOp.sub(stepVec, - // cellPos)); // Print relative step position. - - if (fillHistograms ) { - c2_xz_cell[event.getEventNumber()].fill(cellPos.z(), cellPos.x()); - c2_xz_step[count_b].fill(stepVec.z(), stepVec.x()); - - c2_yz_cell[event.getEventNumber()].fill(cellPos.z(), cellPos.y()); - c2_yz_step[count_b].fill(stepVec.z(), stepVec.y()); - //} - c1_cell_x.fill(cellPos.x()); - c1_cell_y.fill(cellPos.y()); - c1_cell_z.fill(cellPos.z()); - c1_step_x.fill(VecOp.sub(stepVec, cellPos).x()); - c1_step_y.fill(VecOp.sub(stepVec, cellPos).y()); - c1_step_z.fill(VecOp.sub(stepVec, cellPos).z()); - } - - } - } - } - - // System.out.println("before d cut, endcap points = "+count_points); - // System.out.println("after d cut, endcap points = "+used_points.size()); - used_points.clear(); - count_points = 0; - - if ((barrelHits.size() > 0) && (endcapHits.size() > 0)) - count_b++; - - System.out.println("before thr cut, barrel cells with charge = " + br_rpchits.size()); - System.out.println("before thr cut, endcap cells with charge = " + ec_rpchits.size()); - - ArrayList<RPCHit> br_final = new ArrayList<RPCHit>(); // barrel hits after rpc simulation - ArrayList<RPCHit> ec_final = new ArrayList<RPCHit>(); // endcap hits after rpc simulation - - decoder = event.getMetaData(barrelHits).getIDDecoder(); - for (int i = 0; i < br_rpchits.size(); i++) { - if (br_rpchits.get(i).getCharge() > thr) { - br_final.add(br_rpchits.get(i)); - - //if (event.getEventNumber() < 10) { - if (fillHistograms) { - decoder.setID(br_rpchits.get(i).getCid()); - Hep3Vector cPos = decoder.getPositionVector(); - c2_xz_step[event.getEventNumber()].fill(cPos.z(), cPos.x()); - c2_yz_step[event.getEventNumber()].fill(cPos.z(), cPos.y()); - } - //} - } - } - - decoder = event.getMetaData(endcapHits).getIDDecoder(); - for (int i = 0; i < ec_rpchits.size(); i++) { - if (ec_rpchits.get(i).getCharge() > thr) { - ec_final.add(ec_rpchits.get(i)); - if (fillHistograms) { - //if (event.getEventNumber() < 10) { - decoder.setID(ec_rpchits.get(i).getCid()); - Hep3Vector cPos = decoder.getPositionVector(); - c2_xz_step[event.getEventNumber()].fill(cPos.z(), cPos.x()); - c2_yz_step[event.getEventNumber()].fill(cPos.z(), cPos.y()); - //} - } - } - } - - System.out.println("final barrel hits = " + br_final.size()); - System.out.println("final endcap hits = " + ec_final.size()); - - //br_rpchits.clear(); - //br_final.clear(); - //ec_rpchits.clear(); - //ec_final.clear(); - - int flag = 0; - flag = flag | (1 << LCIOConstants.CHBIT_LONG); - flag = flag | (1 << LCIOConstants.RCHBIT_ID1); - - // Store barrel hits in event. - LCMetaData barrelMetaData = event.getMetaData(barrelHits); - event.put( - barrelMetaData.getName() + collectionOutputTag, - br_final, - CalorimeterHit.class, - flag, - barrelMetaData.getName()); - - // Store endcap hits in event. - LCMetaData endcapMetaData = event.getMetaData(endcapHits); - event.put(endcapMetaData.getName() + collectionOutputTag, - ec_final, - CalorimeterHit.class, - 0, - endcapMetaData.getName()); - } -}
\ No newline at end of file
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