+import java.util.*; +import java.io.*; +import org.lcsim.util.aida.AIDA; +import hep.aida.ITree; +import java.util.Random; +import org.lcsim.recon.vertexing.pixsim.*; +import org.lcsim.geometry.FieldMap; +import hep.physics.vec.*; +import hep.aida.*; + +/* + * CarrierPropagatorMakeTable.java + * + * Created on July 31, 2005, 3:03 PM + * + */ +/** + * + * @author Nick Sinev + */ +public class CarrierPropagatorMakeTable +{ + /* + * Here is the example of how to create CarrierPropagator tables. + * The CarrierPropagator can simulate charge collection in the + * pixilated sensor by simulating the movements of every charge + * carreier, generated by ionizing particle on it's path through + * sensor. But this takes a long time, and not sutable for pixilated + * sensor digitization simulation in the event processing framework + * of org.lcsim, where thousands of track segments may be observed in + * the vertex detector sensors in every event. For such application + * it is more rational to do full simulation of charge carrier propagation + * from evry point inside sensitive volume of silicon, and determine + * probability of collecting this charge by central or any of neighboring + * pixels. Such probabilities can be saved in the form of tables in + * special file, and be used for faster simulation of charge collection. + * Creating of such table is very lenghtly process, it depends on how + * fine grid of sample points inside pixel we want to have, and how + * many times we want to propagate particle from the same point to + * get estimation of collection probabilities. + * For 16x16x12 micron pixel, with 0.5 micron spacing between sample + * points (in either direction) and 1000 tests in each point we will + * create grid of 24576 points and generate 24576000 charge carrier tracks. + * On my DELL Optiflex 745 desctop, running 2 CPUs at 2.66 GHz with 3.25 GB + * of memory, it took about 5 hours to make such table. + * Before investing such amount of CPU time into table creation, it makes + * sense to generate more accurate electric field map (not the coarse one, + * used in previous examples). Here I am using map with 0.2 micron spacing. + * I did not put such map into CVS, as it is rather large. So, you will + * need to create your own one, and edit this file for proper map file name. + */ + public static void main(String[] args) throws IOException + { + + /* the name of the output file for CarrierPropagator table: */ + String oname = "Chronopix_20x20x20_hr_B5_EC_propagator.dat"; + + /* The name of electric field map file (converted from TCAD into regular grid map + * see TCADReadSaveExample.java for example of how to do such conversion.) + */ + String efname = "thick3d_20hr_2V_emap"; + String tcmname = efname+".txt"; + /* We are going to use TCAD Field Map for electric field */ + TCADFieldMap tc_fm= new TCADFieldMap(efname); + + + /* Reading in electric field. If there are problems - nothing we can do here. */ + if(!tc_fm.readMap(tcmname)) return; + System.out.println("Found normalized field map file: "+tcmname); + + /* + * Dimensions of field map can tell us about pixel dimensions, + * as correctly saved field map should be exactly of the size of one pixel. + */ + + double[] flims = tc_fm.getMapLimits(); + double pszx = flims[1]-flims[0]; + double pszy = flims[3]-flims[2]; + double pszz = flims[5]-flims[4]; + System.out.println("Field limits: X"+flims[0]+" "+flims[1]+" Y: "+flims[2]+" "+flims[3]+" Z: "+flims[4]+" "+flims[5]); + + /* + * We need to define magnetic field also. + * This should be in the sensor coordinate system, which is same as global for EC + * (inversed B direction for z-minus ec) + * For barrel sensors y axis is directed along global z, so for them Bfield should be (0.,5.,0.) + * Bfield in endcap would not do anything if e-field is always parallel to z axis (as in case of CCD), + * but for complex e-field map, like in active pixel devices, Bfield may have some effect + */ + Hep3Vector B = new BasicHep3Vector(0.,0.,0.); + NamedFieldMap Bfield = new UniformFieldMap(B); + + /* tell them, that this is silicon detector ! */ + Medium bulk=new Silicon(); + + /* PixelConfiguration objec may have name, which the first argument in following call sets. + * Currently this name is not used anywhere, except for interactive info output + */ + PixelConfiguration pcf = new PixelConfiguration("Chrono_20x20x20hr_2V_EC",bulk,tc_fm,Bfield,pszx,pszy,pszz); + /* announce that B field is uniform - to speed up Lorentz angle calculations */ + pcf.setBFieldIsUniform(true); + + /* define pixel inner geometry. + * See CarrierPropagatorTraceExample.java for detailed explanations. + */ + List<SensorRegion> collectors = new ArrayList<SensorRegion>(); + List<SensorRegion> absorbers = new ArrayList<SensorRegion>(); + List<SensorRegion> reflectors = new ArrayList<SensorRegion>(); + SensorRegion ce = new SensorRegion(); + ce.setBox(0.,0.,0.,0.004,0.004,0.002); + collectors.add(ce); + SensorRegion refl = new SensorRegion(); + refl.setBoxLimits(flims[0],flims[1],flims[2],flims[3],flims[5]-0.002,flims[5]+0.002); + reflectors.add(refl); + SensorRegion ab = new SensorRegion(); + ab.setBoxLimits(flims[0],flims[1],flims[2],flims[3],-0.001,0.0001); + absorbers.add(ab); + pcf.setCollectors(collectors); + pcf.setAbsorbers(absorbers); + pcf.setReflectors(reflectors); + + /* Now, when pixel configuration is set, we can create CarrierPropagator for this pixel: */ + CarrierPropagator cpr = new CarrierPropagator(pcf); + + /* setting simulation steps - again, see CarrierPropagatorTraceExample.java for explanations. */ + cpr.setDiffusionSteps(0.0005); + cpr.setDriftSteps(0.0002); +// cpr.setIgnoreDiffusion(true); // for debuging only + + /* In call to makeTable function, we define how many cells in each direction + * our table will be made of. To have cell size 0.5 micron we need 40x40x40 cells for our pixel size. + * Also, the statistics for each point (1000), and number of neighboring pixels (7x7 matrix) are + * passed as arguments. + */ + System.out.println("OK, be partient. It will take long, long time (hours) to finish"); + cpr.makeTable(oname,1000,21,21,21,7,7); + /* + * And finally we will see this message: + */ + System.out.println("Table is recorded !"); + } +}

lcsim/sandbox/NickSinev/Examples

TrackingWithPixSimTest.java added at 1.1

diff -N TrackingWithPixSimTest.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ TrackingWithPixSimTest.java	29 Dec 2010 22:26:56 -0000	1.1
@@ -0,0 +1,255 @@

+import org.lcsim.event.EventHeader;
+import org.lcsim.event.*;
+import org.lcsim.util.Driver;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.recon.vertexing.pixsim.*;
+import org.lcsim.recon.tracking.seedtracker.trackingdrivers.clic_sid.MainTrackingDriver;
+import java.util.*;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+
+public class TrackingWithPixSimTest extends Driver
+{
+   int evnum = 0;
+    AIDA aida = AIDA.defaultInstance();
+    ITree tree = aida.tree();
+
+    IAnalysisFactory af = IAnalysisFactory.create();
+    IHistogramFactory hf = af.createHistogramFactory(tree);
+    IHistogram1D hnhptr = hf.createHistogram1D("Number of hits linked to track",40,0,40);
+    IHistogram1D hnmchpe = hf.createHistogram1D("Number of hits linked to MC particle in event",40,0,40);
+    IHistogram1D hlthr = hf.createHistogram1D("Hit pos. radius of linked to track hits",240,0,1200.);
+    IHistogram1D hltmcr = hf.createHistogram1D("Hit pos. radius of linked to MC particle hits",240,0,1200.);
+    IHistogram2D hlthxy = hf.createHistogram2D("Hit pos. xy of linked to track hits",250,0.,250.,250,0.,250.);
+    IHistogram2D hltmcxy = hf.createHistogram2D("Hit pos. xy of linked to MC particle hits",250,0.,250.,250,0.,250.);
+    IHistogram2D hlthvxy = hf.createHistogram2D("Hit pos. xy of linked to track hits in VTX",160,0.,80.,160,0.,80.);
+    IHistogram2D hltmcvxy = hf.createHistogram2D("Hit pos. xy of linked to MC particle hits in VTX",160,0.,80.,160,0.,80.);
+    IHistogram1D hnarh = hf.createHistogram1D("Number of row hits associated with VTX hit",20,0,20);
+    IHistogram1D hspzr = hf.createHistogram1D("z residual of reconstructed  VTX hit",100,-25.,25.);
+    IHistogram1D hsprpr = hf.createHistogram1D("rphi residual of reconstructed  VTX hit",150,-25.,50.);
+    IHistogram1D hd0res = hf.createHistogram1D("reconstructed track d0 residual, microns",100,-25.,25.);
+    IHistogram1D hz0res = hf.createHistogram1D("reconstructed track z0 residual, microns",100,-25.,25.);
+    IHistogram1D hdptopt = hf.createHistogram1D("reconstructed track deltaPt over Pt",100,-0.02,0.02);
+//    IHistogram2D hrptvspt = hf.createHistogram2D("recon track relative Pt error vs Pt",100,0.,10.,50,-0.02,0.02);
+//    IHistogram2D hrd0vspt = hf.createHistogram2D("recon track d0 (u) vs Pt",100,0.,10.,50,-25,25);
+    IHistogram2D hcovrp = hf.createHistogram2D("rphi cov.matr. error estimate vs res. for  VTX hit",50,0.,25.,50,0.,5.);
+    IHistogram2D hcovz = hf.createHistogram2D("z cov.matr. error estimate vs res. for  VTX hit",50,0.,25.,50,0.,5.);
+    IHistogram1D hcptstm = hf.createHistogram1D("time stamp of central pixel in beam crossing intervals",50,0,200);
+    IHistogram1D hsptstm = hf.createHistogram1D("time stamp of over pixels in cluster in beam crossing intervals",100,0,400);
+
+   public TrackingWithPixSimTest()
+   {
+      // this constructor will build PixilatedSensorManager with sensor option Chronopis 20x20x20 u and use of propagator tables
+      // instead of following track of each charge carrier in the sensor
+        PixilatedSensorManager psm = new PixilatedSensorManager(SensorOption.Chrono20x20x20,true);
+
+      // available options are NOMINAL, SHORT_INT, LONG_INT or CLIC
+        psm.setIlcOption(IlcOption.CLIC);
+      
+      // if option is CLIC we can set time window for accepting clusters. It is set in beam crossings (which is 0.5 ns for CLIC)
+      // and it works only with digiital readout (when chronopix electronics is used). By default for CLIC it is set to 50
+       psm.setTimeWindow(30);
+
+        psm.setNoiseLevel(25.);
+
+       // setting digital readout true or false will assign correct electronics for hit processing
+        psm.setDigitalReadout(true);
+
+      // these settings are relevant only if we are using analog readout. They are ignored if readout is digital
+       psm.setPixelThreshold(3);
+       psm.setClusterThreshold(6);
+
+     // and this  setting only make difference if electronics has analog readout .   
+      psm.setThresholdToNoiseRatio(4.5);
+      // including elrctronic noise hits may increase event processing time by orders of magnitude (depending on noise level and
+      // pixel threshold). So, unless we are interested in the ability of our software to handle huge number of hits, it is
+      // better to keep this option off (set it to false)
+        psm.includeElectronicsNoiseHits(false);
+
+      // we can use pixsim package with any detector, not only clic_sid. To be independent on the fact that detector conditions files
+      // may not have propagator tables for given type of sensors, we can set local option for tables, which will allow us
+      // to read tables not from LCdetectors/detectors/detector_name database, but from local file     
+        psm.useLocalPropagatorTables(false);
+
+      // and the same for resolution tables 
+        psm.useLocalResolutionTables(false);
+
+      // if we defined new sensor configuration and have used CarrierPropagatorMakeTable program to generate propagator tables,
+      // it is easy to build resolution tables, it can be done within pixilated sensor manager
+      // Second argument - number of tracks for each bin of tanlam, tanalpha and sigtonoise. With 1000 it will take about 10 minutes
+      // to build tables with more or less satisfactory precision. 
+
+//        psm.calibrateBarrelResolution(true,1000,"Chronopix_20x20x20_hr_B5_res_bar.dat");  
+//        psm.calibrateEndcapResolution(true,1000,"Chronopix_20x20x20_hr_B5_res_EC.dat","Chronopix_20x20x20_hr_res_EC.dat");
+  
+    add(new MainTrackingDriver(psm));
+   }
+
+   public void process(EventHeader event) 
+   {
+    if(evnum%10 == 0) System.out.println("Processing event: "+evnum);
+    super.process(event);
+     List<List<TrackerHit>> evthits = event.get(TrackerHit.class);
+     List<List<HelicalTrackHit>> helhits = event.get(HelicalTrackHit.class);
+     List<List<Track>> trks = event.get(Track.class);
+     int nmclnk=0;
+     short[] tstmp = new short[50];
+     if(evthits != null)
+     {
+ //     System.out.println("Event has  "+evthits.size()+" TrackerHit collections: ");
+      int narh = 0;
+      double[] mchco=new double[3];
+      for(List<HelicalTrackHit> thits:helhits)
+      {
+        EventHeader.LCMetaData md = event.getMetaData(thits);
+//        System.out.println(md.getName()+" size: "+thits.size());
+        for(HelicalTrackHit tht:thits)
+       {
+        narh=0;
+        List<RawTrackerHit> rhits = tht.getRawHits();
+        boolean lnkMC=false;
+        if(rhits != null)
+        {
+         if(rhits.size() != 0)
+         {
+          for(RawTrackerHit rht:rhits)
+          {
+            short[] adv = rht.getADCValues();
+            if(adv.length == 3)
+            {
+             if(narh<50) tstmp[narh]=adv[2];
+            } 
+            narh++;
+            List<SimTrackerHit> simhits = rht.getSimTrackerHits();
+            if(simhits != null)
+            { 
+             if(simhits.size() != 0)
+             {
+               lnkMC=true;
+               for(SimTrackerHit simhit:simhits)
+              {
+               mchco[0]=(simhit.getPoint())[0];
+               mchco[1]=(simhit.getPoint())[1];
+               mchco[2]=(simhit.getPoint())[2];
+              }
+             }
+            }
+          }
+         }
+        }
+        if(lnkMC)
+        {
+         nmclnk++;
+         double[] hpos=tht.getPosition();
+         double[] covm = tht.getCovMatrix();
+         if(!((covm[0]+covm[2])>0.0000001)) System.out.println("Error in cov. matrix values: "+covm[0]+" , "+covm[2]);
+         double rperr = 1000.*Math.sqrt(covm[0]+covm[2]);
+         double zerr = 1000.*Math.sqrt(covm[5]);
+         double rad=Math.sqrt(hpos[0]*hpos[0]+hpos[1]*hpos[1]);
+         if(rad < 1.) rad=1.;
+         double tl=Math.abs(hpos[2])/rad;
+         if(tl < 0.9) hltmcr.fill(rad);
+         if((hpos[0] > 0.) && (hpos[1] > 0.) && (tl < 0.9) && (rad < 250.)) hltmcxy.fill(hpos[0],hpos[1]);
+         if((hpos[0] > 0.) && (hpos[1] > 0.) && (tl < 0.9) && (rad < 80.)) hltmcvxy.fill(hpos[0],hpos[1]);
+         if((tl<1.0) && (rad < 80.))
+          {
+           hnarh.fill(narh);
+           int mintst=1000;
+           int mintsi=-1;
+           int i;
+           for(i=0; i<narh; i++)
+           {
+             if(tstmp[i] < mintst)
+             {
+               mintst=tstmp[i];
+               mintsi=i;
+             } 
+           }
+           if(mintsi!=-1)
+           {
+            hcptstm.fill(mintst);
+            for(i=0; i<narh; i++)
+             if(i != mintsi) hsptstm.fill(tstmp[i]);
+           }
+           double dz=(hpos[2]-mchco[2])*1000.;
+           double phh = Math.atan2(hpos[1],hpos[0]);
+           double phmc = Math.atan2(mchco[1],mchco[0]);
+           double dx = hpos[0]-mchco[0];
+           double dy= hpos[1]-mchco[1];
+           double drp = Math.sqrt(dx*dx+dy*dy);
+           if(phmc > phh) drp=0.-drp;
+           drp*=1000.;
+           hsprpr.fill(drp);
+           hspzr.fill(dz);
+           hcovrp.fill(Math.abs(drp),rperr);
+           hcovz.fill(Math.abs(dz),zerr);
+          } 
+         }
+       }
+      }
+    }
+    hnmchpe.fill(nmclnk);
+    if(trks != null)
+    {
+//     System.out.println("Seed tracker generated "+trks.size()+" track collections");
+     int cnum=1;
+     for(List<Track> tracks:trks)
+     {
+      EventHeader.LCMetaData md = event.getMetaData(tracks);
+      int tnum=1;
+      for(Track track:tracks)
+      {
+        boolean lnkMC=false;
+      double[] trmom = track.getMomentum();
+      double[] trpars = track.getTrackParameters();
+      double trPt = Math.sqrt(trmom[0]*trmom[0]+trmom[1]*trmom[1]);
+      double trTanL = trmom[2]/trPt;
+      double mcPt=0.;
+      double mcTanL = 0.;
+       List<TrackerHit> trthits=track.getTrackerHits();
+       int nh=trthits.size();
+       hnhptr.fill(nh);
+       tnum++;
+       for(TrackerHit tht:trthits)
+       {
+         double[] hpos=tht.getPosition();
+         double rad=Math.sqrt(hpos[0]*hpos[0]+hpos[1]*hpos[1]);
+         if(rad < 1.) rad=1.;
+         double tl=Math.abs(hpos[2])/rad;
+         if(tl < 0.9) hlthr.fill(rad);
+         if((hpos[0] > 0.) && (hpos[1] > 0.) && (tl < 1.) && (rad < 250.)) hlthxy.fill(hpos[0],hpos[1]);
+         if((hpos[0] > 0.) && (hpos[1] > 0.) && (tl < 1.) && (rad < 80.)) hlthvxy.fill(hpos[0],hpos[1]);
+         List<RawTrackerHit> rhits = tht.getRawHits();
+         for(RawTrackerHit rht:rhits)
+         {
+            List<SimTrackerHit> simhits = rht.getSimTrackerHits();
+            if(simhits != null)
+            { 
+             if(simhits.size() != 0)
+             {
+              MCParticle mcp = simhits.get(0).getMCParticle();
+              lnkMC=true;
+              Hep3Vector mcmom=mcp.getMomentum();
+              mcPt = Math.sqrt(mcmom.x()*mcmom.x()+mcmom.y()*mcmom.y());
+              mcTanL = mcmom.z()/mcPt;
+             }
+           }
+        }
+       }
+       if(lnkMC && (Math.abs(mcTanL) < 1.))
+       {
+         hd0res.fill(trpars[0]*1000.);
+         hz0res.fill(trpars[3]*1000.);
+         hdptopt.fill((trPt-mcPt)/mcPt);
+//         hrptvspt .fill(mcPt,(trPt-mcPt)/mcPt);
+//         hrd0vspt.fill(mcPt,trpars[0]*1000.);
+       }
+
+      }
+     }
+    }
+    evnum++;
+   }
+
+}

\ No newline at end of file

/* And the name for new file, which will be created, and where * transformed map, on regular grid, will be saved: * (I gave it extension txt not to have problems with putting it into CVS.It is really ASCII file, * as all TCAD outputs also). */

Commit in `lcsim/sandbox/NickSinev/Examples` on MAIN
`CarrierPropagatorMakeTable.java`	+142		added 1.1
`TrackingWithPixSimTest.java`	+255		added 1.1
`ChronoPixFieldMove.java`	+7	-6	1.1 -> 1.2
`TCADReadSaveExample.java`	+8	-8	1.1 -> 1.2
	+412	-14