Print

Print


Commit in lcsim/src/org/lcsim/contrib/compile/SteveMagill on MAIN
BTrMipClus.java+191-781.2 -> 1.3


lcsim/src/org/lcsim/contrib/compile/SteveMagill
BTrMipClus.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- BTrMipClus.java	2 Aug 2005 17:18:02 -0000	1.2
+++ BTrMipClus.java	2 Aug 2005 22:33:07 -0000	1.3
@@ -1,3 +1,10 @@
+//  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.jet.JetFinder;
@@ -96,9 +103,13 @@
         }  // 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.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
 //        _sfHB = 0.000025;   //  for now this is tested sampling fraction for SS/RPC HCAL in SiD Detector
-        _sfHB = 0.06;
+        //  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;
 
         _initialized = true;
     }   // end of initialization section
@@ -106,6 +117,7 @@
       // get calorimeter hits by layer - need to calculate hit density for cells
       double EMBESum=0; double HADBESum=0;
       double ESumBTotCAL=0;
+      double hthrshld = 0;
       int nhitsBTotCAL=0;
       double[][] celbtheta = new double[5000][_numbemlayers+_numbhadlayers];
       double[][] celbphi = new double[5000][_numbemlayers+_numbhadlayers];
@@ -135,6 +147,7 @@
 //          {
 //              aida.cloud1D("DIAG : EM B Calib hit E").fill(embhit.getEnergy()/_sfEMB);     
 //          }
+//          EMBESum += embhit.getEnergy();
           EMBESum += embhit.getRawEnergy();
           nhitsEMB++;    
       }
@@ -144,24 +157,33 @@
       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);
       for (SimCalorimeterHit hadbhit : hadbHits)
       {
           long hitID = hadbhit.getCellID();
           hadbdecoder.setID(hitID);
           int hadblayer = hadbdecoder.getLayer();
-          celbtheta[nhitbLayer[_numbemlayers+hadblayer]][_numbemlayers+hadblayer] = hadbdecoder.getTheta();				  
-          celbphi[nhitbLayer[_numbemlayers+hadblayer]][_numbemlayers+hadblayer] = hadbdecoder.getPhi();
-          celbid[nhitbLayer[_numbemlayers+hadblayer]][_numbemlayers+hadblayer] = hitID;
-          nhitbLayer[_numbemlayers+hadblayer]++; 
-          hadbhitmap.put(hitID, hadbhit);          
+//          if (hadbhit.getEnergy() > hthrshld)
+          aida.cloud2D("HAD Cell E vs Time").fill(hadbhit.getTime(),hadbhit.getRawEnergy());
+          if (hadbhit.getRawEnergy() > hthrshld)
+          {
+              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);
 //          }
-          HADBESum += hadbhit.getRawEnergy();
-          nhitsHADB++;
+//          aida.cloud2D("HAD B cell E vs Time").fill(hadbhit.getTime(),hadbhit.getEnergy());
+//              HADBESum += hadbhit.getEnergy();
+              HADBESum += hadbhit.getRawEnergy();
+              nhitsHADB++;
+          }
       }
       aida.cloud1D("Nhits in original HADBhitmap").fill(hadbhitmap.size());
 
@@ -213,11 +235,13 @@
           aida.cloud1D("EM B Calib ESum, mips only").fill(EMBESum/_sfEMB);
           aida.cloud1D("NHits EM B, mips only").fill(nhitsEMB); 
           aida.cloud1D("NHits HAD B, EM B mips only").fill(nhitsHADB);
+          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)
       {
@@ -225,7 +249,7 @@
       }
       aida.cloud1D("Total B CAL NHits").fill(nhitsBTotCAL);     
       aida.cloud1D("Total B CAL Vis ESum").fill(ESumBTotCAL);
-      aida.cloud1D("Total B CAL EM Analog + HAD Digital ESum").fill(EMBESum/_sfEMB+nhitsHADB*0.115*0.2);
+      aida.cloud1D("Total B CAL EM Analog + HAD Digital ESum").fill(EMBESum/_sfEMB+nhitsHADB/_sfHB);
       aida.cloud1D("Total EM B Vis ESum").fill(EMBESum);
       aida.cloud1D("Total EM B Calib ESum").fill(EMBESum/_sfEMB);
       aida.cloud1D("Total HAD B Vis ESum").fill(HADBESum);
