Mattione sim-recon Code Documentation

Plugin-Independent Execution

hd_root: main()

- instantiates MyProcessor and DApplication
  - MyProcessor inherits from JEventProcessor
    - JEventProcessor basically empty/virtual, constructor only initializes a few conrol variables during instantiation
    - MyProcessor constructor basically empty
  - DApplication inherits from JApplication
    - JApplication constructor sets up user signal instructions, mutexes, parses the command line, sets up the parameters in JParameterManager, etc.
    - DApplication constructor instantiates the DEventSourceHDDMGenerator and DFactoryGenerator, and registers them with JApplication
- calls DApplication.Run(JEventProcessor) with the instantiated MyProcessor as the argument (DApplication.Run(JEventProcessor) is JApplication.Run(JEventProcessor))
  - Calls DApplication.Init()
    - Calls JApplication.Init()
      - Attaches plugins (JEventProcessors), adds auto-activated factories (from the JParameterManager), calls JEventProcessor.init() for each processor (MyProcessor.init() creates the output root file) 
    - Checks to see if should use SSE instructions
  - Launches threads (each thread calls JApplication.LaunchThread())
  - Sends the main thread to sleep while the threads execute, wake up occassionally to check status, exits when done
- Exits

JApplication.LaunchThread()

- Create JEventLoop with the instantiated JApplication object as the argument
  - JEventLoop constructor: 
    - registers the JEventLoop with the JApplication object and its member JEvent object "event"
- call JEventLoop.Loop() on the JEventLoop object
  - calls JEventLoop.OneEvent() in an infinite loop until that function tells it to quit the loop (that function loops over each event)
    - calls JEventLoop.Initialize() if not initialized: gets event processors and auto-activated factories from JApplication (which were listed in JParameterManager)
    - calls JEventLoop.ClearFactories(), which calls the .Reset() function of each factory
    - calls JApplication.NextEvent(JEvent), which grabs the next event from the event buffer
    - loops over JEventProcessors (for plugins, functionality is plugin-specific), calling these virtual functions:
      - JEventProcessors.brun() if it's a new run number:
        - if hd_root, calls MyProcessor.brun(): figures out which factories are needed, sets up their information (factory_info_t) in MyProcessor.fac_info
      - JEventProcessors.erun() if the run number changes:
        - if hd_root, calls MyProcessor.erun(): does nothing
      - JEventProcessors.evnt() for each event:
        - if hd_root, calls MyProcessor.evnt(): loops over all factories listed in MyProcessor.fac_info, grabs them from the JEventLoop, for each calls factory->GetNrows()
    - calls JEvent.FreeEvent() to free up the memory from this event
- Exits

JEventLoop.Get(vector<const T*>)

• Call to retrieve data objects of type T with an associated factory (e.g. DChargedTrack, DPhysicsEvent)

- The JEventLoop.Get(vector<const T*>) function retrieves or generates the requested objects in this order:
  - If already generated, returns them.  Else if present in input data file, retrieves them.  Else generates them via the appropriate JFactory.
- These objects are only generated once, and the generating factories retain ownership (only const pointers are returned)
- To check if they were already generated, JEventLoop.Get(vector<const T*>) calls JEventLoop::GetFromFactory()
  - It grabs the associated JFactory.  Note that all of the factories inherit from JFactory, which inherits from JFactory_base, which inherits from JEventProcessor
  - If its JFactory.evnt() method has already been called, then the JFactory.CopyFrom(vector<const T*>) method is called to just copy the already-generated information into the input vector.
- To check if they are located in the source, JFactory.GetCheckSourceFirst() is called.  If so, then JEventLoop.GetFromSource() is called to retrieve the input information. 
- To generate the data, the JFactory.Get(vector<const T*>) method is called.
  - If JFactory.evnt() has already been called, then it calls JFactory.CopyFrom(vector<const T*>) to copy the data into the input vector and returns.  Else the function continues.  
  - It calls JFactory.init() if it has not been called already.  
  - It then calls JFactory.erun() and JFactory.brun() if appropriate.  
  - Then JFactory.evnt() is called to generate the data, then it calls JFactory.CopyFrom(vector<const T*>) to copy the data into the input vector and returns.  

Track Reconstruction

DTrackFitter::FindHitsAndFitTrack()

- call CorrectForELoss (WireBased fit only): assumes momentum from DTrackCandidate is at middle of track (avg), swims backwards to target center & increases momentum based on eloss
- swim reference trajectory
- get DTrackHitSelector_ALT1, all DCDCTrackHit objects, and all DFDCPseudo objects
- call DTrackHitSelector::GetAllHits() with reference trajectory, DTrackHitSelector::fit_type_t (kHelical for kWireBased fit, or kWireBased for kTimeBased fit), DCDCTrackHit's, DFDCPseudo's
  - calls DTrackHitSelector::GetCDCHits
    - calls DTrackHitSelector_ALT1::GetCDCHits (with reference trajectory, DTrackHitSelector::fit_type_t, all DCDCTrackHit objects), sorts hits by ring, then stores them in the fitter
  - calls DTrackHitSelector::GetFDCHits 
    - calls DTrackHitSelector_ALT1::GetFDCHits (with reference trajectory, DTrackHitSelector::fit_type_t, all DFDCPseudo objects), sorts hits by z, then stores them in the fitter
- call FitTrack with position, momentum, q, mass, start time: sets DTrackFitter.input_params (a DKinematicData object) with these info, then calls DTrackFitterKalmanSIMD::FitTrack(void)

DTrackFitterKalmanSIMD::FitTrack(void)

- calls ResetKalmanSIMD to reset data
- push back cdc hits into my_cdchits, fdc hits into my_fdchits
  - for fdc: hit->dE=1e6*fdchit->dE;, hit->xres=hit->yres=1000.;
