lcsim/src/org/lcsim/contrib/compile/SteveMagill
diff -N BTrMipClus.java
--- BTrMipClus.java 9 Aug 2005 23:16:07 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,985 +0,0 @@
-// package org.lcsim.contrib.compile.SteveMagill;
-
-// Overall PFA algorithm includes 1) Track IL-finder in Barrel CAL, creation of EM Mip
-// Clusters, HAD Mip Clusters, CAL Mip Clusters; 2) Nearest-Neighbor clustering of remaining
-// hits, EM NN Clusters and AD NN CLusters; 3) Photon ID of EM NN Clusters, creation of
-// EM Photon Clusters, EM Shower Clusters; 4) Association of EM and HAD Mip Clusters, EM Shower
-// Clusters, and HAD NN Shower Clusters to form Track-CAL Associated Shower Clusters
-// import java.util.List;
-import java.util.*;
-import hep.physics.vec.Hep3Vector;
-// import org.lcsim.event.Cluster;
-// import org.lcsim.event.EventHeader;
-import org.lcsim.event.*;
-import org.lcsim.event.EventHeader.LCMetaData;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-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;
- private double _neupi;
-// 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.0120; // for now this is tested sampling fraction for W/Si EMCAL in SiD Detector
- _sfEMB = 0.0117; // for now, this is effective sf for 20/10 Si/W ECAL from 2 GeV gammas
-// _sfHB = 0.000025; // for now this is tested sampling fraction for SS/RPC HCAL in SiD Detector
- // but for now, use calibration as a sampling fraction (then HAD E in a cluster is OK)
-// _sfHB = 5.3; // 200 MeV per hit in HCAL for SiD SSRPC, 10.0 for pions
- _sfHB = 26.0; // 38 MeV per hit in HCAL for CDC WScin after 1/3 mip thresold cut, 30 for pions
- _neupi = 0.87; // ratio of nhits neu to pi - .87 for CDC, .57 for SiD
-// _sfHB = 0.06; // old SDJan03 SS/Scin analog HCAL sf
-
- _initialized = true;
- } // end of initialization section
-
- // get calorimeter hits by layer - need to calculate hit density for cells
- double EMBESum=0; double EMBESumC=0; double HADBESum=0;
- double EEMlowE=0; double EEMhighE=0;
- double ESumBTotCAL=0;
- double hthrshld = 0;
- double htime = 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);
- EMBESum += embhit.getRawEnergy();
- // for CDC, use this for calibrated B EM ESum
- if (emblayer < 20)
- {
- EMBESumC += embhit.getRawEnergy()/_sfEMB;
- EEMlowE += embhit.getRawEnergy();
- }
- else
- {
- EMBESumC += embhit.getRawEnergy()/(_sfEMB/2);
- EEMhighE += embhit.getRawEnergy();
- }
- 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>();
-
- if (event.getDetectorName().equals("cdcaug05"))
- {
- hthrshld = 0.0005; // 1/3 mip threshold for cdc Scin HCAL
- htime = 100; // 100 ns cut on hit time
-// System.out.println("HAD Threshold = " +hthrshld);
- }
- for (SimCalorimeterHit hadbhit : hadbHits)
- {
- long hitID = hadbhit.getCellID();
- hadbdecoder.setID(hitID);
- int hadblayer = hadbdecoder.getLayer();
- aida.cloud1D("HAD Cell Vis E before thr cut").fill(hadbhit.getRawEnergy());
- aida.cloud1D("HAD Cell Time before timing cut").fill(hadbhit.getTime());
- if (hadbhit.getRawEnergy() > hthrshld && hadbhit.getTime() < 100)
- {
- aida.cloud1D("HAD Cell Vis E after thr cut").fill(hadbhit.getRawEnergy());
- 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);
- HADBESum += hadbhit.getRawEnergy();
- 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(EMBESumC,nhitsEMB);
- aida.cloud2D("NHits HAD B vs HAD B Vis ESum").fill(HADBESum,nhitsHADB);
-// aida.cloud1D("Pho E first 20").fill(EEMlowE);
-// aida.cloud1D("Pho E last 10").fill(EEMhighE);
- if (EMBESumC < 0.40)
- {
- aida.cloud1D("EM B Calib ESum, mips only").fill(EMBESumC);
- aida.cloud1D("NHits EM B, mips only").fill(nhitsEMB);
- aida.cloud1D("NHits HAD B, EM B mips only").fill(nhitsHADB);
- aida.cloud1D("HAD B ESum, EM mips only").fill(nhitsHADB/_sfHB);
- }
- 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);
- aida.cloud1D("HAD B Cal ESum, no EM B Hits").fill(nhitsHADB/_sfHB);
- }
- 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(EMBESumC+nhitsHADB/_sfHB);
- aida.cloud1D("Total EM B Vis ESum").fill(EMBESum);
- aida.cloud1D("Total EM B Calib ESum").fill(EMBESumC);
- aida.cloud1D("Total HAD B Vis ESum").fill(HADBESum);
-
- // Make a new collection that only includes had hits above threshold - for Event display
- List<CalorimeterHit> hlist = new ArrayList(hadbhitmap.values());
- LCMetaData meta = event.getMetaData(hadbHits);
- event.put("ThrHcalBHits",hlist,CalorimeterHit.class,meta.getFlags(),meta.getName());
-
- /** 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 track to mip clusters
- Map<Integer, BasicCluster> trmipclusmap = new HashMap<Integer, BasicCluster>();
- Map<Integer, BasicCluster> tremipmap = new HashMap<Integer, BasicCluster>();
- Map<Integer, BasicCluster> trhmipmap = 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] = _numbemlayers+_numbhadlayers; // 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
- 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())
- {
- if (h == null) continue;
- 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())
- {
- if (eh == null) continue;
- 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())
- {
- if (hh == null) continue;
- hbmclus.addHit(hh);
- bmipclus.addHit(hh);
- }
- }
- hadbmipclusters.add(hbmclus);
- }
- calbmipclusters.add(bmipclus);
- }
- trmipclusmap.put(i, bmipclus);
- tremipmap.put(i,ebmclus);
- trhmipmap.put(i,hbmclus);
- } 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]);
- }
- }
- }
- if (j == 63 && nointeraction)
- {
- if (ebclhitmap.size() > 0)
- {
- for (CalorimeterHit eh : ebclhitmap.values())
- {
- ebmclus.addHit(eh);
- bmipclus.addHit(eh);
- }
- }
- embmipclusters.add(ebmclus);
- if (hbclhitmap.size() > 0)
- {
- for (CalorimeterHit hh : hbclhitmap.values())
- {
- hbmclus.addHit(hh);
- bmipclus.addHit(hh);
- }
- }
- hadbmipclusters.add(hbmclus);
- calbmipclusters.add(bmipclus);
- trmipclusmap.put(i,bmipclus);
- tremipmap.put(i,ebmclus);
- trhmipmap.put(i,hbmclus);
- }
- }
- }
- }
- // 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());
- }
- }
-
- // 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.getRawEnergy()/_sfEMB;
- ebhY += ebhpos[1]*ebhit.getRawEnergy()/_sfEMB;
- ebhZ += ebhpos[2]*ebhit.getRawEnergy()/_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;
- aida.cloud1D("Num Photon B Clusters").fill(phoclusters.size());
- for (Cluster phocluster : phoclusters)
- {
- PhoClusE += phocluster.getEnergy()/_sfEMB;
- }
- aida.cloud1D("Photon Cluster E").fill(PhoClusE);
-
- // 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();
- Map<Integer, BasicCluster> trcaclusmap = new HashMap<Integer, BasicCluster>();
- Map<Integer, Cluster> treshmap = new HashMap<Integer, Cluster>();
- Map<Integer, Cluster> trhshmap = new HashMap<Integer, Cluster>();
-
-// for (int i=0; i<ntrkswim; ++i)
-// System.out.println("Number of Tracks in List = " + itracks.size());
- for (Integer i : itracks)
- {
- BasicCluster tcclus = new BasicCluster();
- double trmipE = 0;
- double treshE = 0;
- double trhshE = 0;
- double trtotE = 0;
- int niter = 1;
-// aida.cloud1D("Track Momentum").fill(trbpp[i]);
- // first, check for track-cell match with a mip cluster at last mip layer (IntLay[i]-1)
-// System.out.println("IL = " + IntLay[i]);
- if (IntLay[i] == _numbemlayers+_numbhadlayers)
- {
- // don't do anything - all mips - no E/P test
- tcclus.addCluster(trmipclusmap.get(i));
- } else if (IntLay[i] < _numbemlayers+_numbhadlayers)
- {
- if (IntLay[i] > 0)
- {
- tcclus.addCluster(trmipclusmap.get(i));
- double emmipE = tremipmap.get(i).getEnergy()/_sfEMB;
- double hmipE = trhmipmap.get(i).getCalorimeterHits().size()/_sfHB;
- trmipE += emmipE+hmipE;
- aida.cloud1D("Mip E per Track").fill(trmipE);
- aida.cloud1D("Track E over P after adding mips").fill(trmipE/trbpp[i]);
- }
- do {
- for (Iterator<Cluster> ieclus = ebshowclusters.iterator(); ieclus.hasNext();)
- {
- Cluster eclus = ieclus.next();
- List<CalorimeterHit> ebshowhits = eclus.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);
- if (tsd < 0.1) aida.cloud1D("Track-Cell Distance in Track EM Shower Cluster").fill(tsd);
- if (tsd<niter*0.01) ntsh++;
- }
- aida.cloud1D("Number of cell matches per EM shower").fill(ntsh);
- if (ntsh > 0)
- {
- tcclus.addCluster(eclus);
- treshE += eclus.getEnergy()/_sfEMB;
- aida.cloud1D("Track E over P after adding EM shower").fill((trmipE+treshE)/trbpp[i]);
- treshmap.put(i, eclus);
- ieclus.remove();
- }
- }
- for (Iterator<Cluster> ihclus = hbclusters.iterator(); ihclus.hasNext();)
- {
- Cluster hclus = ihclus.next();
- List<CalorimeterHit> hbhits = hclus.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+5]);
- double thdelph = Math.abs(hhph-trbphi[i][_numbemlayers+5]);
- if (thdelph > Math.PI) thdelph = 2*Math.PI-thdelph;
- thsd = Math.sqrt(thdelth*thdelth+thdelph*thdelph);
- if (thsd < 0.1) aida.cloud1D("Track-shower Distance in Track HAD Shower Cluster").fill(thsd);
- if (thsd<niter*0.025) nthsh++;
- }
- aida.cloud1D("Number of cell matches per HAD shower").fill(nthsh);
- if (nthsh > 0)
- {
- tcclus.addCluster(hclus);
- trhshE += hclus.getCalorimeterHits().size()*_neupi/_sfHB;
-// System.out.println("Added " + hclus.getCalorimeterHits().size() + " hits to track match");
- aida.cloud1D("Track E over P after adding HAD shower").fill((trmipE+treshE+trhshE)/trbpp[i]);
- trhshmap.put(i, hclus);
- ihclus.remove();
- }
- }
-// System.out.println(" Iteration number = " + niter);
-// System.out.println("E over P ratio = " + (trmipE+treshE+trhshE)/trbpp[i]);
- if (niter == 4) break;
- niter++;
- } while ((trmipE+treshE+trhshE)/trbpp[i] < 0.95);
- }
- trtotE += trmipE+treshE+trhshE;
- if (trtotE > 0) aida.cloud1D("Track E over P").fill(trtotE/trbpp[i]);
- trkcalclusters.add(tcclus);
- trcaclusmap.put(i, tcclus);
- }
- if (trkcalclusters.size() > 0) event.put("TrClMatClusters",trkcalclusters);
-// aida.cloud1D("Number of Track/CAL Matches").fill(trkcalclusters.size());
- // if there are any EM shhower clusters left, put to event
- if (ebshowclusters.size() > 0) event.put("EBShowClusters",ebshowclusters);
- double EMShowClusE = 0;
- aida.cloud1D("Num EM B Shower Clusters").fill(ebshowclusters.size());
- for (Cluster ebshowcluster : ebshowclusters)
- {
- EMShowClusE += ebshowcluster.getEnergy()/_sfEMB;
- }
- aida.cloud1D("EM Shower Cluster E").fill(EMShowClusE);
-
- if (hbclusters.size() > 0) event.put("HADBNNClusters",hbclusters);
-
-
- // Compare E of Clusters to P of Track - this doesn't work - must only
- // be able to map one cluster to one track
-// for (Integer itrack : itracks)
-// {
-// double EMmipE = 0; double HmipE = 0; double EMshowE = 0; double HshowE = 0;
-// if (trcaclusmap.get(itrack) != null)
-// {
-// if (tremipmap.get(itrack) != null) EMmipE = tremipmap.get(itrack).getEnergy()/_sfEMB;
-// if (trhmipmap.get(itrack) != null) HmipE = trhmipmap.get(itrack).getCalorimeterHits().size()/(_sfHB*1.15);
-// if (treshmap.get(itrack) != null) EMshowE = treshmap.get(itrack).getEnergy()/_sfEMB;
-// if (trhshmap.get(itrack) != null) HshowE = trhshmap.get(itrack).getCalorimeterHits().size()/(_sfHB*1.15);
-// aida.cloud1D("EM mip E per track").fill(EMmipE);
-// aida.cloud1D("HAD mip E per track").fill(HmipE);
-// aida.cloud1D("Total CAL mip E per track").fill(EMmipE+HmipE);
-// aida.cloud1D("EM shower E per track").fill(EMshowE);
-// aida.cloud1D("HAD shower E per track").fill(HshowE);
-// aida.cloud1D("Total CAL shower E per track").fill(EMshowE+HshowE);
-// aida.cloud1D("Total CAL E per track").fill(EMshowE+HshowE+EMmipE+HmipE);
-// double trcalE = EMmipE+HmipE+EMshowE+HshowE;
-// double trcalP = trbpp[itrack];
-// if (trcalE/trcalP > 0) aida.cloud1D("E over P of TrackCAL match").fill(trcalE/trcalP);
-// }
-// }
-
- // get some MC info to compare
- int npho=0; int nneu=0; int nchar=0;
- double phoE = 0; double neuE = 0; double charE = 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 mckE = mcob.getEnergy()-mcob.getMass();
- 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);
- if (nhitsEMB < 1) aida.cloud2D("Neutral E vs NhadHits").fill(mckE,nhitsHADB);
- }
- } else
- {
- // its charged
- if (mcppt > 0.9525)
- {
- nchar++;
- charE += mcpE;
- }
- }
- }
- if (neuE+charE+phoE > 75)
- {
- aida.cloud1D("Num MC Photons").fill(npho);
- aida.cloud1D("True Photon E").fill(phoE);
- aida.cloud1D("True Photon E - Photon Cluster E").fill(phoE-PhoClusE);
- aida.cloud1D("Num True Photons - Photon Clusters").fill(npho-phoclusters.size());
- aida.cloud1D("Num Neutrals").fill(nneu);
- aida.cloud1D("Neutral E Sum").fill(neuE);
- if (nhitsEMB < 1 && nhitsHADB > 0)
- {
- aida.cloud1D("NHits HAD B over Neutral ESum").fill(nhitsHADB/(Math.sqrt(neuE*neuE-0.940*0.940)));
- aida.cloud2D("NHits HAD B vs Neutral ESum").fill(Math.sqrt(neuE*neuE-0.940*0.940),nhitsHADB);
- aida.cloud1D("NHits HAD B over Vis ESum").fill(nhitsHADB/HADBESum);
- aida.cloud2D("NHits HAD B vs Vis ESum").fill(HADBESum,nhitsHADB);
- aida.cloud1D("Kinetic ESum of neutrals").fill(Math.sqrt(neuE*neuE-0.940*0.940));
- }
- aida.cloud1D("Total Charged, Neutral and Photon ESum").fill(neuE+charE+phoE);
- aida.cloud1D("Total Barrrel ESum wcuts").fill(EMBESum/_sfEMB+nhitsHADB/_sfHB);
- aida.cloud1D("Num Swum MC Charged Particles").fill(itracks.size());
- aida.cloud1D("Charged Particle E Sum").fill(charE);
- aida.cloud2D("Num MC Photons vs Num Photon Clusters").fill(npho,phoclusters.size());
- aida.cloud2D("Num Swum MC Charged Particles vs Num Trk Cal Clusters").fill(itracks.size(),trkcalclusters.size());
- }
- // end of everything in event
- }
-}
lcsim/src/org/lcsim/contrib/compile/SteveMagill
diff -N Helix.java
--- Helix.java 23 Jul 2005 01:49:56 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,248 +0,0 @@
-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;
-}
-
-
-
-