SoftBWParticleConstraint.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * Forked from https://github.com/iLCSoft/MarlinKinfit                    *
 *                                                                        *
 * Further information about the fit engine and the user interface        *
 * provided in MarlinKinfit can be found at                               *
 * https://www.desy.de/~blist/kinfit/doc/html/                            *
 * and in the LCNotes LC-TOOL-2009-001 and LC-TOOL-2009-004 available     *
 * from http://www-flc.desy.de/lcnotes/                                   *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#ifdef MARLIN_USE_ROOT
 
#ifndef __SOFTBWPARTICLECONSTRAINT_H
#define __SOFTBWPARTICLECONSTRAINT_H
 
#include "analysis/OrcaKinFit/BaseSoftConstraint.h"
#include "analysis/OrcaKinFit/BaseFitObject.h"
 
#include<vector>
#include<cassert>
#include<limits>
 
namespace Belle2 {
  namespace OrcaKinFit {
 
    class ParticleFitObject;
 
//  Class SoftBWParticleConstraint:
 
    class SoftBWParticleConstraint: public BaseSoftConstraint {
    public:
      SoftBWParticleConstraint(double gamma_,                                            
                               double emin_ = -std::numeric_limits<double>::infinity(),   
                               double emax_ =  std::numeric_limits<double>::infinity()    
                              );
      virtual ~SoftBWParticleConstraint() {};
 
      virtual void setFOList(std::vector <ParticleFitObject*>* fitobjects_ 
                            )
      {
        for (int i = 0; i < (int) fitobjects_->size(); i++) {
          fitobjects.push_back((*fitobjects_)[i]);
          flags.push_back(1);
        }
      };
      virtual void addToFOList(ParticleFitObject& fitobject, int flag = 1
                              )
      {
        fitobjects.push_back(&fitobject);
        flags.push_back(flag);
      };
 
      virtual double getValue() const override = 0;
 
      virtual double getChi2() const override;
 
      virtual double getError() const override;
 
      virtual double getGamma() const;
 
      virtual double setGamma(double gamma_     
                             );
 
      virtual void getDerivatives(int idim,      
                                  double der[]   
                                 ) const override = 0;
      virtual void add2ndDerivativesToMatrix(double* M,     
                                             int idim       
                                            ) const override;
 
      virtual void addToGlobalChi2DerVector(double* y,    
                                            int idim     
                                           ) const override;
 
 
      void invalidateCache() const;
 
      virtual void updateCache() const;
 
      virtual bool cacheValid() const;
 
      void test1stDerivatives();
      void test2ndDerivatives();
 
      double num1stDerivative(int ifo,      
                              int ilocal,  
                              double eps   
                             );
      double num2ndDerivative(int ifo1,     
                              int ilocal1, 
                              double eps1, 
                              int ifo2,    
                              int ilocal2, 
                              double eps2  
                             );
 
 
 
      static double erfinv(double x);
 
      static double normal_quantile(double x);
      static double normal_quantile_1stderiv(double x);
      static double normal_quantile_2ndderiv(double x);
      static double normal_pdf(double x);
      static double normal_pdf_deriv(double x);
 
      double penalty(double e) const;
      double penalty1stder(double e) const;
      double penalty2ndder(double e) const;
 
      int getVarBasis() const;
 
    protected:
 
      virtual bool secondDerivatives(int i,                         
                                     int j,                        
                                     double* derivatives           
                                    ) const = 0;
      virtual bool firstDerivatives(int i,                         
                                    double* derivatives           
                                   ) const = 0;
 
 
      typedef std::vector <ParticleFitObject*> FitObjectContainer;
      typedef FitObjectContainer::iterator FitObjectIterator;
      typedef FitObjectContainer::const_iterator ConstFitObjectIterator;
      FitObjectContainer fitobjects;
      std::vector <double> derivatives;
      std::vector <int> flags;
 
      double gamma;
      double emin;   
      double emax;   
 
      mutable bool cachevalid;
      mutable double atanxmin;
      mutable double atanxmax;
      mutable double diffatanx;
 
      enum { VAR_BASIS = BaseDefs::VARBASIS_EPXYZ }; // this means that the constraint knows about E,px,py,pz
 
    };
 
  }// end OrcaKinFit namespace
} // end Belle2 namespace
 
 
#endif // __SOFTBWPARTICLECONSTRAINT_H
 
#endif // MARLIN_USE_ROOT