JFIT::JORCAShowerFit Class Reference

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

#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 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

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

Definition at line 57 of file JORCAShowerFit.hh.

Member Typedef Documentation

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

private

Definition at line 61 of file JORCAShowerFit.hh.

Constructor & Destructor Documentation

JFIT::JORCAShowerFit::JORCAShowerFit ( )

inline

Default constructor.

Definition at line 72 of file JORCAShowerFit.hh.

    {}

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 83 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};
    }

JFIT::JORCAShowerFit::~JORCAShowerFit ( )

inline

Definition at line 97 of file JORCAShowerFit.hh.

{}

Member Function Documentation

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

inlinestatic

Function to sort different shower hypotheses by their likelihood.

Definition at line 117 of file JORCAShowerFit.hh.

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

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 144 of file JORCAShowerFit.hh.

{
  using JFIT::JShower3EZ;
  using JFIT::JPoint4D;
  using JFIT::JVersor3Z;
  using JGEOMETRY3D::JVector3D;

  const double Dmax_m               = parameters_.Dmax_m;
  const double R_Hz                 = parameters_.R_Hz;
  const std::size_t numberOfPrefits = parameters_.numberOfPrefits;

  if ( InPreFits.empty() ) return;

  typedef std::vector<JTRIGGER::JHitL0>  JDataL0_t;
  typedef std::vector<JFIT::JPMTW0> JPMTW0_t;
  JTRIGGER::JBuildL0<JTRIGGER::JHitL0>  buildL0;

  const JDETECTOR::JDetector &detector = moduleRouter_->getReference();

  JEvt::iterator __end = InPreFits.end();
  std::copy(InPreFits.begin(), __end, std::back_inserter(OutFits));
  
  if (numberOfPrefits > 0) {
    std::advance(__end = InPreFits.begin(), std::min(numberOfPrefits, InPreFits.size()));
  }
  
  std::partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
  
  JDataL0_t dataL0;
  buildL0(timeSliceBuildL0, *moduleRouter_, std::back_inserter(dataL0));
  
  for (JEvt::const_iterator shower = InPreFits.begin(); shower != __end; ++shower) {
    const JVector3D pos(getPosition(*shower)); //.rotate(R));  
        
    JFIT::JPoint4D pt(pos, shower->getT());
    JVector3D start_dir(0,0,0);
    
    const JVector3D DetC(0.0, 0.0, 117);
    const double radius = pos.getDistance(DetC);
    
    std::multiset<KM3NETDAQ::JDAQPMTIdentifier> top;
    
    double Nhits = 0;
    
    for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
      
      const JFIT::JHitW0 hit(*i, R_Hz);
      const JGEOMETRY3D::JPosition3D hit_pos(hit.getPosition());

      const double D = hit_pos.getDistance(pt.getPosition());
      const double t_res = hit.getT() - pt.getT(hit_pos);

      const JGEOMETRY3D::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 = exp((Nhits + 26)/26); 
    if(start_E > 100) start_E = 100;
    
    JFIT::JORCAShowerFit::JShowerHypothesis shower_hypothesis;
    std::vector<JFIT::JORCAShowerFit::JShowerHypothesis> shower_hypotheses;
    
    std::vector<JFIT::JPMTW0> buffer;
    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 JGEOMETRY3D::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 * JTOOLS::PI / 180;
          const double phi    = ph * JTOOLS::PI / 180;
          const double scan_E = start_E + E; // [GeV] 

          const JGEOMETRY3D::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 JGEOMETRY3D::JDirection3D conversion(dir_shifted.getX(), 
                                                     dir_shifted.getY(), 
                                                     dir_shifted.getZ());
          const JGEOMETRY3D::JRotation3D R(conversion);
          pt.rotate(R);
          
          for (JDETECTOR::JDetector::const_iterator module = detector.begin(); 
               module != detector.end(); ++module) {

            JGEOMETRY3D::JPosition3D pos(module->getPosition());
            pos.rotate(R);

            if (pt.getDistance(pos) < Dmax_m) {
              for (unsigned int i = 0; i != module->size(); ++i) {
                const KM3NETDAQ::JDAQPMTIdentifier id(module->getID(), i);
                JDETECTOR::JPMT pmt(module->getPMT(i));
                pmt.rotate(R);
                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);
          
        }
      }
    }

    std::vector<JShowerHypothesis>::iterator __end2 = shower_hypotheses.begin() + Nmin;
    std::partial_sort(shower_hypotheses.begin(), __end2, shower_hypotheses.end(), sortLogLik);

    for(std::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), JFIT::JEnergy());
      fit_.step[1] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(0.0, dir_step), JFIT::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 JGEOMETRY3D::JDirection3D conversion(startVersor.getDX(), startVersor.getDY(), startVersor.getDZ());
      JGEOMETRY3D::JRotation3D R(conversion);
      pt.rotate(R);

      for(JDETECTOR::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 KM3NETDAQ::JDAQPMTIdentifier id(module->getID(), i);
            JDETECTOR::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(), JFIT::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 JGEOMETRY3D::JTrack3D td(resPos, resDir, shower->getT());
      const double E_reco_corrected = fit_.value.getE() - 4;
      JGEOMETRY3D::JTrack3E tb(td, E_reco_corrected);
      tb.rotate_back(R);
      
      const JFIT::JFit &outFit = JFIT::getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT), 
                                              tb, -chi2, NDF,
                                              JFIT::JFitStatus_t::OKAY); 

      OutFits.push_back(outFit);
      OutFits.rbegin()->setE(E_reco_corrected);
    }
  }
  std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
}

Member Data Documentation

JRegressor_t JFIT::JORCAShowerFit::fit_

mutableprivate

Definition at line 62 of file JORCAShowerFit.hh.

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

private

Definition at line 64 of file JORCAShowerFit.hh.

JFIT::JShowerFitParameters_t JFIT::JORCAShowerFit::parameters_

private

Definition at line 65 of file JORCAShowerFit.hh.

---

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

• JORCAShowerFit.hh