lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N JetDriverExtended.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ JetDriverExtended.java 8 Jun 2009 18:29:05 -0000 1.1
@@ -0,0 +1,325 @@
+package org.lcsim.contrib.Mbussonn;
+
+import hep.physics.jet.JadeEJetFinder;
+import hep.physics.jet.JetFinder;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.HepLorentzVector;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.LinkedList;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.digisim.MyLCRelation;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Vertex;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * A simple driver which can be used to find jets from ReconstructedParticles.
+ * The resuslting jets are stored in a new collection of ReconstructedParticles.
+ * @author tonyj
+ * @version $Id: JetDriverExtended.java,v 1.1 2009/06/08 18:29:05 mbussonn Exp $
+ */
+public class JetDriverExtended extends Driver
+{
+ private static final Hep3Vector IP = new BasicHep3Vector(0,0,0);
+ private static final String defaultOutputCollectionName = "Jets";
+ private String inputCollectionName;
+ private String outputCollectionName = defaultOutputCollectionName;
+ private JetFinder finder = defaultJetFinder();
+ List<LCRelation> rc2mc = new ArrayList<LCRelation>();
+ List<ReconstructedParticle> collectionOfRCP;
+
+ /** Creates a new instance of JetFinder with the default properties */
+ public JetDriverExtended()
+ {
+ }
+
+
+ public String getInputCollectionName()
+ {
+ return inputCollectionName;
+ }
+
+ /**
+ * The name of the input collection to use. If not set (or set to <code>null</code> uses
+ * the first collection of ReconstructedParticles found in the event.
+ * @param inputCollectionName
+ */
+
+ public void setInputCollectionName(String inputCollectionName)
+ {
+ this.inputCollectionName = inputCollectionName;
+ }
+
+ public String getOutputCollectionName()
+ {
+ return outputCollectionName;
+ }
+
+ /**
+ * The name of the output collection added to the event. If not set, or set to <code>null</code>,
+ * default to "Jets"
+ * @param outputCollectionName
+ */
+
+ public void setOutputCollectionName(String outputCollectionName)
+ {
+ this.outputCollectionName = outputCollectionName == null ? defaultOutputCollectionName : outputCollectionName;
+ }
+
+ public JetFinder getFinder()
+ {
+ return finder;
+ }
+
+ /**
+ * Set the jet finding algorithm to use
+ * @param finder
+ */
+ public void setFinder(JetFinder finder)
+ {
+ this.finder = finder == null ? defaultJetFinder() : finder;
+ }
+
+ protected JetFinder defaultJetFinder()
+ {
+ return new JadeEJetFinder(0.045);
+ }
+
+ @Override
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+ boolean hist = getHistogramLevel() > 0;
+ double totalEnergy=0;
+ // Find the input reconstructed Particles
+ List<ReconstructedParticle> input = null;
+ collectionOfRCP = new ArrayList<ReconstructedParticle>();
+ if (inputCollectionName ==null)
+ {
+ List<List<ReconstructedParticle>> listOfLists = event.get(ReconstructedParticle.class);
+ if (listOfLists.isEmpty()) {
+ input = new LinkedList<ReconstructedParticle>();
+ //List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+
+ for(MCParticle mcp : event.getMCParticles()){
+ if(mcp.getGeneratorStatus()==MCParticle.FINAL_STATE){
+ MCReconstructedParticle rcp = new MCReconstructedParticle(mcp);
+ input.add(rcp);
+ MyLCRelation lcr = new MyLCRelation(rcp, mcp);
+ rc2mc.add(lcr);
+ System.out.println("lcr to string mcp :"+lcr.getTo().toString()+" rcp "+lcr.getFrom().getClass());
+ collectionOfRCP.add(rcp);
+ System.out.println("adding rcp -> "+collectionOfRCP.size());
+ totalEnergy += mcp.getEnergy();
+
+ }
+ }
+ }
+ else{
+ input = listOfLists.get(0);
+ }
+ }
+ else
+ {
+ input = event.get(ReconstructedParticle.class,inputCollectionName);
+ }
+ // Build a list of 4-vectors from the reconstructed particles
+ Map<HepLorentzVector, ReconstructedParticle> map = new HashMap<HepLorentzVector, ReconstructedParticle>();
+ for (ReconstructedParticle p : input)
+ {
+ map.put(p.asFourVector(),p);
+ }
+
+ // Pass the data to the Jet finder
+ finder.setEvent(map.keySet());
+
+
+
+ //if (hist)
+ {
+ AIDA aida = AIDA.defaultInstance();
+ aida.cloud1D("sumE").fill(totalEnergy);
+ double tje =0;
+ for(int i=0;i<finder.njets();i++){
+ HepLorentzVector hlv = finder.jet(i);
+ double px = hlv.v3().x();
+ double py = hlv.v3().y();
+ double pz = hlv.v3().z();
+ double pr = Math.sqrt(px*px+py*py);
+ double pt = hlv.v3().magnitude();
+ double cosTheta = pz/pt;
+ double theta = Math.acos(cosTheta);
+ double eta = Math.log(Math.tan(theta/2.));
+ double phi = Math.atan2(px,py);
+ tje +=hlv.v3().magnitude() ;
+ if(1==finder.nParticlesPerJet(i)){
+ aida.cloud1D("One particle/jet eta ").fill(eta);
+ aida.cloud1D("One particle/jet phi ").fill(phi);
+ aida.cloud1D("One particle/jet pt").fill(pt);
+ aida.cloud1D("One particle/jet |H4|").fill(finder.jet(i).magnitude());
+ aida.cloud1D("One particle/jet pr").fill(pr);
+ aida.cloud1D("One particle/jet pz").fill(pz);
+ aida.cloud2D("One particle/jet direction").fill(pz,pr);
+
+ }else{
+
+ aida.cloud1D("Jet particle/jet eta ").fill(eta);
+ aida.cloud1D("Jet particle/jet phi ").fill(phi);
+ aida.cloud1D("Jet particle/jet pt").fill(pt);
+ aida.cloud1D("Jet particle/jet |H4|").fill(finder.jet(i).magnitude());
+ aida.cloud1D("Jet particle/jet pr").fill(pr);
+ aida.cloud1D("Jet particle/jet pz").fill(pz);
+ aida.cloud2D("Jet particle/jet direction").fill(pz,pr);
+
+ aida.cloud2D("jet direction ").fill(hlv.v3().x(),hlv.v3().y());
+ aida.cloud1D("JetDriver/particlesPerJet").fill(finder.nParticlesPerJet(i));
+ aida.cloud1D("JetDriver/nJets").fill(finder.njets());
+ }
+
+ }
+ aida.cloud1D("total jet energy").fill(tje);
+
+
+ }
+
+ // Loop over the output jets and create a new ReconstructedParticle for each one, pointing back to
+ // the original particles
+
+ List<ReconstructedParticle> output = new ArrayList<ReconstructedParticle>();
+ int nJets = finder.njets();
+ for (int i=0; i<nJets; i++)
+ {
+ Jet jet = new Jet(finder.jet(i));
+ for (HepLorentzVector pj : (List<HepLorentzVector>) finder.particlesInJet(i))
+ {
+ jet.addParticle(map.get(pj));
+ }
+ output.add(jet);
+ }
+
+ // Add the list of jets to the event
+
+ event.put("collectionOfRCP",collectionOfRCP,ReconstructedParticle.class,0);
+ event.put(outputCollectionName, output, ReconstructedParticle.class, 0);
+ event.put("rc2mc", rc2mc, LCRelation.class, 0);
+ }
+ private class Jet implements ReconstructedParticle
+ {
+ private HepLorentzVector fourVector;
+ private List<ReconstructedParticle> particles = new ArrayList<ReconstructedParticle>();
+ private double charge;
+
+ Jet(HepLorentzVector fourVector)
+ {
+ this.fourVector = fourVector;
+ }
+
+ public void addTrack(org.lcsim.event.Track track)
+ {
+ throw new UnsupportedOperationException("Add track to jet");
+ }
+
+ public void addParticleID(org.lcsim.event.ParticleID pid)
+ {
+ throw new UnsupportedOperationException("Add track to jet");
+ }
+
+ public void addParticle(ReconstructedParticle particle)
+ {
+ particles.add(particle);
+ charge += particle.getCharge();
+ }
+
+ public int getType()
+ {
+ // Fixme: what should we return for Jet?
+ return 0;
+ }
+
+ public List<org.lcsim.event.Track> getTracks()
+ {
+ return Collections.EMPTY_LIST;
+ }
+
+ public hep.physics.vec.Hep3Vector getReferencePoint()
+ {
+ return IP;
+ }
+
+ public List<ReconstructedParticle> getParticles()
+ {
+ return particles;
+ }
+
+ public List<org.lcsim.event.ParticleID> getParticleIDs()
+ {
+ return Collections.EMPTY_LIST;
+ }
+
+ public void addCluster(org.lcsim.event.Cluster cluster)
+ {
+ throw new UnsupportedOperationException("Add track to jet");
+ }
+
+ public HepLorentzVector asFourVector()
+ {
+ return fourVector;
+ }
+
+ public double getCharge()
+ {
+ return charge;
+ }
+
+ public List<org.lcsim.event.Cluster> getClusters()
+ {
+ return Collections.EMPTY_LIST;
+ }
+
+ public double[] getCovMatrix()
+ {
+ // Fixme: Does it make sense to return anything here?
+ return null;
+ }
+
+ public double getEnergy()
+ {
+ return fourVector.t();
+ }
+
+ public double getGoodnessOfPID()
+ {
+ return 0;
+ }
+
+ public double getMass()
+ {
+ return fourVector.magnitude();
+ }
+
+ public hep.physics.vec.Hep3Vector getMomentum()
+ {
+ return fourVector.v3();
+ }
+
+ public org.lcsim.event.ParticleID getParticleIDUsed()
+ {
+ return null;
+ }
+ public Vertex getStartVertex()
+ {
+ return null;
+ }
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.1 -r1.2
--- MCParticleExtended.java 3 Jun 2009 00:20:10 -0000 1.1
+++ MCParticleExtended.java 8 Jun 2009 18:29:05 -0000 1.2
@@ -6,24 +6,52 @@
package org.lcsim.contrib.Mbussonn;
import hep.physics.particle.properties.ParticleType;
+import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.HepLorentzVector;
import java.util.List;
+import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
import org.lcsim.event.base.BaseMCParticle;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
/**
*
* @author matthiasbussonnier
*/
public class MCParticleExtended implements MCParticle {
- public MCParticle mcp;
-
-
+ protected MCParticle mcp;
+ protected EventHeader event=null;
+ protected double bfield;
+ protected HelixParamCalculator helix=null;
public MCParticleExtended(MCParticle newMCP){
mcp= newMCP;
}
-
+ public MCParticleExtended(MCParticle newMCP , EventHeader event){
+ this.mcp = newMCP;
+ this.setEvent(event);
+ }
+/**
+ * Should allow hashmap to reconize MCParticleExtends as normal MCPs
+ */
+ @Override
+ public int hashCode(){
+ return mcp.hashCode();
+ }
+ @Override
+ public boolean equals(Object obj) {
+ if (obj == null) {
+ return false;
+ }
+ if (getClass() != obj.getClass()) {
+ return false;
+ }
+ final MCParticleExtended other = (MCParticleExtended) obj;
+ if (this.mcp != other.mcp && (this.mcp == null || !this.mcp.equals(other.mcp))) {
+ return false;
+ }
+ return true;
+ }
/**
* give the absolute value of the radial momentum
* @return
@@ -126,5 +154,24 @@
public HepLorentzVector asFourVector() {
return mcp.asFourVector();
}
+ /**
+ * return the distance of closest approch
+ * need 'event' to be set
+ * @return double dca
+ */
+ public double getDCA(){
+ if(helix == null){throw new IllegalArgumentException("can't get DCA, initialise MCPextend without event ?");}
+ return this.helix.getDCA();
+ }
+ public double getZ0(){
+ return this.helix.getZ0();
+ }
+ public void setEvent(EventHeader eventx){
+ if (eventx == null){ throw new IllegalArgumentException(" MCPE event point to null");}
+ this.event = eventx;
+ this.bfield = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 0.)).z();
+ this.helix = new HelixParamCalculator(mcp, bfield);
+ }
+
}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.1 -r1.2
--- JetFinder.java 3 Jun 2009 00:20:10 -0000 1.1
+++ JetFinder.java 8 Jun 2009 18:29:05 -0000 1.2
@@ -3,11 +3,27 @@
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.util.ArrayList;
import java.util.List;
+import java.util.Set;
+import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
+import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack;
+import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack.Ignore;
+import org.lcsim.contrib.sATLAS.TrackReconstructionDriver;
import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseRelationalTable;
import org.lcsim.event.util.JetDriver;
-import org.lcsim.mc.fast.MCFast;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -22,14 +38,16 @@
public JetFinder()
{
//super.add(new MCFast());
- super.add(new JetDriver());
+ super.add(new JetDriverExtended());
+ super.add(new TrackReconstructionDriver());
+
}
@Override
+ @SuppressWarnings("unchecked")
protected void process(EventHeader event)
{
super.process(event);
-
List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class,"Jets");
aida.cloud1D("nJets").fill(jets.size());
for (ReconstructedParticle jet : jets)
@@ -37,6 +55,150 @@
List<ReconstructedParticle> particlesInJet = jet.getParticles();
aida.cloud1D("nParticles").fill(particlesInJet.size());
}
+ RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+ List<Track> tracklist = event.getTracks();
+ for (LCRelation relation : mcrelations) {
+ hittomc.add(relation.getFrom(), relation.getTo());
+ }
+ //******ajoute une ralation entre reconstructed particle and mcp ********//
+ RelationalTable rc2mc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> rcmcfe = event.get(LCRelation.class, "rc2mc");
+ for (LCRelation relation : rcmcfe) {
+ rc2mc.add(relation.getFrom(), relation.getTo());
+ }
+
+ for (Track track : tracklist) {
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+ MCParticle mcp = tkanal.getMCParticle();
+ tkanal=null;
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double ptotal = Math.sqrt(px*px+py*py+pz*pz);
+ double ptcut = 50.;
+
+
+
+ if(track.getCharge()==0)
+ continue;
+ if( ptotal < ptcut)
+ continue;
+ if (mcp != null)
+ trktomc.add(track, mcp);
+ }
+ String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-JeffMarch26.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(sfile);
+ FindableTrack findable = new FindableTrack(event);
+
+ List<MCParticle> mcplist = event.getMCParticles();
+ System.out.println("mcp in event : "+mcplist.size());
+ //mcplist.clear();
+ for(ReconstructedParticle jet : jets){
+ List<ReconstructedParticle> listOfParticles = jet.getParticles();
+ for(ReconstructedParticle rcpInJet : listOfParticles){
+ MCParticle mcpp = (MCParticle) rc2mc.to(rcpInJet);
+ System.out.println("one more mcp");
+ //MCParticleExtended mcpe = new MCParticleExtended(mcpp,event);
+ //mcplist.add(mcpp);
+ }
+
+ }
+ System.out.println("mcp in jets : "+mcplist.size());
+ for(MCParticle mcpx : mcplist){
+ MCParticleExtended mcp = new MCParticleExtended(mcpx,event);
+ Set<Track> setOfTrack= trktomc.allTo(mcpx);
+ Set<Track> setOfTx = trktomc.allTo(mcp);
+ double ptotal = mcp.getPTotal();
+ double ptcut=50.;
+ double nrcp = rc2mc.allTo(mcpx).size();
+ aida.cloud1D("rcp per mcp").fill(nrcp);
+ if( ptotal < ptcut
+ //||(mcp.getPDGID()!=13 && mcp.getPDGID()!=-13)
+ ||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
+ ||Math.abs(mcp.getEta())>2.5)
+ {continue;}
+ else{
+ //if(!findable.isFindable(mcpx, slist)){
+ if(false){
+ System.out.println("muon non trouvable...");
+ List<Ignore> ignores = new ArrayList<Ignore>();
+ ignores.add(Ignore.NoZ0Cut);
+ if (findable.isFindable(mcpx, slist,ignores)){
+ System.out.println("not findable because of Z0");
+ continue;
+ }
+ ignores.add(Ignore.NoDCACut);
+ if (findable.isFindable(mcpx, slist,ignores)){
+ System.out.println("not findable because of DCA");
+ continue;
+ }
+ ignores.add(Ignore.NoMinHitCut);
+ if (findable.isFindable(mcpx, slist,ignores)){
+ System.out.println("not findable because of min hit");
+ continue;
+ }
+ ignores.add(Ignore.NoSeedCheck);
+ if (findable.isFindable(mcpx, slist,ignores)){
+ System.out.println("not findable because of seed check");
+ continue;
+ }
+ ignores.add(Ignore.NoConfirmCheck);
+ if (findable.isFindable(mcpx, slist,ignores)){
+ System.out.println("not findable because of no confirm check");
+ continue;
+ }
+ ignores=null;
+ System.out.println("not findable at all... why ?");
+ }
+ double weight=0.;
+
+ if(setOfTrack.size()>1){
+ aida.cloud1D("track with more than 1 mcp").fill(setOfTrack.size()); continue;
+ }
+ if(setOfTrack.size()==0) {
+ weight=0.;
+ int nhits = findable.LayersHit(mcpx);
+ aida.cloud1D("mcp/eta for non-finded particle").fill(mcp.getEta());
+ aida.cloud1D("mcp/theta for non-finded particle").fill(mcp.getTheta());
+ aida.cloud1D("nhits for non-finded muons").fill(nhits);
+ }
+ else if (setOfTrack.size()==1){
+ System.out.println("tracktomcp : "+setOfTx.size()+" Vs "+setOfTrack.size());
+ System.out.println("rcp associed with mcp "+rc2mc.allTo(mcpx).size());
+ weight=1;
+ Track track = setOfTrack.iterator().next();
+ double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double pt = Math.sqrt(px*px+py*py);
+ double ps = Math.sqrt(px*px+py*py+pz*pz);
+// double bfield = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 0.)).z();
+// HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
+// double d0mc = helix.getDCA();
+// double z0mc = helix.getZ0();
+ aida.cloud1D("track/d0").fill(d0);
+ aida.cloud1D("track/z0").fill(z0);
+ aida.cloud2D("track/do VS Z0").fill(z0, d0);
+ aida.cloud1D("track/tr momentum").fill(pt);
+ aida.cloud1D("track/scalar momentum").fill(ps);
+ aida.cloud1D("track/residual t momentum").fill(mcp.getPR()-pt);
+ aida.cloud1D("track/residual d0").fill(d0-mcp.getDCA());
+ aida.cloud1D("track/residual z0").fill(z0-mcp.getZ0());
+ aida.cloud1D("track/number of Tracker hits").fill(track.getTrackerHits().size());
+ aida.cloud1D("track/number of mcp hits").fill(findable.LayersHit(mcp));
+ aida.cloud1D("track/chiSquared").fill(track.getChi2());
+ }
+
+
+
+
+ }
+ }//ENDFOR
+ super.printStatistics(System.out);
}
