Commit in lcsim/src/org/lcsim/rpc on MAIN | |||
PadDigiDriver.java | +113 | -342 | 1.2 -> 1.3 |
refactored to use one method for digitization of both barrel and endcap
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)); + //} +} +*/
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