lcsim/sandbox/TimBarklow
diff -N VertexFinderDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ VertexFinderDriver.java 30 Oct 2007 08:10:25 -0000 1.1
@@ -0,0 +1,221 @@
+package org.lcsim.contrib.timb;
+
+import org.lcsim.mc.fast.reconstructedparticle.MCFastReconstructedParticle;
+import org.lcsim.mc.fast.reconstructedparticle.MCFastParticleID;
+import hep.physics.jet.JadeEJetFinder;
+import hep.physics.jet.JetFinder;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.HepLorentzVector;
+import hep.physics.vec.BasicHepLorentzVector;
+import hep.physics.vec.VecOp;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticlePropertyProvider;
+import hep.physics.particle.properties.ParticleType;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Vector;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Track;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import hep.lcd.vertexing.zvtop.ZvTrackList;
+import hep.lcd.vertexing.zvtop.ZvTrack;
+import hep.lcd.vertexing.zvtop.ZvTopVertexer;
+import hep.lcd.vertexing.zvtop.ZvVertexList;
+import hep.lcd.vertexing.zvtop.ZvVertex;
+import hep.lcd.vertexing.zvtop.ZvBField;
+import static java.lang.Math.pow;
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
+
+//
+// Find vertices using ZVTOP and modify MCFastReconstructedParticles accordingly
+//
+public class VertexFinderDriver extends Driver
+{
+ private static int nevent = 0 ;
+ // private static final int nevent_min = 11 ;
+ private static final int nevent_min = 0 ;
+ private static final int nprint_max = 2 ;
+ private static final int n_chg_min = 2 ;
+ private static final boolean convert_mm_to_cm = true;
+ private static final double cmTOmm = 10.;
+ private static final double av = 0.1;
+ private static final double ov = 0.1;
+ private static final Hep3Vector IP = new BasicHep3Vector(0,0,0);
+ private static final String defaultInCollectionName = "MCFastReconstructedParticles";
+ private static final String defaultJetCollectionName = "Jets";
+ private static final String defaultOutCollectionName = "MCFastZVReconstructedParticles";
+ private String inCollectionName = defaultInCollectionName;
+ private String jetCollectionName = defaultJetCollectionName;
+ private String outCollectionName = defaultOutCollectionName;
+ private ZvTopVertexer vertexer = VertexFinderFuncs.initZvTopVertexer();
+ private static final ParticleType k0s;
+ private static final ParticleType lambdab0;
+ private static ParticlePropertyProvider ppp;
+ static {
+ ppp = ParticlePropertyManager.getParticlePropertyProvider();
+ k0s = ppp.get(310);
+ lambdab0 = ppp.get(5122);
+ }
+
+ /** Creates a new instance of VertexFinderDriver */
+ public VertexFinderDriver()
+ {
+ }
+ public VertexFinderDriver(String inCollectionName, String jetCollectionName, String outCollectionName)
+ {
+ this.inCollectionName = inCollectionName;
+ this.jetCollectionName = jetCollectionName;
+ this.outCollectionName = outCollectionName;
+ }
+
+ public String getInCollectionName()
+ {
+ return inCollectionName;
+ }
+
+ public String getJetCollectionName()
+ {
+ return jetCollectionName;
+ }
+
+ public String getOutCollectionName()
+ {
+ return outCollectionName;
+ }
+
+ protected void process(EventHeader event)
+ {
+ nevent++;
+ double bField = event.getDetector().getFieldMap().getField(IP.v())[2];
+ ZvBField.setBZ(bField);
+ boolean zp = nevent <= nprint_max;
+ // Find the input reconstructed Particles
+ List<ReconstructedParticle> list_in;
+ list_in = event.get(ReconstructedParticle.class,inCollectionName);
+ // Find the jet reconstructed Particles
+ List<ReconstructedParticle> list_jet;
+ list_jet = event.get(ReconstructedParticle.class,jetCollectionName);
+
+ // create output list
+ List<ReconstructedParticle> list_out = new ArrayList<ReconstructedParticle>(list_in);
+
+ if(zp) System.out.println(" VertexFinderDriver nevent= "+nevent+" list_in.size= "+list_in.size()+" list_jet.size= "+list_jet.size());
+ for(ReconstructedParticle rpj : list_jet){
+ Map<hep.lcd.mc.fast.tracking.ReconTrack,ReconstructedParticle> m_tr = new HashMap<hep.lcd.mc.fast.tracking.ReconTrack,ReconstructedParticle>();
+ Vector<hep.lcd.mc.fast.tracking.ReconTrack> vector_heprt = new Vector<hep.lcd.mc.fast.tracking.ReconTrack>();
+ List<ReconstructedParticle> list_rp = rpj.getParticles();
+ for(ReconstructedParticle rp : list_rp) {
+ List<Track> list_trk = rp.getTracks();
+ // System.out.println(" VertexFinderDriver list_trk.size= "+list_trk.size()+" list_rp.size= "+list_rp.size());
+ for(Track trk : list_trk){
+ double [] m_parm = new double [5];
+ for(int j=0 ; j<5 ; j++) m_parm[j] = trk.getTrackParameter(j);
+ SymmetricMatrix sym_err = trk.getErrorMatrix();
+ double [][] m_err = new double [5][5];
+ for(int j=0 ; j<5 ; j++) for(int k=0 ; k<5 ; k++) m_err[k][j]=sym_err.e(k,j);
+ if(zp) {
+ for(int j=0; j<5 ; j++) System.out.println(" VertexFinderDriver before any conversion j= "+j+" m_parm[j]= "+m_parm[j]);
+ System.out.println(" VertexFinderDriver trk.getCharge= "+trk.getCharge()+" trk.getMomentum= "+Arrays.toString(trk.getMomentum()));
+ }
+ m_parm[2]=-m_parm[2]; // hep.lcd as opposite sign convention for omega
+ if(convert_mm_to_cm) {
+ m_parm[0]*=av;
+ m_parm[2]/=ov;
+ m_parm[3]*=av;
+ m_err[0][0]*=pow(av,2);
+ m_err[0][1]*=av;
+ m_err[0][2]=av*m_err[0][2]/ov;
+ m_err[0][3]*=pow(av,2);
+ m_err[0][4]*=av;
+ m_err[1][2]/=ov;
+ m_err[1][3]*=av;
+ m_err[2][2]/=pow(ov,2);
+ m_err[2][3]=av*m_err[2][3]/ov;
+ m_err[2][4]/=ov;
+ m_err[3][3]*=pow(av,2);
+ m_err[3][4]*=av;
+ for(int j=0 ; j<4 ; j++) for(int k=j+1 ; k<5 ; k++) m_err[k][j]=m_err[j][k] ;
+ }
+ hep.lcd.mc.fast.tracking.DocaTrackParameters doca = new hep.lcd.mc.fast.tracking.DocaTrackParameters(m_parm,m_err);
+ hep.lcd.mc.fast.tracking.ReconTrack hrt = new hep.lcd.mc.fast.tracking.ReconTrack(doca);
+ vector_heprt.addElement(hrt);
+ m_tr.put(hrt,rp);
+ // System.out.println(" VertexFinderDriver heprt.getCharge= "+((hep.lcd.mc.fast.tracking.ReconTrack)(vector_heprt.lastElement())).getCharge()
+ // +" heprt.getMomentum= "+Arrays.toString(((hep.lcd.mc.fast.tracking.ReconTrack)(vector_heprt.lastElement())).getMomentum()));
+ // System.out.println(" VertexFinderDriver vector_heprt.size= "+vector_heprt.size()+" list_rp.size= "+list_rp.size());
+ }
+ }
+ if(zp) System.out.println(" VertexFinderDriver vector_heprt.size= "+vector_heprt.size()+" list_rp.size= "+list_rp.size());
+ if(vector_heprt.size()>=n_chg_min) {
+ ZvTrackList jetTrackList = new ZvTrackList(vector_heprt.elements());
+ if(zp) for(int j=0; j<jetTrackList.getNTracks() ; j++) System.out.println(" VertexFinderDriver zvtrk.getCharge= "+jetTrackList.trackAt(j).getCharge()
+ +" zvtrk.getMomentum= "+Arrays.toString(jetTrackList.trackAt(j).getMomentum()));
+ if(nevent >= nevent_min) {
+ ZvVertexList vtxList = VertexFinderFuncs.calcSeedVertex(vertexer,jetTrackList);
+ int nVtx = vtxList.getNVertices();
+ if(zp) System.out.println(" From VertexFinderDriver nevent= "+nevent);
+ for(int i=0; i< nVtx ; i++) {
+ ZvVertex vtx = vtxList.vertexAt(i);
+ double [] vpos = vtx.getVertexPos();
+ for(int j=0 ; j<vpos.length ; j++) vpos[j] *= cmTOmm;
+ Hep3Vector bpos = new BasicHep3Vector(vpos);
+ if(zp) System.out.println(" i= "+i+" vtx.getNumTracks= "+vtx.getNumTracks()+" vtx.getSumCharge= "+vtx.getSumCharge()+" vpos[0]= "+vpos[0]+" vpos[1]= "+vpos[1]+" vpos[2]= "+vpos[2]+" mag= "+bpos.magnitude());
+ if(i>0) {
+ ZvTrackList zvtl = vtx.getVtxTracks();
+ HepLorentzVector vtp = new BasicHepLorentzVector();
+ for(int j=0; j<zvtl.getNTracks(); j++){
+ ZvTrack zvt = zvtl.trackAt(j);
+ hep.lcd.mc.fast.tracking.ReconTrack trk = zvt.getOrigReconTrack();
+ int lind = list_in.indexOf(m_tr.get(trk));
+ ReconstructedParticle rp = list_in.get(lind);
+ Hep3Vector zvtmom = new BasicHep3Vector(zvt.getMomentum());
+ vtp = VecOp.add(vtp,new BasicHepLorentzVector(sqrt(zvtmom.magnitudeSquared()+pow(rp.getMass(),2)),zvtmom));
+ }
+ MCFastReconstructedParticle rplambdab0 = new MCFastReconstructedParticle(vpos,(vtp.v3()).v(),vtp.magnitude(),vtx.getSumCharge(),lambdab0);
+ for(int j=0; j<zvtl.getNTracks(); j++){
+ ZvTrack zvt = zvtl.trackAt(j);
+ hep.lcd.mc.fast.tracking.ReconTrack trk = zvt.getOrigReconTrack();
+ int lind = list_in.indexOf(m_tr.get(trk));
+ // for(int k=0; k<5 ; k++) System.out.println(" VertexFinderFuncs j= "+j+" k= "+k+" m_parm[k]= "+(trk.getTrackParameters())[k]);
+ if(zp) System.out.println(" VertexFinderFuncs lind= "+lind+" i= "+i+" j= "+j+" zvt.getCharge= "+zvt.getCharge()+" trk.getMomentum= "+Arrays.toString(trk.getMomentum()));
+ ReconstructedParticle rp = list_in.get(lind);
+ int typeold = rp.getType();
+ int typenew = 9900000+abs(typeold);
+ if(typeold < 0 ) typenew = -typenew;
+ rp.addParticleID(new MCFastParticleID(ppp.get(typenew)));
+ if(zp) System.out.println(" VertexFinderFuncs i= "+i+" j= "+j+" reconst part charge= "+rp.getCharge()+" rp.getMomentum= "+Arrays.toString((rp.getMomentum()).v())
+ +" typeold= "+typeold+" typenew= "+typenew+" rp.getType= "+rp.getType());
+ rplambdab0.addParticle(rp);
+
+ }
+ list_out.add(rplambdab0);
+ List<ReconstructedParticle> list_rpinrplambdab0 = rplambdab0.getParticles();
+ if(zp) System.out.println(" VertexFinderFuncs i= "+i+" rplambdab0.charge = "+rplambdab0.getCharge()+" mass= "+rplambdab0.getMass()+" Momentum= "+Arrays.toString((rplambdab0.getMomentum()).v())
+ +" ref= "+Arrays.toString((rplambdab0.getReferencePoint()).v())+" type= "+rplambdab0.getType()+" list_rpinrplambdab0.size= "+list_rpinrplambdab0.size());
+ }
+ }
+ if(zp) {
+ for(ReconstructedParticle rp : list_rp) {
+ System.out.println(" ");
+ List<Track> list_trk = rp.getTracks();
+ for(Track trk : list_trk){
+ double [] m_parm = trk.getTrackParameters();
+ for(int j=0; j<5 ; j++) System.out.println(" VertexFinderDriver after filling recon particles j= "+j+" m_parm[j]= "+m_parm[j]);
+ }
+ }
+ }
+ }
+ }
+ }
+ event.put(outCollectionName, list_out, ReconstructedParticle.class, 0);
+ }
+}
