PeakAna.h

/*
Author: David Kapukchyan
@[October 20, 2022]
> Added doxygen style comments

@[September 29, 2022]
> Got rid of the virtual painter as it is no longer needed

@[August 15, 2022]
> Small fix to GetPossiblePeaks() so that it does not stop at slope==0 on the decreasing part of the peak

@[June 27, 2022]
> Correctly sets x value of the end point of the found peak. This can result in overlaps when drawing the peak lines.

@[June 26, 2022]
> Changed GetPossiblePeaks() to check if dealing with last point. If there is an active start time and loop reaches last graph point, last graph point becomes the end of that peak window.

@[June 23, 2022]
> Added MergeByProbability() which can be used to merge peaks by probability after GetPossiblePeaks() has been called.

@[June 12, 2022]
> Implemented TObject's "Clone" and "Copy" functions so that "DrawClone" will correctly create a clone of this object for drawing purposes. Also removed setting the found peak line color from 'AnalyzeForPeak()' to *PeakAnaPainter*

@[June 8, 2022]
> Implemented a Gaussian Filter.

@[May 10, 2022]
> Implemented the ability to color and style the graph and peak lines. Changed 'GetPossiblePeaks' to use swap rather than return a vector. 'GetPeak' returns a reference to the PeakWindow since not returning a reference actually returned a copy

@[May 9, 2022]
> Implemented the "Mean" filter from Rtools into this class, and made "filter" variables. Also got rid of 'scale' in 'AnalyzeForPeak' and made it its own variable 'BaselineSigmaScale'

@[March 9, 2022]
> Moved the draw functions to their own separate painter class

@[February 25, 2022]
> Added methods for using the peak "chirality" from peak window

@[February 18, 2022]
> Moved *PeakWindow* to it's own implementation file

@[February 7, 2022]
> Added ConvertToGraph/Ana functions that is useful for super noisy data where each PeakWindow can be thought of as a point in some graph. This graph essentially captures the shape of the data without the noise.

@[February 2, 2022]
> Changed the probability function for peak window from Exp[-x]*Erfc[y] to 1/(x^2+1)*Erfc[y]. This was done as Exp[-x] was going to zero too quickly and thus made it difficult to find good parameters with reasonable probabilites. Also made some modifications to how things are drawn

@[January 29, 2022]
> Changes to how drawing functions work and a compare for PeakWindow

@[January 28, 2022]
> Fixes to peak tunnel method includng adding a tunnel sigma variable since using baseline sigma did not work as expected

@[January 24, 2022]
> Added draw methods

@[January 23, 2022]
> Added the peak tunneling method for skipping over noise data

@[January 19, 2022]
> First instance of this class that is adopted from work on *StFcsPulseFit*

@[January 21, 2022]>Inherit PeakAna from TGraph??
*/
#ifndef PEAKANA_H
#define PEAKANA_H

//C++ Headers
#include <iostream>
#include <vector>
#include <utility> //std::move (Needs C++ 11) and std::pair
#include <iterator>

//ROOT Headers
#include "TObject.h"
#include "TAttLine.h"
#include "TKey.h"
#include "TH1.h"
#include "TPaveText.h"

#include "St_base/StMessMgr.h"

//Custom Headers
#include "PeakWindow.h"
class PeakAnaPainter;

class PeakAna : public TObject
{
public:
  PeakAna();
  PeakAna( int size, double* xvals, double* yvals ); 
  PeakAna( TGraph* Sig ); 
  PeakAna( TH1* hist);    
  PeakAna(const PeakAna &OldAna,TGraph* graph=0); 
  PeakAna& operator=(const PeakAna& rhs);         
  virtual ~PeakAna();
  
  virtual void Copy(TObject& obj) const;              
  virtual TObject* Clone(const char* newname) const;  

  virtual void AddPeakStats(TPaveText* pave, const char* opt="");
  
  static TGraph* ConvertHistToGraph(TH1* graph, UInt_t numavgs=1);
  
  static TGraph* ConvertPeaksToGraph(const std::vector<PeakWindow> &Peaks);
  void ConvertPeaksToGraph();  
  static PeakAna ConvertPeaksToAna(const PeakAna &Ana);  
  PeakAna ConvertPeaksToAna();

