StFmsClusterFitter.h

Go to the documentation of this file.

// $Id: StFmsClusterFitter.h,v 1.11 2020/08/05 16:17:23 akio Exp $
//
// $Log: StFmsClusterFitter.h,v $
// Revision 1.11  2020/08/05 16:17:23  akio
// Avoid floating point exception (from Yuri)
//
// Revision 1.10  2019/06/26 16:49:53  akio
// shower shape scaling for all shapes
//
// Revision 1.9  2018/03/02 20:27:29  akio
// Big update from      Zhanwen Zhu with new shower shape and six z slices
//
// Revision 1.8  2016/06/07 15:51:44  akio
// Making code better based on Coverity reports
//
// Revision 1.7  2015/11/05 17:54:57  akio
// Adding option to scale up shower shape function for large cells
//
// Revision 1.6  2015/10/30 21:33:56  akio
// fix parameter initialization
// adding new cluster categorization method
//
// Revision 1.5  2015/10/29 21:14:55  akio
// increase max number of clusters
// a bug fixes in valley tower association
// removing some debug comments
//
// Revision 1.4  2015/10/23 15:05:07  akio
// fixed constexpr
//
// Revision 1.3  2015/10/21 15:58:04  akio
// Code speed up (~x2) by optimizing minimization fuctions and showershape function
// Add option to merge small cells to large, so that it finds cluster at border
// Add option to perform 1photon fit when 2photon fit faield
// Add option to turn on/off global refit
// Moment analysis done without ECUTOFF when no tower in cluster exceed ECUTOFF=0.5GeV
//
// Revision 1.2  2015/09/02 15:01:32  akio
// Removing StFmsGeometry class, and now it uses StFmsDbMaker to get appropriate parameters.
//
// Revision 1.1  2015/03/10 14:38:53  jeromel
// First version of FmsUtil from Yuxi Pan - reviewd 2015/02
//
#ifndef STROOT_STFMSPOINTMAKER_STFMSCLUSTERFITTER_H_
#define STROOT_STFMSPOINTMAKER_STFMSCLUSTERFITTER_H_

#ifndef __CINT__  // Hide std::array from CINT, as CINT cannot parse C++11
#include <array>
#endif  // __CINT__
#include <list>
#include <vector>

#include "TMinuit.h"
#include "TObject.h"

#include "StFmsFittedPhoton.h"
#include "StFmsTowerCluster.h"

class TF2;
class TString;

namespace FMSCluster {  // $NMSPC
typedef std::list<StFmsFittedPhoton> PhotonList;
  //class StFmsGeometry;
class StFmsTower;
class StFmsClusterFitter : public TObject {
 public:
  StFmsClusterFitter(/*const StFmsGeometry* geometry,*/ Int_t detectorId, Float_t xw, Float_t yw, 
                     Float_t scaleShowerShapeLarge , Float_t scaleShowerShapeSmall,
                     Int_t ShowerShapeWithAngle , Int_t MergeSmallToLarge, double vertexZ);
  StFmsClusterFitter() { }
  virtual ~StFmsClusterFitter();
  TF2* showerShapeFunction();
  void setTowers(StFmsTowerCluster::Towers* towers);
  Double_t fitNPhoton(const std::vector<double>& parameteres,
                      const std::vector<double>& steps,
                      const std::vector<double>& lower,
                      const std::vector<double>& up,
                      PhotonList* photons);
  Double_t fitNPhoton(const std::vector<double>& parameters,
                      const std::vector<double>& lower,
                      const std::vector<double>& up,
                      PhotonList* photons);
#ifndef __CINT__  // Hide std::array from CINT, as CINT cannot parse C++11
  Int_t fit2Photon(const std::array<double, 7>& parameters,
                   const std::array<double, 7>& steps,
                   const std::array<double, 7>& lower,
                   const std::array<double, 7>& upper,
                   PhotonList* photons);
#endif  // __CINT__

  static inline Double_t energyDepositionInTowerOLD(Double_t x, Double_t y, Double_t* par);

  static Double_t energyDepositionInTower(Double_t* x, Double_t* par);
  static inline Double_t energyDepositionInTower(Double_t x, Double_t y,Double_t xun, Double_t yun, Double_t* parmeters,Int_t  MergeSmallToLarge,double vertexz);
  static inline Double_t energyDepositionInTowerSingleLayer(Double_t x, Double_t y, Double_t* parameters);

  static int maxNFittedPhotons();