- start time (mT0) set from input, start time variance (mVarT0) fixed to 0.09
- calls SetSeed to set track parameters from input pos, mom, and charge
  - Forward: x_, y_, z_, tx_ (px/pz), ty_ (py/pz), q_over_p_
  - Central: phi_ (p.phi()), tanl_ (l = pi/2 - p.theta()), q_over_pt_  
- set the mass (m_ratio = ELECTRON_MASS/MASS)
- call DTrackFitterKalmanSIMD::KalmanLoop() to do fit
- set fit results into DTrackFitter base class input_params (for subsequent output)
  - charge, momentum, t0 (from cdc or fdc), position, mass
  - also cov matrix: from "fcov" for forward params, else from "cov"
- sets cdchits_used_in_fit & fdchits_used_in_fit

DTrackFitterKalmanSIMD::KalmanForward()

- loop through track trajectory steps (stored in forward_traj) from the last step (largest z) to the first step
  - transform the state vector to the new step
  - project the error matrix from the previous step to the current step, and add the uncertainties due to eloss and multiple scattering
  - if any fdc hits (should always be true!!):
    - if no fdc hits associated with this step on the trajectory, continue
    - call fdc_y_variance() to calculate the variance on the track position in the coordinate along the wire (hard-coded, only hit energy dependent)
    - compute the residual: difference between the projected hit and the measured hit: 'Mdiff'. 
      - if timebased tracking, use drift distance for Mdiff(0) too, not just doca. Also call fdc_drift_variance() to calculate the variance due to the drift (hard-coded func of drift-time)  
    - compute the transformation matrix 'H' to transform between the hit position coordinates (u, v) to the track parameter coordinates (x, y)
    - if more than one fdc hit at this track step:
      - calculate the variance on the first fdc hit: transform covariance matrix 'C' using 'H' to u, v coordinates, add V
      - calculate the chi-squared of the hit residual
      - calculate the probability of the first fdc hit using the residual chi-squared
      - if chisq is too large, ignore the hit (outlier) (but still look at other fdc hits)
      - loop over other fdc hits, doing the same as the first fdc hit (calc residual, variance, chisq, probability; if chisq too high, ignore it)
      - loop over the hits whose chi-squared wasn't too large, compute the new state vector 'S' and covariance matrix 'C' using the hit probabilities to determine the weights (also set in fdc_updates)
    - else
      - calculate the variance on the hit: transform covariance matrix 'C' using 'H' to u, v coordinates, add V
      - calculate the chi-squared of the hit residual
      - if the chi-squared of the hit residual is too large, ignore the hit (outlier) (continue loop through track trajectory steps)
      - update the track state vector and covariance matrix
      - add the chi-squared of the hit residual to the total chi-squared (starts at zero) (track chi-squared = total chi-squared of the hit residuals)
  - else if any cdc hits:
    - NOT FINISHED!!!! (Can this EVER be called?? The only place that KalmanForward() is called is in KalmanLoop(), and only if there ARE fdc hits)
- set x_, y_, and z_ to those of the final state vector (lowest z)
- multiply chisq by the 'anneal_factor' which is 1 (but the fdc y & drift variances are already multiplied by this!!)

Hit Selection

DTrackHitSelector_ALT1::GetCDCHits()

- loop over DCDCTrackHit objects:
  - calc residual: diff btw "use reference trajectory to find DOCA of track to hit" and "use cell size & drift times (for TimeBased) to find "measured" distance"
  - calc sigma (on residual), chisq, & FOM: sigma is relatively arbitrary
     - sigma = scale factor * multiple scattering factor * momentum factor
        - scale factor is fit-type based, m-s factor is a linear function of the distance along the trajectory, momentum factor is 1 + factor*e^(-p^2)
  - if FOM > 5% then keep the DCDCTrackHit
- sort kept DCDCTrackHit objects by ring & hit probability (reject hits in the same ring that are far away)
- output saved DCDCTrackHit objects

DTrackHitSelector_ALT1::GetFDCHits()

- loop over DFDCPseudo objects:
   - calc track doca to hit
   - calc sigma (on track doca to hit), chisq, & FOM: sigma is relatively arbitrary
   - if FOM > 1% then keep the DCDCTrackHit
- sort kept DFDCPseudo objects by ring & hit probability (reject hits in the same ring that are far away)
- output saved DFDCPseudo objects

DTrackHitSelector_ALT2::GetCDCHits()

- loop over DCDCTrackHit objects:
  - calc residual: diff btw "use reference trajectory to find DOCA of track to hit" and "use cell size & drift times (for TimeBased) to find "measured" distance"
  - calc sigma (on residual), chisq, & FOM: sigma is relatively arbitrary
     - sigma = scale factor * multiple scattering factor * momentum factor
        - scale factor is fit-type based, m-s factor is a linear function of "running sum of: (MULS angle theta_0 squared) * (distance along the trajectory)^2", momentum factor is linear function of momentum
        - "MULS angle theta_0 squared" = multiple scattering theta_rms ^ 2
  - if FOM > 5% then keep the DCDCTrackHit
- sort kept DCDCTrackHit objects by ring & hit probability (reject hits in the same ring that are far away)
- output saved DCDCTrackHit objects

DTrackHitSelector_ALT2::GetFDCHits()

- loop over DFDCPseudo objects:
   - calc track doca to hit
   - calc sigma (on track doca to hit), chisq, & FOM: sigma is relatively arbitrary
   - if FOM > 1% then keep the DCDCTrackHit
- sort kept DFDCPseudo objects by ring & hit probability (reject hits in the same ring that are far away)
- output saved DFDCPseudo objects