  virtual void GetPossiblePeaks();
  
  virtual Int_t SearchForPeak(const std::vector<PeakWindow> &PossiblePeaks );
  Int_t SearchForPeak(const std::vector<PeakWindow> &PossiblePeaks, const PeakWindow& search);
  Int_t SearchForPeak(const std::vector<PeakWindow> &PossiblePeaks, Double_t peak, Double_t width);
  
  //Class methods can be used to store data and preserve found peaks for fitting
  virtual Int_t AnalyzeForPeak();
  virtual Int_t AnalyzeForPeak(Double_t peak, Double_t width);//<! same as #AnalyzeForPeak() but also sets search parameters for peak
  virtual Int_t AnalyzeForNoisyPeak();//<! First calls #ConvertToAna() function and then #AnalyzeForPeak() on the converted data

  virtual void MergeByProbability(std::vector<PeakWindow>& newpeaks) const;
  virtual void MergeByChirality(std::vector<PeakWindow>& newpeaks) const;
  virtual short MergeLeftOrRight(UInt_t index) const;
    
  Double_t Baseline()const{return mBaseline;}                     
  Double_t BaselineSigma()const{return fabs(mBaselineSigma); }    
  Double_t BaselineSigmaScale()const{return mBaselineSigmaScale;} 
  Double_t MinX()const{ return mXRangeMin; }                      
  Double_t MinY()const{ return mYRangeMin; }                      
  Double_t MaxX()const{ return mXRangeMax; }                      
  Double_t MaxY()const{ return mYRangeMax; }                      
  Double_t SearchPeak()const{ return mSearch.mStartX; }           
  Double_t SearchWidth()const{ return mSearch.mEndX; }            
  Double_t TunnelScale()const{ return mTunnelScale; }             
  Double_t TunnelSigma()const{ return mTunnelSigma; }             
  Double_t TunnelThreshold()const{ return mTunnelThreshold; }     

  virtual void Draw(Option_t *opt="");
  virtual void Paint(Option_t *opt=""); //<! see #Draw() for options
  PeakAnaPainter* GetPainter(Option_t *opt=""); //<! Creates a #PeakAnaPainter if one doesn't exist @return #mPainter
  
  void ForceInternal(){ mInternalSignal=true; } //<! Call this when you need to tell this class to delete the internal TGraph object
  PeakAna* GausFilter(Int_t sizeavgs=0, bool copy=true);
  bool GoodWindow(); 
  virtual TGraph* GetData()const{return mG_Data;}
  UInt_t GetDebug()const{return mDEBUG;}         
  const PeakWindow& GetPeak(UInt_t peakidx)const{return mPeaks.at(peakidx);} 
  PeakWindow& GetPeak(UInt_t peakidx){return mPeaks.at(peakidx);}            
  UInt_t GetFilter()const{ return mFilter; }                                 
  Int_t GetFilterScale()const{ return mFilterScale; }                        

  Int_t FoundPeakIndex()const{return mComputedIndex;}                        
  virtual void GetXYMax(Double_t xmin, Double_t xmax);                       
  Double_t MaxXval(){if( mMaxX<mXRangeMax ){GetXYMax(mXRangeMin,mXRangeMax); } return mMaxX; } 
  Double_t MaxYval(){if( mMaxY<mYRangeMin ){GetXYMax(mXRangeMin,mXRangeMax);} return mMaxY; }  
  
  PeakAna* MeanFilter(Int_t sizeavgs=0, bool copy=true);
  int NPeaks()const{return mPeaks.size();}  
  bool ValidPeakIdx() const;                
  
  Double_t PeakStart(){if( mComputedIndex<0 ){this->AnalyzeForPeak();} return mFoundPeak.mStartX;}
  Double_t PeakEnd(){  if( mComputedIndex<0 ){this->AnalyzeForPeak();} return mFoundPeak.mEndX;}  
  Double_t PeakX();  
  Double_t PeakY();  
  void PeakXY(Double_t &xval, Double_t &yval);  

  bool PeakTunnel(const PeakWindow &window) const;
  