 private:
  StFmsClusterFitter(const StFmsClusterFitter&);
  StFmsClusterFitter& operator=(const StFmsClusterFitter&);
  int runMinuitMinimization();
  static inline Double_t energyDepositionDistribution(Double_t x, Double_t y, Double_t* par);
  static Double_t energyDepositionDistribution(Double_t* x, Double_t* par);
  static void minimizationFunctionNPhoton(Int_t& npar, Double_t* grad,
                                          Double_t& fval, Double_t* par,
                                          Int_t iflag);
  static void minimizationFunction2Photon(Int_t& nparam, Double_t* grad,
                                          Double_t& fval, Double_t* param,
                                          Int_t iflag);
  template<class Container>
  int setMinuitParameter(int index, const TString& name,
                         const Container& parameters,
                         const Container& steps,
                         const Container& lower,
                         const Container& upper);
  int readMinuitParameters(std::vector<double>& parameters,
                           std::vector<double>& errors);
  TMinuit mMinuit;  
  static StFmsTowerCluster::Towers* mTowers;  
  static Double_t mEnergySum;
  Float_t mScaleShowerShapeLarge=1.0;
  Float_t mScaleShowerShapeSmall=1.0;
  Int_t mShowerShapeWithAngle;
  Int_t mMergeSmallToLarge;


  ClassDef(StFmsClusterFitter, 0)
};  // class StFmsClusterFitter

inline double asmsqrt(double x) { //slower than std::sqrt()                                                                                                                       
    __asm__ ("fsqrt" : "+t" (x));
    return x;
}

// Returns a * f(x,y,b) as defined here:
// https://drupal.star.bnl.gov/STAR/blog/leun/2010/aug/02/fms-meeting-20100802
inline double showerShapeComponent(double x, double y, double a, double b) {
    if (a == 0.0 || b == 0.0) return 0.0;
    return a * atan(x * y / (b * sqrt(b * b + x * x + y * y)));
}

inline Double_t StFmsClusterFitter::energyDepositionDistribution(double x, double y, double* parameters){
    constexpr double ootwopi = 1.0/2.0/3.14159265358979323846;
    return (     showerShapeComponent(x, y, parameters[1], parameters[4])
               + showerShapeComponent(x, y, parameters[2], parameters[5])
               + showerShapeComponent(x, y, parameters[3], parameters[6]) ) * ootwopi;
}

inline Double_t StFmsClusterFitter::energyDepositionInTowerOLD(Double_t x, Double_t y, Double_t* parameters){
    //old version without incident angle effect
    const double w = parameters[0]/2.0;
    return  energyDepositionDistribution(x-w,y-w,parameters)
          - energyDepositionDistribution(x-w,y+w,parameters)
          - energyDepositionDistribution(x+w,y-w,parameters)
          + energyDepositionDistribution(x+w,y+w,parameters);
}

//x,y are cell position, and xun,yun are photon position, in cell coordinate at Z=ZMAX
inline Double_t StFmsClusterFitter::energyDepositionInTower(Double_t x, Double_t y,Double_t xun, Double_t yun, 
                Double_t* parameters,Int_t mMerge ,double vertexz){
    //Double_t para[60];
    Double_t *Zc;
    Double_t ZcS[6] = {720.45,727.95,735.45,742.95,750.45,757.95};
    Double_t ZcL[6] = {722.98,733.01,743.04,753.07,763.10,773.13};
    
    //for(Int_t i = 0; i < 60; i++){ para[i]=parameters[i]; }
    Double_t Zmax,yoff,xoff;
    if (parameters[0] > 4) { 
        yoff=98.8; Zmax=735.45; xoff=0.3; Zc=ZcL;
    }else{
        yoff=46.5; Zmax=735.45; xoff=0.93; Zc=ZcS;
    }
    if (mMerge>0) yoff=98.8;
    
    Double_t tanx = ( xun + xoff) / (Zmax-vertexz); //recalculate tanx and tany for each iteration
    Double_t tany = ( yoff - yun) / (Zmax-vertexz); //use Zmax instead of lzz, 06/25/2012 Yuxi
    
    Double_t sum = 0.0; 
    for(Int_t i = 0; i < 6; i++){
        Double_t xc = xun + tanx*(Zc[i] - Zmax);
        Double_t yc = yun - tany*(Zc[i] - Zmax); //large (positive) tany corresponds to smaller row # (yc)
        Int_t istart = i*10;
        sum += ( energyDepositionInTowerSingleLayer(x-xc,y-yc,&parameters[istart])  *  parameters[istart+7]) ;
    }
    return sum;
}

inline Double_t StFmsClusterFitter::energyDepositionInTowerSingleLayer(Double_t x, Double_t y, Double_t* parameters){
    const double w = parameters[0];
    double ww1=0;
    double ww2=0;
    if (w < 4) {ww1=3.822/2; ww2=3.875/2; }
    else       {ww1=5.812/2; ww2=5.812/2; }
    
    return  energyDepositionDistribution(x-ww1,y-ww2,parameters)
          - energyDepositionDistribution(x-ww1,y+ww2,parameters)
          - energyDepositionDistribution(x+ww1,y-ww2,parameters)
          + energyDepositionDistribution(x+ww1,y+ww2,parameters);
}


}  // namespace FMSCluster
#endif  // STROOT_STFMSPOINTMAKER_STFMSCLUSTERFITTER_H_