lcsim/src/org/lcsim/contrib/compile/SteveMagill
diff -u -r1.3 -r1.4
--- BTrMipClus.java 2 Aug 2005 22:33:07 -0000 1.3
+++ BTrMipClus.java 9 Aug 2005 01:38:42 -0000 1.4
@@ -24,6 +24,7 @@
// import org.lcsim.recon.ztracking.FoundTrack;
// import org.lcsim.geometry.Detector;
+
public class BTrMipClus extends Driver
{
private AIDA aida = AIDA.defaultInstance();
@@ -36,6 +37,7 @@
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
@@ -103,19 +105,21 @@
} // end loop over had layers
// Get (set) sampling fractions here
-// _sfEMB = 0.01207; // for now this is tested sampling fraction for W/Si EMCAL in SiD Detector
- _sfEMB = 0.01167; // for now, this is effective sf for 20/10 Si/W ECAL from 5 GeV gammas
+// _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 = 6.2; // 161 MeV per hit in HCAL for CDC WRPC
-// _sfHB = 24.0; // 42 MeV per hit in HCAL for CDC WScin after 1/3 mip thresold cut
-// _sfHB = 0.06;
+// _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 HADBESum=0;
+ double EMBESum=0; double EMBESumC=0; double HADBESum=0;
+ double EEMlowE=0; double EEMhighE=0;
double ESumBTotCAL=0;
double hthrshld = 0;
int nhitsBTotCAL=0;
@@ -143,12 +147,18 @@
celbid[nhitbLayer[emblayer]][emblayer] = hitID;
nhitbLayer[emblayer]++;
embhitmap.put(hitID, embhit);
-// if (embhit.getEnergy()<0.00061)
-// {
-// aida.cloud1D("DIAG : EM B Calib hit E").fill(embhit.getEnergy()/_sfEMB);
-// }
-// EMBESum += embhit.getEnergy();
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());
@@ -157,30 +167,23 @@
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.0003; // 1/3 mip threshold for cdc Scin HCAL
- System.out.println("HAD Threshold = " +hthrshld);
+ if (event.getDetectorName().equals("cdcaug05")) hthrshld = 0.0002; // 1/3 mip threshold for cdc Scin HCAL
+// System.out.println("HAD Threshold = " +hthrshld);
for (SimCalorimeterHit hadbhit : hadbHits)
{
long hitID = hadbhit.getCellID();
hadbdecoder.setID(hitID);
int hadblayer = hadbdecoder.getLayer();
// if (hadbhit.getEnergy() > hthrshld)
- aida.cloud2D("HAD Cell E vs Time").fill(hadbhit.getTime(),hadbhit.getRawEnergy());
+ aida.cloud1D("HAD Cell Vis E before thr cut").fill(hadbhit.getRawEnergy());
if (hadbhit.getRawEnergy() > hthrshld)
{
+ 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);
-// aida.cloud1D("HAD B hit vis E").fill(hadbhit.getEnergy());
-// if (hadbhit.getEnergy()<0.000004)
-// {
-// aida.cloud1D("DIAG : HAD B Vis hit E").fill(hadbhit.getEnergy());
-// aida.cloud1D("DIAG : HAD B Calib hit E").fill(hadbhit.getEnergy()/_sfHB);
-// }
-// aida.cloud2D("HAD B cell E vs Time").fill(hadbhit.getTime(),hadbhit.getEnergy());
-// HADBESum += hadbhit.getEnergy();
HADBESum += hadbhit.getRawEnergy();
nhitsHADB++;
}
@@ -228,11 +231,13 @@
}
ESumBTotCAL = EMBESum + HADBESum;
nhitsBTotCAL = nhitsEMB + nhitsHADB;
- aida.cloud2D("NHits EM B vs EM B Calib ESum").fill(EMBESum/_sfEMB,nhitsEMB);
+ 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);
- if (EMBESum/_sfEMB < 0.40)
+// 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(EMBESum/_sfEMB);
+ 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);
@@ -249,9 +254,9 @@
}
aida.cloud1D("Total B CAL NHits").fill(nhitsBTotCAL);
aida.cloud1D("Total B CAL Vis ESum").fill(ESumBTotCAL);
- aida.cloud1D("Total B CAL EM Analog + HAD Digital ESum").fill(EMBESum/_sfEMB+nhitsHADB/_sfHB);
+ 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(EMBESum/_sfEMB);
+ aida.cloud1D("Total EM B Calib ESum").fill(EMBESumC);
aida.cloud1D("Total HAD B Vis ESum").fill(HADBESum);
/** This section finds mip clusters in cal attached to a reconstructed (or true) track.
@@ -638,7 +643,7 @@
// aida.cloud1D("Num of hits in had cluster").fill(nnclus.getCalorimeterHits().size());
}
}
- if (hbclusters.size() > 0) event.put("HADBNNClusters",hbclusters);
+
// map hits in cluster to their cluster for combining later
Map<CalorimeterHit, Cluster> hbnnhitclus = new HashMap<CalorimeterHit, Cluster>();
for (Cluster hbcluster : hbclusters)
@@ -721,7 +726,6 @@
}
}
if (phoclusters.size() > 0) event.put("PhoBClusters",phoclusters);
- if (ebshowclusters.size() > 0) event.put("EBShowClusters",ebshowclusters);
double PhoClusE = 0;
aida.cloud1D("Num Photon B Clusters").fill(phoclusters.size());
@@ -730,15 +734,7 @@
PhoClusE += phocluster.getEnergy()/_sfEMB;
}
aida.cloud1D("Photon Cluster E").fill(PhoClusE);
-
- 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);
-
+
// make a map linking hits to EBShower clusters
Map<CalorimeterHit, Cluster> ebshowhitclus = new HashMap<CalorimeterHit, Cluster>();
for (Cluster ebshowcluster : ebshowclusters)
@@ -769,6 +765,7 @@
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]);
@@ -787,9 +784,11 @@
aida.cloud1D("Mip E per Track").fill(trmipE);
aida.cloud1D("Track E over P after adding mips").fill(trmipE/trbpp[i]);
}
- for (Cluster ebshowcluster : ebshowclusters)
+ do {
+ for (Iterator<Cluster> ieclus = ebshowclusters.iterator(); ieclus.hasNext();)
{
- List<CalorimeterHit> ebshowhits = ebshowcluster.getCalorimeterHits();
+ Cluster eclus = ieclus.next();
+ List<CalorimeterHit> ebshowhits = eclus.getCalorimeterHits();
int ntsh=0;
double tsd=99;
for (CalorimeterHit ebshowhit : ebshowhits)
@@ -805,21 +804,24 @@
double tsdelph = Math.abs(shhph-trbphi[i][IntLay[i]]);
if (tsdelph > Math.PI) tsdelph = 2*Math.PI-tsdelph;
tsd = Math.sqrt(tsdelth*tsdelth+tsdelph*tsdelph);
- aida.cloud1D("Track-Cell Distance in Track EM Shower Cluster").fill(tsd);
- if (tsd<0.075) ntsh++;
+ 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(ebshowcluster);
- treshE += ebshowcluster.getEnergy()/_sfEMB;
+ 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, ebshowcluster);
+ treshmap.put(i, eclus);
+ ieclus.remove();
}
}
- for (Cluster hbcluster : hbclusters)
+ for (Iterator<Cluster> ihclus = hbclusters.iterator(); ihclus.hasNext();)
{
- List<CalorimeterHit> hbhits = hbcluster.getCalorimeterHits();
+ Cluster hclus = ihclus.next();
+ List<CalorimeterHit> hbhits = hclus.getCalorimeterHits();
+
int nthsh=0;
double thsd=99;
for (CalorimeterHit hbhit : hbhits)
@@ -835,46 +837,69 @@
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);
- aida.cloud1D("Track-shower Distance in Track HAD Shower Cluster").fill(thsd);
- if (thsd<0.15) nthsh++;
+ 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(hbcluster);
- trhshE += hbcluster.getCalorimeterHits().size()/_sfHB;
+ 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, hbcluster);
+ 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;
- aida.cloud1D("Track E over P").fill(trtotE/trbpp[i]);
+ 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());
-
- // Compare E of Clusters to P of Track
- for (Integer itrack : itracks)
+ // 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)
{
- 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;
- if (treshmap.get(itrack) != null) EMshowE = treshmap.get(itrack).getEnergy()/_sfEMB;
- if (trhshmap.get(itrack) != null) HshowE = trhshmap.get(itrack).getCalorimeterHits().size()/_sfHB;
- aida.cloud1D("EM mip E per track").fill(EMmipE);
- aida.cloud1D("HAD mip E per track").fill(HmipE);
- aida.cloud1D("EM shower E per track").fill(EMshowE);
- aida.cloud1D("HAD shower E per track").fill(HshowE);
- 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);
- }
- }
+ 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;
@@ -922,8 +947,8 @@
}
}
}
-// if (neuE+charE+phoE > 75)
-// {
+ 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);
@@ -944,7 +969,7 @@
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
}
}
