hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -N TriggerTurnOnAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TriggerTurnOnAnalysis.java 2 Oct 2012 01:56:58 -0000 1.1
@@ -0,0 +1,275 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.hps.users.phansson;
+
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.geometry.Detector;
+import org.lcsim.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.lcsim.hps.evio.TriggerData;
+import org.lcsim.hps.recon.ecal.HPSEcalCluster;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author phansson
+ */
+public class TriggerTurnOnAnalysis extends Driver {
+ private boolean _DEBUG = false;
+ private String ecalClusterCollectionName = "EcalTriggerClusters";
+ private String triggerDataCollectionName = "TriggerStatus";
+ private String outputPlotFileName = "trigturnonanalysis.aida";
+
+ private int totalEvents = 0;
+
+ private List<MCParticle> electrons = null;
+ private List<MCParticle> positrons = null;
+
+ private IHistogram1D hCountTrig;
+ //private IHistogram1D hTrigEffThetaySmallest;
+ private IHistogram1D hTrigThetaySmallestEcut;
+ private IHistogram1D hThetaySmallest;
+ private IHistogram1D hThetaySmallestEcut;
+ private IHistogram2D hThetaySmallestvsE;
+ private IHistogram2D hTrigThetaySmallestvsE;
+ private IHistogram2D hThetayvsThetay;
+ private IHistogram2D hTrigThetayvsThetay;
+ private IHistogram2D hele1vsele2;
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IAnalysisFactory af = aida.analysisFactory();
+ IHistogramFactory hf = aida.histogramFactory();
+ IPlotter plotter_count;
+ IPlotter plotter_count_1;
+ IPlotter plotter_count_2;
+ IPlotter plotter_count_3;
+ IPlotter plotter_count_4;
+
+ @Override
+ public void detectorChanged(Detector detector) {
+
+
+ makePlots();
+
+
+
+ }
+
+ @Override
+ public void process(EventHeader event) {
+
+ if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Process event " + event.getEventNumber());
+
+ List<MCParticle> mcparticles = null;
+ if(event.hasCollection(MCParticle.class)) {
+ mcparticles = event.get(MCParticle.class).get(0);
+ if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Number of MC particles = " + mcparticles.size());
+
+ List<MCParticle> fsParticles = HPSMCParticlePlotsDriver.makeGenFSParticleList(mcparticles);
+ if(_DEBUG) System.out.println(this.getClass().getSimpleName()+": Number of FS MC particles = " + fsParticles.size());
+ electrons = this.getTriggerCandidates(11, -1, -1, fsParticles);
+ positrons = this.getTriggerCandidates(-11, -1, -1, fsParticles);
+ if(_DEBUG) {
+ System.out.println(this.getClass().getSimpleName()+": " + electrons.size() + " electrons");
+ System.out.println(this.getClass().getSimpleName()+": " + positrons.size() + " electrons");
+ }
+
+
+ }
+ else {
+ System.out.println(this.getClass().getSimpleName() + ": no MC particles in this event");
+ }
+
+
+
+
+ List<TriggerData> trigger_data = null;
+
+ if(event.hasCollection(TriggerData.class, triggerDataCollectionName)) {
+ trigger_data = event.get(TriggerData.class, triggerDataCollectionName);
+ if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": event " + event.getRunNumber() + " has trigger data");
+ }
+ else {
+ if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": event " + event.getRunNumber() + " has no trigger data");
+ }
+
+ List<HPSEcalCluster> clusters = null;
+ if( event.hasCollection(HPSEcalCluster.class, ecalClusterCollectionName) ) {
+ clusters = event.get(HPSEcalCluster.class, ecalClusterCollectionName);
+ if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": event " + event.getRunNumber() + " has " + clusters.size() + " ecal clusters");
+ }
+ else {
+ if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": event " + event.getRunNumber() + " has no ecal clusters");
+ }
+
+ //fill counting histogram
+ hCountTrig.fill( trigger_data == null ? 0 : 1);
+
+ //find smallest angle of the proposed pair that fired the trigger
+
+ // Use electron/positron with highest E
+ MCParticle electron = null;
+ for(MCParticle e : electrons) {
+ if (electron==null) {
+ electron = e;
+ } else {
+ if(e.getEnergy() > electron.getEnergy()) {
+ electron = e;
+ }
+ }
+ }
+ if(electrons.size()>=2) {
+ hele1vsele2.fill(electrons.get(0).getEnergy(),electrons.get(1).getEnergy());
+ }
+ MCParticle positron = null;
+ for(MCParticle e : positrons) {
+ if (positron==null) {
+ positron = e;
+ } else {
+ if(e.getEnergy() > positron.getEnergy()) {
+ positron = e;
+ }
+ }
+ }
+ if(electron!=null && positron!=null) {
+ double electron_thetay = Math.abs(Math.atan(electron.getMomentum().y()/electron.getMomentum().z()));
+ double positron_thetay = Math.abs(Math.atan(positron.getMomentum().y()/positron.getMomentum().z()));
+ double thetay_smallest = electron_thetay < positron_thetay ? electron_thetay : positron_thetay;
+ double E_thetay_smallest = electron_thetay < positron_thetay ? electron.getEnergy() : positron.getEnergy();
+ double thetay_larger = electron_thetay < positron_thetay ? positron_thetay : electron_thetay;
+ hThetaySmallest.fill( thetay_smallest );
+ hThetaySmallestvsE.fill(thetay_smallest, E_thetay_smallest);
+ hThetayvsThetay.fill(thetay_smallest, thetay_larger );
+ if(E_thetay_smallest>0.0) {
+ hThetaySmallestEcut.fill( thetay_smallest );
+ if(trigger_data!=null) {
+ hTrigThetaySmallestvsE.fill(thetay_smallest, E_thetay_smallest);
+ hTrigThetayvsThetay.fill(thetay_smallest, thetay_larger );
+ hTrigThetaySmallestEcut.fill( thetay_smallest );
+ }
+ }
+ }
+ if(totalEvents % 500 == 0) updatePlots();
+ totalEvents++;
+
+ }
+
+ MCParticle ele_cand = null;
+ MCParticle pos_cand = null;
+
+ List<MCParticle> getTriggerCandidates(int pdgid, double E_cut, double thetay_cut, List<MCParticle> fsParticles) {
+ List<MCParticle> particles = new ArrayList<MCParticle>();
+ for(MCParticle particle : fsParticles) {
+ int pdgID = particle.getPDGID();
+ if(pdgid!=pdgID) continue;
+ if(E_cut>0 && particle.getEnergy()<E_cut) continue;
+ Hep3Vector p = particle.getMomentum();
+ double thetay = Math.abs(Math.atan(p.y()/p.z()));
+ if(thetay_cut>0 && thetay<thetay_cut) continue;
+ particles.add(particle);
+ }
+ return particles;
+ }
+
+ @Override
+ public void endOfData() {
+ updatePlots();
+ if (!"".equals(outputPlotFileName))
+ try {
+ aida.saveAs(outputPlotFileName);
+ } catch (IOException ex) {
+ Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
+ }
+ }
+
+
+ private void makePlots() {
+ hCountTrig = aida.histogram1D("hCountTrig",2,0,2);
+ //hTrigThetaySmallestE = aida.histogram1D("hTrigThetaySmallest",50,0,0.1);
+ hTrigThetaySmallestEcut = aida.histogram1D("hTrigThetaySmallestEcut",50,0,0.1);
+ hThetaySmallest = aida.histogram1D("hThetaySmallest",50,0,0.1);
+ hThetaySmallestEcut = aida.histogram1D("hThetaySmallestEcut",50,0,0.1);
+ hThetaySmallestvsE = aida.histogram2D("hThetaySmallestvsE",25,0,0.1,20,0,2.);
+ hTrigThetaySmallestvsE = aida.histogram2D("hTrigThetaySmallestvsE",25,0,0.1,20,0,2.);
+ hThetayvsThetay = aida.histogram2D("hThetayvsThetay",50,0,0.1,50,0,0.2);
+ hTrigThetayvsThetay = aida.histogram2D("hTrigThetayvsThetay",50,0,0.1,50,0,0.2);
+
+ hele1vsele2 = aida.histogram2D("hele1vsele2", 50,0,2,50,0,2);
+
+ plotter_count = af.createPlotterFactory().create();
+ plotter_count.createRegions(1,1);
+ plotter_count.setTitle("Trigger Count");
+ //plotter_count.style().statisticsBoxStyle().setVisible(true);
+ plotter_count.region(0).plot(hCountTrig);
+ plotter_count.show();
+
+ plotter_count_1 = af.createPlotterFactory().create();
+ plotter_count_1.createRegions(2,2);
+ plotter_count_1.setTitle("Acceptance vs Thetay");
+ //plotter_count.style().statisticsBoxStyle().setVisible(true);
+ plotter_count_1.region(0).plot(hThetaySmallest);
+ plotter_count_1.region(1).plot(hThetaySmallestEcut);
+ plotter_count_1.region(2).plot(hTrigThetaySmallestEcut);
+ plotter_count_1.region(3).style().statisticsBoxStyle().setVisible(false);
+ plotter_count_1.show();
+
+ plotter_count_2 = af.createPlotterFactory().create();
+ plotter_count_2.createRegions(2,2);
+ plotter_count_2.setTitle("Trigger Count");
+ //plotter_count.style().statisticsBoxStyle().setVisible(true);
+ plotter_count_2.style().setParameter("hist2DStyle", "colorMap");
+ plotter_count_2.style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ plotter_count_2.region(0).plot(hThetaySmallestvsE);
+ plotter_count_2.region(1).plot(hTrigThetaySmallestvsE);
+ plotter_count_2.region(2).style().statisticsBoxStyle().setVisible(false);
+ plotter_count_2.show();
+
+ IPlotter plotter_count_3 = af.createPlotterFactory().create();
+ plotter_count_3.createRegions(1,2);
+ plotter_count_3.setTitle("Trigger Count");
+ //plotter_count.style().statisticsBoxStyle().setVisible(true);
+ plotter_count_3.style().setParameter("hist2DStyle", "colorMap");
+ plotter_count_3.style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ plotter_count_3.region(0).plot(hThetayvsThetay);
+ plotter_count_3.region(1).plot(hTrigThetayvsThetay);
+ plotter_count_3.show();
+
+
+ plotter_count_4 = af.createPlotterFactory().create();
+ plotter_count_4.createRegions(1,1);
+ plotter_count_4.setTitle("e- vs e-");
+ //plotter_count.style().statisticsBoxStyle().setVisible(true);
+ plotter_count_4.style().setParameter("hist2DStyle", "colorMap");
+ plotter_count_4.style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ plotter_count_4.region(0).style().statisticsBoxStyle().setVisible(false);
+ plotter_count_4.region(0).plot(hele1vsele2);
+ plotter_count_4.show();
+
+ }
+
+ private void updatePlots() {
+
+ //hTrigEffThetaySmallest = aida.histogram1D("hTrigEffThetaySmallest",50,0,0.1);
+ IHistogram1D hTrigEffThetaySmallest = hf.divide("hTrigEffThetaySmallest", hTrigThetaySmallestEcut, hThetaySmallestEcut);
+ plotter_count_1.region(3).clear();
+ plotter_count_1.region(3).plot(hTrigEffThetaySmallest);
+ plotter_count_1.region(3).style().statisticsBoxStyle().setVisible(false);
+ IHistogram2D hTrigEffThetaySmallestvsE = hf.divide("hTrigEffThetaySmallestvsE", hTrigThetaySmallestvsE,hThetaySmallestvsE);
+ plotter_count_2.region(2).clear();
+ plotter_count_2.region(2).plot(hTrigEffThetaySmallestvsE);
+ plotter_count_2.region(2).style().statisticsBoxStyle().setVisible(false);
+
+ }
+
+}