@@ -257,8 +281,8 @@
       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; 
+          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();
@@ -330,7 +354,10 @@
       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];
@@ -339,12 +366,12 @@
           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>(); 
+          Map<Long, CalorimeterHit> hbclhitmap = new HashMap<Long, CalorimeterHit>();
           BasicCluster ebmclus = new BasicCluster();
           BasicCluster hbmclus = new BasicCluster();
           BasicCluster bmipclus = new BasicCluster();
           boolean nointeraction = true;
-          IntLay[i] = 64;     // set to maximum for each track to start
+          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)
           {
@@ -410,7 +437,6 @@
 //                      if (j == 0 && prden == 0) aida.cloud1D("Track pT for matches in Layer 0").fill(trbppt[i]);
 //                      if (j == 0 && prden == 0) aida.cloud1D("Track theta for matches in Layer 0").fill(trbtheta[i][j]);
                       if (j == 0) continue;    // if interaction in first layer, just get out
-                      trmipclusmap.put(i, bmipclus);
                       if (j < _numbemlayers)
                       {
                           // only have EM cluster
@@ -419,6 +445,7 @@
                           {
                               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
                               }
@@ -432,6 +459,7 @@
                           {
                               for (CalorimeterHit eh : ebclhitmap.values())
                               {
+                                  if (eh == null) continue;
                                   ebmclus.addHit(eh);
                                   bmipclus.addHit(eh);
                               }
@@ -439,18 +467,22 @@
                           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 (hbclhitmap.size() > 0)
                               {
-                                  hbmclus.addHit(hh);
-                                  bmipclus.addHit(hh);
+                                  for (CalorimeterHit hh : hbclhitmap.values())
+                                  {
+                                      if (hh == null) continue;
+                                      hbmclus.addHit(hh);
+                                      bmipclus.addHit(hh);
+                                  }
                               }
-                          }
-                          hadbmipclusters.add(hbmclus);
+                              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
@@ -473,7 +505,32 @@
                           }
                       }
                   }
-              }
+                  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
@@ -664,14 +721,23 @@
           }
       }
       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());
       for (Cluster phocluster : phoclusters)
       {
           PhoClusE += phocluster.getEnergy()/_sfEMB;
       }
-      aida.cloud1D("Num Photon B Clusters").fill(phoclusters.size());
-      if (ebshowclusters.size() > 0) event.put("EBShowClusters",ebshowclusters);
+      aida.cloud1D("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>();
@@ -690,39 +756,37 @@
       //  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;
+//          aida.cloud1D("Track Momentum").fill(trbpp[i]);
           //  first, check for track-cell match with a mip cluster at last mip layer (IntLay[i]-1)
-          if (IntLay[i] > 0) tcclus.addCluster(trmipclusmap.get(i));
-//              for (BasicCluster calbmipcluster : calbmipclusters)
-//              {
-//                  if (IntLay[i] == 0) continue;
-                  // compare cluster hits with track position
-//                  List<CalorimeterHit> clhits = calbmipcluster.getCalorimeterHits();
-//                  int ntc=0;
-//                  double tcd=99;
-//                  for (CalorimeterHit clhit : clhits)
-//                  {
-                      // compare hit-track distance
-//                      double[] clhpos = clhit.getPosition();
-//                      double clhR = Math.sqrt(clhpos[0]*clhpos[0]+clhpos[1]*clhpos[1]);
-//                      double clhth = Math.atan(clhR/clhpos[2]);
-//                      if(clhth<0) clhth+=Math.PI;
-//                      double clhph = Math.atan2(clhpos[1],clhpos[0]);
-//	              if(clhph<0) clhph+=2*Math.PI;
-//                      double tcdelth = Math.abs(clhth-trbtheta[i][IntLay[i]-1]);
-//                      double tcdelph = Math.abs(clhph-trbphi[i][IntLay[i]-1]);
-//                      if (tcdelph > Math.PI) tcdelph = 2*Math.PI-tcdelph;
-//                      tcd = Math.sqrt(tcdelth*tcdelth+tcdelph*tcdelph);
-//                      aida.cloud1D("Track-Cell Distance in Track Mip Cluster").fill(tcd);
-//                      if (tcd<0.01) ntc++;
-//                  }
-//                  aida.cloud1D("Number of cell matches per cluster").fill(ntc);
-//                  if (ntc > 0) tcclus.addCluster(calbmipcluster);
-//              }
+//          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]);
+              }
               for (Cluster ebshowcluster : ebshowclusters)
               {
                   List<CalorimeterHit> ebshowhits = ebshowcluster.getCalorimeterHits();
@@ -742,10 +806,16 @@
                       if (tsdelph > Math.PI) tsdelph = 2*Math.PI-tsdelph;
                       tsd = Math.sqrt(tsdelth*tsdelth+tsdelph*tsdelph);
                       aida.cloud1D("Track-Cell Distance in Track EM Shower Cluster").fill(tsd);
-                      if (tsd<0.175) ntsh++;
+                      if (tsd<0.075) ntsh++;
                   }
                   aida.cloud1D("Number of cell matches per EM shower").fill(ntsh);
