lcsim-contrib/src/main/java/org/lcsim/contrib/twhite
diff -N MuonAnalysisTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonAnalysisTest.java 12 Aug 2009 08:52:02 -0000 1.1
@@ -0,0 +1,35 @@
+package main.java.org.lcsim.contrib.twhite;
+/**
+ * @version $Id: MuonAnalysisTest.java,v 1.1 2009/08/12 08:52:02 twhite Exp $
+ * @author T White
+ */
+
+import junit.framework.TestCase;
+
+
+public class MuonAnalysisTest extends TestCase{
+
+ MuonAnalysis muonanalysis2;
+
+ protected void setUp() throws Exception {
+
+ muonanalysis2 = new MuonAnalysis();
+
+ }
+
+ protected void tearDown() throws Exception {
+
+ }
+
+ public void testgettype(){
+ int returnedtype;
+ int array[] = {13,13,-13,13,12,9912,-14,9912,-16,9912,2212,211,-321,211,211,211,2112,130,10000000,9999};
+ for(int i=0;i<array.length;i=i+2){
+ returnedtype = muonanalysis2.gettype(array[i]);
+ assertEquals(array[i+1], returnedtype);
+ }
+
+ }
+
+
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/twhite
diff -N MuonAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonAnalysis.java 12 Aug 2009 08:52:02 -0000 1.1
@@ -0,0 +1,340 @@
+package main.java.org.lcsim.contrib.twhite;
+/**
+ * @version $Id: MuonAnalysis.java,v 1.1 2009/08/12 08:52:02 twhite Exp $
+ * @author T White
+ */
+
+import hep.physics.particle.Particle;
+import hep.physics.vec.VecOp;
+
+import java.io.File;
+import java.io.IOException;
+import java.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.lcio.LCIOWriter;
+
+public class MuonAnalysis extends Driver {
+
+ //Size of histogram for number of particles in the event MC and RC
+ int numberparthist[] = {150,0,150};
+ //Size of histogram for number of missing particles in the event
+ int numberpartmissinghist[] = {70,-20,50};
+ //Size of histogram for energy of MC particles
+ int energyhist[] = {500,0,500};
+ //Size of histogram for energy of missing particles
+ int energymissinghist[] = {75,0,75};
+ //Size of cos(angle) histogram
+ int anglehist[] = {200,-1,1};
+
+ String reconname = "ReconstructedParticles"; // Collection name from LCIO file
+ String mcparticlename = "MCParticle"; // Collection name from LCIO file
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ int ievt = 0; // Keeps track of event number
+
+ int beampipe; //Set to one if muon travels down beampipe
+
+ int numparticlesmc; //Number of total mc particles
+ int numparticlesrc; //Number of total rc particles
+
+ Object[] allparticlekeys; //Stores particle ids involved in the event
+
+ // Hash maps to store the number of each particle type
+ HashMap<Integer, Integer> particlenumbersmc;
+ HashMap<Integer, Integer> particlenumbersrc;
+
+ // Ensures collections are defined as global
+ List<MCParticle> mcparticles;
+ List<ReconstructedParticle> rcparticles;
+ // Lists to store the muons involved in the event
+ List<MCParticle> MCMuons;
+ List<ReconstructedParticle> RCMuons;
+
+ // Ensures writer is public
+ LCIOWriter writer;
+
+ public MuonAnalysis() {
+ aida = AIDA.defaultInstance();
+
+ // Opens a new writer
+ try {
+ File output = new File(System.getProperty("user.home"), "output.slcio");
+ writer = new LCIOWriter(output);
+ } catch (IOException x) {
+ throw new RuntimeException("Error writing LCIO file", x);
+ }
+ }
+
+ protected void process(EventHeader event) {
+
+ // Stores the muons in the event
+ MCMuons = new ArrayList<MCParticle>();
+ RCMuons = new ArrayList<ReconstructedParticle>();
+ // Retrieves the data from the LCIO file
+ mcparticles = event.get(MCParticle.class, mcparticlename);
+ rcparticles = event.get(ReconstructedParticle.class, reconname);
+
+ // Counts the number of particles and fills the HashMaps - also outputs data on MC particles for comparison
+ numparticlesmc = countmc();
+ numparticlesrc = countrc();
+
+ //Compile list of particles involved in event, either rc or mc
+ //Add any particles that you will be analysis to ensure data is output on EVERY event
+ //Otherwise only events that contain that particle in either the mc or rc particle lists will be used
+ //This can lead to a problem with analysis the average number of particles created per event
+ Set<Integer> ParticleList = new HashSet<Integer>();
+ ParticleList.addAll(particlenumbersmc.keySet());
+ ParticleList.addAll(particlenumbersrc.keySet());
+ ParticleList.add(130);
+ ParticleList.add(13);
+ ParticleList.add(211);
+ ParticleList.add(9912);
+ ParticleList.add(22);
+ ParticleList.add(11);
+ allparticlekeys = ParticleList.toArray();
+
+ //Fills empty hashmap values with 0 to prevent null pointer error
+ for(Object id : allparticlekeys){
+ if (particlenumbersmc.get(id.hashCode()) == null) {
+ particlenumbersmc.put(id.hashCode(), 0);
+ }
+ if (particlenumbersrc.get(id.hashCode()) == null) {
+ particlenumbersrc.put(id.hashCode(), 0);
+ }
+ }
+
+ //Outputs number of each particle
+ for (int i = 0; i < allparticlekeys.length; i++) {
+ aida.histogram1D("MC Data/" + allparticlekeys[i].toString() + "/Numbers",numberparthist[0],numberparthist[1],numberparthist[2]).fill(particlenumbersmc.get(allparticlekeys[i]));
+ aida.histogram1D("RC Data/" + allparticlekeys[i].toString() + "/Numbers",numberparthist[0],numberparthist[1],numberparthist[2]).fill(particlenumbersrc.get(allparticlekeys[i]));
+ }
+
+ //Output data for all events
+ outputhistograms("Events ( All )");
+
+ //Output data for events with specific number of muons
+ outputhistograms("Events with " + MCMuons.size() + " muons");
+
+ //Output cut data - not all muons reconstructed
+ if(RCMuons.size()<MCMuons.size()){
+ outputhistograms("Events (too few muons)");
+
+// Output event to LCIO file
+// try {
+// writer.write(event);
+// } catch (IOException x) {
+// throw new RuntimeException("Error writing LCIO file", x);
+// }
+
+ }
+
+ //Output cut data - excess particles reconstructed
+ if(numparticlesrc>numparticlesmc){
+ outputhistograms("Events (excess particles)");
+ }
+
+ //Looking for correlations between missing particles and gained particles
+ for(int i=0;i<allparticlekeys.length;i++){
+ for(int j=0; j<i;j++){
+ aida.cloud2D("Missing Particle Correlations/Particle " + allparticlekeys[i].toString() + " vs Particle " + allparticlekeys[j].toString()).fill(particlenumbersmc.get(allparticlekeys[i]) - particlenumbersrc.get(allparticlekeys[i]), particlenumbersmc.get(allparticlekeys[j]) - particlenumbersrc.get(allparticlekeys[j]));
+ }
+ }
+
+ //Event counter
+ if (ievt % 100 == 0) {
+ System.out.println("Event = " + ievt);
+ }
+ ievt++;
+
+ }
+
+ protected void outputhistograms(String cuttype){
+ //Output data on all involved particles for all events AND split by beampipe
+ for (int i = 0; i < allparticlekeys.length; i++) {
+ aida.histogram1D(cuttype + "/" + allparticlekeys[i].toString() + "/Number Part RC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(particlenumbersrc.get(allparticlekeys[i]));
+ aida.histogram1D(cuttype + "/" + allparticlekeys[i].toString() + "/beampipe = " + beampipe + "/Number Part RC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(particlenumbersrc.get(allparticlekeys[i]));
+ aida.histogram1D(cuttype + "/" + allparticlekeys[i].toString() + "/Number Part MC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(particlenumbersmc.get(allparticlekeys[i]));
+ aida.histogram1D(cuttype + "/" + allparticlekeys[i].toString() + "/beampipe = " + beampipe + "/Number Part MC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(particlenumbersmc.get(allparticlekeys[i]));
+ aida.histogram1D(cuttype + "/" + allparticlekeys[i].toString() + "/Missing Part",numberpartmissinghist[0],numberpartmissinghist[1],numberpartmissinghist[2]).fill(particlenumbersmc.get(allparticlekeys[i]) - particlenumbersrc.get(allparticlekeys[i]));
+ aida.histogram1D(cuttype + "/" + allparticlekeys[i].toString() + "/beampipe = " + beampipe + "/Missing Part",numberpartmissinghist[0],numberpartmissinghist[1],numberpartmissinghist[2]).fill(particlenumbersmc.get(allparticlekeys[i]) - particlenumbersrc.get(allparticlekeys[i]));
+
+ //Produces data on specific not found particles rather than just the number of particles used up to this point
+ particlehistograms(cuttype,allparticlekeys[i].hashCode());
+ }
+
+ //Data on all particles not split by type as above
+ aida.histogram1D(cuttype + "/all/beampipe = " + beampipe + "/Number Part RC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(numparticlesrc);
+ aida.histogram1D(cuttype + "/all/beampipe = " + beampipe + "/Number Part MC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(numparticlesmc);
+ aida.histogram1D(cuttype + "/all/beampipe = " + beampipe + "/Missing Part",numberpartmissinghist[0],numberpartmissinghist[1],numberpartmissinghist[2]).fill(numparticlesmc - numparticlesrc);
+ aida.histogram1D(cuttype + "/all/Number Part RC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(numparticlesrc);
+ aida.histogram1D(cuttype + "/all/Number Part MC",numberparthist[0],numberparthist[1],numberparthist[2]).fill(numparticlesmc);
+ aida.histogram1D(cuttype + "/all/Missing Part",numberpartmissinghist[0],numberpartmissinghist[1],numberpartmissinghist[2]).fill(numparticlesmc - numparticlesrc);
+ }
+
+ protected void particlehistograms(String cuttype, int type) {
+ //Loop over MC particles and output details about particles that are MISSING from the reconstruction
+ for (MCParticle mcparticle : mcparticles) {
+ if ((mcparticle.getGeneratorStatus() == Particle.FINAL_STATE) && (mcparticle.getPDGID() == type)) {
+ //Flag for particle found in reconstruction
+ int found = 0;
+ // Check if particle can be found in corresponding bin
+ for (ReconstructedParticle rcparticle : rcparticles) {
+ if ((gettype(rcparticle.getParticleIDUsed().getPDG()) == gettype(mcparticle.getPDGID())) && (mcparticle.getCharge() == rcparticle.getCharge()) && ((VecOp.dot(mcparticle.getMomentum(), rcparticle.getMomentum()) / (mcparticle.getMomentum().magnitude() * rcparticle.getMomentum().magnitude())) > 0.99) && (Math.abs(rcparticle.getEnergy() - mcparticle.getEnergy()) < 5)) {
+ found = 1;
+ }
+ }
+ // If particle cannot be found plot data on mcparticle that is missing
+ if (found == 0) {
+ //Data split by type AND all
+ aida.histogram1D(cuttype + "/" + type + "/Missing/Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(mcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/" + type + "/Missing/Energy",energymissinghist[0],energymissinghist[1],energymissinghist[2]).fill(mcparticle.getEnergy());
+ aida.histogram1D(cuttype + "/all/Missing/Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(mcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/all/Missing/Energy",energymissinghist[0],energymissinghist[1],energymissinghist[2]).fill(mcparticle.getEnergy());
+ if (Math.abs(VecOp.cosTheta(mcparticle.getMomentum())) < 0.99) {
+ //Beampipe data split by type AND all
+ aida.histogram1D(cuttype + "/" + type + "/Missing/Angle Not Beampipe",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(mcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/" + type + "/Missing/Energy Not Beampipe",energymissinghist[0],energymissinghist[1],energymissinghist[2]).fill(mcparticle.getEnergy());
+ aida.histogram1D(cuttype + "/all/Missing/Angle Not Beampipe",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(mcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/all/Missing/Energy Not Beampipe",energymissinghist[0],energymissinghist[1],energymissinghist[2]).fill(mcparticle.getEnergy());
+ }
+ }
+ }
+ }
+ //Loops over the RC particles and output details about particles that are CREATED in the reconstruction
+ for (ReconstructedParticle rcparticle : rcparticles) {
+ if ((rcparticle.getParticleIDUsed().getPDG() == type)) {
+ int found = 0;
+ //Search for corresponding particle in corresponding bin
+ for (MCParticle mcparticle : mcparticles) {
+ if ((gettype(rcparticle.getParticleIDUsed().getPDG()) == gettype(mcparticle.getPDGID())) && (mcparticle.getCharge() == rcparticle.getCharge())
+ && ((VecOp.dot(mcparticle.getMomentum(), rcparticle.getMomentum()) / (mcparticle.getMomentum().magnitude() * rcparticle.getMomentum().magnitude())) > 0.99) && (Math.abs(rcparticle.getEnergy() - mcparticle.getEnergy()) < 5)) {
+ found = 1;
+ }
+ }
+ // If particle cannot be found plot data on mcparticle
+ if (found == 0) {
+ aida.histogram1D(cuttype + "/" + type + "/Created/Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(rcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/" + type + "/Created/Energy",energyhist[0],energyhist[1],energyhist[2]).fill(rcparticle.getEnergy());
+ aida.histogram1D(cuttype + "/all/Created/Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(rcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/all/Created/Energy",energyhist[0],energyhist[1],energyhist[2]).fill(rcparticle.getEnergy());
+ //Loops over muons to compare angle to created particles
+ for(MCParticle mcmuon : MCMuons){
+ aida.histogram1D(cuttype + "/" + type + "/Created/Angle to muons",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.dot(mcmuon.getMomentum(),rcparticle.getMomentum())/(mcmuon.getMomentum().magnitude()*rcparticle.getMomentum().magnitude()));
+ }
+ //Loops over clusters associated with created particle and plot number of hits against energy (perhaps need to sum up hits - or check number of clusters is 1)
+ for(Cluster cluster : rcparticle.getClusters()){
+ aida.cloud2D(cuttype + "/" + type + "/Created/Cal Hits vs energy").fill(cluster.getCalorimeterHits().size(),rcparticle.getEnergy());
+ aida.cloud2D(cuttype + "/all/Created/Cal Hits vs energy").fill(cluster.getCalorimeterHits().size(),rcparticle.getEnergy());
+ }
+ //Output same again without the beampipe events
+ if (Math.abs(VecOp.cosTheta(rcparticle.getMomentum())) < 0.99) {
+ aida.histogram1D(cuttype + "/" + type + "/Created/Angle Not Beampipe",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(rcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/" + type + "/Created/Energy Not Beampipe",energyhist[0],energyhist[1],energyhist[2]).fill(rcparticle.getEnergy());
+ aida.histogram1D(cuttype + "/all/Created/Angle Not Beampipe",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(rcparticle.getMomentum()));
+ aida.histogram1D(cuttype + "/all/Created/Energy Not Beampipe",energyhist[0],energyhist[1],energyhist[2]).fill(rcparticle.getEnergy());
+ for(Cluster cluster : rcparticle.getClusters()){
+ aida.cloud2D(cuttype + "/" + type + "/Created/Cal Hits vs energy Not Beampipe").fill(cluster.getCalorimeterHits().size(),rcparticle.getEnergy());
+ }
+ }
+ }
+ }
+ }
+ }
+
+ protected int gettype(int originaltype) {
+ // Used to change the types of particle to basic types (also removed the sign)
+ int type;
+ switch (Math.abs(originaltype)) {
+ case 12:
+ case 14:
+ case 16:
+ type = 9912; // Neutrinos collected together and called 9912
+ break;
+ case 2112:
+ type = 130; // Neutrons called a K-longs
+ break;
+ case 2212:
+ case 321:
+ type = 211; // Charged particles called pions
+ break;
+ default:
+ type = Math.abs(originaltype);
+ break;
+ }
+ return type;
+ }
+
+ protected int countmc() {
+ // Counts the number of mc particles and adds data to a histogram and particlenumbersmc
+ beampipe = 0;
+ int i = 0;
+ double energy = 0.0;
+ particlenumbersmc = new HashMap<Integer, Integer>();
+ for (MCParticle mcparticle : mcparticles) {
+ if (mcparticle.getGeneratorStatus() == Particle.FINAL_STATE) {
+ if (particlenumbersmc.get(gettype(mcparticle.getPDGID())) == null) {
+ particlenumbersmc.put(gettype(mcparticle.getPDGID()), 1);
+ } else {
+ particlenumbersmc.put(gettype(mcparticle.getPDGID()), particlenumbersmc.get(gettype(mcparticle.getPDGID())) + 1);
+ }
+ i++;
+ energy+=mcparticle.getEnergy();
+ //Outputs data of MC particles for comparison
+ aida.histogram1D("MC Data/" + gettype(mcparticle.getPDGID()) + "/Energy",energyhist[0],energyhist[1],energyhist[2]).fill(mcparticle.getEnergy());
+ aida.histogram1D("MC Data/" + gettype(mcparticle.getPDGID()) + "/Cos Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(mcparticle.getMomentum()));
+ aida.histogram1D("MC Data/all/Energy",energyhist[0],energyhist[1],energyhist[2]).fill(mcparticle.getEnergy());
+ aida.histogram1D("MC Data/all/Cos Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(mcparticle.getMomentum()));
+ //Add muons to specific list
+ if(gettype(mcparticle.getPDGID())==13){
+ MCMuons.add(mcparticle);
+ }
+ //Set if this is a beampipe event or not
+ if(VecOp.cosTheta(mcparticle.getMomentum())>0.99){
+ beampipe = 1;
+ }
+ }
+ }
+ aida.histogram1D("MC Data/Total Energy",energyhist[0],energyhist[1],energyhist[2]).fill(energy);
+ return i;
+ }
+
+ protected int countrc() {
+ // Counts the number of rc particles and adds data to a histogram and particlenumbersrc
+ int i = 0;
+ double energy = 0.0;
+ particlenumbersrc = new HashMap<Integer, Integer>();
+ for (ReconstructedParticle rcparticle : rcparticles) {
+ if (particlenumbersrc.get(gettype(rcparticle.getParticleIDUsed().getPDG())) == null) {
+ particlenumbersrc.put(gettype(rcparticle.getParticleIDUsed().getPDG()), 1);
+ } else {
+ particlenumbersrc.put(gettype(rcparticle.getParticleIDUsed().getPDG()), particlenumbersrc.get(gettype(rcparticle.getParticleIDUsed().getPDG())) + 1);
+ }
+ i++;
+ energy+=rcparticle.getEnergy();
+ //Muons added to specific list for use throughout the program
+ if (rcparticle.getParticleIDUsed().getPDG() == 13) {
+ RCMuons.add(rcparticle);
+ }
+ aida.histogram1D("RC Data/" + gettype(rcparticle.getParticleIDUsed().getPDG()) + "/Energy",energyhist[0],energyhist[1],energyhist[2]).fill(rcparticle.getEnergy());
+ aida.histogram1D("RC Data/" + gettype(rcparticle.getParticleIDUsed().getPDG()) + "/Cos Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(rcparticle.getMomentum()));
+ aida.histogram1D("RC Data/all/Energy",energyhist[0],energyhist[1],energyhist[2]).fill(rcparticle.getEnergy());
+ aida.histogram1D("RC Data/all/Cos Angle",anglehist[0],anglehist[1],anglehist[2]).fill(VecOp.cosTheta(rcparticle.getMomentum()));
+ }
+ aida.histogram1D("RC Data/Total Energy",energyhist[0],energyhist[1],energyhist[2]).fill(energy);
+ return i;
+ }
+
+ protected void endOfData() {
+// It is possible to output values from the HashMaps used to store the data at the end of a run for more precise data analysis
+// Object ParticleList2[] = particlenumbersmc.keySet().toArray();
+// for(int i = 1; i< ParticleList2.length;i++){
+// System.out.println("Particle ID = " + ParticleList2[i].hashCode() + " number of particles = " + particlenumbersmc.get(ParticleList2[i]));
+// }
+ }
+}
