lcsim/src/org/lcsim/contrib/NickSinev/PixSim
diff -u -r1.1 -r1.2
--- CarrierPropagator.java 1 Jul 2008 21:40:03 -0000 1.1
+++ CarrierPropagator.java 25 Aug 2008 23:57:49 -0000 1.2
@@ -7,6 +7,7 @@
import hep.physics.vec.*;
import org.lcsim.geometry.Detector;
import org.lcsim.util.cache.*;
+import org.lcsim.conditions.*;
import java.util.zip.*;
/**
@@ -22,6 +23,8 @@
{
boolean debug = false;
boolean trace = false;
+ boolean init_OK = true;
+ boolean pconf_real = true;
String name="CCD_classic";
Random rnd = new Random();
RandomVector rnv = new RandomVector();
@@ -29,6 +32,7 @@
double ddt = 0.1E-10; // units for collection time in table
boolean use_table = false;
private DecimalFormat df = new DecimalFormat();
+ Interpolation3D inter = new Interpolation3D();
//
// pixel characteristics
//
@@ -152,6 +156,8 @@
double disp = 0.;
double delT=0.;
double expd = 0.;
+ int epOx = 0;
+ int epOy = 0;
/**
*
@@ -166,44 +172,67 @@
this(new PixelConfiguration());
}
+ public CarrierPropagator(String tfnam)
+ {
+ this((IPixelConfiguration) null);
+ use_table = useTable(tfnam);
+ if(!use_table) System.out.println("Could not read Carreir Propagator table "+tfnam+". Will not be able to proceed! ");
+ init_OK = use_table;
+ pconf_real=false;
+ }
+
+ public CarrierPropagator(String tfnam, String dnam)
+ {
+ this((IPixelConfiguration) null);
+ use_table = useTable(tfnam,dnam);
+ if(!use_table) System.out.println("Could not read Carreir Propagator table "+tfnam+" for detector "+dnam+". Will not be able to proceed! ");
+ init_OK = use_table;
+ pconf_real=false;
+ }
+
public CarrierPropagator(IPixelConfiguration pc)
{
- pconf=pc;
- name = pconf.getName();
- silicon = pconf.getBulk();
- psizex=pconf.getPixelSizeX();
- psizey=pconf.getPixelSizeY();
- epidep=pconf.getEpiThickness();
- charge = pconf.getCarrierCharge();
- EField = pconf.getEField();
- BField = pconf.getBField();
- List<SensorRegion> cregs = pconf.getCollectionRegions();
- if(cregs.size()>0)
+ init_OK = false;
+ if(pc != null)
{
- for(SensorRegion sr:cregs)
- collectors.add(sr);
- }
- List<SensorRegion> rregs = pconf.getReflectionRegions();
- if(rregs.size()>0)
- {
- for(SensorRegion sr:rregs)
- reflectors.add(sr);
- }
- List<SensorRegion> aregs = pconf.getAbsorbtionRegions();
- if(aregs.size()>0)
- {
- for(SensorRegion sr:aregs)
- absorbers.add(sr);
+ pconf=pc;
+ name = pconf.getName();
+ silicon = pconf.getBulk();
+ psizex=pconf.getPixelSizeX();
+ psizey=pconf.getPixelSizeY();
+ epidep=pconf.getEpiThickness();
+ charge = pconf.getCarrierCharge();
+ EField = pconf.getEField();
+ BField = pconf.getBField();
+ List<SensorRegion> cregs = pconf.getCollectionRegions();
+ if(cregs.size()>0)
+ {
+ for(SensorRegion sr:cregs)
+ collectors.add(sr);
+ }
+ List<SensorRegion> rregs = pconf.getReflectionRegions();
+ if(rregs.size()>0)
+ {
+ for(SensorRegion sr:rregs)
+ reflectors.add(sr);
+ }
+ List<SensorRegion> aregs = pconf.getAbsorbtionRegions();
+ if(aregs.size()>0)
+ {
+ for(SensorRegion sr:aregs)
+ absorbers.add(sr);
+ }
+ uniformBField=pconf.bFieldIsUniform();
+ unidirectEField = pconf.eFieldIsUnidirect();
+ if(unidirectEField) eFieldDirection=pconf.getEFieldDirection();
+ BField.getField(cp,bfield);
+ lzUnit = silicon.getLorentzVector(bfield,eFieldDirection,charge);
+ lorentzVSet=true;
+ Hep3Vector bFieldDirection = VecOp.unit(bfield);
+ befc = (VecOp.cross(bFieldDirection,eFieldDirection)).magnitude();
+ setLimits();
+ init_OK = true;
}
- uniformBField=pconf.bFieldIsUniform();
- unidirectEField = pconf.eFieldIsUnidirect();
- if(unidirectEField) eFieldDirection=pconf.getEFieldDirection();
- BField.getField(cp,bfield);
- lzUnit = silicon.getLorentzVector(bfield,eFieldDirection,charge);
- lorentzVSet=true;
- Hep3Vector bFieldDirection = VecOp.unit(bfield);
- befc = (VecOp.cross(bFieldDirection,eFieldDirection)).magnitude();
- setLimits();
}
public void setTrace(boolean yes)
@@ -220,6 +249,8 @@
public IPixelConfiguration getPixelConfiguration() { return pconf; }
+ public double[] getPixelDimensions() { double[] pd = new double[3]; pd[0]=psizex; pd[1]=psizey; pd[2]=epidep; return pd; }
+
private void setLimits()
{
x0= -psizex/2.;
@@ -292,55 +323,24 @@
use_table = false;
String dfname = "/detectors/"+dname+".zip";
String enam = "PixilatedSensorTables/"+nam;
- URL url = null;
- try
- {
- url = new URL("http","www.lcsim.org",dfname);
+ ConditionsManager cm = ConditionsManagerImplementation.defaultInstance();
+ try {cm.setDetector(dname, 0); }
+ catch(Exception e) {
+ throw new RuntimeException("Could not find conditions information for detector "+dname);
}
- catch(MalformedURLException e)
- {
- }
- finally
- {
- }
- FileCache ucache = null;
- try {ucache = new FileCache();} catch(IOException e) {} finally{}
- File lroot = new File(ucache.getCacheRoot(),".lcsim");
- File lrootc = new File(lroot,"cache");
- FileCache fcache = null;
- if(lrootc != null){ try {fcache = new FileCache(lrootc);} catch(IOException e) {} finally{} }
- File path = null;
- if(url != null) { try{ path = fcache.getCachedFile(url); } catch(IOException e){} finally {} }
- BufferedReader reader = null;
- if(path != null)
- {
-// System.out.println("Found: "+path.getAbsolutePath());
- ZipFile zipf = null;
- try { zipf = new ZipFile(path);} catch(Exception e) {} finally {}
- Enumeration<?extends ZipEntry> zcont = zipf.entries();
- ZipEntry ztab = null;
- while(zcont.hasMoreElements())
- {
- ZipEntry ze = (ZipEntry) zcont.nextElement();
-// System.out.println("Entry name: "+ze.getName());
- if(ze.getName().equals(enam))
- {
- System.out.println("Found table "+enam+" in detector "+dname+" configuration zip file");
- ztab = ze;
- }
- }
- if(ztab != null)
- {
- InputStream tinps = null;
- try { tinps = zipf.getInputStream(ztab);} catch(IOException e) {} finally {}
- if(tinps != null) reader = new BufferedReader(new InputStreamReader(tinps));
- }
+ BufferedReader reader;
+ try {
+ InputStream is = cm.getRawConditions(enam).getInputStream();
+ reader = new BufferedReader(new InputStreamReader(is));
+ System.out.println("Found table "+enam+" in detector "+dname+" conditions");
+ } catch (Exception e) {
+ throw new RuntimeException(e.getMessage());
}
if(reader != null) use_table = readTable(reader);
return use_table;
}
- public void doNotUseTable() { use_table=false; }
+ public void doNotUseTable() { use_table=false; init_OK = pconf_real;}
public boolean readTable(String fnam)
@@ -465,10 +465,20 @@
tpsz=dpar[2];
System.out.println("Pixel dimensions: "+tpsx+" x "+tpsy+" x "+tpsz);
boolean dmatch=true;
- if(Math.abs(1.-psizex/tpsx) > 0.0001) dmatch=false;
- if(Math.abs(1.-psizey/tpsy) > 0.0001) dmatch=false;
- if(Math.abs(1.-epidep/tpsz) > 0.0001) dmatch=false;
- if(!dmatch) System.out.println("Pixel dimensions in table do not match dimensions in PixelConfiguration!");
+ if(pconf != null)
+ {
+ if(Math.abs(1.-psizex/tpsx) > 0.0001) dmatch=false;
+ if(Math.abs(1.-psizey/tpsy) > 0.0001) dmatch=false;
+ if(Math.abs(1.-epidep/tpsz) > 0.0001) dmatch=false;
+ if(!dmatch) System.out.println("Pixel dimensions in table do not match dimensions in PixelConfiguration!");
+ }
+ else
+ {
+ pconf = new PixelConfiguration(tpsx,tpsy,tpsz);
+ psizex = tpsx;
+ psizey = tpsy;
+ epidep = tpsz;
+ }
if(dmatch)
{
h_table = new float[tnx][tny][tnz][tnpx][tnpy];
@@ -539,9 +549,12 @@
} // end of for loop
} // end of while loop
tab_read=true;
- tdx=psizex/tnx;
- tdy=psizey/tny;
- tdz=epidep/tnz;
+ tdx=psizex/(tnx-1);
+ tdy=psizey/(tny-1);
+ tdz=epidep/(tnz-1);
+ x0 = - psizex/2.;
+ y0 = - psizey/2.;
+ z0 = 0.;
for(int i=0; i<tnx; i++)
{
for(int j=0; j<tny; j++)
@@ -606,9 +619,9 @@
int ofy = (npy-1)/2;
int pxi = 0;
int pyi = 0;
- double ddx = psizex/nx;
- double ddy = psizey/ny;
- double ddz = epidep/nz;
+ double ddx = psizex/(nx-1);
+ double ddy = psizey/(ny-1);
+ double ddz = epidep/(nz-1);
double x0 = -psizex/2;
double y0 = -psizey/2;
double z0 = 0.;
@@ -631,7 +644,7 @@
}
for(int n=0; n<nsmp; n++)
{
- if(n%nrep==0) System.out.println("Generated "+n+" tests in each point");
+ if(n%nrep==0) System.out.println("Generated "+n+" tests in each point (out of "+nsmp+" requested");
for(int i=0; i<nx; i++)
{
x=x0+i*ddx;
@@ -685,18 +698,16 @@
chunks.clear();
startP.setV(sp.x(),sp.y(),sp.z());
equivalentPoint(startP,ep);
- EField.getField(ep,efield);
- E=efield.magnitude();
-// if((!use_table) || (E>0.1))
if(!use_table)
{
+ EField.getField(ep,efield);
+ E=efield.magnitude();
boolean newcc = true;
for(int i=0; i< npairs; i++)
{
newcc=true;
startP.setV(sp.x(),sp.y(),sp.z());
- if(use_table && (i%5 == 0)) transport();
- if(!use_table) transport();
+ transport();
if(status==1)
{
for(ChargeChunk cc:chunks)
@@ -715,7 +726,7 @@
}
}
}
- }
+ } // end if(!use_table)
else
{
int cxi = (tnpx-1)/2;
@@ -728,17 +739,23 @@
double x=ep.x();
double y=ep.y();
double z=ep.z();
+ int tpx = 0;
+ int tpy = 0;
+ int tpz = 0;
tix = (int) Math.floor((x-x0)/tdx);
- if(tix>=tnx) tix=tnx-1;
+ if(tix>=tnx-1) tix=tnx-2;
if(tix<0) tix=0;
+ tpx=tix+1;
tddx = (x-x0)-tix*tdx;
tiy = (int) Math.floor((y-y0)/tdy);
if(tiy<0) tiy=0;
- if(tiy>=tny) tiy=tny-1;
+ if(tiy>=tny-1) tiy=tny-2;
+ tpy=tiy+1;
tddy = (y-y0)-tiy*tdy;
tiz = (int) Math.floor((z-z0)/tdz);
- if(tiz>=tnz) tiz=tnz-1;
+ if(tiz>=tnz-1) tiz=tnz-2;
if(tiz<0) tiz=0;
+ tpz=tiz+1;
tddz = (z-z0)-tiz*tdz;
if((tix != tixp) || (tiy != tiyp) || (tiz != tizp))
{
@@ -760,13 +777,34 @@
}
}
}
+ double[] fpr = new double[8];
+ double[] fct = new double[8];
+ double[] dsp = new double[3];
+ dsp[0]=tddx/tdx;
+ dsp[1]=tddy/tdy;
+ dsp[2]=tddz/tdz;
for(int i=0; i<tnpx; i++)
{
for(int j=0; j<tnpy; j++)
{
- pprob[i][j]=h_tc[i][j]+tddx*hxd_tc[i][j]+tddy*hyd_tc[i][j]+tddz*hzd_tc[i][j];
- ctval[i][j]=t_tc[i][j]+tddx*txd_tc[i][j]+tddy*tyd_tc[i][j]+tddz*tzd_tc[i][j];
- if(ctval[i][j] < 0.) ctval[i][j]=0;
+ fpr[0]=h_table[tix][tiy][tiz][i][j];
+ fpr[1]=h_table[tix][tiy][tpz][i][j];
+ fpr[2]=h_table[tix][tpy][tiz][i][j];
+ fpr[3]=h_table[tix][tpy][tpz][i][j];
+ fpr[4]=h_table[tpx][tiy][tiz][i][j];
+ fpr[5]=h_table[tpx][tiy][tpz][i][j];
+ fpr[6]=h_table[tpx][tpy][tiz][i][j];
+ fpr[7]=h_table[tpx][tpy][tpz][i][j];
+ pprob[i][j]=inter.interpolate(fpr,dsp);
+ fct[0]=t_table[tix][tiy][tiz][i][j];
+ fct[1]=t_table[tix][tiy][tpz][i][j];
+ fct[2]=t_table[tix][tpy][tiz][i][j];
+ fct[3]=t_table[tix][tpy][tpz][i][j];
+ fct[4]=t_table[tpx][tiy][tiz][i][j];
+ fct[5]=t_table[tpx][tiy][tpz][i][j];
+ fct[6]=t_table[tpx][tpy][tiz][i][j];
+ fct[7]=t_table[tpx][tpy][tpz][i][j];
+ ctval[i][j]=inter.interpolate(fct,dsp);
}
}
for(int i=0; i<tnpx; i++)
@@ -800,7 +838,7 @@
}
}
}
- }
+ } // end else
return chunks;
}
@@ -811,17 +849,23 @@
double x=ep.x();
double y=ep.y();
double z=ep.z();
+ int tpx = 0;
+ int tpy = 0;
+ int tpz = 0;
tix = (int) Math.floor((x-x0)/tdx);
- if(tix>=tnx) tix=tnx-1;
+ if(tix>=tnx-1) tix=tnx-2;
if(tix<0) tix=0;
+ tpx=tix+1;
tddx = (x-x0)-tix*tdx;
tiy = (int) Math.floor((y-y0)/tdy);
if(tiy<0) tiy=0;
- if(tiy>=tny) tiy=tny-1;
+ if(tiy>=tny-1) tiy=tny-2;
+ tpy=tiy+1;
tddy = (y-y0)-tiy*tdy;
tiz = (int) Math.floor((z-z0)/tdz);
- if(tiz>=tnz) tiz=tnz-1;
+ if(tiz>=tnz-1) tiz=tnz-2;
if(tiz<0) tiz=0;
+ tpz=tiz+1;
tddz = (z-z0)-tiz*tdz;
if((tix != tixp) || (tiy != tiyp) || (tiz != tizp))
{
@@ -843,12 +887,38 @@
}
}
}
+ double[] fpr = new double[8];
+ double[] fct = new double[8];
+ double[] dsp = new double[3];
+ dsp[0]=tddx/tdx;
+ dsp[1]=tddy/tdy;
+ dsp[2]=tddz/tdz;
for(int i=0; i<tnpx; i++)
{
for(int j=0; j<tnpy; j++)
{
+ fpr[0]=h_table[tix][tiy][tiz][i][j];
+ fpr[1]=h_table[tix][tiy][tpz][i][j];
+ fpr[2]=h_table[tix][tpy][tiz][i][j];
+ fpr[3]=h_table[tix][tpy][tpz][i][j];
+ fpr[4]=h_table[tpx][tiy][tiz][i][j];
+ fpr[5]=h_table[tpx][tiy][tpz][i][j];
+ fpr[6]=h_table[tpx][tpy][tiz][i][j];
+ fpr[7]=h_table[tpx][tpy][tpz][i][j];
+ pprob[i][j]=inter.interpolate(fpr,dsp);
+ fct[0]=t_table[tix][tiy][tiz][i][j];
+ fct[1]=t_table[tix][tiy][tpz][i][j];
+ fct[2]=t_table[tix][tpy][tiz][i][j];
+ fct[3]=t_table[tix][tpy][tpz][i][j];
+ fct[4]=t_table[tpx][tiy][tiz][i][j];
+ fct[5]=t_table[tpx][tiy][tpz][i][j];
+ fct[6]=t_table[tpx][tpy][tiz][i][j];
+ fct[7]=t_table[tpx][tpy][tpz][i][j];
+ ctval[i][j]=inter.interpolate(fct,dsp);
+/*
pprob[i][j]=h_tc[i][j]+tddx*hxd_tc[i][j]+tddy*hyd_tc[i][j]+tddz*hzd_tc[i][j];
- ctval[i][j]=t_tc[i][j]+tddx*txd_tc[i][j]+tddy*tyd_tc[i][j]+tddz*tzd_tc[i][j];
+ ctval[i][j]=t_tc[i][j]+tddx*txd_tc[i][j]+tddy*tyd_tc[i][j]+tddz*tzd_tc[i][j]; */
+
if(ctval[i][j] < 0.) ctval[i][j]=0;
}
}
@@ -922,7 +992,7 @@
{
iy++;
if(iy==tnpy) {iy=0; ix++;}
- }
+ }
if(ix==tnpx) return false;
double xx=cp.x()+psizex*(ix-cxi);
double yy=cp.y()+psizey*(iy-cyi);
@@ -976,6 +1046,8 @@
private void transport()
{
+ if(!init_OK) { finalstat[2]=0; return; }
+ if(use_table) trace=false;
if(trace) p_trace.clear();
cp.setV(startP.x(),startP.y(),startP.z());
if(trace) { Hep3Vector trckp = new BasicHep3Vector(cp.x(),cp.y(),cp.z()); p_trace.add(trckp); }
@@ -990,9 +1062,6 @@
if(!finish)
{
equivalentPoint(cp,ep);
- EField.getField(ep,efield);
- E=efield.magnitude();
-// if((E < 0.1) && use_table)
if(use_table)
{
collected = lookupTable();
@@ -1105,8 +1174,8 @@
}
}
}
- }
- }
+ } // end while(!finish)
+ } // end if(!finish)
offsetX=0;
offsetY=0;
status = 0;
@@ -1119,33 +1188,13 @@
dx=cp.x()+psizex/2.;
dy=cp.y()+psizey/2.;
dz=cp.z();
-/* while(dx > psizex/2.)
- {
- offsetX++;
- dx-=psizex;
- }
- while(dx < -psizex/2.)
- {
- offsetX--;
- dx+=psizex;
- }
- while(dy > psizey/2.)
- {
- offsetY++;
- dy-=psizey;
- }
- while(dy < -psizey/2.)
- {
- offsetY--;
- dy+=psizey;
- } */
offsetX = (int) Math.floor(dx/psizex);
offsetY = (int) Math.floor(dy/psizey);
}
- finalstat[0]=offsetX;
- finalstat[1]=offsetY;
+ finalstat[0]=offsetX+epOx;
+ finalstat[1]=offsetY+epOy;
finalstat[2]=status;
finalstat[3]=nst;
}
@@ -1161,5 +1210,7 @@
public void equivalentPoint(Hep3Vector p, BasicHep3Vector e)
{
e.setV(p.x()-(Math.floor(p.x()/psizex+0.5))*psizex,p.y()-(Math.floor(p.y()/psizey+0.5))*psizey,p.z());
+ epOx = (int) Math.floor((p.x()-ep.x())/psizex);
+ epOy = (int) Math.floor((p.y()-ep.y())/psizey);
}
}
\ No newline at end of file
lcsim/src/org/lcsim/contrib/NickSinev/PixSim
diff -u -r1.2 -r1.3
--- PixilatedSensorManager.java 15 Jul 2008 01:18:29 -0000 1.2
+++ PixilatedSensorManager.java 25 Aug 2008 23:57:49 -0000 1.3
@@ -13,13 +13,14 @@
import org.lcsim.detector.solids.*;
import java.text.*;
import java.util.*;
+import org.lcsim.util.lcio.LCIOConstants;
/**
* This class are designed for creating a list of pixilated sensor
* for each event (only sensors which have hits in them are created)
*
* @author Nick Sinev
- * @version $Id: PixilatedSensorManager.java,v 1.2 2008/07/15 01:18:29 sinev Exp $
+ * @version $Id: PixilatedSensorManager.java,v 1.3 2008/08/25 23:57:49 sinev Exp $
*/
public class PixilatedSensorManager extends Driver
@@ -760,13 +761,18 @@
} // end if(evhits != null)
List<TrackerHit> ebhits = new ArrayList<TrackerHit>();
List<TrackerHit> eehits = new ArrayList<TrackerHit>();
+ List<RawTrackerHit> rawhits = new ArrayList<RawTrackerHit>();
for(PixilatedSensor psn:psbrl)
{
psn.processEvent(phys_bc);
IDetectorElement de = psn.getParent();
IReadout ro = de.getReadout();
List<TrackerHit> sthts = ro.getHits(TrackerHit.class);
- for(TrackerHit ht:sthts) ebhits.add(ht);
+ for(TrackerHit ht:sthts)
+ {
+ ebhits.add(ht);
+ rawhits.addAll(ht.getRawHits());
+ }
}
// System.out.println("Barrel hits processed!");
for(PixilatedSensor psn:psecpl)
@@ -775,7 +781,11 @@
IDetectorElement de = psn.getParent();
IReadout ro = de.getReadout();
List<TrackerHit> sthts = ro.getHits(TrackerHit.class);
- for(TrackerHit ht:sthts) eehits.add(ht);
+ for(TrackerHit ht:sthts)
+ {
+ eehits.add(ht);
+ rawhits.addAll(ht.getRawHits());
+ }
}
// System.out.println("Endcap plus hits processed!");
for(PixilatedSensor psn:psecml)
@@ -784,11 +794,16 @@
IDetectorElement de = psn.getParent();
IReadout ro = de.getReadout();
List<TrackerHit> sthts = ro.getHits(TrackerHit.class);
- for(TrackerHit ht:sthts) eehits.add(ht);
+ for(TrackerHit ht:sthts)
+ {
+ eehits.add(ht);
+ rawhits.addAll(ht.getRawHits());
+ }
}
// System.out.println("Endcap minus hits processed!");
event.put("RecVtxBarrHits",ebhits,TrackerHit.class,0);
event.put("RecVtxEndcapHits",eehits,TrackerHit.class,0);
+ event.put("RecVtxRawHits",rawhits,RawTrackerHit.class,(1 << LCIOConstants.RTHBIT_HITS));
} // end function process