class to handle first step of the shower reconstruction in ORCA: it reconstructs the shower vertex, intended as the shower brightest point, as the barycenter of the hits More...

#include <JShowerPrefit.hh>

Inheritance diagram for JRECONSTRUCTION::JShowerPrefit:

Classes
struct	match_t

Public Member Functions
	JShowerPrefit (const JShowerPrefitParameters_t &parameters, const JModuleRouter &router, const int debug=0)
	Parameterized constructor. More...

JEvt	operator() (const KM3NETDAQ::JDAQEvent &event) const
	Declaration of the operator that performs the reconstruction. More...

void	reset ()
	Reset fit parameters. More...

bool	equals (const JShowerPrefitParameters_t &parameters) const
	Equality. More...

	ClassDef (JShowerPrefitParameters_t, 2)

Public Attributes
const JModuleRouter &	router

int	factoryLimit

size_t	NMaxHits

double	sigma_ns

int	numberOfOutliers

double	TMaxLocal_ns

double	ctMin

double	roadWidth_m

Detailed Description

class to handle first step of the shower reconstruction in ORCA: it reconstructs the shower vertex, intended as the shower brightest point, as the barycenter of the hits

Definition at line 50 of file JShowerPrefit.hh.

Constructor & Destructor Documentation

JRECONSTRUCTION::JShowerPrefit::JShowerPrefit	(	const JShowerPrefitParameters_t &	parameters,
		const JModuleRouter &	router,
		const int	debug = `0`
	)

inline

Parameterized constructor.

Parameters

parameters	struct that holds default-optimized parameters for the reconstruction.
router	module router, this is built via detector file.
debug	debug

Definition at line 63 of file JShowerPrefit.hh.

                                                               :
       JShowerPrefitParameters_t(parameters),
       router(router)
     {}

Member Function Documentation

JEvt JRECONSTRUCTION::JShowerPrefit::operator() ( const KM3NETDAQ::JDAQEvent & event ) const

inline

Declaration of the operator that performs the reconstruction.

Parameters

event which is a JDAQEvent

Definition at line 76 of file JShowerPrefit.hh.

     {
 
       using namespace std;
       using namespace JPP;
 
       const double STANDARD_DEVIATIONS       = 3.0;
 
       typedef JHitR1           hit_type;
       typedef vector<hit_type> buffer_type;
 
       typedef JEstimator<JPoint4D> JEstimator_t;
 
       const JBuildL0<hit_type> buildL0;
       const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));
       const JMatch3G<hit_type> match3G(roadWidth_m, TMaxLocal_ns);  // causality relation for showers 
       const JHitL1Comparator compare;
   
       //define empty vectors to be filled with JHitL0 and JHitL1
       buffer_type dataL0;
       buffer_type dataL1;
   
       JEvt out;   
 
       buffer_type buffer;
 
       buildL0(JDAQTimeslice(event, true),  router, back_inserter(buffer));  // true = snapshot
       buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1));  // false = triggered
       
       copy(dataL1.begin(), dataL1.end(), back_inserter(dataL0));
 
       for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
 
         if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxLocal_ns)) == dataL1.end()) {
           dataL0.push_back(*i); 
         } 
       }
       
       for (buffer_type::const_iterator root = dataL1.begin(); root != dataL1.end(); ++root) {
 
         buffer_type data(1, *root);
       
         JBinder2nd<hit_type> matching = JBind2nd(match3G, *root);
       
         for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
         
           if(( root->getModuleIdentifier() != i->getModuleIdentifier() ) && matching(*i)){
             data.push_back(*i);     
           }
         }
       
         buffer_type::iterator __end1 = clusterizeWeight(data.begin() + 1, data.end(), match3G);
       
         // 4D fit
         JEstimator_t fit;
         JPoint4D  vx;
         double    chi2 = numeric_limits<double>::max();
         int       NDF  = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS;
         int       N    = getCount(data.begin(), __end1);
     
         if(distance(data.begin(), __end1) <= factoryLimit){
       
           double ymin = numeric_limits<double>::max();
       
           buffer_type::iterator __end2 = __end1;
       
           for (int n = 0; n <= numberOfOutliers && distance(data.begin(), __end2) > 
                  JEstimator_t::NUMBER_OF_PARAMETERS; ++n, --__end2) {
         
             sort(data.begin() + 1, __end1, compare);
         
             do {
               try {
             
                 fit(data.begin(), __end2);
 
                 double y = getChi2(fit, data.begin(), __end2, sigma_ns);
             
                 if (y < ymin) {
                   ymin = y;
                   vx   = fit;
                   chi2 = ymin;
                   NDF  = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
                   N    = getCount(data.begin(), __end2);
                 }
               }
               catch(JException& error) { }
           
             } while (next_permutation(data.begin() + 1, __end2, __end1, compare));
         
             ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
           }
       
         } else {
 
           const int number_of_outliers = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS - 1;
 
           buffer_type::iterator __end2 = __end1;
 
           for (int n = 0; n <= number_of_outliers; ++n) {         
 
             try{        
 
               fit(data.begin(), __end2);
               vx   = fit;
               chi2 = getChi2(fit, data.begin(), __end2, sigma_ns);
               NDF  = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
               N    = getCount(data.begin(), __end2);
 
             }
             catch(const JException& error){ }
 
             double                ymax = 0;
             buffer_type::iterator imax = __end2;
 
             for (buffer_type::iterator i = data.begin() + 1; i != __end2; ++i) {
 
               double y = getChi2(fit, *i, sigma_ns);
 
               if (y > ymax) {
                 ymax = y;
                 imax = i;
               }
             }
 
             if (ymax > STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {     
               --__end2;
               swap(*imax, *__end2);       
             } else {
               break;
             }
           }     
         }
 
         if (NDF >= 0) {
         
           const JGEOMETRY3D::JTrack3D fitResult(vx, JGEOMETRY3D::JVersor3Z());
         
           out.push_back(getFit(JHistory(JSHOWERPREFIT), fitResult, getQuality(chi2, N), NDF));
         
           out.rbegin()->setW(13, chi2);
           out.rbegin()->setW(14, N);           
         }
       }
    
       return out;
     
     }

void JRECONSTRUCTION::JShowerPrefitParameters_t::reset ( )

inlineinherited

Reset fit parameters.

Definition at line 35 of file JShowerPrefitParameters_t.hh.

     {
       factoryLimit        = 40;
       NMaxHits            = 20;
       sigma_ns            =  3.0;
       numberOfOutliers    =  4;
       TMaxLocal_ns        = 20.0;
       ctMin               = -0.7;
       roadWidth_m         = 50.0;
     }

bool JRECONSTRUCTION::JShowerPrefitParameters_t::equals ( const JShowerPrefitParameters_t & parameters ) const

inlineinherited

Equality.

Parameters

parameters fit parameters

Returns: true if equals; else false

Definition at line 52 of file JShowerPrefitParameters_t.hh.

     {
       return (this->factoryLimit          == parameters.factoryLimit                &&
               this->NMaxHits              == parameters.NMaxHits                    &&
               this->sigma_ns              == parameters.sigma_ns                    &&
               this->numberOfOutliers      == parameters.numberOfOutliers            &&         
               this->TMaxLocal_ns          == parameters.TMaxLocal_ns                &&
               this->ctMin                 == parameters.ctMin                       &&
               this->roadWidth_m           == parameters.roadWidth_m);
     }