lcsim/src/org/lcsim/contrib/LeiXia
diff -N rpc_PFA_np.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ rpc_PFA_np.java 21 Nov 2005 22:33:44 -0000 1.1
@@ -0,0 +1,2488 @@
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.*;
+import org.lcsim.util.swim.HelixSwim;
+
+/*
+ * Reconstruction: SiD study (rpc active medium) - non projective
+ */
+
+public class rpc_PFA_np extends Driver
+{
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ public void process(EventHeader event)
+ {
+
+// System.out.println();
+// System.out.println();
+ System.out.println();
+ System.out.println("event "+event.getEventNumber());
+
+//------------------
+//get MC truth first
+
+ List<MCParticle> PList = new ArrayList<MCParticle>(); // List of interesting particles
+
+// will save the following particles in this list (excluding neutrinos)
+// 1. ANY particles that appear on the front surface of em calorimeter
+// r = 1260.mm for barrel (real surface at r = 1270.mm)
+// z = 1670.mm for endcap (real surface at z = 1680.mm)
+// 2. or ANY charged FINAL STATE particle with P < 1.2GeV
+// charged particle with P < 0.953GeV will never get to barrel calorimeter
+// and may stop before reaching endcap
+// 3. ANY other particle that generate calorimeter hits
+// some particles have end point before reaching calorimeter (likely to be wrong ending point)
+// some hits were assigned to Final State partices decayed before reaching calorimeter
+// 1 and 2 will be added to the list from mc particle list
+// 3 will be added later when looking at calorimeter hits
+// NOTE: special care is needed to avoid double counting between 1,2 and 3.
+
+ List<MCParticle> EventMCList = event.getMCParticles();
+ double E_PList = 0.; // visible energy of particle 1,2 in PList
+ double E_nu = 0.; // neutrino energy (all neutrinos)
+ double E_nu_fs = 0.; // neutrino energy (final state neutrinos only)
+
+ int fg_barrel = 0; // label for barrel events (uds quark (P > 45 GeV) Pz/P < 0.707)
+
+ for (MCParticle mc : EventMCList)
+ {
+// other: 0 final_state: 1 Intermediate: 2 documentation: 3
+// there's 'unknown' type particle, with PDGID = 0, should get ride of them
+
+ if (((mc.getPDGID() == 1)||(mc.getPDGID() == 2)||(mc.getPDGID() == 3))&&(mc.getEnergy() > 45.))
+ {
+ if (Math.abs(mc.getPZ())/mc.getMomentum().magnitude() < 0.707)
+ {
+ fg_barrel = 1;
+ System.out.println("barrel event");
+ }
+ else
+ {
+ System.out.println("endcap event");
+ }
+ }
+
+ Hep3Vector mc_o = mc.getOrigin();
+ Hep3Vector mc_e = mc.getEndPoint();
+
+ double R_origin = Math.sqrt(mc_o.x()*mc_o.x() + mc_o.y()*mc_o.y());
+ double R_end = Math.sqrt(mc_e.x()*mc_e.x() + mc_e.y()*mc_e.y());
+
+// fill particles in my list
+ if (((R_origin < 1260.)&&(Math.abs(mc_o.z()) < 1670.))&&((R_end > 1260.)||(Math.abs(mc_e.z()) > 1670.)))
+ {
+// System.out.println("particle: "+mc.getType().getName()+" "+mc.getEnergy()+" "+mc.getGeneratorStatus());
+ if (mc.getGeneratorStatus() == 0) aida.cloud1D("R_origin of other particles").fill(R_origin);
+
+ int mcPDGID = mc.getPDGID();
+ if ((Math.abs(mcPDGID) == 12)||(Math.abs(mcPDGID) == 14)||(Math.abs(mcPDGID) == 16))
+ {
+ E_nu += mc.getEnergy();
+ if (mc.getGeneratorStatus() == 1) E_nu_fs += mc.getEnergy();
+ }
+ else if (mc.getPDGID() != 0)
+ {
+ PList.add(mc);
+ E_PList += mc.getEnergy();
+ }
+ }
+ else if ((Math.abs(mc.getCharge()) != 0)&&(mc.getMomentum().magnitude() < 1.2)&&(mc.getGeneratorStatus() == 1)&&(mc.getPDGID() != 0))
+ {
+// System.out.println("particle: "+mc.getType().getName()+" "+mc.getEnergy()+" "+mc.getGeneratorStatus());
+
+ PList.add(mc); // no neutrino in this case
+ E_PList += mc.getEnergy();
+ }
+
+ } // end of mc particle loop
+
+// print out the whole list, if wanted
+// for (MCParticle tt : PList) {System.out.println(" particle # "+PList.indexOf(tt)+" "+tt.getType().getName()+" "+tt.getPDGID()+" "+tt.getEnergy());}
+// System.out.println(" size of PList = "+PList.size()); // size = last index + 1
+
+ aida.cloud1D("E_PList").fill(E_PList);
+ aida.cloud1D("E_PList + E_nu").fill(E_PList + E_nu);
+ aida.cloud1D("E_PList + E_nu_fs").fill(E_PList + E_nu_fs);
+ aida.cloud1D("E_nu").fill(E_nu);
+ aida.cloud1D("E_nu_fs").fill(E_nu_fs);
+
+
+//---------------------------------
+//get EM hits, calculate d1-3, v1-3
+
+ List<SimCalorimeterHit> Barrel_em = event.getSimCalorimeterHits("EcalBarrHits");
+ List<SimCalorimeterHit> Endcap_em = event.getSimCalorimeterHits("EcalEndcapHits");
+// System.out.println("Ecal hits: Barrel, Endcap "+Barrel_em.size()+" "+Endcap_em.size());
+
+ int NHit_em = Barrel_em.size() + Endcap_em.size(); // total number of hits in Ecal = barrel + endcap
+ SimCalorimeterHit[] HitList_em = new SimCalorimeterHit[NHit_em]; // list of Ecal hits
+ int[] ShowerP_em = new int[NHit_em]; // number of contributing particles for each hit (in original data)
+ int[][] SPID_em = new int[NHit_em][10]; // list of contributing particles, ID = index in PList
+
+// fill this later, after ShowerP and SPID are available
+// int[] FinalP_em = new int[NHit_em]; // number of contributing particles for each hit (appropriate contributor for PFA)
+// int[][] FPID_em = new int[NHit_em][10]; // list of contributing particles, ID = index in PList
+ int c_em = 0; //counter of em hits
+
+ double[][] v1_em = new double[NHit_em][3];
+ double[][] v2_em = new double[NHit_em][3];
+ double[][] v3_em = new double[NHit_em][3];
+
+ int[] EM_ClusterID = new int[NHit_em]; // em hit cluster id, need to be cleared (=0) here
+ int[] EM_FinalID = new int[NHit_em]; // em hit cluster id, need to be cleared (=0) here
+
+ for (SimCalorimeterHit hit_em : Barrel_em) // fill barrel hits
+ {
+ HitList_em[c_em] = hit_em;
+ ShowerP_em[c_em] = hit_em.getMCParticleCount();
+// FinalP_em[c_em] = 0;
+
+ double[] x = hit_em.getPosition();
+
+ v3_em[c_em][0] = -x[0]/Math.sqrt(x[0]*x[0]+x[1]*x[1]);
+ v3_em[c_em][1] = -x[1]/Math.sqrt(x[0]*x[0]+x[1]*x[1]);
+ v3_em[c_em][2] = 0;
+
+ v1_em[c_em][0] = 0;
+ v1_em[c_em][1] = 0;
+ v1_em[c_em][2] = 1;
+
+ v2_em[c_em][0] = v3_em[c_em][1];
+ v2_em[c_em][1] = -v3_em[c_em][0];
+ v2_em[c_em][2] = 0;
+// System.out.println("x,y,z = "+v3_em[c_em][0]+" "+v3_em[c_em][1]+" "+v3_em[c_em][2]);
+
+ EM_ClusterID[c_em] = -1;
+ EM_FinalID[c_em] = -1;
+
+ c_em += 1;
+ }
+
+ for (SimCalorimeterHit hit_em : Endcap_em) // fill endcap hits
+ {
+ HitList_em[c_em] = hit_em;
+ ShowerP_em[c_em] = hit_em.getMCParticleCount();
+// FinalP_em[c_em] = 0;
+
+ double[] x = hit_em.getPosition();
+
+ v1_em[c_em][0] = 1;
+ v1_em[c_em][1] = 0;
+ v1_em[c_em][2] = 0;
+
+ v2_em[c_em][0] = 0;
+ v2_em[c_em][1] = 1;
+ v2_em[c_em][2] = 0;
+
+ v3_em[c_em][0] = 0;
+ v3_em[c_em][1] = 0;
+ v3_em[c_em][2] = -x[2]/Math.abs(x[2]);
+
+ EM_ClusterID[c_em] = -1;
+ EM_FinalID[c_em] = -1;
+
+ c_em += 1;
+ }
+// System.out.println("c_em and NHit_em: "+c_em+" "+NHit_em);
+
+ for (int i=0; i < NHit_em; i++) // fill SP information
+ {
+
+ MCParticle loop_mc;
+ for (int j=0; j < ShowerP_em[i]; j++) // loop over contributing particles
+ {
+ loop_mc = HitList_em[i].getMCParticle(j);
+ Hep3Vector p_o = loop_mc.getOrigin();
+ Hep3Vector p_e = loop_mc.getEndPoint();
+ double r_o = Math.sqrt(p_o.x()*p_o.x() + p_o.y()*p_o.y());
+ double r_e = Math.sqrt(p_e.x()*p_e.x() + p_e.y()*p_e.y());
+
+ int Found = 0;
+ int index = 0;
+ for (MCParticle mc : PList)
+ {
+ if ((loop_mc.getPDGID()==mc.getPDGID())&&((Math.abs(loop_mc.getEnergy()-mc.getEnergy())/mc.getEnergy())<0.00001))
+ {
+ if ((Math.abs(loop_mc.getEndPoint().x()-mc.getEndPoint().x())<0.00001)&&(Math.abs(loop_mc.getEndPoint().y()-mc.getEndPoint().y())<0.00001))
+ {
+ Found += 1;
+ index = PList.indexOf(mc);
+ }
+ }
+ }
+
+ if (Found == 1)
+ {
+ SPID_em[i][j] = index;
+ }
+ else if (Found == 0)
+ {
+ if (loop_mc.getGeneratorStatus() == 1)
+ {
+// aida.cloud1D("origin r: FS particle").fill(r_o);
+// aida.cloud1D("end r: FS particle").fill(r_e);
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+
+ int pFound = 0;
+ MCParticle tt = loop_mc;
+ while ((pFound == 0)&&(tt.getOrigin().magnitude() > 20.)) // stop searching while found one or parent is within 20.mm from IP
+ {
+ List<MCParticle> parents = tt.getParents();
+ if (parents.size() > 1) System.out.println("Attn: more than one parents found before finding SP: N = "+parents.size());
+ for (MCParticle p_mc : parents)
+ {
+ tt = p_mc;
+ }
+
+ for (MCParticle mc : PList)
+ {
+ if ((tt.getPDGID()==mc.getPDGID())&&((Math.abs(tt.getEnergy()-mc.getEnergy())/mc.getEnergy())<0.00001))
+ {
+ if ((Math.abs(tt.getEndPoint().x()-mc.getEndPoint().x())<0.00001)&&(Math.abs(tt.getEndPoint().y()-mc.getEndPoint().y())<0.00001))
+ {
+ pFound += 1;
+ SPID_em[i][j] = PList.indexOf(mc);
+ }
+ }
+ }
+ } // searching parents loop
+
+ if (pFound == 0)
+ {
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+ PList.add(loop_mc);
+ SPID_em[i][j] = PList.size() - 1; // last particle in list
+ }
+ }
+ else if (loop_mc.getGeneratorStatus() == 0)
+ {
+// aida.cloud1D("origin r: other").fill(r_o);
+// aida.cloud1D("end r: other").fill(r_e);
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+
+ int pFound = 0;
+ MCParticle tt = loop_mc;
+ while ((pFound == 0)&&(tt.getOrigin().magnitude() > 20.)) // stop searching while found one or parent is within 20.mm from IP
+ {
+ List<MCParticle> parents = tt.getParents();
+ if (parents.size() > 1) System.out.println("Attn: more than one parents found before finding SP: N = "+parents.size());
+ for (MCParticle p_mc : parents)
+ {
+ tt = p_mc;
+ }
+
+ for (MCParticle mc : PList)
+ {
+ if ((tt.getPDGID()==mc.getPDGID())&&((Math.abs(tt.getEnergy()-mc.getEnergy())/mc.getEnergy())<0.00001))
+ {
+ if ((Math.abs(tt.getEndPoint().x()-mc.getEndPoint().x())<0.00001)&&(Math.abs(tt.getEndPoint().y()-mc.getEndPoint().y())<0.00001))
+ {
+ pFound += 1;
+ SPID_em[i][j] = PList.indexOf(mc);
+ }
+ }
+ }
+ } // searching parents loop
+
+ if (pFound == 0)
+ {
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+ if (loop_mc.getPDGID() != 0)
+ {
+ PList.add(loop_mc);
+ SPID_em[i][j] = PList.size() - 1; // last particle in list
+ }
+ else
+ {
+ while (loop_mc.getPDGID() == 0)
+ {
+ loop_mc = loop_mc.getParents().get(0);
+ }
+ PList.add(loop_mc);
+ SPID_em[i][j] = PList.size() - 1; // last particle in list
+ }
+ }
+ }
+ else // so far, no particles found with a type other than 'FS' or 'other'
+ {
+ aida.cloud1D("origin r: rest").fill(r_o);
+ aida.cloud1D("end r: rest").fill(r_e);
+ }
+ }
+ else
+ {
+ System.out.println(" Something went wrong: Found = "+Found+" index = "+index);
+ }
+
+ } // end of contributing particles loop
+
+ } // end of fill SP loop
+
+
+//----------------------------------
+//get HAD hits, calculate d1-3, v1-3
+
+ List<SimCalorimeterHit> Barrel_had = event.getSimCalorimeterHits("HcalBarrHits");
+ List<SimCalorimeterHit> Endcap_had = event.getSimCalorimeterHits("HcalEndcapHits");
+// System.out.println("Hcal hits: Barrel, Endcap "+Barrel_had.size()+" "+Endcap_had.size());
+
+ int NHit_had = Barrel_had.size() + Endcap_had.size(); // total number of hits in Hcal = barrel + endcap
+ SimCalorimeterHit[] HitList_had = new SimCalorimeterHit[NHit_had]; // list of Hcal hits
+ int[] ShowerP_had = new int[NHit_had]; // number of contributing particles for each hit
+ int[][] SPID_had = new int[NHit_had][10]; // list of contributing particles, ID = index in PList
+
+// fill this later, after ShowerP and SPID are available
+// int[] FinalP_had = new int[NHit_had]; // number of contributing FinalState particles for each hit
+// int[][] FPID_had = new int[NHit_had][10]; // list of contributing FPs, ID = FP index in MyMCList
+ int c_had = 0; //counter of had hits
+
+ double[][] v1_had = new double[NHit_had][3];
+ double[][] v2_had = new double[NHit_had][3];
+ double[][] v3_had = new double[NHit_had][3];
+
+ int[] HAD_ClusterID = new int[NHit_had]; // had hit cluster id, need to be cleared (=0) here
+ int[] HAD_FinalID = new int[NHit_had]; // had hit final cluster id, need to be cleared (=0) here
+
+ for (SimCalorimeterHit hit_had : Barrel_had) // fill barrel hits
+ {
+ HitList_had[c_had] = hit_had;
+ ShowerP_had[c_had] = hit_had.getMCParticleCount();
+
+ double[] x = hit_had.getPosition();
+
+ v3_had[c_had][0] = -x[0]/Math.sqrt(x[0]*x[0]+x[1]*x[1]);
+ v3_had[c_had][1] = -x[1]/Math.sqrt(x[0]*x[0]+x[1]*x[1]);
+ v3_had[c_had][2] = 0;
+
+ v1_had[c_had][0] = 0;
+ v1_had[c_had][1] = 0;
+ v1_had[c_had][2] = 1;
+
+ v2_had[c_had][0] = v3_had[c_had][1];
+ v2_had[c_had][1] = -v3_had[c_had][0];
+ v2_had[c_had][2] = 0;
+
+ HAD_ClusterID[c_had] = -1;
+ HAD_FinalID[c_had] = -1;
+
+ c_had += 1;
+ }
+
+ for (SimCalorimeterHit hit_had : Endcap_had) // fill endcap hits
+ {
+ HitList_had[c_had] = hit_had;
+ ShowerP_had[c_had] = hit_had.getMCParticleCount();
+
+ double[] x = hit_had.getPosition();
+
+ v1_had[c_had][0] = 1;
+ v1_had[c_had][1] = 0;
+ v1_had[c_had][2] = 0;
+
+ v2_had[c_had][0] = 0;
+ v2_had[c_had][1] = 1;
+ v2_had[c_had][2] = 0;
+
+ v3_had[c_had][0] = 0;
+ v3_had[c_had][1] = 0;
+ v3_had[c_had][2] = -x[2]/Math.abs(x[2]);
+
+ HAD_ClusterID[c_had] = -1;
+ HAD_FinalID[c_had] = -1;
+
+ c_had += 1;
+ }
+// System.out.println("c_had and NHit_had: "+c_had+" "+NHit_had);
+
+
+ for (int i=0; i < NHit_had; i++) // fill SP information
+ {
+
+ MCParticle loop_mc;
+ for (int j=0; j < ShowerP_had[i]; j++) // loop over contributing particles
+ {
+ loop_mc = HitList_had[i].getMCParticle(j);
+ Hep3Vector p_o = loop_mc.getOrigin();
+ Hep3Vector p_e = loop_mc.getEndPoint();
+ double r_o = Math.sqrt(p_o.x()*p_o.x() + p_o.y()*p_o.y());
+ double r_e = Math.sqrt(p_e.x()*p_e.x() + p_e.y()*p_e.y());
+
+ int Found = 0;
+ int index = 0;
+ for (MCParticle mc : PList)
+ {
+ if ((loop_mc.getPDGID()==mc.getPDGID())&&((Math.abs(loop_mc.getEnergy()-mc.getEnergy())/mc.getEnergy())<0.00001))
+ {
+ if ((Math.abs(loop_mc.getEndPoint().x()-mc.getEndPoint().x())<0.00001)&&(Math.abs(loop_mc.getEndPoint().y()-mc.getEndPoint().y())<0.00001))
+ {
+ Found += 1;
+ index = PList.indexOf(mc);
+ }
+ }
+ }
+
+ if (Found == 1)
+ {
+ SPID_had[i][j] = index;
+ }
+ else if (Found == 0)
+ {
+ if (loop_mc.getGeneratorStatus() == 1)
+ {
+// aida.cloud1D("origin r: FS particle").fill(r_o);
+// aida.cloud1D("end r: FS particle").fill(r_e);
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+
+ int pFound = 0;
+ MCParticle tt = loop_mc;
+ while ((pFound == 0)&&(tt.getOrigin().magnitude() > 20.)) // stop searching while found one or parent is within 20.mm from IP
+ {
+ List<MCParticle> parents = tt.getParents();
+ if (parents.size() > 1) System.out.println("Attn: more than one parents found before finding SP: N = "+parents.size());
+ for (MCParticle p_mc : parents)
+ {
+ tt = p_mc;
+ }
+
+ for (MCParticle mc : PList)
+ {
+ if ((tt.getPDGID()==mc.getPDGID())&&((Math.abs(tt.getEnergy()-mc.getEnergy())/mc.getEnergy())<0.00001))
+ {
+ if ((Math.abs(tt.getEndPoint().x()-mc.getEndPoint().x())<0.00001)&&(Math.abs(tt.getEndPoint().y()-mc.getEndPoint().y())<0.00001))
+ {
+ pFound += 1;
+ SPID_had[i][j] = PList.indexOf(mc);
+ }
+ }
+ }
+ } // searching parents loop
+
+ if (pFound == 0)
+ {
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+ PList.add(loop_mc);
+ SPID_had[i][j] = PList.size() - 1; // last particle in list
+ }
+ }
+ else if (loop_mc.getGeneratorStatus() == 0)
+ {
+// aida.cloud1D("origin r: other").fill(r_o);
+// aida.cloud1D("end r: other").fill(r_e);
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+
+ int pFound = 0;
+ MCParticle tt = loop_mc;
+ while ((pFound == 0)&&(tt.getOrigin().magnitude() > 20.)) // stop searching while found one or parent is within 20.mm from IP
+ {
+ List<MCParticle> parents = tt.getParents();
+ if (parents.size() > 1) System.out.println("Attn: more than one parents found before finding SP: N = "+parents.size());
+ for (MCParticle p_mc : parents)
+ {
+ tt = p_mc;
+ }
+
+ for (MCParticle mc : PList)
+ {
+ if ((tt.getPDGID()==mc.getPDGID())&&((Math.abs(tt.getEnergy()-mc.getEnergy())/mc.getEnergy())<0.00001))
+ {
+ if ((Math.abs(tt.getEndPoint().x()-mc.getEndPoint().x())<0.00001)&&(Math.abs(tt.getEndPoint().y()-mc.getEndPoint().y())<0.00001))
+ {
+ pFound += 1;
+ SPID_had[i][j] = PList.indexOf(mc);
+ }
+ }
+ }
+ } // searching parents loop
+
+ if (pFound == 0)
+ {
+// System.out.println(" particle "+loop_mc.getType().getName()+" "+loop_mc.getEnergy()+" origin, end = "+r_o+" "+p_o.z()+" "+r_e+" "+p_e.z());
+ if (loop_mc.getPDGID() != 0)
+ {
+ PList.add(loop_mc);
+ SPID_em[i][j] = PList.size() - 1; // last particle in list
+ }
+ else
+ {
+ while (loop_mc.getPDGID() == 0)
+ {
+ loop_mc = loop_mc.getParents().get(0);
+ }
+ PList.add(loop_mc);
+ SPID_em[i][j] = PList.size() - 1; // last particle in list
+ }
+ }
+ }
+ else // so far, no particles found with a type other than 'FS' or 'other'
+ {
+ aida.cloud1D("origin r: rest").fill(r_o);
+ aida.cloud1D("end r: rest").fill(r_e);
+ }
+ }
+ else
+ {
+ System.out.println(" Something went wrong: Found = "+Found+" index = "+index);
+ }
+
+ } // end of contributing particles loop
+
+ } // end of fill SP loop
+
+
+//--------------------------------------------------------------------------------------------
+// calculate perfect PFA here: althought PList have double counting in energy, we still try to
+// get an idea how things work
+
+ double E_AD = 0.; // total visible energy from P_PFA
+ double E_p = 0.; // total particle momentum, for all mc particles in the list, to see how bad is the double counting
+
+ double E_photon = 0.; // photon energy
+ double E_neuhad = 0.; // neutral hadron energy
+ double E_endcap = 0.; // endcap neutral energy
+
+ for (MCParticle mc : PList)
+ {
+ double E_em = 0.; // particle partial energy in ecal
+ double N_had = 0.; // particle partial energy (number of hits) in hcal
+
+ E_p += mc.getMomentum().magnitude();
+
+// System.out.println(" loop over particle #"+PList.indexOf(mc)+" "+mc.getType().getName());
+ if (Math.abs(mc.getType().getCharge()) > 0.01)
+ { E_AD += mc.getMomentum().magnitude();}
+ else
+ {
+ for (int j = 0; j < NHit_em; j++)
+ {
+ int f_contributer = 0;
+ for (int k = 0; k < ShowerP_em[j]; k++)
+ { if (SPID_em[j][k] == PList.indexOf(mc)) f_contributer = 1;}
+
+ if (f_contributer == 1) E_em += HitList_em[j].getRawEnergy()/ShowerP_em[j];
+ }
+
+ for (int j = 0; j < NHit_had; j++)
+ {
+ int f_contributer = 0;
+ for (int k = 0; k < ShowerP_had[j]; k++)
+ { if (SPID_had[j][k] == PList.indexOf(mc)) f_contributer = 1;}
+
+ if (f_contributer == 1) N_had += 1./ShowerP_had[j];
+ }
+
+// here is another place that calibration comes in! can choose one of the (two)/three now:
+
+// 1. original, done with all hits of single particle
+// if ((Math.abs(PList.get(CFID[i][0]).getPDGID()) == 11)||(Math.abs(PList.get(CFID[i][0]).getPDGID()) == 22))
+// {EFinal[i] = 82.5*EF_em[i] + 0.11*NF_had[i];}
+// else
+// {EFinal[i] = (EF_em[i]*114.19 + 0.12192*NF_had[i])/(1 + EF_em[i]*1.3359 + NF_had[i]*0.0012218);}
+//
+// 2. after clustering, done with 1 major cluster only
+// if ((Math.abs(PList.get(CFID[i][0]).getPDGID()) == 11)||(Math.abs(PList.get(CFID[i][0]).getPDGID()) == 22))
+// {EFinal[i] = (87.883*EF_em[i] + 0.12*NF_had[i])/(1 + 0.0602*EF_em[i]);}
+// else
+// {EFinal[i] = (EF_em[i]*153.29 + 0.12076*NF_had[i])/(1 + EF_em[i]*2.7321 + NF_had[i]*0.0003390);}
+//
+// 3. new calibration for sidaug05_np non-projective detector!
+// if ((Math.abs(PList.get(CFID[i][0]).getPDGID()) == 11)||(Math.abs(PList.get(CFID[i][0]).getPDGID()) == 22))
+// {EFinal[i] = 82.5*EF_em[i] + 0.11*NF_had[i];}
+// else
+// {EFinal[i] = Math.pow((EF_em[i]*93.54 + 0.12947*NF_had[i]),0.95323);}
+
+ if (mc.getPDGID() == 22) // photons, from RPC_calibration
+ {
+ E_AD += E_em*82.5 + 0.11*N_had;
+ E_photon += mc.getEnergy();
+ if ((mc.getPZ()*mc.getPZ())/(mc.getPX()*mc.getPX()+mc.getPY()*mc.getPY()+mc.getPZ()*mc.getPZ())>0.8) E_endcap += mc.getEnergy();
+ }
+ else // neutral hadrons, from RPC_calibration
+ {
+ E_AD += Math.pow((E_em*93.54 + 0.12947*N_had),0.95323);
+ E_neuhad += mc.getEnergy();
+ if ((mc.getPZ()*mc.getPZ())/(mc.getPX()*mc.getPX()+mc.getPY()*mc.getPY()+mc.getPZ()*mc.getPZ())>0.8) E_endcap += mc.getEnergy();
+// System.out.println("real and measured energy: "+MyEList[i]+" "+(E_em*114.19 + 0.12192*N_had)/(1 + E_em*1.3359 + N_had*0.0012218));
+ }
+ }
+ }
+
+ aida.cloud1D("Perfect PFA (all events)").fill(E_AD);
+ aida.cloud1D("All Particle P's (all events)").fill(E_p);
+ if (E_nu_fs < 0.2)
+ {
+ aida.cloud1D("Perfect PFA (E_nu_fs < 200MeV)").fill(E_AD);
+
+ if (E_photon + E_neuhad < 0.5) aida.cloud1D("Perfect PFA (barrel,0.2)").fill(E_AD);
+ else if (E_endcap/(E_photon + E_neuhad) < 0.05) aida.cloud1D("Perfect PFA (barrel,0.2)").fill(E_AD);
+ else aida.cloud1D("Perfect PFA (endcap,0.2)").fill(E_AD);
+ }
+
+ aida.cloud1D("photon energy").fill(E_photon);
+ aida.cloud1D("neutral hadron energy").fill(E_neuhad);
+ aida.cloud1D("all neutral energy").fill(E_photon + E_neuhad);
+
+
+//------------------------------
+//calculate density for all hits
+
+ double[] density_em = new double[NHit_em];
+ double[] density_had = new double[NHit_had];
+
+// double[] d_em = new double[NHit_em];
+// double[] d_had = new double[NHit_had];
+
+ for (int i=0; i < NHit_em; i++) // calculate density for em hits
+ {
+ density_em[i] = 0;
+// d_em[i] = 0;
+ double[] x0 = HitList_em[i].getPosition();
+
+ double d1_em = 4.0;
+ double d2_em = 4.0;
+ double d3_em = 3.75;
+
+ for (int j=0; j < NHit_em; j++) // loop over other em hits
+ {
+ if (j != i)
+ {
+ double[] x1 = HitList_em[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_em[i][0]+(x1[1]-x0[1])*v1_em[i][1]+(x1[2]-x0[2])*v1_em[i][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_em[i][0]+(x1[1]-x0[1])*v2_em[i][1]+(x1[2]-x0[2])*v2_em[i][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_em[i][0]+(x1[1]-x0[1])*v3_em[i][1]+(x1[2]-x0[2])*v3_em[i][2]);
+
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 3))
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 0.5)) {d_em[i] += 1.;}
+ if (((DotProduct1/d1_em) < 15)&&((DotProduct2/d2_em) < 15)&&((DotProduct3/d3_em) < 15))
+ {
+ density_em[i] += Math.exp(-DotProduct1/d1_em -DotProduct2/d2_em -DotProduct3/d3_em);
+ }
+ }
+ }
+
+ for (int j=0; j < NHit_had; j++) // loop over had hits
+ {
+ double[] x1 = HitList_had[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_em[i][0]+(x1[1]-x0[1])*v1_em[i][1]+(x1[2]-x0[2])*v1_em[i][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_em[i][0]+(x1[1]-x0[1])*v2_em[i][1]+(x1[2]-x0[2])*v2_em[i][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_em[i][0]+(x1[1]-x0[1])*v3_em[i][1]+(x1[2]-x0[2])*v3_em[i][2]);
+
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 3))
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 0.5)) {d_em[i] += 1.;}
+ if (((DotProduct1/d1_em) < 15)&&((DotProduct2/d2_em) < 15)&&((DotProduct3/d3_em) < 15))
+ {
+ density_em[i] += Math.exp(-DotProduct1/d1_em -DotProduct2/d2_em -DotProduct3/d3_em);
+ }
+ }
+
+// System.out.println(i+" x,y,z and density = "+x0[0]+" "+x0[1]+" "+x0[2]+" "+density_em[i]);
+ }
+
+ for (int i=0; i < NHit_had; i++) // calculate density for had hits
+ {
+ density_had[i] = 0;
+// d_had[i] = 0;
+ double[] x0 = HitList_had[i].getPosition();
+
+ double d1_had = 10.;
+ double d2_had = 10.;
+ double d3_had = 28.0;
+
+ for (int j=0; j < NHit_em; j++) // loop over em hits
+ {
+ double[] x1 = HitList_em[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_had[i][0]+(x1[1]-x0[1])*v1_had[i][1]+(x1[2]-x0[2])*v1_had[i][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_had[i][0]+(x1[1]-x0[1])*v2_had[i][1]+(x1[2]-x0[2])*v2_had[i][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_had[i][0]+(x1[1]-x0[1])*v3_had[i][1]+(x1[2]-x0[2])*v3_had[i][2]);
+
+// if (((DotProduct1/d1_had) < 3)&&((DotProduct2/d2_had) < 3)&&((DotProduct3/d3_had) < 3))
+// if (((DotProduct1/d1_had) < 3)&&((DotProduct2/d2_had) < 3)&&((DotProduct3/d3_had) < 0.1)) {d_had[i] += 1.;}
+ if (((DotProduct1/d1_had) < 15)&&((DotProduct2/d2_had) < 15)&&((DotProduct3/d3_had) < 15))
+ {
+ density_had[i] += Math.exp(-DotProduct1/d1_had -DotProduct2/d2_had -DotProduct3/d3_had);
+ }
+ }
+
+ for (int j=0; j < NHit_had; j++) // loop over other had hits
+ {
+ if (j != i)
+ {
+ double[] x1 = HitList_had[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_had[i][0]+(x1[1]-x0[1])*v1_had[i][1]+(x1[2]-x0[2])*v1_had[i][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_had[i][0]+(x1[1]-x0[1])*v2_had[i][1]+(x1[2]-x0[2])*v2_had[i][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_had[i][0]+(x1[1]-x0[1])*v3_had[i][1]+(x1[2]-x0[2])*v3_had[i][2]);
+
+// if (((DotProduct1/d1_had) < 3)&&((DotProduct2/d2_had) < 3)&&((DotProduct3/d3_had) < 3))
+// if (((DotProduct1/d1_had) < 3)&&((DotProduct2/d2_had) < 3)&&((DotProduct3/d3_had) < 0.5)) {d_had[i] += 1.;}
+ if (((DotProduct1/d1_had) < 15)&&((DotProduct2/d2_had) < 15)&&((DotProduct3/d3_had) < 15))
+ {
+ density_had[i] += Math.exp(-DotProduct1/d1_had -DotProduct2/d2_had -DotProduct3/d3_had);
+ }
+ }
+ }
+
+// System.out.println(i+" x,y,z and density = "+x0[0]+" "+x0[1]+" "+x0[2]+" "+density_had[i]);
+ }
+
+
+
+//-------------------------------------------------------
+// try to find all clusters, by looping through seed hits
+
+ int NClusters = 0; // total number of clusters
+ int NFinal = 0; // total number of clusters that survived quality cut
+
+ int NewClusters = 1; // number of new clusters picked up in this loop
+ while (NewClusters > 0)
+ {
+ NewClusters = 0;
+
+ int Seed = -1;
+ double SeedDensity = -1.;
+ int SeedID = 1; // SeedID = 1 for ECAL seed, SeedID = -1 for HCAL seed
+ for (int i = 0; i < NHit_em; i++) // try to find a seed among EM hits first
+ {
+ if ((EM_ClusterID[i] == -1)&&(density_em[i] > SeedDensity))
+ {
+ Seed = i;
+ SeedDensity = density_em[i];
+ }
+ }
+
+ if (SeedDensity >= 0) // found one seed in em hits (>0 will exclude isolated hits, >=0 will include them)
+ {
+ NewClusters = 1;
+ NClusters += 1;
+ NFinal += 1;
+
+ EM_ClusterID[Seed] = NClusters - 1;
+ EM_FinalID[Seed] = NFinal - 1;
+ SeedID = 1;
+ }
+ else // used up em hits, now turn to had hits
+ {
+ Seed = -1;
+ for (int i = 0; i < NHit_had; i++)
+ {
+ if ((HAD_ClusterID[i] == -1)&&(density_had[i] > SeedDensity))
+ {
+ Seed = i;
+ SeedDensity = density_had[i];
+ }
+ }
+
+ if (SeedDensity >= 0)
+ {
+ NewClusters = 1;
+ NClusters += 1;
+ NFinal += 1;
+
+ HAD_ClusterID[Seed] = NClusters - 1;
+ HAD_FinalID[Seed] = NFinal - 1;
+ SeedID = -1;
+ }
+ } // here, one should find one seed, unless used up all hits in em and had
+
+ if (NewClusters > 0) // if there's a seed, pick up nearby hits first, with a loose cut
+ {
+ double LocalDensity = 0.;
+ int NumberofHits = 0; // number of hits picked up in a round
+
+ if (SeedID == 1) // em seed
+ {
+ double[] x0 = HitList_em[Seed].getPosition();
+ double d1_em = 4.0;
+ double d2_em = 4.0;
+ double d3_em = 3.75;
+
+ for (int j=0; j < NHit_em; j++) // loop over em hits
+ {
+ if (EM_ClusterID[j] == -1)
+ {
+ LocalDensity = 0.;
+ double[] x1 = HitList_em[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_em[Seed][0]+(x1[1]-x0[1])*v1_em[Seed][1]+(x1[2]-x0[2])*v1_em[Seed][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_em[Seed][0]+(x1[1]-x0[1])*v2_em[Seed][1]+(x1[2]-x0[2])*v2_em[Seed][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_em[Seed][0]+(x1[1]-x0[1])*v3_em[Seed][1]+(x1[2]-x0[2])*v3_em[Seed][2]);
+
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 0.1)) {LocalDensity = 1.;}
+ if (((DotProduct1/d1_em) < 15)&&((DotProduct2/d2_em) < 15)&&((DotProduct3/d3_em) < 15))
+ {
+ LocalDensity = Math.exp(-DotProduct1/d1_em -DotProduct2/d2_em -DotProduct3/d3_em);
+ }
+
+ if (LocalDensity > 0.010) // add nearby hits to the cluster
+ {
+ EM_ClusterID[j] = NClusters - 1;
+ EM_FinalID[j] = NFinal - 1;
+ NumberofHits += 1;
+
+// System.out.println(" em seed: "+Seed+" pickup em hit: "+j+" LocalDensity = "+LocalDensity);
+ }
+ }
+ }
+
+ for (int j=0; j < NHit_had; j++) // loop over had hits
+ {
+ if (HAD_ClusterID[j] == -1)
+ {
+ LocalDensity = 0.;
+ double[] x1 = HitList_had[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_em[Seed][0]+(x1[1]-x0[1])*v1_em[Seed][1]+(x1[2]-x0[2])*v1_em[Seed][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_em[Seed][0]+(x1[1]-x0[1])*v2_em[Seed][1]+(x1[2]-x0[2])*v2_em[Seed][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_em[Seed][0]+(x1[1]-x0[1])*v3_em[Seed][1]+(x1[2]-x0[2])*v3_em[Seed][2]);
+
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 0.5)) {LocalDensity = 1.;}
+ if (((DotProduct1/d1_em) < 15)&&((DotProduct2/d2_em) < 15)&&((DotProduct3/d3_em) < 15))
+ {
+ LocalDensity = Math.exp(-DotProduct1/d1_em -DotProduct2/d2_em -DotProduct3/d3_em);
+ }
+
+ if (LocalDensity > 0.010) // add nearby hits to the cluster
+ {
+ HAD_ClusterID[j] = NClusters - 1;
+ HAD_FinalID[j] = NFinal - 1;
+ NumberofHits += 1;
+
+// System.out.println(" em seed: "+Seed+" pickup had hit: "+j+" LocalDensity = "+LocalDensity);
+ }
+ }
+ }
+ } // end bracket for em seed
+
+ if (SeedID == -1) // had seed
+ {
+ double[] x0 = HitList_had[Seed].getPosition();
+ double d1_had = 10.;
+ double d2_had = 10.;
+ double d3_had = 28.0;
+
+ for (int j=0; j < NHit_em; j++) // loop over em hits
+ {
+ if (EM_ClusterID[j] == -1)
+ {
+ LocalDensity = 0.;
+ double[] x1 = HitList_em[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_had[Seed][0]+(x1[1]-x0[1])*v1_had[Seed][1]+(x1[2]-x0[2])*v1_had[Seed][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_had[Seed][0]+(x1[1]-x0[1])*v2_had[Seed][1]+(x1[2]-x0[2])*v2_had[Seed][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_had[Seed][0]+(x1[1]-x0[1])*v3_had[Seed][1]+(x1[2]-x0[2])*v3_had[Seed][2]);
+
+// if (((DotProduct1/d1_had) < 3)&&((DotProduct2/d2_had) < 3)&&((DotProduct3/d3_had) < 0.1)) {LocalDensity = 1.;}
+ if (((DotProduct1/d1_had) < 15)&&((DotProduct2/d2_had) < 15)&&((DotProduct3/d3_had) < 15))
+ {
+ LocalDensity = Math.exp(-DotProduct1/d1_had -DotProduct2/d2_had -DotProduct3/d3_had);
+ }
+
+ if (LocalDensity > 0.010) // add nearby hits to the cluster
+ {
+ EM_ClusterID[j] = NClusters - 1;
+ EM_FinalID[j] = NFinal - 1;
+ NumberofHits += 1;
+
+// System.out.println(" had seed: "+Seed+" pickup em hit: "+j+" LocalDensity = "+LocalDensity);
+ }
+ }
+ }
+
+ for (int j=0; j < NHit_had; j++) // loop over had hits
+ {
+ if (HAD_ClusterID[j] == -1)
+ {
+ LocalDensity = 0.;
+ double[] x1 = HitList_had[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_had[Seed][0]+(x1[1]-x0[1])*v1_had[Seed][1]+(x1[2]-x0[2])*v1_had[Seed][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_had[Seed][0]+(x1[1]-x0[1])*v2_had[Seed][1]+(x1[2]-x0[2])*v2_had[Seed][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_had[Seed][0]+(x1[1]-x0[1])*v3_had[Seed][1]+(x1[2]-x0[2])*v3_had[Seed][2]);
+
+// if (((DotProduct1/d1_had) < 3)&&((DotProduct2/d2_had) < 3)&&((DotProduct3/d3_had) < 0.5)) {LocalDensity = 1.;}
+ if (((DotProduct1/d1_had) < 15)&&((DotProduct2/d2_had) < 15)&&((DotProduct3/d3_had) < 15))
+ {
+ LocalDensity = Math.exp(-DotProduct1/d1_had -DotProduct2/d2_had -DotProduct3/d3_had);
+ }
+
+ if (LocalDensity > 0.010) // add nearby hits to the cluster
+ {
+ HAD_ClusterID[j] = NClusters - 1;
+ HAD_FinalID[j] = NFinal - 1;
+ NumberofHits += 1;
+
+// System.out.println(" had seed: "+Seed+" pickup had hit: "+j+" LocalDensity = "+LocalDensity);
+ }
+ }
+ }
+ } // end bracket for had seed
+
+ } // end bracket for picking up nearby hits
+
+
+ if (NewClusters > 0) // if there's a seed at all, continue to find the cluster, otherwise, go to the end of loop
+ {
+ int NumberofHits = 1;
+ while (NumberofHits > 0) // try to pick up more hits from a seed cluster
+ {
+ NumberofHits = 0;
+
+ for (int i=0; i < NHit_em; i++) // see if em hits can be included in cluster
+ {
+ double Product1=0., Product2=0., Product3=0.;
+ if (EM_ClusterID[i] == -1) // if not assigned to any cluster
+ {
+ double LocalDensity = 0;
+ double[] x0 = HitList_em[i].getPosition();
+ double d1_em = 4.0;
+ double d2_em = 4.0;
+ double d3_em = 3.75;
+
+ for (int j = 0; j < NHit_em; j++) // loop over existing em hits of the cluster
+ {
+ if (EM_ClusterID[j] == (NClusters - 1))
+ {
+ double[] x1 = HitList_em[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_em[i][0]+(x1[1]-x0[1])*v1_em[i][1]+(x1[2]-x0[2])*v1_em[i][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_em[i][0]+(x1[1]-x0[1])*v2_em[i][1]+(x1[2]-x0[2])*v2_em[i][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_em[i][0]+(x1[1]-x0[1])*v3_em[i][1]+(x1[2]-x0[2])*v3_em[i][2]);
+
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 0.5)) {LocalDensity += 1.;}
+ if (((DotProduct1/d1_em) < 15)&&((DotProduct2/d2_em) < 15)&&((DotProduct3/d3_em) < 15))
+ {
+ LocalDensity += Math.exp(-DotProduct1/d1_em -DotProduct2/d2_em -DotProduct3/d3_em);
+ }
+ }
+ }
+
+ for (int j = 0; j < NHit_had; j++) // loop over existing had hits of the cluster
+ {
+ if (HAD_ClusterID[j] == (NClusters - 1))
+ {
+ double[] x1 = HitList_had[j].getPosition();
+ double DotProduct1 = Math.abs((x1[0]-x0[0])*v1_em[i][0]+(x1[1]-x0[1])*v1_em[i][1]+(x1[2]-x0[2])*v1_em[i][2]);
+ double DotProduct2 = Math.abs((x1[0]-x0[0])*v2_em[i][0]+(x1[1]-x0[1])*v2_em[i][1]+(x1[2]-x0[2])*v2_em[i][2]);
+ double DotProduct3 = Math.abs((x1[0]-x0[0])*v3_em[i][0]+(x1[1]-x0[1])*v3_em[i][1]+(x1[2]-x0[2])*v3_em[i][2]);
+
+// if (((DotProduct1/d1_em) < 3)&&((DotProduct2/d2_em) < 3)&&((DotProduct3/d3_em) < 0.5)) {LocalDensity += 1.;}
+ if (((DotProduct1/d1_em) < 15)&&((DotProduct2/d2_em) < 15)&&((DotProduct3/d3_em) < 15))
+ {
+ LocalDensity += Math.exp(-DotProduct1/d1_em -DotProduct2/d2_em -DotProduct3/d3_em);
+ }
+ }
+ }
+
+ if (LocalDensity > 0.01)
+ {
+ EM_ClusterID[i] = NClusters - 1;
+ EM_FinalID[i] = NFinal - 1;
+ NumberofHits += 1;
+
+// System.out.println(" pickup em hit: "+i+" LocalDensity = "+LocalDensity);
+ }
+
+ }
+ } // end bracket for looking at new em hits
+
+
+ for (int i=0; i < NHit_had; i++) // see if had hits can be included in cluster
+ {
+ double Product1=0., Product2=0., Product3=0.;
+ if (HAD_ClusterID[i] == -1) // if not assigned to any cluster
+ {
+ double LocalDensity = 0;
+ double[] x0 = HitList_had[i].getPosition();
+ double d1_had = 10.;
[truncated at 1000 lines; 1492 more skipped]