-                  if (ntsh > 0) tcclus.addCluster(ebshowcluster);
+                  if (ntsh > 0) 
+                  {
+                      tcclus.addCluster(ebshowcluster);
+                      treshE += ebshowcluster.getEnergy()/_sfEMB;
+                      aida.cloud1D("Track E over P after adding EM shower").fill((trmipE+treshE)/trbpp[i]);
+                      treshmap.put(i, ebshowcluster);
+                  }
               }
               for (Cluster hbcluster : hbclusters)
               {
@@ -761,30 +831,60 @@
                       if(hhth<0) hhth+=Math.PI;
                       double hhph = Math.atan2(hhpos[1],hhpos[0]);
 	              if(hhph<0) hhph+=2*Math.PI;
-                      double thdelth = Math.abs(hhth-trbtheta[i][_numbemlayers+1]);
-                      double thdelph = Math.abs(hhph-trbphi[i][_numbemlayers+1]);
+                      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);
                       aida.cloud1D("Track-shower Distance in Track HAD Shower Cluster").fill(thsd);
-                      if (thsd<0.175) nthsh++;
+                      if (thsd<0.15) nthsh++;
                   }
                   aida.cloud1D("Number of cell matches per HAD shower").fill(nthsh);
-                  if (nthsh > 0) tcclus.addCluster(hbcluster);
+                  if (nthsh > 0) 
+                  {
+                      tcclus.addCluster(hbcluster);
+                      trhshE += hbcluster.getCalorimeterHits().size()/_sfHB;
+                      aida.cloud1D("Track E over P after adding HAD shower").fill((trmipE+treshE+trhshE)/trbpp[i]);
+                      trhshmap.put(i, hbcluster);
+                  }
               }
+          }
+          trtotE += trmipE+treshE+trhshE;
+          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);      
+      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)
+      {
+          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);
+          }
+      }   
 
       //  get some MC info to compare
       int npho=0;  int nneu=0;  int nchar=0;
-      double phoE = 0;  double neuE = 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; 
+          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();
@@ -815,23 +915,36 @@
           } else
           {
               // its charged
-              if (mcppt > 0.9525) nchar++;
+              if (mcppt > 0.9525)
+              {
+                  nchar++;
+                  charE += mcpE;
+              }
           }
       }
-      aida.cloud1D("Number of MC Photons").fill(npho);
-      aida.cloud1D("True Photon E").fill(phoE);
-      aida.cloud1D("Photon Cluster E").fill(PhoClusE);
-      aida.cloud1D("True Photon E - Photon Cluster E").fill(phoE-PhoClusE);
-      aida.cloud1D("Number of True Photons - Photon Clusters").fill(npho-phoclusters.size());
-      aida.cloud1D("Number of Neutrals").fill(nneu);
-      aida.cloud1D("Neutral E Sum").fill(neuE);
-      aida.cloud1D("Number of charged Particles").fill(nchar);
-      aida.cloud2D("Num MC Photons vs Num Photon Clusters").fill(npho,phoclusters.size());
-      aida.cloud2D("Num MC charged Particles vs Num EM B Mip Clusters").fill(ntrkswim,embmipclusters.size());
-      aida.cloud2D("Num MC charged Particles vs Num EM B Mip Clusters").fill(nchar,embmipclusters.size());
-      aida.cloud2D("Num Swum charged Particles vs Num EM B Mip Clusters").fill(ntrkswim,embmipclusters.size());
-      aida.cloud2D("Num MC charged Particles vs Num Swum charged").fill(nchar,ntrkswim);
-         
+//      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   
    }
 }
CVSspam 0.2.8