View Javadoc

1   /*
2    * HelixFitter.java
3    *
4    * Created on January 22, 2008, 9:25 AM
5    *
6    */
7   
8   package org.lcsim.recon.tracking.seedtracker;
9   
10  import java.util.List;
11  import java.util.Map;
12  
13  import org.lcsim.fit.circle.CircleFit;
14  import org.lcsim.fit.helicaltrack.*;
15  import org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus;
16  import org.lcsim.fit.line.SlopeInterceptLineFit;
17  import org.lcsim.recon.tracking.seedtracker.diagnostic.ISeedTrackerDiagnostics;
18  import org.lcsim.fit.zsegment.ZSegmentFit;
19  
20  /**
21   *
22   * @author Richard Partridge
23   * @version 1.0
24   */
25  public class HelixFitter {
26      private HelicalTrackFitter _fitter = new HelicalTrackFitter();
27      protected MultipleScattering _scattering;
28      private HelicalTrackFit _helix;
29      private MaterialManager _materialmanager;
30      private ConstrainHelix _constrain;
31      private double _bfield = 0.;
32      private CircleFit _circlefit;
33      private SlopeInterceptLineFit _linefit;
34      private ZSegmentFit _zsegmentfit;
35      private FitStatus _status;
36      private ISeedTrackerDiagnostics _diag = null;
37      TrackCheck _trackCheck; // set by SeedTracker
38      private boolean _debug = false;
39  
40      /**
41       * Creates a new instance of HelixFitter
42       */
43      public HelixFitter(MaterialManager materialmanager) {
44          _materialmanager = materialmanager;
45          _scattering = new MultipleScattering(_materialmanager);
46          _constrain = new ConstrainHelix();
47      }
48      
49      public boolean FitCandidate(SeedCandidate seed, SeedStrategy strategy) {
50          if(_debug) System.out.printf("%s: FitCandidate() with stategy: \"%s\"\n",this.getClass().getSimpleName(),strategy.getName());
51          
52          //  Initialize fit results to null objects
53          _helix = null;
54          
55          //  Check that we have set the magnetic field
56          if (_bfield == 0.) throw new RuntimeException("B Field must be set before calling the Helix fitter");
57          
58          //  Set the tolerance in the fitter
59          _fitter.setTolerance(Math.sqrt(strategy.getMaxChisq()));
60          
61          //  Retrieve list of hits to be fit
62          List<HelicalTrackHit> hitlist = seed.getHits();
63          
64          if(_debug) {
65              System.out.println(this.getClass().getSimpleName()+": hitlist size " + hitlist.size() + ":");
66              double z_prev = -99999999.9;
67              for (HelicalTrackHit hit : hitlist) {
68                  System.out.printf("%s: (%.2f,%.2f,%.2f) corrected  %s\n",
69                                      this.getClass().getSimpleName(),hit.getPosition()[0],hit.getPosition()[1],hit.getPosition()[2],
70                                      hit.getCorrectedPosition().toString());
71                  if(z_prev<-99999999.0) z_prev=hit.getPosition()[2];
72                  else {
73                      if(Math.signum(z_prev)!=Math.signum(hit.getPosition()[2])) {
74                          System.out.printf("%s: topBotHits in event\n",this.getClass().getSimpleName());
75                      }
76                      z_prev = hit.getPosition()[2];
77                  }
78              }
79          }
80          //  Retrieve the old helix
81          HelicalTrackFit oldhelix = seed.getHelix();
82          
83          //  If this is the candidate's first helix fit, first do a fit without MS errors
84          if (oldhelix == null) {
85       
86              if(_debug) 
87                  System.out.println(this.getClass().getSimpleName()+": no old helix do the fit without MS scattering map" );
88              
89              //  Reset the stereo hit positions to their nominal value
90              for (HelicalTrackHit hit : hitlist) {
91                  if (hit instanceof HelicalTrackCross) ((HelicalTrackCross) hit).resetTrackDirection();
92              }
93              //  Do the fit
94              _status = _fitter.fit(hitlist);
95              SaveFit();
96              
97              //  Check for unrecoverable fit errors and call appropriate diagnostic
98              if (_status != FitStatus.Success) {
99                  if(_diag!=null) _diag.fireHelixFitFailed(seed, _status, true);
100                 return false;
101             }
102             
103             //  Retrieve the helix parameters from this fit and save them in the seed
104             oldhelix = _fitter.getFit();
105             seed.setHelix(oldhelix);
106 
107             if(_debug) System.out.printf("%s: fit succeeded, will be used as seed, with chi2=%.3f and helix:\n%s \n",this.getClass().getSimpleName(),oldhelix.chisqtot(),oldhelix.toString());
108 
109             //  Calculate the multiple scattering angles for this helix
110             try {
111                seed.setScatterAngles(_scattering.FindScatters(oldhelix)); 
112             } catch (Exception e) {
113                System.err.println(e);
114                if(_debug)
115                {
116                    e.printStackTrace();
117                }
118                return false;
119             }
120             
121             if(_debug) {
122                 System.out.printf("%s: after calculating the MS map it has %d size:\n",this.getClass().getSimpleName(),seed.getMSMap().size());
123                 for(Map.Entry<HelicalTrackHit, MultipleScatter> ms : seed.getMSMap().entrySet()) {
124                     System.out.printf("%s: Hit at layer %d and position %s has MS drpdhi=%f and dz=%f\n",
125                                     this.getClass().getSimpleName(),ms.getKey().Layer(),ms.getKey().getCorrectedPosition().toString(),ms.getValue().drphi(),ms.getValue().dz());
126                 }
127             }
128             
129         }
130         
131         if(_debug) 
132             System.out.printf("%s: update the stereo hit positions with the old helix:\n%s \n",this.getClass().getSimpleName(),oldhelix.toString());
133         
134         //  Update the stereo hit positions and covariance matrices
135         CorrectStereoHits(hitlist, oldhelix);
136         
137         //  Do a helix fit including MS errors
138         if(_debug) {
139             System.out.printf("%s: do the helix fit including MS map this time: \n",this.getClass().getSimpleName());
140             for(Map.Entry<HelicalTrackHit, MultipleScatter> ms : seed.getMSMap().entrySet()) {
141                     System.out.printf("%s: Hit at layer %d and position %s has MS drpdhi=%f and dz=%f\n",
142                                     this.getClass().getSimpleName(),ms.getKey().Layer(),ms.getKey().getCorrectedPosition().toString(),ms.getValue().drphi(),ms.getValue().dz());
143             }
144         }
145         
146         _status = _fitter.fit(hitlist, seed.getMSMap(), oldhelix);
147         SaveFit();
148         
149         //  Check for unrecoverable fit errors and call appropriate diagnostic
150         if (_status != FitStatus.Success) {
151             if(_diag!=null) _diag.fireHelixFitFailed(seed, _status, false);
152             return false;
153         }
154         
155         //  Retrieve and save the new helix fit
156         _helix = _fitter.getFit();
157         seed.setHelix(_helix);
158          if ((_trackCheck != null) && (! _trackCheck.checkSeed(seed))) {
159             return false;
160          }
161 
162         //  Set the non-holonomic constraint chi square
163         _constrain.setConstraintChisq(strategy, _helix, hitlist);
164         
165         //  If the total chi square is below the cut, we have a successful fit
166         boolean success = _helix.chisqtot() <= strategy.getMaxChisq();
167         if (!success)
168             if (_diag != null) _diag.fireFailedChisqCut(seed);
169 
170         //  If fit was successful, set the new multiple scattering angles
171         if (success) {
172                 
173             seed.setScatterAngles(_scattering.FindScatters(_helix));
174             if(_debug) {
175                 System.out.printf("%s: this fit was successful, chi2=%f with resulting helix paramaters:\n%s\n",this.getClass().getSimpleName(),_helix.chisqtot(),_helix.toString());
176                 System.out.printf("%s: updated MS map before returning from fitCandidate():\n",this.getClass().getSimpleName());       
177                 for(Map.Entry<HelicalTrackHit, MultipleScatter> ms : seed.getMSMap().entrySet()) {
178                     System.out.printf("%s: Hit at layer %d and position %s has MS drpdhi=%f and dz=%f\n",
179                                     this.getClass().getSimpleName(),ms.getKey().Layer(),ms.getKey().getCorrectedPosition().toString(),ms.getValue().drphi(),ms.getValue().dz());
180                 }
181             }
182         } else {
183             if(_debug) 
184                 System.out.printf("%s: this fit with chi2=%f failed!\n",this.getClass().getSimpleName(),_helix.chisqtot());
185         }
186         
187         return success;
188     }
189     
190     public void setDiagnostics(ISeedTrackerDiagnostics d) {
191         _diag = d;
192         return;
193     }
194     
195     public HelicalTrackFit getHelix() {
196         return _helix;
197     }
198     
199     public FitStatus getFitStatus() {
200         return _status;
201     }
202     
203     public CircleFit getCircleFit() {
204         return _circlefit;
205     }
206     
207     public SlopeInterceptLineFit getLineFit() {
208         return _linefit;
209     }
210     
211     public ZSegmentFit getZSegmentFit() {
212         return _zsegmentfit;
213     }
214     
215     public void setBField(double bfield) {
216         _bfield = bfield;
217         _scattering.setBField(_bfield);
218         _constrain.setBField(_bfield);
219         return;
220     }
221 
222      public void setReferencePoint(double x,double y) {
223         _fitter.setReferencePoint(x, y);
224         return;
225     }
226 
227     private void SaveFit() {
228         
229         //  Default to no fit results when circle fit fails
230         _circlefit = null;
231         _linefit = null;
232         _zsegmentfit = null;
233         
234         //  Save the circle fit results if they exist
235         if (_status == FitStatus.CircleFitFailed) return;
236         _circlefit = _fitter.getCircleFit();
237         
238         //  If we have a circle fit, try to save the line fit / zsegment fit results
239         if (_status == FitStatus.InconsistentSeed) return;
240         if (_status == FitStatus.LineFitFailed) return;
241         if (_status == FitStatus.ZSegmentFitFailed) return;
242         _linefit = _fitter.getLineFit();
243         _zsegmentfit = _fitter.getZSegmentFit();
244         
245         return;
246     }
247     
248     private void CorrectStereoHits(List<HelicalTrackHit> hitlist, HelicalTrackFit helix) {
249         
250         //  Get the PathMap - used to find the track direction at the hit location
251         Map<HelicalTrackHit, Double> pathmap = helix.PathMap();
252         
253         //  Loop over the hits and look for stereo hits
254         for (HelicalTrackHit hit : hitlist) {
255             if (hit instanceof HelicalTrackCross) {
256                 
257                 //  Found a stereo hit - calculate the track direction and pass it to the hit
258                 ((HelicalTrackCross) hit).setTrackDirection(helix);
259             }
260         }
261         return;
262     }
263 
264     public void setDebug(boolean debug) {
265         this._debug = debug;
266     }
267  }