  Double_t PeakProb(const PeakWindow& window, Double_t scale, Double_t sigma ) const;
  
  virtual void Print(Option_t* opt="") const; 
  void ResetPeak();                           
  
  void SetDebug(UInt_t level){ mDEBUG=level; }   

  virtual void SetData(TGraph* graph);

  virtual void SetData(TH1* hist, UInt_t numavgs=1);
  
  void SetBaseline(Double_t base, Double_t sigma);   
  void SetBaselineSigmaScale(Double_t scale);        
  void SetContinuity(Double_t val){ mDeltaX = val; } 
  
  void SetFilter( UInt_t filter, Int_t scale, Double_t sigma=0 );
  
  void SetRange( Double_t xmin, Double_t ymin, Double_t xmax, Double_t ymax);
  
  void SetSearchWindow(Double_t peak, Double_t width);
  
  void SetTunnelScale(Double_t value);                 
  void SetTunnelSigma(Double_t value);                 
  void SetTunnelPars(Double_t scale, Double_t sigma);  
  void SetTunnelThreshold(Double_t value);             
  
  Double_t ChiralityPeakScale() const { return mChiralityPeakScale; }
  Double_t ChiralityScale() const { return mChiralityScale; }
  Double_t ChiralityProbScale() const { return mChiralityProbScale; }
  Double_t ChiralityThreshold() const { return mChiralityThreshold; }
  void SetChiralityPeakScale(Double_t v) { mChiralityPeakScale=v; }
  void SetChiralityScale(Double_t v) { mChiralityScale=v;}
  void SetChiralityProbScale(Double_t v){ mChiralityProbScale=v;}
  void SetChiralityThreshold(Double_t v){ if(v<1){mChiralityThreshold=v;} }
  
  void SetPeak(const Int_t peakpoint, const Double_t peakx );
  void SetWindow(const Int_t start, const Int_t end );
  void SetWindow(PeakWindow window); 
  
  //Styling for graph
  Color_t GetLineColor() const;
  Style_t GetLineStyle() const;
  Width_t GetLineWidth() const;
  
  Color_t GetFillColor() const;
  Style_t GetFillStyle() const;
  
  Color_t GetMarkerColor() const;
  Size_t GetMarkerSize() const;
  Style_t GetMarkerStyle() const;
  
  void SetLineColor(Color_t color);
  void SetLineColorAlpha(Color_t color, Float_t alpha);
  void SetLineStyle(Style_t style);
  void SetLineWidth(Width_t width);
  
  void SetFillColor(Color_t color);
  void SetFillColorAlpha(Color_t color, Float_t alpha);
  void SetFillStyle(Style_t style);
  
  void SetMarkerColor(Color_t color=1);
  void SetMarkerColorAlpha(Color_t color, Float_t alpha);
  void SetMarkerSize(Size_t size=1);
  void SetMarkerStyle(Style_t style=1);
  
  //Styling for peak painter
  Color_t GetBaseLineColor() const;
  Style_t GetBaseLineStyle() const;
  Width_t GetBaseLineWidth() const;
  
  Color_t GetHitLineColor() const;
  Style_t GetHitLineStyle() const;
  Width_t GetHitLineWidth() const;
  
  void SetBaseLineColor(Color_t color);
  void SetBaseLineColorAlpha(Color_t color,Float_t alpha);
  void SetBaseLineStyle(Style_t style);
  void SetBaseLineWidth(Width_t width);
  
  void SetHitLineColor(Color_t color);
  void SetHitLineColorAlpha(Color_t color,Float_t alpha);
  void SetHitLineStyle(Style_t style);
  void SetHitLineWidth(Width_t width);
  
  Color_t GetPeakLineColorS(UInt_t peakidx) const;
  Color_t GetPeakLineColorE(UInt_t peakidx) const;
  Style_t GetPeakStyleS(UInt_t peakidx) const;
  Style_t GetPeakStyleE(UInt_t peakidx) const;
  Width_t GetPeakWidthS(UInt_t peakidx) const;
  Width_t GetPeakWidthE(UInt_t peakidx) const;
  
