lcsim/src/org/lcsim/rpc
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
lcsim/src/org/lcsim/rpc
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