lcsim/src/org/lcsim/rpc
diff -u -r1.2 -r1.3
--- PadDigiDriver.java 5 Feb 2013 00:25:57 -0000 1.2
+++ PadDigiDriver.java 6 Feb 2013 01:42:38 -0000 1.3
@@ -18,9 +18,6 @@
import org.lcsim.util.Driver;
import org.lcsim.util.lcio.LCIOConstants;
-// FIXME: The code for charge deposition on barrel and endcap is essentially copied and pasted.
-// This should instead be one method called with different parameters.
-// TODO: Move histograms to separate Driver.
public class PadDigiDriver extends Driver
{
// Charge distribution data.
@@ -30,8 +27,11 @@
double minimumCharge = 0.400;
// distance cut in mm
- double distanceCut = 0.90;
-
+ double distanceCut = 0.90;
+
+ // parameter for slope
+ double q0 = -0.062;
+
// Collection names with defaults. These can all be overridden via set methods.
private String barrelInputCollection = "HcalBarrelHits";
private String endcapInputCollection = "HcalEndcapHits";
@@ -81,17 +81,12 @@
{
this.endcapOutputRelationCollection = endcapOutputRelationCollection;
}
-
- public void setup() throws Exception
- {
- // Initialize charge distrubtion tables.
- chargeTables = new ChargeDistributionTables();
- }
-
+
public void startOfData()
{
try {
- setup();
+ // Initialize charge distrubtion tables.
+ chargeTables = new ChargeDistributionTables();
} catch (Exception e) {
throw new RuntimeException(e);
}
@@ -99,51 +94,76 @@
public void process(EventHeader event)
{
- // parameter for slope = 0.084
- // FIXME: What are these parameters? Should be settable on Driver?
- double q0 = -0.062; // ?????
-
+ // Input hits.
List<SimCalorimeterHit> barrelInputHits = event.get(SimCalorimeterHit.class, barrelInputCollection);
- //System.out.println("number of hits in barrel = " + barrelInputHits.size());
-
List<SimCalorimeterHit> endcapInputHits = event.get(SimCalorimeterHit.class, endcapInputCollection);
- //System.out.println("number of hits in endcap = " + endcapInputHits.size());
-
- IDDecoder decoder = event.getMetaData(barrelInputHits).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<BaseCalorimeterHit> br_rpchits = new ArrayList<BaseCalorimeterHit>();
- //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
- ArrayList<BaseCalorimeterHit> ec_rpchits = new ArrayList<BaseCalorimeterHit>();
+ // Output digi hits.
+ ArrayList<BaseCalorimeterHit> barrelRpcHits = new ArrayList<BaseCalorimeterHit>();
+ ArrayList<BaseCalorimeterHit> endcapRpcHits = new ArrayList<BaseCalorimeterHit>();
+ // Relations between digi and sim hits.
ArrayList<LCRelation> barrelRelations = new ArrayList<LCRelation>();
-
ArrayList<LCRelation> endcapRelations = new ArrayList<LCRelation>();
+ // Run barrel digi.
+ IDDecoder decoder = event.getMetaData(barrelInputHits).getIDDecoder();
+ runPadDigi(decoder, barrelInputHits, barrelRpcHits, barrelRelations);
+
+ // Run endcap digi.
+ decoder = event.getMetaData(endcapInputHits).getIDDecoder();
+ runPadDigi(decoder, endcapInputHits, endcapRpcHits, endcapRelations);
+
+ // Collection flags.
+ 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(barrelInputHits);
+ event.put(
+ barrelOutputCollection,
+ barrelRpcHits,
+ CalorimeterHit.class,
+ flag,
+ barrelMetaData.getName());
+
+ // Store final collection of endcap hits in event.
+ LCMetaData endcapMetaData = event.getMetaData(endcapInputHits);
+ event.put(
+ endcapOutputCollection,
+ endcapRpcHits,
+ CalorimeterHit.class,
+ flag,
+ endcapMetaData.getName());
+
+ // Barrel relations.
+ event.put(barrelOutputRelationCollection, barrelRelations, LCRelation.class, 0);
+
+ // Endcap relations.
+ event.put(endcapOutputRelationCollection, endcapRelations, LCRelation.class, 0);
+ }
+
+ private void runPadDigi(
+ IDDecoder decoder,
+ List<SimCalorimeterHit> inputHits,
+ ArrayList<BaseCalorimeterHit> outputHits,
+ ArrayList<LCRelation> barrelRelations) {
+
+ double startTime = System.currentTimeMillis();
+
int count = 0;
- for (SimCalorimeterHit hit : barrelInputHits) {
+
+ ArrayList<Hep3Vector> usedPoints = new ArrayList<Hep3Vector>();
+
+ for (SimCalorimeterHit hit : inputHits) {
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)
- }
+
+ long[] orderedIds = new long[25];
// these are two independent, orthogonal vector bases in the readout pads plane
Hep3Vector vx = new BasicHep3Vector();
@@ -156,9 +176,9 @@
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)
- {
+ // for all neighbors, need to remove cells that are not on the same
+ // detector sector (and thus, not on the same RPC)
+ if (nb_ids.length > 0) {
ArrayList<Hep3Vector> bases = new ArrayList<Hep3Vector>();
cnt = 0;
@@ -202,11 +222,10 @@
int idy = (int) Math.floor(dy / 10. + 0.00001);
int ind = idx + 2 + (idy + 2) * 5;
- ordered_ids[ind] = nb_final_ids[i];
+ orderedIds[ind] = nb_final_ids[i];
}
- ordered_ids[12] = hit.getCellID(); // beside neighboring cells, also need to set id for
+ orderedIds[12] = hit.getCellID(); // beside neighboring cells, also need to set id for
// the hit cell itself.
-
}
count++;
@@ -220,21 +239,21 @@
int ifound = 0;
//count_points++;
- if (used_points.size() > 0) {
- for (Hep3Vector ivec : used_points) {
+ if (usedPoints.size() > 0) {
+ for (Hep3Vector ivec : usedPoints) {
Hep3Vector tmpv = VecOp.sub(ivec, stepVec);
if (Math.sqrt(VecOp.dot(tmpv, tmpv)) < distanceCut)
ifound = 1;
}
if (ifound == 0)
- used_points.add(stepVec);
+ usedPoints.add(stepVec);
} else {
- used_points.add(stepVec);
+ usedPoints.add(stepVec);
}
- if (ifound == 0) // only deal with points that survived dcut
- {
+ // only deal with points that survived dcut
+ if (ifound == 0) {
Hep3Vector pdr = VecOp.sub(stepVec, cellPos);
pdx = VecOp.dot(pdr, vx) / 10. + 5.;
pdy = VecOp.dot(pdr, vy) / 10. + 5.;
@@ -291,16 +310,15 @@
charge[k] = icharge * (mid1 + mid2 * k2) / (1 + k2);
}
- 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
+ // Check if pad exists on the neighber or hit pad list, and there's
+ // charge to deposit.
+ if ((Math.abs(orderedIds[k]) > 0) && (charge[k] > 0))
{
// first hit, just add
- if (br_rpchits.size() == 0) {
+ if (outputHits.size() == 0) {
// Make a new hit.
- BaseCalorimeterHit tmphit = new BaseCalorimeterHit(ordered_ids[k], charge[k], 0.);
- br_rpchits.add(tmphit);
+ BaseCalorimeterHit tmphit = new BaseCalorimeterHit(orderedIds[k], charge[k], 0.);
+ outputHits.add(tmphit);
// Connect digi hit with original hit.
barrelRelations.add(new BaseLCRelation(tmphit, hit));
@@ -310,13 +328,13 @@
int hfound = 0;
int j = 0;
// FIXME: Should use a map here instead. This loop is probably slow.
- while ((hfound == 0) && (j < br_rpchits.size())) {
- if (br_rpchits.get(j).getCellID() == ordered_ids[k]) {
+ while ((hfound == 0) && (j < outputHits.size())) {
+ if (outputHits.get(j).getCellID() == orderedIds[k]) {
hfound = 1;
- br_rpchits.get(j).addCorrectedEnergy(charge[k]);
+ outputHits.get(j).addCorrectedEnergy(charge[k]);
// Connect digi hit with original hit.
- barrelRelations.add(new BaseLCRelation(br_rpchits.get(j), hit));
+ barrelRelations.add(new BaseLCRelation(outputHits.get(j), hit));
}
j++;
@@ -324,8 +342,8 @@
if (hfound == 0) {
// Make a new hit.
- BaseCalorimeterHit tmphit = new BaseCalorimeterHit(ordered_ids[k], charge[k], 0.);
- br_rpchits.add(tmphit);
+ BaseCalorimeterHit tmphit = new BaseCalorimeterHit(orderedIds[k], charge[k], 0.);
+ outputHits.add(tmphit);
// Connect digi hit with original hit.
barrelRelations.add(new BaseLCRelation(tmphit, hit));
@@ -338,279 +356,32 @@
}
}
}
+ double elapsedTime = System.currentTimeMillis() - startTime;
+ System.out.println("digitized " + inputHits.size() + " hits in " + elapsedTime + " ms = " + (elapsedTime / ((double)inputHits.size())) + " ms per hit");
+ }
+}
- used_points.clear();
- //count_points = 0;
-
- decoder = event.getMetaData(endcapInputHits).getIDDecoder();
-
- count = 0;
- for (SimCalorimeterHit hit : endcapInputHits) {
- 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]);
- BaseCalorimeterHit tmphit = new BaseCalorimeterHit(ordered_ids[k], charge[k], 0);
- ec_rpchits.add(tmphit);
-
- // Connect original hit with its digi hit.
- endcapRelations.add(new BaseLCRelation(tmphit, hit));
-
- } 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).getCellID() == ordered_ids[k]) {
- hfound = 1;
- ec_rpchits.get(j).addCorrectedEnergy(charge[k]);
-
- // Connect digi hit with original hit.
- barrelRelations.add(new BaseLCRelation(ec_rpchits.get(j), hit));
- }
-
- j++;
- }
-
- if (hfound == 0) {
- //RPCHit tmphit = new RPCHit(ordered_ids[k], charge[k]);
- BaseCalorimeterHit tmphit = new BaseCalorimeterHit(ordered_ids[k], charge[k], 0);
- ec_rpchits.add(tmphit);
-
- // Connect original hit with its digi hit.
- endcapRelations.add(new BaseLCRelation(tmphit, hit));
- }
- }
- }
-
- } // 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<BaseCalorimeterHit> barrelHitsFinal = new ArrayList<BaseCalorimeterHit>(); // barrel hits after rpc simulation
- ArrayList<BaseCalorimeterHit> endcapHitsFinal = new ArrayList<BaseCalorimeterHit>(); // endcap hits after rpc simulation
- decoder = event.getMetaData(barrelInputHits).getIDDecoder();
- for (int i = 0; i < br_rpchits.size(); i++) {
- // DEBUG: Removed cut for now.
- //if (br_rpchits.get(i).getCorrectedEnergy() > minimumCharge) {
- barrelHitsFinal.add(br_rpchits.get(i));
- //}
- }
+/*
- decoder = event.getMetaData(endcapInputHits).getIDDecoder();
- for (int i = 0; i < ec_rpchits.size(); i++) {
- // DEBUG: Removed cut for now.
- //if (ec_rpchits.get(i).getCorrectedEnergy() > minimumCharge) {
- endcapHitsFinal.add(ec_rpchits.get(i));
- //}
- }
+//// DEBUG: Energy cuts after digi which are removed for now. ////
- //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(barrelInputHits);
- event.put(
- barrelOutputCollection,
- barrelHitsFinal,
- CalorimeterHit.class,
- flag,
- barrelMetaData.getName());
-
- // Store final collection of endcap hits in event.
- LCMetaData endcapMetaData = event.getMetaData(endcapInputHits);
- event.put(
- endcapOutputCollection,
- endcapHitsFinal,
- CalorimeterHit.class,
- flag,
- endcapMetaData.getName());
-
- // Barrel relations.
- event.put(barrelOutputRelationCollection, barrelRelations, LCRelation.class, 0);
-
- // Endcap relations.
- event.put(endcapOutputRelationCollection, endcapRelations, LCRelation.class, 0);
- }
-}
\ No newline at end of file
+ArrayList<BaseCalorimeterHit> barrelHitsFinal = new ArrayList<BaseCalorimeterHit>(); // barrel hits after rpc simulation
+ArrayList<BaseCalorimeterHit> endcapHitsFinal = new ArrayList<BaseCalorimeterHit>(); // endcap hits after rpc simulation
+
+decoder = event.getMetaData(barrelInputHits).getIDDecoder();
+for (int i = 0; i < barrelRpcHits.size(); i++) {
+ // DEBUG: Removed cut for now.
+ //if (br_rpchits.get(i).getCorrectedEnergy() > minimumCharge) {
+ barrelHitsFinal.add(barrelRpcHits.get(i));
+ //}
+}
+
+decoder = event.getMetaData(endcapInputHits).getIDDecoder();
+for (int i = 0; i < endcapRpcHits.size(); i++) {
+ // DEBUG: Removed cut for now.
+ //if (ec_rpchits.get(i).getCorrectedEnergy() > minimumCharge) {
+ endcapHitsFinal.add(endcapRpcHits.get(i));
+ //}
+}
+*/