StFcsWaveformFitMaker.h

// $Id: StFcsWaveformFitMaker.h,v 1.2 2021/05/30 21:26:53 akio Exp $
// $Log: StFcsWaveformFitMaker.h,v $
// Revision 1.2  2021/05/30 21:26:53  akio
// Added mFitDrawOn=2 for resetting mHitIdx for end of page, instead of each event
// Increased accepted tb range as hit, need further tuning
//
// Revision 1.1  2021/03/30 13:40:13  akio
// FCS code after peer review and moved from $CVSROOT/offline/upgrades/akio
//
// Revision 1.13  2021/02/25 21:56:28  akio
// Int_t -> int
//
// Revision 1.12  2021/02/25 19:24:33  akio
// Modified for STAR code review (Hongwei)
//
// Revision 1.11  2021/02/13 13:43:09  akio
// Some cleanups
//
// Revision 1.10  2021/02/13 03:11:44  akio
// Added TGraphAsymmErrorsWithReset and separate resetGraph() and simplify getGraph(idx=-1)
//
// Revision 1.9  2021/02/10 04:11:24  akio
// Code from David on improving TGraphArray housekeeping
//
// Revision 1.8  2021/01/26 18:43:14  akio
// Include David's fitter & holding of TGA for the event. Separate drawFit().
//
// Revision 1.7  2021/01/11 14:46:18  akio
// Adding another tail functional form, adding fit drawing to PDF
//
// Revision 1.6  2021/01/02 21:38:03  akio
// Fix some minor issues
//
// Revision 1.5  2021/01/02 21:04:54  akio
// added 11 for ped sub
//
// Revision 1.4  2020/09/17 19:01:52  akio
// fix bugs David reported
//
// Revision 1.3  2020/09/16 14:52:46  akio
// Use of TGraphAsymmErrors to deal with adc saturation
// Impelment gausFit() for gaussian fit
//
// Revision 1.2  2020/09/03 20:15:49  akio
// Adding res array for peak height/position/sigma and chi2
//
// Revision 1.1  2020/08/19 19:51:08  akio
// Initial version
//

#ifndef STROOT_STFCSWAVEFORMFITMAKER_STFCSWAVEFORMFITMAKER_H_
#define STROOT_STFCSWAVEFORMFITMAKER_STFCSWAVEFORMFITMAKER_H_

#include "StMaker.h"
#include "StFcsPulseAna.h"

#include "TH2.h"

class StFcsCollection;
class StFcsHit;
class StFcsDb;
class StFcsDbPulse;
class TGraphAsymmErrors;
class TGraph;
class TCanvas;

class StFcsWaveformFitMaker : public StMaker {
public:
    
    StFcsWaveformFitMaker(const char* name = "StFcsWaveformFitMaker");
    ~StFcsWaveformFitMaker();
    virtual int Init();
    virtual int InitRun(int runNumber);
    virtual int Make();
    virtual int Finish();
    virtual void Clear(Option_t* option = "");
    
    void setDebug(int v=1)        {SetDebug(v);}
    void setTest(int v);           
    void setEnergySelect(int ecal=13, int hcal=13, int pres=1) {mEnergySelect[0]=ecal; mEnergySelect[1]=hcal; mEnergySelect[2]=pres;}
    void setEnergySumScale(double ecal=1.0, double hcal=1.0, double pres=1.0) {mEnergySumScale[0]=ecal; mEnergySumScale[1]=hcal; mEnergySumScale[2]=pres;}
    void setCenterTimeBins(int v, int min=0, int max=512) {mCenterTB=v; mMinTB=min; mMaxTB=max;}
    void setAdcSaturation(int v)  {mAdcSaturation=(double)v;}
    void setError(double v)       {mError=v;}
    void setErrorSaturated(int v) {mErrorSaturated=v;}
    void setMinAdc(int v)         {mMinAdc=v;}
    void setTail(int v);                            
    void setMaxPeak(int v)        {mMaxPeak=v;}     
    void setPedTimeBins(int min, int max) {mPedMin=min; mPedMax=max;} 
    void setAnaWaveform(bool value=true){ mAnaWaveform=value; }   

    TGraphAsymmErrors* resetGraph();    
  
