lcsim/src/org/lcsim/recon/pfa/filter
diff -N EventFilter.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EventFilter.java 20 Jul 2009 23:59:27 -0000 1.1
@@ -0,0 +1,238 @@
+package org.lcsim.recon.pfa.filter;
+
+import java.util.List;
+import hep.aida.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.*;
+
+/**
+ *
+ * @author rfc
+ *
+ * Accept or reject events based on presence or absence of a particular
+ * particle in the event. Intended for use in studies of leading-particle biasing
+ * effects in calorimeters
+ *
+ * $Id: EventFilter.java,v 1.1 2009/07/20 23:59:27 rfc Exp $
+ *
+ * $Log: EventFilter.java,v $
+ * Revision 1.1 2009/07/20 23:59:27 rfc
+ * Copy org.lcsim.contrib.mit.EventFilter to org.lcsim.recon.pfa.filter.EventFilter to work around present limitations in running from xml files (lcsim-contrib classes cannot be called due to build/dependency structure).
+ *
+ * Revision 1.1 2009/07/18 02:01:47 rfc
+ * Added driver lcsim.contrib.mit.EventFilter.
+ *
+ *
+ */
+public class EventFilter extends Driver {
+
+ String version="0.10";
+
+ /**
+ * AIDA histogram declarations
+ */
+ private static AIDA aida = AIDA.defaultInstance();
+ ICloud1D h1;
+
+ protected int PdgId;
+ protected double EnergyCut;
+ protected int RequireOrRejectFlag;
+ protected int EventCount, AcceptedEventCount, RejectedEventCount;
+ protected int nPrinted, nPrintMax;
+
+ public void setPdgId( int pdgid ) { PdgId = pdgid; }
+ public void setEnergyCut( double cut ) { EnergyCut = cut; }
+ public void setRequireOrRejectFlag( int flag ) { RequireOrRejectFlag = flag; }
+
+ /**
+ * Various constructors
+ *
+ */
+ public EventFilter() {
+ this( 111, 30.0, 1); // pi0, 30 GeV, reject
+ };
+
+ public EventFilter( int pdgid ) {
+ this( pdgid, 30.0, 1); // pdgid, 30 GeV cutoff, reject
+ }
+
+ public EventFilter( int pdgid, double cut ) {
+ this( pdgid, cut, 1 ); // pdgid, cut, reject
+ }
+
+ public EventFilter( int pdgid, double cut, int flag ) {
+ PdgId = pdgid;
+ EnergyCut = cut;
+ RequireOrRejectFlag = flag;
+ EventCount = 0;
+ AcceptedEventCount = 0;
+ RejectedEventCount = 0;
+ nPrinted = 0;
+ nPrintMax = 100;
+ }
+
+ /**
+ * Startup initialization
+ */
+ protected void startOfData() {
+ System.out.println( "EventFilter version " + version +
+ " initialization -------------------");
+ printConfig();
+ /*
+ * Initialize histograms
+ */
+
+
+ System.out.println( "...end EventFilter initialization -------------------");
+ }
+
+ /**
+ * Event loop
+ * @param event The event.
+ */
+ @Override
+ protected void process(EventHeader event)
+ {
+ EventCount++;
+ /**
+ * Look for leading particles of type PdgId with energy > EnergyCut.
+ * If RequireOrRejectFlag == 1, reject events with this particle; if
+ * RequireOrRejectFlag == -1, reject events without this particle; if
+ * RequireOrRejectFlag has any other value, do nothing (which causes
+ * the event to be accepted, but no accept/reject counters to be
+ * incremented.
+ */
+ //if (event.getEventNumber() % 2 != 0)
+ int reject;
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+
+ double Zmass = 0.;
+ double ZE = 0.;
+ double cost = 0.;
+
+ if( RequireOrRejectFlag == 1 ) {
+ reject = 0;
+ } else {
+ reject = 0;
+ }
+
+ aida.cloud1D("evtfilt_all_nTracks").fill(mcl.size());
+ for(MCParticle p:mcl)
+ {
+ aida.cloud1D("evtfilt_all_energy").fill(p.getEnergy());
+ aida.cloud1D("evtfilt_all_cosTheta").fill(VecOp.cosTheta(p.getMomentum()));
+ aida.cloud1D("evtfilt_all_phi").fill(VecOp.phi(p.getMomentum()));
+
+ int id = Math.abs(p.getPDGID());
+ //System.out.println( "___id = " + id );
+
+ /*
+ * Apply event filtering cut
+ */
+ if( id == PdgId ) {
+ // A few histograms for debugging
+ aida.cloud1D("evtfilt_pdgid"+PdgId+"_energy").fill(p.getEnergy());
+ aida.cloud1D("evtfilt_pdgid"+PdgId+"_all_cosTheta").fill(VecOp.cosTheta(p.getMomentum()));
+ aida.cloud1D("evtfilt_pdgid"+PdgId+"_all_phi").fill(VecOp.phi(p.getMomentum()));
+
+ /*if( RequireOrRejectFlag == -1 ) {
+ // Require one or more pdgid's in the event with energy > cut
+ if( p.getEnergy() >= EnergyCut ) {
+ reject = 0;
+ }
+ } else if ( RequireOrRejectFlag == +1 ) {*/
+ // Require zero particles in the event with energy > cut
+ if( p.getEnergy() >= EnergyCut ) {
+ reject = 1;
+ if( RequireOrRejectFlag == 1 ) {
+ if( nPrinted < nPrintMax ) {
+ nPrinted++;
+ System.out.println( "EventFilter rejecting event #"
+ + EventCount +
+ " containing particle with id = " + id +
+ " with energy " + p.getEnergy() + " GeV");
+ }
+ }
+ break; // Out of loop over p:mcl
+ }
+ /*}*/
+ }
+
+ /*if( (id == 1 )||(id == 2 )||(id == 3 ) )
+ {
+ if(p.getParents().get(0).getPDGID() == 23 )
+ {
+ Zmass = p.getParents().get(0).getMass();
+ ZE = p.getParents().get(0).getEnergy();
+ Hep3Vector pp = p.getMomentum();
+ cost = Math.abs(pp.z()/pp.magnitude());
+ aida.cloud1D("evtfilt_Z_mass").fill(Zmass);
+ aida.cloud1D("evtfilt_Z_energy").fill(ZE);
+ aida.cloud1D("evtfile_Z_ptot").fill(pp.magnitude());
+ aida.cloud1D("evtfile_Z_costheta").fill(cost);
+ //System.out.print( "______Zmass, ZE, pp, cost = ");
+ //System.out.println( Zmass + ", " + ZE + ", " + pp + ", " + cost);
+ }
+ }*/
+ }
+
+ if( RequireOrRejectFlag == -1 ) {
+ // Flip reject flag
+ if( reject == 1 ) { reject = 0; } else { reject = 1; }
+ }
+
+ if( ( reject == 1 ) && ( RequireOrRejectFlag == -1 ) ) {
+ if( nPrinted < nPrintMax ) {
+ nPrinted++;
+ System.out.println( "EventFilter rejecting event #" + EventCount +
+ " containing no particle with id = " + PdgId +
+ " with energy > " + EnergyCut + " GeV");
+ }
+ }
+ if( nPrinted == nPrintMax ) {
+ nPrinted++;
+ System.out.println( "EventFilter printout terminated: max print count" +
+ " (" + nPrintMax + " lines) reached");
+ }
+
+ if( reject == 1 )
+ {
+ RejectedEventCount++;
+ throw new Driver.NextEventException();
+ } else {
+ AcceptedEventCount++;
+ // This takes care that the child Drivers are loaded and processed.
+ super.process(event);
+ }
+ }
+
+ /**
+ * Called by the framework when event processing is suspended, usually at
+ * the end of data processing
+ */
+ protected void suspend()
+ {
+ System.out.println( "EventFilter version " + version +
+ " summary -----------------------------" );
+ printConfig();
+ System.out.println( "...EventCount = " + EventCount );
+ System.out.println( "...RejectedEventCount = " + RejectedEventCount );
+ System.out.println( "...AcceptedEventCount = " + AcceptedEventCount );
+ System.out.println( "...end EventFilter summary -----------------------------" );
+ }
+
+ protected void printConfig() {
+ System.out.println( "...PdgId = " + PdgId );
+ System.out.println( "...EnergyCut = " + EnergyCut + " GeV");
+ System.out.println( "...RequireOrRejectFlag = " + RequireOrRejectFlag );
+ if( ( RequireOrRejectFlag != 1 ) && ( RequireOrRejectFlag != -1 ) ) {
+ System.out.println( "...EventFilter warning: RequireOrRejectFlag" +
+ "should be +1 or -1" );
+ }
+ System.out.println( "...nPrintMax = " + nPrintMax + " lines");
+ }
+
+}
