lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill
diff -N ClPointDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ClPointDriver.java 12 Aug 2009 18:35:07 -0000 1.1
@@ -0,0 +1,522 @@
+package org.lcsim.contrib.SteveMagill;
+
+// Takes a cluster list and determines the best option for pointing of each
+// cluster. In this case, clusters are grouped according to whether they
+// point at the IP, the IL Spacepoint of a track mip match, or neither. Later,
+// tests may be added to determine if clusters point to other clusters.
+// input : calorimeter clusters
+// output : separated cluster lists with pointing properties
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.util.Driver;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.aida.*;
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+import org.lcsim.spacegeom.*;
+
+public class ClPointDriver extends Driver
+{
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IAnalysisFactory analysisFactory = IAnalysisFactory.create();
+ private AIDA aida = AIDA.defaultInstance();
+ private double _mindIP;
+ private double _mindILSP;
+ private double _mindcl;
+ private String _clustername;
+ private String[] _oclnames;
+ private boolean clpdb = false;
+
+ public ClPointDriver(double mindIP, double mindILSP, double mindcl)
+ {
+ // add arguments if needed
+ _mindIP = mindIP; // minimum distance at IP
+ _mindILSP = mindILSP; // minimum distance at IL SP
+ _mindcl = mindcl; // minimum distance for IP PhoTrack def
+ }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event); // executes all added drivers
+
+ // get cluster list from event
+ List<BasicCluster> clusters = event.get(BasicCluster.class,_clustername);
+
+ // create lists of clusters with common pointing
+ List<BasicCluster> IPclus = new ArrayList<BasicCluster>(); // points at IP
+ List<BasicCluster> ILSPclus = new ArrayList<BasicCluster>(); // points at ILSP
+ List<BasicCluster> NPclus = new ArrayList<BasicCluster>(); // no pointing
+
+ // Get maps linking IL and mip clusters to tracks, tracks linked to SpacePoints
+ Map<Track, BasicCluster> trmipmap = (Map<Track, BasicCluster>) event.get("TrackMipMap");
+ Map<Track, SpacePoint> trilmap = (Map<Track, SpacePoint>) event.get("TrackILPosMap");
+
+ // Make a new map linking ILSP clusters to tracks
+ Map<Track, BasicCluster> trclmap = new HashMap<Track, BasicCluster>();
+ Map<BasicCluster, Track> cltrmap = new HashMap<BasicCluster, Track>();
+
+ // loop over all clusters, determining pointing properties
+ for (BasicCluster clus : clusters)
+ {
+ int clindex = 0;
+ // first, make sure this cluster has enough hits, if not, put in
+ // nopoint list must be at least 4 hits
+ if (clus.getCalorimeterHits().size()>3)
+ {
+ // continue with analysis, get this cluster's position, Ith, Iph
+ // make into vector of x,y,z
+ double[] clpos = clus.getPosition();
+ double clrxyz = Math.sqrt(clpos[0]*clpos[0]+clpos[1]*clpos[1]+clpos[2]*clpos[2]);
+ double thdir = clus.getITheta();
+ if (thdir<0) thdir+=Math.PI;
+ double phdir = clus.getIPhi();
+ if (phdir<0) phdir+=2*Math.PI;
+ double[] clvec = new double[3];
+ clvec[0] = clrxyz*Math.sin(thdir)*Math.cos(phdir);
+ clvec[1] = clrxyz*Math.sin(thdir)*Math.sin(phdir);
+ clvec[2] = clrxyz*Math.cos(thdir);
+ aida.cloud1D("Cluster ITheta").fill(thdir); // should be 0 to pi
+ aida.cloud1D("Cluster IPhi").fill(phdir); // should be 0 to 2pi
+ // also get IP direction vector using cluster position
+ double IPph = Math.atan2(clpos[1],clpos[0]);
+ if (IPph<0) IPph+=2*Math.PI;
+ double ccr = Math.sqrt(clpos[0]*clpos[0]+clpos[1]*clpos[1]);
+ double IPth = Math.atan(ccr/clpos[2]);
+ if (IPth<0) IPth+=Math.PI;
+ aida.cloud1D("IP dir Theta").fill(IPth);
+ aida.cloud1D("IP dir Phi").fill(IPph);
+ double[] ipvec = new double[3];
+ ipvec[0] = clpos[0];
+ ipvec[1] = clpos[1];
+ ipvec[2] = clpos[2];
+ double dth = Math.abs(thdir-IPth);
+ aida.cloud1D("Diff Ith IPth").fill(dth);
+ double dph = Math.abs(phdir-IPph);
+ if (dph>Math.PI) dph = 2*Math.PI-dph;
+ aida.cloud1D("Diff Iph IPph").fill(dph);
+ aida.cloud2D("Diff theta phi").fill(dph,dth);
+ double dthph = Math.sqrt(dth*dth+dph*dph);
+ aida.cloud1D("Diff dthph").fill(dthph);
+ double magcl = Math.sqrt(clvec[0]*clvec[0]+clvec[1]*clvec[1]+clvec[2]*clvec[2]);
+ double magip = Math.sqrt(ipvec[0]*ipvec[0]+ipvec[1]*ipvec[1]+ipvec[2]*ipvec[2]);
+ // check with ILSPs, keeping best match to track
+// double bestdis = 999.;
+ double bestdis = 0.;
+ for (Track itr : trilmap.keySet())
+ {
+ // get associated IL coordinates from IL map
+ SpacePoint trIL = trilmap.get(itr);
+ double trilx = trIL.x();
+ double delx = trilx-clpos[0];
+ double trily = trIL.y();
+ double dely = trily-clpos[1];
+ double delr = Math.sqrt(delx*delx+dely*dely);
+ double trilz = trIL.z();
+ double delz = trilz-clpos[2];
+// double trxyzil = trIL.rxyz(); // rxyz of track spacepoint at IL
+ // get theta-phi of IL SP
+ double trilth = Math.atan(trIL.rxy()/trIL.z());
+ if (trilth<0) trilth+=Math.PI;
+ double trilph = Math.atan2(trIL.y(),trIL.x());
+ if (trilph<0) trilph+=2*Math.PI;
+ // check if cluster close to this IL in theta-phi
+ double dist = 999;
+ double dccth = Math.abs(trilth-IPth);
+ double dccph = Math.abs(trilph-IPph);
+ if (dccph>Math.PI) dccph = 2*Math.PI-dccph;
+ dist = Math.sqrt(dccth*dccth+dccph*dccph);
+ if (dist>_mindILSP) continue;
+// double trtheta = trIL.theta(); // track theta at IL
+// double trphi = trIL.phi(); // track phi at IL
+ double clilth = Math.atan(delr/delz);
+ if (clilth<0) clilth+=Math.PI;
+ aida.cloud1D("Cl IL Theta").fill(clilth);
+ double clilph = Math.atan2(-dely,-delx);
+ if (clilph<0) clilph+=2.*Math.PI;
+ aida.cloud1D("Cl IL Phi").fill(clilph);
+ // now get difference between cluster I and IL
+ double dilth = Math.abs(thdir-clilth);
+ aida.cloud1D("Diff Ith ILSPth").fill(dilth);
+ double dilph = Math.abs(phdir-clilph);
+ if (dilph>Math.PI) dilph = 2*Math.PI-dilph;
+ aida.cloud1D("Diff Iph ILSPph").fill(dilph);
+ aida.cloud2D("Cl IPhi vs ILSPPhi").fill(clilph,phdir);
+ double[] ilvec = new double[3];
+// ilvec[0] = clrxyz*Math.sin(clilth)*Math.cos(clilph);
+// ilvec[1] = clrxyz*Math.sin(clilth)*Math.sin(clilph);
+// ilvec[2] = clrxyz*Math.cos(clilth);
+ ilvec[0] = delx;
+ ilvec[1] = dely;
+ ilvec[2] = delz;
+ double magil = Math.sqrt(ilvec[0]*ilvec[0]+ilvec[1]*ilvec[1]+ilvec[2]*ilvec[2]);
+ // dot products
+ Hep3Vector cl3v = new BasicHep3Vector(clvec);
+ Hep3Vector il3v = new BasicHep3Vector(ilvec);
+ Hep3Vector ip3v = new BasicHep3Vector(ipvec);
+ double dotiil = VecOp.dot(cl3v, il3v);
+ double dotiip = VecOp.dot(cl3v, ip3v);
+ double cosiil = Math.abs(dotiil/(magcl*magil));
+ double cosiip = Math.abs(dotiip/(magcl*magip));
+ aida.cloud1D("CosDel of IclIL vectors").fill(cosiil);
+ aida.cloud1D("CosDel of IclIP vectors").fill(cosiip);
+ aida.cloud1D("Diff cdil cdip").fill(cosiil-cosiip);
+ aida.cloud2D("Diff ILSP theta phi").fill(dilph,dilth);
+ double dilthph = Math.sqrt(dilth*dilth+dilph*dilph);
+ aida.cloud1D("Diff dilthph").fill(dilthph);
+ aida.cloud1D("Diff IP IL").fill(dthph-dilthph);
+ aida.cloud2D("IP vs IL diff").fill(dilthph,dthph);
+ aida.cloud2D("DiffIIL vs CosDel for all clus").fill(cosiil,dilthph);
+// if (dilthph<dthph && dilthph<bestdis) // points at ILSP more than IP
+ if (cosiil>cosiip && cosiil>bestdis)
+ {
+ aida.cloud1D("Diff I and IL for ILSP Clus").fill(dilthph);
+ aida.cloud1D("CosDel for ILSP Clus").fill(cosiil);
+ aida.cloud2D("DiffIIL vs CosDel for ILSP clus").fill(cosiil,dilthph);
+// bestdis = dilthph;
+ bestdis = cosiil;
+ cltrmap.put(clus, itr);
+ clindex++;
+// ILSPclus.add(clus);
+// trclmap.put(itr, clus);
+ }
+// if (dilthph<dthph) break;
+ } // end of track loop
+ if (clindex>0)
+ {
+ ILSPclus.add(clus);
+ }
+// System.out.println(" AfterTr Cl Index " + clindex);
+// if (clindex>0) continue;
+ if (clindex==0 && dthph<_mindIP) clindex--;
+ }
+ // check index, sort clusters
+// System.out.println(" Cl Index " + clindex);
+ if (clindex==0) NPclus.add(clus);
+// if (clindex>0) ILSPclus.add(clus);
+ if (clindex<0) IPclus.add(clus);
+ } // end of cluster loop
+
+ // have clusters defined into 3 classes - pointing to the interaction layer
+ // spacepoint, pointing to the IP, non-pointing
+ event.put("TrackILSPMap",trclmap);
+ event.put("ClusILSPMap",cltrmap);
+// if (NPclus.size()>0) event.put("NPClusters",NPclus);
+ event.put("NPClusters",NPclus);
+ aida.cloud1D("Number of NoPoint Clus").fill(NPclus.size());
+// if (IPclus.size()>0) event.put("IPClusters",IPclus);
+ event.put("IPClusters",IPclus);
+ aida.cloud1D("Number of IPPoint Clus").fill(IPclus.size());
+// if (ILSPclus.size()>0) event.put("ILSPClusters",ILSPclus);
+ event.put("ILSPClusters",ILSPclus);
+ aida.cloud1D("Number of ILSPPoint Clus").fill(ILSPclus.size());
+
+ // check the IP pointing clusters - divide into close/far from track using
+ // distance measure
+ List<BasicCluster> IPtrclus = new ArrayList<BasicCluster>(); // matched to tracks
+ List<BasicCluster> IPphclus = new ArrayList<BasicCluster>(); // not matched to tracks
+
+ for (BasicCluster ipcl : IPclus)
+ {
+ int ipclind = 0;
+ double[] clpos = ipcl.getPosition();
+ double clph = Math.atan2(clpos[1],clpos[0]);
+ if (clph<0) clph+=2*Math.PI;
+ double ccr = Math.sqrt(clpos[0]*clpos[0]+clpos[1]*clpos[1]);
+ double clth = Math.atan(ccr/clpos[2]);
+ if (clth<0) clth+=Math.PI;
+ // compare with extrapolated track position
+ for (Track itr : trilmap.keySet())
+ {
+ // get associated IL coordinates from IL map
+ SpacePoint trIL = trilmap.get(itr);
+ double trilth = Math.atan(trIL.rxy()/trIL.z());
+ if (trilth<0) trilth+=Math.PI;
+ double trilph = Math.atan2(trIL.y(),trIL.x());
+ if (trilph<0) trilph+=2*Math.PI;
+ // check if cluster close to this IL in theta-phi
+ double dist = 999;
+ double dccth = Math.abs(trilth-clth);
+ double dccph = Math.abs(trilph-clph);
+ if (dccph>Math.PI) dccph = 2*Math.PI-dccph;
+ dist = Math.sqrt(dccth*dccth+dccph*dccph);
+ if (dist<_mindcl) // best combination is 0.07 for ZZ in SiD01
+ {
+ ipclind++;
+ }
+ } // end of track loop
+ // put cluster in proper list
+ if (ipclind>0) IPtrclus.add(ipcl);
+ if (ipclind==0) IPphclus.add(ipcl);
+ } // end of cluster loop
+
+ // now have 2 classes of IP pointing clusters - those that are photon-like
+ // and those that are track-like as determined only by distance to closest track
+// if (IPtrclus.size()>0) event.put("IPTrackClusters",IPtrclus);
+ event.put("IPTrackClusters",IPtrclus);
+ aida.cloud1D("Number of IPPoint Track Clus").fill(IPtrclus.size());
+// if (IPphclus.size()>0) event.put("IPPhoClusters",IPphclus);
+ event.put("IPPhoClusters",IPphclus);
+ aida.cloud1D("Number of IPPoint Pho Clus").fill(IPphclus.size());
+
+ // check pointing of IP and NP clusters to ILSP clusters - try to
+ // link fragments of charged hadrons
+ List<BasicCluster> NPILclus = new ArrayList<BasicCluster>(); // NP clusters pointing to ILSP clusters
+
+ for (BasicCluster npc : NPclus)
+ {
+ int iclilpt = 0;
+ double distclil = 999.;
+ // get this cluster's pointing directions and positions
+ double[] clpos = npc.getPosition();
+// Hep3Vector npcl = new BasicHep3Vector(clpos);
+ double thdir = npc.getITheta();
+ if (thdir<0) thdir+=Math.PI;
+ double phdir = npc.getIPhi();
+ if (phdir<0) phdir+=2*Math.PI;
+ aida.cloud1D("NPCluster ITheta").fill(thdir); // should be 0 to pi
+ aida.cloud1D("NPCluster IPhi").fill(phdir); // should be 0 to 2pi
+ // get x,y,z for local pointing
+ double npclrxyz = Math.sqrt(clpos[0]*clpos[0]+clpos[1]*clpos[1]+clpos[2]*clpos[2]);
+ double[] nppvec = new double[3];
+ nppvec[0] = npclrxyz*Math.sin(thdir)*Math.cos(phdir);
+ nppvec[1] = npclrxyz*Math.sin(thdir)*Math.sin(phdir);
+ nppvec[2] = npclrxyz*Math.cos(thdir);
+ Hep3Vector npp = new BasicHep3Vector(nppvec);
+ double magnpp = Math.sqrt(nppvec[0]*nppvec[0]+nppvec[1]*nppvec[1]+nppvec[2]*nppvec[2]);
+ // get vector of cluster position
+ double[] npclvec = new double[3];
+ npclvec[0] = clpos[0];
+ npclvec[1] = clpos[1];
+ npclvec[2] = clpos[2];
+ Hep3Vector npcl = new BasicHep3Vector(npclvec);
+ double magnpcl = Math.sqrt(npclvec[0]*npclvec[0]+npclvec[1]*npclvec[1]+npclvec[2]*npclvec[2]);
+ // now get ILSP clusters to compare to
+ for (BasicCluster ilspc : ILSPclus)
+ {
+ // get this cluster's position
+ double[] ilclpos = ilspc.getPosition();
+// Hep3Vector ilspcl = new BasicHep3Vector(ilclpos);
+ // calculate Itheta and Iphi for the clus-clus direction
+ double delx = ilclpos[0]-npclvec[0];
+ double dely = ilclpos[1]-npclvec[1];
+ double delz = ilclpos[2]-npclvec[2];
+ double[] npilvec = new double[3];
+ npilvec[0] = delx;
+ npilvec[1] = dely;
+ npilvec[2] = delz;
+ Hep3Vector npil = new BasicHep3Vector(npilvec);
+ double magnpi = Math.sqrt(npilvec[0]*npilvec[0]+npilvec[1]*npilvec[1]+npilvec[2]*npilvec[2]);
+ double Iclph = Math.atan2(dely,delx);
+ if (Iclph<0) Iclph+=2*Math.PI;
+ double ccr = Math.sqrt(delx*delx+dely*dely);
+ double Iclth = Math.atan(ccr/delz);
+ if (Iclth<0) Iclth+=Math.PI;
+ // calculate distance in x,y between clusters
+ double dclxyz = Math.sqrt(delx*delx+dely*dely+delz*delz);
+ aida.cloud1D("Distance in xyz between NP and ILSP clusters").fill(dclxyz);
+ // calculate difference between clus pointing and clus-clus direction
+ double dilth = Math.abs(thdir-Iclth);
+ aida.cloud1D("Diff Iclth ILSPth").fill(dilth);
+ double dilph = Math.abs(phdir-Iclph);
+ if (dilph>Math.PI) dilph = 2*Math.PI-dilph;
+ aida.cloud1D("Diff Iclph ILSPph").fill(dilph);
+ aida.cloud2D("Diff IclILSP theta phi").fill(dilph,dilth);
+ double dilthph = Math.sqrt(dilth*dilth+dilph*dilph);
+ aida.cloud1D("Diff diclilthph").fill(dilthph);
+ // compare dot product of cluster pointing to npil direction
+ double cpdot = VecOp.dot(npp, npil);
+ double cthcp = Math.abs(cpdot/(magnpp*magnpi));
+ double ccdot = VecOp.dot(npp, npcl);
+ double cthcc = Math.abs(ccdot/(magnpp*magnpcl));
+ aida.cloud1D("CosDel np IL").fill(cthcc);
+// if (dclxyz<300. && dilth<1.5 && dilth<distclil)
+ if (cthcp>cthcc)
+ {
+ // this might be a fragment of a charged hadron
+ distclil = dilth;
+ iclilpt++;
+ }
+ } // end of ILSP cluster loop
+ if (iclilpt>0) NPILclus.add(npc);
+ }
+ event.put("NPILClusters",NPILclus);
+
+ // check purity of ILSP clusters - should be linked to charged particles
+ int nhev = 0;
+ int nchev = 0;
+ int nphev = 0;
+ int nnhev = 0;
+ for (BasicCluster icl : ILSPclus)
+ {
+ int numchhits = 0;
+ int numphhits = 0;
+ int numnhhits = 0;
+ int nclhits = icl.getCalorimeterHits().size();
+ nhev += icl.getCalorimeterHits().size();
+ for (CalorimeterHit chit : icl.getCalorimeterHits())
+ {
+ SimCalorimeterHit schit = (SimCalorimeterHit) chit;
+ double hitch = Math.abs(schit.getMCParticle(0).getCharge());
+ double hitm = schit.getMCParticle(0).getMass();
+ if (hitch>0) numchhits++;
+ if (hitch==0 && hitm==0) numphhits++;
+ if (hitch==0 && hitm>0) numnhhits++;
+ if (hitch>0) nchev++;
+ if (hitch==0 && hitm==0) nphev++;
+ if (hitch==0 && hitm>0) nnhev++;
+ }
+ double dnh = nclhits;
+ double dch = numchhits;
+ double dph = numphhits;
+ double dnuh = numnhhits;
+ aida.cloud1D("MC Purity for ILSP Clusters").fill(dch/dnh);
+ aida.cloud1D("Pho cont to ILSP").fill(dph/dnh);
+ aida.cloud1D("NH cont to ILSP").fill(dnuh/dnh);
+ }
+ double dnhev = nhev;
+ double dnchev = nchev;
+ double dnphev = nphev;
+ double dnnhev = nnhev;
+ aida.cloud1D("MC Purity for ILSP Clus per event").fill(dnchev/dnhev);
+ aida.cloud1D("Pho cont to ILSP per event").fill(dnphev/dnhev);
+ aida.cloud1D("NH cont to ILSP per event").fill(dnnhev/dnhev);
+
+ // now check the IP-pointing clusters - divided into the 2 classes
+ // first, the photon-like class
+ int nhev2 = 0;
+ int nchev2 = 0;
+ int nphev2 = 0;
+ int nnhev2 = 0;
+ for (BasicCluster icl : IPphclus)
+ {
+ int numchhits = 0;
+ int numphhits = 0;
+ int numnhhits = 0;
+ int nclhits = icl.getCalorimeterHits().size();
+ nhev2 += icl.getCalorimeterHits().size();
+ for (CalorimeterHit chit : icl.getCalorimeterHits())
+ {
+ SimCalorimeterHit schit = (SimCalorimeterHit) chit;
+ double hitch = Math.abs(schit.getMCParticle(0).getCharge());
+ double hitm = schit.getMCParticle(0).getMass();
+ if (hitch>0) numchhits++;
+ if (hitch==0 && hitm==0) numphhits++;
+ if (hitch==0 && hitm>0) numnhhits++;
+ if (hitch>0) nchev2++;
+ if (hitch==0 && hitm==0) nphev2++;
+ if (hitch==0 && hitm>0) nnhev2++;
+ }
+ double dnh = nclhits;
+ double dch = numchhits;
+ double dph = numphhits;
+ double dnuh = numnhhits;
+ aida.cloud1D("MC Pho Purity for IPPho Clusters").fill(dph/dnh);
+ aida.cloud1D("NH cont to IPPho Clusters").fill(dnuh/dnh);
+ aida.cloud1D("MC Pho+NH Purity for IPPho Clusters").fill((dph+dnuh)/dnh);
+ aida.cloud1D("ChH cont to IPPho CLusters").fill(dch/dnh);
+ }
+ double dnhev2 = nhev2;
+ double dnchev2 = nchev2;
+ double dnphev2 = nphev2;
+ double dnnhev2 = nnhev2;
+ aida.cloud1D("MC Pho Purity for IPPho Clus per event").fill(dnphev2/dnhev2);
+ aida.cloud1D("NH cont to IPPho per event").fill(dnnhev2/dnhev2);
+ aida.cloud1D("MC Pho+NH Purity for IPPho Clus per event").fill((dnphev2+dnnhev2)/dnhev2);
+ aida.cloud1D("ChH cont to IPPho per event").fill(dnchev2/dnhev2);
+
+ // now check purity of IPTrack Clusters
+ int nhev4 = 0;
+ int nchev4 = 0;
+ int nphev4 = 0;
+ int nnhev4 = 0;
+ for (BasicCluster icl : IPtrclus)
+ {
+ int numchhits = 0;
+ int numphhits = 0;
+ int numnhhits = 0;
+ int nclhits = icl.getCalorimeterHits().size();
+ nhev4 += icl.getCalorimeterHits().size();
+ for (CalorimeterHit chit : icl.getCalorimeterHits())
+ {
+ SimCalorimeterHit schit = (SimCalorimeterHit) chit;
+ double hitch = Math.abs(schit.getMCParticle(0).getCharge());
+ double hitm = schit.getMCParticle(0).getMass();
+ if (hitch>0) numchhits++;
+ if (hitch==0 && hitm==0) numphhits++;
+ if (hitch==0 && hitm>0) numnhhits++;
+ if (hitch>0) nchev4++;
+ if (hitch==0 && hitm==0) nphev4++;
+ if (hitch==0 && hitm>0) nnhev4++;
+ }
+ double dnh = nclhits;
+ double dch = numchhits;
+ double dph = numphhits;
+ double dnuh = numnhhits;
+ aida.cloud1D("Pho cont to IPTr Clusters").fill(dph/dnh);
+ aida.cloud1D("NH cont to IPTr Clusters").fill(dnuh/dnh);
+ aida.cloud1D("MC ChH Purity for IPTr CLusters").fill(dch/dnh);
+ }
+ double dnhev4 = nhev4;
+ double dnchev4 = nchev4;
+ double dnphev4 = nphev4;
+ double dnnhev4 = nnhev4;
+ aida.cloud1D("Pho cont to IPTr Clus per event").fill(dnphev4/dnhev4);
+ aida.cloud1D("NH cont to IPTr Clus per event").fill(dnnhev4/dnhev4);
+ aida.cloud1D("MC ChH Purity for IPTr per event").fill(dnchev4/dnhev4);
+
+ // now check remaining Non-Pointing clusters - are these neutral hadrons?
+ int nhev3 = 0;
+ int nchev3 = 0;
+ int nphev3 = 0;
+ int nnhev3 = 0;
+ for (BasicCluster icl : NPclus)
+ {
+ int numchhits = 0;
+ int numphhits = 0;
+ int numnhhits = 0;
+ int nclhits = icl.getCalorimeterHits().size();
+ nhev3 += icl.getCalorimeterHits().size();
+ for (CalorimeterHit chit : icl.getCalorimeterHits())
+ {
+ SimCalorimeterHit schit = (SimCalorimeterHit) chit;
+ double hitch = Math.abs(schit.getMCParticle(0).getCharge());
+ double hitm = schit.getMCParticle(0).getMass();
+ if (hitch>0) numchhits++;
+ if (hitch==0 && hitm==0) numphhits++;
+ if (hitch==0 && hitm>0) numnhhits++;
+ if (hitch>0) nchev3++;
+ if (hitch==0 && hitm==0) nphev3++;
+ if (hitch==0 && hitm>0) nnhev3++;
+ }
+ double dnh = nclhits;
+ double dch = numchhits;
+ double dph = numphhits;
+ double dnuh = numnhhits;
+ aida.cloud1D("MC NH Purity for NP Clusters").fill(dnuh/dnh);
+ aida.cloud1D("ChH cont to NP Clusters").fill(dch/dnh);
+ aida.cloud1D("Pho cont to NP Clusters").fill(dph/dnh);
+ }
+ double dnhev3 = nhev3;
+ double dnchev3 = nchev3;
+ double dnphev3 = nphev3;
+ double dnnhev3 = nnhev3;
+ aida.cloud1D("MC NH Purity for NP Clus per event").fill(dnnhev3/dnhev3);
+ aida.cloud1D("ChH cont to NP Clus per event").fill(dnchev3/dnhev3);
+ aida.cloud1D("Pho cont to NP Clus per event").fill(dnphev3/dnhev3);
+
+ }
+
+ public void setInputClusterName(String name)
+ {
+ _clustername = name;
+ }
+ public void setOutputClusterLists(String[] outclnames)
+ {
+ _oclnames = outclnames;
+ }
+}
+