Mattione sim-recon Code Documentation
From GlueXWiki
Contents
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 Candidate Reconstruction
DTrackCandidate_factory::evnt()
DTrackCandidate_factory_CDC::evnt()
- Call DTrackCandidate_factory_CDC::GetCDCHits()
- call JEventLoop::Get() to get the "raw" hits (DCDCTrackHit) (already have the wire associated with them).
- loop over DCDCTrackHit's and create DCDCTrkHit objects (one per hit wire), initially labeling each hit as noise
- sort the DCDCTrkHit objects into vectors by superlayer number, and sort the vectors by decreasing values of ring# (if same ring, increasing # of straw)
- remove noise hits: hits that are not within some minimum distance (01/23/2012 : 12) of any other hit (using the wire origins as the hit positions along the wire)
- Call DTrackCandidate_factory_CDC::FindSeeds() with the DCDCTrkHit vectors for the axial superlayers (1, 3, 5). Fills objects of vector<DCDCSeed>.
- look for subseeds: all hits are subseeds, but group adjacent hits within a ring into a single subseed. Within each ring there may be several subseeds.
- link subseeds together: loop over subseeds within each ring. If a subseed hasn't been linked yet, call LinkSubSeeds()
- in the next ring, look for subseeds that are close to the last linked subseed found (the input if none found yet). If one is found, link them together and recursively call LinkSubSeeds() with the next ring.
- If no link is found in the next ring, then save the list of linked hits as a DCDCSeed and return into the subseed loop of the previous ring
- Note: if two subseeds in a ring are close enough to be linked, then during the chain of recursive returns a new branch will be started back at its location, re-starting the recursive LinkSubSeeds() calls, and resulting in a new DCDCSeed.
- Call DTrackCandidate_factory_CDC::LinkSeeds() twice: first with the DCDCSeeds from superlayers 3 & 5, then with the DCDCSeeds formed from that combo & from superlayer 1
- Call DTrackCandidate_factory_CDC::PickupUnmatched()
- Check to add lone hits in SL3 to seeds with only SL1 hits
- NOT FINISHED
- Call DTrackCandidate_factory_CDC::DropIncompleteSeeds()
- Drop seeds containing hits from only a single axial superlayer besides SL1
- NOT FINISHED
- Call DTrackCandidate_factory_CDC::FitCircle() to each DCDCSeed
- Call DTrackCandidate_factory_CDC::FilterCloneSeeds()
- Filter out duplicates of seeds by clearing their "valid" flags
- NOT FINISHED
- Loop over DCDCSeeds:
- Call DTrackCandidate_factory_CDC::AddStereoHits() separately for each stereo superlayer (2 then 4)
- Call DTrackCandidate_factory_CDC::FindThetaZRegression() to find tan(lambda) & z_vertex
- NOT FINISHED
- If linear regression fit didn't work, call DTrackCandidate_factory_CDC::FindThetaZ() to find tan(lambda) & z_vertex via the histogramming method
- NOT FINISHED
- If fit doesn't work, reject seed ( continue to next DCDCSeed in loop)
- Call DTrackCandidate_factory_CDC::GetPositionAndMomentum() to extract the track momentum and position from the DCDCSeed.
- Create a new DTrackCandidate object, and set it up with the momentum, position, charge, and DCDCTrackHit objects
DTrackCandidate_factory_CDC::LinkSeeds()
- nested loop over the DCDCSeeds: loop over the first superlayer, and inside of it loop over the other superlayer
- select the wires of the first and last subseeds of the current DCDCSeeds in the loop
- compare the rings of these 4 wires: select the ring in each DCDCSeed that is closest to the other DCDCSeed.
- compare the phi of the wires of these hits: if they are within pi/6 (30 degrees): merge the hits from both DCDCSeeds into a new DCDCSeed (and save it)
- if the hits were merged AND if the length of the chord (of the circle of the innermost ring of the two subseeds) between the two subseeds is less than 12 cm:
- mark each subseed as linked in the input lists as long as the other seed doesn't have too many hits (12 for the superlayer 3 & 5 comparison, 24 for the comparison against superlayer 1)
- else leave as unlinked so that it has the option to be used to create an additional track candidate at the next level
- for each of the input DCDCSeeds that were not marked as linked, create a new DCDCSeed for each (and save it)
- loop over all of the new DCDCSeeds and set the average drift time (from the drift times of the individual hits).
- if the drift time is > 1000, it is labeled as an "OUT_OF_TIME" hit and is excluded from the average
DTrackCandidate_factory_CDC::FitCircle()
- Add the wire origin position of each hit in the seed that is not "OUT_OF_TIME" to hit fit list.
- Add any fdchits (DFDCPseudo) to the hit fit list if they exist
- Fit the hits to a circle using DHelicalFit::FitCircleRiemann(). If it fails, try using a regression fit DHelicalFit::FitCircle(). If it fails, try a straight line fit with DHelicalFit::FitCircleStraightTrack()
- In the DHelicalFit::FitCircleRiemann() function, a fake hit at the beamline is added to the track for the fit (with sigma of 0.1)
- In the DHelicalFit::FitCircle() function, it assumes that the track originates from the beamline.
- Count how many hits on the DCDCSeed are within 4 cm of the circle fit
- if fewer than 2 hits passed, then reject the fit
- If fewer than half of the hits passed, then reject the fit
- if a fit fails or is rejected, return false (then the DCDCSeed is discarded ("valid" flag set to false))
DTrackCandidate_factory_CDC::AddStereoHits()
-
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 (theoretically should enter sigma as sqrt (because factor is sigma^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 residual for both anode and cathode:
- anode is difference btw track doca to hit and measured distance "use cell size & drift times (for TimeBased) to find "measured" distance"
- cathode is difference btw the distances along the strip determined from the reference trajectory and from the hit
- calc sigma (on residual), chisq, & FOM: sigma is relatively arbitrary (NOT FINISHED!!)
- sigma (both anode & cathode, but different parameters for each)= scale factor * multiple scattering factor * momentum factor
- scale factor is fit-type based, m-s factor is e^(linear function of: log of "running sum of: (MULS angle theta_0 squared) * (distance along the trajectory)^2"), momentum factor is linear (1 for now!)
- if FOM > 5% 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
