lcsim/sandbox/Partridge/JetOccupancy
diff -N JetOccupancy.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ JetOccupancy.java 2 Nov 2007 17:49:04 -0000 1.1
@@ -0,0 +1,161 @@
+/*
+ * JetOccupancy.java
+ *
+ * Created on October 21, 2007, 10:2 PM
+ *
+ */
+
+import java.util.ArrayList;
+import org.lcsim.contrib.tracking.TrackerHitCheater;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.base.BaseTrackerHitMC;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
+import java.util.List;
+
+import hep.aida.ICloud1D;
+import hep.physics.particle.properties.UnknownParticleIDException;
+import hep.physics.jet.EventShape;
+import hep.physics.particle.Particle;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+/**
+ * Find the occupancy in the core of a jet
+ *
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class JetOccupancy extends Driver {
+ private TrackerHitCheater cheater = new TrackerHitCheater();
+ private EventShape es = new EventShape();
+ private AIDA aida = AIDA.defaultInstance();
+ private double minEventEnergy = 0.;
+
+ /** Creates a new instance of SeedTracker */
+ public JetOccupancy() {
+ }
+
+ protected void process(EventHeader event) {
+
+ List<MCParticle> fsList = new ArrayList<MCParticle>();
+ double eventEnergy = 0.;
+
+ for(MCParticle p : event.getMCParticles()){
+ if (p.getGeneratorStatus() == Particle.FINAL_STATE) {
+ fsList.add(p);
+ eventEnergy += p.getEnergy();
+ }
+ }
+ List<Hep3Vector> momentaList = momenta(fsList);
+ es.setEvent(momentaList);
+
+ if (jetAccept(es)) {
+ aida.cloud1D("Energy of passed events").fill(eventEnergy);
+ Hep3Vector thrust = es.thrustAxis();
+ double ththeta = Math.acos(thrust.z() / thrust.magnitude());
+ double thphi = Math.atan2(thrust.y(),thrust.x());
+
+ List<TrackerHit> trkhits = new ArrayList<TrackerHit>();
+ // Get the SimTrackerHit collections
+ List<List<SimTrackerHit>> allgenhits = event.get(SimTrackerHit.class);
+ // Loop over the various hit collections
+ for (List<SimTrackerHit> simhitcol : allgenhits) {
+ List<TrackerHit> trkhitcol = cheater.makeTrackerHits(simhitcol);
+ boolean status = trkhits.addAll(trkhitcol);
+ }
+ for (TrackerHit hit : trkhits) {
+ double[] pos = hit.getPosition();
+ BasicHep3Vector posvec = new BasicHep3Vector(pos);
+ double r = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
+ double calpha = VecOp.dot(thrust,new BasicHep3Vector(pos));
+ calpha /= thrust.magnitude() * posvec.magnitude();
+ double xt = r * Math.cos(thphi);
+ double yt = r * Math.sin(thphi);
+ double zt = r * Math.cos(ththeta) / Math.sin(thphi);
+ if (calpha<0) {
+ xt = -xt;
+ yt = -yt;
+ zt = -zt;
+ }
+
+ double dx = pos[0] - xt;
+ double dy = pos[1] - yt;
+ double dz = pos[2] - zt;
+ double ds = Math.sqrt(dx*dx+dy*dy);
+ double dist = Math.sqrt(ds*ds+dz*dz);
+// System.out.println(" ds "+ds+" dz "+dz+"calpha"+calpha);
+ BaseTrackerHitMC hitMC = (BaseTrackerHitMC) hit;
+ List<SimTrackerHit> simlist = hitMC.getSimHits();
+ if (simlist.size()==0) System.out.println(" No SimTrackerHits ");
+ if (simlist.size()>0) {
+ SimTrackerHit simhit = (SimTrackerHit) simlist.get(0);
+ int layer = simhit.getLayer();
+ String name = simhit.getSubdetector().getName();
+ BarrelEndcapFlag flag = simhit.getIDDecoder().getBarrelEndcapFlag();
+ if (flag == BarrelEndcapFlag.BARREL && simhit.getPathLength()>10 {
+ System.out.println(" Long path length: "+simhit.getPathLength());
+ }
+ if (flag == BarrelEndcapFlag.BARREL && simhit.getPathLength()<10) {
+// aida.cloud2D("xt vs yt"+name).fill(xt,yt);
+// aida.cloud2D("xm vs ym"+name).fill(pos[0],pos[1]);
+// aida.cloud2D("dist vs cos(alpha)"+name+layer).fill(dist,calpha);
+ double rt = Math.sqrt(xt*xt+yt*yt);
+ double rm = Math.sqrt(Math.pow(pos[0],2)+Math.pow(pos[1],2));
+ double rs = Math.sqrt(Math.pow(simhit.getPoint()[0],2)+Math.pow(simhit.getPoint()[1],2));
+// aida.cloud2D("rm vs zm"+name).fill(rm,pos[2]);
+// aida.cloud2D("rs vs zs"+name).fill(rs,simhit.getPoint()[2]);
+
+ double dsmax = 0.5;
+ double dzmax = 0.5;
+ if (r>15) {
+ dsmax = 1;
+ dzmax = 1.;
+ }
+ if (r>100) {
+ dsmax = 2;
+ dzmax = 50.;
+ }
+ if (r>300) {
+ dsmax = 5;
+ dzmax = 50.;
+ }
+ if (Math.abs(ds)<dsmax && Math.abs(dz)<dzmax) {
+ aida.cloud1D("ds "+name+" "+layer).fill(ds);
+ aida.cloud1D("dz "+name+" "+layer).fill(dz);
+ aida.cloud1D("r "+name+" "+layer).fill(r);
+ }
+ }
+ }
+ }
+ }
+ }
+
+//=========================== Private Functions ==============================//
+
+ private boolean jetAccept(EventShape es){
+ boolean jA = false;
+ double tat = Math.sqrt(es.thrustAxis().x()*es.thrustAxis().x() + es.thrustAxis().y()*es.thrustAxis().y());
+ double cosTAtheta = Math.abs(Math.cos(Math.atan2(tat, es.thrustAxis().z())));
+ if(cosTAtheta < 0.5 && es.thrust().x() >= 0.94) {
+ jA = true;}
+ aida.cloud1D("Cos(theta) of thrust axis").fill(cosTAtheta);
+ aida.cloud1D("Thrust").fill(es.thrust().x());
+ return jA;
+ }
+
+ private List<Hep3Vector> momenta(List<MCParticle> it){
+ List<Hep3Vector> l = new ArrayList<Hep3Vector>();
+ for(int i=0; i<it.size(); i++){
+ Hep3Vector b = it.get(i).getMomentum();
+ l.add(i, b);
+ }
+ return l;
+ }
+
+}
\ No newline at end of file
