lcsim/src/org/lcsim/contrib/compile/SteveMagill
diff -N BTrMipClus.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ BTrMipClus.java 23 Jul 2005 01:49:56 -0000 1.1
@@ -0,0 +1,837 @@
+// import java.util.List;
+import java.util.*;
+import hep.physics.jet.JetFinder;
+import hep.physics.jet.DurhamJetFinder;
+import hep.physics.vec.HepLorentzVector;
+import hep.physics.vec.Hep3Vector;
+// import org.lcsim.event.Cluster;
+// import org.lcsim.event.EventHeader;
+import org.lcsim.event.*;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.fourvec.Momentum4Vector;
+import org.lcsim.geometry.*;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.nn.NearestNeighborCluster;
+// import org.lcsim.recon.ztracking.FoundTrack;
+// import org.lcsim.geometry.Detector;
+
+public class BTrMipClus extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ private boolean _initialized = false;
+// The number of layers in the EM Calorimeter Barrel
+ private int _numbemlayers;
+// The number of layers in the Hadronic Calorimeter Barrel
+ private int _numbhadlayers;
+// The EM Barrel sampling fraction
+ private double _sfEMB;
+// The HAD Barrel sampling fraction
+ private double _sfHB;
+// The EM calorimeter hits, by layer
+ private List[] _emBLayerHits;
+// The HAD calorimeter hits, by layer
+ private List[] _hadBLayerHits;
+// The radii of the barrel calorimeter layers
+ private double[] _emBRadii;
+ private double[] _hadBRadii;
+// modified Barrel helix swimmer for org.lcsim test
+ private Helix[] _emBSwimmers; //EM
+ private Helix[] _hadBSwimmers; //HAD
+// Z extent of the central barrel calorimeters.
+ private double _embZ; //EM
+ private double _hadbZ; //HAD
+// The central magnetic field strength
+ private double[] _fieldStrength;
+
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+// System.out.println(events);
+
+// Initialize things here
+ if(!_initialized)
+ {
+ System.out.println("Detector is "+event.getDetectorName());
+ if (event.getDetectorName().equals("sidmay05")) System.out.println(" SiD!! ARE YOU SERIOUS?");
+ Detector det = event.getDetector();
+ double[] zero = {0, 0, 0};
+ _fieldStrength = det.getFieldMap().getField(zero);
+ System.out.println("B Field " +_fieldStrength[2]+ " Tesla");
+
+ CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
+ _embZ = emb.getZMin();
+ System.out.println("EM Barrel z " +_embZ);
+// _numbemlayers = 30;
+ _numbemlayers = emb.getLayering().getLayerCount();
+ System.out.println("EM Barrel Layers " +_numbemlayers);
+ _emBLayerHits = new ArrayList[_numbemlayers];
+ _emBRadii = new double[_numbemlayers];
+ _emBSwimmers = new Helix[_numbemlayers];
+ for (int i=0; i<_numbemlayers; ++i)
+ {
+ _emBLayerHits[i] = new ArrayList();
+ _emBRadii[i]=emb.getLayering().getDistanceToLayerSensorMid(i);
+ System.out.println("EM Barrel Layer Number " +i+ " EM Barrel Radius " +_emBRadii[i]);
+ _emBSwimmers[i] = new Helix(_fieldStrength[2], _emBRadii[i], Math.abs(_embZ));
+ _emBSwimmers[i].setCurvatureFlip(true);
+ } // end loop over em layers
+
+ CylindricalCalorimeter hadb = ((CylindricalCalorimeter) det.getSubdetectors().get("HADBarrel"));
+ _hadbZ = hadb.getZMin();
+ System.out.println("HAD Barrel z " +_hadbZ);
+ _numbhadlayers = hadb.getLayering().getLayerCount();
+ System.out.println("HAD Barrel Layers " +_numbhadlayers);
+ _hadBLayerHits = new ArrayList[_numbhadlayers];
+ _hadBRadii = new double[_numbhadlayers];
+ _hadBSwimmers = new Helix[_numbhadlayers];
+ for (int i=0; i<_numbhadlayers; ++i)
+ {
+ _hadBLayerHits[i] = new ArrayList();
+ _hadBRadii[i]=hadb.getLayering().getDistanceToLayerSensorMid(i);
+ System.out.println("HAD Barrel Layer Number " +i+ " HAD Barrel Radius " +_hadBRadii[i]);
+ _hadBSwimmers[i] = new Helix(_fieldStrength[2], _hadBRadii[i], Math.abs(_hadbZ));
+ _hadBSwimmers[i].setCurvatureFlip(true);
+ } // end loop over had layers
+
+ // Get (set) sampling fractions here
+ _sfEMB = 0.01207; // for now this is tested sampling fraction for W/Si EMCAL in SiD Detector
+// _sfHB = 0.000025; // for now this is tested sampling fraction for SS/RPC HCAL in SiD Detector
+ _sfHB = 0.06;
+
+ _initialized = true;
+ } // end of initialization section
+
+ // get calorimeter hits by layer - need to calculate hit density for cells
+ double EMBESum=0; double HADBESum=0;
+ double ESumBTotCAL=0;
+ int nhitsBTotCAL=0;
+ double[][] celbtheta = new double[5000][_numbemlayers+_numbhadlayers];
+ double[][] celbphi = new double[5000][_numbemlayers+_numbhadlayers];
+ double[][] celbden = new double[5000][_numbemlayers+_numbhadlayers];
+ long[][] celbid = new long[5000][_numbemlayers+_numbhadlayers];
+ int[] nhitbLayer = new int[_numbemlayers+_numbhadlayers];
+ int nhitsEMB=0; int nhitsHADB=0;
+
+ // this uses SimCalorimeterHit, but need to know name of collection for this.
+ List<SimCalorimeterHit> embHits = event.getSimCalorimeterHits("EcalBarrHits");
+ CalorimeterIDDecoder embdecoder = (CalorimeterIDDecoder) event.getMetaData(embHits).getIDDecoder();
+ // create a map of cells keyed on their index
+ Map<Long, CalorimeterHit> embhitmap = new HashMap<Long, CalorimeterHit>();
+
+ for (SimCalorimeterHit embhit : embHits)
+ {
+ long hitID = embhit.getCellID();
+ embdecoder.setID(hitID);
+ int emblayer = embdecoder.getLayer();
+// System.out.println("EM layer num = " + emblayer);
+ celbtheta[nhitbLayer[emblayer]][emblayer] = embdecoder.getTheta();
+ celbphi[nhitbLayer[emblayer]][emblayer] = embdecoder.getPhi();
+ celbid[nhitbLayer[emblayer]][emblayer] = hitID;
+ nhitbLayer[emblayer]++;
+ embhitmap.put(hitID, embhit);
+// if (embhit.getEnergy()<0.00061)
+// {
+// aida.cloud1D("DIAG : EM B Calib hit E").fill(embhit.getEnergy()/_sfEMB);
+// }
+ EMBESum += embhit.getEnergy();
+ nhitsEMB++;
+ }
+ aida.cloud1D("Nhits in original EMBhitmap").fill(embhitmap.size());
+
+ List<SimCalorimeterHit> hadbHits = event.getSimCalorimeterHits("HcalBarrHits");
+ CalorimeterIDDecoder hadbdecoder = (CalorimeterIDDecoder) event.getMetaData(hadbHits).getIDDecoder();
+ // create a map of cells keyed on their index
+ Map<Long, CalorimeterHit> hadbhitmap = new HashMap<Long, CalorimeterHit>();
+
+ for (SimCalorimeterHit hadbhit : hadbHits)
+ {
+ long hitID = hadbhit.getCellID();
+ hadbdecoder.setID(hitID);
+ int hadblayer = hadbdecoder.getLayer();
+ celbtheta[nhitbLayer[_numbemlayers+hadblayer]][_numbemlayers+hadblayer] = hadbdecoder.getTheta();
+ celbphi[nhitbLayer[_numbemlayers+hadblayer]][_numbemlayers+hadblayer] = hadbdecoder.getPhi();
+ celbid[nhitbLayer[_numbemlayers+hadblayer]][_numbemlayers+hadblayer] = hitID;
+ nhitbLayer[_numbemlayers+hadblayer]++;
+ hadbhitmap.put(hitID, hadbhit);
+// if (hadbhit.getEnergy()<0.000004)
+// {
+// aida.cloud1D("DIAG : HAD B Vis hit E").fill(hadbhit.getEnergy());
+// aida.cloud1D("DIAG : HAD B Calib hit E").fill(hadbhit.getEnergy()/_sfHB);
+// }
+ HADBESum += hadbhit.getEnergy();
+ nhitsHADB++;
+ }
+ aida.cloud1D("Nhits in original HADBhitmap").fill(hadbhitmap.size());
+
+ for (int i=0; i<_numbemlayers+_numbhadlayers; ++i)
+ {
+ aida.cloud2D("NHits vs layer").fill(i,nhitbLayer[i]);
+ // calculate density for each cell hit in layer - use 5X5 window?
+ for (int j=0; j<nhitbLayer[i]; ++j)
+ {
+ celbden[j][i] = 1;
+ for (int k=0; k<nhitbLayer[i]; ++k)
+ {
+ if (k == j) continue;
+ // in ECAL, for 5 mm X 5 mm cells, each cell is ~ .00315 X .00315 rad
+ // .0075 corresponds to 5X5 without outermost 4 corners, 0.005 is 3X3
+ // in HCAL, for 1 X 1 cm cells, each cell is ~ .00525 X .00525 rad
+ // .018 is 7X7 without corners (3 cells missing at each corner)
+ double dmine = 0.005;
+ // double dminh = 0.0125;
+ double dminh = 3*dmine;
+ double dmin = dmine;
+ double deltheta = Math.abs(celbtheta[k][i]-celbtheta[j][i]);
+ aida.cloud1D("Delta Theta Cell").fill(deltheta);
+ double delphi = Math.abs(celbphi[k][i]-celbphi[j][i]);
+ if (delphi > Math.PI) delphi = 2*Math.PI-delphi;
+ aida.cloud1D("Delta Phi Cell").fill(delphi);
+ double distance = Math.sqrt(deltheta*deltheta+delphi*delphi);
+ if (i >= _numbemlayers) dmin = dminh;
+ if (distance < dmin)
+ {
+ if (i < 30) aida.cloud1D("Distance between EM CAL hits").fill(distance);
+ if (i > 29) aida.cloud1D("Distance between HAD CAL hits").fill(distance);
+ if (i < 30) celbden[j][i] += 0.00315/distance;
+ if (i > 29) celbden[j][i] += 0.00525/distance;
+ }
+ }
+ if (nhitbLayer[i] > 0 && i < 30) aida.cloud1D("EM CAL cell density").fill(celbden[j][i]);
+ if (nhitbLayer[i] > 0 && i > 29) aida.cloud1D("HAD CAL cell density").fill(celbden[j][i]);
+ if (nhitbLayer[i] > 0) aida.cloud2D("CAL cell density vs Layer").fill(i,celbden[j][i]);
+ if (i == 0) aida.cloud1D("Cell Density layer 0").fill(celbden[j][i]);
+ }
+ }
+ ESumBTotCAL = EMBESum + HADBESum;
+ nhitsBTotCAL = nhitsEMB + nhitsHADB;
+ aida.cloud2D("NHits EM B vs EM B Calib ESum").fill(EMBESum/_sfEMB,nhitsEMB);
+ aida.cloud2D("NHits HAD B vs HAD B Vis ESum").fill(HADBESum,nhitsHADB);
+ if (EMBESum/_sfEMB < 0.40)
+ {
+ aida.cloud1D("EM B Calib ESum, mips only").fill(EMBESum/_sfEMB);
+ aida.cloud1D("NHits EM B, mips only").fill(nhitsEMB);
+ aida.cloud1D("NHits HAD B, EM B mips only").fill(nhitsHADB);
+ }
+ if (nhitsEMB < 1 && nhitsHADB > 0)
+ {
+ aida.cloud1D("NHits HAD B, no EM B hits").fill(nhitsHADB);
+ aida.cloud1D("HAD B Vis ESum, no EM B Hits").fill(HADBESum);
+ }
+ if (nhitsHADB < 1 && nhitsEMB > 0)
+ {
+ aida.cloud1D("NHits EM B, no HAD B hits").fill(nhitsEMB);
+ }
+ aida.cloud1D("Total B CAL NHits").fill(nhitsBTotCAL);
+ aida.cloud1D("Total B CAL Vis ESum").fill(ESumBTotCAL);
+ aida.cloud1D("Total B CAL EM Analog + HAD Digital ESum").fill(EMBESum/_sfEMB+nhitsHADB*0.115*0.2);
+ aida.cloud1D("Total EM B Vis ESum").fill(EMBESum);
+ aida.cloud1D("Total EM B Calib ESum").fill(EMBESum/_sfEMB);
+ aida.cloud1D("Total HAD B Vis ESum").fill(HADBESum);
+
+ /** This section finds mip clusters in cal attached to a reconstructed (or true) track.
+ * The starting point is the first layer of the ECAL and the endpoint is the point
+ * of first interaction in the CAL. This end point can be used to start Track/CAL
+ * object association later.
+ * 1) get list of reconstructed tracks or true MC charged particles (determined by
+ * user
+ * 2) loop over all tracks, swim to layer, check for mip hit.
+ * 3) continue until either a) no hit in next layer, or b) large energy or density of hits
+ * is found in next layer.
+ * 4) list of hits forms a cluster.
+ **/
+ // Here is the track list - either use cheater or rec tracks
+// List<List<FoundTrack>> trackCollections = event.get(FoundTrack.class);
+// List<Track> chptrks = event.getTracks();
+
+ // First, this is for true tracks - MC charged final states
+ int ntrkswim = 0;
+ double[][] trbtheta = new double[200][_numbemlayers+_numbhadlayers];
+ double[][] trbphi = new double[200][_numbemlayers+_numbhadlayers];
+ double[] trbpE = new double[200];
+ double[] trbpp = new double[200];
+ double[] trbpct = new double[200];
+ double[] trbppt = new double[200];
+ List<MCParticle> mcps = event.getMCParticles();
+ for (MCParticle mcp : mcps)
+ {
+ // must be final state particle
+// if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) continue;
+ if (mcp.getGeneratorStatus() == mcp.DOCUMENTATION || mcp.getGeneratorStatus() == 0) continue;
+// System.out.println("Generator Status " +mcp.getGeneratorStatus());
+ // must be charged particle
+ int iq = (int)mcp.getCharge();
+ if (iq == 0) continue;
+ Hep3Vector origin = mcp.getOrigin();
+ Hep3Vector momentum = mcp.getMomentum();
+ double chpE = mcp.getEnergy();
+ double chpmom = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()+mcp.getPZ()*mcp.getPZ());
+ double chpct = mcp.getPZ()/chpmom;
+ double chppt = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY());
+ aida.cloud1D("All Ch Particle pT").fill(chppt);
+ // cut on barrel for now - later add endcaps
+ if (Math.abs(chpct)<0.80)
+ {
+ aida.cloud1D("Ch Particle pT in Barrel").fill(chppt);
+ // only check tracks that make it into CAL
+ // 0.9525 GeV to get to inner radius of ECAL, 0.975 GeV to get to shower max of ECAL, 1.04 GeV
+ // to get to inner radius of HCAL, 1.75 GeV to get to outer radius of HCAL
+ if (chppt>0.9525)
+ {
+ aida.cloud1D("Ch Particle pT for swimming").fill(chppt);
+ // swim through the EM Barrel calorimeter
+ for(int i=0; i<_numbemlayers; ++i)
+ {
+ _emBSwimmers[i].swim(momentum, origin, iq);
+ double[] intersection = _emBSwimmers[i].getIntersect();
+ double phi = Math.atan2(intersection[1],intersection[0]);
+ if(phi<0) phi+=2*Math.PI;
+ aida.cloud1D("Track Phi").fill(phi);
+ double theta = Math.atan(_emBRadii[i]/intersection[2]);
+ if(theta<0) theta+=Math.PI;
+ aida.cloud1D("Track Theta").fill(theta);
+ trbtheta[ntrkswim][i] = theta;
+ trbphi[ntrkswim][i] = phi;
+ }
+ // now swim through HAD Barrel Calorimeter
+ for(int i=0; i<_numbhadlayers; ++i)
+ {
+ _hadBSwimmers[i].swim(momentum, origin, iq);
+ double[] intersection = _hadBSwimmers[i].getIntersect();
+ double phi = Math.atan2(intersection[1],intersection[0]);
+ if(phi<0) phi+=2*Math.PI;
+ double theta = Math.atan(_hadBRadii[i]/intersection[2]);
+ if(theta<0) theta+=Math.PI;
+ trbtheta[ntrkswim][i+_numbemlayers] = theta;
+ trbphi[ntrkswim][i+_numbemlayers] = phi;
+ }
+ trbpp[ntrkswim] = chpmom;
+ trbpE[ntrkswim] = chpE;
+ trbpct[ntrkswim] = chpct;
+ trbppt[ntrkswim] = chppt;
+ ntrkswim++;
+ }
+ } else if (Math.abs(chpct)>0.8 && Math.abs(chpct)<0.993)
+ {
+ // endcap tracks
+ aida.cloud1D("End Cap Track pT").fill(chppt);
+ } else
+ {
+ // beampipe tracks
+ aida.cloud1D("Beampipe Track pT").fill(chppt);
+ }
+ }
+ aida.cloud1D("Number of swum Tracks").fill(ntrkswim);
+
+ // Now, can find interaction layer by comparing track position with
+ // CAL hits - first loop over all barrel tracks
+ List<BasicCluster> embmipclusters = new ArrayList();
+ List<BasicCluster> hadbmipclusters = new ArrayList();
+ List<BasicCluster> calbmipclusters = new ArrayList();
+ List<Integer> itracks = new ArrayList();
+ Map<Integer, BasicCluster> trmipclusmap = new HashMap<Integer, BasicCluster>();
+ int[] IntLay = new int[200];
+ int[] nhitMipClus = new int[200];
+ long[][] hitidMipClus = new long[200][1000];
+ for (int i=0; i<ntrkswim; ++i)
+ {
+ itracks.add(i);
+ // map integer into cluster to link tracks and mips
+ Map<Long, CalorimeterHit> ebclhitmap = new HashMap<Long, CalorimeterHit>();
+ Map<Long, CalorimeterHit> hbclhitmap = new HashMap<Long, CalorimeterHit>();
+ BasicCluster ebmclus = new BasicCluster();
+ BasicCluster hbmclus = new BasicCluster();
+ BasicCluster bmipclus = new BasicCluster();
+ boolean nointeraction = true;
+ IntLay[i] = 64; // set to maximum for each track to start
+ // now, loop over all Barrel layers
+ for (int j=0; j<_numbemlayers+_numbhadlayers; ++j)
+ {
+ // continue if no interaction yet, if already have interaction, skip
+ if (nointeraction)
+ {
+ double prden = 0;
+ int ntcmtch = 0;
+ long[] layerids = new long[200];
+ // loop over all hits in this layer
+ for (int k=0; k<nhitbLayer[j]; ++k)
+ {
+ // here is some diagnostic stuff - only need for layer 0, some for first 5 layers
+ if (j == 0)
+ {
+ aida.cloud1D("DIAG - CAL cell theta").fill(celbtheta[k][j]);
+ aida.cloud1D("DIAG - CAL cell cos theta").fill(Math.cos(celbtheta[k][j]));
+ aida.cloud1D("DIAG - CAL cell phi").fill(celbphi[k][j]);
+ aida.cloud2D("DIAG - CAL cell Theta vs Phi").fill(celbphi[k][j],celbtheta[k][j]);
+ aida.cloud1D("DIAG - Track theta").fill(trbtheta[i][j]);
+ aida.cloud1D("DIAG - Track cos theta").fill(Math.cos(trbtheta[i][j]));
+ aida.cloud1D("DIAG - Track phi").fill(trbphi[i][j]);
+ aida.cloud2D("DIAG - Track Theta vs Phi").fill(trbphi[i][j],trbtheta[i][j]);
+ }
+ if (j == 1) aida.cloud2D("DIAG - CAL cell Theta vs Phi 1").fill(celbphi[k][j],celbtheta[k][j]);
+ if (j == 2) aida.cloud2D("DIAG - CAL cell Theta vs Phi 2").fill(celbphi[k][j],celbtheta[k][j]);
+ if (j == 3) aida.cloud2D("DIAG - CAL cell Theta vs Phi 3").fill(celbphi[k][j],celbtheta[k][j]);
+ if (j == 4) aida.cloud2D("DIAG - CAL cell Theta vs Phi 4").fill(celbphi[k][j],celbtheta[k][j]);
+
+ // get cell-track distances and test for matches
+ double trcedmin = 0.0075; // eventually make this an argument for the mip cluster code - make min .0075 for ECAL and 2*min for HCAL
+ if (j >= _numbemlayers) trcedmin = 2*trcedmin;
+ double dtheta = Math.abs(celbtheta[k][j] - trbtheta[i][j]);
+ double dphi = Math.abs(celbphi[k][j] - trbphi[i][j]);
+ if (dphi > Math.PI) dphi = 2*Math.PI-dphi;
+ dphi = Math.abs(dphi);
+ aida.cloud1D("DIAG - delta theta track cell").fill(dtheta);
+ aida.cloud1D("DIAG - delta phi track cell").fill(dphi);
+ double trcedist = Math.sqrt(dtheta*dtheta+dphi*dphi);
+ if (j == 0) aida.cloud1D("DIAG - Track-Cell Distance in Layer 0").fill(trcedist);
+ if (j < 30) aida.cloud1D("DIAG - EM Track-Cell Distance").fill(trcedist);
+ if (j > 29) aida.cloud1D("DIAG - HAD Track-Cell Distance").fill(trcedist);
+ if (trcedist < 2*trcedmin) aida.cloud1D("Near Track-Cell Distance").fill(trcedist);
+ if (trcedist < trcedmin)
+ {
+ if (j < _numbemlayers) aida.cloud1D("DIAG - Matched Track-Cell Distance EM").fill(trcedist);
+ if (j >= _numbemlayers) aida.cloud1D("DIAG - Matched Track-Cell Distance HAD").fill(trcedist);
+ prden += celbden[k][j];
+ layerids[ntcmtch] = celbid[k][j];
+ ntcmtch++;
+ }
+ }
+ aida.cloud1D("Density Sum per layer").fill(prden);
+ if (ntcmtch > 0) prden = prden/ntcmtch;
+ // should probably make this an argument - the density sum limit
+ if (prden == 0 || prden > 2.5)
+ {
+ // have interaction - form mip cluster from hits in cluster hitmap
+ nointeraction = false;
+ IntLay[i] = j;
+ aida.cloud1D("Cell Density Sum for interactions").fill(prden);
+// if (j == 0) aida.cloud1D("Cell Density for matches in layer 0").fill(prden);
+// if (j == 0 && prden == 0) aida.cloud1D("Track pT for matches in Layer 0").fill(trbppt[i]);
+// if (j == 0 && prden == 0) aida.cloud1D("Track theta for matches in Layer 0").fill(trbtheta[i][j]);
+ if (j == 0) continue; // if interaction in first layer, just get out
+ trmipclusmap.put(i, bmipclus);
+ if (j < _numbemlayers)
+ {
+ // only have EM cluster
+ aida.cloud1D("DIAG - Nhits in EMClus hitmap for match").fill(ebclhitmap.size());
+ if (ebclhitmap.size() > 0)
+ {
+ for (CalorimeterHit h : ebclhitmap.values())
+ {
+ ebmclus.addHit(h); // add hit to EM B Mip cluster
+ bmipclus.addHit(h); // add hit to CAL B Mip Cluster
+ }
+ }
+ embmipclusters.add(ebmclus);
+ calbmipclusters.add(bmipclus);
+ } else
+ {
+ // can have both EM and HAD clusters
+ if (ebclhitmap.size() > 0)
+ {
+ for (CalorimeterHit eh : ebclhitmap.values())
+ {
+ ebmclus.addHit(eh);
+ bmipclus.addHit(eh);
+ }
+ }
+ embmipclusters.add(ebmclus);
+ if (j != _numbemlayers) // don't make hadmipcluster if interaction in first layer of HCAL
+ {
+ if (hbclhitmap.size() > 0)
+ {
+ for (CalorimeterHit hh : hbclhitmap.values())
+ {
+ hbmclus.addHit(hh);
+ bmipclus.addHit(hh);
+ }
+ }
+ hadbmipclusters.add(hbmclus);
+ }
+ calbmipclusters.add(bmipclus);
+ }
+ } else
+ {
+ // add hits to temporary mip cluster, remove hits from hitmap
+ aida.cloud1D("Cell Density for Noninteractions").fill(prden);
+ for (int l=0; l<ntcmtch; ++l)
+ {
+ hitidMipClus[i][nhitMipClus[i]] = layerids[l];
+ nhitMipClus[i]++;
+ if (j < _numbemlayers)
+ {
+ CalorimeterHit t = embhitmap.get(layerids[l]);
+ ebclhitmap.put(layerids[l],t);
+ embhitmap.remove(layerids[l]);
+ }
+ if (j >= _numbemlayers)
+ {
+ CalorimeterHit t = hadbhitmap.get(layerids[l]);
+ hbclhitmap.put(layerids[l],t);
+ hadbhitmap.remove(layerids[l]);
+ }
+ }
+ }
+ }
+ }
+ }
+ // Name the mip clusters and add to event
+ if (embmipclusters.size() > 0) event.put("EBMipClusters",embmipclusters);
+ if (hadbmipclusters.size() > 0) event.put("HBMipClusters",hadbmipclusters);
+ if (calbmipclusters.size() > 0) event.put("CALBMipClusters",calbmipclusters);
+ // Done with mip cluster finding
+
+ // some histograms
+ aida.cloud1D("Nhits in EMBhitmap after MipClus").fill(embhitmap.size());
+ aida.cloud1D("Nhits in HADBhitmap after MipClus").fill(hadbhitmap.size());
+ aida.cloud1D("Num EM B Mip Clusters").fill(embmipclusters.size());
+ aida.cloud1D("Num HAD B Mip Clusters").fill(hadbmipclusters.size());
+ aida.cloud1D("Num CAL B Mip Clusters").fill(calbmipclusters.size());
+
+ for (BasicCluster embmipcluster : embmipclusters)
+ {
+ aida.cloud1D("EM B Mip Cluster E").fill(embmipcluster.getEnergy()/_sfEMB);
+ aida.cloud1D("EM B Mip Cluster NHits").fill(embmipcluster.getCalorimeterHits().size());
+ }
+ for (BasicCluster hadbmipcluster : hadbmipclusters)
+ {
+ aida.cloud1D("HAD B Mip Cluster NHits").fill(hadbmipcluster.getCalorimeterHits().size());
+ }
+
+ // make a map linking hits to mip clusters
+ Map<CalorimeterHit, BasicCluster> calmiphitclus = new HashMap<CalorimeterHit, BasicCluster>();
+ for (BasicCluster calbmipcluster : calbmipclusters)
+ {
+ aida.cloud1D("CAL B Mip Cluster NHits").fill(calbmipcluster.getCalorimeterHits().size());
+ for (int i=0; i<calbmipcluster.getCalorimeterHits().size(); ++i)
+ {
+ CalorimeterHit cmiph = calbmipcluster.getCalorimeterHits().get(i);
+ calmiphitclus.put(cmiph,calbmipcluster);
+ }
+ }
+
+ for (int i=0; i<ntrkswim; ++i)
+ {
+ aida.cloud1D("Number of mip hits per Track").fill(nhitMipClus[i]);
+ aida.cloud1D("CAL Interaction Layer per Track").fill(IntLay[i]);
+ aida.cloud2D("Number of mip hits vs Int Layer").fill(IntLay[i],nhitMipClus[i]);
+// System.out.println("Track " +i+ " IL " +IntLay[i]+ " Energy " +trbpE[i]+ " Momentum " +trbpp[i]+ " Phi " +trbphi[i][0]);
+ }
+ // End of Track-Mip section - break here eventually
+
+ // Now, go on to Track-Shower matching using endpoint of matched Track-Mips
+ // first, cluster remaining hits with nearest neighbor clusterer
+ // e.g., 1,1,0,1 gives nn clusters in each layer with min hits of 1
+ // first, make a copy of hitmaps to send to NN clusterer (deletes all hits)
+ Map<Long, CalorimeterHit> ebnnhitmap = new HashMap<Long, CalorimeterHit>();
+ ebnnhitmap.putAll(embhitmap);
+ Map<Long, CalorimeterHit> hbnnhitmap = new HashMap<Long, CalorimeterHit>();
+ hbnnhitmap.putAll(hadbhitmap);
+
+ int _dU = 2;
+ int _dV = 2;
+ int _dLayer = 2;
+ int _minNcells = 4;
+ // get EM clusters
+ List<Cluster> ebclusters = new ArrayList();
+ for(;;)
+ {
+ if (ebnnhitmap.isEmpty()) break;
+ Long k = ebnnhitmap.keySet().iterator().next();
+ CalorimeterHit emhit = ebnnhitmap.get(k);
+ // start a cluster, constructor will aggregate remaining hits unto itself
+ NearestNeighborCluster nnclus = new NearestNeighborCluster(embdecoder, ebnnhitmap, emhit, k, _dU, _dV, _dLayer);
+ if(nnclus.getCalorimeterHits().size()>_minNcells)
+ {
+ ebclusters.add(nnclus);
+ }
+ }
+// String name = event.getMetaData(collection).getName();
+ if (ebclusters.size() > 0) event.put("EMBNNClusters",ebclusters);
+ for (Cluster ebcluster : ebclusters)
+ {
+ for (int i=0; i<ebcluster.getCalorimeterHits().size(); ++i)
+ {
+ // remove cluster hits from EM B hitmap
+ CalorimeterHit ebh = ebcluster.getCalorimeterHits().get(i);
+ long ebhid = ebh.getCellID();
+ embhitmap.remove(ebhid);
+ }
+ }
+ aida.cloud1D("Nhits in EMBhitmap after NN clustering").fill(embhitmap.size());
+ //might be able to ID photons with these clusters? - see below
+
+ _dU = 4;
+ _dV = 4;
+ _dLayer = 4;
+ _minNcells = 2;
+ // get had clusters
+ List<Cluster> hbclusters = new ArrayList();
+ for (;;)
+ {
+ if (hbnnhitmap.isEmpty()) break;
+ Long k = hbnnhitmap.keySet().iterator().next();
+ CalorimeterHit hadhit = hbnnhitmap.get(k);
+ NearestNeighborCluster nnclus = new NearestNeighborCluster(hadbdecoder, hbnnhitmap, hadhit, k, _dU, _dV, _dLayer);
+// aida.cloud1D("Num of hits in had cluster").fill(nnclus.getCalorimeterHits().size());
+ if (nnclus.getCalorimeterHits().size()>_minNcells)
+ {
+ hbclusters.add(nnclus);
+// aida.cloud1D("Num of hits in had cluster").fill(nnclus.getCalorimeterHits().size());
+ }
+ }
+ if (hbclusters.size() > 0) event.put("HADBNNClusters",hbclusters);
+ // map hits in cluster to their cluster for combining later
+ Map<CalorimeterHit, Cluster> hbnnhitclus = new HashMap<CalorimeterHit, Cluster>();
+ for (Cluster hbcluster : hbclusters)
+ {
+ for (int i=0; i<hbcluster.getCalorimeterHits().size(); ++i)
+ {
+ // remove cluster hits from HAD B hitmap
+ CalorimeterHit hbh = hbcluster.getCalorimeterHits().get(i);
+ long hbhid = hbh.getCellID();
+ hadbhitmap.remove(hbhid);
+ // fill cluster-hit map
+ hbnnhitclus.put(hbh,hbcluster);
+ }
+ }
+ aida.cloud1D("Nhits in HADBhitmap after NN clustering").fill(hadbhitmap.size());
+ aida.cloud1D("Num HAD B NNClusters").fill(hbclusters.size());
+ for (Cluster hbcluster : hbclusters)
+ {
+ aida.cloud1D("Num of hits in HAD B NNCluster").fill(hbcluster.getCalorimeterHits().size());
+ aida.cloud1D("Vis Energy of HAD B NNCluster").fill(hbcluster.getEnergy());
+ aida.cloud2D("VisE Clus vs Nhits in HAD B NN Cluster").fill(hbcluster.getCalorimeterHits().size(),hbcluster.getEnergy());
+ }
+
+ // PHOTON-Finder here (use this QAD one for now)
+ // check these clusters and try to ID photons - look for track-cluster matches at shower max
+ // require that EM cluster has hits at this layer.
+ aida.cloud1D("Num EM B NNClusters").fill(ebclusters.size());
+ List<Cluster> phoclusters = new ArrayList();
+ List<Cluster> ebshowclusters = new ArrayList();
+ for (Cluster ebcluster : ebclusters)
+ {
+ aida.cloud1D("Energy of EM B NNCluster").fill(ebcluster.getEnergy()/_sfEMB);
+ aida.cloud1D("Nhits in EM B NNCluster").fill(ebcluster.getCalorimeterHits().size());
+ aida.cloud2D("E Clus vs Nhits in EM B NNCluster").fill(ebcluster.getCalorimeterHits().size(),ebcluster.getEnergy()/_sfEMB);
+ // determine E-weighted position of each cluster
+ double ebhX=0; double ebhY=0; double ebhZ=0;
+ double ebclR=0; double ebclth=0; double ebclph=0;
+ for (int i=0; i<ebcluster.getCalorimeterHits().size(); ++i)
+ {
+ CalorimeterHit ebhit = ebcluster.getCalorimeterHits().get(i);
+ double[] ebhpos = ebhit.getPosition();
+ ebhX += ebhpos[0]*ebhit.getEnergy()/_sfEMB;
+ ebhY += ebhpos[1]*ebhit.getEnergy()/_sfEMB;
+ ebhZ += ebhpos[2]*ebhit.getEnergy()/_sfEMB;
+ }
+ ebhX = ebhX/(ebcluster.getEnergy()/_sfEMB);
+ ebhY = ebhY/(ebcluster.getEnergy()/_sfEMB);
+ ebhZ = ebhZ/(ebcluster.getEnergy()/_sfEMB);
+ ebclR = Math.sqrt(ebhX*ebhX+ebhY*ebhY);
+ aida.cloud1D("R Pos of EM B NNCluster").fill(ebclR/1000);
+ aida.cloud1D("Z Pos of EM B NNCluster").fill(ebhZ/1000);
+ aida.cloud2D("R vs Z of EM B NNCluster").fill(ebhZ/1000,ebclR/1000);
+ ebclth = Math.atan(ebclR/ebhZ);
+ if(ebclth<0) ebclth+=Math.PI;
+ ebclph = Math.atan2(ebhY,ebhX);
+ if(ebclph<0) ebclph+=2*Math.PI;
+ aida.cloud1D("EM B NNCluster Theta").fill(ebclth);
+ aida.cloud1D("EM B NNCluster Phi").fill(ebclph);
+ aida.cloud2D("EM B NNCluster Theta vs Phi").fill(ebclth,ebclph);
+ // Now, check Track-NNCluster distance for match at track swim layer 8 (shower max?)
+ double trcldelth=0; double trcldelph=0; double trcldist=0;
+ int ntrclmtch=0;
+ for (int i=0; i<ntrkswim; ++i)
+ {
+ trcldelth = Math.abs(ebclth-trbtheta[i][7]);
+ aida.cloud1D("Cluster-Track theta diff").fill(trcldelth);
+ trcldelph = Math.abs(ebclph-trbphi[i][7]);
+ if (trcldelph > Math.PI) trcldelph = 2*Math.PI-trcldelph;
+ aida.cloud1D("Cluster-Track phi diff").fill(trcldelph);
+ trcldist = Math.sqrt(trcldelth*trcldelth+trcldelph*trcldelph);
+ aida.cloud1D("Track-Cluster Distance").fill(trcldist);
+ if (trcldist<0.1) ntrclmtch++;
+ }
+ if (ntrclmtch == 0 && ebclR/1000 < 1.32 && ebcluster.getCalorimeterHits().size() > 4)
+ {
+ phoclusters.add(ebcluster);
+ } else
+ {
+ ebshowclusters.add(ebcluster);
+ }
+ }
+ if (phoclusters.size() > 0) event.put("PhoBClusters",phoclusters);
+ double PhoClusE = 0;
+ for (Cluster phocluster : phoclusters)
+ {
+ PhoClusE += phocluster.getEnergy()/_sfEMB;
+ }
+ aida.cloud1D("Num Photon B Clusters").fill(phoclusters.size());
+ if (ebshowclusters.size() > 0) event.put("EBShowClusters",ebshowclusters);
+ aida.cloud1D("Num EM B Shower Clusters").fill(ebshowclusters.size());
+
+ // make a map linking hits to EBShower clusters
+ Map<CalorimeterHit, Cluster> ebshowhitclus = new HashMap<CalorimeterHit, Cluster>();
+ for (Cluster ebshowcluster : ebshowclusters)
+ {
+ for (int i=0; i<ebshowcluster.getCalorimeterHits().size(); ++i)
+ {
+ // get CAL hits and put in map
+ CalorimeterHit ebsh = ebshowcluster.getCalorimeterHits().get(i);
+ ebshowhitclus.put(ebsh,ebshowcluster);
+ }
+ }
+
+ // Now that mips and photons are found, try to link Mip clusters, EBShower Clusters and HADBNN
+ // clusters with tracks to form the complete charged hadron shower.
+ // Start with the tracks
+ List<BasicCluster> trkcalclusters = new ArrayList();
+ List<BasicCluster> neucalclusters = new ArrayList();
+
+// for (int i=0; i<ntrkswim; ++i)
+ for (Integer i : itracks)
+ {
+ BasicCluster tcclus = new BasicCluster();
+ // first, check for track-cell match with a mip cluster at last mip layer (IntLay[i]-1)
+ if (IntLay[i] > 0) tcclus.addCluster(trmipclusmap.get(i));
+// for (BasicCluster calbmipcluster : calbmipclusters)
+// {
+// if (IntLay[i] == 0) continue;
+ // compare cluster hits with track position
+// List<CalorimeterHit> clhits = calbmipcluster.getCalorimeterHits();
+// int ntc=0;
+// double tcd=99;
+// for (CalorimeterHit clhit : clhits)
+// {
+ // compare hit-track distance
+// double[] clhpos = clhit.getPosition();
+// double clhR = Math.sqrt(clhpos[0]*clhpos[0]+clhpos[1]*clhpos[1]);
+// double clhth = Math.atan(clhR/clhpos[2]);
+// if(clhth<0) clhth+=Math.PI;
+// double clhph = Math.atan2(clhpos[1],clhpos[0]);
+// if(clhph<0) clhph+=2*Math.PI;
+// double tcdelth = Math.abs(clhth-trbtheta[i][IntLay[i]-1]);
+// double tcdelph = Math.abs(clhph-trbphi[i][IntLay[i]-1]);
+// if (tcdelph > Math.PI) tcdelph = 2*Math.PI-tcdelph;
+// tcd = Math.sqrt(tcdelth*tcdelth+tcdelph*tcdelph);
+// aida.cloud1D("Track-Cell Distance in Track Mip Cluster").fill(tcd);
+// if (tcd<0.01) ntc++;
+// }
+// aida.cloud1D("Number of cell matches per cluster").fill(ntc);
+// if (ntc > 0) tcclus.addCluster(calbmipcluster);
+// }
+ for (Cluster ebshowcluster : ebshowclusters)
+ {
+ List<CalorimeterHit> ebshowhits = ebshowcluster.getCalorimeterHits();
+ int ntsh=0;
+ double tsd=99;
+ for (CalorimeterHit ebshowhit : ebshowhits)
+ {
+ // compare hit-track distance
+ double[] shhpos = ebshowhit.getPosition();
+ double shR = Math.sqrt(shhpos[0]*shhpos[0]+shhpos[1]*shhpos[1]);
+ double shhth = Math.atan(shR/shhpos[2]);
+ if(shhth<0) shhth+=Math.PI;
+ double shhph = Math.atan2(shhpos[1],shhpos[0]);
+ if(shhph<0) shhph+=2*Math.PI;
+ double tsdelth = Math.abs(shhth-trbtheta[i][IntLay[i]]);
+ double tsdelph = Math.abs(shhph-trbphi[i][IntLay[i]]);
+ if (tsdelph > Math.PI) tsdelph = 2*Math.PI-tsdelph;
+ tsd = Math.sqrt(tsdelth*tsdelth+tsdelph*tsdelph);
+ aida.cloud1D("Track-Cell Distance in Track EM Shower Cluster").fill(tsd);
+ if (tsd<0.175) ntsh++;
+ }
+ aida.cloud1D("Number of cell matches per EM shower").fill(ntsh);
+ if (ntsh > 0) tcclus.addCluster(ebshowcluster);
+ }
+ for (Cluster hbcluster : hbclusters)
+ {
+ List<CalorimeterHit> hbhits = hbcluster.getCalorimeterHits();
+ int nthsh=0;
+ double thsd=99;
+ for (CalorimeterHit hbhit : hbhits)
+ {
+ // compare hit-track distance
+ double[] hhpos = hbhit.getPosition();
+ double hR = Math.sqrt(hhpos[0]*hhpos[0]+hhpos[1]*hhpos[1]);
+ double hhth = Math.atan(hR/hhpos[2]);
+ if(hhth<0) hhth+=Math.PI;
+ double hhph = Math.atan2(hhpos[1],hhpos[0]);
+ if(hhph<0) hhph+=2*Math.PI;
+ double thdelth = Math.abs(hhth-trbtheta[i][_numbemlayers+1]);
+ double thdelph = Math.abs(hhph-trbphi[i][_numbemlayers+1]);
+ if (thdelph > Math.PI) thdelph = 2*Math.PI-thdelph;
+ thsd = Math.sqrt(thdelth*thdelth+thdelph*thdelph);
+ aida.cloud1D("Track-shower Distance in Track HAD Shower Cluster").fill(thsd);
+ if (thsd<0.175) nthsh++;
+ }
+ aida.cloud1D("Number of cell matches per HAD shower").fill(nthsh);
+ if (nthsh > 0) tcclus.addCluster(hbcluster);
+ }
+ trkcalclusters.add(tcclus);
+ }
+ if (trkcalclusters.size() > 0) event.put("TrClMatClusters",trkcalclusters);
+
+
+ // get some MC info to compare
+ int npho=0; int nneu=0; int nchar=0;
+ double phoE = 0; double neuE = 0;
+ List<MCParticle> mcobs = event.getMCParticles();
+ for (MCParticle mcob : mcobs)
+ {
+ // must be final state particle
+// if (mcob.getGeneratorStatus() != mcob.FINAL_STATE) continue;
+ if (mcob.getGeneratorStatus() == mcob.DOCUMENTATION || mcob.getGeneratorStatus() == 0) continue;
+// System.out.println("Generator Status " +mcob.getGeneratorStatus());
+ // must be charged particle
+ int iqob = (int)mcob.getCharge();
+ double mcppx = mcob.getPX();
+ double mcppy = mcob.getPY();
+ double mcppz = mcob.getPZ();
+ double mcpE = mcob.getEnergy();
+ double mcpmt = Math.sqrt(mcppx*mcppx+mcppy*mcppy+mcppz*mcppz);
+ double mcpct = mcppz/mcpmt;
+ double mcppt = Math.sqrt(mcppx*mcppx+mcppy*mcppy);
+ if (Math.abs(mcpct)>0.8) continue;
+ if (iqob == 0)
+ {
+ int pid = Math.abs(mcob.getType().getPDGID());
+ // if photon, count up number, E
+ if (pid == 22)
+ {
+ npho++;
+ phoE += mcpE;
+ }
+ // if neutron or Klong, count up
+ if (pid == 2112 | pid==130 || pid==310 || pid==311)
+ {
+ nneu++;
+ neuE += mcpE;
+ aida.cloud1D("Neutral E Distribution").fill(mcpE);
+ }
+ } else
+ {
+ // its charged
+ if (mcppt > 0.9525) nchar++;
+ }
+ }
+ aida.cloud1D("Number of MC Photons").fill(npho);
+ aida.cloud1D("True Photon E").fill(phoE);
+ aida.cloud1D("Photon Cluster E").fill(PhoClusE);
+ aida.cloud1D("True Photon E - Photon Cluster E").fill(phoE-PhoClusE);
+ aida.cloud1D("Number of True Photons - Photon Clusters").fill(npho-phoclusters.size());
+ aida.cloud1D("Number of Neutrals").fill(nneu);
+ aida.cloud1D("Neutral E Sum").fill(neuE);
+ aida.cloud1D("Number of charged Particles").fill(nchar);
+ aida.cloud2D("Num MC Photons vs Num Photon Clusters").fill(npho,phoclusters.size());
+ aida.cloud2D("Num MC charged Particles vs Num EM B Mip Clusters").fill(ntrkswim,embmipclusters.size());
+ aida.cloud2D("Num MC charged Particles vs Num EM B Mip Clusters").fill(nchar,embmipclusters.size());
+ aida.cloud2D("Num Swum charged Particles vs Num EM B Mip Clusters").fill(ntrkswim,embmipclusters.size());
+ aida.cloud2D("Num MC charged Particles vs Num Swum charged").fill(nchar,ntrkswim);
+
+ // end of everything in event
+ }
+}
lcsim/src/org/lcsim/contrib/compile/SteveMagill
diff -N Helix.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Helix.java 23 Jul 2005 01:49:56 -0000 1.1
@@ -0,0 +1,248 @@
+import java.util.*;
+import hep.physics.vec.*;
+
+/**
+ *
+ * Given a particle of initial momentum, position and
+ * charge (p, r0, and iq), calculates its trajectory in a
+ * uniform magnetic field in the z-direction and determines
+ * its point of intersection (r_hit) with a cylinder centered at
+ * the origin with axis parallel to the z-axis and with ends.
+ * Trajectories are parameterized as fn of path length, s.
+ * There are no approximations made.
+ *
+ * Units are Tesla, mm, GeV/c
+ * @author Ray Frey - Javaized by Tony Johnson
+ * modified by S Magill for org.lcsim
+ * changed from cm to mm (LCIO units)
+ * added swimmer for endcap geometry
+ */
+public class Helix
+{
+
+ /**
+ * Create a Helix swimmer for symmetric (around z=0) Barrel
+ * @param B field strength; uniform, solenoidal, directed along z-axis
+ * @param r_cyl radius of cylinder
+ * @param z_cyl z-coordinate at end of cylinder (assumes symmetric placement at +/- z )
+ */
+ public Helix(double B, double r_cyl, double z_cyl)
+ {
+ this.B = B;
+ this.r_cyl = r_cyl/1000;
+ this.z_cyl = z_cyl/1000;
+ }
+ /**
+ * Create a Helix swimmer for EndCaps (at fixed +/-z and symmetric around r=0)
+ * @param B field strength; uniform, solenoidal, directed along z-axis
+ * @param r_min min radius of plane
+ * @param r_max max radius of plane
+ * @param z_pl z-coordinate of planes (assumes symmetric placement zround r=0)
+ */
+ public Helix(double B, double r_min, double r_max, double z_pl)
+ {
+ this.B = B;
+ this.r_min = r_min/1000;
+ this.r_max = r_max/1000;
+ this.z_pl = z_pl/1000;
+ }
+ /**
+ * @return position at point of intersection (x,y,z)
+ */
+ public double[] getIntersect()
+ {
+ return r_hit;
+ }
+ /**
+ * @return cyl/plane which is first hit
+ */
+ public int getPlane()
+ {
+ return i_hit;
+ }
+ public double[] getCenter()
+ {
+ double[] result = { xc*1000, yc*1000 };
+ return result;
+ }
+ public double getRadius()
+ {
+ return R*1000;
+ }
+ public final static int PROBLEM = 0;
+ public final static int PLANE_R = 1;
+ public final static int PLANE_PLUSZ = 2;
+ public final static int PLANE_MINUSZ = 3;
+ public final static int PLANE_Z = 4;
+ public final static int PLANE_RMIN = 5;
+ public final static int PLANE_RMAX = 6;
+
+ /**
+ * @param p 3-momentum (px,py,pz)
+ * @param r0_in initial position (x0,y0,z0)
+ * @param iq charge iq = q/|e| = +/- 1
+ * @param isFlipped curvature flipping option
+ * @see #setCurvatureFlip
+ */
+ public void swim(Hep3Vector p , Hep3Vector r0_in, int iq,
+ boolean isFlipped)
+ {
+ setCurvatureFlip(isFlipped);
+ swim(p, r0_in, iq);
+ }
+
+ /**
+ * @param p 3-momentum (px,py,pz)
+ * @param r0_in initial position (x0,y0,z0)
+ * @param iq charge iq = q/|e| = +/- 1
+ */
+ public void swim(Hep3Vector p , Hep3Vector r0_in, int iq)
+ {
+
+
+ double pt, s, q, pmom,
+ sin_lambda, cos_lambda, sin_phi0, cos_phi0,
+ temp, szhit, srhit, sin_lam_min;
+ double Rc, phic, phi0, Rq, darg, diff, tdiff;
+ int i;
+ i_hit = PROBLEM;
+ Hep3Vector r0 = new BasicHep3Vector(r0_in.x()/1000., r0_in.y()/1000., r0_in.z()/1000.);
+// r0[0] = r0_in[0]/1000;
+// r0[1] = r0_in[1]/1000;
+// r0[2] = r0_in[2]/1000;
+//
+// calculate useful quantities: lambda is dip angle (w/rt x-y plane;
+// phi0 is initial azimuthal angle in x-y plane w/rt x-axis;
+// R is radius of circle in x-y plane;
+// Rc**2=xc**2 + yc**2 is distance to center of circle in x-y plane
+//
+ temp = p.x()*p.x() + p.y()*p.y();
+ pt = Math.sqrt(temp);
+ pmom = Math.sqrt(temp + p.z()*p.z());
+ sin_lambda = p.z()/pmom;
+ cos_lambda = pt/pmom;
+ cos_phi0 = p.x()/pt;
+ sin_phi0 = p.y()/pt;
+ phi0 = Math.atan2(p.y(),p.x());
+//
+// First, check intercept with planes at z=+/- z_cyl
+// Eqn of motion z(s)=z0 + s*sin(lambda)
+// Check +z_cyl plane:
+//
+ szhit = sbig;
+ sin_lam_min = z_cyl/sbig;
+ if (sin_lambda > sin_lam_min)
+ {
+ i_hit = PLANE_PLUSZ;
+ szhit = ( z_cyl - r0.z() )/sin_lambda;
+ }
+// likewise check -z_cyl plane
+ else if (sin_lambda < -sin_lam_min)
+ {
+ i_hit = PLANE_MINUSZ;
+ szhit = (-z_cyl - r0.z() )/sin_lambda;
+ }
+//
+// Now check intercept of x-y plane path with cylinder
+//
+// First, check that pt is finite
+ if (pt < pt_tiny)
+ {
+ r_hit[0] = 0.;
+ r_hit[1] = 0.;
+ r_hit[2] = 1000*(r0.z() + szhit*sin_lambda);
+ return;
+ }
+//
+// calculate center of circle and circle radius=R
+//
+ q = -iq*m_flip; // this should allow for full flipping
+ R = pt/(0.3*B);
+ Rq = R/q;
+ xc = r0.x() + (Rq)*sin_phi0;
+ yc = r0.y() - (Rq)*cos_phi0;
+ Rc = Math.sqrt( xc*xc + yc*yc );
+ phic = Math.atan2(yc,xc);
+//
+// now calculate path length at intersection point
+//
+ darg = r_cyl*r_cyl/(2.*Rq*Rc) - Rc/(2.*Rq) - Rq/(2.*Rc);
+ if (Math.abs(darg) > 1.0)
+ {
+ srhit = sbig;
+ }
+ else
+ {
+ diff = Math.asin(darg) + phi0 - phic;
+ tdiff = Math.tan(diff);
+ srhit = (Rq/cos_lambda)*Math.atan(tdiff);
+ }
+//
+// Determine if cyl or z-plane hit first (smaller path length)
+//
+ if (Math.abs(srhit) <= szhit)
+ {
+ s = srhit;
+ i_hit = PLANE_R;
+ }
+ else
+ {
+ s = szhit;
+ }
+ //
+ // Calculate position at this intercept point
+ //
+ darg = s*cos_lambda/Rq - phi0;
+ r_hit[0] = 1000*(xc + Rq*Math.sin( darg ));
+ r_hit[1] = 1000*(yc + Rq*Math.cos( darg ));
+ r_hit[2] = 1000*(r0.x() + s*sin_lambda);
+ }
+
+ /**
+ * Flag for flipped curvature. <br>
+ * This enables us to compare this Helix with GISMO hits
+ * for example.
+ * @param isFlipped true if curvatures should be flipped
+ * (default: with flipping)
+ * @see #isFlipped
+ */
+ public void setCurvatureFlip(boolean isFlipped)
+ {
+ m_flip = isFlipped ? -1 : 1;
+ }
+
+ /**
+ * Status of curvature flipping
+ * @return true if curvature flipping option is on
+ * @see #setCurvatureFlip
+ */
+ public boolean isFlipped()
+ {
+ return (m_flip < 0);
+ }
+
+ //
+ //
+ //
+ private double B;
+ private double r_cyl;
+ private double z_cyl;
+ private double r_min;
+ private double r_max;
+ private double z_pl;
+ private double xc, yc;
+ private double R;
+
+ // for curvature flipping
+ private int m_flip = 1; // is NOT flipped.
+
+ private static final double sbig = 99999.;
+ private static final double pt_tiny = 0.010;
+
+ private double[] r_hit = new double[3];
+ private int i_hit;
+}
+
+
+
+