  Color_t GetFoundPeakLineColorS() const;
  Color_t GetFoundPeakLineColorE() const;
  Style_t GetFoundPeakStyleS() const;
  Style_t GetFoundPeakStyleE() const;
  Width_t GetFoundPeakWidthS() const;
  Width_t GetFoundPeakWidthE() const;
  
  void SetFoundPeakLineColorS(Color_t s_color);
  void SetFoundPeakLineColorE(Color_t e_color);
  void SetFoundPeakLineColor(Color_t s_color, Color_t e_color);
  void SetFoundPeakStyleS(Style_t s_style);
  void SetFoundPeakStyleE(Style_t e_style);
  void SetFoundPeakStyle(Style_t s_style,Style_t e_style);
  void SetFoundPeakWidthS(Width_t s_width);
  void SetFoundPeakWidthE(Width_t e_width);
  void SetFoundPeakWidth(Width_t s_width, Width_t e_width);
  
  void SetPeakLineColorS(UInt_t peakidx, Color_t s_color);
  void SetPeakLineColorE(UInt_t peakidx, Color_t e_color);
  void SetPeakLineColor(UInt_t peakidx, Color_t s_color, Color_t e_color);
  void SetPeakLineColorAlphaS(UInt_t peakidx, Color_t s_color,Float_t s_alpha);
  void SetPeakLineColorAlphaE(UInt_t peakidx, Color_t e_color,Float_t e_alpha);
  void SetPeakLineColorAlpha(UInt_t peakidx, Color_t s_color,Float_t s_alpha, Color_t e_color, Float_t e_alpha);
  void SetPeakLineStyleS(UInt_t peakidx, Style_t s_style);
  void SetPeakLineStyleE(UInt_t peakidx, Style_t e_style);
  void SetPeakLineStyle(UInt_t peakidx, Style_t s_style, Style_t e_style);
  void SetPeakLineWidthS(UInt_t peakidx, Width_t s_width);
  void SetPeakLineWidthE(UInt_t peakidx, Width_t e_width);
  void SetPeakLineWidth(UInt_t peakidx, Width_t s_width, Width_t e_width);
  
  void SetAllPeakLineColor(Color_t s_color, Color_t e_color);
  void SetAllPeakLineStyle(Style_t s_style, Color_t e_style);
  void SetAllPeakLineWidth(Width_t s_width, Width_t e_width);
  
protected:
  void Init();

  Int_t mComputedIndex;
  
  std::vector<PeakWindow> mPeaks; 
  
  PeakWindow mFoundPeak;
  
  PeakWindow mSearch;

  Double_t mBaseline;
  
  Double_t mBaselineSigma;
  
  Double_t mBaselineSigmaScale;
  
  Double_t mMaxX;
  Double_t mMaxY;

  Double_t mDeltaX;
  
  TGraph* mG_Data;
  
  //A bit redundant but used for internal checks
  Double_t mXRangeMin;  
  Double_t mXRangeMax;  
  Double_t mYRangeMin;  
  Double_t mYRangeMax;  
  
  Double_t mTunnelScale;      
  Double_t mTunnelSigma;      
  Double_t mTunnelThreshold;  
  
  Double_t mChiralityPeakScale; //how much to scale the slope of the line formed by the start and end points of the window
  Double_t mChiralityScale;     //modify whether peak position should be near start or end of peak window (see PeakWindow)
  Double_t mChiralityProbScale; //Scale for probablity 
  Double_t mChiralityThreshold; //Threshold chirality for peaks (if abs(chirality)>mChiralityThreshold merge +,- chirality)
  
  UInt_t mFilter;
  Double_t* mFilterWeights;  
  Int_t mFilterScale;        
  
  PeakAnaPainter* mPainter;  

  static double* GaussianMatrix2D(int rx,  double sx=0, int ry=0, double sy=0, bool kNorm=true);
  
  virtual bool MergeLeft(Double_t leftchir, Double_t thischir, Double_t rightchir) const;
  virtual std::vector< std::pair<int,int> > MergeIndices(std::vector<short>& vec) const;
  
  TString mOption;      
  
private:
  UInt_t mDEBUG;        
  bool mInternalSignal; 
  
  ClassDef(PeakAna,3);
};

#endif