29 added + 1 removed + 5 modified, total 35 files
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N pi0Analysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ pi0Analysis.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,127 @@
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.VecOp;
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.Cluster;
+import org.lcsim.util.Driver;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.mc.fast.MCFast;
+/*
+
+ */
+
+public class pi0Analysis extends Driver {
+
+ public pi0Analysis() {
+ add(new MCFast());
+ }
+ private AIDA aida = AIDA.defaultInstance();
+
+ @Override
+ protected void process(EventHeader event) {
+ super.process(event);
+ final double piplusmass = 0.13957;
+ // Get the list of MCParticles from the event
+ List<MCParticle> particles = event.get(MCParticle.class, EventHeader.MC_PARTICLES);
+ List<Track> tracks = event.getTracks();
+ List<Cluster> clusters = event.getClusters();
+ // Histogram the number of particles per event
+ // Loop over the particles
+ int found = 0;
+ int ntracks = 0;
+ double sumE = 0.0;
+ Hep3Vector sumP = new BasicHep3Vector(0., 0., 0.);
+ for (MCParticle particle : particles) {
+ if (Math.abs(particle.getPDGID()) == 111) {
+ aida.cloud1D("pi0invmass").fill(particle.getMass());
+ }
+ if (particle.getGeneratorStatus() == MCParticle.FINAL_STATE) {
+ ntracks++;
+ sumE = sumE + particle.getEnergy();
+ sumP = VecOp.add(sumP, particle.getMomentum());
+ }
+ } // end loop over particles
+ double psquare = VecOp.dot(sumP, sumP);
+ double ainvmass = Math.sqrt(sumE * sumE - psquare);
+ aida.cloud1D("ainvmass").fill(ainvmass);
+ aida.cloud1D("nTracks").fill(ntracks);
+
+ double ptotal = 0;
+ double energy;
+ double[] pos;
+ double length;
+
+// double sumRE = 0.0;
+// Hep3Vector sumRP = new BasicHep3Vector(0., 0., 0.);
+ for (int ii = 0; ii < clusters.size(); ii++) {
+ Cluster cluster = clusters.get(ii);
+ List<CalorimeterHit> calorimeterHits = cluster.getCalorimeterHits();
+// System.out.println("Hits:"+calorimeterHits.size());
+ for (CalorimeterHit hit : calorimeterHits) {
+// System.out.println("subdet"+hit.getSubdetector());
+ double[] rvec = hit.getPosition();
+ double tof = hit.getTime();
+ double r = Math.sqrt(rvec[0]*rvec[0] + rvec[1]*rvec[1] + rvec[2]*rvec[2]);
+ double t0 = r/2.99792485e10;
+ double t_corr = tof - t0;
+ aida.cloud1D("cluster hit times").fill(t_corr);
+ aida.cloud1D("cluster edep time").fill(t_corr*hit.getRawEnergy());
+ }
+// energy = cluster.getEnergy();
+// sumRE = sumRE + energy;
+// pos = cluster.getPosition();
+// length = Math.sqrt(pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]);
+// Hep3Vector gammamom = new BasicHep3Vector(pos[0] * energy / length, pos[1] * energy / length, pos[2] * energy / length);
+// sumRP = VecOp.add(sumRP, gammamom);
+ }
+// double rpsquare = VecOp.dot(sumRP, sumRP);
+// double Rinvmass = Math.sqrt(sumRE * sumRE - rpsquare);
+// aida.cloud1D("recinvmass").fill(Rinvmass);
+
+ List<ReconstructedParticle> rcps = event.get(ReconstructedParticle.class, "PandoraPFOCollection");
+ List<ReconstructedParticle> rcgammas = new ArrayList<ReconstructedParticle>();
+ for (ReconstructedParticle rcp:rcps) {
+ aida.cloud1D("pandora/rcp energy").fill(rcp.getEnergy());
+ aida.cloud1D("pandora/rcp PDGID").fill(rcp.getType());
+ if (rcp.getType()==22) {
+ rcgammas.add(rcp);
+ aida.cloud1D("pandora/rc gamma energy").fill(rcp.getEnergy());
+ }
+ }
+ aida.cloud1D("pandora/num rcps").fill(rcps.size());
+ aida.cloud1D("pandora/num rc gammas").fill(rcgammas.size());
+ aida.cloud1D("pandora/num rcgammas over rcps").fill(rcgammas.size()/rcps.size());
+
+ double sumRE = 0.0;
+ Hep3Vector sumRP = new BasicHep3Vector(0.,0.,0.);
+ if (rcgammas.size() == 2) {
+ Hep3Vector p1 = rcgammas.get(0).getMomentum();
+ Hep3Vector p2 = rcgammas.get(1).getMomentum();
+ for (ReconstructedParticle rg:rcgammas) {
+ sumRE += rg.getEnergy();
+ sumRP = VecOp.add(sumRP, rg.getMomentum());
+ aida.cloud1D("pandora/single-gamma energy").fill(rg.getEnergy());
+ }
+ double rpsquare = sumRP.magnitudeSquared();
+ double rinvmass = Math.sqrt(sumRE*sumRE - rpsquare);
+ aida.cloud1D("pandora/two-gamma sum E").fill(sumE);
+ aida.cloud1D("pandora/two-gamma sum P*P").fill(Math.sqrt(rpsquare));
+ aida.cloud1D("pandora/two-gamma invariant mass").fill(rinvmass);
+
+ double costhetagg = VecOp.dot(p1, p2)/(p1.magnitude()*p2.magnitude());
+ aida.cloud1D("pandora/two-gamma invariant mass v2").fill(Math.sqrt(2.0*p1.magnitude()*p2.magnitude()*(1.0-costhetagg)));
+ }
+ if (rcgammas.size() == 1) {
+ ReconstructedParticle g = rcgammas.get(0);
+ aida.cloud1D("pandora/one-found energy").fill(g.getEnergy());
+ }
+
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N TmpDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TmpDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,292 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+//package org.lcsim.mcd.analysis;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.contrib.Cassell.recon.analysis.MakeMCParticlePlots;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.mc.fast.MCFast;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.util.JetDriver;
+import hep.physics.jet.JetFinder;
+import hep.physics.jet.DurhamJetFinder;
+/**
+ *
+ * @author aconway
+ */
+public class TmpDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ MakeMCParticlePlots mcp;
+
+ public TmpDriver(){
+ mcp = new MakeMCParticlePlots();
+ add(new MCFast());
+ JetFinder jFind = new DurhamJetFinder(0.7);
+ JetDriver jDrive = new JetDriver();
+ jDrive.setFinder(jFind);
+ add(jDrive);
+ }
+
+ public void process(EventHeader event) {
+ //super.process(event);
+ this.processChildren(event);
+ MCParticle Boson = findBoson(event);
+ List<MCParticle> BChilds = getBosonChildren(Boson);
+ List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class,"Jets");
+
+ boolean isZgamma = (Boson.getMass()<=124.9);//hasHighEnGamma(event);
+
+ if(isZgamma){
+ aida.cloud1D("myPlots/ZGamma/num Jets").fill(jets.size());
+ aida.cloud1D("myPlots/ZGamma/Z Mass").fill(Boson.getMass());
+ aida.cloud1D("myPlots/ZGamma/Z CosTheta").fill(VecOp.cosTheta(Boson.getMomentum()));
+ Hep3Vector q1_pt = new BasicHep3Vector();
+ Hep3Vector q2_pt = new BasicHep3Vector();
+ boolean isQuark = false;
+ boolean isB = false;
+ for(MCParticle p:BChilds){
+ if (Math.abs(p.getPDGID()) == 5){
+ isB = true;
+ }
+ if (Math.abs(p.getPDGID()) <= 8) {
+ if (isQuark) q2_pt = p.getMomentum();
+ else {
+ q1_pt = p.getMomentum();
+ isQuark = true;
+ }
+ }
+ aida.cloud1D("myPlots/ZGamma/Z Child PDGID's").fill(Math.abs(p.getPDGID()));
+ }
+ if (isQuark) {
+ aida.cloud1D("myPlots/ZGamma/Z->qq: q dot q").fill(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude()));
+ aida.cloud1D("myPlots/ZGamma/Z->qq: q-q angle").fill(
+ Math.acos(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude())));
+ aida.cloud1D("myPlots/ZGamma/Z->qq: q Momenta").fill(q1_pt.magnitude());
+ aida.cloud1D("myPlots/ZGamma/Z->qq: q Momenta").fill(q2_pt.magnitude());
+ aida.cloud1D("myPlots/ZGamma/Z->qq: q CosTheta").fill(VecOp.cosTheta(q1_pt));
+ aida.cloud1D("myPlots/ZGamma/Z->qq: q CosTheta").fill(VecOp.cosTheta(q2_pt));
+ jetPlots("myPlots/ZGamma/Z->qq: jets/",jets);
+ }
+ if (isB) {
+ aida.cloud1D("myPlots/ZGamma/Z->bb: b dot b").fill(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude()));
+ aida.cloud1D("myPlots/ZGamma/Z->bb: b-b angle").fill(
+ Math.acos(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude())));
+ aida.cloud1D("myPlots/ZGamma/Z->bb: b Momenta").fill(q1_pt.magnitude());
+ aida.cloud1D("myPlots/ZGamma/Z->bb: b Momenta").fill(q2_pt.magnitude());
+ aida.cloud1D("myPlots/ZGamma/Z->bb: b CosTheta").fill(VecOp.cosTheta(q1_pt));
+ aida.cloud1D("myPlots/ZGamma/Z->bb: b CosTheta").fill(VecOp.cosTheta(q2_pt));
+ jetPlots("myPlots/ZGamma/Z->bb: jets/",jets);
+ }
+ }
+ else {
+ aida.cloud1D("myPlots/Zqq/Z Mass").fill(Boson.getMass());
+ aida.cloud1D("myPlots/Zqq/Z CosTheta").fill(VecOp.cosTheta(Boson.getMomentum()));
+ Hep3Vector q1_pt = new BasicHep3Vector();
+ Hep3Vector q2_pt = new BasicHep3Vector();
+ boolean isQuark = false;
+ boolean isB = false;
+ for(MCParticle p:BChilds){
+ if (Math.abs(p.getPDGID()) == 5){
+ isB = true;
+ }
+ if (Math.abs(p.getPDGID()) <= 8) {
+ if (isQuark) q2_pt = p.getMomentum();
+ else {
+ q1_pt = p.getMomentum();
+ isQuark = true;
+ }
+ }
+ aida.cloud1D("myPlots/Zqq/Z Child PDGID's").fill(Math.abs(p.getPDGID()));
+ }
+ if (isQuark) {
+ aida.cloud1D("myPlots/Zqq/Z->qq: q dot q").fill(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude()));
+ aida.cloud1D("myPlots/Zqq/Z->qq: q-q angle").fill(
+ Math.acos(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude())));
+ aida.cloud1D("myPlots/Zqq/Z->qq: q Momenta").fill(q1_pt.magnitude());
+ aida.cloud1D("myPlots/Zqq/Z->qq: q Momenta").fill(q2_pt.magnitude());
+ aida.cloud1D("myPlots/Zqq/Z->qq: q CosTheta").fill(VecOp.cosTheta(q1_pt));
+ aida.cloud1D("myPlots/Zqq/Z->qq: q CosTheta").fill(VecOp.cosTheta(q2_pt));
+ aida.histogram1D("myPlots/Zqq/z->qq: q Momenta Hist",50,62.2,62.6).fill(q1_pt.magnitude());
+ aida.histogram1D("myPlots/Zqq/z->qq: q Momenta Hist",50,62.2,62.6).fill(q2_pt.magnitude());
+ jetPlots("myPlots/Zqq/Z->qq: jets/",jets);
+ }
+ if (isB) {
+ aida.cloud1D("myPlots/Zqq/Z->bb: b dot b").fill(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude()));
+ aida.cloud1D("myPlots/Zqq/Z->bb: b-b angle").fill(
+ Math.acos(VecOp.dot(q1_pt,q2_pt)/(q1_pt.magnitude()*q2_pt.magnitude())));
+ aida.cloud1D("myPlots/Zqq/Z->bb: b Momenta").fill(q1_pt.magnitude());
+ aida.cloud1D("myPlots/Zqq/Z->bb: b Momenta").fill(q2_pt.magnitude());
+ aida.cloud1D("myPlots/Zqq/Z->bb: b CosTheta").fill(VecOp.cosTheta(q1_pt));
+ aida.cloud1D("myPlots/Zqq/Z->bb: b CosTheta").fill(VecOp.cosTheta(q2_pt));
+ aida.histogram1D("myPlots/Zqq/z->bb: b Momenta Hist",50,62.2,62.6).fill(q1_pt.magnitude());
+ aida.histogram1D("myPlots/Zqq/z->bb: b Momenta Hist",50,62.2,62.6).fill(q2_pt.magnitude());
+ jetPlots("myPlots/Zqq/Z->bb: jets/",jets);
+ }
+ }
+
+ aida.cloud1D("All boson masses").fill(Boson.getMass());
+ if (Math.abs(BChilds.get(0).getPDGID()) >=0 ) {
+ //if (Boson.getMass() < 110.0) {
+ List<List<MCParticle>> mcll = event.get(MCParticle.class);
+ for(List<MCParticle> mcl:mcll)
+ {
+ String mcs = event.getMetaData(mcl).getName();
+ if(mcs.indexOf("Final") >= 0)
+ {
+ mcp.makePlots(mcs,mcl,0.993);
+ List<MCParticle> mcl2 = new ArrayList<MCParticle>();
+ int hiEgammaCount = 0;
+ for (MCParticle p:mcl){
+ int pdgid = p.getPDGID();
+ double en = p.getEnergy();
+ if( !( (pdgid == 22) && (en >= 20.0)) ) {
+ mcl2.add(p);
+ } else {
+ hiEgammaCount++;
+ }
+ }
+ aida.cloud1D("HighEGammasPerEvent").fill(hiEgammaCount);
+ if (isZgamma) {
+ mcp.makePlots("ZGamma", mcl, 0.993);
+ aida.cloud1D("ZGamma Z Masses").fill(Boson.getMass());
+ for (MCParticle p:mcl) {
+ if((p.getPDGID() == 22) &&
+ (p.getGeneratorStatus() == 1) &&
+ (p.getEnergy() >= 20.0) ) {
+ aida.cloud1D("myPlots/ZGamma/HighEnGamma En").fill(p.getEnergy());
+ aida.cloud1D("myPlots/ZGamma/HighEnGamma CosTheta").fill(VecOp.cosTheta(p.getMomentum()));
+ aida.cloud2D("myPlots/ZGamma/Z-gamma Energy Correlation").fill(p.getEnergy(), Boson.getEnergy());
+ }
+ }
+ }
+ else {
+ mcp.makePlots("Zqq",mcl,0.993);
+ aida.cloud1D("Zqq Z Masses").fill(Boson.getMass());
+ }
+ }
+ }
+ //}
+
+ }
+
+// List<List<MCParticle>> mcll = event.get(MCParticle.class);
+// for(List<MCParticle> mcl:mcll)
+// {
+// String mcs = event.getMetaData(mcl).getName();
+// if(mcs.indexOf("Final") >= 0)
+// {
+// mcp.makePlots(mcs,mcl,0.993);
+// }
+// }
+ }
+
+ public MCParticle findBoson(EventHeader event) {
+ MCParticle bos = event.getMCParticles().get(0);
+ for (MCParticle p:event.getMCParticles()) {
+ if ( (p.getPDGID() == 23) && (p.getGeneratorStatus() == 3) ){
+ bos = p;
+ break;
+ }
+ else if ( (p.getPDGID() == 25) && (p.getGeneratorStatus() == 3) ){
+ bos = p;
+ break;
+ }
+ }
+ return bos;
+ }
+
+ public List<MCParticle> getBosonChildren(MCParticle parent) {
+ List<MCParticle> children = parent.getDaughters();
+ List<MCParticle> rv = new ArrayList<MCParticle>();
+ for (MCParticle p:children) {
+ if (p.getGeneratorStatus() == 3) {
+ rv.add(p);
+ }
+ }
+ return rv;
+ }
+
+ private boolean hasHighEnGamma(EventHeader event) {
+ List<MCParticle> mcpl = event.getMCParticles();
+ int count = 0;
+ for (MCParticle mcp:mcpl) {
+ int pdgid = mcp.getPDGID();
+ int status = mcp.getGeneratorStatus();
+ double en = mcp.getEnergy();
+ if ( (pdgid == 22) && (status == 1) && (en >= 15.0) && (Math.abs(mcp.getParents().get(0).getPDGID()) == 13)) {
+ count++;
+ }
+ }
+ return (count == 1);
+ }
+
+ private void jetPlots(String prefix, List<ReconstructedParticle> jets) {
+ aida.cloud1D(prefix+"num Jets").fill(jets.size());
+ if (jets.size() > 0) {
+ double sumE = 0.0;
+ double sumM = 0.0;
+ double sumP = 0.0;
+ Hep3Vector sumPvec = new BasicHep3Vector(0.0,0.0,0.0);
+ int sumChrg = 0;
+ for(ReconstructedParticle jet:jets) {
+ aida.cloud1D(prefix+"jet energy").fill(jet.getEnergy());
+ aida.cloud1D(prefix+"jet mass").fill(jet.getMass());
+ aida.cloud1D(prefix+"jet momentum").fill(jet.getMomentum().magnitude());
+ aida.cloud1D(prefix+"jet cosTheta").fill(VecOp.cosTheta(jet.getMomentum()));
+ aida.cloud1D(prefix+"jet charge").fill(jet.getCharge());
+ sumE += jet.getEnergy();
+ sumM += jet.getMass();
+ sumP += jet.getMomentum().magnitude();
+ sumChrg += jet.getCharge();
+ VecOp.add(sumPvec,jet.getMomentum());
+ }
+ aida.cloud1D(prefix+"total jet energy").fill(sumE);
+ aida.cloud1D(prefix+"total Jet mass").fill(sumM);
+ aida.cloud1D(prefix+"total jet momentum").fill(sumP);
+ aida.cloud1D(prefix+"total jet charge").fill(sumChrg);
+ aida.cloud1D(prefix+"net momentum vec").fill(sumPvec.magnitude());
+ aida.cloud1D(prefix+"net momentum cosTheta").fill(VecOp.cosTheta(sumPvec));
+ }
+ if (jets.size() == 2) {
+ double sumE = 0.0;
+ double sumM = 0.0;
+ double sumP = 0.0;
+ Hep3Vector p1 = new BasicHep3Vector(0.0,0.0,0.0);
+ Hep3Vector p2 = new BasicHep3Vector(0.0,0.0,0.0);
+ double e1 = 0.0;
+ double e2 = 0.0;
+ int count = 0;
+ for(ReconstructedParticle jet:jets) {
+ aida.cloud1D(prefix+"dijet energy").fill(jet.getEnergy());
+ aida.cloud1D(prefix+"dijet mass").fill(jet.getMass());
+ aida.cloud1D(prefix+"dijet momentum").fill(jet.getMomentum().magnitude());
+ aida.cloud1D(prefix+"dijet cosTheta").fill(VecOp.cosTheta(jet.getMomentum()));
+ aida.cloud1D(prefix+"dijet charge").fill(jet.getCharge());
+ sumE += jet.getEnergy();
+ sumM += jet.getMass();
+ sumP += jet.getMomentum().magnitude();
+ if(count == 0){
+ p1 = jet.getMomentum();
+ e1 = jet.getEnergy();
+ } else {
+ p2 = jet.getMomentum();
+ e2 = jet.getEnergy();
+ }
+ }
+ aida.cloud1D(prefix+"total dijet energy").fill(sumE);
+ aida.cloud1D(prefix+"total dijet mass").fill(sumM);
+ aida.cloud1D(prefix+"total dijet momentum").fill(sumP);
+ aida.cloud1D(prefix+"dijet invariant mass"
+ ).fill(Math.sqrt((e1+e2)*(e1+e2) - VecOp.add(p1, p2).magnitudeSquared()));
+ }
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N EnergyCorrectionDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EnergyCorrectionDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,77 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+//package org.lcsim.mcd.analysis;
+
+import hep.aida.ICloud1D;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class EnergyCorrectionDriver extends Driver {
+
+ private AIDA aida;
+
+ public EnergyCorrectionDriver() {
+ aida = AIDA.defaultInstance();
+ }
+
+ protected void process(EventHeader event) {
+ double E_in = 0.0;
+ double E_kin = 0.0;
+ String Particlename = null;
+ List<MCParticle> particles = event.get(MCParticle.class, "MCParticle");
+ for (MCParticle particle : particles) {
+ if (particle.getGeneratorStatus() == 1) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+ if (particle.getProductionTime() == 0.0) {
+ Particlename = particle.getType().getName();
+ }
+ }
+ }
+ Integer Ein = (int) Math.floor(E_kin + 0.5d);
+ String E_str = Ein.toString();
+ String DirName = Particlename.concat(E_str);
+ ICloud1D Edep = aida.cloud1D(DirName+"/Edep",100000);
+ ICloud1D Eceren = aida.cloud1D(DirName+"/Eceren",100000);
+ List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+ double sumEEdep = 0.0;
+ double sumECeren = 0.0;
+ for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+ String CollectionName = event.getMetaData(simCalorimeterHits).getName();
+
+ if (CollectionName.contains("Edep")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ } // end if Edep
+ if (CollectionName.contains("Opti")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ } // end if Ceren
+ } // end loop over calorimeter hit collections
+
+ Edep.fill(sumEEdep);
+ Eceren.fill(sumECeren);
+ if (E_in != 0) {
+ aida.cloud2D("Ein vs Edep").fill(E_in, sumEEdep / E_in);
+ aida.cloud2D("Ein vs Eceren").fill(E_in, sumECeren / E_in);
+ aida.cloud2D("Ein vs Edep+ECeren").fill(E_in, (sumEEdep + sumECeren) / E_in);
+ if (sumEEdep != 0 && sumECeren != 0) {
+ aida.cloud2D("Edep over Eceren").fill(E_in, sumEEdep / sumECeren);
+ }
+ }
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N OpticalHitsTesterDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ OpticalHitsTesterDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,67 @@
+//package org.lcsim.mcd.analysis;
+
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class OpticalHitsTesterDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ protected void process(EventHeader event) {
+ List<MCParticle> mcpList = event.get(MCParticle.class, event.MC_PARTICLES);
+ List<SimCalorimeterHit> EdepHitList = new ArrayList<SimCalorimeterHit>();
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "ECalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "HCalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "ECalEndcapEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "HCalEndcapEdepHits"));
+ List<SimCalorimeterHit> OptiHitList = new ArrayList<SimCalorimeterHit>();
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "ECalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "HCalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "ECalEndcapOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "HCalEndcapOptiHits"));
+
+ double totMCP = 0.0;
+// for (List<MCParticle> mcps : mcpList) {
+// for (MCParticle mcp : mcps) {
+// if (mcp.getGeneratorStatus() == 1) {
+// totMCP += mcp.getEnergy();
+// }
+// }
+// }
+ for (MCParticle mcp : mcpList) {
+ if (mcp.getGeneratorStatus() == 0) {
+ totMCP += mcp.getEnergy();
+ }
+ }
+ double totEdep = 0.0;
+ double totOpti = 0.0;
+ for (SimCalorimeterHit hit : EdepHitList) {
+ totEdep += hit.getRawEnergy();
+ }
+ for (SimCalorimeterHit hit : OptiHitList) {
+ totOpti += hit.getRawEnergy();
+ }
+
+ System.out.println(totMCP);
+ System.out.println(totEdep);
+ System.out.println(totOpti);
+ aida.cloud2D("total Edep vs Ein").fill(totMCP,totEdep);
+ aida.cloud2D("total Opti vs Ein").fill(totMCP,totOpti);
+ aida.cloud2D("total All vs Ein").fill(totMCP,totEdep+totOpti);
+ aida.cloud2D("ratio Edep vs Ein").fill(totMCP, totEdep/totMCP);
+ aida.cloud2D("ratio Opti vs Ein").fill(totMCP, totOpti/totMCP);
+ aida.cloud2D("ratio All vs Ein").fill(totMCP, (totEdep+totOpti)/totMCP);
+ }
+
+
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N EventShape2.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EventShape2.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,575 @@
+//package hep.physics.jet;
+
+import hep.physics.filter.Predicate;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.HepLorentzVector;
+
+import java.util.Collection;
+import java.util.Iterator;
+
+/**
+ * Event Shape and Thrust utilities.
+ * <p>
+ * This is a transcription of the Jetset thrust and event shape
+ * finders into Java.
+ * <p>
+ * Event shape extracts from the input enumeration 3-momenta which are formed
+ * into a kind of (symmetric) momentum tensor similar to an inertia tensor.
+ * From this tensor the 3-principal axes are determined along with their
+ * associated eigenvalues.
+ * <p>
+ * Traditionally, the nomenclature for the three axes are:
+ * <ul>
+ * <li>Thrust Axis is associated with the largest eigenvalue called the Thrust.
+ * <li>Major Axis is associated with the middle eigenvalue called the Major Thrust.
+ * <li>Minor Axis is associated with the smallest eigenvalue called the Minor Thrust.
+ * </ul>
+ *
+ * @author G.Bower 9/21/98
+ */
+
+public class EventShape2
+{
+ // data: parameters
+ // PARU(41): Power of momentum dependence in sphericity finder.
+ private double m_dSphMomPower = 2.0;
+ // PARU(42): Power of momentum dependence in thrust finder.
+ private double m_dDeltaThPower = 0.0;
+ // MSTU(44): # of initial fastest particles choosen to start search.
+ private int m_iFast = 4;
+ // PARU(48): Convergence criteria for axis maximization.
+ private double m_dConv = 0.0001;
+ // MSTU(45): # different starting configurations that must
+ // converge before axis is accepted as correct.
+ private int m_iGood = 2;
+ // data: results
+ // m_dAxes[0] is the Thrust axis.
+ // m_dAxes[1] is the Major axis.
+ // m_dAxes[2] is the Minor axis.
+ private double[][] m_dAxes;
+ private double[] m_dThrust;
+ private double m_dOblateness;
+ private BasicHep3Vector m_EigenVector1;
+ private BasicHep3Vector m_EigenVector2;
+ private BasicHep3Vector m_EigenVector3;
+ private double m_dEigenValue1;
+ private double m_dEigenValue2;
+ private double m_dEigneValue3;
+
+ private final static int m_maxpart = 1000;
+
+ /**
+ * Call this to input a new event to the event shape routines.
+ * @param e An Enumeration of either HepLorentzVectors (HepLorentzVectors) or Hep3Vectors
+ */
+ public EventShape2()
+ {
+ // The zero element in each array in m_dAxes is ignored. Elements
+ // 1,2 and 3 are the x, y, and z direction cosines of the axis.
+ // Also the zeroth axis and thrust are ignored.
+ m_dAxes = new double[4][4];
+ m_dThrust = new double[4];
+ }
+ /**
+ * Call this to input a new event to the event shape routines.
+ *
+ * @param e An Enumeration of either HepLorentzVectors (HepLorentzVectors) or Hep3Vectors
+ */
+ public void setEvent(Collection data)
+ {
+ setEvent(data,null);
+ }
+ /**
+ * Call this to input a new event to the event shape routines.
+ *
+ * Only elements of the enumeration which are accepted by the predicate will be used
+ * for jet finding.
+ *
+ * @param e An Enumeration of either HepLorentzVectors (HepLorentzVectors) or Hep3Vectors
+ * @param cut A predicate that is applied to each element of e, or null to accept all elements
+ */
+ public void setEvent(Collection data, Predicate cut)
+ {
+
+ //To make this look like normal physics notation the
+ //zeroth element of each array, mom[i][0], will be ignored
+ //and operations will be on elements 1,2,3...
+ double[][] mom = new double[m_maxpart][6];
+ double tmax = 0;
+ double phi = 0.;
+ double the = 0.;
+ double sgn;
+ double[][] fast = new double[ m_iFast + 1 ][6];
+ double[][] work = new double[11][6];
+ double tdi[] = new double[4];
+ double tds;
+ double tpr[] = new double[4];
+ double thp;
+ double thps;
+ double[][] temp = new double[3][5];
+
+ int np = 0;
+ for (Iterator i = data.iterator(); i.hasNext(); )
+ {
+ Object o = i.next();
+ if (cut != null && !cut.accept(o)) continue;
+
+ if (np >= m_maxpart) throw new RuntimeException("Too many particles input to EventShape");
+
+ Hep3Vector v;
+ if (o instanceof Hep3Vector)
+ {
+ v = (Hep3Vector) o;
+ }
+ else if (o instanceof HepLorentzVector)
+ {
+ HepLorentzVector l = (HepLorentzVector) o;
+ v = l.v3();
+ }
+ else throw new RuntimeException("Element input to EventShape is not a Hep3Vector or an HepLorentzVector");
+
+ mom[np][1] = v.x();
+ mom[np][2] = v.y();
+ mom[np][3] = v.z();
+ mom[np][4] = v.magnitude();
+ if ( Math.abs( m_dDeltaThPower ) <= 0.001 )
+ {
+ mom[np][5] = 1.0;
+ }
+ else
+ {
+ mom[np][5] = Math.pow( mom[np][4], m_dDeltaThPower );
+ }
+ tmax = tmax + mom[np][4]*mom[np][5];
+ np++;
+ }
+ if ( np < 2 )
+ {
+ m_dThrust[1] = -1.0;
+ m_dOblateness = -1.0;
+ return;
+ }
+ // for pass = 1: find thrust axis.
+ // for pass = 2: find major axis.
+ for ( int pass=1; pass < 3; pass++)
+ {
+ if ( pass == 2 )
+ {
+ phi = ulAngle( m_dAxes[1][1], m_dAxes[1][2] );
+ ludbrb( mom, 0, -phi, 0., 0., 0. );
+ for ( int i = 0; i < 3; i++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ temp[i][j] = m_dAxes[i+1][j];
+ }
+ temp[i][4] = 0;
+ }
+ ludbrb(temp,0.,-phi,0.,0.,0.);
+ for ( int i = 0; i < 3; i++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ m_dAxes[i+1][j] = temp[i][j];
+ }
+ }
+ the = ulAngle( m_dAxes[1][3], m_dAxes[1][1] );
+ ludbrb( mom, -the, 0., 0., 0., 0. );
+ for ( int i = 0; i < 3; i++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ temp[i][j] = m_dAxes[i+1][j];
+ }
+ temp[i][4] = 0;
+ }
+ ludbrb(temp,-the,0.,0.,0.,0.);
+ for ( int i = 0; i < 3; i++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ m_dAxes[i+1][j] = temp[i][j];
+ }
+ }
+ }
+ for ( int ifas = 0; ifas < m_iFast + 1 ; ifas++ )
+ {
+ fast[ifas][4] = 0.;
+ }
+ // Find the m_iFast highest momentum particles and
+ // put the highest in fast[0], next in fast[1],....fast[m_iFast-1].
+ // fast[m_iFast] is just a workspace.
+ for ( int i = 0; i < np; i++ )
+ {
+ if ( pass == 2 )
+ {
+ mom[i][4] = Math.sqrt( mom[i][1]*mom[i][1]
+ + mom[i][2]*mom[i][2] );
+ }
+ for ( int ifas = m_iFast - 1; ifas > -1; ifas-- )
+ {
+ if ( mom[i][4] > fast[ifas][4] )
+ {
+ for ( int j = 1; j < 6; j++ )
+ {
+ fast[ifas+1][j] = fast[ifas][j];
+ if ( ifas == 0 )
+ {
+ fast[ifas][j] = mom[i][j];
+ }
+ }
+ }
+ else
+ {
+ for ( int j = 1; j < 6; j++ )
+ {
+ fast[ifas+1][j] = mom[i][j];
+ }
+ break;
+ }
+ }
+ }
+ // Find axis with highest thrust (case 1)/ highest major (case 2).
+ for ( int i = 0; i < work.length; i++ )
+ {
+ work[i][4] = 0.;
+ }
+ int p = Math.min( m_iFast, np ) - 1;
+ // Don't trust Math.pow to give right answer always.
+ // Want nc = 2**p.
+ int nc = iPow(2,p);
+ for ( int n = 0; n < nc; n++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ tdi[j] = 0.;
+ }
+ for ( int i = 0; i < Math.min(m_iFast,n); i++ )
+ {
+ sgn = fast[i][5];
+ if (iPow(2,(i+1))*((n+iPow(2,i))/iPow(2,(i+1))) >= i+1)
+ {
+ sgn = -sgn;
+ }
+ for ( int j = 1; j < 5-pass; j++ )
+ {
+ tdi[j] = tdi[j] + sgn*fast[i][j];
+ }
+ }
+ tds = tdi[1]*tdi[1] + tdi[2]*tdi[2] + tdi[3]*tdi[3];
+ for ( int iw = Math.min(n,9); iw > -1; iw--)
+ {
+ if( tds > work[iw][4] )
+ {
+ for ( int j = 1; j < 5; j++ )
+ {
+ work[iw+1][j] = work[iw][j];
+ if ( iw == 0 )
+ {
+ if ( j < 4 )
+ {
+ work[iw][j] = tdi[j];
+ }
+ else
+ {
+ work[iw][j] = tds;
+ }
+ }
+ }
+ }
+ else
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ work[iw+1][j] = tdi[j];
+ }
+ work[iw+1][4] = tds;
+ }
+ }
+ }
+ // Iterate direction of axis until stable maximum.
+ m_dThrust[pass] = 0;
+ thp = -99999.;
+ int nagree = 0;
+ for ( int iw = 0; iw < Math.min(nc,10) && nagree < m_iGood; iw++ )
+ {
+ thp = 0.;
+ thps = -99999.;
+ while ( thp > thps + m_dConv )
+ {
+ thps = thp;
+ for ( int j = 1; j < 4; j++ )
+ {
+ if ( thp <= 1E-10 )
+ {
+ tdi[j] = work[iw][j];
+ }
+ else
+ {
+ tdi[j] = tpr[j];
+ tpr[j] = 0;
+ }
+ }
+ for ( int i = 0; i < np; i++ )
+ {
+ sgn = sign(mom[i][5],
+ tdi[1]*mom[i][1] +
+ tdi[2]*mom[i][2] +
+ tdi[3]*mom[i][3]);
+ for ( int j = 1; j < 5 - pass; j++ )
+ {
+ tpr[j] = tpr[j] + sgn*mom[i][j];
+ }
+ }
+ thp = Math.sqrt( tpr[1]*tpr[1]
+ + tpr[2]*tpr[2]
+ + tpr[3]*tpr[3])/tmax;
+ }
+ // Save good axis. Try new initial axis until enough
+ // tries agree.
+ if ( thp < m_dThrust[pass] - m_dConv )
+ {
+ break;
+ }
+ if ( thp > m_dThrust[pass] + m_dConv )
+ {
+ nagree = 0;
+ sgn = iPow( -1, (int)Math.round(Math.random()) );
+ for ( int j = 1; j < 4; j++ )
+ {
+ m_dAxes[pass][j] = sgn*tpr[j]/(tmax*thp);
+ }
+ m_dThrust[pass] = thp;
+ }
+ nagree = nagree + 1;
+ }
+ }
+ // Find minor axis and value by orthogonality.
+ sgn = iPow( -1, (int)Math.round(Math.random()));
+ m_dAxes[3][1] = -sgn*m_dAxes[2][2];
+ m_dAxes[3][2] = sgn*m_dAxes[2][1];
+ m_dAxes[3][3] = 0.;
+ thp = 0.;
+ for ( int i = 0; i < np; i++ )
+ {
+ thp = thp + mom[i][5]*Math.abs( m_dAxes[3][1]*mom[i][1]
+ + m_dAxes[3][2]*mom[i][2]);
+ }
+ m_dThrust[3] = thp/tmax;
+ // Rotate back to original coordinate system.
+ for ( int i = 0; i < 3; i++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ temp[i][j] = m_dAxes[i+1][j];
+ }
+ temp[i][4] = 0;
+ }
+ ludbrb(temp,the,phi,0.,0.,0.);
+ for ( int i = 0; i < 3; i++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ m_dAxes[i+1][j] = temp[i][j];
+ }
+ }
+ m_dOblateness = m_dThrust[2] - m_dThrust[3];
+ }
+ // Setting and getting parameters.
+ public void setThMomPower(double tp)
+ {
+ // Error if sp not positive.
+ if ( tp > 0. )
+ {
+ m_dDeltaThPower = tp - 1.0;
+ }
+ return;
+ }
+ public double getThMomPower()
+ {
+ return 1.0 + m_dDeltaThPower;
+ }
+ public void setFast(int nf)
+ {
+ // Error if sp not positive.
+ if ( nf > 3 )
+ {
+ m_iFast = nf;
+ }
+ return;
+ }
+ public int getFast()
+ {
+ return m_iFast;
+ }
+
+ // Returning results
+ public BasicHep3Vector thrustAxis()
+ {
+ return new BasicHep3Vector(m_dAxes[1][1],m_dAxes[1][2],m_dAxes[1][3]);
+ }
+ public BasicHep3Vector majorAxis()
+ {
+ return new BasicHep3Vector(m_dAxes[2][1],m_dAxes[2][2],m_dAxes[2][3]);
+ }
+ public BasicHep3Vector minorAxis()
+ {
+ return new BasicHep3Vector(m_dAxes[3][1],m_dAxes[3][2],m_dAxes[3][3]);
+ }
+ /**
+ * Element x = Thrust
+ * Element y = Major Thrust
+ * Element z = Minor Thrust
+ */
+ public BasicHep3Vector thrust()
+ {
+ return new BasicHep3Vector(m_dThrust[1],m_dThrust[2],m_dThrust[3]);
+ }
+ /**
+ * Oblateness = Major Thrust - Minor Thrust
+ */
+ public double oblateness()
+ {
+ return m_dOblateness;
+ }
+ // BasicHep3Vector eigenVector1()
+ // {
+ // return;
+ // }
+ // BasicHep3Vector eigenVector2()
+ // {
+ // return;
+ // }
+ // BasicHep3Vector eigenVector3()
+ // {
+ // return;
+ // }
+ // double eigenValue1()
+ // {
+ // return;
+ // }
+ // double eigenValue2()
+ // {
+ // return;
+ // }
+ // double eigenValue3()
+ // {
+ // return;
+ // }
+ // utilities(from Jetset):
+ private double ulAngle(double x, double y)
+ {
+ double ulangl = 0;
+ double r = Math.sqrt(x*x + y*y);
+ if ( r < 1.0E-20 )
+ {
+ return ulangl;
+ }
+ if ( Math.abs(x)/r < 0.8 )
+ {
+ ulangl = sign(Math.acos(x/r),y);
+ }
+ else
+ {
+ ulangl = Math.asin(y/r);
+ if ( x < 0. && ulangl >= 0. )
+ {
+ ulangl = Math.PI - ulangl;
+ }
+ else if ( x < 0. )
+ {
+ ulangl = -Math.PI - ulangl;
+ }
+ }
+ return ulangl;
+ }
+ private double sign(double a, double b)
+ {
+ if ( b < 0 )
+ {
+ return -Math.abs(a);
+ }
+ else
+ {
+ return Math.abs(a);
+ }
+ }
+ private void ludbrb(double[][] mom,
+ double the,
+ double phi,
+ double bx,
+ double by,
+ double bz)
+ {
+ // Ignore "zeroth" elements in rot,pr,dp.
+ // Trying to use physics-like notation.
+ double rot[][] = new double[4][4];
+ double pr[] = new double[4];
+ double dp[] = new double[5];
+ int np = mom.length;
+ if ( the*the + phi*phi > 1.0E-20 )
+ {
+ rot[1][1] = Math.cos(the)*Math.cos(phi);
+ rot[1][2] = -Math.sin(phi);
+ rot[1][3] = Math.sin(the)*Math.cos(phi);
+ rot[2][1] = Math.cos(the)*Math.sin(phi);
+ rot[2][2] = Math.cos(phi);
+ rot[2][3] = Math.sin(the)*Math.sin(phi);
+ rot[3][1] = -Math.sin(the);
+ rot[3][2] = 0.0;
+ rot[3][3] = Math.cos(the);
+ for ( int i = 0; i < np; i++ )
+ {
+ for ( int j = 1; j < 4; j++ )
+ {
+ pr[j] = mom[i][j];
+ mom[i][j] = 0;
+ }
+ for ( int j = 1; j < 4; j++)
+ {
+ for ( int k = 1; k < 4; k++)
+ {
+ mom[i][j] = mom[i][j] + rot[j][k]*pr[k];
+ }
+ }
+ }
+ double beta = Math.sqrt( bx*bx + by*by + bz*bz );
+ if ( beta*beta > 1.0E-20 )
+ {
+ if ( beta > 0.99999999 )
+ {
+ //send message: boost too large, resetting to <~1.0.
+ bx = bx*(0.99999999/beta);
+ by = by*(0.99999999/beta);
+ bz = bz*(0.99999999/beta);
+ beta = 0.99999999;
+ }
+ double gamma = 1.0/Math.sqrt(1.0 - beta*beta);
+ for ( int i = 0; i < np; i++ )
+ {
+ for ( int j = 1; j < 5; j++ )
+ {
+ dp[j] = mom[i][j];
+ }
+ double bp = bx*dp[1] + by*dp[2] + bz*dp[3];
+ double gbp = gamma*(gamma*bp/(1.0 + gamma) + dp[4]);
+ mom[i][1] = dp[1] + gbp*bx;
+ mom[i][2] = dp[2] + gbp*by;
+ mom[i][3] = dp[3] + gbp*bz;
+ mom[i][4] = gamma*(dp[4] + bp);
+ }
+ }
+ }
+ return;
+ }
+ private int iPow(int man, int exp)
+ {
+ int ans = 1;
+ for( int k = 0; k < exp; k++)
+ {
+ ans = ans*man;
+ }
+ return ans;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N bbbbAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ bbbbAnalysisDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,255 @@
+package org.lcsim.mcd.analysis;
+
+import hep.aida.IHistogram1D;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.mc.fast.MCFast;
+import java.util.List;
+import java.util.ArrayList;
+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.Hep3Vector;
+import hep.physics.vec.VecOp;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.Track;
+
+/**
+ *
+ * @author aconway
+ */
+public class bbbbAnalysisDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ int num_Bs;
+ int num_Bs_accepted_TP;
+ int num_Bs_accepted_CL;
+ int num_Bs_accepted_MC;
+ int nevents;
+ int events_accepted_TP;
+ int events_accepted_CL;
+ int events_accepted_MC;
+
+ public bbbbAnalysisDriver() {
+ add(new MCFast(true,false,true));
+ num_Bs = 0;
+ num_Bs_accepted_TP = 0;
+ num_Bs_accepted_CL = 0;
+ num_Bs_accepted_MC = 0;
+ nevents = 0;
+ events_accepted_TP = 0;
+ events_accepted_CL = 0;
+ events_accepted_MC = 0;
+ }
+
+ public void process(EventHeader event) {
+ nevents++;
+ this.processChildren(event);
+ List<MCParticle> particles = event.get(MCParticle.class,"MCParticle");
+ List<LCRelation> tracks2MCP = event.get(LCRelation.class,"TracksToMCP");
+ List<Track> all_tracks = event.get(Track.class,"Tracks");
+
+ List<MCParticle> allBmesons = new ArrayList();
+
+ for (MCParticle p:particles) {
+ int pdgid = Math.abs(p.getPDGID());
+ if ( (p.getGeneratorStatus()==2)
+ && ( ( (pdgid%10000 > 500) && (pdgid%10000 < 600) )
+ ||( (pdgid > 5000) && (pdgid < 6000) ))
+ && (p.getOrigin().magnitude() == 0.0)
+ && (p.getEndPoint().magnitude() != 0.0) )
+ {
+ aida.cloud1D("Bs/B energy").fill(p.getEnergy());
+ aida.cloud1D("Bs/B cosTheta").fill(Math.abs(VecOp.cosTheta(p.getMomentum())));
+ aida.cloud1D("Bs/B distance").fill(p.getEndPoint().magnitude());
+ aida.histogram1D("Bs/B meson ctau",100,0,2).fill(p.getEndPoint().magnitude()*p.getMass()/p.getMomentum().magnitude());
+ aida.cloud1D("Bs/B meson ctau cloud").fill(p.getEndPoint().magnitude()*p.getMass()/p.getMomentum().magnitude());
+ aida.cloud1D("Bs/B meson p over m").fill(p.getMomentum().magnitude()/p.getMass());
+ allBmesons.add(p);
+ }
+ }
+ int Bs_accepted_TP = 0;
+ int Bs_accepted_CL = 0;
+ int Bs_accepted_MC = 0;
+ for (MCParticle B:allBmesons) {
+ num_Bs++;
+ List<Track> trks = new ArrayList();
+ List<MCParticle> trkmcps = getXChildren(B,5);
+ for (MCParticle p:trkmcps) {
+ Track trk = getTrack(p,tracks2MCP);
+ if (trk != null)
+ trks.add(trk);
+ }
+ int[] ntracks = trackPlots(trks,tracks2MCP,"B tracks/");
+ if (ntracks[0] >= 2) Bs_accepted_TP++;
+ if (ntracks[1] >= 2) Bs_accepted_CL++;
+ if (ntracks[2] >= 2) Bs_accepted_MC++;
+ }
+ aida.cloud1D("Bs accepted per event - TP").fill(Bs_accepted_TP);
+ aida.cloud1D("Bs accepted per event - CL").fill(Bs_accepted_CL);
+ aida.cloud1D("Bs accepted per event - MC").fill(Bs_accepted_MC);
+ if (Bs_accepted_TP >= 4) events_accepted_TP++;
+ if (Bs_accepted_CL >= 4) events_accepted_CL++;
+ if (Bs_accepted_MC >= 4) events_accepted_MC++;
+ //trackPlots(all_tracks, tracks2MCP, "all tracks/");
+ }
+
+ public int[] trackPlots(List<Track> all_tracks, List<LCRelation> tracks2MCP, String lbl) {
+ int tracks_accepted_TP = 0;
+ int tracks_accepted_CL = 0;
+ int tracks_accepted_MC = 0;
+ int tracks_rejected = 0;
+ int tracks = 0;
+ for (Track trk:all_tracks) {
+ tracks++;
+ double [] tpms = trk.getTrackParameters();
+ double dd0 = trk.getErrorMatrix().diagonal(0);
+ aida.histogram1D(lbl+"tracks/d0",50,-1.0,1.0).fill(tpms[0]);
+ aida.histogram1D(lbl+"tracks/d0 over error",200,-10,10).fill(tpms[0]/Math.sqrt(dd0));
+ aida.cloud1D(lbl+"tracks/phi0").fill(tpms[1]);
+ aida.cloud1D(lbl+"tracks/omega").fill(tpms[2]);
+ aida.cloud1D(lbl+"tracks/z0").fill(tpms[3]);
+ aida.cloud1D(lbl+"tracks/s").fill(tpms[4]);
+ aida.histogram1D(lbl+"tracks/costheta",100,0,1).fill(Math.abs(VecOp.cosTheta(new BasicHep3Vector(trk.getMomentum()))));
+
+ MCParticle mcp = null;
+ for (LCRelation lcr:tracks2MCP) {
+ if (trk.equals(lcr.getFrom())) {
+ mcp = (MCParticle) lcr.getTo();
+ break;
+ }
+ }
+ if (mcp != null) {
+ double d0 = tpms[0];
+ double z0 = tpms[3];
+ dd0 = Math.sqrt(trk.getErrorMatrix().diagonal(0));
+ double dz0 = Math.sqrt(trk.getErrorMatrix().diagonal(3));
+ double IP3D = Math.sqrt(d0*d0 + z0*z0);
+ double IP3D_cverr = IP3D*((dd0/d0) + (dz0/z0));
+ double IP3D_CLerr = getIPError(mcp,0.005,0.015);
+ double IP3D_MCerr = getIPError(mcp,0.008,0.025);
+ double z = mcp.getOriginZ();
+ double r = Math.sqrt(mcp.getOriginX()*mcp.getOriginX() + mcp.getOriginY()*mcp.getOriginY());
+ double rmag = mcp.getOrigin().magnitude();
+ double oCostheta = Math.abs(VecOp.cosTheta(mcp.getOrigin()));
+ double pCostheta = Math.abs(VecOp.cosTheta(mcp.getMomentum()));
+ if (rmag == 0.0) {
+ aida.histogram1D(lbl+"origin tracks/d0",50,-1,1).fill(d0);
+ aida.histogram1D(lbl+"origin tracks/dd0",50,0,1).fill(dd0);
+ aida.histogram1D(lbl+"origin tracks/z0",50,-1,1).fill(z0);
+ aida.histogram1D(lbl+"origin tracks/dz0",50,0,1).fill(dz0);
+ aida.histogram1D(lbl+"origin tracks/d0 over dd0",200,-10,10).fill(d0/dd0);
+ aida.histogram1D(lbl+"origin tracks/z0 over dz0",200,-10,10).fill(z0/dz0);
+ aida.histogram1D(lbl+"origin tracks/IP3D",50,0,2).fill(IP3D);
+ aida.histogram1D(lbl+"origin tracks/ID3D TP error",50,0,1).fill(IP3D_cverr);
+ aida.histogram1D(lbl+"origin tracks/IP3D CLIC error",50,0,1).fill(IP3D_CLerr);
+ aida.histogram1D(lbl+"origin tracks/IP3D MuCo error",50,0,1).fill(IP3D_MCerr);
+ aida.histogram1D(lbl+"origin tracks/IP3D over IP3D TP error",200,0,10).fill(IP3D/IP3D_cverr);
+ aida.histogram1D(lbl+"origin tracks/IP3D over IP3D CLIC error",200,0,10).fill(IP3D/IP3D_CLerr);
+ aida.histogram1D(lbl+"origin tracks/IP3D over IP3D MuCo error",200,0,10).fill(IP3D/IP3D_MCerr);
+ aida.histogram1D(lbl+"origin tracks/IP3D TP error over IP3D CLIC error",50,0,2).fill(IP3D_cverr/IP3D_CLerr);
+ aida.histogram1D(lbl+"origin tracks/IP3D TP error over IP3D MuCo error",50,0,2).fill(IP3D_cverr/IP3D_MCerr);
+ aida.histogram1D(lbl+"origin tracks/IP3D MuCo error over IP3D CLIC error",50,1.4,1.8).fill(IP3D_MCerr/IP3D_CLerr);
+ aida.histogram1D(lbl+"origin tracks/track momentum cos theta",50,0,1).fill(pCostheta);
+ aida.histogram1D(lbl+"origin tracks/track momentum",50,0,50).fill(mcp.getMomentum().magnitude());
+ } else {
+ if ( (z < 700.0) && (r < 220.0)
+ && (oCostheta < Math.cos(10.0*Math.PI/180.0))
+ && (pCostheta < Math.cos(10.0*Math.PI/180.0))) {
+ aida.histogram1D(lbl+"B tracks/d0",50,-1,1).fill(d0);
+ aida.histogram1D(lbl+"B tracks/dd0",50,0,1).fill(dd0);
+ aida.histogram1D(lbl+"B tracks/z0",50,-1,1).fill(z0);
+ aida.histogram1D(lbl+"B tracks/dz0",50,0,1).fill(dz0);
+ aida.histogram1D(lbl+"B tracks/d0 over dd0",200,-10,10).fill(d0/dd0);
+ aida.histogram1D(lbl+"B tracks/z0 over dz0",200,-10,10).fill(z0/dz0);
+ aida.histogram1D(lbl+"B tracks/IP3D",50,0,2).fill(IP3D);
+ aida.histogram1D(lbl+"B tracks/ID3D TP error",50,0,1).fill(IP3D_cverr);
+ aida.histogram1D(lbl+"B tracks/IP3D CLIC error",50,0,1).fill(IP3D_CLerr);
+ aida.histogram1D(lbl+"B tracks/IP3D MuCo error",50,0,1).fill(IP3D_MCerr);
+ aida.histogram1D(lbl+"B tracks/IP3D over IP3D TP error",100,0,20).fill(IP3D/IP3D_cverr);
+ aida.histogram2D(lbl+"B tracks/Ip3D over IP3Derr vs Origin",
+ 200,0,25,200,0,10).fill(IP3D/IP3D_cverr, rmag);
+ aida.histogram1D(lbl+"B tracks/IP3D over IP3D CLIC error",100,0,20).fill(IP3D/IP3D_CLerr);
+ aida.histogram1D(lbl+"B tracks/IP3D over IP3D MuCo error",100,0,20).fill(IP3D/IP3D_MCerr);
+ aida.histogram1D(lbl+"B tracks/IP3D TP error over IP3D CLIC error",50,0,2).fill(IP3D_cverr/IP3D_CLerr);
+ aida.histogram1D(lbl+"B tracks/IP3D TP error over IP3D MuCo error",50,0,2).fill(IP3D_cverr/IP3D_MCerr);
+ aida.histogram1D(lbl+"B tracks/IP3D MuCo error over IP3D CLIC error",50,1.4,1.8).fill(IP3D_MCerr/IP3D_CLerr);
+ aida.histogram1D(lbl+"B tracks/track momentum cos theta",50,0,1).fill(pCostheta);
+ aida.histogram1D(lbl+"B tracks/track momentum",50,0,50).fill(mcp.getMomentum().magnitude());
+ aida.histogram1D(lbl+"B tracks/track origin",100,0,200).fill(rmag);
+ if (IP3D/IP3D_cverr >= 3.0) {
+ tracks_accepted_TP++;
+ }
+ if (IP3D/IP3D_CLerr >= 3.0) {
+ tracks_accepted_CL++;
+ }
+ if (IP3D/IP3D_MCerr >= 3.0) {
+ tracks_accepted_MC++;
+ }
+ }
+ }
+ }
+ }
+ aida.cloud1D(lbl+"B tracks/tracks accepted TP").fill(tracks_accepted_TP);
+ aida.cloud1D(lbl+"B tracks/tracks accepted CL").fill(tracks_accepted_CL);
+ aida.cloud1D(lbl+"B tracks/tracks accepted MC").fill(tracks_accepted_MC);
+ aida.cloud1D(lbl+"B tracks/num tracks").fill(tracks);
+ aida.cloud1D(lbl+"B tracks/frac tracks accepted TP").fill((double)tracks_accepted_TP/(double)(tracks));
+ aida.cloud1D(lbl+"B tracks/frac tracks accepted CL").fill((double)tracks_accepted_CL/(double)(tracks));
+ aida.cloud1D(lbl+"B tracks/frac tracks accepted MC").fill((double)tracks_accepted_MC/(double)(tracks));
+ if (tracks_accepted_TP >= 2) {
+ num_Bs_accepted_TP++;
+ }
+ if (tracks_accepted_CL >= 2) {
+ num_Bs_accepted_CL++;
+ }
+ if (tracks_accepted_MC >= 2) {
+ num_Bs_accepted_MC++;
+ }
+ return (new int[]{ tracks_accepted_TP, tracks_accepted_CL, tracks_accepted_MC });
+ }
+
+ @Override
+ public void endOfData() {
+ System.out.println("Single-B acceptance TP: " + ((double)num_Bs_accepted_TP/(double)num_Bs));
+ System.out.println("Single-B acceptance CL: " + ((double)num_Bs_accepted_CL/(double)num_Bs));
+ System.out.println("Single-B acceptance MC: " + ((double)num_Bs_accepted_MC/(double)num_Bs));
+ System.out.println("Four-B acceptance TP: " + ((double)events_accepted_TP/(double)nevents));
+ System.out.println("Four-B acceptance CL: " + ((double)events_accepted_CL/(double)nevents));
+ System.out.println("Four-B acceptance MC: " + ((double)events_accepted_MC/(double)nevents));
+ }
+
+ public List<MCParticle> getXChildren(MCParticle parent, int gens) {
+ List<MCParticle> rv = new ArrayList();
+ for (MCParticle p:parent.getDaughters()) {
+ if ((p.getGeneratorStatus() == 1) && (p.getCharge() != 0))
+ rv.add(p);
+ else if (gens > 1) {
+ rv.addAll(getXChildren(p,gens-1));
+ }
+ }
+ return rv;
+ }
+ Track getTrack(MCParticle mcp, List<LCRelation> trackmap) {
+ Track rv = null;
+ for (LCRelation lc:trackmap) {
+ if (lc.getTo().equals(mcp)) {
+ rv = (Track) lc.getFrom();
+ break;
+ }
+ }
+ return rv;
+ }
+
+ double getIPError(MCParticle mcp, double a, double b) {
+ Hep3Vector p = mcp.getMomentum();
+ Hep3Vector xy = new BasicHep3Vector(1.0,1.0,0.0);
+ double pt = VecOp.dot(p,xy);
+ double theta = Math.acos(Math.abs(VecOp.cosTheta(p)));
+ aida.cloud1D("theta").fill(theta);
+ aida.cloud1D("cos theta").fill(Math.cos(theta));
+ theta = Math.pow(Math.sin(theta), 1.5);
+ double sigsq = (a*a + b*b/(pt*pt*theta*theta))*pt*pt;
+ return Math.sqrt(sigsq);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N visEnDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ visEnDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,65 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+//package org.lcsim.mcd.analysis;
+
+import hep.aida.IHistogram1D;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class visEnDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ IHistogram1D h_visEn;
+
+ double _nevents;
+ protected void startOfData() {
+ h_visEn = aida.histogram1D("Visible Energy", 20, 0, 20);
+ _nevents = 0.0;
+ }
+
+ public void process(EventHeader event) {
+ _nevents += 1.0;
+ List<MCParticle> allMC = event.get(MCParticle.class, "MCParticle");
+ List<MCParticle> allFS = new ArrayList<MCParticle>();
+ double totEn = 0.0;
+ for (MCParticle mcp:allMC) {
+ int pdgid = Math.abs(mcp.getPDGID());
+ if ( (mcp.getGeneratorStatus()==1)
+ && (pdgid!=12)
+ && (pdgid!=14)
+ && (pdgid!=16) ){
+ allFS.add(mcp);
+ totEn += mcp.getEnergy();
+ }
+ }
+
+ for (double i=0.0; i<21.0; i+=1.0) {
+ double angle = i*Math.PI*2.0/180.0;
+ double visEn = 0.0;
+ for (MCParticle mcp:allFS) {
+ aida.cloud1D("cosTheta").fill(Math.abs(VecOp.cosTheta(mcp.getMomentum())));
+ aida.cloud1D("cosTheta weight En").fill(Math.abs(VecOp.cosTheta(mcp.getMomentum())),mcp.getEnergy());
+ if (Math.abs(VecOp.cosTheta(mcp.getMomentum())) < Math.cos(angle)) {
+ visEn += mcp.getEnergy();
+ }
+ }
+ h_visEn.fill(i,visEn/totEn);
+ }
+ }
+
+ protected void endOfData(){
+ h_visEn.scale(1.0/_nevents);
+ }
+
+}
\ No newline at end of file
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N DRCorrectionDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DRCorrectionDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,15 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.mcd.analysis;
+
+/**
+ *
+ * @author aconway
+ */
+public class DRCorrectionDriver {
+
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N TimingAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TimingAnalysisDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,258 @@
+package org.lcsim.mcd.analysis;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.IHistogram3D;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * Exploratory processor
+ * @author aconway
+ */
+public class TimingAnalysisDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ double totEdep;
+ double totOpti;
+ //IHistogram2D rvzEdep;
+ IHistogram2D rvzOpti;
+ IHistogram1D fracEdep;
+ IHistogram1D fracOpti;
+ IHistogram1D fracTotal;
+
+ IHistogram3D trzeEdep;
+ IHistogram3D trzeOpti;
+ IHistogram2D rvzOrigins;
+
+ public void startOfData() {
+ totEdep=0;
+ totOpti=0;
+ fracEdep = aida.histogram1D("Edep/fraction energy deposited vs time",
+ 100,-10,90);
+ fracOpti = aida.histogram1D("Opti/fraction energy deposited vs time",
+ 100,-10,90);
+ fracTotal = aida.histogram1D("Total/fraction energy deposited vs time",
+ 100,-10,90);
+ IHistogram2D rvzEdep = aida.histogram2D("Edep/2D r v z",
+ 100,0,4000,
+ 100,-2000,2000);
+ rvzOpti = aida.histogram2D("Opti/2D r v z",
+ 100,1250,4000,
+ 100,-2000,2000);
+
+ trzeEdep = aida.histogram3D("Edep/3D t v r v z",
+ 100, 0.0, 1000.0, // t in ns
+ 100, 1250.0, 4000.0, // r in mm
+ 100, -2000.0, 2000.0); // z in mm
+ trzeOpti = aida.histogram3D("Opti/3D t v r v z",
+ 100, 0.0, 1000.0, // t in ns
+ 100, 1250.0, 4000.0, // r in mm
+ 100, -2000.0, 2000.0); // z in mm
+
+ rvzOrigins = aida.histogram2D("MCP/r z",
+ 100,-8000,8000,
+ 100,0,2200);
+ }
+ public void process(EventHeader event) {
+ List<CalorimeterHit> EdepHitList = new ArrayList<CalorimeterHit>();
+ EdepHitList.addAll(event.get(CalorimeterHit.class, "ECalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(CalorimeterHit.class, "HCalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(CalorimeterHit.class, "ECalEndcapEdepHits"));
+ EdepHitList.addAll(event.get(CalorimeterHit.class, "HCalEndcapEdepHits"));
+
+ List<CalorimeterHit> OptiHitList = new ArrayList<CalorimeterHit>();
+ OptiHitList.addAll(event.get(CalorimeterHit.class, "ECalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(CalorimeterHit.class, "HCalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(CalorimeterHit.class, "ECalEndcapOptiHits"));
+ OptiHitList.addAll(event.get(CalorimeterHit.class, "HCalEndcapOptiHits"));
+ //List<List<SimCalorimeterHit>> calHitList = event.get(SimCalorimeterHit.class);
+ List<SimTrackerHit> trkrHitList = new ArrayList<SimTrackerHit>();
+ trkrHitList.addAll(event.get(SimTrackerHit.class, "TkrBarrHits"));
+ trkrHitList.addAll(event.get(SimTrackerHit.class, "TkrEndcapHits"));
+ trkrHitList.addAll(event.get(SimTrackerHit.class, "VtxBarrHits"));
+ trkrHitList.addAll(event.get(SimTrackerHit.class, "VtxEndcapHits"));
+ List<MCParticle> mcpList = event.get(MCParticle.class, "MCParticle");
+
+
+ for (CalorimeterHit hit : EdepHitList) {
+ totEdep += hit.getRawEnergy();
+ double tofl = hit.getTime() - hit.getPositionVec().magnitude()/299.79;
+ fracEdep.fill(tofl, hit.getRawEnergy());
+ fracTotal.fill(tofl, hit.getRawEnergy());
+ }
+ for (CalorimeterHit hit : OptiHitList) {
+ totOpti += hit.getRawEnergy();
+ double tofl = hit.getTime() - hit.getPositionVec().magnitude()/299.79;
+ fracOpti.fill(tofl, hit.getRawEnergy());
+ fracTotal.fill(tofl, hit.getRawEnergy());
+ }
+
+ makeHists(EdepHitList, "Edep/");
+ makeHists(OptiHitList, "Opti/");
+// trackerHists(trkrHitList);
+// MCPHists(mcpList, "MCP/");
+ }
+
+ @Override
+ public void endOfData() {
+ fracEdep.scale(1.0/fracEdep.sumAllBinHeights());
+ fracOpti.scale(1.0/fracOpti.sumAllBinHeights());
+ fracTotal.scale(1.0/fracTotal.sumAllBinHeights());
+
+
+ for (double t = 0.0; t < 10; t += 1.0) {
+ double sumEdep = 0.0;
+ double sumOpti = 0.0;
+ double sumTotal = 0.0;
+ for (int i=fracEdep.coordToIndex(-1.0); i<=fracEdep.coordToIndex(t); i++) {
+ sumEdep += fracEdep.binHeight(i);
+ }
+ for (int i=fracEdep.coordToIndex(-1.0); i<=fracOpti.coordToIndex(t); i++) {
+ sumOpti += fracOpti.binHeight(i);
+ }
+ for (int i=fracTotal.coordToIndex(-1.0); i<=fracTotal.coordToIndex(t); i++) {
+ sumTotal += fracTotal.binHeight(i);
+ }
+ System.out.println("-1 to "+String.valueOf(t)+"ns:");
+ System.out.println("Edep: "+String.valueOf(sumEdep));
+ System.out.println("Opti: "+String.valueOf(sumOpti));
+ System.out.println("Total: "+String.valueOf(sumTotal));
+ }
+ //save3DHist(trzeEdep,"3DHist-trz-EdepBarr");
+ //save3DHist(trzeOpti,"3DHist-trz-OptiBarr");
+ }
+
+ private void MCPHists (List<MCParticle> mcpList, String prefix) {
+ for (MCParticle mcp : mcpList) {
+ Hep3Vector start = mcp.getOrigin();
+ Hep3Vector pvec = mcp.getMomentum();
+ double r = Math.sqrt(start.x()*start.x() + start.y()*start.y());
+ double z = start.z();
+ double phi = VecOp.phi(start);
+ double t = mcp.getProductionTime();
+ double p = mcp.getMomentum().magnitude();
+ double pr = Math.sqrt(pvec.x()*pvec.x()+pvec.y()*pvec.y());
+
+ Hep3Vector stop = mcp.getEndPoint();
+ double rstop = Math.sqrt(stop.x()*stop.x() + stop.y()+stop.y());
+ double zstop = stop.z();
+
+ rvzOrigins.fill(z,r,pr/p);
+ aida.cloud2D(prefix+"r phi").fill(r,phi,p);
+ aida.cloud2D(prefix+"z phi").fill(z,phi, p);
+ aida.cloud1D(prefix+"t").fill(t,p);
+ aida.cloud2D(prefix+"r v z stop").fill(z,r,p);
+ int pdgid = Math.abs(mcp.getPDGID());
+ String name;
+ if (pdgid == 22) {
+ name = "gamma/";
+ } else if (pdgid == 11) {
+ name = "electron/";
+ } else if (pdgid == 2112) {
+ name = "neutron/";
+ } else {
+ name = "other/";
+ }
+ aida.cloud1D(prefix+name+"momentum").fill(p);
+ if (t<100) {
+ aida.cloud1D(prefix+name+"time").fill(t,p);
+ if (mcp.getCharge() != 0) {
+ aida.cloud1D(prefix+"charged/time").fill(t,p);
+ } else {
+ aida.cloud1D(prefix+"uncharged/time").fill(t,p);
+ }
+ }
+
+ }
+ }
+
+ private void makeHists (List<CalorimeterHit> hitList, String prefix) {
+ for (CalorimeterHit hit : hitList) {
+ Hep3Vector pos = hit.getPositionVec();
+ double r = Math.sqrt(pos.x()*pos.x() + pos.y()*pos.y());
+ double z = pos.z();
+ double t = hit.getTime();
+ double tof = t - pos.magnitude()/299.79;
+ double e = hit.getRawEnergy();
+ aida.histogram2D(prefix + "2D r v z").fill(r, z, e);
+ aida.histogram3D(prefix + "3D t v r v z").fill(t, r, z, e);
+ timingHists(t,z,r,e,prefix+"/abs t/");
+ timingHists(tof,z,r,e,prefix+"/tof/");
+ }
+ }
+ private void trackerHists(List<SimTrackerHit> hitList) {
+ for (SimTrackerHit hit : hitList) {
+ Hep3Vector pos = hit.getPositionVec();
+ double r = Math.sqrt(pos.x()*pos.x() + pos.y()*pos.y());
+ double z = pos.z();
+ double t = hit.getTime();
+ double tof = t - pos.magnitude()/299.79;
+ double[] pvec = hit.getMomentum();
+ double p = Math.sqrt(pvec[0]*pvec[0]+pvec[1]*pvec[1]+pvec[2]*pvec[2]);
+ timingHists(t, z, r, p, "tracking/abs t/");
+ timingHists(tof, z, r, p, "tracking/tof/");
+ }
+ }
+
+ private void timingHists(double t, double z, double r, double e, String prefix) {
+ if (Math.abs(t) < 10.0) {
+ double timebinwidth = 1;
+ double binnedtime = t - t % timebinwidth;
+ int bin = (int) (binnedtime / timebinwidth);
+ aida.cloud2D(prefix + "r v z/" + String.valueOf(bin * timebinwidth) + "ns").fill(z, r, e);
+ aida.cloud1D(prefix + "t/" + String.valueOf(bin * timebinwidth) + "ns").fill(t,e);
+ } else if (Math.abs(t) < 20.0) {
+ double timebinwidth = 2.0;
+ double binnedtime = t - t % timebinwidth;
+ int bin = (int) (binnedtime / timebinwidth);
+ aida.cloud2D(prefix + "r v z/" + String.valueOf(bin * timebinwidth) + "ns").fill(z, r, e);
+ aida.cloud1D(prefix + "t/" + String.valueOf(bin * timebinwidth) + "ns").fill(t,e);
+ } else if (Math.abs(t) < 100.0) {
+ double timebinwidth = 20.0;
+ double binnedtime = t - t % timebinwidth;
+ int bin = (int) (binnedtime / timebinwidth);
+ aida.cloud2D(prefix + "r v z/" + String.valueOf(bin * timebinwidth) + "ns").fill(z, r, e);
+ aida.cloud1D(prefix + "t/" + String.valueOf(bin * timebinwidth) + "ns").fill(t,e);
+ }
+ }
+
+ private void save3DHist(IHistogram3D hist, String filename) {
+ try {
+ PrintWriter writer = new PrintWriter("/home/aconway/fermi/3Dhists/"+filename+".txt", "US-ASCII");
+ int nbinsx = hist.xAxis().bins();
+ int nbinsy = hist.yAxis().bins();
+ int nbinsz = hist.zAxis().bins();
+ for (int i = 0; i < nbinsx; i++) {
+ double t = hist.xAxis().binLowerEdge(i);
+ for (int j = 0; j < nbinsy; j++) {
+ double r = hist.yAxis().binLowerEdge(j);
+ for (int k = 0; k < nbinsz; k++) {
+ double z = hist.zAxis().binLowerEdge(k);
+ double h = hist.binHeight(i, j, k);
+ writer.println(
+ String.valueOf(t) + " " +
+ String.valueOf(r) + " " +
+ String.valueOf(z) + " " +
+ String.valueOf(h));
+ }
+ }
+ }
+ writer.close();
+ System.out.println("Successfully saved"+filename+".txt!");
+ } catch (IOException ex) {
+ }
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N trilinearHiggsDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ trilinearHiggsDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,644 @@
+//package org.lcsim.mcd.analysis;
+
+import hep.aida.IHistogram1D;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.mc.fast.MCFast;
+import java.util.List;
+import java.util.ArrayList;
+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.Hep3Vector;
+import hep.physics.vec.VecOp;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.Track;
+import org.lcsim.mc.fast.tracking.DocaTrackParameters;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+/**
+ * Simple driver to do generator-level study on mu+ mu- -> h h vm vm~ via W-fusion
+ * Events generated in MadGraph with process mu+ mu- > h vm vm~ / z , h > h h / z
+ * @author aconway
+ */
+public class trilinearHiggsDriver extends Driver {
+
+ public trilinearHiggsDriver()
+ {
+ add(new MCFast());
+ //add(new MCFastTracking(false,true));
+ }
+ private AIDA aida = AIDA.defaultInstance();
+ protected String mcParticlesSkimmedName = "MCParticle";
+ List<LCRelation> tracks2MCP;
+ List<Track> all_tracks;
+ IHistogram1D cone;
+ IHistogram1D conexef;
+ IHistogram1D costhetastar;
+ IHistogram1D fourBacceptanceCLIC;
+ IHistogram1D fourBacceptanceMuC;
+
+ double MuC_cone = 10.0*Math.PI/180.0;
+ double CLIC_cone = 20.0*Math.PI/180.0;
+
+ double evt_wt = 1.732;
+
+ int numPassMuCTest = 0;
+ int numPassCLICTest = 0;
+ int[] passConeTest = new int[21];
+ int[] passMuCoTest = new int[21];
+ int[] passCLICTest = new int[21];
+ int bbbarevts = 0;
+
+ protected void startOfData() {
+ cone = aida.histogram1D("cone angle",40,0,40);
+ conexef = aida.histogram1D("cone angle weighted by fraction of visible energy",40,0,40);
+ costhetastar = aida.histogram1D("Cosine of Decay Angle",20,0.0,1.0);
+ }
+
+ public void process(EventHeader event) {
+ this.processChildren(event);
+ List<MCParticle> particles = event.get(MCParticle.class,mcParticlesSkimmedName);
+ tracks2MCP = event.get(LCRelation.class,"TracksToMCP");
+ all_tracks = event.get(Track.class,"Tracks");
+ List<MCParticle> children = new ArrayList<MCParticle>();
+
+ Hep3Vector xy = new BasicHep3Vector(1.0,1.0,0.0);
+ List<MCParticle> higgses = new ArrayList<MCParticle>();
+ List<MCParticle> allFSPs = new ArrayList<MCParticle>();
+ List<MCParticle> allFSPs_nucut = new ArrayList<MCParticle>();
+ List<MCParticle> allBmesons = new ArrayList<MCParticle>();
+ List<MCParticle> allBquarks = new ArrayList<MCParticle>();
+
+ for (Track trk:all_tracks) {
+ double [] tpms = trk.getTrackParameters();
+ double dd0 = trk.getErrorMatrix().diagonal(0);
+ aida.histogram1D("tracks/d0",50,-1.0,1.0).fill(tpms[0]);
+ aida.histogram1D("tracks/d0 over error",200,-10,10).fill(tpms[0]/Math.sqrt(dd0));
+ aida.cloud1D("tracks/phi0").fill(tpms[1]);
+ aida.cloud1D("tracks/omega").fill(tpms[2]);
+ aida.cloud1D("tracks/z0").fill(tpms[3]);
+ aida.cloud1D("tracks/s").fill(tpms[4]);
+ aida.histogram1D("tracks/costheta",100,0,1).fill(Math.abs(VecOp.cosTheta(new BasicHep3Vector(trk.getMomentum()))));
+
+ MCParticle mcp = null;
+ for (LCRelation lcr:tracks2MCP) {
+ if (trk.equals(lcr.getFrom())) {
+ mcp = (MCParticle) lcr.getTo();
+ }
+ }
+ if (mcp != null) {
+ if (mcp.getOrigin().magnitude() == 0.0) {
+ aida.histogram1D("origin tracks/d0",50,-1,1).fill(tpms[0]);
+ aida.histogram1D("origin tracks/d0 over error",200,-10,10).fill(tpms[0]/Math.sqrt(dd0));
+ }
+ }
+ }
+
+ for (MCParticle p:particles) {
+ int pdgid = Math.abs(p.getPDGID());
+ // get Higgs particles in documentation status
+ // MCParticle list does not show the original Higgs
+ //if ((p.getPDGID()==25) ) {//&& (p.getGeneratorStatus()==3)) {
+ if ((p.getPDGID()==25) && (p.getGeneratorStatus()==3)) {
+ higgses.add(p);
+ // Basic plots
+ aida.cloud1D("higgs/masses").fill(p.getMass());
+ aida.cloud1D("higgs/momentum").fill(p.getMomentum().magnitude());
+ aida.cloud1D("higgs/transverse momentum").fill(VecOp.dot(p.getMomentum(), xy));
+ aida.cloud1D("higgs/cos theta (lab)").fill(VecOp.cosTheta(p.getMomentum()));
+ } else
+ // add all the FS particles except for the IS neutrinos
+ // Have to check the muon neutrinos if their parent is same
+ if ( (p.getGeneratorStatus()==1) ) {
+ if (pdgid!=14) {
+ allFSPs.add(p);
+ // if it is a muon neutrino, check if its parent is a muon neutrino
+ // only the IS neutrinos will have muon neutrino parents.
+ } else if (Math.abs(p.getParents().get(0).getPDGID())!=14) {
+ allFSPs.add(p);
+ }
+ // cut all neutrinos
+ if ( (pdgid!=12) && (pdgid!=14) && (pdgid!=16) ){
+ allFSPs_nucut.add(p);
+ }
+ }
+ /*
+ * Here we find the mesons from b-decays and their children and
+ * make some plots.
+ */
+ else if ( (p.getGeneratorStatus()==2)
+ && (
+ ( (pdgid%10000 > 500) && (pdgid%10000 < 600) )
+ ||
+ ( (pdgid > 5000) && (pdgid < 6000) )
+// pdgid == 511
+ )
+ && (p.getOrigin().magnitude() == 0.0)
+ && (p.getEndPoint().magnitude() != 0.0) )
+ {
+ aida.cloud1D("newBs/B energy").fill(p.getEnergy());
+ aida.cloud1D("newBs/B cosTheta").fill(Math.abs(VecOp.cosTheta(p.getMomentum())));
+ aida.cloud1D("newBs/B distance").fill(p.getEndPoint().magnitude());
+ List<MCParticle> ch = p.getDaughters();
+ aida.cloud1D("newBs/B decay multiplicity").fill(ch.size());
+ List<MCParticle> allCh = getAllChildren(p);
+ aida.cloud1D("newBs/all B children").fill(allCh.size());
+ double endPointMag = p.getEndPoint().magnitude();
+ boolean addIt = true;
+ for (MCParticle mcp:allCh) {
+ if ( (mcp.getOrigin().magnitude() == endPointMag) && (mcp.getGeneratorStatus()==1)) {
+ aida.cloud1D("newBs/B decay pdgids").fill(Math.abs(mcp.getPDGID()));
+
+ }
+ aida.cloud1D("newBs/all B decay pdgids").fill(Math.abs(mcp.getPDGID()));
+ }
+ allBmesons.add(p);
+ }
+ }
+ IPCutPlots(allBmesons,"IP Cut Plots/");
+ aida.cloud1D("num B-mesons").fill(allBmesons.size());
+ if (allBmesons.size()==4) {
+ bbbarevts++;
+ }
+
+ }
+
+ void IPCutPlots(List<MCParticle> bMesons, String lbl) {
+ double mint = 0.0;
+ double maxt = 20.0;
+ double nsteps = 20.0;
+ double step = (maxt-mint)/nsteps;
+ for (double i=mint; i<=maxt; i+=step) {
+ int numvisCLIC = 0;
+ int numvisMuC = 0;
+ int numvisBs = 0;
+ double toten = 0.0;
+ double visen = 0.0;
+ double enlep = 0.0;
+ double enchr = 0.0;
+ double cosCone = Math.abs(Math.cos(i*Math.PI/180.0));
+ for (MCParticle B:bMesons) {
+ int numlepCLIC = 0;
+ int numlepMuC = 0;
+ int numchrCLIC = 0;
+ int numchrMuC = 0;
+ Hep3Vector bVtex = B.getEndPoint();
+ // Only bother if B doesn't hit cone
+ if (Math.abs(VecOp.cosTheta(B.getEndPoint())) < cosCone) {
+ aida.cloud1D(lbl+"Single B Acceptance (cone cut)").fill(i);
+ aida.histogram1D(lbl+"B meson ctau",100,0,2).fill(B.getEndPoint().magnitude()*B.getMass()/B.getMomentum().magnitude());
+ aida.cloud1D(lbl+"B meson ctau cloud").fill(B.getEndPoint().magnitude()*B.getMass()/B.getMomentum().magnitude());
+ aida.cloud1D(lbl+"B meson p over m").fill(B.getMomentum().magnitude()/B.getMass());
+ aida.cloud1D(lbl+"B meson vertex").fill(B.getEndPoint().magnitude());
+ numvisBs++;
+ List<MCParticle> allCh = getXChildren(B,3);
+ //List<MCParticle> allCh = getAllChildren(B);
+ //List<MCParticle> allCh = B.getDaughters();
+ for (MCParticle ch:allCh) {
+ int pdgid = Math.abs(ch.getPDGID());
+ double cosCh = Math.abs(VecOp.cosTheta(ch.getMomentum()));
+ Hep3Vector orig_xy = new BasicHep3Vector(ch.getOriginX(),ch.getOriginY(),0.0);
+ Hep3Vector orig_z = new BasicHep3Vector(0.0,0.0,ch.getOriginZ());
+ // requirements for good track:
+ if ( (ch.getGeneratorStatus() == 1) &&
+ (cosCh < cosCone) &&
+ (orig_xy.magnitude() < 220.0) &&
+ (orig_z.magnitude() < 300.0) &&
+ (Math.abs(VecOp.cosTheta(ch.getOrigin())) < cosCone) &&
+ (ch.getCharge() != 0) ) {
+ Track trk = getTrack(ch,tracks2MCP);
+ if (trk != null) {
+ double d0 = trk.getTrackParameter(0);
+ double dd0 = trk.getErrorMatrix().diagonal(0);
+ double z0 = trk.getTrackParameter(3);
+ double IP = Math.sqrt(d0*d0 + z0*z0);
+ // double IP = getImpactParameter(ch);
+ double IPeCLIC = getIPError(ch,0.005,0.015);
+ //double IPeMuC = getIPError(ch,0.008,0.025);
+ double IPeMuC = Math.sqrt(dd0 + trk.getErrorMatrix().diagonal(3));
+ IPeMuC = Math.sqrt(IPeMuC*IPeMuC + 0.006*0.006);
+ if (i==1.0) {
+ aida.cloud1D(lbl+"num FS Children per B").fill(allCh.size());
+ aida.histogram1D(lbl+"d0 of B tracks",50,-1,1).fill(d0);
+ aida.histogram1D(lbl+"d0 over dd0",200,-10,10).fill(d0/Math.sqrt(dd0));
+ aida.cloud2D(lbl+"IP vs IPeCLIC").fill(IP, IPeCLIC);
+ aida.cloud2D(lbl+"IP vs IPeMuC").fill(IP,IPeMuC);
+ aida.cloud1D(lbl+"IP over IPeCLIC").fill(IP/IPeCLIC);
+ aida.histogram1D(lbl+"IP over IPeCLIC hist",
+ 20,0.0,20.0).fill(IP/IPeCLIC);
+ aida.cloud1D(lbl+"track origin").fill(ch.getOrigin().magnitude());
+ aida.histogram1D(lbl+"IP",50,0,1).fill(IP);
+ aida.histogram1D(lbl+"track origin hist",100,0.0,200.0).fill(ch.getOrigin().magnitude());
+ aida.cloud1D(lbl+"Impact parameter").fill(IP);
+ aida.cloud1D(lbl+"Impact parameter error CLIC").fill(IPeCLIC);
+ aida.cloud1D(lbl+"Impact parameter error MuC").fill(IPeMuC);
+ if (IP > 200) {
+ System.out.println("IP > 200: PDGID:" + pdgid + ", Parent PDGID:" + ch.getParents().get(0).getPDGID());
+ }
+ }
+ if (IP > 3.0*IPeCLIC) {
+ aida.cloud1D(lbl+"Track acceptance at CLIC").fill(i);
+ if ( (pdgid == 11) || (pdgid == 13) || (pdgid==15) ) {
+ numchrCLIC++;
+ aida.cloud1D(lbl+"Energy of Leptons").fill(ch.getEnergy());
+ aida.cloud1D(lbl+"|p| of Charged tracks").fill(ch.getMomentum().magnitude());
+ } else if ( (pdgid!=14) && (pdgid!=12) && (pdgid!=16) ){
+ numchrCLIC++;
+ aida.cloud1D(lbl+"Energy of Charged tracks").fill(ch.getEnergy());
+ aida.cloud1D(lbl+"|p| of Charged tracks").fill(ch.getMomentum().magnitude());
+ }
+ }
+ if (IP > 3.0*IPeMuC) {
+ aida.cloud1D(lbl+"Track acceptance at MuC").fill(i);
+ if ( (pdgid == 11) || (pdgid == 13) || (pdgid==15) ) {
+ numchrMuC++;
+ } else if ( (pdgid!=14) && (pdgid!=12) && (pdgid!=16) ){
+ numchrMuC++;
+ }
+ }
+ }
+ }
+ }
+ if ((i==1) && (allCh.size() != 0))
+ aida.cloud1D(lbl+"Frac of Tracks accepted per B (MuC)").fill((double)numchrMuC/(double)allCh.size());
+ }
+ if ( (numlepCLIC >= 1) || (numchrCLIC >= 2) ) {
+ numvisCLIC++;
+ aida.cloud1D(lbl+"Single B Acceptance (CLIC)").fill(i);
+ }
+ if ( (numlepMuC >= 1) || (numchrMuC >= 2) ) {
+ numvisMuC++;
+ aida.cloud1D(lbl+"Single B Acceptance (MuC)").fill(i);
+ }
+ if (i == 1) {
+ aida.cloud1D(lbl+"num Tracks per B (MuC)").fill(numvisMuC);
+ }
+ }
+ aida.cloud1D(lbl+"visible B's").fill(numvisBs);
+ if (numvisBs == 4) {
+ aida.cloud1D(lbl+"4-B acceptance").fill(i);
+
+ }
+ aida.cloud1D(lbl+"numvisCLIC hist").fill(numvisCLIC);
+ if (numvisCLIC == 4) {
+ passCLICTest[(int)i] += 1;
+ aida.cloud1D(lbl+"CLIC 4-b acceptance").fill(i);
+ }
+ if (numvisMuC == 4) {
+ passMuCoTest[(int)i] += 1;
+ aida.cloud1D(lbl+"MuC 4-b acceptance").fill(i);
+ }
+ }
+ }
+
+ double getImpactParameter(MCParticle p) {
+ double bField = 5.0;
+ DocaTrackParameters dtp = new DocaTrackParameters(bField,p.getMomentum(),p.getOrigin(),p.getCharge());
+ /**
+ * The track parameters for LCD are defined as follows
+ * <table>
+ * <tr><th>Index</th><th>Meaning</th></tr>
+ * <tr><td> 0 </td><td> d0 = XY impact parameter </td><tr>
+ * <tr><td> 1 </td><td> phi0 </td><tr> </td><tr>
+ * <tr><td> 2 </td><td> omega = 1/curv.radius (negative for negative tracks) </td><tr>
+ * <tr><td> 3 </td><td> z0 = z of track (z impact parameter) </td><tr>
+ * <tr><td> 4 </td><td> s = tan lambda </td><tr>
+ * </table>
+ * **/
+ double[] tps = dtp.getTrackParameters();
+ //return Math.sqrt(tps[0]*tps[0] + tps[3]*tps[3]);
+ return tps[0];
+ }
+ double getIPError(MCParticle mcp, double a, double b) {
+ Hep3Vector p = mcp.getMomentum();
+ double pt = p.magnitude();
+ double theta = Math.acos(Math.abs(VecOp.cosTheta(mcp.getMomentum())));
+ aida.cloud1D("theta").fill(theta);
+ theta = Math.pow(Math.sin(theta), 1.5);
+ double sigsq = a*a + b*b/(pt*pt*theta*theta);
+ return Math.sqrt(sigsq);
+ }
+
+ Track getTrack(MCParticle mcp, List<LCRelation> trackmap) {
+ Track rv = null;
+ for (LCRelation lc:trackmap) {
+ if (lc.getTo().equals(mcp)) {
+ rv = (Track) lc.getFrom();
+ break;
+ }
+ }
+ return rv;
+ }
+ /*
+ * Plot visible energy as function of cone angle;
+ * assume cone is straight cone to IP
+ * lbl string is plot subfolder, eg. "bbbar"
+ */
+// void coneAnglePlots(List<MCParticle> mcps, String lbl) {
+//// double mint = 0.0;
+//// double maxt = Math.PI/4.0;
+//// double nsteps = 50.0;
+//// double step = (maxt-mint)/nsteps;
+// bbbarevts++;
+// double mint = 0.0;
+// double maxt = 20.0;
+// double nsteps = 20.0;
+// double step = (maxt-mint)/nsteps;
+// for (double i=mint; i<=maxt; i+=step) {
+// double visen = 0.0;
+// double toten = 0.0;
+// int numvis = 0;
+// for (MCParticle p:mcps) {
+// if (VecOp.cosTheta(p.getMomentum()) < Math.cos(i*Math.PI/180.0)) {
+// int pdgid = Math.abs(p.getPDGID());
+// if ( (pdgid!=14) && (pdgid!=12) && (pdgid!=16) ){
+// if (childCut(p,i)) {
+// visen += p.getEnergy();
+// numvis++;
+// }
+// }
+// }
+// toten += p.getEnergy();
+// }
+// aida.histogram1D(lbl+"/number of b's visible",20,0,20).fill(numvis);
+// if (numvis==4) {
+// aida.histogram1D(lbl+"/events with all b's",20,0,20).fill(i);
+// passConeTest[(int)i] += 1;
+// }
+// aida.cloud1D(lbl+"/weighted").fill(i, visen/toten);
+// aida.cloud1D(lbl+"/unweighted").fill(i);
+// aida.cloud2D(lbl+"/visible energy").fill(i, visen);
+// aida.cloud2D(lbl+"/total energy").fill(i,toten);
+// aida.cloud1D(lbl+"/fraction visible").fill(i, visen/toten);
+// aida.cloud2D(lbl+"/fraction visible (deg)").fill(i*180.0/Math.PI, visen/toten);
+// aida.histogram2D(lbl+"/fraction visible (deg) hist",
+// (int) nsteps,mint,maxt,
+// 50,0,1.001
+// ).fill(i, visen/toten);
+// //aida.profile1D(lbl+"/fraction vis en profile").fill(i,visen/toten);
+// conexef.fill(i,visen/toten);
+// cone.fill(i);
+// }
+//
+// }
+
+
+ protected void endOfData() {
+ conexef.scale(1.0/10000.0);
+// System.out.println(Integer.toString(numPassMuCTest));
+// System.out.println(Integer.toString(numPassCLICTest));
+ System.out.println("CLIC:");
+ for(int i=0; i<passCLICTest.length; i++) {
+ double frac = ((double)passCLICTest[i])/((double)bbbarevts);
+ System.out.print(Double.toString(frac)+",");
+ }
+ System.out.println();
+ System.out.println("MuCo:");
+ for(int i=0; i<passMuCoTest.length; i++) {
+ double frac = ((double)passMuCoTest[i])/((double)bbbarevts);
+ System.out.print(Double.toString(frac)+",");
+ }
+ costhetastar.scale(evt_wt/costhetastar.entries());
+ aida.cloud1D("IP Cut Plots/"+"Single B Acceptance (MuC)").scale(1.0/40000.0);
+ }
+
+ /*
+ * Helper functions
+ */
+ private boolean childCut(MCParticle b, double angle) {
+ boolean rv = false;
+ int nlep = 0;
+ int ntrk = 0;
+ List<MCParticle> daughters = getAllFSChildren(b);
+ for (MCParticle p:daughters) {
+ if (Math.abs(VecOp.cosTheta(p.getMomentum())) < Math.cos(angle*Math.PI/180.0)) {
+ int pdgid = Math.abs(p.getPDGID());
+ if ( (pdgid == 11) || (pdgid == 13) || (pdgid==15) ) {
+ nlep++;
+ } else if (Math.abs(p.getCharge())==1.0) {
+ ntrk++;
+ }
+ }
+ }
+ rv = ( (nlep>=1) || (ntrk>=2));
+ return rv;
+ }
+ public List<MCParticle> getXChildren(MCParticle parent, int gens) {
+ List<MCParticle> rv = new ArrayList();
+ for (MCParticle p:parent.getDaughters()) {
+ if ((p.getGeneratorStatus() == 1) && (p.getCharge() != 0))
+ rv.add(p);
+ else if (gens > 1) {
+ rv.addAll(getXChildren(p,gens-1));
+ }
+ }
+ return rv;
+ }
+ public List<MCParticle> getAllFSChildren(MCParticle parent) {
+ List<MCParticle> AllChildren = new ArrayList<MCParticle>();
+ List<MCParticle> children = parent.getDaughters();
+ for (MCParticle child:children) {
+ if (child.getGeneratorStatus() == 1) {
+ AllChildren.add(child);
+ }
+ else {
+ AllChildren.addAll(getAllFSChildren(child));
+ }
+ }
+ return AllChildren;
+ }
+ public List<MCParticle> getAllChildren(MCParticle parent) {
+ List<MCParticle> AllChildren = new ArrayList<MCParticle>();
+ List<MCParticle> children = parent.getDaughters();
+ if (children.size() == 0) {
+ return AllChildren;
+ } else {
+ for (MCParticle child:children) {
+ AllChildren.add(child);
+ AllChildren.addAll(getAllFSChildren(child));
+ }
+ }
+ return AllChildren;
+ }
+ /*
+ * Basic Higgs plots
+ */
+ /*
+ MCParticle h1 = higgses.get(0);
+ MCParticle h2 = higgses.get(1);
+ Hep3Vector p1 = h1.getMomentum();
+ Hep3Vector p2 = h2.getMomentum();
+ List<Hep3Vector> hs = new ArrayList<Hep3Vector>();
+ hs.add(p1);
+ hs.add(p2);
+ Hep3Vector cm = VecOp.centerOfMass(hs);
+ aida.cloud1D("com/momentum magnitude").fill(cm.magnitude());
+ aida.cloud1D("com/momentum cos theta").fill(VecOp.cosTheta(cm));
+ aida.cloud1D("com/momentum theta").fill(Math.acos(VecOp.cosTheta(cm)));
+ aida.cloud1D("com/transverse momentum").fill(VecOp.dot(cm, xy));
+ aida.cloud1D("higgs/cos opening angle (lab)").fill(VecOp.dot(p1,p2)/(p1.magnitude()*p2.magnitude()));
+ aida.cloud1D("higgs/opening angle (lab)").fill(Math.acos(VecOp.dot(p1,p2)/(p1.magnitude()*p2.magnitude())));
+ Hep3Vector diff = VecOp.sub(p1, p2);
+ double hinvmass = Math.sqrt( (h1.getEnergy()+h2.getEnergy())*(h1.getEnergy()+h2.getEnergy()) - diff.magnitudeSquared());
+ aida.cloud1D("higgs/two-higgs invariant mass").fill(hinvmass);
+ aida.cloud1D("higgs/decay angle (com)").fill(Math.acos(VecOp.dot(diff,cm)/(diff.magnitude()*cm.magnitude())));
+ aida.cloud1D("higgs/cos decay angle (com)").fill(VecOp.dot(diff,cm)/(diff.magnitude()*cm.magnitude()));
+ costhetastar.fill(Math.abs(VecOp.dot(diff,cm)/(diff.magnitude()*cm.magnitude())));
+ */
+ /*
+ * Higgs-decay determination and plots
+ */
+ /*
+ // Get Higgs decay products
+ boolean first_bbdecay = false;
+ List<MCParticle> h1decay = new ArrayList<MCParticle>();
+ for (MCParticle p:h1.getDaughters()) {
+ if ( (Math.abs(p.getPDGID())!=14) && (p.getPDGID()!=25) ) {
+ h1decay.add(p);
+ aida.cloud1D("decay products/PDGID").fill(Math.abs(p.getPDGID()));
+ aida.cloud1D("decay products/momentum").fill(p.getMomentum().magnitude());
+ aida.cloud1D("decay products/transverse momentum").fill(VecOp.dot(p.getMomentum(), xy));
+ aida.cloud1D("decay products/cos theta").fill(VecOp.cosTheta(p.getMomentum()));
+ aida.cloud1D("decay products/theta").fill(Math.acos(VecOp.cosTheta(p.getMomentum())));
+ if (Math.abs(p.getPDGID())==5) {
+ first_bbdecay = true;
+ }
+ }
+ }
+ boolean second_bbdecay = false;
+ List<MCParticle> h2decay = new ArrayList<MCParticle>();
+ for (MCParticle p:h2.getDaughters()) {
+ if ( (Math.abs(p.getPDGID())!=14) && (p.getPDGID()!=25) ) {
+ h2decay.add(p);
+ aida.cloud1D("decay products/PDGID").fill(Math.abs(p.getPDGID()));
+ aida.cloud1D("decay products/momentum").fill(p.getMomentum().magnitude());
+ aida.cloud1D("decay products/transverse momentum").fill(VecOp.dot(p.getMomentum(), xy));
+ aida.cloud1D("decay products/cos theta").fill(VecOp.cosTheta(p.getMomentum()));
+ aida.cloud1D("decay products/theta").fill(Math.acos(VecOp.cosTheta(p.getMomentum())));
+ if (Math.abs(p.getPDGID())==5) {
+ second_bbdecay = true;
+ }
+ }
+ }
+ boolean four_bs = first_bbdecay && second_bbdecay;
+ if (four_bs) {
+ }
+ // True if all b's have either 1+ lepton or 2+ tracks outside cone.
+ boolean decaysInMuCCone = true;
+ boolean decaysInCLICCone = true;
+ if (false) {
+ List<MCParticle> bs = new ArrayList<MCParticle>();
+ bs.add(0,h1decay.get(0));
+ bs.add(1,h1decay.get(1));
+ bs.add(2,h2decay.get(0));
+ bs.add(3,h2decay.get(1));
+ double totEn = 0.0;
+ double totEnCLIC = 0.0;
+ double totEnMuC = 0.0;
+ for(MCParticle b:bs) {
+ boolean thisInMuCCone = false;
+ boolean thisInCLICCone = false;
+ int leps_muc = 0;
+ int leps_clic = 0;
+ int trks_muc = 0;
+ int trks_clic = 0;
+ Hep3Vector pb = b.getMomentum();
+ aida.cloud1D("bs/cos theta (lab)").fill(VecOp.cosTheta(pb));
+ aida.cloud1D("bs/momentum ").fill(pb.magnitude());
+ for (MCParticle ch:getAllFSChildren(b)) {
+ int pdgid = Math.abs(ch.getPDGID());
+ Hep3Vector pch = ch.getMomentum();
+ double chCosT = Math.abs(VecOp.cosTheta(pch));
+ double chen = ch.getEnergy();
+ double en_thresh = 0.0;
+ aida.cloud1D("bs/child pdgids").fill(pdgid);
+ aida.cloud1D("bs/child cos theta (lab)").fill(chCosT);
+ aida.cloud1D("bs/child cos theta wtd by en").fill(chCosT, chen/126.0);
+ aida.cloud1D("bs/child momentum").fill(pch.magnitude());
+ if ( (pdgid==11) || (pdgid==13) || (pdgid==15) ) {
+ if ( (chen > en_thresh) && (chCosT < Math.cos(MuC_cone)) ) {
+ leps_muc++;
+ }
+ if ( (chen > en_thresh) && (chCosT < Math.cos(CLIC_cone)) ) {
+ leps_clic++;
+ }
+ } else if (ch.getCharge() != 0.0) {
+ if ( (chen > en_thresh) && chCosT < Math.cos(MuC_cone)) {
+ trks_muc++;
+ }
+ if ( (chen > en_thresh) && chCosT < Math.cos(CLIC_cone)) {
+ trks_clic++;
+ }
+ }
+ if ( (pdgid!=12) && (pdgid!=14) && (pdgid!=16)) {
+ totEn += chen;
+ if (chCosT < Math.cos(MuC_cone)) {
+ totEnMuC += chen;
+ }
+ if (chCosT < Math.cos(CLIC_cone)) {
+ totEnCLIC += chen;
+ }
+
+ }
+ }
+ if ( (leps_muc < 1) && (trks_muc < 2) ) {
+ decaysInMuCCone = false;
+ }
+ if ( (leps_clic < 1) && (trks_clic < 2) ) {
+ decaysInCLICCone = false;
+ }
+ aida.cloud1D("bs/leps_muc").fill(leps_muc);
+ aida.cloud1D("bs/leps_clic").fill(leps_clic);
+ aida.cloud1D("bs/trks_muc").fill(trks_muc);
+ aida.cloud1D("bs/trks_clic").fill(trks_clic);
+ }
+ if (decaysInMuCCone) {
+ aida.cloud1D("bs/all bs detected (0->MuC, 1->CLIC)").fill(0);
+ numPassMuCTest += 1;
+ }
+ if (decaysInCLICCone) {
+ aida.cloud1D("bs/all bs detected (0->MuC, 1->CLIC)").fill(1);
+ numPassCLICTest += 1;
+ }
+ aida.cloud1D("bs/fraction visible energy - CLIC").fill(totEnCLIC/totEn);
+ aida.cloud1D("bs/fraction visible energy - MuC").fill(totEnMuC/totEn);
+ }
+ * */
+ /*
+ * Final state particle-based analysis (except for IS neutrinos!)
+ */
+ /*
+ double totalEnergy = 0.0;
+ double totalEnergy_nucut = 0.0;
+ for (MCParticle p:allFSPs) {
+ totalEnergy += p.getEnergy();
+ int pdgid = Math.abs(p.getPDGID());
+ // neutrino cut
+ if ( (pdgid!=12) && (pdgid!=14) && (pdgid!=16)) {
+ totalEnergy_nucut += p.getEnergy();
+ }
+ }
+ aida.cloud1D("fs/total energy").fill(totalEnergy);
+ aida.cloud1D("fs/total energy - nu cut").fill(totalEnergy_nucut);
+ aida.cloud1D("fs/total energy - fraction nu cut").fill(totalEnergy_nucut/totalEnergy);
+ */
+ /*
+ * Cone angle cuts
+ */
+// coneAnglePlots(allBmesons, "B-mesons");
+// aida.cloud1D("num B-mesons").fill(allBmesons.size());
+// coneAnglePlots(allFSPs,"cone cuts - all");
+ //coneAnglePlots(allFSPs_nucut,"cone cuts - all - nucut");
+// if (
+// (Math.abs(h1decay.get(0).getPDGID())==5) &&
+// (Math.abs(h2decay.get(0).getPDGID())==5) ) {
+// coneAnglePlots(allFSPs,"cone cuts - b-bbar");
+// coneAnglePlots(allFSPs_nucut,"cone cuts - b-bbar - nucut");
+// }
+// else if (
+// (Math.abs(h1decay.get(0).getPDGID())==24) &&
+// (Math.abs(h2decay.get(0).getPDGID())==24) ) {
+// coneAnglePlots(allFSPs,"cone cuts - WW*");
+// coneAnglePlots(allFSPs_nucut,"cone cuts - WW* - nucut");
+// }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N EvtShapeCutsDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EvtShapeCutsDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,173 @@
+package org.lcsim.mcd.analysis;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.mc.fast.MCFast;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.jet.EventShape;
+import hep.physics.vec.BasicHepLorentzVector;
+import hep.physics.vec.HepLorentzVector;
+import org.lcsim.contrib.Cassell.recon.Cheat.MakePerfectReconParticles;
+
+/**
+ *
+ * @author aconway
+ */
+public class EvtShapeCutsDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ MCFast MCFast;
+ double _ctcut = 0.922;
+
+ public EvtShapeCutsDriver(){
+ MCFast = new MCFast();
+ }
+
+ public void process(EventHeader event) {
+ this.add(MCFast);
+ this.processChildren(event);
+ this.remove(MCFast);
+
+ List<MCParticle> FSChildren = event.get(MCParticle.class,"GenFinalStateParticles");
+ List<MCParticle> ctCutList = cosThetaCut(FSChildren, _ctcut);
+ List<MCParticle> nuCutList = neutrinoCut(FSChildren);
+ List<MCParticle> ctnuCutList = neutrinoCut(ctCutList);
+ event.put("noCut", FSChildren, MCParticle.class,0);
+ event.put("ctCut", ctCutList, MCParticle.class, 0);
+ event.put("nuCut", nuCutList, MCParticle.class, 0);
+ event.put("ctnuCut", ctnuCutList, MCParticle.class, 0);
+
+ MakePerfectReconParticles mprp1 = new MakePerfectReconParticles("noCut","rcp noCut");
+ MakePerfectReconParticles mprp2 = new MakePerfectReconParticles("ctCut","rcp ctCut");
+ MakePerfectReconParticles mprp3 = new MakePerfectReconParticles("nuCut","rcp nuCut");
+ MakePerfectReconParticles mprp4 = new MakePerfectReconParticles("ctnuCut","rcp ctnuCut");
+ this.add(mprp1);
+ this.add(mprp2);
+ this.add(mprp3);
+ this.add(mprp4);
+ this.processChildren(event);
+ this.remove(mprp1);
+ this.remove(mprp2);
+ this.remove(mprp3);
+ this.remove(mprp4);
+
+ MCParticle boson = getBoson(event);
+ String MCTag = getMCTag(boson);
+ }
+
+ private List<MCParticle> cosThetaCut(List<MCParticle> mcpl, double ctcut) {
+ List<MCParticle> rv = new ArrayList<MCParticle>();
+ for (MCParticle p:mcpl) {
+ double cosT = Math.abs(VecOp.cosTheta(p.getMomentum()));
+ if (cosT < ctcut) {
+ rv.add(p);
+ }
+ }
+ return rv;
+ }
+
+ private List<MCParticle> neutrinoCut(List<MCParticle> mcpl) {
+ List<MCParticle> rv = new ArrayList<MCParticle>();
+ for(MCParticle p:mcpl) {
+ if( !( (p.getPDGID() == 12)
+ || (p.getPDGID() == 14)
+ || (p.getPDGID() == 16)
+ || (p.getPDGID() == 18) ) )
+ rv.add(p);
+ }
+ return rv;
+ }
+ private MCParticle getBoson(EventHeader event) {
+ MCParticle boson = event.getMCParticles().get(0);
+ for (MCParticle p:event.getMCParticles()) {
+ if ( ((p.getPDGID() == 25) || (p.getPDGID() == 23))
+ && (p.getGeneratorStatus() == 3)
+ && (Math.abs(p.getParents().get(0).getPDGID()) == 13)) {
+ //&& (p.getMass() > 40.0) ){
+ boson = p;
+ }
+ }
+ return boson;
+ }
+
+ private String getMCTag(MCParticle boson) {
+ List<MCParticle> children = boson.getDaughters();
+
+ String chan = "Unknown";
+ if (!children.isEmpty()) {
+ int pdgid1 = 0;
+ int pdgid2 = 0;
+ boolean ffound = false;
+ boolean sfound = false;
+ for (MCParticle p:children) {
+ if (p.getGeneratorStatus()==3) { // Documentation state
+ int pdgid = p.getPDGID();
+
+ if ( !ffound && !sfound ) {
+ pdgid1 = pdgid;
+ ffound = true;
+ }
+ if ( ffound && !sfound ) {
+ pdgid2 = pdgid;
+ sfound = true;
+ }
+ }
+ }
+ if (ffound && sfound) {
+ if (Math.abs(pdgid1)==Math.abs(pdgid2)) {
+ aida.cloud1D("Decay Channels").fill(Math.abs(pdgid1));
+ switch (Math.abs(pdgid1)) {
+ case 1 :
+ chan = "UDS";
+ break;
+ case 2 :
+ chan = "UDS";
+ break;
+ case 3 :
+ chan = "UDS";
+ break;
+ case 4 :
+ chan = "CCBAR";
+ break;
+ case 5 :
+ chan = "BBBAR";
+ break;
+ case 6 :
+ chan = "TTBAR";
+ break;
+ case 15 :
+ chan = "TAUTAU";
+ break;
+ case 21 :
+ chan = "GLUGLU";
+ break;
+ case 22 :
+ chan = "GAMGAM";
+ break;
+ case 23 :
+ chan = "ZZ";
+ break;
+ case 24 :
+ chan = "WW";
+ break;
+ default :
+ chan = "OTHER";
+ break;
+ }
+ }
+ else {
+ aida.cloud1D("Unknown decay channels by child PDGID").fill(pdgid1);
+ aida.cloud1D("Unknown decay channels by child PDGID").fill(pdgid2);
+ }
+ }
+ }
+ return chan;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N MCPContribTestDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCPContribTestDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,232 @@
+package org.lcsim.mcd.analysis;
+
+import hep.aida.IAxis;
+import hep.aida.IHistogram1D;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.base.BaseCalorimeterHit;
+import org.lcsim.event.base.BaseSimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+/**
+ * Driver to test whether we lose anything by combining hit contributions
+ * from the same MCParticles.
+ * If not, then we won't need to run SLIC with the /lcio/PDGFlag option,
+ * saving simulation time and file size.
+ * Want to look at the time delay of the edeps compared to the time of the
+ * first hit, and compared to the first hit from that MCP.
+ * Want to look at the effect of timing cuts on energy response.
+ * Also want to determine if we lose much by combining all the contributions.
+ * @author aconway
+ */
+public class MCPContribTestDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ String[] pfixes = {"Edep/raw/","Edep/comb/","Opti/raw/","Opti/comb/"};
+ String[] collnames = {"Edep","Opti"};
+ double[] timebins = {-0.1,-0.05,0.0,0.1, 0.25, 0.5, 1.0, 2.5, 5.0, 10.0, 25.0, 50.0, 100.0};
+
+ double nevents;
+ Hep3Vector shower_center;
+ double tot_edep;
+ double tot_opti;
+
+ @Override
+ protected void startOfData(){
+ nevents = 0.0;
+ for (String coll : collnames) {
+ aida.histogram1D("contrib plots/"+coll+"/delta t 0-10",
+ 100,0,10);
+ aida.histogram1D("contrib plots/"+coll+"/delta t 0-100",
+ 100,0,100);
+ }
+ for (String pfix:pfixes) {
+ aida.histogram1D(pfix+"time of flight: 0 to 1",
+ 50,-1.0,1.0);
+ aida.histogram1D(pfix+"time of flight: 0 to 5",
+ 50,-1.0,5.0);
+ aida.histogram1D(pfix+"time of flight: 0 to 25",
+ 50,-1.0,25.0);
+ aida.histogram1D(pfix+"time of flight: 0 to 100",
+ 50,-1.0,100.0);
+ aida.histogram1D(pfix+"time of flight: 0 to 100",
+ 50,-1.0,10000.0);
+ }
+ }
+ @Override
+ protected void process(EventHeader event) {
+ tot_edep = 0.0;
+ tot_opti = 0.0;
+ nevents += 1.0;
+
+ List<SimCalorimeterHit> EdepHitList = new ArrayList<SimCalorimeterHit>();
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "ECalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "HCalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "ECalEndcapEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "HCalEndcapEdepHits"));
+ List<SimCalorimeterHit> OptiHitList = new ArrayList<SimCalorimeterHit>();
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "ECalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "HCalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "ECalEndcapOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "HCalEndcapOptiHits"));
+
+ double maxe = 0.0;
+ SimCalorimeterHit centerhit = EdepHitList.get(0);
+ for (SimCalorimeterHit hit:EdepHitList) {
+ double en = hit.getRawEnergy();
+ tot_edep += en;
+ if (en>maxe) {
+ maxe = en;
+ centerhit = hit;
+ }
+ }
+ for (SimCalorimeterHit hit:OptiHitList) {
+ tot_opti += hit.getRawEnergy();
+ }
+ shower_center = centerhit.getPositionVec();
+ aida.cloud1D("shower center radius").fill(shower_center.magnitude());
+
+ List<SimCalorimeterHit> comb_EdepHitList = combineMCPs(EdepHitList, "Edep");
+ List<SimCalorimeterHit> comb_OptiHitList = combineMCPs(OptiHitList, "Opti");
+
+ makeHists(EdepHitList,tot_edep,"Edep/raw/");
+ makeHists(comb_EdepHitList,tot_edep,"Edep/comb/");
+ makeHists(OptiHitList,tot_opti,"Opti/raw/");
+ makeHists(comb_OptiHitList,tot_opti,"Opti/comb/");
+ }
+
+ private void makeHists(List<SimCalorimeterHit> hits, double tot_e, String pfix) {
+ for (SimCalorimeterHit hit : hits) {
+ Hep3Vector pos = hit.getPositionVec();
+ double r = Math.sqrt(pos.x()*pos.x() + pos.y()*pos.y());
+ double dr = VecOp.sub(pos, shower_center).magnitude();
+ double z = pos.z();
+ double t0 = hit.getTime();
+ double tof0 = t0 - pos.magnitude()/299.79;
+ double etot = hit.getRawEnergy();
+ aida.cloud1D(pfix+"tof0 wt by etot").fill(tof0,etot);
+ int ncontribs = hit.getMCParticleCount();
+ aida.cloud1D(pfix+"ncontribs").fill(ncontribs);
+ for (int i=0; i<ncontribs; i++) {
+ MCParticle mcp = hit.getMCParticle(i);
+ double ei = hit.getContributedEnergy(i);
+ double ti = hit.getContributedTime(i);
+ double dt = ti - t0;
+ double tofi = tof0 + dt;
+ aida.cloud1D(pfix+"delta t").fill(dt,ei);
+ aida.cloud2D(pfix+"delta t (x) vs dist from impact (y)").
+ fill(dt, dr, ei);
+ aida.histogram1D(pfix+"time of flight: 0 to 1").fill(tofi,ei/tot_e);
+ aida.histogram1D(pfix+"time of flight: 0 to 5").fill(tofi,ei/tot_e);
+ aida.histogram1D(pfix+"time of flight: 0 to 25").fill(tofi,ei/tot_e);
+ aida.histogram1D(pfix+"time of flight: 0 to 100").fill(tofi,ei/tot_e);
+ }
+ }
+ }
+ @Override
+ protected void endOfData() {
+ for (String pfix:pfixes) {
+ aida.histogram1D(pfix+"time of flight: 0 to 1").scale(1.0/nevents);
+ aida.histogram1D(pfix+"time of flight: 0 to 5").scale(1.0/nevents);
+ aida.histogram1D(pfix+"time of flight: 0 to 25").scale(1.0/nevents);
+ aida.histogram1D(pfix+"time of flight: 0 to 100").scale(1.0/nevents);
+ integrate_histogram(
+ aida.histogram1D(pfix+"time of flight: 0 to 1"),
+ pfix+"integral time of flight: 0 to 1");
+ integrate_histogram(
+ aida.histogram1D(pfix+"time of flight: 0 to 5"),
+ pfix+"integral time of flight: 0 to 5");
+ integrate_histogram(
+ aida.histogram1D(pfix+"time of flight: 0 to 25"),
+ pfix+"integral time of flight: 0 to 25");
+ integrate_histogram(
+ aida.histogram1D(pfix+"time of flight: 0 to 100"),
+ pfix+"integral time of flight: 0 to 100");
+ }
+ }
+
+ /*
+ Take a list of hits from the calorimeter and create new hits with the
+ contributions from each MCP combined.
+ */
+ private List<SimCalorimeterHit> combineMCPs(List<SimCalorimeterHit> hits, String collname) {
+ List<SimCalorimeterHit> rvlist = new ArrayList<SimCalorimeterHit>();
+
+ for (SimCalorimeterHit hit:hits) {
+ List<MCParticle> newmcps = new ArrayList<MCParticle>();
+ List<Float> newens = new ArrayList<Float>();
+ List<Float> newts = new ArrayList<Float>();
+ List<Integer> newpdgs = new ArrayList<Integer>();
+ int ncontribs = hit.getMCParticleCount();
+ // loop over contribs, check if we've seen a contrib from this MCP
+ // in this hit's contribs, if so, increment energy entry, set
+ // new lower time, etc
+ for (int i=0; i<ncontribs; i++) {
+ MCParticle mcp = hit.getMCParticle(i);
+ if (newmcps.contains(mcp)) {
+ int index = newmcps.indexOf(mcp);
+ newens.set(index, (float) hit.getContributedEnergy(i) + newens.get(index));
+ newts.set(index, Math.min((float) hit.getContributedTime(i), newts.get(index)));
+ } else {
+ newmcps.add(mcp);
+ newens.add((float) hit.getContributedEnergy(i));
+ newts.add((float) hit.getContributedTime(i));
+ newpdgs.add(hit.getPDG(i));
+ }
+ }
+ // now, make plots comparing edep timing for individual mcps
+ for (int i=0; i<ncontribs; i++) {
+ MCParticle mcp = hit.getMCParticle(i);
+ double ei = hit.getContributedEnergy(i);
+ double ti = hit.getContributedTime(i);
+ assert newmcps.contains(mcp);
+ int newindex = newmcps.indexOf(mcp);
+ double etot = newens.get(newindex);
+ double t0 = newts.get(newindex);
+ aida.cloud1D("contrib plots/"+collname+"/delta t wt en").
+ fill(ti-t0,ei);
+ aida.histogram1D("contrib plots/"+collname+"/delta t 0-10")
+ .fill(ti-t0,ei);
+ aida.histogram1D("contrib plots/"+collname+"/delta t 0-100")
+ .fill(ti-t0,ei);
+ }
+ float[] ensarray = new float[newens.size()];
+ for (int i=0; i<newens.size(); i++) {
+ ensarray[i] = newens.get(i);
+ }
+ float[] tsarray = new float[newts.size()];
+ for (int i=0; i<newts.size(); i++) {
+ tsarray[i] = newts.get(i);
+ }
+ int[] pdgarray = new int[newpdgs.size()];
+ for (int i=0; i<newpdgs.size(); i++) {
+ pdgarray[i] = newpdgs.get(i);
+ }
+ BaseCalorimeterHit newbasehit = new BaseSimCalorimeterHit(
+ hit.getCellID(),hit.getRawEnergy(), hit.getTime(),
+ newmcps.toArray(), ensarray, tsarray, pdgarray);
+ newbasehit.setPosition(hit.getPosition());
+ SimCalorimeterHit newhit = (SimCalorimeterHit) newbasehit;
+ rvlist.add(newhit);
+ }
+ return rvlist;
+ }
+ private IHistogram1D integrate_histogram(IHistogram1D hist, String newname) {
+ IAxis ax = hist.axis();
+ IHistogram1D rv = aida.histogram1D(newname, ax.bins(), ax.lowerEdge(), ax.upperEdge());
+ for (int i=0; i<ax.bins(); i++) {
+ double sum = 0.0;
+ for (int j=0; j<=i; j++) {
+ sum += hist.binHeight(j);
+ }
+ rv.fill(ax.binLowerEdge(i),sum);
+ }
+ return rv;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N mcpHitContribAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ mcpHitContribAnalysisDriver.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,78 @@
+
+package org.lcsim.mcd.analysis;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+/**
+ *
+ * @author aconway
+ */
+public class mcpHitContribAnalysisDriver extends Driver {
+ private AIDA aida = AIDA.defaultInstance();
+
+ @Override
+ protected void process(EventHeader event) {
+ List<List<MCParticle>> mcpList = event.get(MCParticle.class);
+ List<SimCalorimeterHit> EdepHitList = new ArrayList<SimCalorimeterHit>();
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "ECalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "HCalBarrelEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "ECalEndcapEdepHits"));
+ EdepHitList.addAll(event.get(SimCalorimeterHit.class, "HCalEndcapEdepHits"));
+ List<SimCalorimeterHit> OptiHitList = new ArrayList<SimCalorimeterHit>();
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "ECalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "HCalBarrelOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "ECalEndcapOptiHits"));
+ OptiHitList.addAll(event.get(SimCalorimeterHit.class, "HCalEndcapOptiHits"));
+ double sumEdep = 0.0;
+ for (SimCalorimeterHit hit : EdepHitList) {
+ double t0 = hit.getTime();
+ double tofl0 = getToFl(hit);
+ double e0 = hit.getRawEnergy();
+ aida.cloud1D("Edep/hit/time").fill(tofl0,e0);
+ int nmc = hit.getMCParticleCount();
+ aida.cloud1D("Edep/num MCP contribs").fill(nmc);
+ for (int i=0; i<nmc; i++) {
+ double e = hit.getContributedEnergy(i);
+ sumEdep += e;
+ double t = hit.getContributedTime(i);
+ double tofl = getToFl(hit,t);
+ aida.cloud1D("Edep/contrib/time").fill(tofl,e);
+ aida.cloud1D("Edep/contrib/delta time").fill(t-t0,e);
+ int pdgid = hit.getMCParticle(i).getPDGID();
+ aida.cloud1D("Edep/contrib/pdgid "+String.valueOf(pdgid)+"/time").fill(tofl,e);
+ aida.cloud1D("Edep/contrib/pdgid "+String.valueOf(pdgid)+"/delta time").fill(tofl,e);
+ }
+ }
+ for (SimCalorimeterHit hit : OptiHitList) {
+ double t0 = hit.getTime();
+ double tofl0 = getToFl(hit);
+ double e0 = hit.getRawEnergy();
+ aida.cloud1D("Opti/hit/time").fill(tofl0,e0);
+ int nmc = hit.getMCParticleCount();
+ aida.cloud1D("Opti/num MCP contribs").fill(nmc);
+ for (int i=0; i<nmc; i++) {
+ double e = hit.getContributedEnergy(i);
+ double t = hit.getContributedTime(i);
+ double tofl = getToFl(hit,t);
+ aida.cloud1D("Opti/contrib/time").fill(tofl,e);
+ aida.cloud1D("Opti/contrib/delta time").fill(t-t0,e);
+ int pdgid = hit.getMCParticle(i).getPDGID();
+ aida.cloud1D("Opti/contrib/pdgid "+String.valueOf(pdgid)+"/time").fill(tofl,e);
+ aida.cloud1D("Opti/contrib/pdgid "+String.valueOf(pdgid)+"/delta time").fill(tofl,e);
+ }
+ }
+ aida.cloud1D("Sum Edep from MCPContribs").fill(sumEdep);
+ }
+
+ private double getToFl(SimCalorimeterHit hit) {
+ return hit.getTime() - hit.getPositionVec().magnitude()/299.79;
+ }
+ private double getToFl(SimCalorimeterHit hit, double time) {
+ return time - hit.getPositionVec().magnitude()/299.79;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N ElectronCorrection.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ElectronCorrection.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,348 @@
+//package org.lcsim.mcd.analysis;
+
+/**
+ *
+ * @author wenzel
+ * This Driver is used to calibrate the correction of a dual readout calorimeter to electrons.
+ * (Determining the energy scale of the calorimeter response.)
+ */
+import hep.aida.ICloud1D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import java.io.BufferedWriter;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.aida.AIDA;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+
+public class ElectronCorrection extends Driver {
+
+ private AIDA aida;
+ String AIDAFile;
+ String file_name;
+ FileWriter fstream;
+ BufferedWriter out;
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsFactory;
+ IFunction gauss;
+ boolean first;
+ boolean firstEvent;
+ double E_in;
+ Integer Ein;
+ Integer Ein_prev;
+ double E_kin;
+ String Particlename;
+ ICloud1D Edep;
+ ICloud1D Eceren;
+ double nsigmas;
+ int nbins;
+ IDataPointSet dps_emean = null;
+ IDataPointSet dps_cerenemean = null;
+ IDataPointSet[] dps_cerene = null;
+ int point_cerene[];
+ int point;
+ int point_cerenemean;
+ double[] result;
+ double errors[];
+ String[] Fitters = {"Chi2", "leastsquares", "bml", "cleverchi2"};
+ IFitter jminuit;
+ IFitResult jminuitResult;
+
+ public ElectronCorrection() {
+ file_name = "ElectronCorrection.txt";
+ AIDAFile = "ElectronCorrection.aida";
+ aida = AIDA.defaultInstance();
+ dps_cerene = new IDataPointSet[Fitters.length];
+ point_cerene = new int[Fitters.length];
+ fitFactory = aida.analysisFactory().createFitFactory();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_cerene_" + Fitters[i];
+ point_cerene[i] = 0;
+ dps_cerene[i] = dpsFactory.create(dpsname, "electron response mean", 2);
+ }
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+ // Create a two dimensional IDataPointSet.
+ dps_emean = dpsFactory.create("dps_emean", "electron response mean", 2);
+ dps_cerenemean = dpsFactory.create("dps_cerenemean", "electron response mean", 2);
+ point = 0;
+ point_cerenemean = 0;
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ }
+
+ @Override
+ protected void process(EventHeader event) {
+ E_in = 0.0;
+ E_kin = 0.0;
+ String DirName = null;
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ if (particle.getGeneratorStatus() == 0) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+ if (particle.getProductionTime() == 0.0) {
+ Particlename = particle.getType().getName();
+ }
+ }
+ break;
+ }
+ Ein = (int) Math.floor(E_kin + 0.5d);
+ String E_str = Ein.toString();
+ DirName = Particlename.concat(E_str);
+ if (Ein != Ein_prev) {
+ if (first) {
+ first = false;
+ } else {
+ System.out.println("ElectronCorrection:First Event");
+ convertandfit();
+ }
+ Ein_prev = Ein;
+ System.out.println("First Event:");
+ System.out.println("E_in: " + E_in);
+ System.out.println("E_kin: " + E_kin);
+ System.out.println("Name: " + Particlename);
+ aida.tree().cd("/");
+ aida.tree().mkdir(DirName);
+ aida.tree().cd(DirName);
+ Edep = aida.histogramFactory().createCloud1D("Edep", "Energy Cloud", 100000, "autoConvert = false");
+ Eceren = aida.histogramFactory().createCloud1D("Eceren", "Cerenkov Cloud", 100000, "autoConvert = false");
+ System.out.println("DirName: " + DirName);
+ }
+ List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+ double sumEEdep = 0.0;
+ double sumECeren = 0.0;
+ for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+ String CollectionName = event.getMetaData(simCalorimeterHits).getName();
+
+ if (CollectionName.contains("Edep")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ } // end if Edep
+ if (CollectionName.contains("Opti")) {
+// if (CollectionName.startsWith("Ceren_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ } // end if Ceren
+ } // end loop over calorimeter hit collections
+
+ Edep.fill(sumEEdep);
+ Eceren.fill(sumECeren);
+ }
+
+ @Override
+ protected void endOfData() {
+ System.out.println("ElectronCorrection:End of Data:");
+ convertandfit();
+ fstream = null;
+ try {
+ fstream = new FileWriter(file_name);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ out = new BufferedWriter(fstream);
+ aida.tree().cd("/");
+ IFunction line = functionFactory.createFunctionByName("line", "p1");
+ IFitter linefit = fitFactory.createFitter("Chi2", "jminuit");
+ IFitResult resultline = linefit.fit(dps_emean, line);
+ functionFactory.cloneFunction("e mean fitted line ", resultline.fittedFunction());
+ double[] fresult = resultline.fittedParameters();
+ try {
+ out.write("Ionization scale results:\n");
+ out.write("Chi2 = " + resultline.quality() + "\n");
+ out.write(fresult[0] + " , " + fresult[1] + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ //
+ // now deal with the cerenkov stuff:
+ //
+ try {
+ out.write("Cerenkov scale results:\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ for (int i = 0; i < Fitters.length; i++) {
+ System.out.println("Fitter: " + Fitters[i]);
+ resultline = linefit.fit(dps_cerene[i], line);
+ String fname = "cerenkov " + Fitters[i] + " fitted line";
+ System.out.println(fname);
+ functionFactory.cloneFunction(fname, resultline.fittedFunction());
+ System.out.println(Fitters[i] + ": " + resultline.quality());
+ fresult = resultline.fittedParameters();
+ try {
+ out.write("Fitter: " + Fitters[i]+"\n");
+ out.write("Chi2 = " + resultline.quality() + "\n");
+ out.write(fresult[0] + " , " + fresult[1] + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ try {
+ out.close();
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ @Override
+ protected void resume() {
+ System.out.println("ElectronCorrection:resume");
+ firstEvent =true;
+ aida.cloud1D("c_Edep_energy").reset();
+ }
+
+ protected void convertandfit() {
+ System.out.println("ElectronCorrection:convertandfit");
+ IHistogram1D Edep_conv;
+ if (Edep.isConverted()) {
+ Edep_conv = Edep.histogram();
+ } else {
+ System.out.println("Converting EDep");
+ double meanc = Edep.mean();
+ double rmsc = Edep.rms();
+ nsigmas = 3.;
+ nbins = 100;
+
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep.setConversionParameters(nbins, minx, maxx);
+ Edep.convertToHistogram();
+ Edep_conv = Edep.histogram();
+ }
+
+ gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_conv.mean());
+ gauss.setParameter("sigma", Edep_conv.rms());
+
+ dps_emean.addPoint();
+ IDataPoint dp = dps_emean.point(point);
+ dp.coordinate(1).setValue(Ein_prev);
+ dp.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+ dp.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+ dp.coordinate(0).setValue(Edep_conv.mean());
+ dp.coordinate(0).setErrorPlus(Edep_conv.rms());
+ dp.coordinate(0).setErrorMinus(Edep_conv.rms());
+ point++;
+// chi 2 fit:
+ System.out.println("chi2 fit:");
+ /*jminuit = fitFactory.createFitter("Chi2", "jminuit");
+ jminuitResult = jminuit.fit(Edep_conv, gauss);
+ System.out.println(Edep_conv.maxBinHeight() + " , " + Edep_conv.mean() + " , " + Edep_conv.rms());
+ System.out.println("jminuit Chi2=" + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ functionFactory.cloneFunction("fitted gauss (jminuitchi2)", jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+
+ // least squares fit:
+ System.out.println("least squares fit:");
+ gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_conv.mean());
+ gauss.setParameter("sigma", Edep_conv.rms());
+ IFitter jminuitls = fitFactory.createFitter("leastsquares", "jminuit");
+ IFitResult jminuitResultls = jminuitls.fit(Edep_conv, gauss);
+ System.out.println("jminuit ls=" + jminuitResultls.quality());
+ double[] resultls = jminuitResultls.fittedParameters();
+ functionFactory.cloneFunction("fitted gauss (jminuitls)", jminuitResultls.fittedFunction());
+ System.out.println(resultls[0] + "," + resultls[1] + "," + resultls[2]);
+*/
+ //
+ // now deal with cerenkov response:
+ //
+ IHistogram1D Eceren_conv;
+
+ if (Eceren.isConverted()) {
+ Eceren_conv = Eceren.histogram();
+ } else {
+ System.out.println("Converting Eceren");
+ double meanc = Eceren.mean();
+ double rmsc = Eceren.rms();
+ nsigmas = 5.;
+ nbins = 100;
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Eceren.setConversionParameters(nbins, minx, maxx);
+ Eceren.convertToHistogram();
+ Eceren_conv = Eceren.histogram();
+ }
+
+ gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+ gauss.setParameter("mean", Eceren_conv.mean());
+ gauss.setParameter("sigma", Eceren_conv.rms());
+ dps_cerenemean.addPoint();
+ IDataPoint dp_cerenemean = dps_cerenemean.point(point_cerenemean);
+ dp_cerenemean.coordinate(1).setValue(Ein_prev);
+ dp_cerenemean.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+ dp_cerenemean.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+ dp_cerenemean.coordinate(0).setValue(Eceren_conv.mean());
+ dp_cerenemean.coordinate(0).setErrorPlus(Eceren_conv.rms());
+ dp_cerenemean.coordinate(0).setErrorMinus(Eceren_conv.rms());
+ point_cerenemean++;
+
+ for (int i = 0; i < Fitters.length; i++) {
+ System.out.println("Fitter: " + Fitters[i]);
+ gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+ gauss.setParameter("mean", Eceren_conv.mean());
+ gauss.setParameter("sigma", Eceren_conv.rms());
+ jminuit = fitFactory.createFitter(Fitters[i], "jminuit");
+ jminuitResult = jminuit.fit(Eceren_conv, gauss);
+ System.out.println("jminuit " + Fitters[i] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ String functionname = "ceren fitted gauss " + Fitters[i];
+ System.out.println(functionname);
+ functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+ IDataPoint dp_cerene = dps_cerene[i].addPoint();
+ dp_cerene.coordinate(0).setValue(result[1]);
+ dp_cerene.coordinate(0).setErrorPlus(Math.abs(result[2]));
+ dp_cerene.coordinate(0).setErrorMinus(Math.abs(result[2]));
+ dp_cerene.coordinate(1).setValue(Ein_prev);
+ dp_cerene.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+ dp_cerene.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+ point_cerene[i]++;
+ }
+
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("ElectronCalibrationDriver:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+
+ public void setMyFilename(String file_name) {
+ System.out.println("ElectronCalibrationDriver:setMyFilename");
+ this.file_name = file_name;
+ System.out.println(file_name);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N evtToTxt.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ evtToTxt.java 20 Mar 2014 20:24:46 -0000 1.1
@@ -0,0 +1,70 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+//package org.lcsim.mcd.analysis;
+
+import org.lcsim.util.Driver;
+import java.io.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import java.util.ArrayList;
+import java.util.List;
+
+/**
+ *
+ * @author aconway
+ */
+public class evtToTxt extends Driver {
+ String ofFolder;
+ FileWriter fstream;
+ BufferedWriter out;
+ int evtNum;
+
+ public evtToTxt() {
+ evtNum = 0;
+ ofFolder = "/home/aconway/fermi/jetFinder/events/";
+ }
+
+ public void process(EventHeader event) {
+ evtNum += 1;
+ String ofName = ofFolder+evtNum+".txt";
+ try {
+ fstream = new FileWriter(ofName);
+ out = new BufferedWriter(fstream);
+ } catch (Exception e) {
+ System.err.println("Error opening output file: \n"+ofName+"\n"+e.getMessage());
+ }
+ int numParticles = 0;
+
+ List<MCParticle> FSList = new ArrayList<MCParticle>();
+ List<MCParticle> MCPList = event.getMCParticles();
+
+ for (MCParticle p:MCPList) {
+ if (p.getGeneratorStatus() == 1) {
+ FSList.add(p);
+ }
+ }
+
+ for (MCParticle p:FSList) {
+ double px = p.getMomentum().x();
+ double py = p.getMomentum().y();
+ double pz = p.getMomentum().z();
+ double en = p.getEnergy();
+ try {
+ out.write(px+","+py+","+pz+","+en+"\n");
+ } catch (Exception e) {
+ System.err.println("Error writing to output file: \n"+ofName+"\n"+e.getMessage());
+ }
+ }
+
+ try {
+ out.close();
+ } catch (Exception e) {
+ System.err.println("Error closing output file: \n"+ofName+"\n"+e.getMessage());
+ }
+ }
+
+ public void endOfData() {
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -u -r1.2 -r1.3
--- HiggsAnalysisDriver.java 19 Mar 2013 20:27:06 -0000 1.2
+++ HiggsAnalysisDriver.java 20 Mar 2014 20:24:46 -0000 1.3
@@ -1,4 +1,4 @@
-package org.lcsim.mcd.analysis;
+//package org.lcsim.mcd.analysis;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
@@ -27,7 +27,7 @@
MCFast MCFast;
JetDriver jDrive;
JetFinder jFind;
- double _ctcut = 0.9;
+ double _ctcut = 0.966;
double _jetcut = 0.2;
public HiggsAnalysisDriver() {
@@ -48,64 +48,170 @@
aida.cloud1D("boson invariant masses")
.fill(bosonInvMass);
+ boolean makespcut = partonCut(boson, 15.0);
+ if (makespcut){
+ aida.cloud1D(MCTag+" makes parton cos(t) cut - cos(t)").fill(VecOp.cosTheta(boson.getDaughters().get(0).getMomentum()));
- List<MCParticle> FSChildren = event.get(MCParticle.class,"GenFinalStateParticles");
-
- List<MCParticle> ctCutList = cosThetaCut(FSChildren, _ctcut);
- List<MCParticle> nuCutList = neutrinoCut(FSChildren);
- List<MCParticle> ctnuCutList = neutrinoCut(ctCutList);
+ List<MCParticle> FSChildren = event.get(MCParticle.class,"GenFinalStateParticles");
+
+ List<MCParticle> ctCutList = cosThetaCut(FSChildren, _ctcut);
+ List<MCParticle> nuCutList = neutrinoCut(FSChildren);
+ List<MCParticle> ctnuCutList = neutrinoCut(ctCutList);
+
+ event.put("noCut", FSChildren, MCParticle.class,0);
+ event.put("ctCut", ctCutList, MCParticle.class, 0);
+ event.put("nuCut", nuCutList, MCParticle.class, 0);
+ event.put("ctnuCut", ctnuCutList, MCParticle.class, 0);
+
+ MakePerfectReconParticles mprp1 = new MakePerfectReconParticles("noCut","rcp noCut");
+ MakePerfectReconParticles mprp2 = new MakePerfectReconParticles("ctCut","rcp ctCut");
+ MakePerfectReconParticles mprp3 = new MakePerfectReconParticles("nuCut","rcp nuCut");
+ MakePerfectReconParticles mprp4 = new MakePerfectReconParticles("ctnuCut","rcp ctnuCut");
+
+ this.add(mprp1);
+ this.add(mprp2);
+ this.add(mprp3);
+ this.add(mprp4);
+ this.processChildren(event);
+ this.remove(mprp1);
+ this.remove(mprp2);
+ this.remove(mprp3);
+ this.remove(mprp4);
+
+ String evt;
+ boolean isWW = false;
+ if (MCTag.equals("BBBAR")) {
+ evt = "b-bbar/";
+ } else if (MCTag.equals("WW")) {
+ isWW = true;
+ String type = getWWDecay(boson.getDaughters());
+ if (type.equals("jjlnu")) {
+ evt = "jjlnu/";
+ } else { evt = "rest/"; }
+ } else { evt = "rest/"; }
+
+ EventShape es = getEventShape(ctnuCutList);
+ double tot_en = 0.0;
+ for (MCParticle p:ctnuCutList) {
+ tot_en += p.getEnergy();
+ }
+ boolean makesCut = evtShapeCut(es, tot_en, 98.0,
+ 0.94,1.0,
+ 0.00,0.20,
+ 0.0,1.0,
+ 0.0,1.0);
+
+ if (boson.getMass() >= 124.9) {
+ evt = evt+"Z*/";
+ }
+ // else if (boson.getMass() >= 110) {
+ // evt = "High Z*/";
+ // }
+ else if (boson.getMass() >= 70) {
+ evt = evt+"Z+gamma/";
+ } else {
+ evt = evt+"Low Z*/";
+ }
+ makePlots(event, FSChildren, "noCut", MCTag, evt);
+ // makePlots(event, ctCutList, "ctCut", MCTag, evt);
+ makePlots(event, nuCutList, "nuCut", MCTag, evt);
+ makePlots(event, nuCutList, "nuCut", MCTag, "all/");
+ makePlots(event, ctnuCutList, "ctnuCut", MCTag, evt);
+ makePlots(event, ctnuCutList, "ctnuCut", MCTag, "all/");
+ if(makesCut){
+ makePlots(event,ctnuCutList, "ctnuCut", MCTag,"all/cuts/");
+ }
- event.put("noCut", FSChildren, MCParticle.class,0);
- event.put("ctCut", ctCutList, MCParticle.class, 0);
- event.put("nuCut", nuCutList, MCParticle.class, 0);
- event.put("ctnuCut", ctnuCutList, MCParticle.class, 0);
-
- MakePerfectReconParticles mprp1 = new MakePerfectReconParticles("noCut","rcp noCut");
- MakePerfectReconParticles mprp2 = new MakePerfectReconParticles("ctCut","rcp ctCut");
- MakePerfectReconParticles mprp3 = new MakePerfectReconParticles("nuCut","rcp nuCut");
- MakePerfectReconParticles mprp4 = new MakePerfectReconParticles("ctnuCut","rcp ctnuCut");
-
- this.add(mprp1);
- this.add(mprp2);
- this.add(mprp3);
- this.add(mprp4);
- this.processChildren(event);
- this.remove(mprp1);
- this.remove(mprp2);
- this.remove(mprp3);
- this.remove(mprp4);
-
- String evt;
- if (boson.getMass() >= 124.9) {
- evt = "Z*/";
- }
-// else if (boson.getMass() >= 110) {
-// evt = "High Z*/";
-// }
- else if (boson.getMass() >= 70) {
- evt = "Z+gamma/";
} else {
- evt = "Low Z*/";
+ aida.cloud1D(MCTag+" misses parton cos(t) cut - cos(t)").fill(VecOp.cosTheta(boson.getDaughters().get(0).getMomentum()));
+ }
+ //makeCutPlots(ctnuCutList, makesWWCut, isWW, "WW cut", MCTag, evt);
+ }
+
+ boolean partonCut(MCParticle boson, double cone_angle) {
+ List<MCParticle> ds = boson.getDaughters();
+ boolean rv = true;
+ for (MCParticle mpc:ds) {
+ if ( (mpc.getGeneratorStatus()==3)
+ && (VecOp.cosTheta(mpc.getMomentum()) > Math.cos(Math.toRadians(cone_angle)))) {
+ rv = false;
+ }
}
- makePlots(event, FSChildren, "noCut", MCTag, evt);
- makePlots(event, ctCutList, "ctCut", MCTag, evt);
- makePlots(event, nuCutList, "nuCut", MCTag, evt);
- makePlots(event, ctnuCutList, "ctnuCut", MCTag, evt);
+ return rv;
+ }
+
+ private void makeCutPlots(List<MCParticle> mcpl, boolean makesCut, boolean mctruth, String cutName, String MCTag, String evt) {
+ String prefix;
+ if (makesCut) {
+ prefix = "makes "+cutName;
+ } else {
+ prefix = "fails "+cutName;
+ }
+ String truth;
+ if (mctruth) {
+ truth = "reals";
+ } else {
+ truth = "fakes";
+ }
+ makeThrustPlots(mcpl,evt+prefix+"/ALL");
+ makeThrustPlots(mcpl,evt+prefix+"/"+MCTag);
+ makeThrustPlots(mcpl,evt+prefix+"/"+truth);
+ makeEnergyPlots(mcpl,evt+prefix+"/ALL");
+ makeEnergyPlots(mcpl,evt+prefix+"/"+MCTag);
+ makeEnergyPlots(mcpl,evt+prefix+"/"+truth);
}
private void makePlots(EventHeader event, List<MCParticle> mcpl, String prefix, String decay, String evt) {
- makeJetPlots(event,prefix,"/ALL",evt);
+ //makeJetPlots(event,prefix,"/ALL",evt);
makeThrustPlots(mcpl,evt+prefix+"/ALL");
- makeGammaPlots(mcpl,evt+prefix+"/ALL");
+ //makeGammaPlots(mcpl,evt+prefix+"/ALL");
makeBosonPlots(mcpl,evt+prefix+"/ALL", event);
+ makeEnergyPlots(mcpl,evt+prefix+"/ALL");
- makeJetPlots(event,prefix,decay,evt);
+ //makeJetPlots(event,prefix,decay,evt);
makeThrustPlots(mcpl,evt+prefix+"/"+decay);
- makeGammaPlots(mcpl,evt+prefix+"/"+decay);
+ //makeGammaPlots(mcpl,evt+prefix+"/"+decay);
+ makeEnergyPlots(mcpl,evt+prefix+"/"+decay);
if(decay.equals("WW")) makeWWPlots(mcpl,evt+prefix+"/"+decay, event);
}
+ private void makeEnergyPlots(List<MCParticle> mcpl, String prefix) {
+ prefix = prefix + "/Energy/";
+ double sum_en = 0.0;
+ double sum_et = 0.0;
+ double sum_mom = 0.0;
+ Hep3Vector net_mom = new BasicHep3Vector(0.0,0.0,0.0);
+ for(MCParticle p:mcpl) {
+ double en = p.getEnergy();
+ Hep3Vector mom = p.getMomentum();
+ double cosT = VecOp.cosTheta(mom);
+ sum_en += en;
+ net_mom = VecOp.add(net_mom, mom);
+ sum_mom += mom.magnitude();
+ sum_et += Math.abs(en*(mom.x()*mom.x()+mom.y()*mom.y())/mom.magnitudeSquared());
+ int pdgid = Math.abs(p.getPDGID());
+ if((( pdgid == 12)
+ || (pdgid == 14)
+ || (pdgid == 16)
+ || (pdgid == 18) ) ){
+ aida.cloud1D(prefix+"neutrinos/energy").fill(en);
+ aida.cloud1D(prefix+"neutrinos/cosTheta").fill(cosT);
+ }
+ }
+ double tot_inv_mass = Math.sqrt(sum_en*sum_en - net_mom.magnitude()*net_mom.magnitude());
+ aida.cloud1D(prefix+"Total Energy").fill(sum_en);
+ aida.cloud1D(prefix+"Total Inv mass").fill(tot_inv_mass);
+ aida.cloud1D(prefix+"Net Momentum - Mag.").fill(net_mom.magnitude());
+ aida.cloud1D(prefix+"Net Momentum - Cos Theta").fill(VecOp.cosTheta(net_mom));
+ aida.cloud1D(prefix+"Net Transverse Momentum").fill(Math.sqrt(net_mom.x()*net_mom.x() + net_mom.y()*net_mom.y()));
+ aida.cloud1D(prefix+"Net Colinear Momentum").fill(Math.abs(net_mom.z()));
+ aida.cloud1D(prefix+"Total Transverse Energy").fill(sum_et);
+ if (sum_en > 98.0) {
+ aida.cloud1D(prefix+"Pass Energy Cut").fill(sum_en);
+ }
+ }
+
private void makeWWPlots(List<MCParticle> mcpl, String prefix, EventHeader event) {
List<MCParticle> Ws = new ArrayList<MCParticle>();
MCParticle boson = getBoson(event);
@@ -151,9 +257,10 @@
aida.cloud1D(prefix+"boson cos theta").fill(VecOp.cosTheta(p));
Hep3Vector p_t = new BasicHep3Vector(p.x(),p.y(),0);
aida.cloud1D(prefix+"boson transverse momentum").fill(p_t.magnitude());
+ aida.cloud1D(prefix+"boson colinear momentum").fill(p.z());
for (MCParticle mcp:mcpl) {
- if( (mcp.getPDGID()==22) && (mcp.getEnergy()>20.0)) {
+ if( (mcp.getPDGID()==22) && (mcp.getEnergy()>1.0)) {
Hep3Vector pg = mcp.getMomentum();
double imass = Math.sqrt(
(mcp.getEnergy()+boson.getEnergy())*(mcp.getEnergy()+boson.getEnergy())
@@ -164,6 +271,9 @@
aida.cloud2D(prefix+"boson+gamma momentum correl.").fill(p.magnitude(),pg.magnitude());
aida.cloud2D(prefix+"boson+gamma cos theta correl.").fill(VecOp.cosTheta(p), VecOp.cosTheta(pg));
aida.cloud2D(prefix+"boson+gamma mass-energy correl.").fill(boson.getMass(),mcp.getEnergy());
+ if (mcp.getOrigin().magnitude() == 0.0 ) {
+ aida.cloud2D("Gamma Energy vs. Boson Mass").fill(boson.getMass(),mcp.getEnergy());
+ }
aida.cloud1D(prefix+"boson inv mass ct cut").fill(invMass);
}
}
@@ -175,7 +285,8 @@
private void makeGammaPlots(List<MCParticle> mcpl, String prefix) {
for (MCParticle p:mcpl) {
- if ( (p.getPDGID()==22) && (p.getEnergy()>20.0)) {
+ Hep3Vector orig = p.getOrigin();
+ if ( (p.getPDGID()==22) && (p.getEnergy()>0.1) && (orig.magnitude()==0.0)) {
aida.cloud2D(prefix+"/gamma/costheta vs en"
).fill(p.getEnergy(),VecOp.cosTheta(p.getMomentum()));
aida.cloud1D(prefix+"/gamma/costheta").fill(VecOp.cosTheta(p.getMomentum()));
@@ -191,16 +302,19 @@
Hep3Vector thrust = es.thrust();
Hep3Vector ThrustAxis = es.thrustAxis();
double ob = es.oblateness();
-
- if(thrust.x()>0) aida.cloud1D(prefix+"/Thrst").fill(thrust.x());
+ prefix = prefix+"/EvtShape";
+ if(thrust.x()>-1) {
+ aida.cloud1D(prefix+"/Thrust").fill(thrust.x());
+ aida.cloud2D(prefix+"/Correl Thrust vs Thrust Axis CosT"
+ ).fill(thrust.x(), VecOp.cosTheta(ThrustAxis));
+ aida.cloud2D(prefix+"/Correl Thrust Major").fill(thrust.x(), thrust.y());
+ aida.cloud2D(prefix+"/Correl Thrust Oblateness").fill(thrust.x(), ob);
+ aida.cloud2D(prefix+"/Thrust minus Maj vs Obl").fill(thrust.x() - thrust.y(), ob);
+ }
aida.cloud1D(prefix+"/Major Axis").fill(thrust.y());
aida.cloud1D(prefix+"/Minor Axis").fill(thrust.z());
aida.cloud1D(prefix+"/Oblateness").fill(ob);
- aida.cloud2D(prefix+"/Correl Thrust Major").fill(thrust.x(), thrust.y());
- aida.cloud2D(prefix+"/Correl Thrust Oblateness").fill(thrust.x(), ob);
aida.cloud1D(prefix+"/Thrust Axis CosTheta").fill(VecOp.cosTheta(ThrustAxis));
- aida.cloud2D(prefix+"/Correl Thrust vs Thrust Axis CosT"
- ).fill(thrust.x(), VecOp.cosTheta(ThrustAxis));
}
private void makeJetPlots(EventHeader event, String prefix, String decay,String evt) {
@@ -291,6 +405,29 @@
return rv;
}
+ // return true if event makes the cuts
+ private boolean evtShapeCut(EventShape es, double en, double min_en,
+ double min_thr, double max_thr,
+ double min_maj, double max_maj,
+ double min_min, double max_min,
+ double min_obl, double max_obl) {
+ Hep3Vector thrust_vec = es.thrust();
+ double thrust = thrust_vec.x();
+ double major = thrust_vec.y();
+ double minor = thrust_vec.z();
+ double oblate = es.oblateness();
+
+ boolean rv = false;
+ if ( (thrust > min_thr) && (thrust < max_thr) &&
+ (major > min_maj) && (major < max_maj) &&
+ (minor > min_min) && (minor < max_min) &&
+ (oblate > min_obl) && (oblate < max_obl) &&
+ (en > min_en)) {
+ rv = true;
+ }
+ return rv;
+ }
+
private List<MCParticle> cosThetaCut(List<MCParticle> mcpl, double ctcut) {
List<MCParticle> rv = new ArrayList<MCParticle>();
for (MCParticle p:mcpl) {
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -u -r1.1 -r1.2
--- MCHiggsDriver.java 8 Feb 2013 16:43:36 -0000 1.1
+++ MCHiggsDriver.java 20 Mar 2014 20:24:46 -0000 1.2
@@ -1,4 +1,4 @@
-package org.lcsim.mcd.analysis;
+//package org.lcsim.mcd.analysis;
import java.util.List;
import java.util.ArrayList;
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -u -r1.1 -r1.2
--- MCFastClusterAnalysis.java 19 Dec 2012 17:12:43 -0000 1.1
+++ MCFastClusterAnalysis.java 20 Mar 2014 20:24:46 -0000 1.2
@@ -1,4 +1,4 @@
-package org.lcsim.mcd.analysis;
+//package org.lcsim.mcd.analysis;
import hep.aida.ICloud1D;
import java.util.List;
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -u -r1.1 -r1.2
--- MCWTaggingDriver.java 8 Feb 2013 16:44:39 -0000 1.1
+++ MCWTaggingDriver.java 20 Mar 2014 20:24:46 -0000 1.2
@@ -1,11 +1,8 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-package org.lcsim.mcd.analysis;
+//package org.lcsim.mcd.analysis;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
+import hep.physics.jet.EventShape;
import java.util.ArrayList;
import java.util.List;
import org.lcsim.event.EventHeader;
@@ -24,7 +21,7 @@
*
*/
-enum DecayChannel {
+enum WDecayChannel {
UDS, BBAR, CCBAR, TTBAR, // quarks
GLUGLU, GAMGAM, WW, ZZ, ZGAM, // bosons
TAUTAUBAR, // leptons
@@ -41,6 +38,7 @@
// Array of W decay products by PDGID <-> index
+ // 0-index is total
int[] W_DECAY_COUNTS = new int[20];
private AIDA aida = AIDA.defaultInstance();
@@ -48,9 +46,9 @@
public void process(EventHeader event) {
this.processChildren(event);
MCParticle Higgs = findHiggs(event);
- DecayChannel dChan = getMCTag(Higgs);
+ WDecayChannel dChan = getMCTag(Higgs);
- if (dChan == DecayChannel.WW) {
+ if (dChan == WDecayChannel.WW) {
List<MCParticle> HiggsChildren = new ArrayList<MCParticle>(Higgs.getDaughters());
for (MCParticle p:Higgs.getDaughters()) {
if (Math.abs(p.getPDGID()) != 24) {
@@ -62,17 +60,100 @@
MCParticle Wp = HiggsChildren.get(0);
MCParticle Wm = HiggsChildren.get(1);
+ aida.cloud1D("W masses").fill(Wp.getMass());
+ aida.cloud1D("W masses").fill(Wm.getMass());
+
+ for(MCParticle p:Wp.getDaughters()){
+ if (p.getPDGID() != 24) {
+ aida.cloud1D("Wp decay channels").fill(p.getPDGID());
+ }
+ }
+ for(MCParticle p:Wm.getDaughters()){
+ if (p.getPDGID() != -24) {
+ aida.cloud1D("Wm decay channels").fill(p.getPDGID());
+ }
+ }
+
List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class,"Jets");
+ aida.cloud1D("Num jets (total)").fill(jets.size());
- List<ReconstructedParticle> Wpjets = dijetRPMatch(jets,Wp);
- List<ReconstructedParticle> Wmjets = dijetRPMatch(jets,Wm);
-
- List<LCRelation> lcrlist = new ArrayList<LCRelation>();
- lcrlist.add(new MyLCRelation((List<ReconstructedParticle>) Wpjets, (MCParticle) Wp));
- lcrlist.add(new MyLCRelation((List<ReconstructedParticle>) Wmjets, (MCParticle) Wm));
- event.put("DijetToMCP", lcrlist, ReconstructedParticle.class,0);
+ if (jets.size() >= 2) {
+ List<ReconstructedParticle> Wpjets = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> Wmjets = new ArrayList<ReconstructedParticle>();
+ if (Wp.getMass() > Wm.getMass()) {
+ if (Math.abs(Wp.getDaughters().get(0).getPDGID()) < 10) {
+ Wpjets = jetRPMatch(jets,Wp);
+ aida.cloud1D("Real W Masses").fill(Wp.getMass());
+ }
+ }
+ else {
+ if (Math.abs(Wm.getDaughters().get(0).getPDGID()) < 10) {
+ Wmjets = jetRPMatch(jets,Wm);
+ aida.cloud1D("Real W Masses").fill(Wm.getMass());
+ }
+ }
+
+ List<LCRelation> lcrlist = new ArrayList<LCRelation>();
+ lcrlist.add(new MyLCRelation((List<ReconstructedParticle>) Wpjets, (MCParticle) Wp));
+ lcrlist.add(new MyLCRelation((List<ReconstructedParticle>) Wmjets, (MCParticle) Wm));
+ event.put("jetToMCP", lcrlist, LCRelation.class,0);
+ }
}
}
+
+ // Given a list of jets and a particle, find the combination of jets that best reconstruct the particle's mass
+ public List<ReconstructedParticle> jetRPMatch(List<ReconstructedParticle> jets, MCParticle p) {
+ List<ReconstructedParticle> rvlist = new ArrayList<ReconstructedParticle>();
+ List<List<ReconstructedParticle>> jet_combos = new ArrayList<List<ReconstructedParticle>>();
+
+ for(ReconstructedParticle jet:jets) {
+ List<ReconstructedParticle> tmp = new ArrayList<ReconstructedParticle>();
+ tmp.add(jet);
+ jet_combos.add(tmp);
+ }
+ List<List<ReconstructedParticle>> new_jet_combos = jetRPMatchHelper(jet_combos);
+ double closest_mass = 1000000.0;
+ for (List<ReconstructedParticle> jetset:new_jet_combos) {
+ double inv_mass_sum = 0;
+ for (ReconstructedParticle jet:jetset) {
+ inv_mass_sum += Math.sqrt(
+ jet.getEnergy()*jet.getEnergy()
+ - jet.getMomentum().magnitudeSquared());
+ }
+ aida.cloud1D("All jet combo masses").fill(inv_mass_sum);
+ if (Math.abs(p.getMass() - inv_mass_sum) < Math.abs(p.getMass() - closest_mass)){
+ rvlist = jetset;
+ closest_mass = inv_mass_sum;
+ }
+ }
+ aida.cloud1D("Best jet combo masses").fill(closest_mass);
+ aida.cloud1D("Best jet combo differences").fill(p.getMass() - closest_mass);
+ aida.cloud1D("Best jet combo num jets").fill(rvlist.size());
+ aida.cloud1D("Num total jets - combo jets").fill(rvlist.size()-jets.size());
+
+ return rvlist;
+ }
+ // Returns a list with every possible combination of RP's in the original list,
+ // provided that the original list items each contain one RP
+ // Uses head recursion, which is okay since there should be fewer than 10 jets.
+ public List<List<ReconstructedParticle>> jetRPMatchHelper(List<List<ReconstructedParticle>> jet_combos) {
+ if (jet_combos.size() == 1) {
+ return jet_combos;
+ }
+ else {
+ List<List<ReconstructedParticle>> sublist = new ArrayList<List<ReconstructedParticle>>(jet_combos);
+ List<ReconstructedParticle> first = new ArrayList<ReconstructedParticle>(sublist.remove(0));
+ List<List<ReconstructedParticle>> newlist = jetRPMatchHelper(sublist);
+ int nl_size = newlist.size();
+ for(int i=0; i<nl_size; i++) {
+ List<ReconstructedParticle> tmplist = new ArrayList<ReconstructedParticle>(newlist.get(i));
+ tmplist.add(first.get(0));
+ newlist.add(tmplist);
+ }
+ newlist.add(first);
+ return newlist;
+ }
+ }
// Given a list of jets and a particle, find the two jets that best reconstruct the particle
public List<ReconstructedParticle> dijetRPMatch(List<ReconstructedParticle> jets, MCParticle p){
@@ -87,9 +168,9 @@
double p1Charge = p.getCharge();
double j1Charge = jet1.getCharge();
double j2Charge = jet2.getCharge();
- if (j1Charge + j2Charge != p1Charge) {
- break;
- }
+// if (j1Charge + j2Charge != p1Charge) {
+// break;
+// }
double mj1 = jet1.getMass();
double mj2 = jet2.getMass();
double ej1 = jet1.getEnergy();
@@ -107,10 +188,16 @@
rvlist.add(jet2);
best_mass_diff = diff_dj_mass;
}
+
+ aida.cloud1D("All dijet combination masses").fill(dijet_inv_mass);
}
}
- aida.cloud1D("W-matched dijet masses").fill(dijet_inv_mass);
- aida.cloud1D("W-matched dijet mass minus MC W mass").fill(best_mass_diff);
+
+ if(best_mass_diff < 10000000.0) {
+ aida.cloud1D("W-matched dijet mass minus MC W mass").fill(best_mass_diff);
+ aida.cloud1D("W-matched dijet masses").fill(dijet_inv_mass);
+ }
+
return rvlist;
}
@@ -126,11 +213,11 @@
}
// Determine the type of higgs decay by identifying its children by PDGID
- // Return as DecayChannel enum
- public DecayChannel getMCTag(MCParticle higgs) {
+ // Return as WDecayChannel enum
+ public WDecayChannel getMCTag(MCParticle higgs) {
List<MCParticle> HiggsChildren = higgs.getDaughters();
- DecayChannel chan = DecayChannel.UNKNOWN;
+ WDecayChannel chan = WDecayChannel.UNKNOWN;
if (!HiggsChildren.isEmpty()) {
int pdgid1 = 0;
int pdgid2 = 0;
@@ -155,37 +242,37 @@
aida.cloud1D("Decay Channels").fill(Math.abs(pdgid1));
switch (Math.abs(pdgid1)) {
case 1 :
- chan = DecayChannel.UDS;
+ chan = WDecayChannel.UDS;
break;
case 2 :
- chan = DecayChannel.UDS;
+ chan = WDecayChannel.UDS;
break;
case 3 :
- chan = DecayChannel.UDS;
+ chan = WDecayChannel.UDS;
break;
case 4 :
- chan = DecayChannel.CCBAR;
+ chan = WDecayChannel.CCBAR;
break;
case 5 :
- chan = DecayChannel.BBAR;
+ chan = WDecayChannel.BBAR;
break;
case 6 :
- chan = DecayChannel.TTBAR;
+ chan = WDecayChannel.TTBAR;
break;
case 21 :
- chan = DecayChannel.GLUGLU;
+ chan = WDecayChannel.GLUGLU;
break;
case 22 :
- chan = DecayChannel.GAMGAM;
+ chan = WDecayChannel.GAMGAM;
break;
case 23 :
- chan = DecayChannel.ZZ;
+ chan = WDecayChannel.ZZ;
break;
case 24 :
- chan = DecayChannel.WW;
+ chan = WDecayChannel.WW;
break;
default :
- chan = DecayChannel.OTHER;
+ chan = WDecayChannel.OTHER;
break;
}
}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -N MCHiggs.java
--- MCHiggs.java 8 Feb 2013 16:39:52 -0000 1.10
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,133 +0,0 @@
-package org.lcsim.mcd.analysis;
-
-import java.util.List;
-import java.util.ArrayList;
-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 hep.physics.vec.Hep3Vector;
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.VecOp;
-import org.lcsim.mc.fast.MCFast;
-
-/**
- *
- * @author agias
- */
-public class MCHiggs extends Driver{
-
- public MCHiggs()
- {
- add(new MCFast());
- }
-
- private AIDA aida = AIDA.defaultInstance();
-
- protected String mcParticlesSkimmedName = "MCParticle";
-
- public void process(EventHeader event) {
- List<MCParticle> particles = event.get(MCParticle.class,mcParticlesSkimmedName);
- List<MCParticle> children = new ArrayList<MCParticle>();
-
- boolean hfound = false;
- for (MCParticle particle : particles) {
- // 25 == h0/H01
- if ((particle.getPDGID()==25)&&(!hfound)) {
- children.addAll(particle.getDaughters());
- hfound=true;
- }
- }
- MCParticle p1 = children.get(0);
- MCParticle p2 = children.get(0);
-
- for (MCParticle particle : children) {
- /*
- * Only count positive PDGID's because negatives
- * are just antiparticles
- * 4 == c
- * 5 == b
- * 15 == tau-
- * 21 == gluon
- * 22 == gamma
- * 23 == Zo
- * 24 == W
- * 25 == h0/H01
- */
- if ((particle.getPDGID() > 0) && (particle.getPDGID() != 25)) {
- aida.cloud1D("decay modes").fill(particle.getPDGID());
- p1 = particle;
- }
- if ((particle.getPDGID() < 0)) {
- p2 = particle;
- }
- }
-
- if (p1.getPDGID() == 5) {
- Hep3Vector p = p1.getMomentum();
- Hep3Vector p_bar = p2.getMomentum();
-
- double mass1 = p1.getMass();
- double en1 = p1.getEnergy();
- double mass2 = p2.getMass();
- double en2 = p2.getEnergy();
-
- Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
-
- double p1_m = p.magnitude();
- double p2_m = p_bar.magnitude();
-
- double pt = VecOp.dot(p, z_unit);
- double pt_bar = VecOp.dot(p_bar, z_unit);
-
- aida.cloud1D("b momentum").fill(p1_m);
- aida.cloud1D("b_bar momentum").fill(p2_m);
- aida.histogram1D("b momentum hist",100,62.3108,62.3126).fill(p1_m);
- aida.histogram1D("b_bar momentum hist",100,62.3,62.4).fill(p2_m);
-
- aida.cloud1D("b transverse momentum").fill(pt);
- aida.cloud1D("b_bar transverse momentum").fill(pt_bar);
- aida.histogram1D("b transverse momentum hist", 100,0,63).fill(pt);
- aida.histogram1D("b_bar transverse momentum hist", 100,0,63).fill(pt_bar);
-
-
- double inv_mass = Math.sqrt(mass1*mass1 + mass2*mass2 + 2.*(en1*en2-VecOp.dot(p, p_bar)));
- aida.cloud1D("H invariant mass").fill(inv_mass);
- aida.histogram1D("H inv mass hist",100,124.990,124.995).fill(inv_mass);
- aida.histogram1D("H inv mass hist big",200,100,130).fill(inv_mass);
-
- aida.cloud1D("b phi").fill(VecOp.phi(p));
- aida.cloud1D("bbar phi").fill(VecOp.phi(p_bar));
- aida.cloud1D("b cos theta").fill(VecOp.cosTheta(p));
- aida.cloud1D("bbar cos theta").fill(VecOp.cosTheta(p_bar));
- }
- if (p1.getPDGID() == 22) {
- double mass1 = p1.getMass();
- double ener1 = p1.getEnergy();
- double mass2 = p2.getMass();
- double ener2 = p2.getEnergy();
-
- Hep3Vector pt1 = p1.getMomentum();
- Hep3Vector pt2 = p2.getMomentum();
- double p1_m = pt1.magnitude();
- double p2_m = pt2.magnitude();
-
- Hep3Vector z_unit = new BasicHep3Vector(0,0,1.);
-
- double p1trans = VecOp.dot(pt1, z_unit);
- double p2trans = VecOp.dot(pt2,z_unit);
-
-
- double inv_mass = Math.sqrt(2*(ener1*ener2 - VecOp.dot(pt1, pt2)));
-
- aida.cloud1D("gamma1 momentum").fill(p1_m);
- aida.cloud1D("gamma2 momentum").fill(p2_m);
- aida.cloud1D("gamma1 transverse momentum").fill(p1trans);
- aida.cloud1D("gamma2 transverse momentum").fill(p2trans);
- aida.cloud1D("gamma higgs inv mass").fill(inv_mass);
- }
- aida.cloud1D("num children").fill(children.size());
- }
-
-}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N DRCalibrationDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DRCalibrationDriver.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,110 @@
+/*
+ * * @author wenzel
+ * This Driver is used to obtain the dualreadout correction of a dual readout calorimeter.
+ * The Determination of the energy scale of the calorimeter response is done in a previous step
+ * using electrons (see: ElectronCalibrationDriver)
+ * first it runs the Digitiser and then the Dual readout calibration module.
+ */
+package org.lcsim.mcd.analysis.DRCorrection;
+
+//import org.lcsim.contrib.HansWenzel.DualCorrection.DigiSim.CSClusdMJetDriver;
+//import org.lcsim.mcd.analysis.DRCorrection.DualCorrection;
+import org.lcsim.util.Driver;
+
+public class DRCalibrationDriver extends Driver {
+ //
+ // default detector configuration:
+ //
+
+ private String Detector = "mcdrcal01";
+ private String Material = "PbW04";
+ private Double Density = 8.28;
+ private Double rindex = 2.21;
+ private String CollectionName = "NoCuts";
+ private Double CerenkovThres = 0.02;
+ private Double IonizationThres = 0.02;
+ private String PhysicsList = "FTFP_BERT";
+ String AIDAFile = null;
+ String file_name = null;
+ DualCorrection dual;
+
+ public DRCalibrationDriver() {
+ dual = new DualCorrection();
+ add(dual);
+ }
+
+ @Override
+ public void startOfData() {
+ System.out.println("DRCalibrationDriver:startOfData");
+ System.out.println(AIDAFile);
+ if (AIDAFile != null) {
+ dual.setMyAIDAFilename(AIDAFile);
+ } else {
+ System.err.println("DRCalibrationDriver: AIDAFile variable must be set");
+ System.exit(1); // exit if variable is not set
+ }
+ if (file_name != null) {
+ dual.setMyFilename(file_name);
+ } else {
+ System.err.println("DRCalibrationDriver: file_name variable must be set");
+ System.exit(1); // exit if variable is not set
+ }
+ dual.setMyPhysicsList(PhysicsList);
+ dual.setMyDetector(Detector);
+ dual.setMyMaterial(Material);
+ dual.setMyDensity(Density);
+ dual.setMyrindex(rindex);
+ dual.setMyCollectionName(CollectionName);
+ dual.setMyCerenkovThres(CerenkovThres);
+ dual.setMyIonizationThres(IonizationThres);
+ dual.startOfData();
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("DRCalibrationDriver:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+
+ public void setMyFilename(String file_name) {
+ System.out.println("DRCalibrationDriver:setMyFilename");
+ this.file_name = file_name;
+ System.out.println(file_name);
+ }
+
+ public void setMyPhysicsList(String list) {
+ System.out.println("setMyPhysicsList");
+ this.PhysicsList = list;
+ System.out.println(list);
+ }
+
+ public void setMyDetector(String Detector) {
+ this.Detector = Detector;
+ }
+
+ public void setMyMaterial(String Material) {
+ this.Material = Material;
+ }
+
+ public void setMyDensity(double Density) {
+ System.out.println("DRCalibrationDriver:setMyDensity");
+ this.Density = Density;
+ System.out.println(Density);
+ }
+
+ public void setMyrindex(double rindex) {
+ this.rindex = rindex;
+ }
+
+ public void setMyCollectionName(String CollectionName) {
+ this.CollectionName = CollectionName;
+ }
+
+ public void setMyCerenkovThres(double CerenkovThres) {
+ this.CerenkovThres = CerenkovThres;
+ }
+
+ public void setMyIonizationThres(double IonizationThres) {
+ this.IonizationThres = IonizationThres;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N Resolution.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Resolution.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,391 @@
+package org.lcsim.mcd.analysis.DRCorrection;
+
+/**
+ *
+ * @author wenzel
+ * @date
+ */
+import hep.aida.ICloud1D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import java.io.BufferedWriter;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.aida.AIDA;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+
+public class Resolution extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ String AIDAFile;
+ String file_name;
+ FileWriter fstream;
+ BufferedWriter out;
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsFactory;
+ IFunction gauss;
+ static boolean first;
+ static boolean firstEvent;
+ static double E_in;
+ Integer Ein;
+ Integer Ein_prev;
+ double E_kin;
+ String Particlename;
+ ICloud1D Edep;
+ ICloud1D ECeren;
+ ICloud1D Edep_cor;
+ ICloud1D Edep_dcor;
+ ICloud1D Eceren_cor;
+ double nsigmas=5.0;
+ int nbins=100;
+ IDataPointSet[] dps_Correctede = null;
+ IDataPointSet[] dps_Resolution = null;
+ IDataPointSet dps_ECorrFrac = null;
+ IDataPointSet dps_CCorrFrac = null;
+ IDataPointSet[] dps_DCorrFrac = null;
+ int point_Correctede[];
+ double[] result;
+ double errors[];
+ String[] Fitters = {"Chi2", "leastsquares"};
+ // available are "Chi2", "leastsquares", "bml", "cleverchi2"
+ // but "bml", "cleverchi2" seem to produce rubish
+
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ IFunction E_cor;
+ IFunction C_cor;
+ IFunction D_cor;
+ double xval[] = {10.};
+
+ double t_thresh;
+ double e_thresh;
+
+ //
+ // default detector configuration:
+ //
+ private String Detector = "mcdrcal01";
+ private String Material = "PbW04";
+ private Double Density = 8.28;
+ private Double rindex = 2.21;
+ private String CollectionName = "NoCuts";
+ private Double CerenkovThres = 0.02;
+ private Double IonizationThres = 0.02;
+ private String PhysicsList = "FTFP_BERT";
+
+ @Override
+ protected void startOfData() {
+ System.out.println("Start of Data:");
+ t_thresh = -1.0;
+ e_thresh = -1.0;
+ dps_Correctede = new IDataPointSet[Fitters.length];
+ dps_Resolution = new IDataPointSet[Fitters.length];
+ dps_DCorrFrac = new IDataPointSet[Fitters.length];
+ point_Correctede = new int[Fitters.length];
+ fitFactory = aida.analysisFactory().createFitFactory();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ dps_ECorrFrac = dpsFactory.create("dps_ECorrFrac_", 2);
+ dps_CCorrFrac = dpsFactory.create("dps_CCorrFrac_", 2);
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_Correctede_" + Fitters[i];
+ point_Correctede[i] = 0;
+ dps_Correctede[i] = dpsFactory.create(dpsname, "electron response mean", 2);
+ dps_Resolution[i] = dpsFactory.create("dps_Resolution_"+Fitters[i], 2);
+ dps_DCorrFrac[i] = dpsFactory.create("dps_DCorrFrac_"+Fitters[i], 2);
+ }
+ DRFunctionFactory mapoff = DRFunctionFactory.getInstance();
+ DetectorConfiguration dc = new DetectorConfiguration(Detector, Material,
+ Density, rindex,
+ CollectionName, CerenkovThres,
+ IonizationThres, PhysicsList);
+ if (mapoff.checkdc(dc)) {
+ E_cor = mapoff.getecFunction(dc);
+ C_cor = mapoff.getccFunction(dc);
+ D_cor = mapoff.getdcFunction(dc);
+ } else {
+ System.exit(0);
+ }
+
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ fstream = null;
+ try {
+ fstream = new FileWriter(file_name);
+ } catch (IOException ex) {
+ Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ out = new BufferedWriter(fstream);
+ }
+
+ @Override
+ protected void process(EventHeader event) {
+// System.out.println("Process event #"+String.valueOf(event.getEventNumber()));
+ E_in = 0.0;
+ E_kin = 0.0;
+ String DirName = null;
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ if (particle.getParents().isEmpty()) {
+// if (particle.getGeneratorStatus() == 0) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+ if (particle.getProductionTime() == 0.0) {
+ Particlename = particle.getType().getName();
+ }
+ }
+ break;
+ }
+ Ein = (int) Math.floor(E_kin + 0.5d);
+ String E_str = Ein.toString();
+// System.out.println(E_str+" GeV "+Particlename);
+ DirName = Particlename.concat(E_str);
+ if (Ein != Ein_prev) {
+ if (first) {
+ first = false;
+ } else {
+ System.out.println("E_in: " + E_in);
+ System.out.println("Ein_prev: " + Ein_prev);
+ convertandfit();
+ }
+ Ein_prev = Ein;
+ System.out.println("First Event:");
+ System.out.println("E_in: " + E_in);
+ System.out.println("E_kin: " + E_kin);
+ System.out.println("Name: " + Particlename);
+ aida.tree().cd("/");
+ aida.tree().mkdir(DirName);
+ aida.tree().cd(DirName);
+ Edep = aida.histogramFactory().createCloud1D("Edep", "uncorrected Energy Cloud", 100000, "autoConvert = false");
+ ECeren = aida.histogramFactory().createCloud1D("ECeren", "uncorrected Cerenkov Cloud", 100000, "autoConvert = false");
+ Edep = aida.histogramFactory().createCloud1D("Edep", "Energy Cloud", 100000, "autoConvert = false");
+ Edep_cor = aida.histogramFactory().createCloud1D("Edep_cor", "corrected Energy Cloud", 100000, "autoConvert = false");
+ Edep_dcor = aida.histogramFactory().createCloud1D("Edep_dcor", "dual readout corrected Energy Cloud", 100000, "autoConvert = false");
+ Eceren_cor = aida.histogramFactory().createCloud1D("Eceren_cor", "corrected Cerenkov Cloud", 100000, "autoConvert = false");
+ System.out.println("DirName: " + DirName);
+ }
+
+ List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+ double sumEEdep = 0.0;
+ double sumECeren = 0.0;
+ for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+ String colName = event.getMetaData(simCalorimeterHits).getName();
+ if (colName.contains("Edep")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (t_thresh > 0.0 || e_thresh > 0.0) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - calorimeterHit.getPositionVec().magnitude()/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < t_thresh) && (edep > e_thresh) ) {
+ sumEEdep += edep;
+ }
+ }
+ }
+ else {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ }
+ } // end if Edep
+ if (colName.contains("Opti")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (t_thresh > 0.0 || e_thresh > 0.0) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - calorimeterHit.getPositionVec().magnitude()/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < t_thresh) && (edep > e_thresh) ) {
+ sumECeren += edep;
+ }
+ }
+ }
+ else {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ }
+ }
+ }
+
+ Edep.fill(sumEEdep);
+ xval[0] = sumEEdep;
+ double sumEEdep_cor = E_cor.value(xval);
+ Edep_cor.fill(sumEEdep_cor);
+ ECeren.fill(sumECeren);
+ xval[0] = sumECeren;
+ double sumECeren_cor = C_cor.value(xval);
+ Eceren_cor.fill(sumECeren_cor);
+ double ratio = sumECeren_cor / sumEEdep_cor;
+ double fraction = sumEEdep_cor / E_in;
+ double cfac;
+ if (ratio < 1.) {
+ xval[0] = ratio;
+ cfac = D_cor.value(xval);
+ } else {
+ cfac = 1.0;
+ }
+ aida.cloud2D("/alex/Ecor over Ein vs Ccor over Ein").fill(ratio, fraction/cfac);
+ IDataPoint dp_ECorrFrac = dps_ECorrFrac.addPoint();
+ Edep_dcor.fill(sumEEdep_cor / cfac);
+ }
+
+ @Override
+ protected void endOfData() {
+ System.out.println("End of Data:");
+ convertandfit();
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ @Override
+ protected void resume() {
+ System.out.println("resume:");
+ firstEvent = true;
+ aida.cloud1D("c_Edep_energy").reset();
+ }
+
+ protected void convertandfit() {
+ System.out.println("convert and fit:");
+ IHistogram1D Edep_dcor_conv;
+
+ if (Edep_dcor.isConverted()) {
+ Edep_dcor_conv = Edep_dcor.histogram();
+ } else {
+ double meanc = Edep_dcor.mean();
+ double rmsc = Edep_dcor.rms();
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep_dcor.setConversionParameters(nbins, minx, maxx);
+ Edep_dcor.convertToHistogram();
+ Edep_dcor_conv = Edep_dcor.histogram();
+ }
+ gauss.setParameter("amplitude", Edep_dcor_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_dcor_conv.mean());
+ gauss.setParameter("sigma", Edep_dcor_conv.rms());
+ System.out.print(Particlename + " " + Ein_prev + " GeV\n");
+ System.out.print(Edep_dcor_conv.maxBinHeight() + "," + Edep_dcor_conv.mean() + "," + Edep_dcor_conv.rms() + "\n");
+ try {
+ out.write(Particlename + " " + Ein_prev + " GeV\n");
+ out.write(Edep_dcor_conv.maxBinHeight() + "," + Edep_dcor_conv.mean() + "," + Edep_dcor_conv.rms() + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+ for (int i = 0; i < Fitters.length; i++) {
+ System.out.println("Fitter: " + Fitters[i]);
+ gauss.setParameter("amplitude", Edep_dcor_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_dcor_conv.mean());
+ gauss.setParameter("sigma", Edep_dcor_conv.rms());
+ jminuit = fitFactory.createFitter(Fitters[i], "jminuit");
+ jminuitResult = jminuit.fit(Edep_dcor_conv, gauss);
+ System.out.println("jminuit " + Fitters[i] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ double[] errsplus = jminuitResult.errorsPlus();
+ double[] errsminu = jminuitResult.errorsMinus();
+ String functionname = "Corrected fitted gauss " + Fitters[i];
+ System.out.println(functionname);
+ functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+ try {
+ out.write(result[0] + "," + result[1] + "," + result[2] + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ IDataPoint dp_Correctede = dps_Correctede[i].addPoint();
+ dp_Correctede.coordinate(0).setValue(result[1]);
+ dp_Correctede.coordinate(0).setErrorPlus(Math.abs(result[2]));
+ dp_Correctede.coordinate(0).setErrorMinus(Math.abs(result[2]));
+ dp_Correctede.coordinate(1).setValue(Ein_prev);
+ dp_Correctede.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+ dp_Correctede.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+ point_Correctede[i]++;
+ IDataPoint dp_Resolution = dps_Resolution[i].addPoint();
+ double sigma = Math.abs(result[2]);
+ dp_Resolution.coordinate(0).setValue(1.0/Math.sqrt(Ein_prev));
+ dp_Resolution.coordinate(1).setValue(Math.abs(sigma)/Ein_prev);
+ dp_Resolution.coordinate(1).setErrorMinus(errsplus[2]);
+ dp_Resolution.coordinate(1).setErrorPlus(errsminu[2]);
+ IDataPoint dp_DCorrFrac = dps_DCorrFrac[i].addPoint();
+ double mean = result[1];
+ dp_DCorrFrac.coordinate(0).setValue(Ein_prev);
+ dp_DCorrFrac.coordinate(1).setValue(mean/Ein_prev);
+ dp_DCorrFrac.coordinate(1).setErrorMinus(sigma);
+ dp_DCorrFrac.coordinate(1).setErrorPlus(sigma);
+ }
+
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setMyFilename(String filename) {
+ file_name = filename;
+ }
+
+ public void setMyAIDAFilename(String aidafile) {
+ AIDAFile = aidafile;
+ }
+ public void setMyPhysicsList(String list) {
+ System.out.println("setMyPhysicsList");
+ this.PhysicsList = list;
+ System.out.println(list);
+ }
+
+ public void setMyDetector(String Detector) {
+ this.Detector = Detector;
+ }
+
+ public void setMyMaterial(String Material) {
+ this.Material = Material;
+ }
+
+ public void setMyDensity(Double Density) {
+ this.Density = Density;
+ }
+
+ public void setMyrindex(Double rindex) {
+ this.rindex = rindex;
+ }
+
+ public void setMyCollectionName(String CollectionName) {
+ this.CollectionName = CollectionName;
+ }
+
+ public void setMyCerenkovThres(Double CerenkovThres) {
+ this.CerenkovThres = CerenkovThres;
+ }
+
+ public void setMyIonizationThres(Double IonizationThres) {
+ this.IonizationThres = IonizationThres;
+ }
+ public void setMynsigmas(double nsigmas) {
+ this.nsigmas = nsigmas;
+ }
+ public void setMynbins(int nbins) {
+ this.nbins = nbins;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N pion_lowen_calib_steering.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ pion_lowen_calib_steering.xml 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,37 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+ </classpath>
+
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-0.2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-0.5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-0.75gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-2gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+
+ <execute>
+ <driver name="DRCalibrationDriver"/>
+ </execute>
+ <drivers>
+ <driver name="DRCalibrationDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.DRCalibrationDriver">
+ <myAIDAFilename>pipl_calibration-LowEnNoCuts.aida</myAIDAFilename>
+ <myFilename>pipl_calibration-LowEnNoCuts.txt</myFilename>
+ <myPhysicsList>FTFP_BERT</myPhysicsList>
+ <myDetector>mcdrcal01</myDetector>
+ <myMaterial>PbW04</myMaterial>
+ <myDensity>8.28</myDensity>
+ <myrindex>2.21</myrindex>
+ <myCollectionName>LowEnNoCuts</myCollectionName>
+ <myCerenkovThres>0.02</myCerenkovThres>
+ <myIonizationThres>0.02</myIonizationThres>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N DetectorConfiguration.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DetectorConfiguration.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,110 @@
+package org.lcsim.mcd.analysis.DRCorrection;
+
+/**
+ *
+ * @author wenzel
+ *
+ * This class is used to describe a specifc detector configuration as it is used to study
+ * properties of dual read out calorimetry.
+ * Calibration constants like energy scale for the ionization and cerenkov response as well
+ * as dual readout correction functions are associated with a unique detector configuration
+ * nd have to be obtained for different materials and physics lists.
+ * All correction function can be accessed via the DRFunctionFactory class.
+ */
+public class DetectorConfiguration {
+
+ private String Detector;
+ private String Material;
+ private Double Density;
+ private Double rindex;
+ private String CollectionName;
+ private Double CerenkovThres;
+ private Double IonizationThres;
+ private String PhysicsList;
+
+ public DetectorConfiguration(String d, String m, Double den, Double r, String n, Double Ct, Double It, String pl) {
+ Detector = d;
+ Material = m;
+ Density = den;
+ rindex = r;
+ CollectionName = n;
+ CerenkovThres = Ct;
+ IonizationThres = It;
+ PhysicsList = pl;
+ }
+
+ public void print() {
+ System.out.println("Detector: " + Detector);
+ System.out.println("Material: " + Material);
+ System.out.println("Density: " + Density);
+ System.out.println("Refraction Index: " + rindex);
+ System.out.println("Collection Name: " + CollectionName);
+ System.out.println("Cerenkov Threshold: " + CerenkovThres);
+ System.out.println("Ionization Theshold: " + IonizationThres);
+ System.out.println("Physics List: " + PhysicsList);
+ }
+
+ @Override
+ public boolean equals(Object o) {
+ if ((o instanceof DetectorConfiguration) &&
+ (((DetectorConfiguration) o).getDetector().equals(Detector)) &&
+ (((DetectorConfiguration) o).getMaterial().equals(Material)) &&
+ (((DetectorConfiguration) o).getDensity().compareTo(Density) == 0) &&
+ (((DetectorConfiguration) o).getrindex().compareTo(rindex) == 0) &&
+ (((DetectorConfiguration) o).getCollectionName().equals(CollectionName)) &&
+ (((DetectorConfiguration) o).getCerenkovThres().compareTo(CerenkovThres) == 0) &&
+ (((DetectorConfiguration) o).getIonizationThres().compareTo(IonizationThres) == 0) &&
+ (((DetectorConfiguration) o).getPhysicsList().equals(PhysicsList))) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ @Override
+ public int hashCode() {
+ int hash = 7;
+ hash = 53 * hash + (this.Detector != null ? this.Detector.hashCode() : 0);
+ hash = 53 * hash + (this.Material != null ? this.Material.hashCode() : 0);
+ hash = 53 * hash + (this.Density != null ? this.Density.hashCode() : 0);
+ hash = 53 * hash + (this.rindex != null ? this.rindex.hashCode() : 0);
+ hash = 53 * hash + (this.CollectionName != null ? this.CollectionName.hashCode() : 0);
+ hash = 53 * hash + (this.CerenkovThres != null ? this.CerenkovThres.hashCode() : 0);
+ hash = 53 * hash + (this.IonizationThres != null ? this.IonizationThres.hashCode() : 0);
+ hash = 53 * hash + (this.PhysicsList != null ? this.PhysicsList.hashCode() : 0);
+ return hash;
+ }
+
+ public String getDetector() {
+ return Detector;
+ }
+
+ public String getMaterial() {
+ return Material;
+ }
+
+ public Double getDensity() {
+ return Density;
+ }
+
+ public Double getrindex() {
+ return rindex;
+ }
+
+ public String getCollectionName() {
+ return CollectionName;
+ }
+
+ public Double getCerenkovThres() {
+ return CerenkovThres;
+ }
+
+ public Double getIonizationThres() {
+ return IonizationThres;
+ }
+
+ public String getPhysicsList() {
+ return PhysicsList;
+ }
+}
+
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N ElectronCalibrationDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ElectronCalibrationDriver.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,54 @@
+/*
+ * This Driver is used to calibrate the correction of a dual readout calorimeter to electrons.
+ * (Determining the energy scale of the calorimeter response.)
+ * first it runs the Digitiser and then the Electron calibration module.
+ */
+package org.lcsim.mcd.analysis.DRCorrection;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author wenzel
+ */
+public class ElectronCalibrationDriver extends Driver {
+
+ String AIDAFile = null;
+ String file_name = null;
+ ElectronCorrection elec;
+
+ public ElectronCalibrationDriver() {
+ elec = new ElectronCorrection();
+ add(elec);
+ }
+
+ @Override
+ public void startOfData() {
+ System.out.println("ElectronCalibrationDriver:startOfData");
+ if (AIDAFile != null) {
+ elec.setMyAIDAFilename(AIDAFile);
+ } else {
+ System.err.println("ElectronCalibrationDriver: AIDAFile variable must be set");
+ System.exit(1); // exit if variable is not set
+ }
+ if (file_name != null) {
+ elec.setMyFilename(file_name);
+ } else {
+ System.err.println("ElectronCalibrationDriver: file_name variable must be set");
+ System.exit(1); // exit if variable is not set
+ }
+
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("ElectronCalibrationDriver:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+
+ public void setMyFilename(String file_name) {
+ System.out.println("ElectronCalibrationDriver:setMyFilename");
+ this.file_name = file_name;
+ System.out.println(file_name);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N DualCorrection.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DualCorrection.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,507 @@
+package org.lcsim.mcd.analysis.DRCorrection;
+
+/*
+ * @author wenzel
+ * @date
+ * This Driver is used to obtain the dualreadout correction of a dual readout calorimeter.
+ * The Determination of the energy scale of the calorimeter response is done in a previous step
+ * using electrons (see: ElectronCalibrationDriver)
+ */
+import hep.aida.ICloud1D;
+import hep.aida.ICloud2D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IProfile1D;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.io.BufferedWriter;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.text.DecimalFormat;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.aida.AIDA;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+
+public class DualCorrection extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ String AIDAFile;
+ String file_name;
+ FileWriter fstream;
+ BufferedWriter out;
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsFactory;
+ boolean first;
+ boolean firstEvent;
+ double E_in;
+ Integer Ein;
+ Integer Ein_prev;
+ String Part_prev;
+ double E_kin;
+ String Particlename;
+ ICloud1D Edep;
+ ICloud1D Eceren;
+ ICloud1D Edep_cor;
+ ICloud1D Eceren_cor;
+ ICloud1D dlength;
+ ICloud2D c_efrac_ratio;
+ ICloud2D c_Ceren_vs_Edep;
+ ICloud1D[] slice_comb;
+ ICloud1D[] slice;
+ IHistogram1D[] conv_slice_comb;
+ IHistogram1D[] conv_slice;
+ IFunction gauss;
+ IProfile1D prof_combined;
+ IProfile1D ratio;
+ double nsigmas;
+ int nbins;
+ IDataPointSet[] dps_ratio;
+ int point_ratio[];
+ double[] result;
+ double errors[];
+ String[] Fitters = {"Chi2", "leastsquares"};
+ // Available are: "Chi2", "leastsquares", "bml", "cleverchi2";
+ // but bml and cleverchi2 results don't make too much sense
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ IFunction E_cor; // electron Ionisation correction function
+ IFunction C_cor; // electron cerenkov correction function
+ IFunction poly;
+ double xval[] = {10.};
+ double binwidth = 0.04;
+ int maxbin = (int) Math.rint(1.0 / binwidth);
+
+ double t_thresh;
+ double e_thresh;
+ //
+ // default detector configuration:
+ //
+ private String Detector = "mcdrcal01";
+ private String Material = "PbW04";
+ private Double Density = 8.28;
+ private Double rindex = 2.21;
+ private String CollectionName = "LowEnNoCuts";
+ private Double CerenkovThres = 0.02;
+ private Double IonizationThres = 0.02;
+ private String PhysicsList = "FTFP_BERT";
+
+ @Override
+ protected void startOfData() {
+ // if e_thresh > 0 then t_thresh MUST be > 0 as well.
+ t_thresh = -1.0;
+ e_thresh = -1.0;
+ System.out.println("DualCorrection:Start of Data");
+ fitFactory = aida.analysisFactory().createFitFactory();
+ dps_ratio = new IDataPointSet[Fitters.length];
+ point_ratio = new int[Fitters.length];
+ dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ poly = functionFactory.createFunctionByName("poly", "p4");
+ jminuit = fitFactory.createFitter("Chi2", "jminuit");
+ poly.setParameter("p0", 0.63);
+ poly.setParameter("p1", 0.54);
+ poly.setParameter("p2", -0.58);
+ poly.setParameter("p3", 0.5);
+ poly.setParameter("p4", 0.0);
+ DRFunctionFactory mapoff = DRFunctionFactory.getInstance();
+ System.out.println(Density);
+ DetectorConfiguration dc = new DetectorConfiguration(Detector, Material,
+ Density, rindex,
+ CollectionName, CerenkovThres,
+ IonizationThres, PhysicsList);
+ if (mapoff.checkdc(dc)) {
+ E_cor = mapoff.getecFunction(dc);
+ C_cor = mapoff.getccFunction(dc);
+ } else {
+ System.exit(0);
+ }
+// E_cor = functionFactory.createFunctionByName("ec_mcdrcal01_NoCuts_FTFP_BERT", "p1");
+// E_cor.setParameter("p0", 4.026e-3);
+// E_cor.setParameter("p1", 0.9998);
+// C_cor = functionFactory.createFunctionByName("cc_mcdrcal01_NoCuts_FTFP_BERT", "p1");
+// C_cor.setParameter("p0", 3.6605e-3);
+// C_cor.setParameter("p1", 1.87977e-5);
+ dlength = aida.histogramFactory().createCloud1D("dlength", "decay length combined ", 500000, "autoConvert = false");
+ prof_combined = aida.histogramFactory().createProfile1D("ratio_comb", "ratio combined", 30, 0.2, 1.);
+ c_efrac_ratio = aida.histogramFactory().createCloud2D("c_efrac_ratio", "c_efrac_ratio combined ", 500000, "autoConvert = false");
+ slice_comb = new ICloud1D[maxbin + 1];
+ conv_slice_comb = new IHistogram1D[maxbin + 1];
+ for (int ibin = 0; ibin < maxbin + 1; ibin++) {
+ String cloudy = "ratio_" + ibin;
+ slice_comb[ibin] = aida.histogramFactory().createCloud1D(cloudy, cloudy, 500000, "autoConvert = false");
+ }
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+ fstream = null;
+ try {
+ fstream = new FileWriter(file_name);
+ } catch (IOException ex) {
+ Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ out = new BufferedWriter(fstream);
+ }
+
+ @Override
+ protected void process(EventHeader event) {
+ E_in = 0.0;
+ E_kin = 0.0;
+ String DirName = null;
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ if (particle.getParents().isEmpty()) {
+// if (particle.getGeneratorStatus() == 0) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+// Hep3Vector end_point = particle.getEndPoint();
+// Hep3Vector start_point = particle.getOrigin();
+// Hep3Vector diff = VecOp.sub(end_point, start_point);
+// double decaylength = diff.magnitude();
+// dlength.fill(decaylength);
+// if (decaylength > 500.) {
+// return;
+// }
+ if (particle.getProductionTime() == 0.0) {
+ Particlename = particle.getType().getName();
+ }
+ }
+ break;
+ }
+ Ein = (int) Math.floor(E_kin + 0.5d);
+ String E_str = Ein.toString();
+ DirName = Particlename.concat(E_str);
+ if (Ein != Ein_prev || Particlename.equals(Part_prev) != true) {
+ if (first) {
+ first = false;
+ } else {
+ System.out.println("Process event: ");
+ fitprofile(aida.profile1D("ratio"), "dodo");
+ convertandfit(slice, conv_slice);
+ }
+ Ein_prev = Ein;
+ Part_prev = Particlename;
+ System.out.println("First Event:");
+ System.out.println("E_in: " + E_in);
+ System.out.println("E_kin: " + E_kin);
+ System.out.println("Name: " + Particlename);
+ aida.tree().cd("/");
+ aida.tree().mkdir(DirName);
+ aida.tree().cd(DirName);
+ Edep = aida.histogramFactory().createCloud1D("Edep", "Energy Cloud", 100000, "autoConvert = false");
+ Eceren = aida.histogramFactory().createCloud1D("Eceren", "Cerenkov Cloud", 100000, "autoConvert = false");
+ Edep_cor = aida.histogramFactory().createCloud1D("Edep_cor", "corrected Energy Cloud", 100000, "autoConvert = false");
+ Eceren_cor = aida.histogramFactory().createCloud1D("Eceren_cor", "corrected Cerenkov Cloud", 100000, "autoConvert = false");
+ c_Ceren_vs_Edep = aida.histogramFactory().createCloud2D("c_Ceren_vs_Edep", "ceren vs Edep", 100000, "autoConvert = false");
+ slice = new ICloud1D[maxbin + 1];
+ conv_slice = new IHistogram1D[maxbin + 1];
+ for (int ibin = 0; ibin < maxbin + 1; ibin++) {
+ String cloudy = "ratio_" + ibin;
+ slice[ibin] = aida.histogramFactory().createCloud1D(cloudy, cloudy, 500000, "autoConvert = false");
+ }
+ System.out.println("DirName: " + DirName);
+ }
+
+ List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+ double sumEEdep = 0.0;
+ double sumECeren = 0.0;
+ aida.profile1D("ratio", 30, 0.2, 1.);
+ for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+ String colName = event.getMetaData(simCalorimeterHits).getName();
+ if (colName.contains("Edep")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (t_thresh > 0.0 || e_thresh > 0.0) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - calorimeterHit.getPositionVec().magnitude()/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < t_thresh) && (edep > e_thresh) ) {
+ sumEEdep += edep;
+ }
+ }
+ }
+ else {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ }
+ } // end if Edep
+ if (colName.contains("Opti")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (t_thresh > 0.0 || e_thresh > 0.0) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - calorimeterHit.getPositionVec().magnitude()/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < t_thresh) && (edep > e_thresh) ) {
+ sumECeren += edep;
+ }
+ }
+ }
+ else {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ }
+ }
+ }
+
+ Edep.fill(sumEEdep);
+ xval[0] = sumEEdep;
+ double sumEEdep_cor = E_cor.value(xval);
+ Edep_cor.fill(sumEEdep_cor);
+ Eceren.fill(sumECeren);
+ xval[0] = sumECeren;
+ double sumECeren_cor = C_cor.value(xval);
+ Eceren_cor.fill(sumECeren_cor);
+ double ratio = sumECeren_cor / sumEEdep_cor;
+ double fraction = sumEEdep_cor / E_in;
+ aida.cloud1D("frac").fill(fraction);
+ aida.cloud1D("c_CerenEdep_ratio").fill(ratio);
+ c_Ceren_vs_Edep.fill(sumECeren_cor, sumEEdep_cor);
+ aida.cloud2D("c_efrac_ratio").fill(ratio, fraction);
+ aida.profile1D("ratio").fill(ratio, fraction);
+ xval[0] = ratio;
+ prof_combined.fill(ratio, fraction);
+ c_efrac_ratio.fill(ratio, fraction);
+ highratioplots(event,simCalorimeterHitCollections,ratio, fraction);
+ int bin = (int) Math.rint(ratio / binwidth);
+ if (bin < 0) {
+ bin = 0;
+ }
+ if (bin > maxbin - 1) {
+ bin = maxbin;
+ }
+ slice_comb[bin].fill(fraction);
+ slice[bin].fill(fraction);
+ }
+
+ @Override
+ protected void endOfData() {
+ System.out.println("DualCorrection: endOfData");
+ convertandfit(slice, conv_slice);
+ aida.tree().cd("/");
+ fitprofile(prof_combined, "profcombined");
+ convertandfit(slice_comb, conv_slice_comb);
+ try {
+ out.close();
+ } catch (IOException ex) {
+ Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ @Override
+ protected void resume() {
+ System.out.println("resume:");
+ firstEvent = true;
+ aida.cloud1D("c_Edep_energy").reset();
+ }
+
+ protected void highratioplots(EventHeader event, List<List<SimCalorimeterHit>> hitlists, double ratio, double fraction) {
+ double ratio_thresh = 1.25;
+ double frac_upper = 0.91;
+ double frac_lower = 0.85;
+ boolean higher = ratio > ratio_thresh;
+ boolean inbox = higher && fraction < frac_upper && fraction > frac_lower;
+ for (List<SimCalorimeterHit> simCalorimeterHits : hitlists) {
+ String colName = event.getMetaData(simCalorimeterHits).getName();
+ if (colName.contains("Edep")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ int pdgid = calorimeterHit.getMCParticle(i).getPDGID();
+ double en = calorimeterHit.getContributedEnergy(i);
+ if (inbox) {
+ aida.cloud1D("ratioplots/edep/higher/mcp pdgids").fill(pdgid);
+ aida.cloud1D("ratioplots/edep/higher/mcp pdgids by en").fill(pdgid, en);
+ } else {
+ aida.cloud1D("ratioplots/edep/lower/mcp pdgids").fill(pdgid);
+ aida.cloud1D("ratioplots/edep/lower/mcp pdgids by en").fill(pdgid, en);
+ }
+ }
+ }
+ } // end if Edep
+ if (colName.contains("Opti")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ int pdgid = calorimeterHit.getMCParticle(i).getPDGID();
+ double en = calorimeterHit.getContributedEnergy(i);
+ if (inbox) {
+ aida.cloud1D("ratioplots/opti/higher/mcp pdgids").fill(pdgid);
+ aida.cloud1D("ratioplots/opti/higher/mcp pdgids by en").fill(pdgid, en);
+ } else {
+ aida.cloud1D("ratioplots/opti/lower/mcp pdgids").fill(pdgid);
+ aida.cloud1D("ratioplots/opti/lower/mcp pdgids by en").fill(pdgid, en);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ protected void fitprofile(IProfile1D profile, String Functionname) {
+
+ System.out.println("DualCorrection: fitprofile");
+
+ int profbins = profile.axis().bins();
+ double profxmin = profile.axis().lowerEdge();
+ double profxmax = profile.axis().upperEdge();
+
+ jminuitResult = jminuit.fit(profile, poly);
+ System.out.println("jminuit Chi2=" + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ functionFactory.cloneFunction(Functionname, jminuitResult.fittedFunction());
+ System.out.println("correction function: " + result[0] + "," + result[1] + "," + result[2] + "," + result[3]);
+ try {
+ out.write(Particlename + " " + Ein_prev + " GeV\n");
+ out.write(result[0] + "," + result[1] + "," + result[2] + "," + result[3]);
+ } catch (IOException ex) {
+ Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ protected void convertandfit(ICloud1D[] slices, IHistogram1D[] conv_slices) {
+ System.out.println("DualCorrection: convert and fit:");
+ try {
+ out.write(Particlename + " " + Ein_prev + " GeV\n");
+ } catch (IOException ex) {
+ Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_ratio_" + Fitters[i];
+ point_ratio[i] = 0;
+ dps_ratio[i] = dpsFactory.create(dpsname, "ratio", 2);
+ }
+ for (int i = 0; i < slices.length; i++) {
+ if (slices[i].isConverted()) {
+ System.out.println("convert and fit already converted");
+ conv_slices[i] = slices[i].histogram();
+ } else {
+ double meanc = slices[i].mean();
+ double rmsc = slices[i].rms();
+ nsigmas = 3.;
+ nbins = 100;
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ slices[i].setConversionParameters(nbins, minx, maxx);
+ slices[i].convertToHistogram();
+ conv_slices[i] = slices[i].histogram();
+ }
+ int entries = conv_slices[i].entries();
+ if (entries > 300) {
+ for (int ii = 0; ii < Fitters.length; ii++) {
+ System.out.println("Fitter: " + Fitters[ii]);
+ gauss.setParameter("amplitude", conv_slices[i].maxBinHeight());
+ gauss.setParameter("mean", conv_slices[i].mean());
+ gauss.setParameter("sigma", conv_slices[i].rms());
+ jminuit = fitFactory.createFitter(Fitters[ii], "jminuit");
+ jminuitResult = jminuit.fit(conv_slices[i], gauss);
+ System.out.println("jminuit " + Fitters[ii] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ String functionname = "ratio fitted gauss slices: " + i + " " + Fitters[ii];
+ System.out.println(functionname);
+ functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+ try {
+ out.write(result[0] + "," + result[1] + "," + result[2] + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ IDataPoint dp_cerene = dps_ratio[ii].addPoint();
+ double ratio = (double) (i * binwidth) + 0.5 * binwidth;
+ dp_cerene.coordinate(0).setValue(ratio);
+ dp_cerene.coordinate(0).setErrorPlus(binwidth * 0.5);
+ dp_cerene.coordinate(0).setErrorMinus(binwidth * 0.5);
+ dp_cerene.coordinate(1).setValue(result[1]);
+ dp_cerene.coordinate(1).setErrorPlus(Math.abs(result[2]));
+ dp_cerene.coordinate(1).setErrorMinus(Math.abs(result[2]));
+ point_ratio[ii]++;
+ }
+ }
+ }
+ //Now fit the data point sets:
+ for (int ii = 0; ii < Fitters.length; ii++) {
+ System.out.println("Fitter: " + Fitters[ii]);
+ jminuitResult = jminuit.fit(dps_ratio[ii], poly);
+ String fname = "correction " + Fitters[ii] + " fitted poly";
+ System.out.println(fname);
+ functionFactory.cloneFunction(fname, jminuitResult.fittedFunction());
+ System.out.println(Fitters[ii] + ": " + jminuitResult.quality());
+ }
+ }
+
+ public void setMyAIDAFilename(String jpgfile) {
+ System.out.println("setMyVariable");
+ AIDAFile = jpgfile;
+ System.out.println(AIDAFile);
+ }
+
+ public void setMyFilename(String filename) {
+ System.out.println("setMyFilename");
+ file_name = filename;
+ System.out.println(file_name);
+ }
+
+ public void setMyPhysicsList(String list) {
+ System.out.println("setMyPhysicsList");
+ this.PhysicsList = list;
+ System.out.println(list);
+ }
+
+ public void setMyDetector(String Detector) {
+ this.Detector = Detector;
+ }
+
+ public void setMyMaterial(String Material) {
+ this.Material = Material;
+ }
+
+ public void setMyDensity(Double Density) {
+ this.Density = Density;
+ }
+
+ public void setMyrindex(Double rindex) {
+ this.rindex = rindex;
+ }
+
+ public void setMyCollectionName(String CollectionName) {
+ this.CollectionName = CollectionName;
+ }
+
+ public void setMyCerenkovThres(Double CerenkovThres) {
+ this.CerenkovThres = CerenkovThres;
+ }
+
+ public void setMyIonizationThres(Double IonizationThres) {
+ this.IonizationThres = IonizationThres;
+ }
+
+ double roundTwoDecimals(double d) {
+ DecimalFormat twoDForm = new DecimalFormat("#.##");
+ return Double.valueOf(twoDForm.format(d));
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N DRFunctionFactory.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DRFunctionFactory.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,199 @@
+/*
+ * this class returns the electron, erenkov and dualreadout correction for
+ * a specific detector configuration
+ */
+package org.lcsim.mcd.analysis.DRCorrection;
+
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.Iterator;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author wenzel
+ */
+final public class DRFunctionFactory {
+
+ private AIDA aida;
+ private IFunctionFactory functionFactory;
+ static DRFunctionFactory mof;
+ private static Map<DetectorConfiguration, ArrayList<IFunction>> mp;
+
+ /**
+ use very simple singleton pattern to make
+ sure that only one instance of this class.
+ */
+ public static DRFunctionFactory getInstance() {
+ if (mof == null) {
+ mof = new DRFunctionFactory();
+ }
+ return mof;
+ }
+
+ private DRFunctionFactory() {
+ aida = AIDA.defaultInstance();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ mp = new HashMap<DetectorConfiguration, ArrayList<IFunction>>();
+ ArrayList<IFunction> al = new ArrayList<IFunction>();
+ al = new ArrayList<IFunction>();
+ DetectorConfiguration dualcor = new DetectorConfiguration("mcdrcal01", "PbW04", 8.28, 2.21, "NoCuts", 0.02, 0.02, "FTFP_BERT");
+ IFunction Corfu = functionFactory.createFunctionByName("mcdrcal01_NoCuts_FTFP_BERT", "p4");
+ Corfu.setParameter("p0", 0.794055);
+ Corfu.setParameter("p1", 0.539570);
+ Corfu.setParameter("p2", -0.68976);
+ Corfu.setParameter("p3", 0.383136);
+ Corfu.setParameter("p4", 0.0);
+ al.add(Corfu);
+ // adding or set elements in Map by put method key and value pair
+ IFunction efu = functionFactory.createFunctionByName("ec_mcdrcal01_NoCuts_FTFPP_BERT", "p1");
+ efu.setParameter("p0", 4.026e-3);
+ efu.setParameter("p1", 0.9998);
+ al.add(efu);
+ IFunction Cerfu = functionFactory.createFunctionByName("cc_mcdrcal01_NoCuts_FTFP_BERT", "p1");
+ Cerfu.setParameter("p0", 3.6605e-3);
+ Cerfu.setParameter("p1", 1.87977e-5);
+ al.add(Cerfu);
+ mp.put(dualcor, al);
+ dualcor = new DetectorConfiguration("mcdrcal01", "PbW04", 8.28, 2.21, "10ns", 0.02, 0.02, "FTFP_BERT");
+ Corfu = functionFactory.createFunctionByName("mcdrcal01_10ns_FTFP_BERT", "p4");
+ Corfu.setParameter("p0", 0.435559);
+ Corfu.setParameter("p1", 0.552895);
+ Corfu.setParameter("p2", -0.48328);
+ Corfu.setParameter("p3", 0.518759);
+ Corfu.setParameter("p4", 0.0);
+ al.add(Corfu);
+ // adding or set elements in Map by put method key and value pair
+ efu = functionFactory.createFunctionByName("ec_mcdrcal01_10ns_FTFPP_BERT", "p1");
+ efu.setParameter("p0", 1.5401e-2);
+ efu.setParameter("p1", 1.0008);
+ al.add(efu);
+ Cerfu = functionFactory.createFunctionByName("cc_mcdrcal01_10ns_FTFP_BERT", "p1");
+ Cerfu.setParameter("p0", 0.00734336);
+ Cerfu.setParameter("p1", 1.881998E-5);
+ al.add(Cerfu);
+ mp.put(dualcor, al);
+ dualcor = new DetectorConfiguration("mcdrcal01", "PbW04", 8.28, 2.21, "proton-NoCuts", 0.02, 0.02, "FTFP_BERT");
+ Corfu = functionFactory.createFunctionByName("mcdrcal01_protNoCuts_FTFP_BERT", "p4");
+ Corfu.setParameter("p0", 0.794055);
+ Corfu.setParameter("p1", 0.539570);
+ Corfu.setParameter("p2", -0.68976);
+ Corfu.setParameter("p3", 0.383136);
+ Corfu.setParameter("p4", 0.0);
+ al.add(Corfu);
+ // adding or set elements in Map by put method key and value pair
+ efu = functionFactory.createFunctionByName("ec_mcdrcal01_protNoCuts_FTFPP_BERT", "p1");
+ efu.setParameter("p0", 4.026e-3);
+ efu.setParameter("p1", 0.9998);
+ al.add(efu);
+ Cerfu = functionFactory.createFunctionByName("cc_mcdrcal01_protNoCuts_FTFP_BERT", "p1");
+ Cerfu.setParameter("p0", 3.6605e-3);
+ Cerfu.setParameter("p1", 1.87977e-5);
+ al.add(Cerfu);
+ mp.put(dualcor, al);
+ dualcor = new DetectorConfiguration("mcdrcal01", "PbW04", 8.28, 2.21, "LowEnNoCuts", 0.02, 0.02, "FTFP_BERT");
+ Corfu = functionFactory.createFunctionByName("mcdrcal01_LowEnNoCuts_FTFP_BERT", "p4");
+ Corfu.setParameter("p0", 0.794055);
+ Corfu.setParameter("p1", 0.539570);
+ Corfu.setParameter("p2", -0.68976);
+ Corfu.setParameter("p3", 0.383136);
+ Corfu.setParameter("p4", 0.0);
+ al.add(Corfu);
+ // adding or set elements in Map by put method key and value pair
+ efu = functionFactory.createFunctionByName("ec_mcdrcal01_LowEnNoCuts_FTFPP_BERT", "p1");
+ efu.setParameter("p0", 4.026e-3);
+ efu.setParameter("p1", 0.9998);
+ al.add(efu);
+ Cerfu = functionFactory.createFunctionByName("cc_mcdrcal01_LowEnNoCuts_FTFP_BERT", "p1");
+ Cerfu.setParameter("p0", 3.6605e-3);
+ Cerfu.setParameter("p1", 1.87977e-5);
+ al.add(Cerfu);
+ mp.put(dualcor, al);
+ }
+
+ public void print() {
+ System.out.println("Currently " + mp.size() + " Detector Configurations are available");
+ //Get Map in Set interface to get key and value
+ Set s = mp.entrySet();
+ //Move next key and value of Map by iterator
+ Iterator it = s.iterator();
+ while (it.hasNext()) {
+ // key=value separator this by Map.Entry to get key and value
+ Map.Entry m = (Map.Entry) it.next();
+ // getKey is used to get key of Map
+ DetectorConfiguration key = (DetectorConfiguration) m.getKey();
+ key.print();
+ // getValue is used to get value of key in Map
+ ArrayList<IFunction> map = (ArrayList<IFunction>) m.getValue();
+ IFunction value = map.get(0);
+ System.out.println("Nr of Parameters: " + value.numberOfParameters());
+ System.out.println("Name of Function: " + value.title());
+ double[] params = value.parameters();
+ for (int i = 0; i < params.length; i++) {
+ System.out.println(params[i]);
+ }
+ value = map.get(1);
+ System.out.println("Nr of Parameters: " + value.numberOfParameters());
+ System.out.println("Name of Function: " + value.title());
+ double[] params1 = value.parameters();
+ for (int i = 0; i < params1.length; i++) {
+ System.out.println(params1[i]);
+ }
+ value = map.get(2);
+ System.out.println("Nr of Parameters: " + value.numberOfParameters());
+ System.out.println("Name of Function: " + value.title());
+ double[] params2 = value.parameters();
+ for (int i = 0; i < params2.length; i++) {
+ System.out.println(params2[i]);
+ }
+ }
+ }
+// get dual readout correction function
+
+ public IFunction getdcFunction(DetectorConfiguration dc) {
+ IFunction refunc = null;
+ if (mp.containsKey(dc)) {
+ ArrayList<IFunction> map = (ArrayList<IFunction>) mp.get(dc);
+ refunc = (IFunction) map.get(0);
+ }
+ int n = refunc.parameters().length;
+ System.out.println("dcfunc");
+ for (int j = 0; j < n; j++) {
+ System.out.println(refunc.parameters()[j]);
+ }
+ return refunc;
+ }
+// get cerenkov scale
+
+ public IFunction getccFunction(DetectorConfiguration dc) {
+ IFunction refunc = null;
+ if (mp.containsKey(dc)) {
+ ArrayList<IFunction> map = (ArrayList<IFunction>) mp.get(dc);
+ refunc = (IFunction) map.get(2);
+ }
+ return refunc;
+ }
+// get electron scale:
+
+ public IFunction getecFunction(DetectorConfiguration dc) {
+ IFunction refunc = null;
+ if (mp.containsKey(dc)) {
+ ArrayList<IFunction> map = (ArrayList<IFunction>) mp.get(dc);
+ refunc = (IFunction) map.get(1);
+ }
+ return refunc;
+ }
+/// check if there is an entry for a given DetectorConfiguration
+
+ boolean checkdc(DetectorConfiguration dc) {
+ if (mp.containsKey(dc)) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N electron_calib_steering.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ electron_calib_steering.xml 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,29 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+ </classpath>
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-10gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-25gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-50gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-100gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+ <execute>
+ <driver name="ElectronCalibrationDriver"/>
+ </execute>
+ <drivers>
+ <driver name="ElectronCalibrationDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.ElectronCalibrationDriver">
+ <myAIDAFilename>elec_calibration.aida</myAIDAFilename>
+ <myFilename>elec_calibration.txt</myFilename>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N DRResolutionDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DRResolutionDriver.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,114 @@
+/*
+ *
+ * This Driver is used to run the module that checks that all the calibration using
+ * electrons and pions makes sense (e.g. peaks show up at the right spot).
+ * In addition the module estimate the energy resolution of the dual readout calorimeter
+ * for single particles.
+ */
+package org.lcsim.mcd.analysis.DRCorrection;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author wenzel
+ */
+public class DRResolutionDriver extends Driver {
+ //
+ // default detector configuration:
+ //
+
+ private String Detector = "mcdrcal01";
+ private String Material = "PbW04";
+ private Double Density = 8.28;
+ private Double rindex = 2.21;
+ private String CollectionName = "10ns";
+ private Double CerenkovThres = 0.02;
+ private Double IonizationThres = 0.02;
+ private String PhysicsList = "FTFP_BERT";
+ String AIDAFile = null;
+ String file_name = null;
+ Resolution resol;
+
+ public DRResolutionDriver() {
+
+ resol = new Resolution();
+ add(resol);
+ }
+
+ @Override
+ public void startOfData() {
+ System.out.println("DRCalibrationDriver:startOfData");
+ System.out.println(AIDAFile);
+ if (AIDAFile != null) {
+ resol.setMyAIDAFilename(AIDAFile);
+ } else {
+ System.err.println("DRCalibrationDriver: AIDAFile variable must be set");
+ System.exit(1); // exit if variable is not set
+ }
+ if (file_name != null) {
+ resol.setMyFilename(file_name);
+ } else {
+ System.err.println("DRCalibrationDriver: file_name variable must be set");
+ System.exit(1); // exit if variable is not set
+ }
+ resol.setMyPhysicsList(PhysicsList);
+ resol.setMyDetector(Detector);
+ resol.setMyMaterial(Material);
+ resol.setMyDensity(Density);
+ resol.setMyrindex(rindex);
+ resol.setMyCollectionName(CollectionName);
+ resol.setMyCerenkovThres(CerenkovThres);
+ resol.setMyIonizationThres(IonizationThres);
+ resol.startOfData();
+
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("DRCalibrationDriver:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+
+ public void setMyFilename(String file_name) {
+ System.out.println("DRCalibrationDriver:setMyFilename");
+ this.file_name = file_name;
+ System.out.println(file_name);
+ }
+
+ public void setMyPhysicsList(String list) {
+ System.out.println("setMyPhysicsList");
+ this.PhysicsList = list;
+ System.out.println(list);
+ }
+
+ public void setMyDetector(String Detector) {
+ this.Detector = Detector;
+ }
+
+ public void setMyMaterial(String Material) {
+ this.Material = Material;
+ }
+
+ public void setMyDensity(double Density) {
+ System.out.println("DRCalibrationDriver:setMyDensity");
+ this.Density = Density;
+ System.out.println(Density);
+ }
+
+ public void setMyrindex(double rindex) {
+ this.rindex = rindex;
+ }
+
+ public void setMyCollectionName(String CollectionName) {
+ this.CollectionName = CollectionName;
+ }
+
+ public void setMyCerenkovThres(double CerenkovThres) {
+ this.CerenkovThres = CerenkovThres;
+ }
+
+ public void setMyIonizationThres(double IonizationThres) {
+ this.IonizationThres = IonizationThres;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N ElectronCorrection.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ElectronCorrection.java 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,380 @@
+package org.lcsim.mcd.analysis.DRCorrection;
+
+//package org.lcsim.mcd.analysis;
+
+/**
+ *
+ * @author wenzel
+ * This Driver is used to calibrate the correction of a dual readout calorimeter to electrons.
+ * (Determining the energy scale of the calorimeter response.)
+ */
+import org.lcsim.mcd.analysis.*;
+import hep.aida.ICloud1D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import java.io.BufferedWriter;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.aida.AIDA;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+
+public class ElectronCorrection extends Driver {
+
+ private AIDA aida;
+ String AIDAFile;
+ String file_name;
+ FileWriter fstream;
+ BufferedWriter out;
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsFactory;
+ IFunction gauss;
+ boolean first;
+ boolean firstEvent;
+ double E_in;
+ Integer Ein;
+ Integer Ein_prev;
+ double E_kin;
+ String Particlename;
+ ICloud1D Edep;
+ ICloud1D Eceren;
+ double nsigmas;
+ int nbins;
+ IDataPointSet dps_emean = null;
+ IDataPointSet dps_cerenemean = null;
+ IDataPointSet[] dps_cerene = null;
+ int point_cerene[];
+ int point;
+ int point_cerenemean;
+ double[] result;
+ double errors[];
+ String[] Fitters = {"Chi2", "leastsquares", "bml", "cleverchi2"};
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ double t_thresh;
+ double e_thresh;
+
+ public ElectronCorrection() {
+ // if e_thresh > 0 then t_thresh MUST be > 0 as well.
+ t_thresh = 10.0;
+ e_thresh = -1.0;
+ file_name = "ElectronCorrection.txt";
+ AIDAFile = "ElectronCorrection.aida";
+ aida = AIDA.defaultInstance();
+ dps_cerene = new IDataPointSet[Fitters.length];
+ point_cerene = new int[Fitters.length];
+ fitFactory = aida.analysisFactory().createFitFactory();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_cerene_" + Fitters[i];
+ point_cerene[i] = 0;
+ dps_cerene[i] = dpsFactory.create(dpsname, "electron response mean", 2);
+ }
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+ // Create a two dimensional IDataPointSet.
+ dps_emean = dpsFactory.create("dps_emean", "electron response mean", 2);
+ dps_cerenemean = dpsFactory.create("dps_cerenemean", "electron response mean", 2);
+ point = 0;
+ point_cerenemean = 0;
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ }
+
+ @Override
+ protected void process(EventHeader event) {
+ E_in = 0.0;
+ E_kin = 0.0;
+ String DirName = null;
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ if (particle.getGeneratorStatus() == 0) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+ if (particle.getProductionTime() == 0.0) {
+ Particlename = particle.getType().getName();
+ }
+ }
+ break;
+ }
+ Ein = (int) Math.floor(E_kin + 0.5d);
+ String E_str = Ein.toString();
+ DirName = Particlename.concat(E_str);
+ if (Ein != Ein_prev) {
+ if (first) {
+ first = false;
+ } else {
+ System.out.println("ElectronCorrection:First Event");
+ convertandfit();
+ }
+ Ein_prev = Ein;
+ System.out.println("First Event:");
+ System.out.println("E_in: " + E_in);
+ System.out.println("E_kin: " + E_kin);
+ System.out.println("Name: " + Particlename);
+ aida.tree().cd("/");
+ aida.tree().mkdir(DirName);
+ aida.tree().cd(DirName);
+ Edep = aida.histogramFactory().createCloud1D("Edep", "Energy Cloud", 100000, "autoConvert = false");
+ Eceren = aida.histogramFactory().createCloud1D("Eceren", "Cerenkov Cloud", 100000, "autoConvert = false");
+ System.out.println("DirName: " + DirName);
+ }
+ List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+ double sumEEdep = 0.0;
+ double sumECeren = 0.0;
+ for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+ String CollectionName = event.getMetaData(simCalorimeterHits).getName();
+ if (CollectionName.contains("Edep")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (t_thresh > 0.0 || e_thresh > 0.0) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - calorimeterHit.getPositionVec().magnitude()/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < t_thresh) && (edep > e_thresh) ) {
+ sumEEdep += edep;
+ }
+ }
+ }
+ else {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ }
+ } // end if Edep
+ if (CollectionName.contains("Opti")) {
+// if (CollectionName.startsWith("Edep_") && CollectionName.endsWith("DigiHits")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (t_thresh > 0.0 || e_thresh > 0.0) {
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - calorimeterHit.getPositionVec().magnitude()/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < t_thresh) && (edep > e_thresh) ) {
+ sumECeren += edep;
+ }
+ }
+ }
+ else {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ }
+ }
+ }
+ Edep.fill(sumEEdep);
+ Eceren.fill(sumECeren);
+ }
+
+ @Override
+ protected void endOfData() {
+ System.out.println("ElectronCorrection:End of Data:");
+ convertandfit();
+ fstream = null;
+ try {
+ fstream = new FileWriter(file_name);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ out = new BufferedWriter(fstream);
+ aida.tree().cd("/");
+ IFunction line = functionFactory.createFunctionByName("line", "p1");
+ IFitter linefit = fitFactory.createFitter("Chi2", "jminuit");
+ IFitResult resultline = linefit.fit(dps_emean, line);
+ functionFactory.cloneFunction("e mean fitted line ", resultline.fittedFunction());
+ double[] fresult = resultline.fittedParameters();
+ try {
+ out.write("Ionization scale results:\n");
+ out.write("Chi2 = " + resultline.quality() + "\n");
+ out.write(fresult[0] + " , " + fresult[1] + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ //
+ // now deal with the cerenkov stuff:
+ //
+ try {
+ out.write("Cerenkov scale results:\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ for (int i = 0; i < Fitters.length; i++) {
+ System.out.println("Fitter: " + Fitters[i]);
+ resultline = linefit.fit(dps_cerene[i], line);
+ String fname = "cerenkov " + Fitters[i] + " fitted line";
+ System.out.println(fname);
+ functionFactory.cloneFunction(fname, resultline.fittedFunction());
+ System.out.println(Fitters[i] + ": " + resultline.quality());
+ fresult = resultline.fittedParameters();
+ try {
+ out.write("Fitter: " + Fitters[i]+"\n");
+ out.write("Chi2 = " + resultline.quality() + "\n");
+ out.write(fresult[0] + " , " + fresult[1] + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ try {
+ out.close();
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ @Override
+ protected void resume() {
+ System.out.println("ElectronCorrection:resume");
+ firstEvent =true;
+ aida.cloud1D("c_Edep_energy").reset();
+ }
+
+ protected void convertandfit() {
+ System.out.println("ElectronCorrection:convertandfit");
+ IHistogram1D Edep_conv;
+ if (Edep.isConverted()) {
+ Edep_conv = Edep.histogram();
+ } else {
+ System.out.println("Converting EDep");
+ double meanc = Edep.mean();
+ double rmsc = Edep.rms();
+ nsigmas = 3.;
+ nbins = 100;
+
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep.setConversionParameters(nbins, minx, maxx);
+ Edep.convertToHistogram();
+ Edep_conv = Edep.histogram();
+ }
+
+ gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_conv.mean());
+ gauss.setParameter("sigma", Edep_conv.rms());
+
+ dps_emean.addPoint();
+ IDataPoint dp = dps_emean.point(point);
+ dp.coordinate(1).setValue(Ein_prev);
+ dp.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+ dp.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+ dp.coordinate(0).setValue(Edep_conv.mean());
+ dp.coordinate(0).setErrorPlus(Edep_conv.rms());
+ dp.coordinate(0).setErrorMinus(Edep_conv.rms());
+ point++;
+// chi 2 fit:
+ System.out.println("chi2 fit:");
+ /*jminuit = fitFactory.createFitter("Chi2", "jminuit");
+ jminuitResult = jminuit.fit(Edep_conv, gauss);
+ System.out.println(Edep_conv.maxBinHeight() + " , " + Edep_conv.mean() + " , " + Edep_conv.rms());
+ System.out.println("jminuit Chi2=" + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ functionFactory.cloneFunction("fitted gauss (jminuitchi2)", jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+
+ // least squares fit:
+ System.out.println("least squares fit:");
+ gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_conv.mean());
+ gauss.setParameter("sigma", Edep_conv.rms());
+ IFitter jminuitls = fitFactory.createFitter("leastsquares", "jminuit");
+ IFitResult jminuitResultls = jminuitls.fit(Edep_conv, gauss);
+ System.out.println("jminuit ls=" + jminuitResultls.quality());
+ double[] resultls = jminuitResultls.fittedParameters();
+ functionFactory.cloneFunction("fitted gauss (jminuitls)", jminuitResultls.fittedFunction());
+ System.out.println(resultls[0] + "," + resultls[1] + "," + resultls[2]);
+*/
+ //
+ // now deal with cerenkov response:
+ //
+ IHistogram1D Eceren_conv;
+
+ if (Eceren.isConverted()) {
+ Eceren_conv = Eceren.histogram();
+ } else {
+ System.out.println("Converting Eceren");
+ double meanc = Eceren.mean();
+ double rmsc = Eceren.rms();
+ nsigmas = 5.;
+ nbins = 100;
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Eceren.setConversionParameters(nbins, minx, maxx);
+ Eceren.convertToHistogram();
+ Eceren_conv = Eceren.histogram();
+ }
+
+ gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+ gauss.setParameter("mean", Eceren_conv.mean());
+ gauss.setParameter("sigma", Eceren_conv.rms());
+ dps_cerenemean.addPoint();
+ IDataPoint dp_cerenemean = dps_cerenemean.point(point_cerenemean);
+ dp_cerenemean.coordinate(1).setValue(Ein_prev);
+ dp_cerenemean.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+ dp_cerenemean.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+ dp_cerenemean.coordinate(0).setValue(Eceren_conv.mean());
+ dp_cerenemean.coordinate(0).setErrorPlus(Eceren_conv.rms());
+ dp_cerenemean.coordinate(0).setErrorMinus(Eceren_conv.rms());
+ point_cerenemean++;
+
+ for (int i = 0; i < Fitters.length; i++) {
+ System.out.println("Fitter: " + Fitters[i]);
+ gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+ gauss.setParameter("mean", Eceren_conv.mean());
+ gauss.setParameter("sigma", Eceren_conv.rms());
+ jminuit = fitFactory.createFitter(Fitters[i], "jminuit");
+ jminuitResult = jminuit.fit(Eceren_conv, gauss);
+ System.out.println("jminuit " + Fitters[i] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ String functionname = "ceren fitted gauss " + Fitters[i];
+ System.out.println(functionname);
+ functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+ IDataPoint dp_cerene = dps_cerene[i].addPoint();
+ dp_cerene.coordinate(0).setValue(result[1]);
+ dp_cerene.coordinate(0).setErrorPlus(Math.abs(result[2]));
+ dp_cerene.coordinate(0).setErrorMinus(Math.abs(result[2]));
+ dp_cerene.coordinate(1).setValue(Ein_prev);
+ dp_cerene.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+ dp_cerene.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+ point_cerene[i]++;
+ }
+
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("ElectronCalibrationDriver:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+
+ public void setMyFilename(String file_name) {
+ System.out.println("ElectronCalibrationDriver:setMyFilename");
+ this.file_name = file_name;
+ System.out.println(file_name);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N pion_calib_steering.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ pion_calib_steering.xml 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,38 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+ </classpath>
+
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-10gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-25gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-50gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+
+ <execute>
+ <driver name="DRCalibrationDriver"/>
+ </execute>
+ <drivers>
+ <driver name="DRCalibrationDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.DRCalibrationDriver">
+ <myAIDAFilename>pipl_calibration-NoCuts.aida</myAIDAFilename>
+ <myFilename>pipl_calibration-NoCuts.txt</myFilename>
+ <myPhysicsList>FTFP_BERT</myPhysicsList>
+ <myDetector>mcdrcal01</myDetector>
+ <myMaterial>PbW04</myMaterial>
+ <myDensity>8.28</myDensity>
+ <myrindex>2.21</myrindex>
+ <myCollectionName>NoCuts</myCollectionName>
+ <myCerenkovThres>0.02</myCerenkovThres>
+ <myIonizationThres>0.02</myIonizationThres>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N resolution_pions.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ resolution_pions.xml 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,35 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+</classpath>
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-10gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-25gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/pipl/pipl-90deg-50gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+ <execute>
+ <driver name="DRResolutionDriver"/>
+ </execute>
+ <drivers>
+ <driver name="DRResolutionDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.DRResolutionDriver">
+ <myAIDAFilename>resolution_pipl-NoCuts.aida</myAIDAFilename>
+ <myFilename>resolution_pipl-NoCuts.txt</myFilename>
+ <myPhysicsList>FTFP_BERT</myPhysicsList>
+ <myDetector>mcdrcal01</myDetector>
+ <myMaterial>PbW04</myMaterial>
+ <myDensity>8.28</myDensity>
+ <myrindex>2.21</myrindex>
+ <myCollectionName>NoCuts</myCollectionName>
+ <myCerenkovThres>0.02</myCerenkovThres>
+ <myIonizationThres>0.02</myIonizationThres>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -N proton_calib_steering.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ proton_calib_steering.xml 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,38 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+ </classpath>
+
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-10gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-25gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-50gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+
+ <execute>
+ <driver name="DRCalibrationDriver"/>
+ </execute>
+ <drivers>
+ <driver name="DRCalibrationDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.DRCalibrationDriver">
+ <myAIDAFilename>proton_calibration-NoCuts.aida</myAIDAFilename>
+ <myFilename>proton_calibration-NoCuts.txt</myFilename>
+ <myPhysicsList>FTFP_BERT</myPhysicsList>
+ <myDetector>mcdrcal01</myDetector>
+ <myMaterial>PbW04</myMaterial>
+ <myDensity>8.28</myDensity>
+ <myrindex>2.21</myrindex>
+ <myCollectionName>proton-NoCuts</myCollectionName>
+ <myCerenkovThres>0.02</myCerenkovThres>
+ <myIonizationThres>0.02</myIonizationThres>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis
diff -N nbactions.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ nbactions.xml 20 Mar 2014 20:24:47 -0000 1.1
@@ -0,0 +1,22 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<actions>
+ <action>
+ <actionName>build</actionName>
+ <goals>
+ <goal>install</goal>
+ </goals>
+ <properties>
+ <skipTests>true</skipTests>
+ </properties>
+ </action>
+ <action>
+ <actionName>rebuild</actionName>
+ <goals>
+ <goal>clean</goal>
+ <goal>install</goal>
+ </goals>
+ <properties>
+ <skipTests>true</skipTests>
+ </properties>
+ </action>
+ </actions>
mcd-analysis
diff -u -r1.10 -r1.11
--- pom.xml 23 Sep 2013 23:51:12 -0000 1.10
+++ pom.xml 20 Mar 2014 20:24:47 -0000 1.11
@@ -77,5 +77,10 @@
<artifactId>lcsim-contrib</artifactId>
<version>2.0-SNAPSHOT</version>
</dependency>
+ <dependency>
+ <groupId>${project.groupId}</groupId>
+ <artifactId>lcsim-contrib</artifactId>
+ <version>2.0-SNAPSHOT</version>
+ </dependency>
</dependencies>
</project>
CVSspam 0.2.12