lcsim/src/org/lcsim/contrib/SteveMagill
diff -N PhotonFinderDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PhotonFinderDriver.java 23 Apr 2007 20:04:31 -0000 1.1
@@ -0,0 +1,530 @@
+package org.lcsim.contrib.SteveMagill;
+
+// This driver extrapolates tracks into CAL, associating mips and determining interaction layer
+// output : Mip Clusters and IL per track
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.util.Driver;
+import org.lcsim.util.swim.*;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.ITree;
+import org.lcsim.geometry.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.recon.ztracking.cheater.*;
+import org.lcsim.recon.cluster.fixedcone.*;
+import org.lcsim.recon.cluster.analysis.*;
+import org.lcsim.recon.emid.hmatrix.HMatrix;
+import org.lcsim.recon.emid.hmatrix.HMatrixTask;
+import org.lcsim.recon.emid.hmatrix.HMatrixBuilder;
+import org.lcsim.recon.emid.hmatrix.HMatrixConditionsConverter;
+import org.lcsim.math.chisq.ChisqProb;
+import org.lcsim.recon.cluster.nn.NearestNeighborCluster;
+import org.lcsim.spacegeom.*;
+import org.lcsim.util.hitmap.HitMap;
+
+public class PhotonFinderDriver extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ private int _mincells;
+ private double _dTrcl;
+ private HMatrixTask _task = HMatrixTask.ANALYZE;
+ private int _nLayersB;
+ private int _nLayersEC;
+ private HMatrixBuilder _hmb;
+ private HMatrix _hmx;
+ private ITree _tree;
+ private boolean _initialized = false;
+ private int _nmeasB;
+ private int _nmeasEC;
+ private static final double _log10inv = 1./Math.log(10.0);
+ private int _logEIndexB;
+ private int _logEIndexEC;
+ double[] _valsB;
+ double[] _valsEC;
+ private double[] BRadii = new double[100];
+ private double[] ECZs = new double[100];
+ private boolean phoD = false;
+ private boolean phoRes = false;
+
+ public PhotonFinderDriver(int mincells, double dTrcl)
+ {
+ _mincells = mincells;
+ _dTrcl = dTrcl;
+ // The HMatrix could possibly change, so be sensitive to this.
+ getConditionsManager().registerConditionsConverter(new HMatrixConditionsConverter());
+ }
+
+// public PhotonFinderDriver()
+// {
+// this(HMatrixTask.ANALYZE);
+// }
+
+// public PhotonFinderDriver(HMatrixTask task)
+// {
+// _tree = aida.tree();
+// _task = task;
+// if(_task==HMatrixTask.ANALYZE)
+// {
+// // The HMatrix could possibly change, so be sensitive to this.
+// getConditionsManager().registerConditionsConverter(new HMatrixConditionsConverter());
+// }
+// }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event); // executes all added drivers
+
+ boolean mcfltr = true;
+// mcfltr = (Boolean)event.get("MCFilt");
+// System.out.println("MCFilt is " +mcfltr);
+ if (mcfltr)
+ {
+ if(!_initialized)
+ {
+ CylindricalCalorimeter calsubB = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMBarrel");
+ _nLayersB = calsubB.getLayering().getLayerCount();
+// System.out.println("found "+_nLayersB+" layers in the EMBarrel");
+ _nmeasB = _nLayersB;
+ _logEIndexB = _nmeasB;
+ _nmeasB+=1;
+ _valsB = new double[_nmeasB];
+ for (int i=0; i<_nLayersB; ++i)
+ {
+ BRadii[i]=calsubB.getLayering().getDistanceToLayerSensorMid(i);
+ }
+
+ CylindricalCalorimeter calsubEC = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMEndcap");
+ _nLayersEC = calsubEC.getLayering().getLayerCount();
+// System.out.println("found "+_nLayersEC+" layers in the EMEndcap");
+ _nmeasEC = _nLayersEC;
+ _logEIndexEC = _nmeasEC;
+ _nmeasEC+=1;
+ _valsEC = new double[_nmeasEC];
+ for (int i=0; i<_nLayersEC; ++i)
+ {
+ ECZs[i]=calsubEC.getLayering().getDistanceToLayerSensorMid(i);
+ }
+
+ //FIXME key needs to be better defined
+ int key = 0;
+ if(_task==HMatrixTask.ANALYZE)
+ {
+ //FIXME need to fetch name of HMAtrix file to use from a conditions file
+// System.out.println("Looking for HMatrix...");
+ _hmx = getConditionsManager().getCachedConditions(HMatrix.class,
+ "LongitudinalHMatrix.hmx").getCachedData();
+ }
+ _initialized = true;
+ }
+ // First, need to make list of track positions at EM shower max for each track
+ List<PerfectTrack> evtracks = event.get(PerfectTrack.class,"PerfectTracks");
+ double[] trth = new double[200];
+ double[] trph = new double[200];
+ double TotTrP = 0;
+ int ntr=0;
+ for (Track ictrk : evtracks)
+ {
+ Detector det = event.getDetector();
+ double[] zero = {0, 0, 0};
+ double[] _fieldStrength = det.getFieldMap().getField(zero);
+ HelixSwimmer tswim = new HelixSwimmer(_fieldStrength[2]);
+ double[] trp = ictrk.getMomentum();
+ Hep3Vector trp3 = new BasicHep3Vector(trp);
+ double[] tror = ictrk.getReferencePoint();
+ Hep3Vector tror3 = new BasicHep3Vector(tror);
+ int trq = ictrk.getCharge();
+ double cttrp = Math.sqrt(trp[0]*trp[0]+trp[1]*trp[1]+trp[2]*trp[2]);
+ if (phoD) aida.cloud1D("Track momentum in PhoFinder").fill(cttrp);
+ TotTrP += cttrp;
+ tswim.setTrack(trp3, tror3, trq);
+ double torad = tswim.getDistanceToRadius(BRadii[3]); // this is pre-shower max (sm=7) in EM Barrel Cal
+ double toz = tswim.getDistanceToZ(ECZs[3]); // this is pre-shower max (sm=7) in EM Endcap Cal
+ if (torad < toz) // in Barrel
+ {
+ SpacePoint trackSP = tswim.getPointAtDistance(torad);
+// System.out.println("SpacePoint at EMB Shower Max " + trackSP);
+ trth[ntr] = trackSP.theta();
+ if (phoD) aida.cloud1D("Track Theta at EMB Shower Max").fill(trth[ntr]);
+ trph[ntr] = trackSP.phi();
+ if (phoD) aida.cloud1D("Track Phi at EMB Shower Max").fill(trph[ntr]);
+ if (phoD) aida.cloud2D("Theta vs Phi for tracks at EM Shower Max").fill(trph[ntr],trth[ntr]);
+ } else // in Endcap
+ {
+ SpacePoint trackSP = tswim.getPointAtDistance(toz);
+// System.out.println("SpacePoint at EMEC Shower Max " + trackSP);
+ trth[ntr] = trackSP.theta();
+ if (phoD) aida.cloud1D("Track Theta at EMEC Shower Max").fill(trth[ntr]);
+ trph[ntr] = trackSP.phi();
+ if (phoD) aida.cloud1D("Track Phi at EMEC Shower Max").fill(trph[ntr]);
+ if (phoD) aida.cloud2D("Theta vs Phi for tracks at EM Shower Max").fill(trph[ntr],trth[ntr]);
+ }
+ ntr++;
+ } // done getting theta and phi of tracks
+
+ // loop over all EM clusters, check starting layer, check to see if no tracks nearby,
+ // make a new NN1118 cluster out of the hits,
+ // evaluate 1118 cluster with the HMatrix
+ List<BasicCluster> CandClusters = new ArrayList<BasicCluster>();
+ List<BasicCluster> TotNNclusters = new ArrayList<BasicCluster>();
+ List<BasicCluster> phoBclusters = new ArrayList<BasicCluster>();
+ List<BasicCluster> phoECclusters = new ArrayList<BasicCluster>();
+ try
+ {
+ List<BasicCluster> EBclusters = event.get(BasicCluster.class,"EBTMHitsFCClus");
+
+ for (BasicCluster ebclus : EBclusters)
+ {
+ List<BasicCluster> ennclusters = new ArrayList<BasicCluster>();
+ BasicCluster candclus = new BasicCluster();
+ // check to see if track is close to this cluster before HMatrix evaluation
+// System.out.println("Barrel FC Cluster");
+ double ecp[] = ebclus.getPosition();
+ double ecpx = ecp[0];
+ double ecpy = ecp[1];
+ double ecpz = ecp[2];
+ double ecR = Math.sqrt(ecpx*ecpx+ecpy*ecpy);
+ double ebclth = Math.atan(ecR/ecpz);
+ if (ebclth<0) ebclth+=Math.PI;
+ double ebclph = Math.atan2(ecpy,ecpx);
+ if (ebclph<0) ebclph+=2*Math.PI;
+ if (phoD) aida.cloud2D("Theta vs Phi of EB Clusters").fill(ebclth,ebclph);
+ // check for tracks near clusters - some of these are interacting hadrons
+// System.out.println("Number of tracks to check with cands " + ntr);
+ int ntrclmatch = 0;
+ for (int i=0; i<ntr; ++i)
+ {
+ double trcldelth = Math.abs(ebclth-trth[i]);
+ if (phoD) aida.cloud1D("Clus-Trk theta diff in PhoFinder").fill(trcldelth);
+ double trcldelph = Math.abs(ebclph-trph[i]);
+ if (trcldelph > Math.PI) trcldelph = 2*Math.PI-trcldelph;
+ if (phoD) aida.cloud1D("Clus-Trk phi diff in PhoFinder").fill(trcldelph);
+ double trcldist = Math.sqrt(trcldelth*trcldelth+trcldelph*trcldelph);
+ if (phoD) aida.cloud1D("Trk-Clus Distance in PhoFinder").fill(trcldist);
+ // check for distance from a track - used 0.03, but .025 looks better
+ if (trcldist<_dTrcl) ntrclmatch++;
+ }
+ // if any matches, don't test this cluster - its not a photon
+ if (ntrclmatch>0) continue;
+ // Make a hitmap of cluser hits to re-cluster in NN clusterer
+ HitMap ebhitmap = new HitMap();
+ for (CalorimeterHit ebhit : ebclus.getCalorimeterHits())
+ {
+ long cID = ebhit.getCellID();
+ ebhitmap.put(cID, ebhit);
+ }
+ int _dU = 1;
+ int _dV = 1;
+ int _dLayer = 1;
+ double ethresh = 0.;
+
+ for (;;)
+ {
+ if (ebhitmap.isEmpty()) break;
+ long hID = ebhitmap.keySet().iterator().next();
+ CalorimeterHit hit = ebhitmap.get(hID);
+ NearestNeighborCluster nnclus = new NearestNeighborCluster(ebhitmap, hit, hID, _dU, _dV, _dLayer, ethresh);
+ if(nnclus.getCalorimeterHits().size()>_mincells)
+ {
+ ennclusters.add(nnclus);
+// System.out.println("Formed NN cluster from FC cluster");
+ TotNNclusters.add(nnclus);
+ }
+ }
+ // if number of NN clusters is = 1, use it to test with H-Matrix, else use original FC cluster
+// System.out.println("Number of NN clusters " + ennclusters.size());
+ if (ennclusters.size() == 0)
+ {
+ candclus = ebclus;
+ } else if (ennclusters.size() == 1)
+ {
+ candclus = ennclusters.get(0);
+ } else if (ennclusters.size()>1)
+ {
+ candclus = ebclus;
+ }
+ CandClusters.add(candclus);
+ // check if this cluster passes H-Matrix test, if so, add candclus to list of Photon Clusters
+ double energy = candclus.getRawEnergy();
+ double[] layerE = layerEnergiesB(candclus);
+ // accumulate the longitudinal energy fractions...
+ // tried to kludge to correct for photons interacting late
+ int layfirst = 0;
+ for(int j=0; j<layerE.length; ++j)
+ {
+ layfirst = j;
+ if (layerE[j] > 0.0) break;
+ }
+// if (layfirst > 7) continue;
+ if (phoD) aida.cloud1D("Layer of first photon cand interaction B").fill(layfirst);
+ // test with H-Matrix if cluster has enough cells
+ if (candclus.getCalorimeterHits().size()>=_mincells && layfirst<9)
+ {
+ for(int i=0; i<layerE.length-layfirst; ++i)
+ {
+ _valsB[i] = layerE[i+layfirst];
+ if (phoD) aida.cloud2D("Fractional Energy vs Layer").fill(i,layerE[i+layfirst]);
+ }
+ // Add the correlation to the log of the cluster energy
+ _valsB[_logEIndexB]=Math.log10(energy);
+ // have now filled the vector of measurements. need to either accumulate the HMatrix or apply it
+ if (_task==HMatrixTask.ANALYZE)
+ {
+ if (phoD) aida.cloud1D("nmeas").fill(_nmeasB);
+ double chisq = _hmx.chisquared(_valsB);
+ if (phoD) aida.cloud1D("Chisq").fill(chisq);
+ if (phoD) aida.cloud2D("Chisq vs energy").fill(energy,chisq);
+ if (phoD) aida.cloud1D("Chisq Probability").fill(ChisqProb.gammq(_nmeasB,chisq));
+ double chisqD = _hmx.chisquaredDiagonal(_valsB);
+ if (phoD) aida.cloud1D("ChisqD").fill(chisqD);
+ if (phoD) aida.cloud2D("ChisqD vs energy").fill(energy,chisqD);
+ double chiprobD = ChisqProb.gammq(_nmeasB,chisqD);
+ if(chiprobD<0.0000000001) chiprobD=0.0000000001;
+ if (phoD) aida.cloud1D("ChisqD Probability").fill(chiprobD);
+ double logchiprobD = Math.log10(chiprobD);
+ if (phoD) aida.cloud1D("Log ChisqD Probability").fill(logchiprobD);
+ if (phoD) aida.cloud2D("Chisq vs LogChisqD Prob").fill(logchiprobD,chisq);
+ double chisqND = chisq - chisqD;
+ if (phoD) aida.cloud1D("ChisqND").fill(chisqND);
+ if (logchiprobD > -9.9)
+ {
+ phoBclusters.add(ebclus);
+// System.out.println("Pho Cand Passed H-Matrix");
+ } else
+ {
+// System.out.println("Pho Cand Failed H-Matrix test");
+ int clsize = ebclus.getCalorimeterHits().size();
+ if (clsize>1 && clsize<_mincells && layfirst<5) phoBclusters.add(ebclus);
+ }
+ }
+ } else
+ {
+// System.out.println("No H-Matrix test try layfirst");
+// System.out.println("Num of hits " + candclus.getCalorimeterHits().size() + " First Layer " + layfirst);
+ int cansize = candclus.getCalorimeterHits().size();
+ if (cansize>1 && cansize<_mincells && layfirst<5) phoBclusters.add(ebclus);
+ }
+ } // end of barrel part
+ }
+ catch (java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No barrel clusters in photon finder");
+ }
+ try
+ {
+ List<BasicCluster> EECclusters = event.get(BasicCluster.class,"EECTMHitsFCClus");
+
+ for (BasicCluster eecclus : EECclusters)
+ {
+ List<BasicCluster> ennclusters = new ArrayList<BasicCluster>();
+ BasicCluster candclus = new BasicCluster();
+ double ecp[] = eecclus.getPosition();
+ double ecpx = ecp[0];
+ double ecpy = ecp[1];
+ double ecpz = ecp[2];
+ double ecR = Math.sqrt(ecpx*ecpx+ecpy*ecpy);
+ double ebclth = Math.atan(ecR/ecpz);
+ if (ebclth<0) ebclth+=Math.PI;
+ double ebclph = Math.atan2(ecpy,ecpx);
+ if (ebclph<0) ebclph+=2*Math.PI;
+ if (phoD) aida.cloud2D("Theta vs Phi of EEC Clusters").fill(ebclth,ebclph);
+ // check for tracks near clusters - some of these are interacting hadrons
+// System.out.println("Number of tracks to check with cands " + ntr);
+ int ntrclmatch = 0;
+ for (int i=0; i<ntr; ++i)
+ {
+ double trcldelth = Math.abs(ebclth-trth[i]);
+ if (phoD) aida.cloud1D("Clus-Trk theta diff in PhoFinder").fill(trcldelth);
+ double trcldelph = Math.abs(ebclph-trph[i]);
+ if (trcldelph > Math.PI) trcldelph = 2*Math.PI-trcldelph;
+ if (phoD) aida.cloud1D("Clus-Trk phi diff in PhoFinder").fill(trcldelph);
+ double trcldist = Math.sqrt(trcldelth*trcldelth+trcldelph*trcldelph);
+ if (phoD) aida.cloud1D("Trk-Clus Distance in PhoFinder").fill(trcldist);
+ // check for distance from a track - used 0.03, but .025 looks better
+ if (trcldist<_dTrcl) ntrclmatch++;
+ }
+ // if any matches, don't test this cluster - its not a photon
+ if (ntrclmatch>0) continue;
+ // Make a hitmap of cluser hits to re-cluster in NN clusterer
+ HitMap eechitmap = new HitMap();
+ for (CalorimeterHit eechit : eecclus.getCalorimeterHits())
+ {
+ long cID = eechit.getCellID();
+ eechitmap.put(cID, eechit);
+ }
+ int _dU = 1;
+ int _dV = 1;
+ int _dLayer = 1;
+ double ethresh = 0.;
+
+ for (;;)
+ {
+ if (eechitmap.isEmpty()) break;
+ long hID = eechitmap.keySet().iterator().next();
+ CalorimeterHit hit = eechitmap.get(hID);
+ NearestNeighborCluster nnclus = new NearestNeighborCluster(eechitmap, hit, hID, _dU, _dV, _dLayer, ethresh);
+ if(nnclus.getCalorimeterHits().size()>_mincells)
+ {
+ ennclusters.add(nnclus);
+// System.out.println("Formed NN cluster from FC cluster");
+ TotNNclusters.add(nnclus);
+ }
+ }
+ // if number of NN clusters is = 1, use it to test with H-Matrix, else use original FC cluster
+// System.out.println("Number of NN clusters " + ennclusters.size());
+ if (ennclusters.size() == 0)
+ {
+ candclus = eecclus;
+ } else if (ennclusters.size() == 1)
+ {
+ candclus = ennclusters.get(0);
+ } else if (ennclusters.size()>1)
+ {
+ candclus = eecclus;
+ }
+ CandClusters.add(candclus);
+ // check if this cluster passes H-Matrix test, if so, add candclus to list of Photon Clusters
+ double energy = candclus.getRawEnergy();
+ double[] layerE = layerEnergiesEC(candclus);
+ // accumulate the longitudinal energy fractions...
+ int layfirst = 0;
+ for(int j=0; j<layerE.length; ++j)
+ {
+ layfirst = j;
+ if (layerE[j] > 0.0) break;
+ }
+ if (phoD) aida.cloud1D("Layer of first photon cand interaction EC").fill(layfirst);
+ if (candclus.getCalorimeterHits().size()>=_mincells && layfirst<9)
+ {
+ for(int i=0; i<layerE.length-layfirst; ++i)
+ {
+ _valsEC[i] = layerE[i+layfirst];
+ if (phoD) aida.cloud2D("Fractional Energy vs Layer").fill(i,layerE[i+layfirst]);
+ }
+ // Add the correlation to the log of the cluster energy
+ _valsEC[_logEIndexEC]=Math.log10(energy);
+ // have now filled the vector of measurements. need to either accumulate the HMatrix or apply it
+ if (_task==HMatrixTask.ANALYZE)
+ {
+ if (phoD) aida.cloud1D("nmeas").fill(_nmeasEC);
+ double chisq = _hmx.chisquared(_valsEC);
+ if (phoD) aida.cloud1D("Chisq").fill(chisq);
+ if (phoD) aida.cloud2D("Chisq vs energy").fill(energy,chisq);
+ if (phoD) aida.cloud1D("Chisq Probability").fill(ChisqProb.gammq(_nmeasEC,chisq));
+ double chisqD = _hmx.chisquaredDiagonal(_valsEC);
+ if (phoD) aida.cloud1D("ChisqD").fill(chisqD);
+ if (phoD) aida.cloud2D("ChisqD vs energy").fill(energy,chisqD);
+ double chiprobD = ChisqProb.gammq(_nmeasEC,chisqD);
+ if(chiprobD<0.0000000001) chiprobD=0.0000000001;
+ if (phoD) aida.cloud1D("ChisqD Probability").fill(chiprobD);
+ double logchiprobD = Math.log10(chiprobD);
+ if (phoD) aida.cloud1D("Log ChisqD Probability").fill(logchiprobD);
+ if (phoD) aida.cloud2D("Chisq vs LogChisqD Prob").fill(logchiprobD,chisq);
+ double chisqND = chisq - chisqD;
+ if (phoD) aida.cloud1D("ChisqND").fill(chisqND);
+ if (logchiprobD > -9.9)
+ {
+ phoECclusters.add(eecclus);
+ } else
+ {
+ int clsize = eecclus.getCalorimeterHits().size();
+ if (clsize>1 && clsize<_mincells && layfirst<5) phoECclusters.add(eecclus);
+ }
+ }
+ } else // can't use H-Matrix, try something else
+ {
+ int cansize = candclus.getCalorimeterHits().size();
+ if (cansize>1 && cansize<_mincells && layfirst<5) phoECclusters.add(eecclus);
+ }
+ } // end of Endcap part
+ }
+ catch(java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("No endcap clusters in photon finder");
+ }
+ event.put("AllEMNNClusters",TotNNclusters);
+ event.put("PhoCandClusters",CandClusters);
+ event.put("PhotonBClusters",phoBclusters);
+ event.put("PhotonECClusters",phoECclusters);
+ double PFAPhoE = 0;
+ for (BasicCluster phoB : phoBclusters)
+ {
+ PFAPhoE += phoB.getEnergy()*1.015*1.013;
+// PFAPhoE += phoB.getEnergy();
+ aida.cloud1D("Number of hits in Photon B Clusters").fill(phoB.getCalorimeterHits().size());
+ }
+ for (BasicCluster phoEC : phoECclusters)
+ {
+ PFAPhoE += phoEC.getEnergy()*1.015*1.013;
+// PFAPhoE += phoEC.getEnergy();
+ aida.cloud1D("Number of hits in Photon EC Clusters").fill(phoEC.getCalorimeterHits().size());
+ }
+ if (phoRes && phoBclusters.size()+phoECclusters.size()==1) aida.cloud1D("Single Photon ESum").fill(PFAPhoE);
+ if (phoRes) aida.cloud1D("Total Photon ESum").fill(PFAPhoE);
+ if (phoRes) aida.cloud1D("Total Charged P + Photon E").fill(PFAPhoE+TotTrP);
+ } // end of filter loop
+ } // end of event loop
+ private double[] layerEnergiesB(BasicCluster clusB)
+ {
+ //FIXME could reuse this array
+ double[] layerEnergiesB = new double[_nLayersB];
+ double clusterEnergy = 0.;
+ List<CalorimeterHit> hits = clusB.getCalorimeterHits();
+// CalorimeterIDDecoder _decoderB = hit
+// CalorimeterIDDecoder _decoderB = (CalorimeterIDDecoder) event.getMetaData(chitsB).getIDDecoder();
+ //System.out.println("New cluster with "+hits.size()+ " hits"); // and energy "+clus.getEnergy());
+ for(CalorimeterHit hit : hits)
+ {
+ IDDecoder _decoderB = hit.getIDDecoder();
+ _decoderB.setID(hit.getCellID());
+ double e = hit.getRawEnergy();
+ int layer =_decoderB.getLayer();
+ //System.out.println("layer "+layer+" energy "+e);
+ clusterEnergy+=e;
+ layerEnergiesB[layer]+=e;
+ }
+
+ //System.out.println("Cluster energy = " + clusterEnergy + "\n");
+
+ for(int i=0; i<_nLayersB; ++i)
+ {
+ layerEnergiesB[i]/=clusterEnergy;
+ }
+ //FIXME sum of cell energies does not equal cluster energy!
+// System.out.println(clusterEnergy+" "+clus.getEnergy());
+ return layerEnergiesB;
+ }
+ private double[] layerEnergiesEC(BasicCluster clusEC)
+ {
+ //FIXME could reuse this array
+ double[] layerEnergiesEC = new double[_nLayersEC];
+ double clusterEnergy = 0.;
+ List<CalorimeterHit> hits = clusEC.getCalorimeterHits();
+ //System.out.println("New cluster with "+hits.size()+ " hits"); // and energy "+clus.getEnergy());
+ for(CalorimeterHit hit : hits)
+ {
+ IDDecoder _decoderEC = hit.getIDDecoder();
+ _decoderEC.setID(hit.getCellID());
+ double e = hit.getRawEnergy();
+ int layer =_decoderEC.getLayer();
+ //System.out.println("layer "+layer+" energy "+e);
+ clusterEnergy+=e;
+ layerEnergiesEC[layer]+=e;
+ }
+
+ //System.out.println("Cluster energy = " + clusterEnergy + "\n");
+
+ for(int i=0; i<_nLayersEC; ++i)
+ {
+ layerEnergiesEC[i]/=clusterEnergy;
+ }
+ //FIXME sum of cell energies does not equal cluster energy!
+// System.out.println(clusterEnergy+" "+clus.getEnergy());
+ return layerEnergiesEC;
+ }
+}
+
