class to handle the third step of the shower reconstruction, mainly dedicated for ORCA More...

#include <JShowerEnergyPrefit.hh>

Inheritance diagram for JRECONSTRUCTION::JShowerEnergyPrefit:

Classes
struct	JResult
	Auxiliary class for energy estimation. More...

Public Member Functions
	JShowerEnergyPrefit (const JShowerEnergyPrefitParameters_t &parameters, const JModuleRouter &router, const JSummaryRouter &summary, const std::string pdfFile, const int debug=0)
	Parameterized constructor. More...

JEvt	operator() (const KM3NETDAQ::JDAQEvent &event, const JFIT::JEvt &in)
	Declaration of the member function that actually performs the reconstruction. More...

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

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

	ClassDef (JShowerEnergyPrefitParameters_t, 1)

Public Attributes
const JModuleRouter &	router

const JSummaryRouter &	summary

double	TMax_ns

double	TMin_ns

double	roadWidth_m

double	Emin_GeV

double	Emax_GeV

double	resolution

double	R_Hz

int	mestimator

size_t	numberOfPrefits

size_t	numberOfOutfits

Private Types
typedef JRegressor< JEnergy, JSimplex >	JRegressor_t

Detailed Description

class to handle the third step of the shower reconstruction, mainly dedicated for ORCA

Definition at line 73 of file JShowerEnergyPrefit.hh.

Member Typedef Documentation

typedef JRegressor<JEnergy, JSimplex> JRECONSTRUCTION::JShowerEnergyPrefit::JRegressor_t

private

Definition at line 78 of file JShowerEnergyPrefit.hh.

Constructor & Destructor Documentation

JRECONSTRUCTION::JShowerEnergyPrefit::JShowerEnergyPrefit	(	const JShowerEnergyPrefitParameters_t &	parameters,
		const JModuleRouter &	router,
		const JSummaryRouter &	summary,
		const std::string	pdfFile,
		const int	debug = `0`
	)

inline

Parameterized constructor.

Parameters

parameters	struct that holds default-optimized parameters for the reconstruction, available in $JPP_DATA.
router	module router, this is built via detector file.
summary	summary router
pdfFile	PDF file
debug	debug

Definition at line 92 of file JShowerEnergyPrefit.hh.

                                                                  :
       JShowerEnergyPrefitParameters_t(parameters),
       JRegressor_t(pdfFile),
       router(router),
       summary(summary)
     {
       using namespace JPP;
       
       JRegressor_t::debug  = 1;
       JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
       JRegressor_t::Vmax_npe = 20.0;
       JRegressor_t::MAXIMUM_ITERATIONS = 1000;
       
       this->estimator.reset(getMEstimator(parameters.mestimator));
     }

Member Function Documentation

JEvt JRECONSTRUCTION::JShowerEnergyPrefit::operator()	(	const KM3NETDAQ::JDAQEvent &	event,
		const JFIT::JEvt &	in
	)

inline

Declaration of the member function that actually performs the reconstruction.

Parameters

event	which is a JDAQEvent
in	input fits, which should contain a vertex fit