    TGraphAsymmErrors* getGraph(int idx=-1);
    TGraphAsymmErrors* getGraph(int det, int id);   

    //measuring fit time
    void setMeasureTime(char* file) {mMeasureTime=file;}
    TH1F* mTime = 0;

    //stage0 peak algo study
    TH2F* mTimeIntg[4];

    TGraphAsymmErrors* makeTGraphAsymmErrors(int n, double* t, double* adc);    
    TGraphAsymmErrors* makeTGraphAsymmErrors(TGraph* g);                 
    TGraphAsymmErrors* makeTGraphAsymmErrors(StFcsHit* hit);             
    TGraphAsymmErrors* makeTGraphAsymmErrors(TH1* hist);                 

    float analyzeWaveform(int select, TGraphAsymmErrors* g, float* res, TF1*& f, float ped=0.0); 
    float analyzeWaveform(int select, int n, double* t, double* adc, float* res, TF1*& f, float ped=0.0); 
    float analyzeWaveform(int select, TGraph* g, float* res, TF1*& f,float ped=0.0);                      
    float analyzeWaveform(int select, StFcsHit* hit, float* res, TF1*& f,float ped=0.0);                  

    //Pedestal analysis methods
    float AnaPed( TGraphAsymmErrors* g, float& ped, float& pedstd );
    float AnaPed( StFcsHit* hit, float& ped, float& pedstd );
    //actual analysis methods
    float sum8    (TGraphAsymmErrors* g, float *res);                             
    float sum16   (TGraphAsymmErrors* g, float *res);                             
    float highest (TGraphAsymmErrors* g, float *res);                             
    float highest3(TGraphAsymmErrors* g, float *res);                             
    float gausFit (TGraphAsymmErrors* g, float *res, TF1*& f, float ped=0.0);     
    float gausFitWithPed (TGraphAsymmErrors* g, float *res, TF1*& f);             
    float PulseFit1(TGraphAsymmErrors* g, float* res, TF1*& f, float ped=0.0);    
    float PulseFit2(TGraphAsymmErrors* g, float* res, TF1*& f, float ped=0.0);    
    float PulseFitAll(TGraphAsymmErrors* g, float* res, TF1*& f, float ped=0.0);  
    float PulseFit2WithPed( TGraphAsymmErrors* g, float* res, TF1*& f);           
    float PulseFitAllWithPed( TGraphAsymmErrors* g, float* res, TF1*& f);         
    float PedFitPulseFit( TGraphAsymmErrors* g, float* res, TF1*& f);             
  
    //Draw fits    
    void setMaxPage(int v){mMaxPage=v;}         
    void setSkip(int v){mSkip=v;}               

    void setFileName(char* file, int maxpage=25, int skip=5){mFilename=file; mMaxPage=maxpage; mSkip=skip;} 
    void writeFile(std::string filename);       

    void setFitDrawOn(int v=1) {mFitDrawOn=v;}
    void setFitFilter(char* filter) {mFilter=filter; mFitDrawOn=2;}

    //Draw from David
    int centerTB()const{return mCenterTB;}                                
    void setDavidFitter(StFcsPulseAna* v){if(mPulseFit==0){mPulseFit=v;}} 
    StFcsPulseAna* davidFitter(){return mPulseFit;}                       
    StFcsPulseAna* InitFitter(Double_t ped=0);     

    static int GenericPadPos(int value, int Nvals, int PadNums );
    static int PadNum4x4(int det, int col, int row);
    void drawRegion(int det, int col_low, int row_low, int col_high, int row_high, int event=0);
    void drawEvent(int det, int event=0);    
    void drawFitter(Option_t* opt){ if(mPulseFit!=0){mPulseFit->Draw(opt);} }  

    void drawCh(UInt_t detid, UInt_t ch) const;
    void drawDualFit(UInt_t detid, UInt_t ch);

    void printArray() const;                 
    void printSetup() const;                 

 protected:
    TClonesArray mChWaveData;  
    void drawFit(TGraphAsymmErrors* gg, TF1* func); 
    StFcsPulseAna* mPulseFit;   

    int mTest = 0;

