Modified Hans' DRCorrection stuff for use with mcd. Also some random experimental analysis drivers.
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()); + } + + } +}
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())); + } + } +}
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); + } + } + } +}
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); + } + + +}
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; + } +}
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); + } +}
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
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 { + +}
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) { + } + } +}
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"); +// } +}
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; + } +}
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; + } +}
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; + } +}
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); + } +}
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() { + } +}
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) {
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;
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;
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; } }
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()); - } - -}
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; + } +}
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; + } +}
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>
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; + } +} +
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); + } +}
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)); + } +}
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; + } + } +}
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>
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; + } +}
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); + } +}
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>
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>
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>
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>
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>
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1