Definition at line 151 of file JShowerEnergyPrefit.hh.

     {
       using namespace std;
       using namespace JPP;
 
       typedef vector<JHitL0>  JDataL0_t;
       JBuildL0<JHitL0>        buildL0;
 
       JEvt out;
 
       const JDetector &detector = router.getReference();
 
       JDataL0_t dataL0;
       buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
   
       for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
 
         JShower3E sh(getPosition(*shower), getDirection(*shower), shower->getT(), shower->getE());
         
         multiset<JDAQPMTIdentifier> top;
 
         vector<JShowerNPEHit> data;
         vector<JShowerNPEHit> buffer;
 
         const JModel<JPoint4D> match(JPoint4D(sh.getPosition(), sh.getT()), roadWidth_m, JRegressor_t::T_ns);
 
         for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
 
           if (match(*i)) {
             top.insert(i->getPMTIdentifier());
           }
         }
 
         JDetectorSubset_t subdetector(detector, sh.getPosition(), roadWidth_m);
         
         for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
 
           for (size_t i = 0; i != module->size(); ++i) {
             
             const size_t count  = top.count(JDAQPMTIdentifier(module->getID(), i));
 
             const double rate_Hz = summary.getRate(JDAQPMTIdentifier(module->getID(), i));
 
             data.push_back(JShowerNPEHit(this->getNPE(module->getPMT(i), rate_Hz), count));
           }
         }
 
         const int NDF = data.size() - 1;
 
         if (NDF >= 0) {
           
           // 5-point search between given limits
           const int      N  =  5;
 
           JResult result[N];
 
           for (int i = 0; i != N; ++i) {
 
             result[i].x = log10(Emin_GeV + i * (Emax_GeV - Emin_GeV) / (N-1));
 
           }
           
           do{
   
             int j = 0;
 
             for (int i = 0; i != N; ++i) {
               
               if (!result[i]) {
 
                 result[i].chi2 = (*this)(result[i].x, data.begin(), data.end());
               
               }
               
               if (result[i].chi2 < result[j].chi2) {
                 j = i;
               }
             }       
 
             // squeeze range
 
            switch (j) {
 
            case 0:
              result[4] = result[1];
              result[2] = JResult(0.5 * (result[0].x + result[4].x));
              break;
 
            case 1:
              result[4] = result[2];
              result[2] = result[1];
              break;
 
            case 2:
              result[0] = result[1];
              result[4] = result[3];
              break;
 
            case 3:
              result[0] = result[2];
              result[2] = result[3];
              break;
 
            case 4:
              result[0] = result[3];
              result[2] = JResult(0.5 * (result[0].x + result[4].x));
              break;
            }
            
            result[1] = JResult(0.5 * (result[0].x + result[2].x));
            result[3] = JResult(0.5 * (result[2].x + result[4].x));
            
           } while (result[4].x - result[0].x > resolution);
       
         
           if (result[1].chi2 != result[3].chi2) {
             
             result[2].x += 0.25 * (result[3].x - result[1].x) * (result[1].chi2 - result[3].chi2) / (result[1].chi2 + result[3].chi2 - 2*result[2].chi2);
 
             result[2].chi2  = (*this)(result[2].x, data.begin(), data.end());
           
           }
          
           //      const double chi2 = result[2].chi2;  // this is not used because fits are sorted wrt the position fit
           const double E    = result[2].x.getE();
 
           JShower3E sh_fit(sh.getPosition(), sh.getDirection(), sh.getT(), E);
         
           // the fits of this reco step are sorted wrt the previous reco step 
           // because otherwise it tends to degrade the position reco performances
           out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERENERGYPREFIT), sh_fit, shower->getQ(), 
                                shower->getNDF(), sh_fit.getE())); 
 
         }
       }
       
       return out;
 
     }

void JRECONSTRUCTION::JShowerEnergyPrefitParameters_t::reset ( )

inlineinherited

Reset fit parameters.

Definition at line 37 of file JShowerEnergyPrefitParameters_t.hh.

     {
       TMax_ns                 =  500;
       TMin_ns                 = -500;
       roadWidth_m             =   80;
       Emin_GeV                =    1;
       Emax_GeV                =  100;
       resolution              =    0.01; 
       R_Hz                    =   10.0e3;
       mestimator              = JFIT::EM_NULL;
       numberOfPrefits         =    0;
       numberOfOutfits         =    1;
     }

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

inlineinherited

Equality.

Parameters

parameters fit parameters

Returns: true if equals; else false

Definition at line 57 of file JShowerEnergyPrefitParameters_t.hh.

     {
       return (this->TMax_ns                  == parameters.TMax_ns               &&         
               this->TMin_ns                  == parameters.TMin_ns               &&         
               this->numberOfPrefits          == parameters.numberOfPrefits       &&         
               this->numberOfOutfits          == parameters.numberOfOutfits       &&         
               this->roadWidth_m              == parameters.roadWidth_m           &&
               this->Emin_GeV                 == parameters.Emin_GeV              &&
               this->Emax_GeV                 == parameters.Emax_GeV              &&
               this->resolution               == parameters.resolution            &&
               this->R_Hz                     == parameters.R_Hz                  && 
               this->mestimator               == parameters.mestimator);
     }