    TFile* mOutFile;  
    //For testing Dep0 algo (mTest==1)
    TH2F* mH2_Dep0DepMod[3];
    TH2F* mH2_Sum8Dep0[3];
    TH2F* mH2_Sum8DepMod[3];
    //For testing number of peaks finding (mTest==2||mTest==3||mTest==6||mTest==8)
    TH1F* mH1_NPeaksAkio = 0;                       
    TH1F* mH1_NPeaksFilteredAkio = 0;               
    TH2F* mH2F_AdcTbAkio[6];                        
    TH2F* mH2F_AdcTbMine[6];                        
    TH2F* mH2F_AdcTbValidPeak[7];                   
    TH2F* mH2F_SumFitvSumWin[6];                    
    TH2F* mH2F_APeakvMPeak[6];                      
    TH1F* mH1F_PeakStart[6];                        
    TH1F* mH1F_PeakEnd[6];                          
    TH1F* mH1_PeakTiming = 0;                       

    TH1F* mH1F_NPeaks[7];                           
    TH1F* mH1F_NPeaksFiltered[7];                   
    TH2F* mH2_NPeakvsPeakXdiff = 0;                 
    TH2F* mH2_NPeakvsPeakYratio = 0;                
    TH1F* mH1_VOverlap = 0;                         
    TH2F* mH2_NOvsNPeaks = 0;                       
    TH2F* mH2_VvsNOverlap = 0;                      
    TH2F* mH2F_NOvsId[6];                           
    TH1F* mH1F_Res0[7];                             
    TH1F* mH1F_Res0Zoom[7];                         
    TH1F* mH1F_Sum8Res0[7];                         
    TH1F* mH1F_Sum8Res0Zoom[7];                     
    TH1F* mH1F_FitRes0[7];                          
    TH1F* mH1F_FitRes0Zoom[7];                      
    TH2F* mH2F_Sum8vFit[7];                         
    TH1F* mH1_TimeFitPulse = 0;                     

    //For testing peak height vs. sigma
    TH2F* mH2_HeightvsSigma = 0;                    
    TH2F* mH2_HeightvsSigmaTrig = 0;                
    TH1F* mH1_ChiNdf = 0;                           
    TH2F* mH2_HeightvsChiNdf = 0;                   
    TH2F* mH2_MeanvsChiNdf = 0;                     
    TH2F* mH2_SigmavsChiNdf = 0;                    

    TH1F* mH1_PeakTimingGaus = 0;                   
    TH1F* mH1_PeakTimingPuls = 0;                   
    TH2F* mH2_PeakTimingCompare = 0;                

    void SetupDavidFitterMay2022(Double_t ped=0);   

    int PeakCompare(const PeakWindow& pwin1, const PeakWindow& pwin2 );
    std::vector<int> NPeaksPrePost(int& trigidx, Double_t& xmin, Double_t& xmax) const;

 private:
    StFcsDb* mDb=0;                    
    StFcsDbPulse* mDbPulse = 0;        
    StFcsCollection* mFcsCollection=0; 

    unsigned int mHitIdx = 0; 

    //Figures out error to set on the TGraphAsymmErrors for a given point and adc
    char *mMeasureTime=0;                

    int mEnergySelect[3];                
    double mEnergySumScale[3];           
    bool mAnaWaveform;                   
    int mCenterTB=50;                    
    int mMinTB=0;                        
    int mMaxTB=512;                      

    double mError=1.0;                   
    double mErrorSaturated=1000.0;       
    double mAdcSaturation=4000.0;        
 
    int mPedMin=-1;                      
    int mPedMax=-1;                      

    //gausFit
    int mMinAdc=5;                       
    int mTail=2;                         
    int mMaxPeak=5;                      

    //Drawing fits
    TCanvas* mCanvas=0;                  
    int mPage=0;                         
    int mPad=0;                          
    int mMaxPage=25;                     
    int mSkip=2;                         
    char* mFilename=0;                   
    char* mFilter=0;
    char mDetName[100];                  
    int mFitDrawOn=0;                    

    virtual const Char_t *GetCVS() const {static const Char_t cvs[]="Tag " __DATE__ " " __TIME__ ; return cvs;}

    ClassDef(StFcsWaveformFitMaker, 4)
};
#endif  // STROOT_STFCSWAVEFORMFITMAKER_STFCSWAVEFORMFITMAKER_H_