JFIT::JORCAShowerFit Class Reference

#include <JORCAShowerFit.hh>

Classes
struct	JShowerHypothesis
	Container of shower hypothesis: During the fit the algorithm creates a grid in direction and energy space, for each of these points a shower hypothesis is saved: direction, energy and LogLikelihood value. More...

Public Member Functions
	JORCAShowerFit ()
	Default constructor. More...
	JORCAShowerFit (const JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > &moduleRouter, const JFIT::JShowerFitParameters_t &parameters, const std::string pdfFile)
	Parameterized constructor. More...
	~JORCAShowerFit ()
void	getJEvt (const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0, JFIT::JEvt &InPreFits, JFIT::JEvt &OutFits) const
	Declaration of the member function that actually performs the reconstruction. More...

Static Public Member Functions
static double	estimateInitialEnergy (int nHits)
	Function to estimate the initial em shower energy. More...
static bool	sortLogLik (JShowerHypothesis A, JShowerHypothesis B)
	Function to sort different shower hypotheses by their likelihood. More...

Private Types
typedef JFIT::JRegressor< JFIT::JShower3EZ, JFIT::JSimplex >	JRegressor_t

Private Attributes
JRegressor_t	fit_
JLANG::JSharedPointer< const JDETECTOR::JModuleRouter >	moduleRouter_
JFIT::JShowerFitParameters_t	parameters_

Detailed Description

Definition at line 64 of file JORCAShowerFit.hh.

Member Typedef Documentation

JRegressor_t

typedef JFIT::JRegressor<JFIT::JShower3EZ, JFIT::JSimplex> JFIT::JORCAShowerFit::JRegressor_t

private

Definition at line 68 of file JORCAShowerFit.hh.

Constructor & Destructor Documentation

JORCAShowerFit() [1/2]

JFIT::JORCAShowerFit::JORCAShowerFit ( )

inline

Default constructor.

Definition at line 79 of file JORCAShowerFit.hh.

    {}

JORCAShowerFit() [2/2]

JFIT::JORCAShowerFit::JORCAShowerFit ( const JLANG::JSharedPointer< const JDETECTOR::JModuleRouter > &  moduleRouter, const JFIT::JShowerFitParameters_t &  parameters, const std::string  pdfFile  )

inline

Parameterized constructor.

Parameters

moduleRouter	JSharedPointer of JModuleRouter, this is built via detector file.
parameters	struct that holds default-optimized parameters for the reconstruction, which can be found in $JPP_DATA.
pdfFile	PDF file

Definition at line 90 of file JORCAShowerFit.hh.

                                           :
      moduleRouter_(moduleRouter),
      parameters_(parameters)
    {
      JRegressor_t::debug  =1;
      JRegressor_t::T_ns.setRange(-parameters.Tmax_ns,parameters.Tmax_ns);
      JRegressor_t::Vmax_npe = 20.0;
      JRegressor_t::MAXIMUM_ITERATIONS = 1000;
      
      fit_ = {pdfFile};
      fit_.estimator.reset(new JMEstimatorNull());
    }

~JORCAShowerFit()

JFIT::JORCAShowerFit::~JORCAShowerFit ( )

inline

Definition at line 105 of file JORCAShowerFit.hh.

{}

Member Function Documentation

estimateInitialEnergy()

static double JFIT::JORCAShowerFit::estimateInitialEnergy ( int  nHits )

inlinestatic

Function to estimate the initial em shower energy.

Definition at line 125 of file JORCAShowerFit.hh.

                                                  {
 
      double Energy = exp((nHits + 26)/26);
      if(Energy > 100) Energy = 100;
 
      return Energy;
    }

sortLogLik()

static bool JFIT::JORCAShowerFit::sortLogLik ( JShowerHypothesis  A, JShowerHypothesis  B  )

inlinestatic

Function to sort different shower hypotheses by their likelihood.

Definition at line 138 of file JORCAShowerFit.hh.

    {  
      return(A.LogLik < B.LogLik);
    };    

getJEvt()

void JFIT::JORCAShowerFit::getJEvt	(	const KM3NETDAQ::JDAQTimeslice &	timeSliceBuildL0,
		JFIT::JEvt &	InPreFits,
		JFIT::JEvt &	OutFits
	)		const

Declaration of the member function that actually performs the reconstruction.

member function definition

Parameters

