View Javadoc

1   package org.lcsim.recon.cat;
2   
3   import java.util.ArrayList;
4   import java.util.List;
5   
6   import org.lcsim.event.CalorimeterHit;
7   import org.lcsim.event.Cluster;
8   import org.lcsim.event.MCParticle;
9   
10  /**
11   * Simple data format for Minimum ionizing particles in the ECAL.
12   * MipStub is used within the Garfield track finding package for
13   * calorimeter assisted tracking with the SiD
14   *
15   * @see EmcalMipStubs
16   * @see GarfieldTrackFinder
17   *
18   * @author  E. von Toerne
19   * @version $Id: MipStub.java,v 1.3 2012/07/23 10:56:43 grefe Exp $
20   */
21  final public class MipStub implements Cluster {
22    
23    // -- Data :  ----------------------------------------------------------------
24    
25    public double[] base;
26    public double[] dir;
27    public double kappa;  // curvature estimate
28    private double baseError;
29    private double dirXYError;
30    private double dirZError;
31    private double kappaError;
32    private boolean isEndcap = false;
33    public double energy;  // energy estimate
34    private double energyError;
35    public double angleBaseDir;     // angle FROM BASE to DIR in xy
36    //static double maxKappaValue = 1./50.;  // maximum allowed curvature estimate
37    static final double k_cm = 10.;
38    public int debugLevel;
39    private int nHits=0;
40    private int minLayer=-1;
41    private MCParticle[] mcParticles = null;
42    
43    private Cluster cluster = null;
44    
45    // -- Constructors :  --------------------------------------------------------
46    
47    /**
48     * Default constructor
49     */
50    public MipStub() {
51      this.base = new double[]{0,0,0};
52      this.dir = new double[]{0,0,0};
53      this.kappa = 0.;
54      this.energy = -1.;
55      this.angleBaseDir=0.;
56      nHits = 0;
57      normalizeDir();
58      setBaseError();
59      setDirError();
60      setKappaError();
61    }
62    
63    /**
64     * standard constructor
65     * @param p         base, which is also the entrance point into the ECAL
66     * @param d         direction of MIP at entrance point
67     * @param ka        curvature (one over curvature radius)
68     * @param nHit      number of hits in cluster
69     * @param isNdCap   boolean, true if MipStub and its cluster is in Endcap
70     * @param minLay    Minimum layer in ECAL (approx. 0 for MIPS)
71     *
72     */
73    public MipStub(double[] p, double[] d, double ka, int nHit, boolean isNdCap, int minLay) {
74      this.base = new double[]{p[0],p[1],p[2]};
75      this.dir = new double[]{d[0],d[1],d[2]};
76      this.kappa = ka;
77      this.energy = -1.;
78      this.nHits = nHit; //
79      this.isEndcap = isNdCap;
80      this.minLayer = minLay;
81      normalizeDir();
82      setAngleBaseDir();
83      setBaseError();
84      setDirError();
85      setKappaError();
86    }
87    
88    /**
89     * Standard constructor that also takes a reference to the original cluster
90     * (temporary, D.O.)
91     *
92     * @param p         base, which is also the entrance point into the ECAL
93     * @param d         direction of MIP at entrance point
94     * @param ka        curvature (one over curvature radius)
95     * @param nHit      number of hits in cluster
96     * @param isNdCap   boolean, true if MipStub and its cluster is in Endcap
97     * @param minLay    Minimum layer in ECAL (approx. 0 for MIPS)
98     * @param cluster   Cluster from which this MIP stub was created
99     *
100    */
101   public MipStub(double[] p, double[] d, double ka, int nHit, boolean isNdCap, int minLay, Cluster cluster) {
102     this.base = new double[]{p[0],p[1],p[2]};
103     this.dir = new double[]{d[0],d[1],d[2]};
104     this.kappa = ka;
105     this.energy = -1.;
106     this.nHits = nHit; //
107     this.isEndcap = isNdCap;
108     this.minLayer = minLay;
109     normalizeDir();
110     setAngleBaseDir();
111     setBaseError();
112     setDirError();
113     setKappaError();
114     this.cluster = cluster;
115   }
116   
117   // -- Getters :  -------------------------------------------------------------
118   
119   public double kappa(){ return kappa;}
120   public int getnHits(){ return nHits;}
121   public double getPhi(){return Math.atan2(base[1],base[0]);}
122   public double getAngleBaseDir(){ return angleBaseDir;}
123   public double getDirXYError(){ return dirXYError;}
124   public double getDirZError(){ return dirZError;}
125   public double getBaseError(){return baseError;}
126   public double getKappaError(){return kappaError;}
127   public int getMinLayer(){return minLayer;}
128   public boolean isEndcap(){return isEndcap;}
129   public int charge() { return (kappa<0.) ? -1 : 1 ;}
130   public MCParticle[] getMCParticles() {return mcParticles;}
131   
132   public double getEnergyError()
133   {
134       return energyError;
135   }
136   
137   public void setEnergyError(double energyError)
138   {
139       this.energyError = energyError;
140   }
141   
142   // -- Public setters :  ------------------------------------------------------
143   
144   public void setMCParticles(MCParticle[] mcParticles) {this.mcParticles = mcParticles;}
145   public void setKappa(double ka){kappa=ka;}
146   
147   
148   // -- Private helper functions :  ---------------------------------------------
149   
150   /**
151    * normalizes the direction vector. convention here follows {@link GarfieldHelix}
152    * xy-projection of dir is of length one.
153    *
154    * @see             GarfieldHelix
155    */
156   private void normalizeDir(){
157     // normalization: xy direction has length one.
158     double d=Math.sqrt(dir[0]*dir[0]+dir[1]*dir[1]);
159     if (d>1.E-20){
160       dir[0]=dir[0]/d;
161       dir[1]=dir[1]/d;
162       dir[2]=dir[2]/d;
163     } else {
164       dir[0]=1;
165       dir[1]=0.;
166       dir[2]=dir[2]*1.E10;
167     }
168   }
169   
170   /**
171    * calculates the angle from base to the dir vector in the xy plane.
172    */
173   private void setAngleBaseDir() {
174     double lBase = Math.sqrt(base[0]*base[0]+base[1]*base[1]);
175     if (lBase>1.E-9){
176       double c=(base[0]*dir[0]+base[1]*dir[1])/(lBase);  // please note that dir-xy is of norm one
177       if (Math.abs(c)>1.){ angleBaseDir = 0.;} else {angleBaseDir = Math.acos(c);}
178       if (base[0]*dir[1]-base[1]*dir[0]<0.) angleBaseDir=-angleBaseDir;
179     } else angleBaseDir = 0.;
180   }
181   
182   /**
183    * calculates a rough error estimate for direction vector in XY plane
184    * based on <code>nHits</code> only
185    */
186   private void setDirError(){
187     if (nHits>=10) dirXYError=20.;
188     else if (nHits >=5) dirXYError=25.;
189     else dirXYError= 45.;
190     dirXYError = Math.sqrt(dirXYError*dirXYError+900.*angleBaseDir*angleBaseDir);
191     dirZError = dirXYError / 5. * k_cm; // set dirZError
192     // conversion from degrees into rads
193     dirXYError = dirXYError*Math.PI/180.;
194     if (isEndcap()){
195       dirZError = dirZError * 10.;
196       dirXYError = dirXYError * 10.;
197     }
198   }
199   
200   private void setBaseError(){ baseError = 1.0 * k_cm;}
201   
202   private void setKappaError(){
203     //kappaError = 0.001/k_cm;
204     //kappaError = Math.sqrt(kappaError*kappaError+0.3*0.3*kappa*kappa);
205     if (nHits>=15) kappaError=0.03;
206     else if (nHits>=10) kappaError=0.035;
207     else if (nHits>=5) kappaError=0.05;
208     else kappaError=0.1;
209     if (isEndcap()) kappaError = kappaError * 10.;
210   }
211   
212   // -- Debugging :  -----------------------------------------------------------
213   
214   // plotting, testing and debugging
215   void debug(){
216     //if (debugLevel <1) return;
217     System.out.println("MipStub base="+this.base[0]+" "+this.base[1]+" "+this.base[2]);
218     System.out.println("MipStub dir="+this.dir[0]+" "+this.dir[1]+" "+this.dir[2]+" kappa="+this.kappa+" angle="+angleBaseDir);
219   }
220   
221   /**
222    * prepares graphic presentation of a mipstub by providing a list of positions
223    * this plot routine uses quite a bit of memory.
224    *
225    */
226   public List plot(){
227     // creates a list of 3-D point that represent the MipStub and that that can be used in plotting of histograms
228     List lst = new ArrayList();
229     double s = 2.*k_cm; // step size for the plot
230     for (int j=0;j<=5;j++){
231       double[] pix=new double[]{0,0,0};
232       pix[0]=base[0]+s*dir[0]*j;
233       pix[1]=base[1]+s*dir[1]*j;
234       pix[2]=base[2]+s*dir[2]*j;
235       lst.add(pix);
236     }
237     GarfieldHelix hh = new GarfieldHelix(base,dir,kappa);
238     for (int j=0;j<=10;j++){
239       double[] pix=new double[]{0,0,0};
240       hh.setPointOnHelix(-8.*k_cm*j);
241       pix[0]=hh.getPointOnHelix(0);
242       pix[1]=hh.getPointOnHelix(1);
243       pix[2]=hh.getPointOnHelix(2);
244       lst.add(pix);
245     }
246     return lst;
247   }
248   
249   // -- Implementing org.lcsim.event.Cluster :  --------------------------------
250   // -- (dummy implementations for some methods - this should be fixed later) --
251 
252   public int getType() {return 0;}  // FixMe
253 
254   public double[] getSubdetectorEnergies() {return new double[] {getEnergy()};}
255 
256   public double[] getShape() {return new double[6];}  // Meaning is undefined in LCIO
257 
258   // TODO return meaningful position error
259   public double[] getPositionError() {return new double[] {baseError, baseError, baseError, baseError, baseError, baseError};}
260 
261   public double[] getPosition() {return base;}
262 
263   public double getITheta() {return 0.;}
264 
265   public double getIPhi() {return 0.;}
266 
267   public double[] getHitContributions() {
268     if (cluster == null) {
269       return new double[] {0.};
270     } else {
271       return cluster.getHitContributions();
272     }
273   }
274 
275   public double getEnergy() {return energy;}
276 
277   public double[] getDirectionError() {return new double[] {dirXYError, dirXYError, dirZError};}
278 
279   public List<Cluster> getClusters() {
280     List<Cluster> list = new ArrayList<Cluster>(1);
281     if (cluster != null) {
282       list.add(cluster);
283     }
284     return list;
285   }
286 
287   public List<CalorimeterHit> getCalorimeterHits() {
288     if (cluster == null) {
289       return new ArrayList<CalorimeterHit>(1);
290     } else {
291       return cluster.getCalorimeterHits();
292     }
293   }
294 
295   public int getSize() {return nHits;}
296 
297   public int getParticleId() {
298       return 0;
299   }
300   
301 }
302