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
	Auxiliary class to match hit to root hit. More...

Public Types
typedef JTRIGGER::JHitR1	hit_type

typedef std::vector< hit_type >	buffer_type

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, 3)

Public Attributes
const JModuleRouter &	router

size_t	numberOfPrefits
	number of prefits More...

int	factoryLimit
	factory limit for combinatorics More...

double	sigma_ns
	time resolution [ns] More...

int	numberOfOutliers
	maximum number of outliers More...

double	TMaxLocal_ns
	time window for local coincidences [ns] More...

double	TMaxExtra_ns
	time window for extra coincidences [ns] More...

double	ctMin
	minimal cosine space angle between PMT axes More...

double	DMax_m
	maximal distance to optical module [m] More...

size_t	numberOfL1
	minimal number of L1 More...

int	pos_grid_m
	edge [m] of the position grid More...

int	pos_step_m
	step in [m] of position grid More...

int	time_grid_ns
	edge [ns] of the time grid More...

int	time_step_ns
	step in [ns] of time grid More...

size_t	numberOfGrids
	number of prefits to be used to build a grid around More...

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 52 of file JShowerPrefit.hh.

Member Typedef Documentation

typedef JTRIGGER::JHitR1 JRECONSTRUCTION::JShowerPrefit::hit_type

Definition at line 57 of file JShowerPrefit.hh.

typedef std::vector<hit_type> JRECONSTRUCTION::JShowerPrefit::buffer_type

Definition at line 58 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 68 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 81 of file JShowerPrefit.hh.

     {
 
       using namespace std;
       using namespace JPP;
 
       const double STANDARD_DEVIATIONS       = 3.0;
 
       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(DMax_m, TMaxExtra_ns);  // causality relation for showers 
   
       //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));
 
       if(dataL1.size() <= numberOfL1){
         for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
           if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxExtra_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(NDF > 0){
           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, hit_type::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, hit_type::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;
               }
             }   
           }
   
           out.push_back(getFit(JHistory(JSHOWERPREFIT), JTrack3D(vx, JGEOMETRY3D::JVersor3Z()), getQuality(chi2, N), NDF));
         
         }
       }
 
       out.select(numberOfGrids, qualitySorter);
 
       size_t solutions = out.size();
 
       for(size_t i=0; i < solutions; i++){
         for(int x = -pos_grid_m; x < pos_grid_m + pos_step_m/2.; x += pos_step_m){
           for(int y = -pos_grid_m; y < pos_grid_m + pos_step_m/2.; y += pos_step_m){
             for(int z = -pos_grid_m; z < pos_grid_m + pos_step_m/2.; z += pos_step_m){
               for(int t = -time_grid_ns; t < time_grid_ns + time_step_ns/2.; t += time_step_ns){
                 if (x != 0 || y != 0 || z != 0 || t != 0) {
                   
                   out.push_back(getFit(JHistory(JSHOWERPREFIT),
                                        JTrack3D(JPoint4D(getPosition(out[i]) + JPosition3D(x,y,z), out[i].getT()+t), JVersor3Z()),
                                        0,
                                        0));
                    
                 }
               }
             }
           }
         }
       }
       return out;
     
     }

void JRECONSTRUCTION::JShowerPrefitParameters_t::reset ( )

inlineinherited

Reset fit parameters.

Definition at line 35 of file JShowerPrefitParameters_t.hh.

     {
       factoryLimit        = 40;
       sigma_ns            =  3.0;
       numberOfOutliers    =  4;
       TMaxLocal_ns        = 20.0;
       TMaxExtra_ns        = 60.0;
       ctMin               = 0.0;
       DMax_m              = 50.0;
       numberOfPrefits     = 100;
       numberOfGrids       = 2;
       numberOfL1          = 10;
       pos_grid_m          = 10;
       pos_step_m          = 10;
       time_grid_ns        = 60;
       time_step_ns        = 60;
     }

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 59 of file JShowerPrefitParameters_t.hh.

     {
       return (this->factoryLimit          == parameters.factoryLimit                &&
               this->sigma_ns              == parameters.sigma_ns                    &&
               this->numberOfOutliers      == parameters.numberOfOutliers            &&
               this->TMaxExtra_ns          == parameters.TMaxExtra_ns                &&
               this->TMaxLocal_ns          == parameters.TMaxLocal_ns                &&
               this->ctMin                 == parameters.ctMin                       &&
               this->DMax_m                == parameters.DMax_m                      &&
               this->numberOfPrefits       == parameters.numberOfPrefits             &&
               this->numberOfGrids         == parameters.numberOfGrids               &&
               this->pos_grid_m            == parameters.pos_grid_m                  &&
               this->pos_step_m            == parameters.pos_step_m                  &&
               this->time_grid_ns          == parameters.time_grid_ns                &&
               this->time_step_ns          == parameters.time_step_ns                &&
               this->numberOfL1            == parameters.numberOfL1);
     }