timeSliceBuildL0	which is contructed via a JDAQEvent
InPreFits	input fits, normally produced by JORCAShowerPositionFit
OutFits	output fits

Definition at line 165 of file JORCAShowerFit.hh.

{
 
  const double Dmax_m               = parameters_.Dmax_m;
  const double R_Hz                 = parameters_.R_Hz;
  const size_t numberOfPrefits = parameters_.numberOfPrefits;
 
  if ( InPreFits.empty() ) return;
 
  typedef vector<JTRIGGER::JHitL0>  JDataL0_t;
  typedef vector<JFIT::JPMTW0> JPMTW0_t;
  JBuildL0<JTRIGGER::JHitL0>  buildL0;
 
  const JDETECTOR::JDetector &detector = moduleRouter_->getReference();
 
  JEvt::iterator __end = InPreFits.end();
  std::copy(InPreFits.begin(), __end, back_inserter(OutFits));
  
  if (numberOfPrefits > 0) {
    advance(__end = InPreFits.begin(), min(numberOfPrefits, InPreFits.size()));
  }
  
  partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
  
  JDataL0_t dataL0;
  buildL0(timeSliceBuildL0, *moduleRouter_, back_inserter(dataL0));
  
  for (JEvt::const_iterator shower = InPreFits.begin(); shower != __end; ++shower) {
    const JVector3D pos(getPosition(*shower));   
        
    JDetectorSubset_t subdetector(detector, pos, Dmax_m);
 
    JPoint4D pt(pos, shower->getT());
    JVector3D start_dir(0,0,0);
 
    // get the detector center and find the distance between the reco vertex and it.
    const JCenter3D center(detector.begin(), detector.end());
    const double radius = pos.getDistance(center);
 
    multiset<KM3NETDAQ::JDAQPMTIdentifier> top;
    
    double Nhits = 0;
    
    for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
      
      const JHitW0 hit(*i, R_Hz);
      const JPosition3D hit_pos(hit.getPosition());
 
      const double D = hit_pos.getDistance(pt.getPosition());
      const double t_res = hit.getT() - pt.getT(hit_pos);
 
      const JDirection3D photonDir(hit_pos - pt.getPosition());
      const double cosT = photonDir.getDot(hit.getDirection());
      
      if(D < Dmax_m && (t_res >= -25 && t_res <= 25) && (cosT >= -1 && cosT <= 0.1)){
        top.insert(hit.getPMTIdentifier());
        const JVector3D d(photonDir.getDX(), photonDir.getDY(), photonDir.getDZ());
        start_dir.add(d);
      }
      if(D < Dmax_m && (t_res > -30 && t_res < 20)){
        Nhits++;
      }
    }
    
    double start_E = estimateInitialEnergy(Nhits);
 
    // rotate subdetector wrt start_dir, subdetector is already subtracted by pos
 
    const JDirection3D conversion(start_dir.getX(), start_dir.getY(), start_dir.getZ());
    const JRotation3D R(conversion);
 
    vector<JPMTW0> buffer;
 
    JFIT::JORCAShowerFit::JShowerHypothesis shower_hypothesis;
    vector<JFIT::JORCAShowerFit::JShowerHypothesis> shower_hypotheses;
    
    double chi2 = 0.0;
    double max_theta = 0, max_phi = 0, scan_step = 0;
    double Nmin = 0;
    
    if(start_E > 10 && radius < 100){
      max_theta = 30;
      max_phi = 30;
      scan_step = 5;
      Nmin = 10;
    } else {
      max_theta = 35;
      max_phi = 35;
      scan_step = 5;
      Nmin = 20;
    }
 
    const JAngle3D start_dir_angles(start_dir.getX(), start_dir.getY(), start_dir.getZ());
    
    for(double E = -15; E <= 5; E += 7.5){  // scan in energy
      for(double th = -max_theta; th <= max_theta; th += scan_step){
        for(double ph = -max_phi; ph <= max_phi; ph += scan_step){
          
          buffer.clear();
          chi2 = 0.0;
          
          //todo: use precise value of pi
          const double theta  = th * PI / 180;
          const double phi    = ph * PI / 180;
          const double scan_E = start_E + E; // [GeV] 
 
          const JAngle3D dir_shifted_angles(start_dir_angles.getTheta() + theta, 
                                                         start_dir_angles.getPhi() + phi);
          const JVector3D dir_shifted(dir_shifted_angles.getDX(), 
                                      dir_shifted_angles.getDY(), 
                                      dir_shifted_angles.getDZ());
          const JDirection3D conversion(dir_shifted.getX(), 
                                                     dir_shifted.getY(), 
                                                     dir_shifted.getZ());
          const JRotation3D R(conversion);
          pt.rotate(R);
          
          for (JDetectorSubset_t::iterator module = subdetector.begin();
               module != subdetector.end(); ++module) {
      
            module->rotate(R);
            JPosition3D module_pos(module->getPosition());
      
            for (unsigned int i = 0; i != module->size(); ++i) {
              const JDAQPMTIdentifier id(module->getID(), i);
              JPMT pmt(module->getPMT(i));
              buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
            }
          }
          
          for(JPMTW0_t::iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt ){
            chi2 += fit_(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), *pmt);
          }
          
          shower_hypothesis.LogLik = chi2;
          shower_hypothesis.Energy = scan_E;
          shower_hypothesis.Direction = dir_shifted;
          
          shower_hypotheses.push_back(shower_hypothesis);
 
          pt.rotate_back(R);
          
        }
      }
    }
 
    vector<JShowerHypothesis>::iterator __end2 = shower_hypotheses.begin() + Nmin;
    partial_sort(shower_hypotheses.begin(), __end2, shower_hypotheses.end(), sortLogLik);
 
    for(vector<JShowerHypothesis>::iterator hp = shower_hypotheses.begin(); hp != __end2; ++hp){
 
      fit_.step.resize(3);
      const double dir_step = 0.05;
      fit_.step[0] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(dir_step, 0.0), JEnergy());
      fit_.step[1] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(0.0, dir_step), JEnergy());
      fit_.step[2] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(), 0.05);
      
      const double scan_E = hp->Energy;
 
      buffer.clear();
      
      const JDirection3D startVersor(hp->Direction);
      
      const JDirection3D conversion(startVersor.getDX(), startVersor.getDY(), startVersor.getDZ());
      JRotation3D R(conversion);
      pt.rotate(R);
 
      for(JDetector::const_iterator module = detector.begin(); 
          module != detector.end(); ++module) {
        JPosition3D pos(module->getPosition());
        pos.rotate(R);
 
        if (pt.getDistance(pos) < Dmax_m) {
          for (unsigned int i = 0; i != module->size(); ++i) {
            const JDAQPMTIdentifier id(module->getID(), i);
            JPMT pmt(module->getPMT(i));
            pmt.rotate(R);
            buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
          }
        }
      }
      
      const int NDF = buffer.size() - fit_.step.size();
      
      chi2 = fit_(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), buffer.begin(), buffer.end());
      
      pt.rotate_back(R);
      
      const JVector3D resPos(fit_.value);
      const JVersor3D resDir(fit_.value.getDX(), fit_.value.getDY(), fit_.value.getDZ());
      
      const JTrack3D td(resPos, resDir, shower->getT());
      double E_reco_corrected = abs(fit_.value.getE() - 4);
 
      // Now that JShowerFit uses a Param PDF, the energy has to be corrected with a linear function obtained
      // from a fit to the reconstruction error: this can be improved.
 
      const double p0 = -2.028;
      const double p1 = 0.7683;
      if(E_reco_corrected > 6){
        E_reco_corrected = (E_reco_corrected - p0) / (1 + p1);
      }
 
      JTrack3E tb(td, E_reco_corrected);
      tb.rotate_back(R);
      
      const JFit &outFit = getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT), 
                                  tb, -chi2, NDF, E_reco_corrected,
                                  JFitStatus_t::OKAY); 
 
      OutFits.push_back(outFit);
    }
  }
  sort(OutFits.begin(), OutFits.end(), qualitySorter);
}

Member Data Documentation

fit_

JRegressor_t JFIT::JORCAShowerFit::fit_

mutableprivate

Definition at line 69 of file JORCAShowerFit.hh.

moduleRouter_

JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> JFIT::JORCAShowerFit::moduleRouter_

private

Definition at line 71 of file JORCAShowerFit.hh.

parameters_

JFIT::JShowerFitParameters_t JFIT::JORCAShowerFit::parameters_

private

Definition at line 72 of file JORCAShowerFit.hh.

---

The documentation for this class was generated from the following file:

• JORCAShowerFit.hh