diff --git a/StRoot/StFcsDbMaker/StFcsDb.cxx b/StRoot/StFcsDbMaker/StFcsDb.cxx index bcd98ae3349..b666a7c2f09 100644 --- a/StRoot/StFcsDbMaker/StFcsDb.cxx +++ b/StRoot/StFcsDbMaker/StFcsDb.cxx @@ -523,6 +523,15 @@ int StFcsDb::getDetFromName(const std::string& detname){ } } +unsigned short StFcsDb::getKey(unsigned short detid, unsigned short id){ return ( (detid & 0x7)<<12 | (id & 0xffff) ); } +void StFcsDb::getDetIdFromKey(unsigned short key, unsigned short& detid, unsigned short& id) +{ + detid = (key >> 12) & 0x0007; + id = (key & 0x0fff); +} +unsigned short StFcsDb::getDetFromKey(unsigned short key){ return (key >> 12) & 0x0007; } +unsigned short StFcsDb::getIdFromKey(unsigned short key) { return (key & 0x0fff); } + StThreeVectorD StFcsDb::getDetectorOffset(int det) const{ if(mRun19>0){ const float bOffY=-(17.0*5.81); //40in=101.6cm and 17*5.81=98.76 so I will leave this unchanged @@ -2048,4 +2057,4 @@ const g2t_track_st* StFcsDb::getG2tTrack(StFcsCluster* c, g2t_track_st* g2ttrk, } } return &g2ttrk[parents[order].first-1]; -} \ No newline at end of file +} diff --git a/StRoot/StFcsDbMaker/StFcsDb.h b/StRoot/StFcsDbMaker/StFcsDb.h index c851d24704b..19e81e3e59f 100644 --- a/StRoot/StFcsDbMaker/StFcsDb.h +++ b/StRoot/StFcsDbMaker/StFcsDb.h @@ -151,6 +151,10 @@ class StFcsDb : public TDataSet { void getName(int ehp, int ns, int dep, int ch, char name[]); //! Get Name of a channel static void getFromName(const char name[], int& det, int& id); //! Get det/id from name static int getDetFromName(const std::string& detname); //! Get det from name + static unsigned short getKey(unsigned short detid, unsigned short id);//This key matches the 'mDetId' in StFcsHit but without "zs" value (@[June 21, 2022](David Kapukchyan)>Move static method to StFcsHit?) + static void getDetIdFromKey(unsigned short key, unsigned short& detid, unsigned short& id); + static unsigned short getDetFromKey(unsigned short key); + static unsigned short getIdFromKey(unsigned short key); //! Utility functions related to DetectorPosition StThreeVectorD getDetectorOffset(int det) const; //! get the offset of the detector diff --git a/StRoot/StFcsDbMaker/StFcsDbMaker.cxx b/StRoot/StFcsDbMaker/StFcsDbMaker.cxx index e65c63e0ddf..7478d95a523 100644 --- a/StRoot/StFcsDbMaker/StFcsDbMaker.cxx +++ b/StRoot/StFcsDbMaker/StFcsDbMaker.cxx @@ -118,6 +118,7 @@ #include "StFcsDbMaker.h" #include "StFcsDb.h" +#include "StFcsDbPulse.h" #include "St_db_Maker/St_db_Maker.h" #include "StMessMgr.h" @@ -136,14 +137,18 @@ StFcsDbMaker::StFcsDbMaker(const char *name) : StMaker(name){ LOG_INFO << "******** StFcsDbMaker::StFcsDbMaker = "<Init(); + mFcsDbPulse->Init(); return kStOK; return StMaker::Init(); } diff --git a/StRoot/StFcsDbMaker/StFcsDbMaker.h b/StRoot/StFcsDbMaker/StFcsDbMaker.h index b592a5b5b7e..77de38eddbd 100644 --- a/StRoot/StFcsDbMaker/StFcsDbMaker.h +++ b/StRoot/StFcsDbMaker/StFcsDbMaker.h @@ -93,6 +93,7 @@ #endif class StFcsDb; +class StFcsDbPulse; class StFcsDbMaker : public StMaker { public: @@ -107,9 +108,10 @@ class StFcsDbMaker : public StMaker { private: StFcsDb *mFcsDb; int mDbAccess=1; + StFcsDbPulse* mFcsDbPulse; virtual const Char_t *GetCVS() const {static const Char_t cvs[]="Tag " __DATE__ " " __TIME__ ; return cvs;} - ClassDef(StFcsDbMaker,4) //StAF chain virtual base class for Makers + ClassDef(StFcsDbMaker,0) //StAF chain virtual base class for Makers }; #endif diff --git a/StRoot/StFcsDbMaker/StFcsDbPulse.cxx b/StRoot/StFcsDbMaker/StFcsDbPulse.cxx new file mode 100644 index 00000000000..31f5fecc1fc --- /dev/null +++ b/StRoot/StFcsDbMaker/StFcsDbPulse.cxx @@ -0,0 +1,186 @@ +#include "StFcsDbPulse.h" + +ClassImp(StFcsDbPulse); + +StFcsDbPulse::StFcsDbPulse(const char* name) : TDataSet(name) +{ + mTail = 0; + mTBPerRC = 8; + + mGSigma = 0; + mA1 = 0; + mA2 = 0; + mXoff1 = 0; + mXoff2 = 0; + mTau1 = 0; + mTau2 = 0; + mP1 = 0; + mP2 = 0; +} + +StFcsDbPulse::~StFcsDbPulse(){} + +int StFcsDbPulse::Init() +{ + + mGSigma = 24.5; // signal nominal sigma + //GSigmaSigma = 0; // distibution of signal sigma + mA1 = 1.0/sqrtpi()/0.155/129; + mA2 = 0.2/sqrtpi()/0.155/129; + mXoff1 = 70 - 129; + mXoff2 = 220 - 129; + mTau1 = 200.0; + mTau2 = 40.0; + mP1 = 1.0; + mP2 = 1.0; + + return 0; +} + +void StFcsDbPulse::setTail(int tail) +{ + mTail = tail; + if( tail == 0 ){return;} + if( tail == 1 ) + { + //Data from Gerard 2020 summer + mGSigma = 24.5/nsecPerTB(); + mA1 = 1.0/0.154/mGSigma; + mA2 = 0.2/0.154/mGSigma; + mXoff1 = (70 - 129)/nsecPerTB(); + mXoff2 = (220 - 129)/nsecPerTB(); + mTau1 = 200.0/nsecPerTB(); + mTau2 = 40.0/nsecPerTB(); + mP1 = 1.0; + mP2 = 1.0; + } + else if( tail == 2 ) + { + //Data from WAH with real detector/LED system 2021 Jan + mGSigma = 2.347; + mA1 = 2543.0/854.0/mGSigma; + mA2 = 0.0; + mXoff1 = 211.3-215.7; + mXoff2 = 0.0; + mTau1 = 4.375; + mTau2 = 0.0; + mP1 = 1.0; + mP2 = 0.0; + } + else{ LOG_WARN << "StFcsDbPulse::setTail - Invalid Tail value " << tail << endm; } +} + +void StFcsDbPulse::setTGraphAsymmErrors(TGraphAsymmErrors* gae, const int &i, const double &adc, double Yerr, double YerrSat){ + if( gae == 0 ){LOG_ERROR << "StFcsDbPulse::setTGraphAsymmErrors - Graph object cannot be zero" << endm; return;} + double HighY = 0; + if(adcSetPointError(i,0,0,Yerr,HighY); +} + +//this one is just shower shape function +double StFcsDbPulse::pulseShape(double* x, double* p) { + double ret = p[0]*exp(-0.5*pow((x[0]-p[1])/p[2],2)); + if(mTail>0){ + double x1 = x[0] - p[1] - mXoff1; + if(x1>0){ + double a0 = p[0] * p[2]; + ret += a0*mA1/mTau1/mTau1*pow(x1,mP1)*exp(-x1/mTau1); + if(mA2>0){ + double x2 = x[0] - p[1] - mXoff2; + if(x2>0){ + ret += a0*mA2/mTau2/mTau2*pow(x2,mP2)*exp(-x2/mTau2); + } + } + } + } + return ret; +} + +double StFcsDbPulse::multiPulseShape(double* x, double* p) { + int npulse = p[0]; + double ret = p[1]; + for(int i=0; i=xhigh ){LOG_ERROR << "StFcsDbPulse::createPulse - Invalid range" <(value+(Nvals*PadNums))/static_cast(Nvals) );} + else{ return GenericPadPos(value-(Nvals*PadNums), Nvals, PadNums); } +} + +int StFcsDbPulse::PadNum4x4(int det, int col, int row) +{ + int ncol = 0; + int nrow = 0; + if( det<=1 ) + { + ncol = 2; + nrow = 3; + } + else if( det<=3 ) + { + ncol = 2; + nrow = 2; + } + else{ LOG_ERROR << "This only works for Ecal and Hcal" << endm; return 0;} + int padcol = GenericPadPos(col,ncol,4); + int padrow = GenericPadPos(row,nrow,4); + return 4*(padrow-1)+padcol; +} + +Int_t StFcsDbPulse::getYMinMax(TGraphAsymmErrors* gae, Double_t &Ymin, Double_t &Ymax, Double_t xmin, Double_t xmax) +{ + if( gae==0 ){ return -1; } + Int_t index = -1; + //Start with invalid values + Double_t MinY = Ymax; + Double_t MaxY = Ymin; + for( int i=0; iGetN(); ++i ) + { + Double_t X; Double_t Y; + gae->GetPoint(i,X,Y); + if( Xxmax ){continue;} + if( Y>MaxY ){ MaxY=Y; index=i;} + if( Y + +#include "TDataSet.h" +#include "TMath.h" +#include "TGraphAsymmErrors.h" +#include "TF1.h" +#include "St_base/StMessMgr.h" + +class StFcsDbPulse : public TDataSet { + + public: + StFcsDbPulse(const char* name = "fcsPulse"); //!< Constructor + virtual ~StFcsDbPulse(); //!< Destructor + int Init(); //!< Initialize object + + //Basic contants + static double sqrtpi(){return sqrt(TMath::Pi());} //!< sqrt(TMath::Pi) + static double sqrt2pi(){return sqrt(2.0*TMath::Pi());} //!< sqrt(2*TMath::Pi) + + //Constants for simulating pulses + void setTBPerRC(double v){mTBPerRC = v;} //!< @param v set #mTBPerRC + double TBPerRC()const{return mTBPerRC;} //!< @return #mTBPerRC + double nsecPerTB()const{return 107.0/TBPerRC();} //!< nanoseconds per timebin + double BeamLengthSig()const{return 10.0/nsecPerTB();} //!< beam length sigma + + /**@brief Sets the variables needed by the sum of xexp functions that describe the tail of the pulse shape + + It will set the constants related to the tail of the pulse shape to known values based on previous data analysis. + @param tail set values for pulse shape tail based on predefined values\n + 0 = no tail\n + 1 = tail shape based on LED data from Lab\n + 2 = tail shape based on Run 22 LED data + */ + void setTail(int tail); + + void setGSigma(double v){mGSigma = v;} //!< @param v set #mGSigma + void setA1(double v){mA1 = v;} //!< @param v set #mA1 + void setA2(double v){mA2 = v;} //!< @param v set #mA2 + void setXoff1(double v){mXoff1 = v;} //!< @param v set #mXoff1 + void setXoff2(double v){mXoff2 = v;} //!< @param v set #mXoff2 + void setTau1(double v){mTau1 = v;} //!< @param v set #mTau1 + void setTau2(double v){mTau2 = v;} //!< @param v set #mTau2 + void setP1(double v){mP1 = v;} //!< @param v set #mP1 + void setP2(double v){mP2 = v;} //!< @param v set #mP2 + + //Variables related to the tail function + double GSigma()const{return mGSigma;} //!< @return #mGSigma + double A1()const{return mA1;} //!< @return #mA1 + double A2()const{return mA2;} //!< @return #mA2 + double Xoff1()const{return mXoff1;} //!< @return #mXoff1 + double Xoff2()const{return mXoff2;} //!< @return #mXoff2 + double Tau1()const{return mTau1;} //!< @return #mTau1 + double Tau2()const{return mTau2;} //!< @return #mTau2 + double P1()const{return mP1;} //!< @return #mP1 + double P2()const{return mP2;} //!< @return #mP2 + + /**@brief Figure out and set the errors on FCS pulse data stored in a TGraphAsymmErrors object + + Assign the error on each point of the TGraphAsymmErrors object that either holds the pulse data or is being filled with the pulse data. The error on X is assumed to be zero. ADC saturation is fixed at 4K because it is 12 bit. + @param gae graph object that holds the pulse data + @param i point on the graph object to modify + @param adc ADC value to be checked for error assignment + @param Yerr Y error on ADC for values below ADC saturation + @param YerrSat Y error on ADC for values above ADC saturation + */ + static void setTGraphAsymmErrors(TGraphAsymmErrors* gae, const int &i, const double &adc, double Yerr, double YerrSat); + + /**@brief Single pulse shape gaus+xexp+xexp + + 1-dimensional function to describe the pulse shape which is a Gaussian for the main part plus a tail\n + gaus = Gaussian\n + xexp = A/tau/tau*(x-x0)^n*exp(-(x-x0)/tau))\n + parameters of xexp are fixed constants in this class and describe the tail of the pulse shape + @param x one dimension array as x-value input of function + @param p array of parameters for Gaussian + */ + double pulseShape(double* x, double* p); + + /**@brief Multi-pulse shape function constant+gaus+xexp+xexp for many pulses + + Extends #pulseShape() to be able to generate many pulses. The values associated with xexp are fixed for all pulses with the constants defined in this class + @param x one dimension array as x-value input of function + @param p parameters to describe number of pulses and the Gaussian for those pulses\n + p[0] = number of pulses\n + p[1] = constant offset\n + p[2] = height of Gaussian pulse 1\n + p[3] = mean of Gaussian pulse 1\n + p[4] = sigma of Gaussian pulse 1\n + p[5] = height of Gaussian pulse 2\n + ...\n + */ + double multiPulseShape(double* x, double* p); + + /**@brief Function to create pulse shape for FCS, 5 parameters is minimum + + Creates a ROOT TF1 for #multipulseShape() + @param xlow lowest x-value for function + @param xhigh highest x-value for function + @param npars number of parameters needed by TF1 (2+3*number of pulses) + */ + TF1* createPulse(double xlow=0, double xhigh=1,int npars=5 ); + + /**@brief Function to tell you pad number when drawing multiple objects on the same pad + + This function is mostly used by #PadNum4x4() for drawing FCS Ecal and Hcal channels on the same canvas + @param value is the row or column number of an Ecal or Hcal channel + @param Nvals is the number of stuff per column or row to put on same pad + @param PadNums is the total number of pads in the column or row + */ + static int GenericPadPos(int value, int Nvals, int PadNums ); + + /**@brief Function that gives pad number when drawing a specific detector id + + This function is used to know where to draw a specific Ecal or Hcal channels on utilizing a canvas with 4 rows and 4 columns + @param det Ecal or Hcal detector ID + @param col column number of Ecal or Hcal channel to draw + @param row row number of Ecal or Hcal channel to draw + */ + static int PadNum4x4(int det, int col, int row); + + /**@brief Finds minimum and maximum y-values in a TGraph and returns index for max y + + This function is used to scan a TGraph and find the global minimum and maximum y-values inside a specified x-range. + @param gae TGraph to scan + @param Ymin lowest possible y-value to check and then store result of found minimum y-value + @param Ymax highest possible y-value to check and then store result of found maximum y-value + @param xmin minimum x-value to start searching for y-min/max + @param xmax minimum x-value to start seraching for y-min/max + */ + static Int_t getYMinMax(TGraphAsymmErrors* gae, Double_t &Ymin, Double_t &Ymax, Double_t xmin=-5, Double_t xmax=2000); + + virtual void Print(Option_t* opt = "") const; //!< Print all the constants associated with this class + + protected: + double mTBPerRC; //!< number of timebins in one RHIC crossing + double mGSigma; //!< pulse shape nominal sigma of Gaussian part + //Pulse shape constants + double mA1; //!< pulse shape tail: height of first xexp function + double mA2; //!< pulse shape tail: height of second xexp function + double mXoff1; //!< pulse shape tail: x offset of first xexp function + double mXoff2; //!< pulse shape tail: x offset of second xexp function + double mTau1; //!< pulse shape tail: scale of first xexp function + double mTau2; //!< pulse shape tail: scale of second xexp function + double mP1; //!< pulse shape tail: power of first xexp function + double mP2; //!< pulse shape tail: power of second xexp function + + private: + int mTail; //!< pulse shape tail, turn off or on to known values + static const int mAdcSaturation = 4000; //!< ADC is 12 bit number so it saturates at 4K + + ClassDef(StFcsDbPulse, 1); +}; + +#endif diff --git a/StRoot/StFcsWaveformFitMaker/PeakAna.cxx b/StRoot/StFcsWaveformFitMaker/PeakAna.cxx new file mode 100644 index 00000000000..3c56ee6291d --- /dev/null +++ b/StRoot/StFcsWaveformFitMaker/PeakAna.cxx @@ -0,0 +1,1280 @@ +#include "PeakAna.h" +#include "PeakAnaPainter.h" + +ClassImp(PeakAna) + +PeakAna::PeakAna() +{ + mG_Data = 0; + Init(); +} + +PeakAna::PeakAna(TGraph* data ) +{ + mG_Data = data;//Set graph object now so doesn't get created in 'Init' + Init();//This sets internal signal to true so set it false after call + mInternalSignal=false; +} + +PeakAna::PeakAna( int size, double* xvals, double* yvals ) +{ + mG_Data = new TGraph(size,xvals,yvals); + Init(); +} + +PeakAna::PeakAna( TH1* hist ) +{ + mG_Data = ConvertHistToGraph(hist); + Init(); +} + +void PeakAna::Init() +{ + mDEBUG = 0; + mSearch.mStartX = 0; mSearch.mEndX=1; + + mComputedIndex = -1;//This will keep from having to determine peaks every time + mInternalSignal = true; + mBaseline = -5.0; + mBaselineSigma = 0.75; + mBaselineSigmaScale = 4.0; + + mMaxX = -5.0; + mMaxY = -5.0; + + if( mG_Data==0 ){ mG_Data = new TGraph(); } + + mXRangeMin = 0; + mXRangeMax = 1; + mYRangeMin = 0; + mYRangeMax = 1; + mFoundPeak.SetWindow(mXRangeMin,mXRangeMax); + + mTunnelScale = 1;//1 is default + mTunnelSigma = 1;//1 is default + mTunnelThreshold = -1;//Turn off by default + + mChiralityPeakScale = 1;//default is scale by 1 + mChiralityScale = 1;//default is peaks must be centered in window + mChiralityProbScale = 1; + mChiralityThreshold = -1;//Turn off by default + + mDeltaX = -1.0;//Turn off by default + mFilter = 0; + mFilterScale = 1; + mFilterWeights = 0; + + mPainter = 0; +} + +PeakAna::PeakAna(const PeakAna &OldAna,TGraph* graph) +{ + //Note that the graph is not cloned because PeakAna should not be changing the contents of the data + if( graph!=0 ){ mG_Data = graph; }//Set graph object now so doesn't get created in 'Init' + else{ mG_Data = OldAna.GetData(); }//If no graph object given then clone the graph + if( mG_Data!=0 ){ mG_Data = (TGraph*)mG_Data->Clone(); }//In case 'OldAna' graph is zero then don't call clone, since if it is zero Init() will create a new graph object + Init(); + if( graph!=0 ){mInternalSignal=false;}//Init() sets internal signal to true so set it false after call if external graph object given + + ((PeakAna&)OldAna).Copy(*this); + +} + +PeakAna& PeakAna::operator=(const PeakAna& rhs) +{ + if( this == &rhs ){ return *this; } + //Note that the graph is not cloned because PeakAna should not be changing the contents of the data + mG_Data = rhs.GetData();//Set graph object now so doesn't get created in 'Init' + Init();//This sets internal signal to true so set it false after call + mInternalSignal=false;//Since graph is not copied don't treat as an internal signal + + ((PeakAna&)rhs).Copy(*this); + + return *this; +} + +//Note that this function doesn't copy the graph object or a copies the vector of peaks. This is done on purpose becuase the idea is that the copy should run its own analysis and the algorithm is robust enough to generate the same results. To copy all objects use TObject::Clone +void PeakAna::Copy(TObject& obj) const +{ + TObject::Copy(obj);//Copy TObject bits + + ((PeakAna&)obj).mMaxX = mMaxX; + ((PeakAna&)obj).mMaxY = mMaxY; + + ((PeakAna&)obj).mXRangeMin = mXRangeMin; + ((PeakAna&)obj).mYRangeMin = mYRangeMin; + ((PeakAna&)obj).mXRangeMax = mXRangeMax; + ((PeakAna&)obj).mYRangeMax = mYRangeMax; + ((PeakAna&)obj).mBaseline = mBaseline; + ((PeakAna&)obj).mBaselineSigma = mBaselineSigma; + ((PeakAna&)obj).mBaselineSigmaScale = mBaselineSigmaScale; + ((PeakAna&)obj).mSearch.SetWindow(SearchPeak(),SearchWidth()); + ((PeakAna&)obj).mDeltaX = mDeltaX; + ((PeakAna&)obj).mFilter = mFilter; + ((PeakAna&)obj).mFilterScale = mFilterScale; + //((PeakAna&)obj).mFilterSigma = mFilterSigma; + if( ((PeakAna&)obj).mFilterWeights!=0 ){ delete [] ((PeakAna&)obj).mFilterWeights; ((PeakAna&)obj).mFilterWeights = 0; } + if( mFilterWeights!=0 && mFilterScale>0 ){ + ((PeakAna&)obj).mFilterWeights = new double[2*mFilterScale+1]; + std::copy( mFilterWeights,mFilterWeights+(2*mFilterScale+1), ((PeakAna&)obj).mFilterWeights ); + } + else{ ((PeakAna&)obj).mFilterWeights = 0; } + ((PeakAna&)obj).mTunnelScale = mTunnelScale; + ((PeakAna&)obj).mTunnelSigma = mTunnelSigma; + ((PeakAna&)obj).mTunnelThreshold = mTunnelThreshold; + + ((PeakAna&)obj).mPainter = 0;//Dont' copy painter +} + +TObject* PeakAna::Clone(const char* newname) const +{ + TGraph* cloneg = (TGraph*)GetData()->Clone(newname); + PeakAna* ana = new PeakAna(cloneg); + Copy(*ana); + ana->ForceInternal();//Cloned graphs should always be internal + if( FoundPeakIndex()>=0 ){ + //If mComputedIndex>=0 this means AnalyzeForPeak() was called and need to copy the peaks. This is preffered over re-running the analysis as copy should be faster than re-running the data + ana->mComputedIndex = mComputedIndex; + for( UInt_t i=0; imPeaks.push_back( mPeaks.at(i) ); + } + ana->mFoundPeak = mFoundPeak; + } + return ana; +} + +PeakAna::~PeakAna() +{ + if(mInternalSignal ){delete mG_Data;} + delete [] mFilterWeights; + delete mPainter; +} + +Int_t PeakAna::AnalyzeForPeak() +{ + ResetPeak(); + if( mFilter == 1 ){ MeanFilter(mFilterScale,false); } + if( mFilter == 2 ){ GausFilter(mFilterScale,false); } + this->GetPossiblePeaks(); + //Doesn't make sense to merge for number of peaks <= 1 + if( mChiralityThreshold>=0 && mPeaks.size()>1 ){ + std::vector mergedpeaks; + this->MergeByChirality(mergedpeaks); + mPeaks.swap( mergedpeaks ); + }; + mComputedIndex = this->SearchForPeak( mPeaks );//Returns a valid index for the possible peak vector or the size of the vector for an invalid peak + if( mComputedIndex(mPeaks.size()) ){ mFoundPeak = mPeaks.at(mComputedIndex); } + else{ + mFoundPeak.SetWindow(mXRangeMax+1,mXRangeMax+1);//If no valid index then set found peak to greater than max X values so it is obvious something's wrong + } + //@[May 10, 2022]>mFoundPeak is a copy of the one in the vector. When drawing the windows the one in the vector is accessed not 'mFoundPeak' which is why the color for the one in the vector is changed and not mFoundPeak. + return mComputedIndex; +} + +Int_t PeakAna::AnalyzeForPeak(Double_t peak, Double_t width) +{ + SetSearchWindow(peak,width); + return this->AnalyzeForPeak(); +} + +Int_t PeakAna::AnalyzeForNoisyPeak() +{ + ResetPeak(); + this->GetPossiblePeaks(); + PeakAna noisyana = PeakAna::ConvertPeaksToAna(*this); + noisyana.SetTunnelThreshold(-1);//Don't use peak tunneling for this kind of peak finding + noisyana.AnalyzeForPeak(); + //Replace old peak values with new one (maybe use swap in future??) + mPeaks.clear(); + for( Int_t inoisy=0; inoisySearchForPeak( mPeaks );//Returns a valid index for the possible peak vector or the size of the vector for an invalid peak + if( mComputedIndex(mPeaks.size()) ){ + mFoundPeak = mPeaks.at(mComputedIndex); //No offset for now rethink in the future?? + //@[May 10, 2022]>mFoundPeak is a copy of the one in the vector. When drawing the windows the one in the vector is accessed not 'mFoundPeak' which is why the color for the one in the vector is changed and not mFoundPeak. + //mPeaks.at(mComputedIndex).SetStartLineColor(kRed); + //mPeaks.at(mComputedIndex).SetEndLineColor(kOrange); + } + else{mFoundPeak.SetWindow(mXRangeMax+1,mXRangeMax+1);}//If no valid index then set found peak to greater than max X values so it is obvious something's wrong + return mComputedIndex; +} + +bool PeakAna::ValidPeakIdx() const +{ + if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){return true;} + else{return false;} +} + +bool PeakAna::GoodWindow() +{ + if( mComputedIndex<0 ){this->AnalyzeForPeak();} + //First check if start and end times are within our max timebin window of 0-1023 + if( mFoundPeak.mStartXmXRangeMax ){return false;} + else if( mFoundPeak.mEndX mXRangeMax ){return false;} + //Next check if values make physical sense + else if( mFoundPeak.mStartX==mFoundPeak.mEndX ){return false;} + else if( mFoundPeak.mStartX > mFoundPeak.mEndX ){return false;} + else{return true;} +} + +void PeakAna::GetXYMax(Double_t xmin, Double_t xmax) +{ + if( mG_Data==0 ){ return; } + for( int i=0; iGetN(); ++i ){ + Double_t X; Double_t Y; + mG_Data->GetPoint(i,X,Y); + if( Xxmax ){continue;} + if( Y>mMaxY ){ mMaxY=Y; mMaxX=X; } + } + if( mMaxY<0 && mMaxX<0 ){ + LOG_WARN << "Unable to find a maximum adc\nSetting to impossible values" << endm; + mMaxY = 5000; + mMaxX = 5000; + } +} + +Double_t PeakAna::PeakX() +{ + if( mComputedIndex<0 ){ this->AnalyzeForPeak(); } + if( mFoundPeak.mP_Peak<0){ return 0; } + else{ Double_t x,y; mG_Data->GetPoint(mFoundPeak.mP_Peak,x,y); return x;} +} + +Double_t PeakAna::PeakY() +{ + if( mComputedIndex<0 ){this->AnalyzeForPeak();} + if( mFoundPeak.mP_Peak<0){ return 0; } + else{ Double_t x,y; mG_Data->GetPoint(mFoundPeak.mP_Peak,x,y); return y;} +} + +void PeakAna::PeakXY(Double_t &xval, Double_t &yval) +{ + if( mComputedIndex<0 ){this->AnalyzeForPeak();} + if( mFoundPeak.mP_Peak<0 ){ xval=0; yval=0; } + else{ mG_Data->GetPoint(mFoundPeak.mP_Peak,xval,yval); } +} + +void PeakAna::Print(Option_t* opt) const +{ + std::cout << "|P::Graph:"<MergeLeftOrRight(i) + << std::endl; + } +} + +void PeakAna::ResetPeak() +{ + //If computed index is negative then analyzeforpeak was never called and it is not needed to reset anything + if( mComputedIndex>=0 ){ + mComputedIndex = -1; + mPeaks.clear(); + mFoundPeak.Reset(mXRangeMin,mXRangeMax); + } +} + +void PeakAna::SetData(TGraph* graph) +{ + if( graph==0 ){LOG_ERROR << "PeaAna::SetData - Graph cannot be 0!" << endm; return;} + ResetPeak(); + if( mInternalSignal && mG_Data!=0 ){delete mG_Data;} + mG_Data = graph; + mInternalSignal = false; +} + +void PeakAna::SetData(TH1* hist, UInt_t numavgs) +{ + if( hist==0 ){LOG_ERROR << "PeakAna::SetData - Histogram cannot be 0" << endm; return;} + ResetPeak(); + if( mInternalSignal ){delete mG_Data;} + mG_Data = PeakAna::ConvertHistToGraph(hist,numavgs); + mInternalSignal = true; +} + +void PeakAna::SetFilter( UInt_t filter, Int_t scale, Double_t sigma ) +{ + mFilter = filter; + mFilterScale=scale; + if( scale<=0 ){ return; } + if( sigma==0 ){ sigma=scale/2.0; }//if sigma is zero use scale/2 + if( mFilterWeights!=0 ){ delete [] mFilterWeights; mFilterWeights=0; } + if( filter==2 ){ mFilterWeights = GaussianMatrix2D(scale,sigma); } +} + +void PeakAna::SetBaseline(Double_t value, Double_t sigma) +{ + if( sigma<0){ LOG_WARN << "PeakAna::SetBaseline - Baseline sigma should not be less than zero\n" << endm; } + mBaseline = value; + mBaselineSigma = sigma; +} + +void PeakAna::SetBaselineSigmaScale(Double_t scale) +{ + mBaselineSigmaScale = fabs(scale); +} + +void PeakAna::SetRange( Double_t xmin, Double_t ymin, Double_t xmax, Double_t ymax) +{ + mXRangeMin = xmin; + mYRangeMin = ymin; + mXRangeMax = xmax; + mYRangeMax = ymax; +} +void PeakAna::SetSearchWindow(Double_t peak, Double_t width) +{ + mSearch.SetWindow(peak,width); + ResetPeak();//reset since old found peak should be invalid +} + +void PeakAna::SetPeak(const Int_t peakpoint, const Double_t peakx ) +{ mFoundPeak.mP_Peak=peakpoint; mFoundPeak.mPeakX=peakx; } +void PeakAna::SetWindow(const Int_t start, const Int_t end ) +{ mFoundPeak.SetWindow(start,end); } + +void PeakAna::SetWindow(PeakWindow window) +{ mFoundPeak=window; } + +void PeakAna::SetTunnelScale(Double_t value) +{ + if( value>=0 ){mTunnelScale = value;} +} +void PeakAna::SetTunnelSigma(Double_t value) +{ + if( value>0 ){mTunnelSigma = value;} +} +void PeakAna::SetTunnelPars(Double_t scale, Double_t sigma) +{ + SetTunnelScale(scale); + SetTunnelSigma(sigma); +} +void PeakAna::SetTunnelThreshold(Double_t value) +{ + if( value<=1 ){ mTunnelThreshold=value; } +} + + +TGraph* PeakAna::ConvertHistToGraph(TH1* hist, UInt_t numavgs) +{ + Int_t nbins = hist->GetNbinsX(); + //Only consider averaging when number of bins to average are greater than or equal to 1 or less than the number of bins + if( numavgs<1 || static_cast(numavgs)>nbins ){return 0;}//strictly greater than since GetNbinsX()+1 is overflow + TGraph* AvgGr = new TGraph(); + for( Int_t ibin=1; ibin<=nbins; ibin+=numavgs ){ + Double_t sum=0; + UInt_t counter=0;//In case number of averages is not same as how many sums were performed + for(Int_t i=ibin; i(numavgs) && i<=nbins; ++i ){ sum += hist->GetBinContent(i); ++counter; } + Double_t avg = sum/static_cast(counter); + //Set graph's point to center of x range + Double_t xlow = hist->GetBinLowEdge(ibin); + Double_t xhigh = hist->GetBinLowEdge(ibin+counter); + AvgGr->SetPoint(AvgGr->GetN(),(xlow+xhigh)/2.0,avg);//This is for knowing the x-value + } + return AvgGr; +} + +TGraph* PeakAna::ConvertPeaksToGraph(const std::vector &Peaks) +{ + if( Peaks.size()==0 ){return 0;} + TGraph* graph = new TGraph(); + for( UInt_t ipeak=0; ipeakSetPoint( ipeak,Peaks.at(ipeak).mPeakX,Peaks.at(ipeak).MidPoint() ); + } + return graph; +} + +void PeakAna::ConvertPeaksToGraph() +{ + TGraph* graph = new TGraph(); + for( UInt_t ipeak=0; ipeakSetPoint( ipeak,mPeaks.at(ipeak).mPeakX,mPeaks.at(ipeak).MidPoint() ); + } + SetData(graph); + ForceInternal();//new graph replaces old one so needs to be deleted + return; +} + +PeakAna PeakAna::ConvertPeaksToAna(const PeakAna &Ana) +{ + TGraph* graph = new TGraph(); + for( Int_t ipeak=0; ipeakSetPoint(ipeak, Ana.GetPeak(ipeak).mPeakX, Ana.GetPeak(ipeak).MidPoint() ); + } + PeakAna NewAna(Ana,graph); + NewAna.ForceInternal();//new graph so needs to be deleted + return NewAna; +} + +PeakAna* PeakAna::MeanFilter( Int_t sizeavgs, bool copy ) +{ + if( GetData()==0 ){return this;} + const int npoints = GetData()->GetN(); + if( sizeavgs==0 || npoints<=1 ){ return this; } + if( sizeavgs<0 ){ sizeavgs = -sizeavgs; } + double ynew[npoints]; + double* xdata = GetData()->GetX(); + double* ydata = GetData()->GetY(); + for( int ipoint=0; ipointAlso need to check minimum and maximum x-values?? + Double_t delxerr = 0.001*mDeltaX;//Since comparing doubles add 0.1% tolerance to mDeltaX + if( (nextpoint_m-1)>=0 ){ + if( mDeltaX<=0 || ((mDeltaX-delxerr)<=fabs(xdata[nextpoint_m-1]-lastx_m) && fabs(xdata[nextpoint_m-1]-lastx_m)<=(mDeltaX+delxerr)) ){ + //Next point is mDeltaX away so valid point + nextpoint_m -= 1; + lastx_m = xdata[nextpoint_m]; + if( mFilterWeights!=0 ){ sumvalues += mFilterWeights[sizeavgs-sizecounter-1]*ydata[nextpoint_m]; } + else{ sumvalues += ydata[nextpoint_m]; } + } + else{ + //Next point is not mDeltaX away so invalid point and add baseline and decrement lastx value + lastx_m -= mDeltaX; + if( mFilterWeights!=0 ){ sumvalues += mFilterWeights[sizeavgs-sizecounter-1]*Baseline(); } + else{ sumvalues += Baseline(); } + } + //Some y-value was added so increment number of points/weights + if( mFilterWeights!=0 ){ sumweights += mFilterWeights[sizeavgs-sizecounter-1]; } + else{ sumweights += 1; } + } + if( (nextpoint_p+1)ForceInternal(); + return ana; + } + if( mInternalSignal ){ + //@[July 6, 2022] > (Copy arrays in c++ )[https://stackoverflow.com/questions/16137953/is-there-a-function-to-copy-an-array-in-c-c] + std::copy( ynew, ynew+npoints, ydata); + //[July 3, 2022]>Taken from TGraph CtorAllocate(). Setting minimum and maximum to -1111 effectivley resets the minimum and maximum. + GetData()->SetMinimum(-1111); + GetData()->SetMaximum(-1111); + } + else{ + TGraph* graph = new TGraph(npoints,xdata,ynew ); + SetData(graph); + ForceInternal(); + } + ResetPeak(); + return this; +} + +PeakAna* PeakAna::GausFilter( Int_t sizeavgs, bool copy ) +{ + if( GetData()==0 ){return this;} + const int npoints = GetData()->GetN(); + if( sizeavgs==0 || npoints<=1 ){ return this; } + if( sizeavgs<0 ){ sizeavgs = -sizeavgs; } + double ynew[npoints]; + double* xdata = GetData()->GetX(); + double* ydata = GetData()->GetY(); + for( int ipoint=0; ipoint=0 ){ + if( mDeltaX<=0 || ((mDeltaX-delxerr)<=fabs(xdata[nextpoint_m-1]-lastx_m) && fabs(xdata[nextpoint_m-1]-lastx_m)<=(mDeltaX+delxerr)) ){ + //Next point is mDeltaX away so valid point + nextpoint_m -= 1; + lastx_m = xdata[nextpoint_m]; + if( mFilterWeights!=0 ){ sumvalues += mFilterWeights[sizeavgs-sizecounter-1]*ydata[nextpoint_m]; } + else{ sumvalues+=ydata[nextpoint_m]; } + } + else{ + //Next point is not mDeltaX away so invalid point and add baseline and decrement lastx value + //I should add baselines even if nextpoint+-1 is outside array (is this padding)? + lastx_m -= mDeltaX; + if( mFilterWeights!=0 ){ sumvalues += mFilterWeights[sizeavgs-sizecounter-1]*Baseline(); } + else{ sumvalues += Baseline(); } + } + } + else{ + //nextpoint_m is now negative so add padding by copying first point + if( mFilterWeights!=0 ){ sumvalues += mFilterWeights[sizeavgs-sizecounter-1]*ydata[0]; } + else{ sumvalues += ydata[0]; } + } + if( (nextpoint_p+1)=npoints so add padding by copying last point + if( mFilterWeights!=0 ){ sumvalues += mFilterWeights[sizeavgs+sizecounter+1]*ydata[npoints-1]; } + else{ sumvalues += ydata[npoints-1]; } + } + //Some y-value was added for either case so increment number of points/weights accordingly + if( mFilterWeights!=0 ){ + sumweights += mFilterWeights[sizeavgs-sizecounter-1]; + sumweights += mFilterWeights[sizeavgs+sizecounter+1]; + } + else{ sumweights += 2; } + } + ynew[ipoint] = sumvalues/sumweights; + } + + + if( copy ){ + TGraph* filtered = new TGraph(npoints,xdata,ynew); + PeakAna* ana = new PeakAna(*this,filtered); + ana->ForceInternal(); + return ana; + } + if( mInternalSignal ){ + //@[July 6, 2022] > (Copy arrays in c++ )[https://stackoverflow.com/questions/16137953/is-there-a-function-to-copy-an-array-in-c-c] + std::copy( ynew, ynew+npoints, ydata); + //[July 3, 2022]>Taken from TGraph CtorAllocate(). Setting minimum and maximum to -1111 effectivley resets the minimum and maximum. + GetData()->SetMinimum(-1111); + GetData()->SetMaximum(-1111); + } + else{ + TGraph* graph = new TGraph(npoints,xdata,ynew ); + SetData(graph); + ForceInternal(); + } + ResetPeak(); + return this; +} + +PeakAna PeakAna::ConvertPeaksToAna() +{ + return PeakAna::ConvertPeaksToAna(*this); +} + +bool PeakAna::PeakTunnel(const PeakWindow &window) const +{ + if( mTunnelThreshold>=0 ){ + Double_t prob = window.PeakTunnelProb(mG_Data,mTunnelScale,mTunnelSigma); + if( GetDebug()>1 ){ + LOG_DEBUG << "PeakAna::PeakTunnel - |prob:"<2 ){ LOG_DEBUG << "|scale:"<mTunnelThreshold ){return true; } + else{ return false; } + + } + return false; +} + +Double_t PeakAna::PeakProb(const PeakWindow &window, Double_t scale, Double_t sigma ) const +{ + return window.PeakTunnelProb(mG_Data,scale,sigma); +} + +void PeakAna::GetPossiblePeaks() +{ + std::vector PossiblePeak;//Gather all the possible occurances when signal is larger than Sigma + if( GetDebug()>1 ){LOG_DEBUG << "PeakAna - In GetPossiblePeaks" << endm; } + + Double_t baseline = Baseline();//Already checked baseline above + Double_t slopecutoff = BaselineSigma()*BaselineSigmaScale(); + if( GetDebug() > 2){ + LOG_DEBUG << "|baseline:"<1 ){LOG_DEBUG << "Finding Peak for cutoff " << slopecutoff << endm;} + //The idea is that you start with all things negative and loop through all the ADC values. When the slope is greater than the cutoff you save it and then the signal will rise and then fall so slope will eventually go negative and this is when you keep track of the end values. Once signal goes positive again or you no longer pass the slope cutoff stop and save the start and stop values. + if( GetDebug()>2 ){LOG_DEBUG << "PeakAna::GetPossiblePeaks:Start graph reading loop" << endm;} + Int_t npoints = mG_Data->GetN(); + for( Int_t ismp=0; ismpGetPoint(ismp,LX,LY); + mG_Data->GetPoint(ismp+1,RX,RY); + if( GetDebug()>1 ){LOG_DEBUG << " - |P:"<2 ){LOG_DEBUG << "'LX' too small"< mXRangeMax ){if( GetDebug()>2 ){LOG_DEBUG << "'LX' too large"< mYRangeMax){if( GetDebug()>2 ){LOG_DEBUG << "'LY' outside y range"< mYRangeMax){if( GetDebug()>2 ){LOG_DEBUG << "'RY' outside y range"<0 ){ + Double_t delxerr = 0.001*mDeltaX;//Since comparing doubles add 0.1% tolerance to mDeltaX + if( !((mDeltaX-delxerr)<=(RX-LX) && (RX-LX)<=(mDeltaX+delxerr)) ){ + if( GetDebug()>2 ){ LOG_DEBUG << " - |DeltaX:"<2 ){LOG_DEBUG << " - |Slope:"<0; however since I dynamically change the "end" time until the slope changes to positive this statement ensures that 'continue' only gets called on positive slope results and not negative slopes when it is trying to find the correct end time. + if( GetDebug()>2 ){LOG_DEBUG << " - peak|ismp:"<1 ){LOG_DEBUG << " + No StartTime" << endm;} + if( LY>baseline+slopecutoff && Slope>0 ){ + //Needs to be checked sequentially since we need to reject any points below the baseline+cutoff that may have a large slope + if( GetDebug()>1 ){LOG_DEBUG << " + Passed Slope and baselineCutoff setting as start time" << endm;} + peak.mStartX=LX;//Set start x-value + peak.mStartY=LY;//Set start y-value + LocalMax=LX;//Start checking local maximums + } + //If didn't pass thresholds for start time then continue to next point + } + else{//Found a suitable start time which indicates positve slope above some thershold + //The statements below will now check subsequent values and if the slope is >= 0 then we have wrong local max so set local max to next value. Note this means that localmax is not greater than End which is mXRangeMax+1. Once slope goes below 0 End will take that value and local max will be less than End. Once the slope changes again LocalMax will now be greater than End + if( GetDebug()>1 ){LOG_DEBUG << " + Found StartTime:"<< peak.mStartX << "|LocalMax:"<= 0 ){ + LocalMax=LX; //keep moving local max as long slope>=0 + if( Slope==0 && peak.mP_Peak>=0 ){ peak.mEndX=LX; } //if slope==0 and a peak was found i.e. mP_Peak>=0 then don't end finding but keep moving the end point i.e. peak.mEndX=LX. This condition is needed because sometimes a slope is zero and then will start decreasing again later. The importance of finding a peak is because it establishes the fact that the curve is convave and zero slope is ambiguous in which direction it will go. If slopes keep decreasing then expands peak window, if slope changes sign at next point algorithm will naturally stop. This is the desired behavior. Using the condition `Slope>0` does not produce the same the effect because sometimes this kind of zero slope behavior happens on the increasing slope side of the data and this should not be stopping the search because the next one could be positive. + if( GetDebug()>2){LOG_DEBUG<<" + Slope>=0|"<GetPoint(ismp,peak.mPeakX,peak.mPeakY); + } + if(GetDebug()>2){LOG_DEBUG<<" + Slope<0"<peak.mEndX || LY2 ){LOG_DEBUG << " + Slope positive again|StartX:"<0 ){ + Double_t leftprob = PossiblePeak.back().PeakTunnelProb(mG_Data,mTunnelScale,mTunnelSigma); + Double_t rightprob = peak.PeakTunnelProb(mG_Data,mTunnelScale,mTunnelSigma); + PeakWindow newpeak = PeakWindow::Combine(PossiblePeak.back(),peak, leftprob<=rightprob?true:false);//Take peak with highest probability + PossiblePeak[n-1] = newpeak; + } + } + else{ PossiblePeak.push_back(peak); }//Store windows + PossiblePeak.back().SetPeak(mG_Data); + if( GetDebug()>1){LOG_DEBUG << " * |Found:"; PossiblePeak.back().Print("debug"); LOG_DEBUG<2){LOG_DEBUG << " * |peakreset:"; peak.Print("debug"); LOG_DEBUG<0 ){ LOG_DEBUG << "PeakAna::GetPossiblePeaks:End graph reading loop|SizePeaks:" < &PossiblePeaks) +{ + if( PossiblePeaks.size()==0 ) + { + if( GetDebug()>0 ){LOG_DEBUG << "PeakAna::SearchForPeak - Error:Unable to find a valid peak\nReturning impossible index" << endm; } + return PossiblePeaks.size(); + } + else + { + if( GetDebug()>0 ) + { + LOG_DEBUG << "|Size of Possible peaks:"<1 ) + { + LOG_DEBUG << " - "; + LOG_DEBUG << "|Index:"<1 ) + { + LOG_DEBUG << " + "; + LOG_DEBUG << "|TrueIndex:"<=0 && mSearch.mStartX>0 && PossiblePeaks.at(peakindex).mP_Peak>0 ) + { + if(GetDebug()>1){PossiblePeaks.at(peakindex).Print("debug");LOG_DEBUG< &PossiblePeaks, const PeakWindow& search) +{ + mSearch = search; + return this->SearchForPeak(PossiblePeaks); +} + +Int_t PeakAna::SearchForPeak(const std::vector &PossiblePeaks, Double_t peak, Double_t width) +{ + SetSearchWindow(peak, width); + return this->SearchForPeak(PossiblePeaks); +} + +void PeakAna::MergeByProbability(std::vector& newpeaks) const +{ + if( !newpeaks.empty() ){ newpeaks.clear(); } + for( UInt_t i=0; i& newpeaks) const +{ + //In order to check peaks to merge globally make a vector that keeps the order of the peaks vector and whose value will be +1 to merge with right peak or -1 to merge with left peak, 0 is neutral. + //The idea is consecutive +1 would merge with consecutive -1. The zero is special since a +1 0 -1 pattern should be merged into a single peak but a 1 0 0 -1 pattern should not be + //std::cout << "PeakAna::MergeByChirality - Start" << std::endl; + std::vector mergelist; + for( UInt_t i=0; i mChiralityThreshold ){ mergelist.push_back(0); } + //else{ mergelist.push_back( this->MergeLeftOrRight(i) ); } + mergelist.push_back( this->MergeLeftOrRight(i) ); + } + + //std::cout << "PeakAna::MergeByChirality - peaksize:" << mPeaks.size() << std::endl; + //std::cout << "PeakAna::MergeByChirality - mergesize:" << mergelist.size() << std::endl; + std::vector< std::pair > coupledindexs = this->MergeIndices(mergelist); + + newpeaks.clear();//Clear any data in the peak + for( UInt_t i=0; i(i)==coupledindexs.at(j).first){ + PeakWindow combpeak=mPeaks.at( i ); + for( int k=coupledindexs.at(j).first+1; k<=coupledindexs.at(j).second; ++k ){ + Double_t thisprob = combpeak.PeakTunnelProb(mG_Data,mTunnelScale,mTunnelSigma); + Double_t otherprob = mPeaks.at(k).PeakTunnelProb(mG_Data,mTunnelScale,mTunnelSigma); + combpeak.Combine(mPeaks.at(k),thisprob<=otherprob?true:false); + }//k loop + newpeaks.push_back(combpeak); + i=coupledindexs.at(j).second; + combined = true; + } + } + if( !combined ){ newpeaks.push_back( mPeaks.at(i) ); } + } + //std::cout << "PeakAna::MergeByChirality - newpeaksize:" << newpeaks.size() << std::endl; +} + + +short PeakAna::MergeLeftOrRight(UInt_t index) const +{ + if( index>=mPeaks.size() ){return 0;}//invalid index + if( mPeaks.size()==1 ){return 0;}//Only 1 peak don't need to merge + //PeakWindow combined = mPeaks.at(index); + //if( mPeaks.at(index).PeakChiralityProb(mChiralityProbScale,mChiralityScale) > mChiralityThreshold ){return 0;} + Double_t leftchir = 0; + Double_t thischir = mPeaks.at(index).PeakChirality(mChiralityPeakScale,mChiralityScale); + Double_t rightchir = 0; + if( thischir==0 || thischir==1 ){return 0;}//Sanity check [Feb, 28 2022]>(If slope exactly zero then chirality is 1 for now, maybe use small fluctation around 1 since 0 should mean nonexistent exclusively at least when chirality only depnds on mse??) + + if( index==0 ){ + rightchir = mPeaks.at(index+1).PeakChirality(mChiralityPeakScale,mChiralityScale); + if( mPeaks.at(index+1).PeakChiralityProb(mChiralityProbScale,rightchir) > mChiralityThreshold ){ + rightchir = 0; } + } + else if( index+1==mPeaks.size() ){ + leftchir = mPeaks.at(index-1).PeakChirality(mChiralityPeakScale,mChiralityScale); + if( mPeaks.at(index-1).PeakChiralityProb(mChiralityProbScale,leftchir) > mChiralityThreshold ){ + leftchir = 0; } + } + else if( 0 mChiralityThreshold ){ + leftchir = 0; + } + if( mPeaks.at(index+1).PeakChiralityProb(mChiralityProbScale,rightchir) > mChiralityThreshold ){ + rightchir = 0; + } + } + else{ std::cout << "PeakAna::MergeLeftOrRight - WARNING:Invalid index" << std::endl; } + + //Check the 4 different cases + if( leftchir==0 && rightchir==0 ){ return 0; }//don't merge since left and right did not pass thresholds or do not exist + else if( leftchir!=0 && rightchir==0 ){//ignore right chirality + if( thischir*leftchir<0 ){return -1;}//opposite signs so merge with left + else{ return 0; } + } + else if( + leftchir==0 && rightchir!=0 ){//ignore left chirality + if( thischir*rightchir<0 ){return 1;}//opposite signs so merge with right + else{ return 0; } + } + else{//leftchir!=0 && rightchir!=0 + bool leftoppsign=false;//is left and this chirality opposite signs + if( (leftchir*rightchir)<0 ){leftoppsign=true;} + bool rightoppsign=false;//is right and this chirality opposite signs + if( (rightchir*thischir)<0 ){rightoppsign=true;} + + Double_t leftchirprob = mPeaks.at(index-1).PeakChiralityProb(mChiralityProbScale,leftchir); + Double_t rightchirprob = mPeaks.at(index+1).PeakChiralityProb(mChiralityProbScale,rightchir); + + if( leftoppsign ){ + if( rightoppsign ){//both left and right have opposite signs + //Merge with the lower chirality probability favoring left for equality so they will get merged?? + //or maybe need to keep going recursively by calling this function in the direction of the chirality, this way it will stop when chiralities are on longer going back and forth?? + if( leftchirprob0 ){return 1;} + else{ return -1; } + //if( leftchirprob > PeakAna::MergeIndices(std::vector& vec) const +{ + int firstpos_index = -1; + int lastneg_index = -1; + int zero_count = 0;//counting number of consecutive zeros when firstpos_index is not + std::vector< std::pair > mergeindexs; + for( int i=0; i(vec.size()); ++i ){ + //std::cout << "entrance|("<=0 ){ + if( firstpos_index<0 ){ if( vec.at(i)>0 ){firstpos_index = i;} } + else{ + if( lastneg_index>firstpos_index ){ + //store + std::pair mergepair(firstpos_index,lastneg_index); + mergeindexs.push_back(mergepair); + //std::cout << "store|firstpos:"<=0 ){ + if( vec.at(i)<0 ){ + lastneg_index = i; + } + } + } + return mergeindexs; +} + +double* PeakAna::GaussianMatrix2D(int rx, double sx, int ry, double sy, bool kNorm) +{ + if( rx<0 ){rx = -rx;} + if( ry<0 ){ry = -ry;} + if( sx==0 ){ sx = static_cast(rx)*0.5; } + if( sy==0 ){ sy = static_cast(ry)*0.5; } + const int nsizex = (2*rx+1); + const int nsizey = (2*ry+1); + + double GM_sum = 0; + //@[June 7, 2022]>[How to create 2D dynamic arrays](https://stackoverflow.com/questions/936687/how-do-i-declare-a-2d-array-in-c-using-new) + double* GM_2D = new double[nsizex*nsizey]; + for( int ix=0; ix 0){ mPainter->Paint(opt); } + else{ mPainter->Paint(mOption.Data()); } + } +} + +PeakAnaPainter* PeakAna::GetPainter(Option_t *opt) +{ + if( mPainter==0 ){ + mPainter = new PeakAnaPainter(); + mPainter->SetPeakAna(this); + } + mOption = opt; + return mPainter; +} + +void PeakAna::AddPeakStats(TPaveText* pave, const char* opt) +{ + if( pave==0 ){ return; } + TString option(opt); + bool statsalloption = false; + bool statsdetailoption = false; + + option.ToLower(); + if( option.Contains("a") ){ statsalloption = true; } + if( option.Contains("d") ){ statsdetailoption = true; } + + pave->AddText("--Found Peaks--"); + + if( mComputedIndex<0 ){ this->AnalyzeForPeak(); } + if( mComputedIndex == static_cast(mPeaks.size()) ){if( !statsalloption ){return;}}//if only writing stats for found peak don't process no found peak case + + //In future base the precision of the values based on the x/y range, and tunnel thershold value ?? + if( !statsalloption ){ + TText* t = pave->AddText( Form("I:%u|S:[%2.2f,%2.2f]|E[%2.2f,%2.2f]|P(%d,%2.2f,%2.2f)", + mComputedIndex, + (GetPeak(mComputedIndex)).mStartX, + (GetPeak(mComputedIndex)).mStartY, + (GetPeak(mComputedIndex)).mEndX, + (GetPeak(mComputedIndex)).mEndY, + (GetPeak(mComputedIndex)).mP_Peak, + (GetPeak(mComputedIndex)).mPeakX, + (GetPeak(mComputedIndex)).mPeakY + ) ); + t->SetTextAlign(11); + if( statsdetailoption ){ + TText* t2 = pave->AddText( Form(" + |Prob:%1.4f|Chir:%3.1f|ChirProb:%1.4f", + (GetPeak(mComputedIndex)).PeakTunnelProb(GetData(),TunnelScale(),TunnelSigma()), + (GetPeak(mComputedIndex)).PeakChirality(ChiralityPeakScale(),ChiralityScale()), + (GetPeak(mComputedIndex)).PeakChiralityProb(ChiralityProbScale(),ChiralityPeakScale(),ChiralityScale()) + ) ); + t2->SetTextAlign(11); + } + } + else{ + //pave->AddText( Form("|X:[%2.2f,%2.2f]|Y:[%2.2f,%2.2f]|BBS:(%d,%2.2f)|ProbSST:()") ); + for( Int_t ipeak=0; ipeakAddText( Form("|I:%u|S[%2.2f,%2.2f]|E[%2.2f,%2.2f]", + ipeak, + (GetPeak(ipeak)).mStartX, + (GetPeak(ipeak)).mStartY, + (GetPeak(ipeak)).mEndX, + (GetPeak(ipeak)).mEndY + ) ); + t->SetTextAlign(11); + if( mComputedIndex==ipeak ){ t->SetTextColor(kRed+1); } + if( statsdetailoption ){ + //|Chir:%3.1f|ChirProb:%1.4f + //(GetPeak(ipeak)).PeakChirality(ChiralityPeakScale(),ChiralityScale()), + //(GetPeak(ipeak)).PeakChiralityProb(ChiralityProbScale(),ChiralityPeakScale(),ChiralityScale()) + TText* t2 = pave->AddText( Form(" + |P(%d,%2.2f,%2.2f)|Prob:%1.4f", + (GetPeak(ipeak)).mP_Peak, + (GetPeak(ipeak)).mPeakX, + (GetPeak(ipeak)).mPeakY, + (GetPeak(ipeak)).PeakTunnelProb(GetData(),TunnelScale(),TunnelSigma()) + ) ); + t2->SetTextAlign(11); + if( mComputedIndex==ipeak ){ t2->SetTextColor(kRed+1); } + } + } + } + +} + +//Styling for graph +Color_t PeakAna::GetLineColor() const +{ if(mG_Data!=0 ){return mG_Data->GetLineColor();} return 0; } +Style_t PeakAna::GetLineStyle() const +{ if(mG_Data!=0 ){return mG_Data->GetLineStyle();} return 0; } +Width_t PeakAna::GetLineWidth() const +{ if(mG_Data!=0 ){return mG_Data->GetLineWidth();} return 0; } + +Color_t PeakAna::GetFillColor() const +{ if(mG_Data!=0 ){return mG_Data->GetFillColor();} return 0; } +Style_t PeakAna::GetFillStyle() const +{ if(mG_Data!=0 ){return mG_Data->GetFillStyle();} return 0; } + +Color_t PeakAna::GetMarkerColor() const +{ if(mG_Data!=0 ){return mG_Data->GetMarkerColor();} return 0; } +Size_t PeakAna::GetMarkerSize() const +{ if(mG_Data!=0 ){return mG_Data->GetMarkerSize();} return 0; } +Style_t PeakAna::GetMarkerStyle() const +{ if(mG_Data!=0 ){return mG_Data->GetMarkerStyle();} return 0; } + +void PeakAna::SetLineColor(Color_t color) +{ if(mG_Data!=0 ){mG_Data->SetLineColor(color);} } +void PeakAna::SetLineColorAlpha(Color_t color, Float_t alpha) +{ if(mG_Data!=0 ){mG_Data->SetLineColorAlpha(color,alpha);} } +void PeakAna::SetLineStyle(Style_t style) +{ if(mG_Data!=0 ){mG_Data->SetLineStyle(style);} } +void PeakAna::SetLineWidth(Width_t width) +{ if(mG_Data!=0 ){mG_Data->SetLineWidth(width);} } + +void PeakAna::SetFillColor(Color_t color) +{ if(mG_Data!=0 ){mG_Data->SetFillColor(color);} } +void PeakAna::SetFillColorAlpha(Color_t color, Float_t alpha) +{ if(mG_Data!=0 ){mG_Data->SetFillColorAlpha(color,alpha);} } +void PeakAna::SetFillStyle(Style_t style) +{ if(mG_Data!=0 ){mG_Data->SetFillStyle(style);} } + +void PeakAna::SetMarkerColor(Color_t color) +{ if(mG_Data!=0 ){mG_Data->SetMarkerColor(color);} } +void PeakAna::SetMarkerColorAlpha(Color_t color, Float_t alpha) +{ if(mG_Data!=0 ){mG_Data->SetMarkerColorAlpha(color,alpha);} } +void PeakAna::SetMarkerSize(Size_t size) +{ if(mG_Data!=0 ){mG_Data->SetMarkerSize(size);} } +void PeakAna::SetMarkerStyle(Style_t style) +{ if(mG_Data!=0 ){mG_Data->SetMarkerStyle(style);} } + +//Styling for peak painter +Color_t PeakAna::GetBaseLineColor() const +{ if( mPainter!=0 ){ return mPainter->GetBaseLineColor(); } return 0;} +Style_t PeakAna::GetBaseLineStyle() const +{ if( mPainter!=0 ){ return mPainter->GetBaseLineStyle(); } return 0;} +Width_t PeakAna::GetBaseLineWidth() const +{ if( mPainter!=0 ){ return mPainter->GetBaseLineWidth(); } return 0;} + +Color_t PeakAna::GetHitLineColor() const +{ if( mPainter!=0 ){ return mPainter->GetHitLineColor(); } return 0;} +Style_t PeakAna::GetHitLineStyle() const +{ if( mPainter!=0 ){ return mPainter->GetHitLineStyle(); } return 0;} +Width_t PeakAna::GetHitLineWidth() const +{ if( mPainter!=0 ){ return mPainter->GetHitLineWidth(); } return 0;} + +void PeakAna::SetBaseLineColor(Color_t color) +{GetPainter()->SetBaseLineColor(color);} +void PeakAna::SetBaseLineColorAlpha(Color_t color,Float_t alpha) +{GetPainter()->SetBaseLineColorAlpha(color,alpha);} +void PeakAna::SetBaseLineStyle(Style_t style) +{GetPainter()->SetBaseLineStyle(style);} +void PeakAna::SetBaseLineWidth(Width_t width) +{GetPainter()->SetBaseLineWidth(width);} + +void PeakAna::SetHitLineColor(Color_t color) +{GetPainter()->SetHitLineColor(color);} +void PeakAna::SetHitLineColorAlpha(Color_t color,Float_t alpha) +{GetPainter()->SetHitLineColorAlpha(color,alpha);} +void PeakAna::SetHitLineStyle(Style_t style) +{GetPainter()->SetHitLineStyle(style);} +void PeakAna::SetHitLineWidth(Width_t width) +{GetPainter()->SetHitLineWidth(width);} + +Color_t PeakAna::GetPeakLineColorS(UInt_t peakidx) const +{return GetPeak(peakidx).GetStartLineColor();} +Color_t PeakAna::GetPeakLineColorE(UInt_t peakidx) const +{return GetPeak(peakidx).GetEndLineColor();} +Style_t PeakAna::GetPeakStyleS(UInt_t peakidx) const +{return GetPeak(peakidx).GetStartLineStyle();} +Style_t PeakAna::GetPeakStyleE(UInt_t peakidx) const +{return GetPeak(peakidx).GetEndLineStyle();} +Width_t PeakAna::GetPeakWidthS(UInt_t peakidx) const +{return GetPeak(peakidx).GetStartLineWidth();} +Width_t PeakAna::GetPeakWidthE(UInt_t peakidx) const +{return GetPeak(peakidx).GetEndLineWidth();} + +Color_t PeakAna::GetFoundPeakLineColorS() const +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ return GetPeak(mComputedIndex).GetStartLineColor(); } return 0; } +Color_t PeakAna::GetFoundPeakLineColorE() const +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ return GetPeak(mComputedIndex).GetEndLineColor(); } return 0; } +Style_t PeakAna::GetFoundPeakStyleS() const +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ return GetPeak(mComputedIndex).GetStartLineStyle(); } return 0; } +Style_t PeakAna::GetFoundPeakStyleE() const +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ return GetPeak(mComputedIndex).GetEndLineStyle(); } return 0; } +Width_t PeakAna::GetFoundPeakWidthS() const +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ return GetPeak(mComputedIndex).GetStartLineWidth(); } return 0; } +Width_t PeakAna::GetFoundPeakWidthE() const +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ return GetPeak(mComputedIndex).GetEndLineWidth(); } return 0; } + +void PeakAna::SetFoundPeakLineColorS(Color_t s_color) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetStartLineColor(s_color); } } +void PeakAna::SetFoundPeakLineColorE(Color_t e_color) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetEndLineColor(e_color); } } +void PeakAna::SetFoundPeakLineColor(Color_t s_color, Color_t e_color) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetStartLineColor(s_color); GetPeak(mComputedIndex).SetEndLineColor(e_color); } } +void PeakAna::SetFoundPeakStyleS(Style_t s_style) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetStartLineStyle(s_style); } } +void PeakAna::SetFoundPeakStyleE(Style_t e_style) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetEndLineStyle(e_style); } } +void PeakAna::SetFoundPeakStyle(Style_t s_style,Style_t e_style) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetStartLineStyle(s_style); GetPeak(mComputedIndex).SetEndLineStyle(e_style); } } +void PeakAna::SetFoundPeakWidthS(Width_t s_width) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetStartLineWidth(s_width); } } +void PeakAna::SetFoundPeakWidthE(Width_t e_width) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetEndLineWidth(e_width); } } +void PeakAna::SetFoundPeakWidth(Width_t s_width, Width_t e_width) +{ if( 0<=mComputedIndex && mComputedIndex(mPeaks.size()) ){ GetPeak(mComputedIndex).SetStartLineWidth(s_width); GetPeak(mComputedIndex).SetEndLineWidth(e_width); } } + +void PeakAna::SetPeakLineColorS(UInt_t peakidx, Color_t s_color) +{ GetPeak(peakidx).SetStartLineColor(s_color); } +void PeakAna::SetPeakLineColorE(UInt_t peakidx, Color_t e_color) +{ GetPeak(peakidx).SetStartLineColor(e_color); } +void PeakAna::SetPeakLineColor(UInt_t peakidx, Color_t s_color, Color_t e_color) +{ GetPeak(peakidx).SetStartLineColor(s_color); GetPeak(peakidx).SetStartLineColor(e_color); } +void PeakAna::SetPeakLineColorAlphaS(UInt_t peakidx, Color_t s_color,Float_t s_alpha) +{ GetPeak(peakidx).SetStartLineColorAlpha( s_color, s_alpha); } +void PeakAna::SetPeakLineColorAlphaE(UInt_t peakidx, Color_t e_color,Float_t e_alpha) +{ GetPeak(peakidx).SetEndLineColorAlpha( e_color, e_alpha); } +void PeakAna::SetPeakLineColorAlpha(UInt_t peakidx, Color_t s_color,Float_t s_alpha, Color_t e_color, Float_t e_alpha) +{ GetPeak(peakidx).SetStartLineColorAlpha( s_color, s_alpha); GetPeak(peakidx).SetEndLineColorAlpha( e_color, e_alpha); } +void PeakAna::SetPeakLineStyleS(UInt_t peakidx, Style_t s_style) +{ GetPeak(peakidx).SetStartLineStyle(s_style); } +void PeakAna::SetPeakLineStyleE(UInt_t peakidx, Style_t e_style) +{ GetPeak(peakidx).SetEndLineStyle(e_style); } +void PeakAna::SetPeakLineStyle(UInt_t peakidx, Style_t s_style, Style_t e_style) +{ GetPeak(peakidx).SetStartLineStyle(s_style); GetPeak(peakidx).SetEndLineStyle(e_style); } +void PeakAna::SetPeakLineWidthS(UInt_t peakidx, Width_t s_width) +{ GetPeak(peakidx).SetStartLineWidth(s_width); } +void PeakAna::SetPeakLineWidthE(UInt_t peakidx, Width_t e_width) +{ GetPeak(peakidx).SetEndLineWidth(e_width); } +void PeakAna::SetPeakLineWidth(UInt_t peakidx, Width_t s_width, Width_t e_width) +{ GetPeak(peakidx).SetStartLineWidth(s_width); GetPeak(peakidx).SetEndLineWidth(e_width); } + +void PeakAna::SetAllPeakLineColor(Color_t s_color, Color_t e_color) +{ + for( UInt_t ipeak=0; ipeak 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?? +*/ +/*! + The purpose of this class is to find peaks in a sample data and the range of those peaks; i.e. where the first derivative is zero and second derivate changes sign from positive to negative and back to positive as this indicates a local max. + +There is a helper data class #PeakWindow which holds the found peak data. +*/ + + +#ifndef PEAKANA_H +#define PEAKANA_H + +//C++ Headers +#include +#include +#include //std::move (Needs C++ 11) and std::pair +#include + +//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" +//#include "PeakAnaVirtualPainter.h" +class PeakAnaPainter; + +class PeakAna : public TObject +{ +public: + PeakAna();//!< Default constructor that always creates a new TGraph + PeakAna( int size, double* xvals, double* yvals ); //!< Constructor using array of x and y values with known size + PeakAna( TGraph* Sig ); //!< Constructor with a known TGraph + PeakAna( TH1* hist); //!< Constructor with a histogram, histogram gets converted to graph object + PeakAna(const PeakAna &OldAna,TGraph* graph=0); //!< Copy constructor can be called with a new graph + PeakAna& operator=(const PeakAna& rhs); //!< Assignment operator doesn't clone graph + virtual ~PeakAna();//!< Destructor + + virtual void Copy(TObject& obj) const; //!< Internal method use Clone instead. + virtual TObject* Clone(const char* newname) const; //!< Clones internal graph as opposed to just copying the pointer + + /**@brief Add peak information to a "statistics" box + + "statistics" box can be a TPaveText or a TPaveStats since the latter inherits from the former. + @param pave TPaveText to write peak information + @param opt options for what information to add to stats box\n + "" (No option) means just show found peak + "a" means show all found peaks + "d" means show all detail information for the peak + */ + virtual void AddPeakStats(TPaveText* pave, const char* opt=""); + + /** @brief Converts a histogram to a graph + + Treats each bin in the histogram as a point in a TGraph and the bin content as the y-value. Also, can be used to do a running average over "numavgs" bins + + @param graph histogram to be converted to a graph + @param numavgs number of bins to average over when creating the graph (running average), 1 means no running average + @return TGraph with bin data as points + */ + static TGraph* ConvertHistToGraph(TH1* graph, UInt_t numavgs=1); + + /**@brief Convert each peak position into a TGraph + + The purpose of this function is to help with really noisy data where the "peaks" are really just noise. Essentially the idea is that each "peak" should actually be a point on a kind of "running average" in this way you capture the overall shape of the data without noise. + @param Peaks vector of PeakWindows to convert to a TGraph + @return TGraph with each peak in "Peaks" as a point + */ + static TGraph* ConvertPeaksToGraph(const std::vector &Peaks); + void ConvertPeaksToGraph(); //!< same as #ConvertPeaksToGraph() but creates a new graph that replaces old one + static PeakAna ConvertPeaksToAna(const PeakAna &Ana); //!< same as #ConvertPeaksToGraph() but returns a new PeakAna with same settings without modifying old one + PeakAna ConvertPeaksToAna();//!< same as #ConvertPeaksToAna() but returns new PeakAna with same settings without modifying current one + + /**@brief Finds all possible Peaks in signal with some criteria + + Uses a second derivative test to find all peaks in the data. Found peaks must be above the baseline, and the slope must change by #mBaseline+#mBaselineSigma*#mBaselineSigmaScale. These values can be set with the appropriate function calls. Also, checks each peak's probability and merges peaks below some threshold if set.\n + Called by #AnalyzeForPeak(); + */ + virtual void GetPossiblePeaks(); + + /**@brief searches a vector of #PeakWindow's for a specific peak based on input search criteria + + Search criteria is stored in #mSearch where the #mStartX is the peak location to search for and #mEndX is the tolerance (width) around which the peak should be found.\n + Called by #AnalyzeForPeak() on #mPeaks + + @param PossiblePeaks vector of peaks to search + @return index in #mPeaks of found peak that matched search criteria + */ + virtual Int_t SearchForPeak(const std::vector &PossiblePeaks ); + Int_t SearchForPeak(const std::vector &PossiblePeaks, const PeakWindow& search);//!< same as #SearchForPeak() but and set search criteria using #PeakWindow + Int_t SearchForPeak(const std::vector &PossiblePeaks, Double_t peak, Double_t width);//!< same as #SearchForPeak() and set search window by peak location and width + + //Class methods can be used to store data and preserve found peaks for fitting + /**@brief Main analysis method for finding peaks + + First filters the data if the #mFilter is set, then calls #GetPossiblePeaks(), then calls #SearchForPeak(). + @return index of found peak, if no found peak then returns size of #mPeaks + */ + virtual Int_t AnalyzeForPeak(); + virtual Int_t AnalyzeForPeak(Double_t peak, Double_t width);//& newpeaks) const; + virtual void MergeByChirality(std::vector& newpeaks) const; + virtual short MergeLeftOrRight(UInt_t index) const; + + Double_t Baseline()const{return mBaseline;} //!< @return #mBaseline + Double_t BaselineSigma()const{return fabs(mBaselineSigma); } //!< @return fabs(#mBaselineSigma) + Double_t BaselineSigmaScale()const{return mBaselineSigmaScale;} //!< @return #mBaselineSigmaScale + Double_t MinX()const{ return mXRangeMin; } //!< @return #mXRangeMin + Double_t MinY()const{ return mYRangeMin; } //!< @return #mYRangeMin + Double_t MaxX()const{ return mXRangeMax; } //!< @return #mXRangeMax + Double_t MaxY()const{ return mYRangeMax; } //!< @return #mYRangeMax + Double_t SearchPeak()const{ return mSearch.mStartX; } //!< @return expected peak location for searching + Double_t SearchWidth()const{ return mSearch.mEndX; } //!< @return width from expected peak location to search + Double_t TunnelScale()const{ return mTunnelScale; } //!< @return #mTunnelScale + Double_t TunnelSigma()const{ return mTunnelSigma; } //!< @return #mTunnelSigma + Double_t TunnelThreshold()const{ return mTunnelThreshold; } //!< @return #mTunnelThreshold + + /**@brief Draw method for #PeakAna + + Options are semicolon seperated list with graph options first then the peak qa drawing options and finally the stats box drawing option.\n + peak qa drawing options(case insensitive):\n + - "R" for range\n + - "B" for baselines\n + - "F" found peak qa\n + - "P" for full peak qa\n + - "A" for all, which combines "P" and "B"\n + stats drawing options(case insensitive):\n + - "S" is to show just the found peak\n + - "A" is to show for all peaks\n + - "D" is whether to show detailed printout or not (works with option "S" or "A") + + Example1: "APL;FB;S" means draw graph using option "APL", draw only found peak qa and baselines, and create a stats box showing only basic information about the found peak. + Example2: ";A;AD" means don't draw graph but draw all peaks and put detailed information of all peaks in a stats box + Example3: "PL;A" means draw graph with option "PL" and draw all peaks but don't draw a stats box + */ + virtual void Draw(Option_t *opt=""); + virtual void Paint(Option_t *opt=""); //AnalyzeForPeak();} return mFoundPeak.mStartX;}//!< Found Signal starting x-value + Double_t PeakEnd(){ if( mComputedIndex<0 ){this->AnalyzeForPeak();} return mFoundPeak.mEndX;} //!< Found Signal ending x-value + Double_t PeakX(); //!< x-value of found signal peak + Double_t PeakY(); //!< y-value of found signal peak + void PeakXY(Double_t &xval, Double_t &yval); //!< get peak x, and y value directly by reference + /**@brief test whether a given peak satisifies peak tunnel parameters + + Evaluates a #PeakWindow's probability with function parameters stored in this class, and the internal graph object. + @param window #PeakWindow to be evaluated + @return true if peak probability>#mTunnelThreshold, false otherwise + */ + bool PeakTunnel(const PeakWindow &window) const; + + /**@brief compute a given #PeakWindow probability using internal graph + + Evaluates a given #PeakWindow's prbability using internal graph and external parameters for the probability function + @param window #PeakWindow to compute probability + @param scale x-value scale for probability, see #mTunnelScale + @param sigma y-value scale for probability, see #mTunnelSigma + @return computed probability + */ + Double_t PeakProb(const PeakWindow& window, Double_t scale, Double_t sigma ) const; + + virtual void Print(Option_t* opt="") const; //!< Print peak information + void ResetPeak(); //!< Resets values associated with peak finding + + void SetDebug(UInt_t level){ mDEBUG=level; } //!< @param level sets debug level #mDEBUG + /**@brief sets new data for #PeakAna + + This will change #mG_Data to the graph object passed in and calls #ResetPeak(). Since graph is passed as external parameter this means this object should not delete it in the future so #mInternalSignal is set to false. If you want to change this behavior so this class deletes #mG_Data call #ForceInternal() + @param graph sets #mG_Data + */ + virtual void SetData(TGraph* graph); + + /**@brief sets new data for #PeakAna using histogram object + + Works like #SetData() but first calls #ConvertHistToGraph() to convert a histogram to a graph. Since it creates a new graph #mInternalSignal is set to true so that the new graph gets deleted properly. The second argument can be used to do a running average over 'numavgs' number of bins. + @param hist histogram to use as graph data + @param numavgs controls how many bins you want to average over when creating graph (1 means no average) + */ + virtual void SetData(TH1* hist, UInt_t numavgs=1); + + void SetBaseline(Double_t base, Double_t sigma); //!< @param base sets #mBaseline @param sigma sets #mBaselineSigma + void SetBaselineSigmaScale(Double_t scale); //!< @param scale sets #mBaselineSigmaScale + void SetContinuity(Double_t val){ mDeltaX = val; } //!< @param val sets #mDeltaX + + /**@brief Set the filter to use when peak finding + + There are two filters that can be used: a Mean filter which just averages over small samples of the data at a time, and a Gaussian filter which does a weighted average over the data samples, where the weights resemble a Gaussian function see #GaussianMatrix2D(). The data sample averaging size is determined by #mFilterScale. Note that a filter does not alter the size of the data. + @param filter sets #mFilter (0=none, 1=Mean, 2=Gaussian) + @param scale sets #mFilterScale (number of points to group together when filtering) + @param sigma sigma of Gaussian to use to compute #mFilterWeights (0 means sigma=mFilterScale/2) + */ + void SetFilter( UInt_t filter, Int_t scale, Double_t sigma=0 ); + + /**@brief Sets the absolute range of the data + + This sets the varaibles that tell the finder the range of x,y values possible for the data. This kind of range is needed for the algorithm to work properly since it only checks for peaks inside that range. + @param xmin sets #mXRangeMin + @param ymin sets #mYRangeMin + @param xmax sets #mXRangeMax + @param ymax sets #mYRangeMax + */ + void SetRange( Double_t xmin, Double_t ymin, Double_t xmax, Double_t ymax); + + /**@brief Sets peak search parameters + + This function is needed if you want to pick/check out a particular peak in the vector of #PeakWindows. It allows you to set the x-value of where you expect or want a peak and how far from that ideal position (width) you would like to check for a peak. + @param peak x-location where you expect or want a peak + @param width how far in x you expect or want a peak from the ideal position of "peak" + */ + void SetSearchWindow(Double_t peak, Double_t width); + + void SetTunnelScale(Double_t value); //!< @param value sets #mTunnelScale + void SetTunnelSigma(Double_t value); //!< @param value sets #mTunnelSigma + void SetTunnelPars(Double_t scale, Double_t sigma); //!< @param scale sets #mTunnelScale @param sigma sets #mTunnelSigma + void SetTunnelThreshold(Double_t value); //!< @param value sets #mTunnelThreshold + + 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 );//!< Overwrite peak location in #mFoundpeak. WARNING doesn't set start and end values nor overwrite #mComputedIndex + void SetWindow(const Int_t start, const Int_t end );//!< Overwrite peak start and end values in #mFoundPeak. WARNING doesn't set peak values nor overwrite #mComputedIndex + void SetWindow(PeakWindow window); //!< Overwrite #mFoundPeak with a different #PeakWindow. WARNING doesn't overwrite #mComputedIndex + + //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();//!< Initialize internal varaibles + + /**@brief Index in #mPeaks where a peak was found + + This variables is negative if #AnalyzeForPeak() was not called.\n + It is equal to size of #mPeaks if #AnalyzeForPeak() was called but no peak was found.\n + All other times it equals the index in #mPeaks where the peak was found. The found peak is also stored in #mFoundPeak.\n + This makes it easy to tell whether #AnalyzeForPeak() was called and whether or not a peak was found. + */ + Int_t mComputedIndex; + + std::vector mPeaks; //!< vector that stores all the found peaks in the data + + /**@brief Copy of found peak in #mPeaks + + This gets set after #SearchForPeak() is called and finds a peak with given search parameters. It is separate for ease of access and so it can be overwritten for testing and debugging purposes + */ + PeakWindow mFoundPeak; + + /**@brief Variable that defines peak search parameters + + The values are stored in a #PeakWindow object but only #PeakWindow::mStartX is used for the peak location and #PeakWindow::mEndX is used for the width. The width represents the +- window from #PeakWindow::mStartX to look for a peak. This is used by #SearchForPeak(). Default is peak at 0 with a width of 1. + */ + PeakWindow mSearch; + + /**@brief Minimum threshold to search for a peak + + Threshold for peak start and end values meaning y-value must exceed this variable to even check if it contains a peak. y-values below this variable will be ignored. This defines the noise level + */ + Double_t mBaseline; + + /**@brief Error on #PeakAna::mBaseline + + This is used to characterize how much noise there is. It will be used with #PeakAna::mBaseline to determine the threshold y-value for the data. It effectively sets the resolution of the finder. + This is the standard deviation (sigma) around the baseline for ADCs (Also used as a threshold for the hit) error in baseline (fabs in case sigma is negative), //The sigma will determine the threshold value above baseline('value') i.e. the resolution of the finder + */ + Double_t mBaselineSigma; + + /**@brief scale on #PeakAna::mBaselineSigma to determine final threshold + + This is a multiplicative factor on #PeakAna::mBaselineSigma that is used to determine the final thresholds on the y-value and slope to start the second derivative test in #GetPossiblePeaks(). This "hit" threshold is #PeakAna::mBaseline + #PeakAna::mBaselineSigma * #PeakAna::mBaselineSigmaScale. The default is 4. This threshold also applies to the first positive slope to prevent finding data with small changes as potential peaks. + */ + Double_t mBaselineSigmaScale; + + Double_t mMaxX;//!< x-value where global maximum occurs + Double_t mMaxY;//!< y-value of global maximum + + /**@brief graph known delta-x + + Set this value if you know how far apart each x point in the data/graph needs to be. This way the finder can handle discontinuities (i.e. each adjacent point is not a fixed x-distance) in the actual data. Default is -1, which means don't check continuity. + */ + Double_t mDeltaX; + + /**@brief TGraph that stores the x,y data + + This graph object will not be deleted by this class unless #PeakAna::mInternalSiganl=true + A TGraph is used because it is easier to process but may want to change to dynamic arrays for speed?? + */ + TGraph* mG_Data; + + //A bit redundant but used for internal checks + Double_t mXRangeMin; //!< Absolute possible x-value minimum of data + Double_t mXRangeMax; //!< Absolute possible x-value maximum of data + Double_t mYRangeMin; //!< Absolute possible y-value minimum of data + Double_t mYRangeMax; //!< Absolute possible y-value maximum of data + + Double_t mTunnelScale; //!< Scale on exponential for determining tunneling probability (default is 1) see #PeakWindow::PeakTunelProb() + Double_t mTunnelSigma; //!< Sigma for Gaussian for determining tunneling probability (default is 1) see #PeakWindow::PeakTunelProb() + Double_t mTunnelThreshold; //!< Cutoff probability for peak tunneling method. If threshold less than 0 (default) then skip peak tunnel method, see #PeakTunnel() + + 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) + + /**@brief Filter to use in #AnalyzeForPeak() + + see #SetFilter(). Note that filters do not change the size of the data\n + - 0 = No Filter (default)\n + - 1 = Mean Filter\n + - 2 = Gauss Filter/Blur\n + */ + UInt_t mFilter; + Double_t* mFilterWeights; //!< Array of weights to use in filtering, size is 2*mFilterScale+1 + Int_t mFilterScale; //!< How many points to group together when applying filter + + PeakAnaPainter* mPainter; //!< Painter for this class + + /**@brief Generates up to a 2D matrix of Gaussian weights + + Generates Gaussian weights to use in Gaussian filtering in a flattened 2D matrix. + @param rx number of columns to generate weights + @param sx sigma of Gaussian to use in column weights + @param ry number of rows to generate weights + @param sy sigma of Gaussian to use in row weights + @param kNorm if true use a normalized Gaussian i.e. sum of all weights=1 + @return 1D array containing weights + */ + 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 > MergeIndices(std::vector& vec) const; + + TString mOption; //!< Drawing option + +private: + UInt_t mDEBUG; //!< debug level: 0 is none, 1 is info, 2 is verbose, 3 is deep + bool mInternalSignal; //!< If true will call delete on #mG_Data + + ClassDef(PeakAna,3); +}; + +#endif diff --git a/StRoot/StFcsWaveformFitMaker/PeakAnaPainter.cxx b/StRoot/StFcsWaveformFitMaker/PeakAnaPainter.cxx new file mode 100644 index 00000000000..91922437e12 --- /dev/null +++ b/StRoot/StFcsWaveformFitMaker/PeakAnaPainter.cxx @@ -0,0 +1,243 @@ +#include "PeakAnaPainter.h" +#include "PeakAna.h" + +ClassImp(PeakAnaPainter) + +PeakAnaPainter::PeakAnaPainter() +{ + Init(); +} + +void PeakAnaPainter::Init() +{ + mPA = 0; + + mTheBaseLine = 0; + mTheHitLine = 0; + + mPaveT_PA = 0; + mGraphOption = ""; + mPeakOption = ""; + mStatsOption = ""; +} + +PeakAnaPainter::~PeakAnaPainter() +{ + delete mTheBaseLine; + delete mTheHitLine; mTheHitLine=0; + delete mPaveT_PA; mPaveT_PA=0; +} + +void PeakAnaPainter::SetPeakAna(PeakAna* ana) +{ + if( ana==0 ){return;} + mPA = ana; +} + +void PeakAnaPainter::CleanPainter() +{ + if( mPaveT_PA!=0 ){ mPaveT_PA->Clear(); } +} + +Color_t PeakAnaPainter::GetBaseLineColor() +{if(mTheBaseLine==0){return 0;}else{return mTheBaseLine->GetLineColor();} } +Style_t PeakAnaPainter::GetBaseLineStyle() +{if(mTheBaseLine==0){return 0;}else{return mTheBaseLine->GetLineStyle();} } +Width_t PeakAnaPainter::GetBaseLineWidth() +{if(mTheBaseLine==0){return 0;}else{return mTheBaseLine->GetLineWidth();} } + +Color_t PeakAnaPainter::GetHitLineColor() +{if(mTheHitLine==0){return 0;}else{return mTheHitLine->GetLineColor();} } +Style_t PeakAnaPainter::GetHitLineStyle() +{if(mTheHitLine==0){return 0;}else{return mTheHitLine->GetLineStyle();} } +Width_t PeakAnaPainter::GetHitLineWidth() +{if(mTheHitLine==0){return 0;}else{return mTheHitLine->GetLineWidth();} } + +void PeakAnaPainter::SetBaseLineColor(Color_t color) +{if(mTheBaseLine==0){mTheBaseLine=new TLine();} mTheBaseLine->SetLineColor(color); } +void PeakAnaPainter::SetBaseLineColorAlpha(Color_t color,Float_t alpha) +{if(mTheBaseLine==0){mTheBaseLine=new TLine();} mTheBaseLine->SetLineColorAlpha(color,alpha); } +void PeakAnaPainter::SetBaseLineStyle(Style_t style) +{if(mTheBaseLine==0){mTheBaseLine=new TLine();} mTheBaseLine->SetLineStyle(style); } +void PeakAnaPainter::SetBaseLineWidth(Width_t width) +{if(mTheBaseLine==0){mTheBaseLine=new TLine();} mTheBaseLine->SetLineWidth(width); } + +void PeakAnaPainter::SetHitLineColor(Color_t color) +{if(mTheHitLine==0){mTheHitLine=new TLine();} mTheHitLine->SetLineColor(color); } +void PeakAnaPainter::SetHitLineColorAlpha(Color_t color,Float_t alpha) +{if(mTheHitLine==0){mTheHitLine=new TLine();} mTheHitLine->SetLineColorAlpha(color,alpha); } +void PeakAnaPainter::SetHitLineStyle(Style_t style) +{if(mTheHitLine==0){mTheHitLine=new TLine();} mTheHitLine->SetLineStyle(style); } +void PeakAnaPainter::SetHitLineWidth(Width_t width) +{if(mTheHitLine==0){mTheHitLine=new TLine();} mTheHitLine->SetLineWidth(width); } + +void PeakAnaPainter::Paint(Option_t* opt) +{ + if( mPA==0 ){ return; } + std::string option(opt); + + CleanPainter(); + bool drawgraph = true; + bool region = false; + bool drawbaselines = false; + bool drawfoundpeakqa = false; + bool drawfullpeakqa = false; + + std::size_t firstcolon = option.find(";"); + if( firstcolon!=std::string::npos ){ + mGraphOption = option.substr(0,firstcolon); + std::size_t secondcolon = option.find(";",firstcolon+1); + if( secondcolon!=std::string::npos ){ + mPeakOption = option.substr(firstcolon+1,secondcolon-firstcolon-1); + mStatsOption = option.substr(secondcolon+1); + std::size_t thirdcolon = option.find(";",secondcolon+1); + if( thirdcolon!=std::string::npos ){ LOG_WARN << "PeakAnaPainter::Paint(): Too many semicolons in draw option" << endm; } + } + else{ mPeakOption = option.substr(firstcolon+1); } + } + else{ mGraphOption = option; } + + mPeakOption.ToLower(); + if( mGraphOption.Length()==0 ){ drawgraph=false; } + //if( peakopts.Contains("e") ){ drawgraph=false; } + if( mPeakOption.Contains("r") ){ region=true; } + if( mPeakOption.Contains("b") ){ drawbaselines=true; } + if( mPeakOption.Contains("f") ){ drawfoundpeakqa=true; } + if( mPeakOption.Contains("p") ){ drawfullpeakqa=true; } + if( mPeakOption.Contains("a") ){ drawbaselines=true; drawfullpeakqa=true; } + + if( drawfoundpeakqa && drawfullpeakqa ){ drawfoundpeakqa=false; }//Redundant to do both found and all peaks + if( drawgraph ) { this->PaintRawData(); } + if( region && drawgraph ) { this->PaintFoundPeak(); } //Doesn't make sense to restrict graph range if graph is not being drawn + if( drawbaselines ) { this->PaintBaselines(); } + if( drawfoundpeakqa ) { this->PaintFoundPeakQa(); } + if( drawfullpeakqa ) { this->PaintPeakRanges(); } + if( mStatsOption.Length()!=0 ) { this->PaintStats(); } +} + + +TPaveText* PeakAnaPainter::MakePaveText(Double_t xmin, Double_t ymin, Double_t xmax, Double_t ymax) +{ + if( mPaveT_PA==0 ){ mPaveT_PA = new TPaveText(xmin,ymin,xmax,ymax,"NB NDC"); } + else{ + mPaveT_PA->Clear(); + } + mPaveT_PA->SetFillColorAlpha(kWhite,1);//Fully opaque by default + //mPaveT_PA->SetFillStyle(4000); + return mPaveT_PA; +} + +bool PeakAnaPainter::ValidGraph() +{ + if( mPA==0 || mPA->GetData()==0 || mPA->GetData()->GetN()==0 ){return false;} + else{ return true; } +} + +void PeakAnaPainter::PaintRawData( ) +{ + if( !ValidGraph() ){ return;} + mPA->GetData()->Paint(mGraphOption.Data()); + return; +} + +void PeakAnaPainter::PaintFoundPeak( ) +{ + if( !ValidGraph()){return;} + + Double_t Xmin = mPA->PeakStart(); + Double_t Xmax = mPA->PeakEnd(); + if( !(Xmin>mPA->MaxX() || Xmax>mPA->MaxX()) ){mPA->GetData()->GetXaxis()->SetRangeUser(Xmin-5,Xmax+5);} +} + +void PeakAnaPainter::PaintFoundPeakQa() +{ + if( !ValidGraph() ){return;} + + Int_t computedindex = mPA->FoundPeakIndex(); + if( computedindex<0 ){ computedindex = mPA->AnalyzeForPeak(); } + this->PaintBaselines(); + this->PaintFoundRange(); + + return; +} + +void PeakAnaPainter::PaintPeakQa( ) +{ + if( !ValidGraph() ){return;} + + Int_t computedindex = mPA->FoundPeakIndex(); + if( computedindex<0 ){computedindex = mPA->AnalyzeForPeak();} + + this->PaintBaselines(); + this->PaintPeakRanges(); + + return; +} + +void PeakAnaPainter::PaintBaselines() +{ + Double_t base = mPA->Baseline(); + Double_t hitline = base+mPA->BaselineSigmaScale()*mPA->BaselineSigma(); + if( mTheBaseLine==0 ){mTheBaseLine = new TLine(mPA->MinX(),base,mPA->MaxX(),base);} + else{ mTheBaseLine->SetY1(base); mTheBaseLine->SetY2(base);} + if( mTheHitLine==0 ){mTheHitLine = new TLine(mPA->MinX(),hitline,mPA->MaxX(),hitline);} + else{ mTheHitLine->SetY1(hitline); mTheHitLine->SetY2(hitline);} + mTheBaseLine->SetLineColor(kBlack); + mTheHitLine->SetLineColor(kViolet); + + mTheBaseLine->Paint(); + mTheHitLine->Paint(); +} + +void PeakAnaPainter::PaintFoundRange() +{ + Int_t computedindex = mPA->FoundPeakIndex(); + if( computedindex<0 ){computedindex = mPA->AnalyzeForPeak();} + if( computedindex == mPA->NPeaks() ){return;}//If no peak found then computed index is equal to number of peaks in peak vector + + TLine* sl = mPA->GetPeak(computedindex).GetStartLine(mPA->MinY(),mPA->MaxY()); + sl->SetLineColor(kRed); + sl->Paint(); + TMarker* mp = mPA->GetPeak(computedindex).GetPeakMarker(); + mp->SetMarkerColor(kViolet); + mp->SetMarkerSize(mPA->GetMarkerSize()*2.0); + mp->Paint(); + TLine* el = mPA->GetPeak(computedindex).GetEndLine(mPA->MinY(),mPA->MaxY()); + el->SetLineColor(kOrange); + el->Paint(); +} + +void PeakAnaPainter::PaintPeakRanges( ) +{ + Int_t computedindex = mPA->FoundPeakIndex(); + if(computedindex<0 ){computedindex = mPA->AnalyzeForPeak();} + + if( mPA->GetDebug() > 1){ + LOG_DEBUG<< "|SizePeaks:"<NPeaks() << "|FoundPeak:"<Baseline() << "|Hit:"<Baseline()+mPA->BaselineSigmaScale()*mPA->BaselineSigma(); + if( mPA->NPeaks()!=0 && computedindexNPeaks() ){(mPA->GetPeak(computedindex)).Print("debug");} + LOG_DEBUG << endm; + } + + for( UShort_t ipeak = 0; ipeakNPeaks(); ++ipeak ){ + if( ipeak==computedindex ){continue;} + mPA->GetPeak(ipeak).GetStartLine( mPA->MinY(), mPA->MaxY() )->Paint(); + TMarker* mp = mPA->GetPeak(ipeak).GetPeakMarker(); + mp->SetMarkerSize(mPA->GetMarkerSize()*2.0); + mp->Paint(); + mPA->GetPeak(ipeak).GetEndLine( mPA->MinY(), mPA->MaxY() )->Paint(); + } + this->PaintFoundRange(); + + return; +} + +void PeakAnaPainter::PaintStats() +{ + //@[April 8, 2022]>Need a way to determine size of pavetext from number of peaks as well as information to be added for each peak (this could be as simple as counting the number of lines and scaling the size accordingly) + this->MakePaveText()->SetTextSize(0.025); + mPA->AddPeakStats( mPaveT_PA,mStatsOption.Data() ); + //mPaveT_PA->Draw("same"); + mPaveT_PA->Paint( mPaveT_PA->GetOption() ); + +} + diff --git a/StRoot/StFcsWaveformFitMaker/PeakAnaPainter.h b/StRoot/StFcsWaveformFitMaker/PeakAnaPainter.h new file mode 100644 index 00000000000..97a1451f4e1 --- /dev/null +++ b/StRoot/StFcsWaveformFitMaker/PeakAnaPainter.h @@ -0,0 +1,108 @@ +/* +Author: David Kapukchyan +@[September 29, 2022] +> Got rid of the virtual painter + +@[June 12, 2022] +> Painter now sets the found peak line color + +@[May 10, 2022] +> Changed how options work. Graph options now supported. Use semicolon separated list for graph, peak and stat options respectively. Added style options for the various lines. Moved lines for the peak windows to the 'PeakWindow' class. + +@[March 23, 2022] +> Found a bug where calling "Draw" on ROOT objects is bad because it adds the TLine objects created by this class to the pad. This means the pad can delete those objects which causes seg fault errors. The proper method to call is "Paint". This is because "Draw" for most ROOT objects only adds objects to the TPad's list so that later the pad can call the "Paint" method of the added class. This means when "Draw" is called on a *PeakAna* object it gets added to the pad so that "Paint" can be called. The "Paint" method will now draw all the lines needed. However, if "Draw" is called the lines are also added to the TPad's list. This is not needed since it is sufficient that only the *PeakAna* class be in the list so that the lines are not deleted by TPad. + + If that paint method involves a "Draw" call then it only adds objects to the pad without calling their paint method + +@[March 14, 2022] +> Separated the peak stats box painting method and added options for it + +@[March 9, 2022] +> This class will be used to draw PeakAna. It is a first try after reading ROOT code and and seeing that by having a separate painter class is better due to how TCanvas objects work and drawing in ROOT in general works. Also after splitting up the various classes I realized that it really is easier to have a separate painter and analysis class. The painter class essentially assumes a canvas/pad exists and only draws on it the things it needs independent of what else is on the canvas/pad. This is the philosophy to keep in mind when writing painter classes. +*/ + +/*! +Painter class for PeakAna +*/ + +#ifndef PEAKANAPAINTER_H +#define PEAKANAPAINTER_H + +//C++ Headers +#include +#include //std::move (Needs C++ 11) and std::pair + +//ROOT Headers +#include "TObject.h" +#include "TKey.h" +#include "TH1.h" +#include "TPaveText.h" + +//Custom Headers +#include "PeakWindow.h" + +class PeakAna; + +class PeakAnaPainter +{ +public: + PeakAnaPainter(); + virtual ~PeakAnaPainter(); + + virtual void Paint(Option_t *opt=""); + virtual void PaintRawData(); //!< Raw data with no modifications + virtual void PaintFoundPeak(); //!< Raw data inside zoomed in on found signal region + virtual void PaintFoundPeakQa(); //!< Draw signal and found signal window + virtual void PaintPeakQa(); //!< Show all found signal windows and signal + virtual void PaintBaselines(); //!< Just draw the baseline and hitlines + virtual void PaintFoundRange(); //!< Just draw the found peak on the current pad + virtual void PaintPeakRanges(); //!< Draw all found peaks on the current pad + virtual void PaintStats(); //!< Draw Stats box for peak finding + + virtual void CleanPainter(); //!< Clean up internal objects + + virtual void SetPeakAna(PeakAna* ana); //!< @param ana set the #PeakAna object to paint + + virtual Color_t GetBaseLineColor(); + virtual Style_t GetBaseLineStyle(); + virtual Width_t GetBaseLineWidth(); + + virtual Color_t GetHitLineColor(); + virtual Style_t GetHitLineStyle(); + virtual Width_t GetHitLineWidth(); + + virtual void SetBaseLineColor(Color_t color); + virtual void SetBaseLineColorAlpha(Color_t color,Float_t alpha); + virtual void SetBaseLineStyle(Style_t style); + virtual void SetBaseLineWidth(Width_t width); + + virtual void SetHitLineColor(Color_t color); + virtual void SetHitLineColorAlpha(Color_t color,Float_t alpha); + virtual void SetHitLineStyle(Style_t style); + virtual void SetHitLineWidth(Width_t width); + + bool ValidGraph(); //!< Check if #PeakAna object loaded and has a non-zero TGraph + +protected: + void Init();//!< Initialize internal variables to null + + TLine* mTheBaseLine; //!< line for the #PeakAna::mBaseline + TLine* mTheHitLine; //!< threshold for a peak #PeakAna::mBaseline + #PeakAna::mBaselineSigma*#PeakAna::mBaselineSigmaScale + + virtual TPaveText* MakePaveText(Double_t xmin=0.7,Double_t ymin=0.5, Double_t xmax=1.0, Double_t ymax=1.0); //!< Makes the stats box to show peak infomation + + PeakAna* mPA; //!< pointer to #PeakAna for drawing (PA=PeakAna) + TPaveText* mPaveT_PA; //!< for custom stats box + TString mGraphOption; //!< option for drawing the TGraph + TString mPeakOption; //!< option for drawing the peaks + TString mStatsOption; //!< option for what to put in stats box + + //[March 23, 2022]>Need a variable to know if something needs to repainted or not this will prevent multiple calls to delete?? + +private: + UInt_t mDEBUG; //!< debug level, 0=none + + ClassDef(PeakAnaPainter,2); +}; + +#endif diff --git a/StRoot/StFcsWaveformFitMaker/PeakWindow.cxx b/StRoot/StFcsWaveformFitMaker/PeakWindow.cxx new file mode 100644 index 00000000000..3159b32982d --- /dev/null +++ b/StRoot/StFcsWaveformFitMaker/PeakWindow.cxx @@ -0,0 +1,401 @@ +#include "PeakWindow.h" + +ClassImp(PeakWindow) + +PeakWindow::PeakWindow():mStartX(0),mEndX(1),mStartY(0),mEndY(0),mP_Peak(-1),mPeakX(-1),mPeakY(-1),mStartLine(0),mPeakMarker(0),mEndLine(0) +{ +} + +PeakWindow::PeakWindow(Double_t start, Double_t end):PeakWindow() +{ + if( startGetStartLine( mStartLine->GetY1(), mStartLine->GetY2()); } + if( mPeakMarker!=0 ){ window->GetPeakMarker(); } + if( mEndLine!=0 ) { window->GetEndLine( mEndLine->GetY1(), mStartLine->GetY2()); } + + return window; +} + +PeakWindow::~PeakWindow() +{ + delete mStartLine; + delete mPeakMarker; + delete mEndLine; +} + +void PeakWindow::SetWindow(Double_t s, Double_t e) +{ + mStartX=s; + mEndX=e; +} + +void PeakWindow::GetWindow(Double_t &s, Double_t &e) const +{s=mStartX;e=mEndX;} + +void PeakWindow::Reset(Double_t start, Double_t end) +{ + if( start=gdata->GetN()){return;} + Double_t x2=0; Double_t y2=0; + gdata->GetPoint(mP_Peak,x2,y2); + if( mP_Peak==0 || mP_Peak==gdata->GetN()-1 ){ + mPeakX = x2; + mPeakY = y2; + return; + } + Double_t x1=0; Double_t y1=0; + Double_t x3=0; Double_t y3=0; + gdata->GetPoint(mP_Peak-1,x1,y1); + gdata->GetPoint(mP_Peak+1,x3,y3); + Double_t x00 = (x1+x2)/2.0; + Double_t x11 = (x2+x3)/2.0; + Double_t PosSlope = (y2-y1)/(x2-x1); + Double_t NegSlope = (y3-y2)/(x3-x2); + mPeakX = (x00*NegSlope-x11*PosSlope)/(NegSlope-PosSlope);//Formula for x-intercept of a line with points (x00,PosSlope) and (x11,NegSlope) + mPeakY = y2; + + return; +} + +void PeakWindow::Combine( const PeakWindow &other, bool keepthis ) +{ + //Assuming mStartXother.mStartX ){ + mStartX = other.mStartX; + mStartY = other.mStartY; + } + //Same assumption as above then rightmost end value should be the greater of the two end values + if( mEndXGetPoint(mP_Peak,peakx,peaky); + return mslope*peakx+yint; + } + else{return mslope*mPeakX+yint;} +} + +Double_t PeakWindow::SlopeChirality(Double_t scale) const +{ + //[Feb 28, 2022]>Use small fluctation around zero to return 0?? + Double_t slope = StartEndLineSlope(); + if( slope<0 ){ return -1.0*cosh(scale*slope); } + else{ return cosh(scale*slope); } +} + +Double_t PeakWindow::PeakChirality(Double_t slopescale, Double_t scale) const +{ + if( scale==0 ){scale=1;} + if( scale<0 ){ scale = fabs(scale); } + Double_t startendslope = StartEndLineSlope(); + Double_t peakstartslope = (mPeakY-mStartY)/(mPeakX-mStartX); + Double_t peakendslope = (mPeakY-mEndY)/(mPeakX-mEndX); + + Double_t deltaps = peakstartslope-startendslope; + Double_t deltape = peakendslope-startendslope; + //scale==1 peak center is exactly halfway between StartX and EndX + //if scale>1 peak center i.e. Chirality==0 is towards StartX + //if scale<1 peak center i.e. Chirality==0 is towards EndX + return this->SlopeChirality(slopescale*startendslope)*((deltape*deltape)-(deltaps*deltaps)/(scale*scale)); + //return this->SlopeChirality(slopescale*startendslope); +} + +Double_t PeakWindow::PeakChiralityProb(Double_t probscale,Double_t chirality) const +{ + if( probscale<0 ){ probscale = fabs(probscale); } + //Probability is 1/(chir^2+1). This was chosen because 1/x^2+1 "falls softer" than e^-|x|. This returns a probabilty to NOT merge a peak meaning if the chirality is 0 then should return 1 and when chirality is +-inf should return 0. + return static_cast(1)/(probscale*chirality*chirality+static_cast(1)); +} + +Double_t PeakWindow::PeakChiralityProb(Double_t probscale, Double_t peakscale, Double_t chirscale ) const +{ + Double_t chir = PeakChirality(peakscale, chirscale); + return this->PeakChiralityProb(peakscale,chir); +} + +Double_t PeakWindow::PeakTunnelProb(TGraph* graph, Double_t scale, Double_t sigma) const +{ + if( sigma<0 ){sigma = fabs(sigma);} + else if(sigma==0){LOG_ERROR << "PeakWindow::PeakTunnelProb - ERROR:sigma cannot be 0 - Returning probability of 0" << endm; return 0;} + Double_t xdiff = mEndX-mStartX; + Double_t ydiff = mEndY-mStartY; + Double_t mslope = ydiff/xdiff; //slope of the line formed by the point (mStartX,mStartY) and (mEndX,mEndY) + Double_t yint = mStartY-mslope*mStartX;//y-intercept of the line above "^" + Double_t peakx=0; Double_t peaky=0; + graph->GetPoint(mP_Peak,peakx,peaky); + Double_t yline = mslope*peakx+yint;//y-value at peak x value, use mPeak instead(may not be set so check and pick)?? + Double_t heightdiff = peaky-yline; + //return 1.0/scale*exp(-1.0*scale*heightdiff*xdiff);//Alternative probablity that is not so easy to normalize + if( heightdiff<=0 ){ return 1; }//if peak value is less than the line value automatically tunnel + else{ + //return TMath::Exp(-1.0*fabs(scale)*fabs(xdiff)) * (TMath::Erfc((heightdiff)/(TMath::Sqrt2()*sigma))); + return (1.0/(scale*xdiff*xdiff+1.0))*(TMath::Erfc((heightdiff)/(TMath::Sqrt2()*sigma))); + } + //The idea for this probablity distribution (not probability density function) comes from two assumptions about the underlying data + //1. As the distance between the start and end x-values increases the probability of tunneling through will decrease exponentially with some scale that must be determined based on the input data. One-tenth of the difference in x-values may work. + //2. There is some noise in the y-value of the data that follows a Normal distribution. The mean ('yline') will be the y-value of the line connecting the start and end points evaluated at the x-value of the peak. The sigma should be the noise level of the data. The cumulative distribution function (which is the actual probability) for a Normal Distribution is the complimentary error function (Erfc). Note that the normalization constant of 1/2 has been removed since this is a probability distribution not a probability density function and hence does not need to be normalized as long as the value returned is between 0 and 1. This is true for Erfc as long as 'peaky'>='yline' which is the case here. + //The total proability of tunneling is then (@1)*(@2) with two free parameters the scale and the sigma. Note the scale and sigma must be positive +} + +UShort_t PeakWindow::CompareTo(const PeakWindow& other) const +{ + UShort_t check = 0; + if( mStartX == other.mStartX ){ + if( mEndX == other.mEndX ){ check = 1; } + } + if( check>0 && mP_Peak == other.mP_Peak ){check=2;} + if( check>1 && mStartY == other.mStartY ){ + if( mEndY == other.mEndY ){ + check = 3; + } + } + if( check>2 && mPeakX == other.mPeakX ){ check = 4;} + if( check>3 && mPeakY == other.mPeakY ){ check = 5;} + return check; +} + +void PeakWindow::Draw(Option_t* opt) +{ + AppendPad(opt); +} + +void PeakWindow::Paint(Option_t* opt) +{ + TString option(opt); + option.ToLower(); + bool drawstart = false; + bool drawpeak = false; + bool drawend = false; + + if( option.Length()==0 ){ + drawstart=true; + drawpeak=true; + drawend=true; + } + else{ + if( option.Contains("s") ){ drawstart = true; } + if( option.Contains("p") ){ drawpeak = true; } + if( option.Contains("e") ){ drawend = true; } + } + + if( drawstart ){GetStartLine()->Paint();} + if( drawpeak ){GetPeakMarker()->Paint();} + if( drawend ){GetEndLine()->Paint();} +} + +TLine* PeakWindow::GetStartLine(Double_t ymin, Double_t ymax) +{ + if( mStartLine==0 ){ + if( ymin==ymax ){ + if( mStartYSetLineColor(kGreen); + } + else{ + mStartLine->SetX1(mStartX); + mStartLine->SetX2(mStartX); + if( ymin==ymax ){ return mStartLine; } + mStartLine->SetY1(ymin); + mStartLine->SetY2(ymax); + } + return mStartLine; +} +Color_t PeakWindow::GetStartLineColor() const +{ if( mStartLine!=0 ){return mStartLine->GetLineColor(); } return 0; } +Style_t PeakWindow::GetStartLineStyle() const +{ if( mStartLine!=0 ){return mStartLine->GetLineStyle(); } return 0; } +Width_t PeakWindow::GetStartLineWidth() const +{ if( mStartLine!=0 ){return mStartLine->GetLineWidth(); } return 0; } +void PeakWindow::SetStartLineColor(Color_t color){ GetStartLine()->SetLineColor(color); } +void PeakWindow::SetStartLineColorAlpha(Color_t color,Float_t alpha){ GetStartLine()->SetLineColorAlpha(color,alpha); } +void PeakWindow::SetStartLineStyle(Style_t style){ GetStartLine()->SetLineStyle(style); } +void PeakWindow::SetStartLineWidth(Width_t width){ GetStartLine()->SetLineWidth(width); } + +TMarker* PeakWindow::GetPeakMarker() +{ + if( mPeakMarker==0 ){ + mPeakMarker = new TMarker(mPeakX,mPeakY,4);//Default style is open circle + mPeakMarker->SetMarkerColor(kGreen+1);//Default color is the green color that is between the green colors of start and end + } + else{ + mPeakMarker->SetX(mPeakX); + mPeakMarker->SetY(mPeakY); + } + return mPeakMarker; +} +Color_t PeakWindow::GetPeakMarkerColor() const +{ if( mPeakMarker!=0 ){return mPeakMarker->GetMarkerColor(); } return 0; } +Style_t PeakWindow::GetPeakMarkerStyle() const +{ if( mPeakMarker!=0 ){return mPeakMarker->GetMarkerStyle(); } return 0; } +Size_t PeakWindow::GetPeakMarkerSize() const +{ if( mPeakMarker!=0 ){return mPeakMarker->GetMarkerSize(); } return 0; } +void PeakWindow::SetPeakMarkerColor(Color_t color){ GetPeakMarker()->SetMarkerColor(color); } +void PeakWindow::SetPeakMarkerColorAlpha(Color_t color, Float_t alpha){ GetPeakMarker()->SetMarkerColorAlpha(color,alpha); } +void PeakWindow::SetPeakMarkerStyle(Style_t style){ GetPeakMarker()->SetMarkerStyle(style); } +void PeakWindow::SetPeakMarkerSize(Size_t size){ GetPeakMarker()->SetMarkerSize(size); } + +TLine* PeakWindow::GetEndLine(Double_t ymin, Double_t ymax) +{ + if( mEndLine==0 ){ + if( ymin==ymax ){ + if( mStartYSetLineColor(kGreen+2); + } + else{ + mEndLine->SetX1(mEndX); + mEndLine->SetX2(mEndX); + if( ymin==ymax ){ return mEndLine; } + mEndLine->SetY1(ymin); + mEndLine->SetY2(ymax); + } + return mEndLine; +} +Color_t PeakWindow::GetEndLineColor() const +{ if( mEndLine!=0 ){ return mEndLine->GetLineColor(); } return 0; } +Style_t PeakWindow::GetEndLineStyle() const +{ if( mEndLine!=0 ){ return mEndLine->GetLineStyle(); } return 0; } +Width_t PeakWindow::GetEndLineWidth() const +{ if( mEndLine!=0 ){ return mEndLine->GetLineWidth();} return 0; } +void PeakWindow::SetEndLineColor(Color_t color){ GetEndLine()->SetLineColor(color); } +void PeakWindow::SetEndLineColorAlpha(Color_t color,Float_t alpha){ GetEndLine()->SetLineColorAlpha(color,alpha); } +void PeakWindow::SetEndLineStyle(Style_t style){ GetEndLine()->SetLineStyle(style); } +void PeakWindow::SetEndLineWidth(Width_t width){ GetEndLine()->SetLineWidth(width); } + diff --git a/StRoot/StFcsWaveformFitMaker/PeakWindow.h b/StRoot/StFcsWaveformFitMaker/PeakWindow.h new file mode 100644 index 00000000000..f63ef6afceb --- /dev/null +++ b/StRoot/StFcsWaveformFitMaker/PeakWindow.h @@ -0,0 +1,195 @@ +/* +Author: David Kapukchyan +@[June 27, 2022] +> Added a TMarker for the peak position in drawing PeakWindow objects. Also added options to the Draw() and Paint() methods for this object. + +@[June 26, 2022] +> Fixed SetPeak to work for the extreme cases where the peak point is the first or last point + +@[June 12, 2022] +> Implemented TObject's "Clone" and "Copy" methods in line with *PeakAna* for ease of copying these objects. + +@[May 10, 2022] +> Added TLine objects here for drawing. As well they're accessor functions and styling options. + +@[February 28, 2022] +> Added Sinh(*mse*) into chirality function as both a multiplicative factor and addition to the "hyperbola" as described in Feb. 25, 2022. It did not work as expected since the "hyperbola" form goes to zero but so does the Sinh(*mse*) and so as a multiplicative factor that means whether the slope is zero or the "hyperbola" is zero you get low chirality since it should not be a one or the either case but a mixture. The sum also did not work as this meant the probablity function maxes at 0.5 instead of 1. The bigger issue is that the "hyperbola" term was dominating the chirality and thus the Sinh(*mse*) was not contributing at all to the chirality. The main issue that arises when one takes away the "hyperbola" part, is that the slope is still not sufficient as noisy data can easily cause the slope to be zero since the start and end points will match up. This is a "false" matching in the sense that the overall behavior of the curve could still be increasing but a large noise level causes the next point to be as low as the start point. In order to solve this one may need to take into account the noise level (e.g. "sigma" in tunnel scale) and treat it as some kind of error on the start and end points and then propagate that to the uncertainty in slope and use that as a part of the chirality ?? + +@[February 25, 2022] +> Testing of "peak" chirality function of the form (mpe - mse)^2 - (mps - mse)^2/(scale)^2; essentially a hyperbola where *mpe* is slope of line containing points (mPeakX,mPeakY) and (mEndX,mEndY), *mps* is slope of line containing points (mPeakX,mPeakY) and (mStartX,mStartY), *mse* is slope of line containing (mEndY,mEndX) and (mStartX,mStartY); in all cases order of points mattter; only one scale is needed since scale>1 shifts to the "start" and 0 Added a "peak" chirality function and functions to check if a given PeakWindow should be merged based on chirality + +@[February 18, 2022] +> Moved *PeakWindow* to it's implementation file as it's become a more composite structure +*/ +/*! + Data class to hold properties of a found peak like the peak position and its start and end points. Mostly used as a helper class for #PeakAna. + + A "peak window" consists of 3 points: where a peak starts (first positive slope after negative slopes leading up to peak), where it plateaus (slope is zero i.e. the peak position), and where the peak ends (last negative slope before changing to positive slopes) +*/ + +#ifndef PEAKWINDOW_H +#define PEAKWINDOW_H + +//C++ Headers +#include +#include +#include //std::move (Needs C++ 11) + +//ROOT Headers +#include "TObject.h" +#include "TKey.h" +#include "TH1.h" +#include "TMarker.h" +#include "TGraph.h" +#include "TLine.h" +#include "TMath.h" + +#include "St_base/StMessMgr.h" + +class PeakWindow : public TObject{ +public: + /**@brief Default Constructor + + Start and end points of peak are set to a range of 0 to 1 and #mP_Peak, #mPeakX, and #mPeakY equal to -1 (impossible point) + */ + PeakWindow(); + PeakWindow(Double_t start, Double_t end); //!< Construct with known start and end points, peak gets set to imposible values + PeakWindow(const PeakWindow& oldpeak); //!< Copy Constructor + PeakWindow& operator=(const PeakWindow& rhs); //!< Assignment operator + virtual ~PeakWindow(); //!< Destructor + + virtual void Copy(TObject& obj) const; //!< Only copies variables, to copy TLines use #Clone() + virtual TObject* Clone(const char* newname="") const; //!< Clone whole object, name is irrelevant + + //To ease computation time I store the x,y values of the peak window/range rather than the graph point + Double_t mStartX; //!< x value for start of the peak window + Double_t mEndX; //!< x value for end of the peak window + Double_t mStartY; //!< y value associated with #mStartX + Double_t mEndY; //!< y value associated with #mEndX + Int_t mP_Peak; //!< Point Number of peak in a TGraph object (P for point), point is such that slope with previous point will be positive and next point will be negative + Double_t mPeakX; //!< x-value of peak position as determined by #SetPeak() + Double_t mPeakY; //!< y-value at #mP_Peak + + void SetWindow(Double_t s, Double_t e); //!< @param s set x-value for start of peak @param e set x-value for end of peak + void GetWindow(Double_t &s, Double_t &e) const; //!< @param s get x-value for start of peak @param e get x-value for end of peak + /**@brief sets #mPeakX based on #mP_Peak using line of slopes from points (#mP_Peak-1,#mP_Peak) and (#mP_Peak,#mP_Peak+1) + + This function is used to set #mPeakX to a value from the left and right slopes of #mP_Peak to correct for any discretization coming from points on a graph. Requires #mP_Peak has been set correctly + @param gdata TGraph where data is stored + */ + void SetPeak(TGraph* gdata); + + /**@brief combine two #PeakWindow objects + + @param leftpeak one of the peaks to be combined and assumed to have the lower x-value + @param rightpeak one of the peaks to be combined and assumed to have the higher x-value + @param keepleft boolean to determine which peak position should be kept in the combined peak, true means "leftpeak", false means "rightpeak" + @return combined #PeakWindow object is returned + */ + static PeakWindow Combine(const PeakWindow &leftpeak, const PeakWindow &rightpeak, bool keepleft=true); + /**@brief merges one #PeakWindow into another + + Changes "this" peak to include the "other" peak. + @param other #PeakWindow to merge with this one + @param keepthis true means keep "this" peak's position, false means keep "other" peak's position + */ + virtual void Combine( const PeakWindow &other, bool keepthis=true ); + + /**@brief compare two #PeakWindow objects + + Compares two #PeakWindow's and returns a value based on how different they are + @return comparison flag\n + 0 = different #PeakWindow's\n + 1 = same #mStartX and #mStartY only\n + 2 = "1" + same #mP_Peak only\n + 3 = "1" + "2" + same #mStartY and #mEndY only\n + 4 = "1" + "2" + "3" + same #mPeakX only\n + 5 = same #PeakWindow's ("1"+"2"+"3"+"4"+same #mPeakY) + */ + virtual UShort_t CompareTo( const PeakWindow& other ) const; + Double_t StartEndLineSlope() const; //!< Computes the slope of the line formed by the points (#mStartX,#mStartY) and (#mEndX,#mEndY) + Double_t StartEndSlopeUncertainty(Double_t sigma) const; //!< Uncertainty int the slope of the line formed by the points (#mStartX,#mStartY) and (#mEndX,#mEndY) + Double_t StartEndLineYint() const; //!< Computes the y-intercept of the line formed by the points (#mStartX,#mStartY) and (#mEndX,#mEndY) + + /**@brief Computes the the line formed by the points (#mStartX,#mStartY) and (#mEndX,#mEndY) and evaluates that line at #mPeakX + + The function is mostly needed in computing the probablity for peak tunneling since the difference of #MidPoint() and #mPeakY is used in the probability formula. + @param graph if graph given use that graph's x-value at #mP_Peak rather than #mPeakX + @return y-value evaluated at #mPeakX of line formed by points (#mStartX,#mStartY) and (#mEndX,#mEndY) + */ + Double_t MidPoint(TGraph* graph=0) const; + virtual Double_t SlopeChirality(Double_t scale) const;//function to compute chirality factor for startendslope + virtual Double_t PeakChirality(Double_t slopescale, Double_t peakscale) const;//for peakscale 1 means peak is centered in window, peakscale>1 peak center shifted towards start, peakscale<1 peak center shifted towards end (peakscale>0) + virtual Double_t PeakChiralityProb(Double_t probscale, Double_t chirality) const; + virtual Double_t PeakChiralityProb(Double_t probscale, Double_t peakscale, Double_t chirscale) const;//Returns probabilty not to merge or probabilty is a "real" peak + + /**@brief Compute probablity that a given #PeakWindow is a real peak + + \f{equation}{ Probability = 1/(scale*StartEndDiff^2+1)*Erfc(PeakHeightDiff/(sqrt(2)*sigma)) \f} or 1 if PeakHeightDiff<=0\n + Erfc = complimentary error function\n + StartEndDiff = #mEndX-#mStartX\n + PeakHeightDiff = #mPeakY - #MidPoint()\n + + This function will work even if #SetPeak() is not called since it requires a TGraph and will read the values from there + + @param graph TGraph of data points + @param scale StartEndDiff scale in formula of probability + @param sigma sigma of Erfc to use in formula of probability + */ + virtual Double_t PeakTunnelProb(TGraph* graph, Double_t scale=1.0, Double_t sigma=1.0 ) const; + //bool ChirMerge(const Rtools::PeakWindow& other, Double_t scalechir=1) const;//returns true if "this" should be merged with "other" PeakWindow according to PeakChirality + //bool ChirMergeLeft(const Rtools::PeakWindow& left, const Rtools::PeakWindow& right, Double_t scalechir=1) const;//based on chirality - returns true if should be merged with "left" PeakWindow; false if should be merged with "right" PeakWindow + virtual void Reset(Double_t start, Double_t end); //!< Reset #PeakWindow to constructor state + virtual void Print(Option_t* opt="") const; //!< Prints information about PeakWindow + + /**@brief Draw the PeakWindow + + @param opt options for drawing: + - "" (empty string) means draw start line, peak marker and end line.\n + - "s" means draw start line\n + - "p" means draw peak marker\n + - "e" means draw end line + */ + virtual void Draw(Option_t* opt=""); + virtual void Paint(Option_t* opt=""); //!< paint method see #Draw() for options + + //Options for drawing + TLine* GetStartLine(Double_t ymin=0, Double_t ymax=0); //!< Create and return a TLine for the start of the peak window + Color_t GetStartLineColor() const; + Style_t GetStartLineStyle() const; + Width_t GetStartLineWidth() const; + void SetStartLineColor(Color_t color); + void SetStartLineColorAlpha(Color_t color,Float_t alpha); + void SetStartLineStyle(Style_t style); + void SetStartLineWidth(Width_t width); + + TMarker* GetPeakMarker(); //!< Create and return a TMarker to mark the location of the peak + Color_t GetPeakMarkerColor() const; + Style_t GetPeakMarkerStyle() const; + Size_t GetPeakMarkerSize() const; + void SetPeakMarkerColor(Color_t color); + void SetPeakMarkerColorAlpha(Color_t color, Float_t alpha); + void SetPeakMarkerStyle(Style_t style); + void SetPeakMarkerSize(Size_t size); + + TLine* GetEndLine(Double_t ymin=0, Double_t ymax=0); //!< Create and return a TLine for the end of the peak window + Color_t GetEndLineColor() const; + Style_t GetEndLineStyle() const; + Width_t GetEndLineWidth() const; + void SetEndLineColor(Color_t color); + void SetEndLineColorAlpha(Color_t color,Float_t alpha); + void SetEndLineStyle(Style_t style); + void SetEndLineWidth(Width_t width); + +protected: + TLine* mStartLine; //!< TLine for drawing the start of the peak window + TMarker* mPeakMarker; //!< TMarker for drawing the peak location + TLine* mEndLine; //!< TLine for drawing the end of the peak window + + ClassDef(PeakWindow,3); +}; + +#endif diff --git a/StRoot/StFcsWaveformFitMaker/StFcsPulseAna.cxx b/StRoot/StFcsWaveformFitMaker/StFcsPulseAna.cxx new file mode 100644 index 00000000000..d6aedabb3d9 --- /dev/null +++ b/StRoot/StFcsWaveformFitMaker/StFcsPulseAna.cxx @@ -0,0 +1,700 @@ +#include "StFcsPulseAna.h" + +ClassImp(StFcsPulseAna) + +double StFcsPulseAna::MaxwellBoltzmannDist(double* x, double* p) +{ + if( x[0]Clone( (mName+"_H1_Basline").c_str() );} + if( old.BaselineFit()!=0 ){mF1_BaselineFit = (TF1*)old.BaselineFit()->Clone( (mName+"_F1_BaselineFit").c_str() );} + if( old.SignalFit()!=0 ){mF1_SignalFit = (TF1*)old.SignalFit()->Clone( (mName+"_F1_SignalFit").c_str() );} +} + +StFcsPulseAna& StFcsPulseAna::operator=(const StFcsPulseAna& rhs) +{ + if( this == &rhs ){ return *this; } + //Note that the graph is not cloned because PeakAna should not be changing the contents of the data + PeakAna::operator = (rhs); + mName = rhs.Name()+"_copy"; + + if( rhs.BaselineHist()!=0 ){mH1_Baseline = (TH1F*)rhs.BaselineHist()->Clone( (mName+"_H1_Basline").c_str() );} + if( rhs.BaselineFit()!=0 ){mF1_BaselineFit = (TF1*)rhs.BaselineFit()->Clone( (mName+"_F1_BaselineFit").c_str() );} + if( rhs.SignalFit()!=0 ){ mF1_SignalFit = (TF1*)rhs.SignalFit()->Clone( (mName+"_F1_SignalFit").c_str() ); } + + return *this; +} + +void StFcsPulseAna::Init() +{ + mDbPulse = 0; + + mWindowSum = false; mSumAdc=0; + mFitSum = false; mSumAdcFit=0.0; + + mH1_Baseline = 0; + mF1_BaselineFit = 0; + + mF1_SignalFit = 0; +} + +void StFcsPulseAna::Copy(TObject& obj) const +{ + PeakAna::Copy(obj); + + ((StFcsPulseAna&)obj).mDbPulse = mDbPulse; + ((StFcsPulseAna&)obj).mWindowSum = mWindowSum; + ((StFcsPulseAna&)obj).mSumAdc = mSumAdc; + ((StFcsPulseAna&)obj).mFitSum = mFitSum; + ((StFcsPulseAna&)obj).mSumAdcFit = mSumAdcFit; + + ((StFcsPulseAna&)obj).mH1_Baseline = 0; + ((StFcsPulseAna&)obj).mF1_BaselineFit = 0; + ((StFcsPulseAna&)obj).mF1_SignalFit = 0; +} + +TObject* StFcsPulseAna::Clone(const char* newname) const +{ + TGraph* cloneg = (TGraph*)this->GetData()->Clone(); + StFcsPulseAna* ana = 0; + + if( !strlen(newname) ){ ana = new StFcsPulseAna(cloneg,mName+"_copy"); } + else{ ana = new StFcsPulseAna(cloneg,newname); } + Copy(*ana); + ana->ForceInternal();//Cloned graph should be deleted by cloned object + + if( FoundPeakIndex()>=0 ){ + //If mComputedIndex>=0 this means AnalyzeForPeak() was called and need to copy the peaks. This is preffered over re-running the analysis as copy should be faster than re-running the data + ana->mComputedIndex = mComputedIndex; + for( UInt_t i=0; imPeaks.push_back( mPeaks.at(i) ); } + ana->mFoundPeak = mFoundPeak; + } + + if( mH1_Baseline!=0 ){ ana->mH1_Baseline = (TH1F*)mH1_Baseline->Clone((ana->Name()+"_H1_Baseline").c_str()); } + if( mF1_BaselineFit!=0 ){ ana->mF1_BaselineFit = (TF1*)mF1_BaselineFit->Clone((ana->Name()+"_F1_Baseline").c_str()); } + if( mF1_SignalFit!=0 ){ ana->mF1_SignalFit = (TF1*)mF1_SignalFit->Clone((ana->Name()+"_F1_SignalFit").c_str()); } + + return ana; +} + +StFcsPulseAna::~StFcsPulseAna() +{ + delete mH1_Baseline; + delete mF1_SignalFit; + delete mF1_BaselineFit; + +} + +void StFcsPulseAna::ResetFinder() +{ + PeakAna::SetData((TGraph*)0);//Effectively deletes and sets graph to zero and destroys found peaks + ResetBaseline(); + ResetSum(); +} + +void StFcsPulseAna::ResetBaseline() +{ + mBaseline = -5.0; + mBaselineSigma = 0.75; + if( mH1_Baseline!=0 ){ + mH1_Baseline->Reset(); + mH1_Baseline->GetXaxis()->SetRangeUser(-0.5,4095.5);//Needs to go back to default since calling reset doesn't revert axes to original values + } + if( mF1_BaselineFit!=0 ){mF1_BaselineFit->ResetAttLine();}//Values from old fit will be replaced when 'AnalyzeForBaseline' gets called +} + +void StFcsPulseAna::ResetSum() +{ + if( mF1_SignalFit!=0 ){ delete mF1_SignalFit; mF1_SignalFit=0; }//This was for plotting but still good to do just in case + + mWindowSum=false; mSumAdc = 0; + mFitSum=false; mSumAdcFit = 0.0; +} + +Int_t StFcsPulseAna::AnalyzeForPeak() +{ + mComputedIndex = PeakAna::AnalyzeForPeak(); + if( mTunnelThreshold<0 ){ + mTunnelThreshold *= -1.0; + if( mTunnelThreshold<=1.0 ){ + std::vector merged; + this->MergeByProbability(merged); + mPeaks.swap( merged ); + mTunnelThreshold *= -1.0;//Restore to old value + mComputedIndex = this->SearchForPeak( mPeaks ); + } + } + if( mComputedIndex(mPeaks.size()) ){ mFoundPeak = mPeaks.at(mComputedIndex); } + else{mFoundPeak.SetWindow(mXRangeMax+1,mXRangeMax+1);}//If no valid index then set found peak to greater than max X values so it is obvious something's wrong + + return mComputedIndex; +} + +void StFcsPulseAna::MergeByProbability(std::vector& merged) const +{ + if( !merged.empty() ){ merged.clear(); } + for( UInt_t i=0; i 2.1 ){ merged.push_back(win); continue;}//Skip merging peaks that are more than 2 timebins away + + if( PeakTunnel(win) ){ + if( merged.empty() ){ merged.push_back(win); continue; }//vector is empty so add first found peak + Double_t thisprob = merged.back().PeakTunnelProb(mG_Data,mTunnelScale,mTunnelSigma); + Double_t nextprob = win.PeakTunnelProb(mG_Data,mTunnelScale,mTunnelSigma); + merged.back().Combine(win,thisprob<=nextprob?true:false);//Take peak with highest probability + } + else{ merged.push_back(win); } + } +} + +void StFcsPulseAna::AnalyzeForPedestal() +{ + ResetBaseline();//Resets histogram and function if they exist + if( mH1_Baseline==0){mH1_Baseline = new TH1F( (mName+"_"+"H1_Baseline").c_str(),"",4096,-0.5,4095.5); mH1_Baseline->Sumw2(); mH1_Baseline->SetTitle("Baseline Histogram");} + //Fill histogram of ADC values + for( int ipoint=0; ipointGetData()->GetN(); ++ipoint) + { + Double_t tb; Double_t adc; + Int_t check = this->GetData()->GetPoint(ipoint,tb,adc); + if( check<0 ){ LOG_WARN << "Unable to find point " << ipoint << " in Signal graph" << endm; continue;} + if( adc<0.1 ){continue;} //skip adc=0 from histogram since this is most likely missing data in non-pedestal subtracted data + mH1_Baseline->Fill(adc); + } + /*May 19, 2021: Fitting with too few data points gives bad baselines so just use max of histogram for now*/ + mBaseline = mH1_Baseline->GetMaximumBin()-1;//Offset by one to account for bin counting + mBaselineSigma = 0.75; + + //Create a function to fit to histogram of ADC values + if( mF1_BaselineFit==0 ){ mF1_BaselineFit = new TF1((mName+"_"+"F1_BaselineFit").c_str(),"gaus(0)", + mH1_Baseline->GetXaxis()->GetXmin(), + mH1_Baseline->GetXaxis()->GetXmax() + ); + } + if( FindBaseline() ){return;} + else + { + LOG_WARN << "Unable to find a proper baseline" << endm; + if( GetDebug()>2 ){ LOG_DEBUG << "->";Print();} + else{LOG_INFO << "->Found baseline:"<< mBaseline << " Found BaselineSigma:"<Setting baseline to 0 and sigma to 0.75 to prevent T0 algorithm from failing" << endm;} + mBaseline=0; + mBaselineSigma=0.75; + return; + } +} + +bool StFcsPulseAna::FindBaseline() +{ + if( mBaseline > -0.001 ) + { + if( GetDebug()>2 ) + { + LOG_DEBUG << "|BaselineValue:"<GetEntries()==0 ){mBaseline = -5.0; return false;} + //For nonpedestal data use full range for pedestal data find largest bin and use a +-3 bin buffer for the fit + Int_t XStartVal = mH1_Baseline->GetMaximumBin()-1;//Offset by one to make the value correct + if( GetDebug()>2 ){ LOG_DEBUG << "|Height:"<GetBinContent(XStartVal+1) << "|XStartVal:"<< XStartVal << "|Range:"<<3 <<"|StartSigma:"<<0.9 << endm; } + mF1_BaselineFit->SetRange(XStartVal-3,XStartVal+3);//Only do within this range + mF1_BaselineFit->SetParameter(0,mH1_Baseline->GetBinContent(XStartVal+1)); + mF1_BaselineFit->SetParameter(1,XStartVal);//Use maximum instead of mean for more accurate center + mF1_BaselineFit->SetParameter(2,0.9);//Right now fixed but may need changing in the future? + + //Base line is mean of the gaussian fit + char opt[5] = "NRQ"; + if( GetDebug()>2 ){ opt[2] = '\0';}//Hack to to not quiet fit output + Int_t FitStatus = mH1_Baseline->Fit(mF1_BaselineFit,opt); + if( GetDebug()>1 ){ LOG_DEBUG << "|FitStatus:"<< FitStatus << endm; } + if( FitStatus == 0 ) + { + Double_t ratio = mF1_BaselineFit->GetChisquare()/static_cast(mF1_BaselineFit->GetNDF()); + if( ratio > 10.0 ){ mBaseline = XStartVal; return true; }//Take care of bad baseline fits by setting baseline to maximum + else + { + mBaseline = static_cast(mF1_BaselineFit->GetParameter(1)); + mBaselineSigma = static_cast( fabs(mF1_BaselineFit->GetParameter(2)) );//fabs in case sigma is negative + return true; + } + } + else if( FitStatus > 0 ) + { + //Max works better + mBaseline = mH1_Baseline->GetMaximumBin()-1;//Offset by one to account for bin counting + mBaselineSigma = 0.75; + if(GetDebug()>1){LOG_WARN << "Fit failed to converge setting baseline to mean of histogram and of sigma to 0.75" << endm;} + return true; + } + else{return false;} +} + +void StFcsPulseAna::FillAdc(TGraphAsymmErrors* g, unsigned short& counter, int Start, unsigned short* adcdata) +{ + //adcdata must be size 8 array + if( counter != 0 ){return;}//counter must start with zero size + unsigned int startpoint = 0; + double tb = 0; + double adc = 0; + for( unsigned int ipoint=0; static_cast(ipoint)GetN() && counter<8; ++ipoint ){ + g->GetPoint(ipoint,tb,adc); + if( int(tb)>=Start ){ startpoint=ipoint; break; } + } + + if( startpoint==0 || static_cast(startpoint+8)>=g->GetN() ){return;} + while( int(tb-(Start+counter)) < int(8-counter) ){ + counter = static_cast(tb-Start); + adcdata[counter]=adc; + g->GetPoint(++startpoint,tb,adc); + } +} + +//Copied on Nov. 01, 2021 +int StFcsPulseAna::SumDep0(TGraphAsymmErrors* gdata, int Start, int ped) +{ + unsigned short Size = 0;//Must start with zero size + unsigned short AdcData[8] = {0};//Must be 8 timebin array + StFcsPulseAna::FillAdc(gdata,Size,Start,AdcData); + if( Size==0 ){return 0;} + + //double ped = Baseline(); + int sum = 0 ; + int peak = 0 ; + int last = 0 ; + + for(int tb=0;tb<8;tb++) { + + unsigned short radc = AdcData[tb] ; + + switch(tb) { + case 0 : + last = radc ; + sum = radc ; + peak = 0 ; + break ; + case 1 : + if(radc>sum) peak = 1 ; + sum += radc ; + last = radc ; + break ; + case 2 : + sum += radc ; + last = radc ; + break ; + case 3 : + sum += radc ; + last = radc ; + break ; + case 4 : + sum += radc ; + last = radc ; + break ; + case 5 : + sum += radc ; + last = radc ; + break ; + case 6 : + sum += radc ; + last = radc ; + break ; + case 7 : + default : + //printf("radc %d, last %d, peak %d\n",radc,last,peak) ; + + if(radc>=last && peak==0) { + sum = 0 ; + } + else { + sum += radc ; + + sum -= ped*8 ; // ped is now only 3*ch_ped! + if(sum < 0) sum = 0 ; + } + + break ; + } + } + return sum; + +} + +int StFcsPulseAna::SumDep0Mod(TGraphAsymmErrors* gdata, int Start, int ped){ + unsigned short Size = 0;//Must start with zero size + unsigned short AdcData[8] = {0};//Must be 8 timebin array + StFcsPulseAna::FillAdc(gdata,Size,Start,AdcData); + if( Size==0 ){return 0;} + + //double ped = Baseline(); + int sum = 0 ; + int peak = 0 ; + int last = 0 ; + + for(int tb=0;tb<8;tb++) { + + unsigned short radc = AdcData[tb] ; + + switch(tb) { + case 0 : + last = radc ; + sum = radc ; + peak = 0 ; + break ; + case 1 : + if(radc>sum) peak = 1 ; + sum += radc ; + last = radc ; + break ; + case 2 : + sum += radc ; + last = radc ; + break ; + case 3 : + sum += radc ; + last = radc ; + break ; + case 4 : + sum += radc ; + last = radc ; + break ; + case 5 : + sum += radc ; + last = radc ; + break ; + case 6 : + if( peak!=0 || radc0 ){ sum=0; } + else { + sum += radc ; + + sum -= ped*8 ; // ped is now only 3*ch_ped! + if(sum < 0) sum = 0 ; + } + + break ; + } + } + return sum; + +} + +Int_t StFcsPulseAna::Sum(Int_t Start, Int_t End) +{ + Int_t SumAdc = 0; + Double_t base = Baseline();//Necessary to check if pedestal value has been set. + if( base < -4.0 ){return SumAdc;} + if( mG_Data==0 || mG_Data->GetN()==0 ){return SumAdc;} + for( Int_t ismp=0; ismpGetData()->GetN(); ++ismp ) + { + Double_t tb; Double_t adc; + GetData()->GetPoint(ismp,tb,adc); + if( tb < Start ){continue;} + if( tb > End ){ break; } + SumAdc += adc; + } + if( GetDebug()>1 ){LOG_DEBUG << "|Start:"<1 ){LOG_DEBUG << "|SumAdcA:"<< SumAdc << "|BaseArea:"<GetN()==0 ){return mSumAdcFit;} + if( Start==0 && End==0 ){Start=PeakStart(); End=PeakEnd();}//Default case + //Upon failure to find a window 'Start' and 'End' value will be equal. In this case just set some default window. For run 2019 it will be 35-60 and triggered crossing was at 50 + //For 2019 Ecal Cosmic running signal window was between 110-140 + //if( Start==End ){ SetWindow(110,140); Start=115; End=135; }//40 and 55 here since I add to range later + //Sanity Checks + if( !GoodWindow() ){ return mSumAdcFit; } + if( mF1_SignalFit!=0 ){ delete mF1_SignalFit; mF1_SignalFit=0; } + mF1_SignalFit = new TF1( (mName+"_Gaus_"+"mF1_SignalFit").c_str(),"gaus(0)+[3]",Start,End); + + //mG_Data->GetXaxis()->SetRangeUser(Start-5,End+5); + mF1_SignalFit->SetParameter( 0, MaxY() ); + mF1_SignalFit->SetParameter( 1, MaxX() ); + mF1_SignalFit->SetParameter( 2, static_cast(End-Start)*0.2 );//Seemed to work most of the time with 0.3, 0.5 seemed too much + mF1_SignalFit->SetParameter( 3, base ); + char opt[5] = "NRQ"; + if( GetDebug()>2 ){ opt[2] = '\0';}//Hack to to not quiet fit output + Int_t fitstatus = this->GetData()->Fit( mF1_SignalFit, opt ); + if( fitstatus >= 0 ) + { + //Fix to 4 sigma for now in the future make it a variable? + Double_t FitStart = mF1_SignalFit->GetParameter(1) - 4.0*fabs(mF1_SignalFit->GetParameter(2)); + Double_t FitEnd = mF1_SignalFit->GetParameter(1) + 4.0*fabs(mF1_SignalFit->GetParameter(2)); + mSumAdcFit = mF1_SignalFit->Integral(FitStart,FitEnd);//Sum only in the range and subtract the area under the baseline + if( GetDebug()>2 ){LOG_DEBUG << "|FitStart:"<2 ){LOG_DEBUG << "|FitSumA:"<GetN()==0 ){return mSumAdcFit;} + if( Start==0 && End==0 ){Start=PeakStart(); End=PeakEnd();}//Default case + //Upon failure to find a window 'Start' and 'End' value will be equal. In this case just set some default window. For run 2019 it will be 35-60 and triggered crossing was at 50 + //For 2019 Ecal Cosmic running signal window was between 110-140 + //if( Start==End ){ SetWindow(110,140); Start=115; End=135; }//40 and 55 here since I add to range later + //Sanity Checks + if( !GoodWindow() ){ return mSumAdcFit; } + if( mF1_SignalFit!=0 ){ delete mF1_SignalFit; mF1_SignalFit=0; } + mF1_SignalFit = new TF1( (mName+"_MB_"+"mF1_SignalFit").c_str(),MaxwellBoltzmannDist,Start,End,4); + //mG_Data->GetXaxis()->SetRangeUser(Start-5,End+5); + Double_t height=1; Double_t scale=1; + SignalMBPars(height,scale); + mF1_SignalFit->SetParameter( 0, height ); + mF1_SignalFit->SetParameter( 1, mFoundPeak.mStartX ); + mF1_SignalFit->SetParameter( 2, scale ); + mF1_SignalFit->SetParameter( 3, base ); + mF1_SignalFit->SetParName(0,"Amplitude");//Amplitude + mF1_SignalFit->SetParName(1,"X-Offset"); //x-offset (peak rise location) + mF1_SignalFit->SetParName(2,"Scale"); //scale parameter + mF1_SignalFit->SetParName(3,"Y-Offset"); //y-offset + mF1_SignalFit->FixParameter(3,base);//Works best when y-offset not allowed to change + + char opt[10] = "BNRQ"; + if( GetDebug()>2 ){ opt[3] = '\0';}//Hack to to not quiet fit output + Int_t fitstatus = this->GetData()->Fit( mF1_SignalFit, opt ); + if( fitstatus >= 0 ) + { + mSumAdcFit = mF1_SignalFit->Integral(Start,End);//Sum only in the range and subtract the area under the baseline + if( GetDebug()>2 ){LOG_DEBUG << "|FitStart:"<2 ){LOG_DEBUG << "|FitSumA:"<GetN()==0 ){return mSumAdcFit;} + if( Start==0 && End==0 ){Start=PeakStart(); End=PeakEnd();}//Default case + //Upon failure to find a window 'Start' and 'End' value will be equal. In this case just set some default window. + //Sanity Checks + if( !GoodWindow() ){ return mSumAdcFit; } + if( mF1_SignalFit!=0 ){ delete mF1_SignalFit; mF1_SignalFit=0; } + //Only care about one found peak for now in future can change to include more peaks ?? + int npeak = 1; + double para[5] = {static_cast(npeak),base,PeakY(),PeakX(),mDbPulse->GSigma()}; + mF1_SignalFit = new TF1( (mName+"_Pulse_"+"mF1_SignalFit").c_str(), mDbPulse,&StFcsDbPulse::multiPulseShape,Start,End,2+npeak*3); + mF1_SignalFit->SetParameters(para); + mF1_SignalFit->FixParameter(0,npeak); + mF1_SignalFit->FixParameter(1,base); + for(int i=0; iSetParLimits(1+i*3+1,0.0,40000.0); //limit peak not to go negative + int j=1+i*3+2; + mF1_SignalFit->SetParLimits(j,para[j]-2.0,para[j]+2.0); //limit peak position to +- 2TB + mF1_SignalFit->SetParLimits(1+i*3+3,0.01,10.0); //limit sigma to go too narrow or wide (Note(May 5, 2021):orig 0.5 to 10.0) + } + + char opt[10] = "BNRQ"; + if( GetDebug()>2 ){ opt[3] = '\0';}//Hack to to not quiet fit output + Int_t fitstatus = this->GetData()->Fit( mF1_SignalFit, opt ); + //Int_t fitstatus = 0;//For testing + if( fitstatus >= 0 ) + { + mSumAdcFit = mF1_SignalFit->Integral(Start,End);//Sum only in the range and subtract the area under the baseline + if( GetDebug()>2 ){LOG_DEBUG << "|FitStart:"<2 ){LOG_DEBUG << "|FitSumA:"<FixParameter(0,0); + func->FixParameter(1,0); + } + return; + } + if( FoundPeakIndex()<0 ){ this->AnalyzeForPeak(); } + int npeaks = NPeaks(); + if( func==0 ){ + func = mDbPulse->createPulse(mXRangeMin,mXRangeMax,2+npeaks*3); + } + if( npeaks>33 ){return;}//Since parameter array has 101 this is (101-2)/3 max peaks + double para[101] = {0}; + para[0] = npeaks; + para[1] = Baseline(); + func->SetParName(0,"NPeaks"); + func->SetParName(1,"Ped"); + for( Int_t i=0; iSetParName(j,name); + para[j++] = GetPeak(i).mPeakY; + sprintf(name,"P%d_M",i); + func->SetParName(j,name); + para[j++] = GetPeak(i).mPeakX; + sprintf(name,"P%d_S",i); + func->SetParName(j,name); + para[j] = mDbPulse->GSigma(); + } + func->SetParameters(para); + func->FixParameter(0,npeaks); + func->FixParameter(1,Baseline()); + for(int i=0; iSetParLimits(j++,0.0,50000.0); //limit peak not to go negative + func->SetParLimits(j,para[j]-2.0,para[j]+2.0); //limit peak position to +- 2TB + func->SetParLimits(++j,0.5,10.0); //limit sigma to go too narrow or wide + } +} + +StFcsPulseAna* StFcsPulseAna::DrawCopy(Option_t* opt, const char* name_postfix, TGraph* graph) const +{ + StFcsPulseAna* ana = new StFcsPulseAna(*this,name_postfix,graph);//Don't want to call "Clone" since that will clone the graph which you may not want to do in some cases. If graph==0 then it will be cloned + ana->SetBit(kCanDelete); + ana->AnalyzeForPeak();//Since vector of peaks is not copied re-run analysis; this is in case a filter was used which may replace the graph or a different graph is used + ana->AppendPad(opt);//Only append to pad don't draw yet + + return ana; +} + +void StFcsPulseAna::Print(Option_t* opt) const +{ + TString option(opt); + option.ToLower(); + std::cout << "|Name:"<GetData(); + std::cout << "|P::SignalFit:"< Added Refinements to drawing and copying. Also made some changes to how the tunnel probability is used to merge the peaks. + +@[March 3, 2022] + > Copied StFcsPulseFit so that I can inherit from and use the PeakAna class methods from MyTools +*/ + +/*! + This class is an extension of #PeakAna to analyze specifically STAR data from the DEP boards. The DEP boards are used in the Forward Calorimeter System (FCS) and comprises of an electromagnetic, and hadronic calorimeters as well as a presower. +*/ + +#ifndef STFCSPULSEANA_H +#define STFCSPULSEANA_H + +//ROOT Headers +#include "TMath.h" +#include "TLegend.h" +#include "TLegendEntry.h" +#include "TF1.h" +#include "TGraphAsymmErrors.h" +#include "TMultiGraph.h" +#include "TH1.h" + +#include "TLine.h" +#include "TPaveText.h" +#include "TPaveStats.h" + +#include "StFcsDbMaker/StFcsDbPulse.h" +#include "PeakAna.h" + +class StFcsPulseAna : public PeakAna +{ + public: + StFcsPulseAna(); //!< Default constructor + StFcsPulseAna( std::string name ); //!< Constructor with name + explicit StFcsPulseAna( TGraph* Sig, std::string name = "StFcsPulseAna"); //!< Construct with a known TGraph + StFcsPulseAna(const StFcsPulseAna& old, const char* post_name="_copy", TGraph* graph=0); //!< Copy Constructor can be called with a new graph and post-fix to the name + StFcsPulseAna& operator=(const StFcsPulseAna& rhs); //!< Assignment operator + + virtual ~StFcsPulseAna(); //!< Destructor + + virtual void Copy(TObject& obj) const; //!< Internal method use Clone instead. + virtual TObject* Clone(const char* newname) const; //!< Clones internal graph as opposed to just copying the pointer + + virtual Int_t AnalyzeForPeak(); //!< Overwritten from #PeakAna to process peak tunneling after finding all peaks + virtual Int_t AnalyzeForPeak(Double_t peak, Double_t width){ return PeakAna::AnalyzeForPeak(peak,width); }//& merged) const; //!< Overwritten from #PeakAna::MergeByProbability() to change merge criteria + + protected: + void Init(); //!< Initialize everything to zero except signal and background histograms + + StFcsDbPulse* mDbPulse; //!< Pointer to #StFcsDbPulse + + TH1F* mH1_Baseline; //!< Histogram that holds projection of mG_Data for baseline determination + //Main fit functions + TF1* mF1_SignalFit; //!< Function to fit to the data + TF1* mF1_BaselineFit; //!< Gaussian function to fit to #mH1_Baseline to determine baseline + + bool FindBaseline(); //!< Does Gaussian fitting to #mH1_Baseline to determine baseline + + private: + + std::string mName; //!< Name of class + + bool mWindowSum; //!< true if a sum using #SumWindow() was called + bool mFitSum; //!< true if any sum method that uses fitting was called + + //Variables to hold sum values. This is to avoid tedious recomputation + Int_t mSumAdc; //!< holds the value from #SumWindow() when it is called + Double_t mSumAdcFit; //!< holds the value from a fitting sum method when one is called + + ClassDef(StFcsPulseAna,1) + +}; + +#endif diff --git a/StRoot/StFcsWaveformFitMaker/StFcsPulseFit.cxx b/StRoot/StFcsWaveformFitMaker/StFcsPulseFit.cxx deleted file mode 100644 index 6a65fcfc68e..00000000000 --- a/StRoot/StFcsWaveformFitMaker/StFcsPulseFit.cxx +++ /dev/null @@ -1,523 +0,0 @@ -//Author:David Kapukchyan -//Oct 23, 2019 -//Class implmentation for signal finding in FCS - -#include "StFcsPulseFit.h" - -ClassImp(StFcsPulseFit) - -StFcsPulseFit::StFcsPulseFit() -{ - mName = "StFcsPulseFit"; - Initialize(); -} - -StFcsPulseFit::StFcsPulseFit(std::string name) -{ - mName = name; - mGAE_Signal = 0;//This is so 'Initialize' will create the TGraphAsymmErrors object - Initialize(); -} - -StFcsPulseFit::StFcsPulseFit(TGraphAsymmErrors* Sig, std::string name) -{ - mName = name; - mGAE_Signal = Sig;//Set graph object now so doesn't get created in 'Initialize' - Initialize();//This sets internal signal to true so set it false after call - mInternalSignal=false; -} - -/* -StFcsPulseFit::StFcsPulseFit(int ntb, double* tb, double* adc, std::string name) -{ - mName = name; - mGAE_Signal = new TGraphAsymmErrors(ntb, tb, adc); - Initialize(); -} -*/ - -void StFcsPulseFit::Initialize() -{ - mDEBUG = false; - mWindowSum = false; mSumAdc=0; - mFitSum = false; mSumAdcFit=0.0; - mSearch.Start = 50; mSearch.End=5;//Chosen to match "WaveFormFitMaker" where 50 is the peak - - mT0Computed = false;//This will keep from having to determine T0 every time - mInternalSignal = true; - mSignalT0.Reset(); - mBaseline = -5.0; - mBaselineSigma = 0.75; - - mMaxAdc = -5.0; - mMaxTb = -5.0; - - if( mGAE_Signal==0 ) - { - mGAE_Signal = new TGraphAsymmErrors(); - mGAE_Signal->SetNameTitle((mName+"_"+"G_Signal").c_str(),"Signal Pulse"); - } - mH1_Baseline = 0; - mF1_BaselineFit = 0; - - mF1_SignalFit = 0; - -} - -void StFcsPulseFit::SetSignal(TGraphAsymmErrors* SigGraph) -{ - if( SigGraph==0 ){std::cout << "Signal graph cannot be 0" << std::endl; exit(0);} - ResetFinder(); - if( mInternalSignal ){delete mGAE_Signal;} - mGAE_Signal = SigGraph; - mInternalSignal = false; -} - -StFcsPulseFit::~StFcsPulseFit() -{ - if(mGAE_Signal!=0 && mInternalSignal ){delete mGAE_Signal;} - if(mH1_Baseline!=0) { delete mH1_Baseline; } - if(mF1_SignalFit!=0) { delete mF1_SignalFit; } - if(mF1_BaselineFit!=0) { delete mF1_BaselineFit; } - -} - -void StFcsPulseFit::ResetFinder() -{ - if( mInternalSignal ){mGAE_Signal->Set(0);}//Effectively deletes all points in the graph only for internal graph objects - //delete mGAE_Signal; mGAE_Signal = new TGraphAsymmErrors(); - ResetBaseline(); - ResetT0(); - ResetSum(); -} - -void StFcsPulseFit::ResetBaseline() -{ - mBaseline = -5.0; - mBaselineSigma = 0.75; - if( mH1_Baseline!=0 ){mH1_Baseline->Reset();} - if( mF1_BaselineFit!=0 ){mF1_BaselineFit->ResetAttLine();}//Values from old fit will be replaced when 'AnalyzeForBaseline' gets called -} - -void StFcsPulseFit::ResetT0() -{ - mT0Computed = false; - mSignalT0.Reset(); -} - -void StFcsPulseFit::ResetSum() -{ - if( mF1_SignalFit!=0 ){ mF1_SignalFit->ResetAttLine(); }//This was for plotting but still good to do just in case - - mWindowSum=false; mSumAdc = 0; - mFitSum=false; mSumAdcFit = 0.0; -} - -void StFcsPulseFit::SetBaseline(Double_t value, Double_t sigma) -{ - if( value<=-0.0001 || sigma<=-0.0001){ std::cout << "WARNING:Baseline cannot be less than zero\nKeeping default behavior to search for baseline" << std::endl; return; } - mBaseline = value; - mBaselineSigma = sigma; - -} - -void StFcsPulseFit::AnalyzeForPedestal() -{ - ResetBaseline();//Resets histogram and function if they exist - if( mH1_Baseline==0){mH1_Baseline = new TH1F( (mName+"_"+"H1_Baseline").c_str(),"",4096,-0.5,4095.5); mH1_Baseline->Sumw2(); mH1_Baseline->SetTitle("Baseline Histogram");} - //Fill histogram of ADC values - for( int ipoint=0; ipointGetN(); ++ipoint) - { - Double_t tb; Double_t adc; - Int_t check = mGAE_Signal->GetPoint(ipoint,tb,adc); - if( check<0 ){ std::cout << "WARNING:Unable to find point " << ipoint << " in Signal graph" << std::endl; continue;} - if( adc<0.1 ){continue;} //skip adc=0 from histogram since this is most likely missing data in non-pedestal subtracted data - mH1_Baseline->Fill(adc); - } - //Create a function to fit to histogram of ADC values - if( mF1_BaselineFit==0 ){ mF1_BaselineFit = new TF1((mName+"_"+"F1_BaselineFit").c_str(),"gaus(0)", - mH1_Baseline->GetXaxis()->GetXmin(), - mH1_Baseline->GetXaxis()->GetXmax() - ); - } - if( FindBaseline() ){return;} - else - { - std::cout << "WARNING:Unable to find a proper baseline" << std::endl; - if(mDEBUG){std::cout << "->";PrintInfo();} - else{std::cout << "->Found baseline:"<< mBaseline << " Found BaselineSigma:"<Setting baseline to 0 and sigma to 0.75 to prevent T0 algorithm from failing" << std::endl; - mBaseline=0; - mBaselineSigma=0.75; - return; - } -} - -//In the future make range around XStartVal as arguments to the function -bool StFcsPulseFit::FindBaseline() -{ - //std::cout << "In FindBaseline:" << std::endl; - if( mBaseline > -0.001 ) - { - if( mDEBUG) - { - std::cout << "|BaselineValue:"<GetEntries()==0 ){mBaseline = -5.0; return false;} - //For nonpedestal data use full range for pedestal data find largest bin and use a +-5 bin buffer for the fit - Int_t XStartVal = mH1_Baseline->GetMaximumBin()-1;//Offset by one to make the value correct - if( mDEBUG ){ std::cout << "|Height:"<GetBinContent(XStartVal) << "|XStartVal:"<< XStartVal << "|Range:"<<5 <<"|StartSigma:"<<0.9 <SetRange(XStartVal-5,XStartVal+5);//Only do within this range - mF1_BaselineFit->SetParameter(0,mH1_Baseline->GetBinContent(XStartVal)); - mF1_BaselineFit->SetParameter(1,XStartVal);//Use maximum instead of mean for more accurate center - mF1_BaselineFit->SetParameter(2,0.9);//Right now fixed but may need changing in the future? - /*Change this to use set function - if( mPed==1 ){mF1_BaselineFit->SetParameter( 2,mH1_Baseline->GetRMS() );}//For pedestal RMS is good enough start point - else{mF1_BaselineFit->SetParameter(2,0.6);}//0.6 is a good start value RMS is too large for non-pedestal data - if( mPed==0 )//For physics fix the baseline to 0.0 with a sigma of 0.75 with tb range for 2019 in future change hard coded values - { - mF1_BaselineFit->SetRange(0,350); - mF1_BaselineFit->SetParameter(0,1); - mF1_BaselineFit->SetParameter(1,0.0); - mF1_BaselineFit->SetParameter(2,0.75);//0.75 range so T0 finder will look for ADC > 3 - } - else{mH1_Baseline->Fit(mF1_BaselineFit,"QNR");}//For anything other than physics data fit to find baseline - */ - //Base line is mean of the gaussian fit - Int_t FitStatus = mH1_Baseline->Fit(mF1_BaselineFit,"QNR"); - if( FitStatus >= 0 ) - { - Double_t ratio = mF1_BaselineFit->GetChisquare()/static_cast(mF1_BaselineFit->GetNDF()); - if( ratio > 10.0 ){ mBaseline = XStartVal; return true; }//Take care of bad baseline fits by setting baseline to maximum - else - { - mBaseline = static_cast(mF1_BaselineFit->GetParameter(1)); - mBaselineSigma = static_cast( fabs(mF1_BaselineFit->GetParameter(2)) );//fabs in case sigma is negative - return true; - } - } - else{return false;} -} - -//Vector will hold bin where possible T0 occurred -std::vector StFcsPulseFit::GetPossibleT0(Float_t SigmaScale) -{ - std::vector PossibleT0;//Gather all the possible occurances when signal is larger than Sigma - //std::cout << "In GetPossible T0" << std::endl; - if( !FindBaseline() ){ std::cout << "ERROR:No valid baseline\nPlease either run \"AnalyzeForPedestal\" or call \"SetBaseline\" with values greater than or equal to zero" << std::endl; return PossibleT0; } - Double_t baseline = mBaseline;//Already checked baseline above - Double_t SlopeCutoff = SigmaScale*fabs(mBaselineSigma); - //Ctools::PrintDebug(std::cout, "baseline", baseline ); - //Ctools::PrintDebug(std::cout, "SlopeCutoff", SlopeCutoff ); - //std::cout << std::endl; - SigWindow sig; - Int_t LocalMax=-5;//Variable to help keep track of when slope changes - //std::cout << "Finding T0 for cutoff " << SlopeCutoff << std::endl; - //The idea is that you start with all things negative and loop through all the ADC values. When the slope is greater than the cutoff you save it and then the signal will rise and then fall so slope will eventually go negative and this is when you keep track of the end values. Once signal goes positive again or you no longer pass the slope cutoff stop and save the start and stop values - if( mDEBUG ){std::cout << "StFcsPulseFit::GetPossibleT0:Start graph reading loop" << std::endl;} - for( Int_t ismp=0; ismpGetN()-1; ++ismp ) - { - Double_t Ltb; Double_t Ladc;//L for left (actually current point) - Double_t Rtb; Double_t Radc;//R for right (actually next point) - mGAE_Signal->GetPoint(ismp,Ltb,Ladc); - mGAE_Signal->GetPoint(ismp+1,Rtb,Radc); - //std::cout << " - |P:"< 4095){std::cout << "Error:Invalid Data 'Ladc'"< 4095){std::cout << "Error:Invalid Data 'Radc'"<0; however since I dynamically change the "end" time until the slope changes to positive this statement ensures that 'continue' only gets called on positive slope results and not negative slopes when it is trying to find the correct end time. - //if( (sig.Start>=0 && sig.End!=2000 && Ladc>0) ){continue;} - - //if( fabs(Slope) > SlopeCutoff*10.0 ){ continue; }//Avoid large sudden changes in slope (Hard code to 10 times the slope cutoff for now in future think of better method - //std::cout << "Val|ismp:"< SlopeCutoff || Ladc>baseline+SlopeCutoff ) - { - //std::cout << " + Passed Slope and baselineCutoff setting as start time" << std::endl; - sig.Start=Ltb;//Set start time - LocalMax=Ltb;//Start checking local maximums - } - //If didn't pass thresholds for start time then continue to next point - } - else//Found a suitable start time which indicates positve slope above some thershold - { - //The statements below will now check subsequent values and if the slope is >= 0 then we have wrong local max so set local max to next value. Note this means that localmax is not greater than End which is 2000. Once slope goes below 0 End will take that value and local max will be less than End. Once the slope changes again LocalMax will now be greater than End - //std::cout << " + Found StartTime:"<< sig.Start << std::endl; - if( Slope >= 0 ){ LocalMax=Ltb;/* std::cout<<" + Slope>=0"<sig.End || Ladc 0 ){ sig.End=ismp; } - //if( LBinValue>baseline+SlopeCutoff ){sig.End=ismp;} - //else{PossibleT0.push_back(sig); sig.Reset();} - } - - } - //PossibleT0.push_back(sig); - if( mDEBUG ){ std::cout << "StFcsPulseFit::GetPossibleT0:End graph reading loop|SizeT0:" < &PossibleT0s) -{ - if( PossibleT0s.size()==0 ) - { - if( mDEBUG ){std::cout << "Error:Unable to find a valid T0\nReturning 2000" << std::endl; } - return 2000; - } - /* - else if( PossibleT0s.size()==1 ) - { - if( mDEBUG ){std::cout << "Only found one T0" << std::endl;} - return 0; - } - */ - else - { - //Add this check even if only one entry to get rid of a false positive somewhere outside the range? - //If multiple T0s are found only pick the one with a specific start time - //Int_t StartTime=120;//Also good for pedestal which shouldn't find a T0 - //The values chosen are to match the trigger for 2019 for which the physics trigger fired around timebin ~40 and the pulser fires at timebin ~210 - /*Change this to use set function - if( mPed==0 ){ StartTime = 40; }//physics run - if( mPed==1 ){return 2000;} //For obvious reasons don't do pedestal runs - if( mPed==2 ){ StartTime = 210;}//pulser run - */ - if( mDEBUG ) - { - std::cout << "|Checking Name:"<GetPoint(PossibleT0s.at(it0).P_Peak,PeakLoc,temp); - if( mSearch.Start-mSearch.End<=PeakLoc && PeakLoc <= mSearch.Start+mSearch.End ) - { - T0index=it0; - if( mDEBUG ) - { - std::cout << " + "; - std::cout << "|TrueIndex:"<0 && PossibleT0s.at(T0index).Start>0 ) - if( T0index>=0 && mSearch.Start>0 && PossibleT0s.at(T0index).P_Peak>0 ) - { - if(mDEBUG){PossibleT0s.at(T0index).Print();std::cout< MyT0s = GetPossibleT0(Sigma);//This is the sigma to compute the threshold - UShort_t FoundT0 = SearchForT0(MyT0s);//Returns either a valid index for the possible T0 vector or 2000 - //1024 samples max so that should be the upper bound on the returned index - if( FoundT0 < 1024 ){mSignalT0 = MyT0s.at(FoundT0);}//No offset for now rethink in the future - else{mSignalT0.SetWindow(FoundT0,FoundT0);}//If no valid index then set SignalT0 to 2000 so it is obvious something's wrong - mT0Computed=true; - return; -} - -bool StFcsPulseFit::GoodWindow() -{ - if( !mT0Computed ){AnalyzeForT0();} - //First check if start and end times are within our max timebin window of 0-1023 - if( mSignalT0.Start<0 || mSignalT0.Start>1023 ){/*std::cout<<"StartOut"< 1023 ){/*std::cout<<"EndOut"< mSignalT0.End ){/*std::cout<<"Start>End"<GetN()==0 ){return SumAdc;} - //Ctools::PrintDebug(std::cout,"SumAdc",SumAdc); - for( Int_t ismp=0; ismpGetN(); ++ismp ) - { - Double_t tb; Double_t adc; - mGAE_Signal->GetPoint(ismp,tb,adc); - if( tb < Start ){continue;} - if( tb > End ){ break; } - //Ctools::PrintDebug(std::cout,"SumAdc",SumAdc); - //Ctools::PrintDebug(std::cout,"tb",tb); - //Ctools::PrintDebug(std::cout,"adc",adc); - //std::cout << std::endl; - SumAdc += adc; - } - if( mDEBUG ){std::cout << "|Start:"<GetN()==0 ){return mSumAdcFit;} - if( Start==0 && End==0 ){Start=SignalStart(); End=SignalEnd();}//Default case - //Upon failure to find a window 'Start' and 'End' value will be equal. In this case just set some default window. For run 2019 it will be 35-60 and triggered crossing was at 50 - //For 2019 Ecal Cosmic running signal window was between 110-140 - //if( Start==End ){ SetWindow(110,140); Start=115; End=135; }//40 and 55 here since I add to range later - //Sanity Checks - if( !GoodWindow() ){ return mSumAdcFit; } - //Always make range a little bigger - if( mF1_SignalFit==0 ){ mF1_SignalFit = new TF1( (mName+"_"+"mF1_SignalFit").c_str(),"gaus(0)+[3]",Start,End); } - else{ mF1_SignalFit->SetRange(Start,End); } - //mGAE_Signal->GetXaxis()->SetRangeUser(Start-5,End+5); - mF1_SignalFit->SetParameter( 0, MaxAdc() ); - mF1_SignalFit->SetParameter( 1, MaxTb() ); - mF1_SignalFit->SetParameter( 2, static_cast(End-Start)*0.2 );//Seemed to work most of the time with 0.3, 0.5 seemed too much - mF1_SignalFit->SetParameter( 3, base ); - Int_t fitstatus = mGAE_Signal->Fit( mF1_SignalFit, "QNR" ); - if( fitstatus >= 0 ) - { - //Fix to 4 sigma for now in the future make it a variable? - Double_t FitStart = mF1_SignalFit->GetParameter(1) - 4.0*fabs(mF1_SignalFit->GetParameter(2)); - Double_t FitEnd = mF1_SignalFit->GetParameter(1) + 4.0*fabs(mF1_SignalFit->GetParameter(2)); - mSumAdcFit = mF1_SignalFit->Integral(FitStart,FitEnd);//Sum only in the range and subtract the area under the baseline - if( mDEBUG ){std::cout << "|FitStart:"<GetN(); ++i ) - { - Double_t X; Double_t Y; - mGAE_Signal->GetPoint(i,X,Y); - if( Xxmax ){continue;} - if( Y>mMaxAdc ){ mMaxAdc=Y; mMaxTb=X; } - } - if( mMaxAdc<0 && mMaxTb<0 ) - { - std::cout << "Unable to find a maximum adc\nSetting to impossible values" << std::endl; - mMaxAdc = 5000; - mMaxTb = 5000; - } -} - -Double_t StFcsPulseFit::SignalPeakTb() -{ - if( !mT0Computed ){ AnalyzeForT0(); } - if( mSignalT0.P_Peak<0){ return 0; } - else{ Double_t x,y; mGAE_Signal->GetPoint(mSignalT0.P_Peak,x,y); return x;} -} - -Double_t StFcsPulseFit::SignalPeakAdc() -{ - if( !mT0Computed ){AnalyzeForT0();} - if( mSignalT0.P_Peak<0){ return 0; } - else{ Double_t x,y; mGAE_Signal->GetPoint(mSignalT0.P_Peak,x,y); return y;} -} - -void StFcsPulseFit::SignalPeak(Double_t &tb, Double_t &adc) -{ - if( !mT0Computed ){AnalyzeForT0();} - if( mSignalT0.P_Peak<0 ){ tb=0; adc=0; } - else{ mGAE_Signal->GetPoint(mSignalT0.P_Peak,tb,adc); } -} - -void StFcsPulseFit::PrintInfo() const -{ - std::cout << "|Name:"<=0. The sigma set will determine the threshold ADC. ThresholdADC=sigma*4. The 4 is default but can be changed - See "AnalyzeForT0". This means default ADC threshold is 3 for zs-pedestal subtracted data - //void SetBaselineRange(Int_t Nbins, Double_t Xmin, Double_t Xmax);//Need to implment? - void SetWindow(const Int_t &start, const Int_t &end ){mSignalT0.SetWindow(start,end);mT0Computed=true;}//Set signal window to user defined values. - void SetSearchWindow(Int_t peak, Int_t width){ mSearch.SetWindow(peak,width); }//Set start time to search for and the width of the start time - void SetZS(){SetBaseline(0,0.6);}//Call this for ZS data which uses thresholds that are relevant for that. (like 0 baseline and 0.5 sigma so thereshold is at 2 since ZS is pedestal subtracted. Maybe even 0.3 so it is above one. - - void AnalyzeForPedestal(); //Analyze graph data to determine baseline internally - - void AnalyzeForT0(Float_t Sigma=4.0);//Signal cutoff is 4(sigma)*BaselineSigma above baseline - Int_t SignalStart(){if( !mT0Computed ){AnalyzeForT0();} return mSignalT0.Start;}//Found Signal Start time - Int_t SignalEnd(){ if( !mT0Computed ){AnalyzeForT0();} return mSignalT0.End;}//Found Signal End time - Double_t SignalPeakTb();//Timebin of found signal peak - Double_t SignalPeakAdc();//ADC of found signal peak - void SignalPeak(Double_t &tb, Double_t &adc); - - Double_t Baseline(){ if( FindBaseline() ){return mBaseline;}else{return -5.0;} }//If baseline is found will only return the value - Double_t BaselineSigma(){ if( FindBaseline() ){return fabs(mBaselineSigma);}else{return 0.75;} }//fabs in case the found sigma is negative - //Need checks for fitting future? - //Double_t SignalFitMean(){return mF1_SignalFit->GetParameter(1);} - //Double_t SignalFitSigma(){return fabs(mF1_SignalFit->GetParameter(2));} - - void GetXYMax(Double_t xmin=-5, Double_t xmax=2000); - Double_t MaxAdc(){if( mMaxAdc<0 ){GetXYMax();} return mMaxAdc; } - Double_t MaxTb(){if( mMaxTb<0 ){GetXYMax(); } return mMaxTb; } - - bool GoodWindow(); - - TH1F* BaselineHist(){return mH1_Baseline;} - TF1* SignalFit(const char option ='n');//Replace with a function that actually calls "fit" and then returns the function?Arguments can be values and range?A little more complicated since hard to tell when user would want to call a fit or just get the function itself?Maybe no arguments means give me function, arguments means do a refit? - TF1* BaselineFit(){return mF1_BaselineFit;}//see above "^"? - - //void ReFitBaseline(){mG_Signal->Fit(mF1_BaselineFit,"R");} - - virtual void PrintInfo() const; - - protected: - - bool mDEBUG; - virtual void Initialize(); //Sets everything to zero except signal and background histograms - - bool mT0Computed; //This will keep from having to determine T0 every time - //bool mFindBase; //This will be used to know whether baseline needs to be found or not - - //This struct is for holding the signal window - struct SigWindow - { - Int_t Start; Int_t End; - Int_t P_Peak;//Point Number of peak (P for point) - //Default values are two extremes (Start>End otherwise algorithm won't work properly) - SigWindow(){Start=-5;End=2000;P_Peak=-5;}//2000 since largest timebin is 1024 - void SetWindow(Int_t s, Int_t e){Start=s; End=e;} - void GetWindow(Int_t &s, Int_t &e){s=Start;e=End;} - void Reset(){Start=-5.0; End=2000;P_Peak=-5.0;} - void Print()const{std::cout << "|Start:"< GetPossibleT0(Float_t SigmaScale);//Finds all possible T0s in signal with some criteria - UShort_t SearchForT0(const std::vector &PossibleT0s);//Finds signal start time by checking which start time fits inside the trigger window - - private: - - std::string mName; - - bool mInternalSignal; - bool mWindowSum; - bool mFitSum; - - //Variables to hold sum values. This is to avoid tedious recomputation - Int_t mSumAdc; - Double_t mSumAdcFit; - - Double_t mMaxAdc;//Maximum value for Adc - Double_t mMaxTb;//Timebin where maximum value of Adc occurs - - - ClassDef(StFcsPulseFit,1) - -}; - -#endif diff --git a/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.cxx b/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.cxx index de3edc69fec..e7dadb48226 100644 --- a/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.cxx +++ b/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.cxx @@ -71,6 +71,7 @@ ClassImp(StFcsWaveformFitMaker) #include "StEventTypes.h" #include "StEvent/StFcsHit.h" #include "StFcsDbMaker/StFcsDb.h" +#include "StFcsDbMaker/StFcsDbPulse.h" #include #include @@ -86,64 +87,253 @@ ClassImp(StFcsWaveformFitMaker) #include "TStyle.h" #include "TROOT.h" -namespace { - static const double sqrtpi = sqrt(3.141592654); - static const double sqrt2pi = sqrt(2.0 * 3.141592654); - static const double TBPerRC = 8; - static const double nsecPerTB = 107.0/TBPerRC; - - //pulse shape - double GSigma=24.5/nsecPerTB; - double A1=0; - double A2=0; - double Xoff1=0; - double Xoff2=0; - double Tau1=1; - double Tau2=1; - double P1=1; - double P2=1; - - //Data from Gerard 2020 summer - static const double GSigma_1 = 24.5/nsecPerTB; - static const double A1_1 = 1.0/0.154/GSigma_1; - static const double A2_1 = 0.2/0.154/GSigma_1; - static const double Xoff1_1 = (70 - 129)/nsecPerTB; - static const double Xoff2_1 = (220 - 129)/nsecPerTB; - static const double Tau1_1 = 200.0/nsecPerTB; - static const double Tau2_1 = 40.0/nsecPerTB; - static const double P1_1 = 1.0; - static const double P2_1 = 1.0; - - //Data from WAH with real detector/LED system 2021 Jan - static const double GSigma_2 = 2.347; - static const double A1_2 = 2543.0/854.0/GSigma_2; - static const double A2_2 = 0.0; - static const double Xoff1_2 = 211.3-215.7; - static const double Xoff2_2 = 0.0; - static const double Tau1_2 = 4.375; - static const double Tau2_2 = 0.0; - static const double P1_2 = 1.0; - static const double P2_2 = 0.0; -} StFcsWaveformFitMaker::StFcsWaveformFitMaker(const char* name) : StMaker(name) { mChWaveData.SetClass("TGraphAsymmErrors"); //Initialize with only one graph at first - mPulseFit = new StFcsPulseFit( (TGraphAsymmErrors*)mChWaveData.ConstructedAt(0) ); - mEnergySelect[0]=10; //default gaus fit for Ecal - mEnergySelect[1]=10; //default gaus fit for Hcal + mPulseFit = 0; + + mOutFile = 0; + + mEnergySelect[0]=13; //default PulseFit2 for Ecal + mEnergySelect[1]=13; //default PulseFit2 for Hcal mEnergySelect[2]=1; //default sum8 for Pres + + for( UShort_t i=0; i<7; ++i ){ + if( i<3 ){ + mH2_Dep0DepMod[i]=0; + mH2_Sum8Dep0[i]=0; + mH2_Sum8DepMod[i]=0; + } + if( i<6 ){ + mH2F_AdcTbAkio[i] = 0; + mH2F_AdcTbMine[i] = 0; + mH2F_SumFitvSumWin[i] = 0; + mH2F_APeakvMPeak[i] = 0; + mH1F_PeakStart[i] = 0; + mH1F_PeakEnd[i] = 0; + mH2F_NOvsId[i] = 0; + } + mH2F_AdcTbValidPeak[i] = 0; + mH1F_NPeaks[i] = 0; + mH1F_NPeaksFiltered[i] = 0; + mH1F_Res0[i] = 0; + mH1F_Res0Zoom[i] = 0; + mH1F_Sum8Res0[i] = 0; + mH1F_Sum8Res0Zoom[i] = 0; + mH1F_FitRes0[i] = 0; + mH1F_FitRes0Zoom[i] = 0; + mH2F_Sum8vFit[i] = 0; + } +} + +StFcsWaveformFitMaker::~StFcsWaveformFitMaker() { + mChWaveData.Delete(); + delete mPulseFit; + for( UShort_t i=0; i<7; ++i ){ + if( i<3 ){ + delete mH2_Dep0DepMod[i]; + delete mH2_Sum8Dep0[i]; + delete mH2_Sum8DepMod[i]; + } + if( i<6 ){ + delete mH2F_AdcTbAkio[i]; + delete mH2F_AdcTbMine[i]; + delete mH2F_SumFitvSumWin[i]; + delete mH2F_APeakvMPeak[i]; + delete mH1F_PeakStart[i]; + delete mH1F_PeakEnd[i]; + delete mH2F_NOvsId[i]; + } + delete mH2F_AdcTbValidPeak[i]; + delete mH1F_NPeaks[i]; + delete mH1F_NPeaksFiltered[i]; + + delete mH1F_Res0[i]; + delete mH1F_Res0Zoom[i]; + delete mH1F_Sum8Res0[i]; + delete mH1F_Sum8Res0Zoom[i]; + delete mH1F_FitRes0[i]; + delete mH1F_FitRes0Zoom[i]; + delete mH2F_Sum8vFit[i]; + } + delete mTime; + + delete mH1_NPeaksAkio; + delete mH1_NPeaksFilteredAkio; + + delete mH1_PeakTiming; + + delete mH2_NPeakvsPeakXdiff; + delete mH2_NPeakvsPeakYratio; + delete mH1_VOverlap; + delete mH2_NOvsNPeaks; + delete mH2_VvsNOverlap; + + delete mH1_TimeFitPulse; + + delete mH2_HeightvsSigma; + delete mH2_HeightvsSigmaTrig; + delete mH1_ChiNdf; + delete mH2_HeightvsChiNdf; + delete mH2_MeanvsChiNdf; + delete mH2_SigmavsChiNdf; + + delete mH1_PeakTimingGaus; + delete mH1_PeakTimingPuls; + delete mH2_PeakTimingCompare; + + if( mOutFile!=0 ){ mOutFile->Close(); delete mOutFile; } } -StFcsWaveformFitMaker::~StFcsWaveformFitMaker() { - mChWaveData.Delete();//Completely clear all graphs (Put in Finish??) - delete mPulseFit; +void StFcsWaveformFitMaker::writeFile(std::string filename){ + if( filename.size()==0 ){return;} + if( mOutFile!=0 ){ mOutFile->Close(); delete mOutFile; mOutFile=0; } + mOutFile = new TFile(filename.c_str(),"RECREATE"); } +void StFcsWaveformFitMaker::setTest(int v) +{ + mTest=v; + if( mOutFile==0 ){ writeFile("test.root"); }//All test levels need to have a file output to write to + if( mTest==1 ){//Test level 1 is for testing the DEP triggering algorithm to currently selected energy method. If method is zero switch on most basic which is sum 8 + if( mEnergySelect[0]==0 ){mEnergySelect[0]=1; mEnergySelect[1]=1; mEnergySelect[2]=1; } + } + if( mTest==2 ){//Test level 2 is for testing *StFcsPulseAna* to gausFit() which requires either energy select 10 or 11. + if( mEnergySelect[0]!=10 || mEnergySelect[0]!=11 ){ mEnergySelect[0]=10;} + if( mEnergySelect[1]!=10 || mEnergySelect[1]!=11 ){ mEnergySelect[1]=10;} + if( mEnergySelect[2]!=10 || mEnergySelect[2]!=11 ){ mEnergySelect[2]=10;} + } + if( mTest==3 ){//Test level 3 is for testing the steering code in PulseFit1() which is energy select 12 + if( mEnergySelect[0]!=12 ){ mEnergySelect[0]=12; } + if( mEnergySelect[1]!=12 ){ mEnergySelect[1]=12; } + if( mEnergySelect[2]!=12 ){ mEnergySelect[2]=12; } + } + if( mTest==4 ){//Test level 4 is for testing fitting all signals using my pulsefit code + if( mEnergySelect[0]!=12 ){ mEnergySelect[0]=12; } + if( mEnergySelect[1]!=12 ){ mEnergySelect[1]=12; } + if( mEnergySelect[2]!=12 ){ mEnergySelect[2]=12; } + } + if( mTest==5 ){//Test level 5 is for testing timing on gausFit() and PulseFit1() so all is needed is that the energy select is not zero and drawing is off otherwise it will double up the hits in mChWaveData + if( mEnergySelect[0]==0 ){ mEnergySelect[0]=1; mEnergySelect[1]=1; mEnergySelect[2]=1; mFitDrawOn=0; } + } + if( mTest==6 ){//Test level 6 is for testing the steering code in PulseFit2() + if( mEnergySelect[0]!=13 ){ mEnergySelect[0]=13; } + if( mEnergySelect[1]!=13 ){ mEnergySelect[1]=13; } + if( mEnergySelect[2]!=13 ){ mEnergySelect[2]=13; } + } + if( mTest==7 ){//Test level 6 is for testing the steering code in PulseFit2() + if( mEnergySelect[0]!=13 ){ mEnergySelect[0]=13; } + if( mEnergySelect[1]!=13 ){ mEnergySelect[1]=13; } + if( mEnergySelect[2]!=1 ){ mEnergySelect[2]=1; } + } +} + +void StFcsWaveformFitMaker::setTail(int v){mTail=v;} + void StFcsWaveformFitMaker::Clear(Option_t* option) { - if(mFitDrawOn==1) mHitIdx=0; + if(mFitDrawOn>=1) mHitIdx=0; StMaker::Clear(option); } +int StFcsWaveformFitMaker::Init() +{ + if(mFilename){ + mCanvas=new TCanvas("FCSWaveFormFit","FCSWaveFormFit",10,10,2000,2000); + gStyle->SetOptStat(0); + char file[100]; + if( mFitDrawOn>=0 ){ + sprintf(file,"%s.pdf[",mFilename);//Opens fileA but does not save anything + mCanvas->Print(file); + } + } + if(mMeasureTime){ + mTime=new TH1F("FitTime","FitTime;msec",1000,0,1000); + } + if(mFilter && mFilter[0]=='0'){ + mTimeIntg[0]=new TH2F("Noboth", "No both; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); + mTimeIntg[1]=new TH2F("NoFall", "No Fall; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); + mTimeIntg[2]=new TH2F("NoRise", "No Rise; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); + mTimeIntg[3]=new TH2F("Accepted","Accepted; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); + } + if( mOutFile!=0 ){ + for( UShort_t i=0; i<7; ++i ){ + std::stringstream ss; + ss << i; + if( i<3 ){ + if( mTest==1 ){ + mH2_Dep0DepMod[i]=new TH2F( ("H2_Dep0DepMod_"+ss.str()).c_str(), ("Sum Dep0 vs DepMod "+ss.str()).c_str(),100,0,1000,100,0,1000); + mH2_Sum8Dep0[i]=new TH2F( ("H2_Sum8Dep0_"+ss.str()).c_str(), ("Sum 8 vs Dep0 "+ss.str()).c_str(),100,0,1000,100,0,1000); + mH2_Sum8DepMod[i]=new TH2F( ("H2_Sum8DepMod_"+ss.str()).c_str(), ("Sum 8 vs DepMod "+ss.str()).c_str(),100,0,1000,100,0,1000); + } + } + if( i<6 ){ + if( mTest==2 ){ + mH2F_AdcTbAkio[i] = new TH2F( ("H2_AdcTbAkio_"+ss.str()).c_str(),"Akio Peaks Adc vs. Tb",102,-1.5,100.5,4097,-0.5,4096.5); + mH2F_AdcTbMine[i] = new TH2F( ("H2_AdcTbMine_"+ss.str()).c_str(),"My Peaks Adc vs. Tb",102,-1.5,100.5,4097,-0.5,4096.5); + mH2F_SumFitvSumWin[i] = new TH2F( ("H2_SumFitvSumWin_"+ss.str()).c_str(),"Fitted Sum vs. Peak Window Sum",100,0,2000,100,0,2000); + mH2F_APeakvMPeak[i] = new TH2F( ("H2_APeakvMPeak_"+ss.str()).c_str(),"Peak Locations Akio v. Mine",80,-0.5,79.5,80,-0.5,79.5); + mH1F_PeakStart[i] = new TH1F( ("H1_PeakStart_"+ss.str()).c_str(),"PeakWindow start times",102,-1.5,100.5); + mH1F_PeakEnd[i] = new TH1F( ("H1_PeakEnd_"+ss.str()).c_str(),"PeakWindow end time",102,-1.5,100.5); + } + if( mTest==3 ){ + mH2F_NOvsId[i] = new TH2F( ("H2_NOvsId_"+ss.str()).c_str(),"Number of overlaps vs. Ch Id", 748,-0.5,747.5, 11,-0.5,10.5); + } + } + //i<7 histograms + if( mTest==2 ){ + if( i<1 ){ //Test=2 only needs 1 histogram + mH1F_NPeaks[i] = new TH1F( ("H1_NPeaks_"+ss.str()).c_str(),"Number of peaks from finder", 11,-0.5,10.5); + mH1F_NPeaksFiltered[i] = new TH1F( ("H1_NPeaksFiltered_"+ss.str()).c_str(),"Number of peaks from finder when a valid peak was found", 11,-0.5,10.5); + } + mH2F_AdcTbValidPeak[i] = new TH2F( ("H2_AdcTbValidPeak_"+ss.str()).c_str(),"Valid Peaks Adc vs. Tb",102,-1.5,100.5,4097,-0.5,4096.5); + } + if( mTest==3 || mTest==6 || mTest==7 ){ + mH1F_NPeaks[i] = new TH1F( ("H1_NPeaks_"+ss.str()).c_str(),"Number of peaks from finder", 11,-0.5,10.5); + mH1F_NPeaksFiltered[i] = new TH1F( ("H1_NPeaksFiltered_"+ss.str()).c_str(),"Number of peaks from finder when a valid peak was found", 11,-0.5,10.5); + mH1F_Res0[i] = new TH1F( ("H1_Res0_"+ss.str()).c_str(),"All ADC sums", 100,0,2000 ); + mH1F_Res0Zoom[i] = new TH1F( ("H1_Res0Zoom_"+ss.str()).c_str(),"All ADC sums", 201,-0.5,200.5 ); + } + if( mTest==3 || mTest==6 ){ + mH1F_Sum8Res0[i] = new TH1F( ("H1_Sum8Res0_"+ss.str()).c_str(),"All ADC sums using sum 8", 100,0,2000 ); + mH1F_Sum8Res0Zoom[i] = new TH1F( ("H1_Sum8Res0Zoom_"+ss.str()).c_str(),"All ADC sums using sum 8", 201,-0.5,200.5 ); + mH1F_FitRes0[i] = new TH1F( ("H1_FitRes0_"+ss.str()).c_str(),"All ADC sums using fitting", 100,0,2000 ); + mH1F_FitRes0Zoom[i] = new TH1F( ("H1_FitRes0Zoom_"+ss.str()).c_str(),"All ADC sums using fitting", 201,-0.5,200.5 ); + mH2F_Sum8vFit[i] = new TH2F(("H2_Sum8vFit_"+ss.str()).c_str(),"sum 8 vs fit sum", 100,0,2000, 100,0,2000); + } + } + + if( mTest==2 ){ + mH1_NPeaksAkio = new TH1F("H1_NPeaksAkio","Number of peaks from current method", 11,-0.5,10.5); + mH1_NPeaksFilteredAkio = new TH1F("H1_NPeaksFilteredAkio","Number of peaks from current method when a valid peak is found", 11,-0.5,10.5); + mH1_PeakTiming = new TH1F("H1_PeakTiming","Timing to just find peak",200,0,5); + } + if( mTest==3 ){ + mH2_NPeakvsPeakXdiff = new TH2F("H2_NPeakvsPeakXdiff","NPeak vs. Peak X difference", 41,-20.5,20.5, 11,-0.5,10.5); + mH2_NPeakvsPeakYratio = new TH2F("H2_NPeakvsPeakYratio","NPeak vs. Peak Y ratio", 100,0,10, 11,-0.5,10.5); + mH1_VOverlap = new TH1F("H1_VOverlap","Peak Compare values",4,-0.5,3.5 ); + mH2_NOvsNPeaks = new TH2F("H2_NOvsNPeaks","Number of Overlapping peaks vs. Number of Peaks", 11,-0.5,10.5, 11,-0.5,10.5); + mH2_VvsNOverlap = new TH2F("H2_VvsNOverlap","Value of Comparison vs. Peak index", 21,-10.5,10.5, 4,-0.5,3.5); + } + if( mTest==3 || mTest==6 ){ + mH1_TimeFitPulse = new TH1F("H1_TimeFitPulse","FitTime;msec",1000,0,1000); + } + if( mTest==4 ){ + mH2_HeightvsSigma = new TH2F("H2_HeightvsSigma","Fitted Peak Height vs. Sigma;Sigma;Height", 31,-0.5,30.5, 4000,-0.5,3999.5); + mH2_HeightvsSigmaTrig = new TH2F("H2_HeightvsSigmaTrig","Fitted Peak Height vs. Sigma Triggered Xing;Sigma;Height", 31,-0.5,30.5, 4000,-0.5,3999.5); + mH1_ChiNdf = new TH1F("H1_ChiNdf","Chi^2/NDF for all fits;Chi^2/NDF", 50,-0.5,99.5); + mH2_HeightvsChiNdf = new TH2F("H2_HeightvsChiNdf", "Peak Height TXing vs. Chi^2/NDF;Chi^2/NDF;Height", 50,-0.5,99.5, 500,-0.5,499.5); + mH2_MeanvsChiNdf = new TH2F("H2_MeanvsChiNdf", "Peak Mean TXing vs. Chi^2/NDF;Chi^2/NDF;Mean", 50,-0.5,99.5, 21,39.5,60.5 ); + mH2_SigmavsChiNdf = new TH2F("H2_SigmavsChiNdf", "Peak Sigma TXing vs. Chi^2/NDF;Chi^2/NDF;Sigma", 50,-0.5,99.5, 21,-0.5,20.5 ); + } + if( mTest==5 ){ + mH1_PeakTimingGaus = new TH1F("H1_PeakTimingGaus","gausFit timing;msec", 1000,0,1000 ); + mH1_PeakTimingPuls = new TH1F("H1_PeakTimingPuls","PulseFit timing;msec", 1000,0,1000 ); + mH2_PeakTimingCompare = new TH2F("H2_PeakTimingCompare","PulseFit vs. gausFit;ms;ms", 200,0,6, 200,0,6 ); + } + } + return StMaker::Init(); +} + int StFcsWaveformFitMaker::InitRun(int runNumber) { LOG_DEBUG << "StFcsWaveformFitMaker initializing run" << endm; mDb = static_cast(GetDataSet("fcsDb")); @@ -151,51 +341,23 @@ int StFcsWaveformFitMaker::InitRun(int runNumber) { LOG_ERROR << "StFcsWaveformFitMaker initializing failed due to no StFcsDb" << endm; return kStErr; } - - if(mFilename){ - mCanvas=new TCanvas("FCSWaveFormFit","FCSWaveFormFit",10,10,2000,2000); - gStyle->SetOptStat(0); - } - if(mMeasureTime){ - mTime=new TH1F("FitTime","FitTime; msec",200,0,6); - } - if(mFilter && mFilter[0]=='0'){ - mTimeIntg[0]=new TH2F("Noboth", "No both; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); - mTimeIntg[1]=new TH2F("NoFall", "No Fall; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); - mTimeIntg[2]=new TH2F("NoRise", "No Rise; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); - mTimeIntg[3]=new TH2F("Accepted","Accepted; PeakTB; Integral",100,47.0,54.0,400,0.0,2000); + mDbPulse = static_cast(GetDataSet("fcsPulse")); + if (!mDbPulse) { + LOG_ERROR << "StFcsWaveformFitMaker initializing failed due to no StFcsDbPulse" << endm; + return kStErr; } + mDbPulse->setTail(mTail); + if( mPulseFit==0 ){ mPulseFit = new StFcsPulseAna(); SetupDavidFitterMay2022(); mPulseFit->setDbPulse(mDbPulse); } + mPulseFit->setDbPulse(mDbPulse); - if(mTail==1){ - GSigma = GSigma_1; - A1 = A1_1; - A2 = A2_1; - Xoff1 = Xoff1_1; - Xoff2 = Xoff2_1; - Tau1 = Tau1_1; - Tau2 = Tau2_1; - P1 = P1_1; - P2 = P2_1; - }else if(mTail==2){ - GSigma = GSigma_2; - A1 = A1_2; - A2 = A2_2; - Xoff1 = Xoff1_2; - Xoff2 = Xoff2_2; - Tau1 = Tau1_2; - Tau2 = Tau2_2; - P1 = P1_2; - P2 = P2_2; - } - return StMaker::InitRun(runNumber); } int StFcsWaveformFitMaker::Finish(){ if(mFilename && mPad>=0){ char file[200]; - sprintf(file,"%s.pdf)",mFilename); - mCanvas->Print(file,"pdf"); + sprintf(file,"%s.pdf]",mFilename); + mCanvas->Print(file); } if(mMeasureTime){ printf("WFFit Time Mean=%f RMS=%f\n",mTime->GetMean(),mTime->GetRMS()); @@ -205,6 +367,7 @@ int StFcsWaveformFitMaker::Finish(){ c->cd(1)->SetLogy(); mTime->Draw(); c->SaveAs(mMeasureTime); + delete c; } if(mFilter && mFilter[0]=='0'){ TCanvas *c= new TCanvas("Stage0","Stage0",10,10,800,800); @@ -216,6 +379,67 @@ int StFcsWaveformFitMaker::Finish(){ c->cd(3)->SetLogy(); mTimeIntg[2]->Draw("colz"); c->cd(4)->SetLogy(); mTimeIntg[3]->Draw("colz"); c->SaveAs("stage0.png"); + delete c; + } + if( mOutFile!=0 ){ + mOutFile->cd(); + + for( UShort_t i=0; i<7; ++i ){ + if( i<3 ){ + if( mH2_Dep0DepMod[i]!=0 ){mH2_Dep0DepMod[i]->Write();} + if( mH2_Sum8Dep0[i]!=0 ) {mH2_Sum8Dep0[i]->Write();} + if( mH2_Sum8DepMod[i]!=0 ){mH2_Sum8DepMod[i]->Write();} + } + if( i<6 ){ + if( mH2F_AdcTbAkio[i]!=0 ) { mH2F_AdcTbAkio[i]->Write();} + if( mH2F_AdcTbMine[i]!=0 ) { mH2F_AdcTbMine[i]->Write();} + if( mH2F_SumFitvSumWin[i]!=0 ) { mH2F_SumFitvSumWin[i]->Write();} + if( mH2F_APeakvMPeak[i]!=0 ) { mH2F_APeakvMPeak[i]->Write();} + if( mH1F_PeakStart[i]!=0 ) { mH1F_PeakStart[i]->Write();} + if( mH1F_PeakEnd[i]!=0 ) { mH1F_PeakEnd[i]->Write();} + if( mH2F_NOvsId[i]!=0 ) { mH2F_NOvsId[i]->Write();} + } + if( mH2F_AdcTbValidPeak[i]!=0 ){ mH2F_AdcTbValidPeak[i]->Write();} + if( mH1F_NPeaks[i]!=0 ) { mH1F_NPeaks[i]->Write();} + if( mH1F_NPeaksFiltered[i]!=0 ){ mH1F_NPeaksFiltered[i]->Write();} + + if( mH1F_Res0[i]!=0 ) { mH1F_Res0[i]->Write();} + if( mH1F_Res0Zoom[i]!=0 ) { mH1F_Res0Zoom[i]->Write();} + if( mH1F_Sum8Res0[i]!=0 ) { mH1F_Sum8Res0[i]->Write();} + if( mH1F_Sum8Res0Zoom[i]!=0 ) { mH1F_Sum8Res0Zoom[i]->Write();} + if( mH1F_FitRes0[i]!=0 ) { mH1F_FitRes0[i]->Write();} + if( mH1F_FitRes0Zoom[i]!=0 ) { mH1F_FitRes0Zoom[i]->Write();} + if( mH2F_Sum8vFit[i]!=0 ) { mH2F_Sum8vFit[i]->Write();} + } + + if( mTime!=0 ){ mTime->Write(); } + + if( mH1_NPeaksAkio!=0 ){ mH1_NPeaksAkio->Write(); } + if( mH1_NPeaksFilteredAkio!=0 ){ mH1_NPeaksFilteredAkio->Write(); } + if( mH1_PeakTiming!=0 ){mH1_PeakTiming->Write();} + + if( mH2_NPeakvsPeakXdiff!=0 ) {mH2_NPeakvsPeakXdiff->Write();} + if( mH2_NPeakvsPeakYratio!=0) {mH2_NPeakvsPeakYratio->Write();} + if( mH1_VOverlap!=0 ) {mH1_VOverlap->Write();} + if( mH2_NOvsNPeaks!=0 ) {mH2_NOvsNPeaks->Write();} + if( mH2_VvsNOverlap!=0 ) {mH2_VvsNOverlap->Write();} + if( mH1_TimeFitPulse!=0 ) {mH1_TimeFitPulse->Write();} + + if( mH2_HeightvsSigma!=0 ) {mH2_HeightvsSigma->Write();} + if( mH2_HeightvsSigmaTrig!=0) {mH2_HeightvsSigmaTrig->Write();} + if( mH1_ChiNdf!=0 ) {mH1_ChiNdf->Write();} + if( mH2_HeightvsChiNdf!=0 ) {mH2_HeightvsChiNdf->Write();} + if( mH2_MeanvsChiNdf!=0 ) {mH2_MeanvsChiNdf->Write();} + if( mH2_SigmavsChiNdf!=0 ) {mH2_SigmavsChiNdf->Write();} + + if( mH1_PeakTimingGaus!=0 ) {mH1_PeakTimingGaus->Write();} + if( mH1_PeakTimingPuls!=0 ) {mH1_PeakTimingPuls->Write();} + if( mH2_PeakTimingCompare!=0) {mH2_PeakTimingCompare->Write();} + + //In case you want to save a particular graph needs mFitDrawOn==1 + //TGraphAsymmErrors* g = getGraph(0,161); + //std::cout << "found graph:"<Write("GAE_det0id161"); } } return kStOK; } @@ -285,7 +509,20 @@ int StFcsWaveformFitMaker::Make() { //printf("Fit Time = %lld usec\n",usec); mTime->Fill(float(usec)/1000.0); } - + if( mTest==5 ){ + auto startg = std::chrono::high_resolution_clock::now(); + integral = analyzeWaveform(10,hits[i],res,func,ped); + auto stopg = std::chrono::high_resolution_clock::now(); + long long usecg = chrono::duration_cast(stopg-startg).count(); + auto startp = std::chrono::high_resolution_clock::now(); + integral = analyzeWaveform(12,hits[i],res,func,ped); + auto stopp = std::chrono::high_resolution_clock::now(); + long long usecp = chrono::duration_cast(stopp-startp).count(); + mH1_PeakTimingGaus->Fill(float(usecg)/1000.0); + mH1_PeakTimingPuls->Fill(float(usecp)/1000.0); + mH2_PeakTimingCompare->Fill(float(usecg)/1000.0,float(usecp)/1000.0); + if( usecp>200000.0 ){ std::cout << getGraph()->GetName() << std::endl; } + } } } return kStOk; @@ -307,18 +544,39 @@ TGraphAsymmErrors* StFcsWaveformFitMaker::getGraph(int idx){ return gae; } +TGraphAsymmErrors* StFcsWaveformFitMaker::getGraph(int det, int id) +{ + //Since mHitIdx may have been reset by call to clear use size of array instead this way you can find graphs even if Clear() has been called + for( Int_t i=0; iGetName(),det0,id0); + //std::cout << "|det0:"<=0 && det0==det ){ + if( id0>=0 && id0==id ){ + return gae; + } + } + } + LOG_ERROR << "Unable to find det:"<SetPointError(i,0,0,mError,HighY); } +*/ TGraphAsymmErrors* StFcsWaveformFitMaker::makeTGraphAsymmErrors(int n, double* tb, double* adc){ TGraphAsymmErrors* gae = resetGraph(); for(int i=0; iSetPoint(i,tb[i],adc[i]); - setTGraphAsymmErrors(gae, i, adc[i]); + StFcsDbPulse::setTGraphAsymmErrors(gae, i, adc[i],mError,mErrorSaturated); } return gae; } @@ -336,7 +594,7 @@ TGraphAsymmErrors* StFcsWaveformFitMaker::makeTGraphAsymmErrors(TH1* hist){ for(int i=1; i<=NValues; ++i ){ double adc = hist->GetBinContent(i); gae->SetPoint(i-1,hist->GetBinCenter(i),adc); - setTGraphAsymmErrors(gae,i-1,adc); + StFcsDbPulse::setTGraphAsymmErrors(gae,i-1,adc,mError,mErrorSaturated); } return gae; } @@ -353,7 +611,7 @@ TGraphAsymmErrors* StFcsWaveformFitMaker::makeTGraphAsymmErrors(StFcsHit* hit){ if(tb>=mMinTB){ int adc=hit->adc(i); gae->SetPoint(j,tb,adc); - setTGraphAsymmErrors(gae,j,adc); + StFcsDbPulse::setTGraphAsymmErrors(gae,j,adc,mError,mErrorSaturated); j++; } if(tb>=mMaxTB) break; @@ -373,7 +631,8 @@ float StFcsWaveformFitMaker::analyzeWaveform(int select, TGraphAsymmErrors* g, f case 4: integral = highest3(g, res); break; case 10: integral = gausFit(g, res, func, ped); break; case 11: integral = gausFitWithPed(g, res, func); break; - case 12: integral = PulseFit(g,res,func); break; + case 12: integral = PulseFit1(g,res,func); break; + case 13: integral = PulseFit2(g,res,func); break; case 21: integral = PedFit(g, res, func); break; case 31: integral = LedFit(g, res, func); break; default: @@ -421,10 +680,31 @@ float StFcsWaveformFitMaker::analyzeWaveform(int select, TGraphAsymmErrors* g, f if(integral>50 && dname.Contains(mFilter)) flag=1; } if(mFitDrawOn && flag && mFilename && mPage<=mMaxPage) { - printf("det=%s func=%d mFitDrawOn=%d mFilter=%s mFilename=%s mPage=%d mMaxPage=%d mHitIdx=%d integral=%f\n", - mDetName,func,mFitDrawOn,mFilter,mFilename,mPage,mMaxPage,mHitIdx,integral); + printf("hit:%u det=%s func=%p mFitDrawOn=%d mFilter=%s mFilename=%s mPage=%d mMaxPage=%d integral=%f\n", + mHitIdx,mDetName,(void*)func,mFitDrawOn,mFilter,mFilename,mPage,mMaxPage,integral); drawFit(g,func); } + + if( mOutFile!=0 ){ + if( mTest==1 ){ + int det = 0; int ch=0; + mDb->getFromName( g->GetName(), det,ch ); + int sumdep0 = StFcsPulseAna::SumDep0(g,centerTB()-3); + int sumdepmod = StFcsPulseAna::SumDep0Mod(g,centerTB()-3); + if( sumdep0>sumdepmod ){std::cout << " - StFcsWaveformFitMaker::analyzeWaveform|SumDep0:"<Fill(sumdepmod,sumdep0); + mH2_Sum8Dep0[det/2]->Fill(sumdep0,integral); + mH2_Sum8DepMod[det/2]->Fill(sumdepmod,integral); + } + if( mTest==7 ){ + int det0 = 0; int ch0=0; + mDb->getFromName( g->GetName(), det0,ch0 ); + mH1F_Res0[det0]->Fill(res[0]); + mH1F_Res0Zoom[det0]->Fill(res[0]); + mH1F_Res0[6]->Fill(res[0]); + mH1F_Res0Zoom[6]->Fill(res[0]); + } + } return integral; } @@ -546,34 +826,6 @@ float StFcsWaveformFitMaker::highest3(TGraphAsymmErrors* g, float* res){ return sum; //this is the 3 timebin sum } -//this one is just shower shape function -double StFcsWaveformFitMaker::pulseShape(double* x, double* p) { - double ret = p[0]*exp(-0.5*pow((x[0]-p[1])/p[2],2)); - if(mTail>0){ - double x1 = x[0] - p[1] - Xoff1; - if(x1>0){ - double a0 = p[0] * p[2]; - ret += a0*A1/Tau1/Tau1*pow(x1,P1)*exp(-x1/Tau1); - if(A2>0){ - double x2 = x[0] - p[1] - Xoff2; - if(x2>0){ - ret += a0*A2/Tau2/Tau2*pow(x2,P2)*exp(-x2/Tau2); - } - } - } - } - return ret; -} - -double StFcsWaveformFitMaker::multiPulseShape(double* x, double* p) { - int npulse = p[0]; - double ret = p[1]; - for(int i=0; i2) sprintf(Opt," "); @@ -621,7 +873,7 @@ float StFcsWaveformFitMaker::gausFit(TGraphAsymmErrors* g, float* res, TF1*& fun if( (adc1>adc0 || tb1==mMinTB) && adc1>=adc2){ //if falling from starting point, make it peak para[1+npeak*3+1] = adc1; para[1+npeak*3+2] = tb1; - para[1+npeak*3+3] = GSigma; + para[1+npeak*3+3] = mDbPulse->GSigma(); double dt = fabs(tb1 - mCenterTB); if(trgminSetData(g); } + Int_t compidx = mPulseFit->FoundPeakIndex(); + int myNpeaks = -1; + int npeakidx = -1; + int mypeakidx = -1; + if( mTest==2 ){ + mH1_NPeaksAkio->Fill(npeak); + if( trgx>=0 ){ mH1_NPeaksFilteredAkio->Fill(npeak); }//peaks inside triggered crossing + auto start=std::chrono::high_resolution_clock::now(); + if( compidx<0 ){ compidx = mPulseFit->AnalyzeForPeak(); } + auto stop=std::chrono::high_resolution_clock::now(); + long long usec = chrono::duration_cast(stop-start).count(); + //printf("Fit Time = %lld usec\n",usec); + mH1_PeakTiming->Fill(float(usec)/1000.0); + + myNpeaks = mPulseFit->NPeaks(); + mH1F_NPeaks[0]->Fill(myNpeaks); + if( 0<=compidx && compidxFill(myNpeaks); } + npeakidx = npeak<=4?npeak:5; + mypeakidx = myNpeaks<=4?myNpeaks:5; + //std::cout << "|npeakidx:"<< npeakidx << "|mypeakidx:"<< mypeakidx << std::endl; + //bool checkpeak = false; //This is to check a zero peak with a large adc value in triggered crossing + for( int i=0; iGetN(); ++i ){ + mH2F_AdcTbAkio[npeakidx]->Fill(t[i],a[i]); + mH2F_AdcTbMine[mypeakidx]->Fill(t[i],a[i]); + //if( mypeakidx==0 ){ if( 4020 ){ checkpeak=true; } } } + if( mPulseFit->ValidPeakIdx() ){mH2F_AdcTbValidPeak[mypeakidx]->Fill(t[i],a[i]);} + else{ mH2F_AdcTbValidPeak[6]->Fill(t[i],a[i]); } + //std::cout << "- |i:"<0 && npeakcreatePulse(mMinTB,mMaxTB,2+npeak*3);//Needs to have name "waveform" func->SetLineColor(6); func->SetParameters(para); func->FixParameter(0,npeak); @@ -645,14 +934,15 @@ float StFcsWaveformFitMaker::gausFit(TGraphAsymmErrors* g, float* res, TF1*& fun func->SetParLimits(j,para[j]-2.0,para[j]+2.0); //limit peak position to +- 2TB func->SetParLimits(1+i*3+3,0.5,10.0); //limit sigma to go too narrow or wide } - TFitResultPtr result = g->Fit("waveform",Opt,"",mMinTB,mMaxTB); + //TFitResultPtr result = g->Fit("waveform",Opt,"",mMinTB,mMaxTB); + TFitResultPtr result = g->Fit(func,Opt,"",mMinTB,mMaxTB); if(trgx>=0){ // return pulse closest to center of triggered xing res[1] = func->GetParameter(trgx*3 + 2); res[2] = func->GetParameter(trgx*3 + 3); res[3] = func->GetParameter(trgx*3 + 4); res[4] = func->GetChisquare()/func->GetNDF(); res[5] = npeak; - res[0] = res[1]*res[3]*sqrt2pi; + res[0] = res[1]*res[3]*StFcsDbPulse::sqrt2pi(); }else{ // no peak found in triggered xing //res[1] = 0.0; //res[2] = 0.0; @@ -664,10 +954,22 @@ float StFcsWaveformFitMaker::gausFit(TGraphAsymmErrors* g, float* res, TF1*& fun }else if(npeak>=mMaxPeak){ res[5] = npeak; sum8(g, res); // too many peak, taking sum8 - res[0]/=1.21; // normarized by 1/1.21 + res[0]/=1.21; // normalized by 1/1.21 LOG_INFO << Form("%s Finding too many peaks npeak=%d. Skip fitting. Taking Sum8",mDetName,npeak)<Fill(mPulseFit->SumWindow(),res[0]);} + if( npeakidx==mypeakidx){ + Double_t peakloc = 0; + if( 0<=compidx && compidxGetPeak(compidx).mPeakX; } + mH2F_APeakvMPeak[npeakidx]->Fill(peakloc,res[2]); + } + mH1F_PeakStart[mypeakidx]->Fill(mPulseFit->PeakStart()); + mH1F_PeakEnd[mypeakidx]->Fill(mPulseFit->PeakEnd()); + } return res[0]; } @@ -693,16 +995,35 @@ void StFcsWaveformFitMaker::drawFit(TGraphAsymmErrors* g, TF1* func){ mCanvas->cd(mPad); gg->SetTitle(mDetName); gg->Draw("ALP"); + if( func!=0 ){ + func->SetLineColor(kRed); + func->SetLineWidth(2); + func->DrawCopy("LSAME"); + } + if( mPulseFit!=0 ){ + int det=-1; + int ch=-1; + StFcsDb::getFromName(gg->GetName(),det,ch); + char post[50]; + sprintf(post,"_D%dC%d",det,ch); + StFcsPulseAna* ana = (StFcsPulseAna*)mPulseFit->DrawCopy("LP;A",post,gg);//Sets 'kCanDelete' so canvas will delete the copy + ana->GetData()->SetLineColor(kBlue); + ana->GetData()->SetMarkerStyle(kBlue); + ana->GetData()->SetMarkerStyle(4); + ana->GetData()->SetMarkerSize(0.5); + } + if(mPad==MAXPAD){ char file[100]; - if(mMaxPage==0) {sprintf(file,"%s.pdf",mFilename);} - else if(mPage==0) {sprintf(file,"%s.pdf(",mFilename);} - else if(mPage==mMaxPage) {sprintf(file,"%s.pdf)",mFilename); mPad=-1; mFitDrawOn=0; mHitIdx=0;} - else {sprintf(file,"%s.pdf",mFilename);} + sprintf(file,"%s.pdf",mFilename); + //if(mMaxPage==0) {sprintf(file,"%s.pdf",mFilename);} + //else if(mPage==0) {sprintf(file,"%s.pdf(",mFilename);} + //else if(mPage==mMaxPage) {sprintf(file,"%s.pdf",mFilename); mPad=-1; mFitDrawOn=0; mHitIdx=0;} + //else {sprintf(file,"%s.pdf",mFilename);} printf("Saving pdf with [%s] mPage=%d mPad=%d\n",file,mPage,mPad); - mCanvas->Print(file,"pdf"); + mCanvas->Print(file); // for(int i=0; iGetName(),det,id); + if( det>=0 && det==static_cast(detid) ){ + if( id>=0 && id==static_cast(ch) ){ + ggdraw->Draw("APL"); + ggdraw->SetLineColor(kBlack); + ggdraw->SetMarkerStyle(kBlack); + ggdraw->SetMarkerStyle(4); + ggdraw->SetMarkerSize(0.5); + ggdraw->GetXaxis()->SetTitle("timebin"); + ggdraw->GetYaxis()->SetTitle("ADC"); + if( mPulseFit!=0 ){ + int det=-1; + int ch=-1; + StFcsDb::getFromName(ggdraw->GetName(),det,ch); + char post[50]; + sprintf(post,"_D%dC%d",det,ch); + StFcsPulseAna* ana = mPulseFit->DrawCopy("LP;P",post,ggdraw);//Sets 'kCanDelete' so an external canvas will delete this object when "Clear" is called + ana->GetData()->SetLineColor(kBlue); + ana->GetData()->SetMarkerStyle(kBlue); + ana->GetData()->SetMarkerStyle(4); + ana->GetData()->SetMarkerSize(0.5); + //ana->ResetPeak(); + //ana->SetDebug(5); + //mPulseFit->Print(); + //ana->Print(); + } + } + } + } + +} + float StFcsWaveformFitMaker::gausFitWithPed(TGraphAsymmErrors* g, float* res, TF1*& func){ int n = g->GetN(); double *t = g->GetX(); double *a = g->GetY(); int p=0; + double sumsq = 0;//For variance for(int i=0; i=mPedMin && tb1<=mPedMax){ p+=a[i]; + sumsq += a[i]*a[i]; } } float ped = float(p)/(mPedMax-mPedMin+1.0); - printf("Pedestal=%6.2f/(%6.2f-%6.2f+1)=%6.2f\n",float(p),float(mPedMax),float(mPedMin),ped); + res[6] = ped; + res[7] = sqrt( ( sumsq-((double(p)*double(p))/(mPedMax-mPedMin+1.0)) )/(mPedMax-mPedMin) );//Variance/StdDev using naive algorithm from https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance + printf("Pedestal=%6.2f/(%6.2f-%6.2f+1)=%6.2f +- %6.2f\n",float(p),float(mPedMax),float(mPedMin),ped,res[7]); return gausFit(g, res, func, ped); } void StFcsWaveformFitMaker::drawRegion(int det, int col_low, int row_low, int col_high, int row_high, int event){ const int MAXPAD = 16; int NumDrawnPad[MAXPAD]={0};//To keep track of how much was drawn on each pad + //Keep track of minimum and maximum y for each pad x should be same so leave as is + Double_t MinY[MAXPAD] = {0}; + Double_t MaxY[MAXPAD] = {0}; int MaxDrawPad=0; //int MaxCol = mDb->nColumn(det); //int MaxRow = mDb->nRow(det); if( det<2 ){MaxDrawPad=6;}//Maximum to draw on one pad for Ecal else if( det>1 && det<4 ){MaxDrawPad=4;}//Maximum to draw on one pad for Hcal - else{ std::cout<< "ERROR:This function only works for det<4"<SetOptStat(0); @@ -745,7 +1111,6 @@ void StFcsWaveformFitMaker::drawRegion(int det, int col_low, int row_low, int co } mCanvas->Divide(4,4); - //int Nentries = mChWaveData.GetEntries();//Just to prevent multiple evaluation since ROOT loops through array to count elements (not needed Feb 05, 2021) //Need to loop to mHitIdx since that is the most hits for a given event and this function should be only used after all data is read in. This prevents reading old event that may have had more hits. for( unsigned int iCh=0; iChgetRowNumber(Chdet,Chid); if( !(col_low<=col && col<=col_high) || !(row_low<=row && row<=row_high) ){continue;}//Skip columns and rows not in range - mPad = PadNum4x4(det,col,row);//Won't return zero unless invalid column/row - TPad* padTemp = (TPad*)mCanvas->cd( mPad ); - + mPad = StFcsDbPulse::PadNum4x4(det,col,row);//Won't return zero unless invalid column/row + TPad* padTemp = (TPad*)mCanvas->cd( mPad-- );//Post-incrment should return old value //Since mPad is a number from 1-16 need to offset by one to get 0-15 for array - int color = 100-((static_cast(NumDrawnPad[mPad-1])/static_cast(MaxDrawPad))*(100-51));//Some suitable rainbow root colors (100=red, 51=purple) + int color = 100-((static_cast(NumDrawnPad[mPad])/static_cast(MaxDrawPad))*(100-51));//Some suitable rainbow root colors (100=red, 51=purple) gTemp->SetLineColor(color); gTemp->SetMarkerColor(color); - gTemp->SetMarkerStyle(NumDrawnPad[mPad-1]+24);//Some suitable root styles + gTemp->SetMarkerStyle(NumDrawnPad[mPad]+24);//Some suitable root styles gTemp->SetMarkerSize(0.5); - - //std::cout << "DEBUG::drawRegion|Pad:"<Draw("APL");} - else{ gTemp->Draw("PLsame"); } - if( (++NumDrawnPad[mPad-1]) > MaxDrawPad ){std::cout << "WARNING:Drawing too many on one pad"<cd(); + gTemp->Draw("APL"); + MinY[mPad]=padTemp->GetUymin(); + MaxY[mPad]=padTemp->GetUymax(); + } + else{ + gTemp->Draw("PL"); + Double_t YminNew=5000; Double_t YmaxNew=-5; + StFcsDbPulse::getYMinMax(gTemp,YminNew,YmaxNew); + if( YminNewMaxY[mPad] ){MaxY[mPad]=YmaxNew;} + ((TGraphAsymmErrors*)padTemp->GetListOfPrimitives()->At(1))->GetYaxis()->SetRangeUser(MinY[mPad],MaxY[mPad]);//For this function first graph in list of primitives is the graph with axis object (First object of list is TFrame) + padTemp->Draw();//Update axes + } + if( (++NumDrawnPad[mPad]) > MaxDrawPad ){LOG_WARN << "StFcsWaveformFitMaker::drawRegion(): Drawing too many on one pad"<SetSignal(gae);//Resets finder - mPulseFit->SetSearchWindow(mCenterTB,5);//Some suitable starting search parameter may need to be adjusted based on data - mPulseFit->SetZS();//This will set baseline to zero with some "good" guess for the pedestal sigma - //mPulseFit->SetBaseline(0,0.75);//Need to set baseline to zero, or greater than zero, to prevent baseline finding. - //This is good for zero-suppressed data but should be left unset otherwise. Also can be used to adjust adc threshold - //for acceptance by changing the 0.75 to some other number (default ADC acceptance = 4.0*0.75) - mPulseFit->SetSearchWindow(mCenterTB,3);//Check +- 3tb around triggered crossing - //Many more options to set but excluded for now to see how default values perform - res[0]=mPulseFit->SumAdcFit(); //this is full integral (pedestal subtracted) - //int xlow = mPulseFit->SignalStart(); - //int xhigh = mPulseFit->SignalEnd(); - func = mPulseFit->SignalFit(); //Fitted function - //mPulseFit->GetXYMax(xlow,xhigh); //Get x and y max from graph not fit - res[1]=func->GetParameter(0); //this is peak height from signal fit - //res[1] = mPulseFit->MaxAdc(); - res[2]=func->GetParameter(1); //this is peak position from signal fit [timebin] - //res[2] = mPulseFit->MaxTb(); - res[3]=func->GetParameter(2); //this is width from signal fit - res[4]=func->GetChisquare()/func->GetNDF(); //chi2 from signal fit - return res[0]; +float StFcsWaveformFitMaker::PulseFit1(TGraphAsymmErrors* gae, float* res, TF1*& func){ + if( mPulseFit==0 ){mPulseFit = new StFcsPulseAna(gae); SetupDavidFitterMay2022();} + else{mPulseFit->SetData(gae);}//Resets finder + if( GetDebug()>2 ){mPulseFit->SetDebug(1);} + + Int_t compidx = mPulseFit->FoundPeakIndex(); + if( compidx<0 ){ compidx = mPulseFit->AnalyzeForPeak(); } + int npeaks = mPulseFit->NPeaks(); + + int det0 = 0; int ch0=0; + mDb->getFromName( gae->GetName(), det0,ch0 ); + + if( mTest==4 ){ + if( compidx==npeaks ){ res[0] = 0; return res[0]; } + func = mDbPulse->createPulse(mMinTB,mMaxTB,2+npeaks*3); + mPulseFit->SetFitPars(func); + TFitResultPtr result = gae->Fit(func,"BNRQ"); + for( int i=0; iFill( func->GetParameter(1+i*3+3),func->GetParameter(1+i*3+1) ); + } + res[5] = npeaks; + res[1] = func->GetParameter(compidx*3 + 2); + res[2] = func->GetParameter(compidx*3 + 3); + res[3] = func->GetParameter(compidx*3 + 4); + res[4] = func->GetChisquare()/func->GetNDF(); + res[0] = res[1]*res[3]*StFcsDbPulse::sqrt2pi(); + mH2_HeightvsSigmaTrig->Fill(res[3],res[1]); + mH1_ChiNdf->Fill(res[4]); + mH2_HeightvsChiNdf->Fill(res[4],res[1]); + mH2_MeanvsChiNdf->Fill(res[4],res[2]); + mH2_SigmavsChiNdf->Fill(res[4],res[3]); + + return res[0]; + } + + //std::cout << "|compidx:"<Fill(npeaks); + mH1F_NPeaks[6]->Fill(npeaks); + } + if( compidx==npeaks ){ res[0] = 0; }//no found peak case sum is 0 + else if( npeaks<=1 ){ //0 or 1 peak case with a valid peak; do sum 8 + if( mTest==3 ){ + mH1F_NPeaksFiltered[det0]->Fill(npeaks); + mH1F_NPeaksFiltered[6]->Fill(npeaks); + } + res[0] = sum8(gae,res); + //Scale sum8 to match fitted sum (These may need to be confirmed by data year by year) + if( det0==0 || det0==1 ){res[0]/=1.226;} + if( det0==2 || det0==3 ){res[0]/=1.195;} + if( det0==4 || det0==5 ){res[0]/=1.29;} + res[1] = mPulseFit->GetPeak(compidx).mPeakY; + res[2] = mPulseFit->GetPeak(compidx).mPeakX; + res[3] = mDbPulse->GSigma(); + res[4] = -1; + } + else{ //More than 1 peak and valid peak found + if( mTest==3 ){ + mH1F_NPeaksFiltered[det0]->Fill(npeaks); + mH1F_NPeaksFiltered[6]->Fill(npeaks); + //int mypeakidx = (npeaks<=5)?(npeaks-2):3; + double xfound = mPulseFit->GetPeak(compidx).mPeakX; + double yfound = mPulseFit->GetPeak(compidx).mPeakY; + for( Int_t i=0; iGetPeak(i).mPeakX; + double y = mPulseFit->GetPeak(i).mPeakY; + mH2_NPeakvsPeakXdiff->Fill(xfound-x,npeaks); + mH2_NPeakvsPeakYratio->Fill(yfound/y,npeaks); + } + } + int numoverlaps = 0; + for( Int_t i=0; iGetPeak(compidx),mPulseFit->GetPeak(i) ); + mH1_VOverlap->Fill(comp); + mH2_VvsNOverlap->Fill(i-compidx,comp); + if( comp!=0 ){ ++numoverlaps; } + } + else{ if(PeakCompare( mPulseFit->GetPeak(compidx),mPulseFit->GetPeak(i) )!=0 ){ ++numoverlaps; } } + } + if( numoverlaps>0 ){ + //res[0] = mPulseFit->FitSignal(mMinTb,mMaxTb); + auto start=std::chrono::high_resolution_clock::now();//for timing studies can be commented out as long as not testing + func = mDbPulse->createPulse(mMinTB,mMaxTB,2+npeaks*3); + mPulseFit->SetFitPars(func); + //mPulseFit->SetSignal(func); + TFitResultPtr result = gae->Fit(func,"BNRQ"); + res[1] = func->GetParameter(compidx*3 + 2); + res[2] = func->GetParameter(compidx*3 + 3); + res[3] = func->GetParameter(compidx*3 + 4); + res[4] = func->GetChisquare()/func->GetNDF(); + res[0] = res[1]*res[3]*StFcsDbPulse::sqrt2pi(); + if( mTest==3 ){ + auto stop=std::chrono::high_resolution_clock::now(); + long long usec = chrono::duration_cast(stop-start).count(); + mH1_TimeFitPulse->Fill(float(usec)/1000.0); + mH2F_NOvsId[det0]->Fill(ch0,numoverlaps); + mH2_NOvsNPeaks->Fill(npeaks,numoverlaps); + //printf("Fit Time = %lld usec\n",usec); + } + } + else{//Don't need to fit as other peaks don't contribute to found peak + res[0] = sum8(gae,res); + //Scale sum8 to match fitted sum (These may need to be confirmed by data year by year) + if( det0==0 || det0==1 ){res[0]/=1.226;} + if( det0==2 || det0==3 ){res[0]/=1.195;} + if( det0==4 || det0==5 ){res[0]/=1.29;} + res[1] = mPulseFit->GetPeak(compidx).mPeakY; + res[2] = mPulseFit->GetPeak(compidx).mPeakX; + res[3] = mDbPulse->GSigma(); + res[4] = -1; + } + } + res[5] = npeaks;//change to number of peaks in pre2 and post1 (number of peaks fitted) + if( mTest==3 ){ + mH1F_Res0[det0]->Fill(res[0]); + mH1F_Res0Zoom[det0]->Fill(res[0]); + mH1F_Res0[6]->Fill(res[0]); + mH1F_Res0Zoom[6]->Fill(res[0]); + float sum8res[8] = {0}; + sum8( gae ,sum8res ); + mH1F_Sum8Res0[det0]->Fill(sum8res[0]); + mH1F_Sum8Res0Zoom[det0]->Fill(sum8res[0]); + mH1F_Sum8Res0[6]->Fill(sum8res[0]); + mH1F_Sum8Res0Zoom[6]->Fill(sum8res[0]); + float savesum8 = sum8res[0]; + //Scale sum8 to match fitted sum + //if( det0==0 || det0==1 ){savesum8/=1.226;} + //if( det0==2 || det0==3 ){savesum8/=1.195;} + //if( det0==4 || det0==5 ){savesum8/=1.29;} + + if( npeaks>=1 ){//No found peak skip fitting + if( func==0 ){//No fit performed above so do a fit now + //Here I will reuse the sum8res array + func = mDbPulse->createPulse(mMinTB,mMaxTB,2+npeaks*3); + mPulseFit->SetFitPars(func); + TFitResultPtr result = gae->Fit(func,"BNRQ"); + sum8res[5] = npeaks; + sum8res[1] = func->GetParameter(compidx*3 + 2); + sum8res[2] = func->GetParameter(compidx*3 + 3); + sum8res[3] = func->GetParameter(compidx*3 + 4); + sum8res[4] = func->GetChisquare()/func->GetNDF(); + sum8res[0] = sum8res[1]*sum8res[3]*StFcsDbPulse::sqrt2pi(); + mH1F_FitRes0[det0]->Fill(sum8res[0]); + mH1F_FitRes0Zoom[det0]->Fill(sum8res[0]); + mH1F_FitRes0[6]->Fill(sum8res[0]); + mH1F_FitRes0Zoom[6]->Fill(sum8res[0]); + + mH2F_Sum8vFit[det0]->Fill(sum8res[0],savesum8); + mH2F_Sum8vFit[6]->Fill(sum8res[0],savesum8); + } + else{//Fit was done above so just save that result + mH1F_FitRes0[det0]->Fill(res[0]); + mH1F_FitRes0Zoom[det0]->Fill(res[0]); + mH1F_FitRes0[6]->Fill(res[0]); + mH1F_FitRes0Zoom[6]->Fill(res[0]); + mH2F_Sum8vFit[det0]->Fill(res[0],savesum8); + mH2F_Sum8vFit[6]->Fill(res[0],savesum8); + } + } + } + + return res[0]; +} + +float StFcsWaveformFitMaker::PulseFit2(TGraphAsymmErrors* gae, float* res, TF1*& func){ + if( mPulseFit==0 ){mPulseFit = new StFcsPulseAna(gae); SetupDavidFitterMay2022();} + else{mPulseFit->SetData(gae);}//Resets finder + if( GetDebug()>2 ){mPulseFit->SetDebug(1);} + + Int_t compidx = mPulseFit->FoundPeakIndex(); + if( compidx<0 ){ compidx = mPulseFit->AnalyzeForPeak(); } + int npeaks = mPulseFit->NPeaks(); + + int det0 = 0; int ch0=0; + mDb->getFromName( gae->GetName(), det0,ch0 ); + + if( mTest==6 ){ + mH1F_NPeaks[det0]->Fill(npeaks); + mH1F_NPeaks[6]->Fill(npeaks); + } + if( compidx==npeaks ){ res[0] = 0; }//no found peak case sum is 0 + else if( npeaks<=1 ){ //0 or 1 peak case with a valid peak; do sum 8 + res[0] = sum8(gae,res); + //Scale sum8 to match fitted sum (These may need to be confirmed by data year by year) + if( det0==0 || det0==1 ){res[0]/=1.226;} + if( det0==2 || det0==3 ){res[0]/=1.195;} + if( det0==4 || det0==5 ){res[0]/=1.29;} + res[1] = mPulseFit->GetPeak(compidx).mPeakY; + res[2] = mPulseFit->GetPeak(compidx).mPeakX; + res[3] = mDbPulse->GSigma(); + res[4] = -1; + } + else{ //More than 1 peak and valid peak found + Double_t xmin, xmax; + int trigfitidx = compidx;//this is starting condition is needed to pick out the right index + npeaks = NPeaksPre2Post1(trigfitidx,xmin,xmax);//Check for peaks near triggered crossing (counts triggered crosing peak) + if( mTest==6 || mTest==7 ){ + mH1F_NPeaksFiltered[det0]->Fill(npeaks); + mH1F_NPeaksFiltered[6]->Fill(npeaks); + } + if( npeaks>1 ){//If equal to one then still don't need to fit + //res[0] = mPulseFit->FitSignal(mMinTb,mMaxTb); + auto start=std::chrono::high_resolution_clock::now();//for timing studies can be commented out as long as not testing + func = mDbPulse->createPulse(xmin,xmax,2+npeaks*3);//only fit inside the range of valid peaks + mPulseFit->SetFitPars(func); + //mPulseFit->SetSignal(func); + TFitResultPtr result = gae->Fit(func,"BNRQ"); + //compidx no longer equal to index to triggered crossing peak + res[1] = func->GetParameter(trigfitidx*3 + 2); + res[2] = func->GetParameter(trigfitidx*3 + 3); + res[3] = func->GetParameter(trigfitidx*3 + 4); + res[4] = func->GetChisquare()/func->GetNDF(); + res[0] = res[1]*res[3]*StFcsDbPulse::sqrt2pi(); + if( mTest==6 ){ + auto stop=std::chrono::high_resolution_clock::now(); + long long usec = chrono::duration_cast(stop-start).count(); + mH1_TimeFitPulse->Fill(float(usec)/1000.0); + //printf("Fit Time = %lld usec\n",usec); + } + } + else{//Don't need to fit as other peaks don't contribute to found peak + res[0] = sum8(gae,res); + //Scale sum8 to match fitted sum (These may need to be confirmed by data year by year) + if( det0==0 || det0==1 ){res[0]/=1.226;} + if( det0==2 || det0==3 ){res[0]/=1.195;} + if( det0==4 || det0==5 ){res[0]/=1.29;} + res[1] = mPulseFit->GetPeak(compidx).mPeakY; + res[2] = mPulseFit->GetPeak(compidx).mPeakX; + res[3] = mDbPulse->GSigma(); + res[4] = -1; + } + } + res[5] = npeaks;//number of peaks in pre2 and post1 (number of fitted peaks) + + if( mTest==6 ){ + mH1F_Res0[det0]->Fill(res[0]); + mH1F_Res0Zoom[det0]->Fill(res[0]); + mH1F_Res0[6]->Fill(res[0]); + mH1F_Res0Zoom[6]->Fill(res[0]); + float sum8res[8] = {0}; + sum8( gae ,sum8res ); + mH1F_Sum8Res0[det0]->Fill(sum8res[0]); + mH1F_Sum8Res0Zoom[det0]->Fill(sum8res[0]); + mH1F_Sum8Res0[6]->Fill(sum8res[0]); + mH1F_Sum8Res0Zoom[6]->Fill(sum8res[0]); + float savesum8 = sum8res[0]; + //if( det0==0 || det0==1 ){savesum8/=1.226;} + //if( det0==2 || det0==3 ){savesum8/=1.195;} + //if( det0==4 || det0==5 ){savesum8/=1.29;} + + if( mPulseFit->NPeaks()>=1 ){//No found peak skip fitting + TF1* testfunc = mDbPulse->createPulse(mMinTB,mMaxTB,2+(mPulseFit->NPeaks())*3);//Want to compare when doing fit to all peaks + //Here I will reuse the sum8res array + mPulseFit->SetFitPars(testfunc); + TFitResultPtr result = gae->Fit(testfunc,"BNRQ"); + sum8res[5] = mPulseFit->NPeaks(); + sum8res[1] = testfunc->GetParameter(compidx*3 + 2); + sum8res[2] = testfunc->GetParameter(compidx*3 + 3); + sum8res[3] = testfunc->GetParameter(compidx*3 + 4); + sum8res[4] = testfunc->GetChisquare()/testfunc->GetNDF(); + sum8res[0] = sum8res[1]*sum8res[3]*StFcsDbPulse::sqrt2pi(); + mH1F_FitRes0[det0]->Fill(sum8res[0]); + mH1F_FitRes0Zoom[det0]->Fill(sum8res[0]); + mH1F_FitRes0[6]->Fill(sum8res[0]); + mH1F_FitRes0Zoom[6]->Fill(sum8res[0]); + + mH2F_Sum8vFit[det0]->Fill(sum8res[0],savesum8); + mH2F_Sum8vFit[6]->Fill(sum8res[0],savesum8); + } + } + + return res[0]; } float StFcsWaveformFitMaker::PedFit(TGraphAsymmErrors* gae, float* res, TF1*& func){ - //StFcsPulseFit Fitter(gae); - mPulseFit->SetSignal(gae); + if( mPulseFit==0 ){mPulseFit = new StFcsPulseAna(gae);} + else{mPulseFit->SetData(gae);}//Resets finder + if( GetDebug()>2 ){mPulseFit->SetDebug(1);} + + mPulseFit->SetRange(-4,0,2000,5000); + mPulseFit->SetBaselineFit(func);//Need to do this since func gets deleted outside + mPulseFit->SetData(gae); mPulseFit->AnalyzeForPedestal();//Must be called before anything else otherwise crash func = mPulseFit->BaselineFit(); res[0] = mPulseFit->Baseline(); res[1] = mPulseFit->BaselineSigma(); + if( func==0 ){ res[2]=0.0; return res[0];} res[2] = func->GetChisquare()/func->GetNDF(); //chi2 from signal fit return res[0]; } float StFcsWaveformFitMaker::LedFit(TGraphAsymmErrors* gae, float* res, TF1*& func) { - //StFcsPulseFit Fitter(gae); - mPulseFit->SetSignal(gae); - mPulseFit->SetSearchWindow(mCenterTB,5);//Some suitable starting search parameter may need to be adjusted based on data + if( mPulseFit==0 ){mPulseFit = new StFcsPulseAna(gae);} + else{mPulseFit->SetData(gae);}//Resets finder + if( GetDebug()>2){mPulseFit->SetDebug(1);} + + mPulseFit->SetBaselineSigmaScale(5); + mPulseFit->SetRange(-4,0,2000,5000); + mPulseFit->SetSearchWindow(centerTB(),4);//Check +- 4tb around triggered crossing. Mmay need to be adjusted based on data //Since LED data is raw (not pedestal subtracted) determine baseline first mPulseFit->AnalyzeForPedestal();//Must be called before anything else otherwise crash - func = mPulseFit->BaselineFit(); - res[5] = mPulseFit->Baseline(); - res[6] = mPulseFit->BaselineSigma(); - res[7] = func->GetChisquare()/func->GetNDF(); //chi2 from signal fit - - mPulseFit->SetSearchWindow(mCenterTB,4);//Check +- 4tb around triggered crossing - //Many more options to set but excluded for now to see how default values perform - res[0]=mPulseFit->SumAdcFit(); //this is full integral (pedestal subtracted) - func = mPulseFit->SignalFit(); //Fitted function - res[1]=func->GetParameter(0); //this is peak height from signal fit - res[2]=func->GetParameter(1); //this is peak position from signal fit [timebin] - res[3]=func->GetParameter(2); //this is width from signal fit - res[4]=func->GetChisquare()/func->GetNDF(); //chi2 from signal fit - + + Int_t compidx = mPulseFit->FoundPeakIndex(); + if( compidx<0 ){ compidx = mPulseFit->AnalyzeForPeak(); } + int npeaks = mPulseFit->NPeaks(); + + if( compidx==npeaks ){ res[0] = 0; }//no found peak case sum is 0 + else{//At least 1 peak found + func = mDbPulse->createPulse(mMinTB,mMaxTB,2+npeaks*3); + mPulseFit->SetFitPars(func); + TFitResultPtr result = gae->Fit(func,"BNRQ"); + res[5] = npeaks; + res[1] = func->GetParameter(compidx*3 + 2); + res[2] = func->GetParameter(compidx*3 + 3); + res[3] = func->GetParameter(compidx*3 + 4); + res[4] = func->GetChisquare()/func->GetNDF(); + res[0] = res[1]*res[3]*StFcsDbPulse::sqrt2pi(); + } return res[0]; } @@ -894,9 +1536,80 @@ int StFcsWaveformFitMaker::PadNum4x4(int det, int col, int row) ncol = 2; nrow = 2; } - else{ std::cout << "This only works for Ecal and Hcal" << std::endl; return 0;} + else{ LOG_ERROR << "StFcsWaveformFitMaker::PadNum4x4: This only works for Ecal and Hcal" << endm; return 0;} int padcol = GenericPadPos(col,ncol,4); int padrow = GenericPadPos(row,nrow,4); return 4*(padrow-1)+padcol; } +int StFcsWaveformFitMaker::PeakCompare(const PeakWindow& pwin1, const PeakWindow& pwin2 ) +{ + //return true; + double xdiff = fabs(pwin1.mPeakX-pwin2.mPeakX); + double yratio = pwin1.mPeakY/pwin2.mPeakY; + unsigned int result = 0;//Bit vector to store results of checks + // + First bit xdiff cut + // + Second bit yratio cut + //@[June 22, 2022](David Kapukchyan)>After looking at some plots from Test==3 an xdiff>10 seems like a good cutoff for no overlap between peaks and a y-ratio>2.0 + //check peak-y difference if greater than some value no overlap + //check sum8 of both and if both small no overlap + if( xdiff<10.0 ){ result |= 1<<0; } + if( yratio<2.0 ){ result |= 1<<1; } + return result; +} + +int StFcsWaveformFitMaker::NPeaksPre2Post1(int& trigidx,Double_t& xmin, Double_t &xmax) const +{ + //make initial window 1 RHIC crossing + xmin = mCenterTB - mDbPulse->TBPerRC()/2; + xmax = mCenterTB + mDbPulse->TBPerRC()/2; + int npeaksxing=0; + //@[August 25, 2022](David Kapukchyan)>In future use varaiables as opposed to constants?? + int prexing = 3*mDbPulse->TBPerRC();//Pre-crossing more important than post-crossing so use wider range for pre-crossing + int postxing = 2*mDbPulse->TBPerRC(); + bool foundtrigidx = false; + for( int i=0; iNPeaks(); ++i ){ + if( (mCenterTB-prexing) < mPulseFit->GetPeak(i).mPeakX && mPulseFit->GetPeak(i).mPeakX < (mCenterTB+postxing) ){ + //if pre crossing peak hits zero then don't need to fit pre-crossing (this can be done by checking if pre-endx==trig-startx?? + if( !foundtrigidx && trigidx==i ){ foundtrigidx = true; trigidx = npeaksxing; }//Should be before ++npeaksxing + ++npeaksxing; + if( xmin > mPulseFit->GetPeak(i).mStartX ){ xmin = mPulseFit->GetPeak(i).mStartX; } + if( xmax < mPulseFit->GetPeak(i).mEndX ){ xmax = mPulseFit->GetPeak(i).mEndX; } + } + } + if( !foundtrigidx ){ LOG_ERROR << "StFcsWaveformFitMaker::NPeakPre2Post1: Unable to find a matching triggered crossing possibly due to improper use of function" << endm; } + return npeaksxing; +} + +StFcsPulseAna* StFcsWaveformFitMaker::InitFitter(Double_t ped) +{ + if( mPulseFit==0 ){ mPulseFit = new StFcsPulseAna(); } + //Need to set baseline to zero, or greater than zero, to prevent baseline finding. + //This is good for zero-suppressed data but should be left unset otherwise. Also can be used to adjust adc threshold + //for acceptance by changing the 0.75 to some other number (default ADC acceptance = 4.0*0.75) + mPulseFit->SetBaseline(ped,0.39); + mPulseFit->SetBaselineSigmaScale(5); + mPulseFit->SetRange(-4,0,2000,5000); + mPulseFit->SetSearchWindow(centerTB(),4);//Check +- 4tb around triggered crossing; some suitable starting search parameter may need to be adjusted based on data + mPulseFit->SetContinuity(1.0); + return mPulseFit; +} + +void StFcsWaveformFitMaker::SetupDavidFitterMay2022(Double_t ped) +{ + if( mPulseFit==0 ){return;} + InitFitter(ped); + //mPulseFit->SetFilter(1,1);//Mean filter with "radius" 1 + //mPulseFit->SetFilter(1,2);//Mean filter with "radius" 2 + //mPulseFit->SetFilter(2,1,0.5);//Gaus filter with "radius" 1 and sigma 0.5. + //mPulseFit->SetFilter(2,3,1);//Gaus filter with "radius" 2 and sigma 1. + //Parameters for peak tunneling method, + mPulseFit->SetTunnelScale(1.0); + mPulseFit->SetTunnelSigma(1.5); + mPulseFit->SetTunnelThreshold(-0.001);//Negative value <=-1 means do merging after peak finding, positive value <=1 means do merging wiht peak finding + //Below are some other thresholds for testing + //mPulseFit->SetTunnelThreshold(-0.9);//Negative value >=-1 means do merging after peak finding + //mPulseFit->SetTunnelThreshold(0.1);//Negative value >=-1 means do merging after peak finding + //mPulseFit->SetTunnelThreshold(-2.0);//Negative value >=-1 means do merging after peak finding +} + diff --git a/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.h b/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.h index e50f164d041..93405c8d91f 100644 --- a/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.h +++ b/StRoot/StFcsWaveformFitMaker/StFcsWaveformFitMaker.h @@ -52,57 +52,60 @@ #define STROOT_STFCSWAVEFORMFITMAKER_STFCSWAVEFORMFITMAKER_H_ #include "StMaker.h" -#include "StFcsPulseFit.h" +#include "StFcsPulseAna.h" + +#include "TH2.h" class StFcsCollection; class StFcsHit; class StFcsDb; +class StFcsDbPulse; class TGraphAsymmErrors; class TGraph; class TCanvas; -class TH1F; -class TH2F; class StFcsWaveformFitMaker : public StMaker { public: StFcsWaveformFitMaker(const char* name = "StFcsWaveformFitMaker"); ~StFcsWaveformFitMaker(); - int InitRun(int runNumber); - int Make(); - int Finish(); - void Clear(Option_t* option = ""); + 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); //!< Set test level. Intended to be used for single files. Output file name can be changed with #writeFile() void setEnergySelect(int ecal=10, int hcal=10, int pres=1) {mEnergySelect[0]=ecal; mEnergySelect[1]=hcal; mEnergySelect[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) {mTail=v;} + void setTail(int v); void setMaxPeak(int v) {mMaxPeak=v;} void setPedTimeBins(int min, int max) {mPedMin=min; mPedMax=max;} - // Create or Reset TGraphAsymmErrors (at mHitIdx) - TGraphAsymmErrors* resetGraph(); - + TGraphAsymmErrors* resetGraph(); //! Create or Reset TGraphAsymmErrors (at mHitIdx) + // Return TGraphAsymmErrors at idx // if idx<0 (default) it returns "current" which is idx=mHitIdx-1 // makeTGraphAsymmErrors() makes all at idx=mHitIdx=0 except when // mFitDrawOn=1, then keep incrementing mHitIdx to hold them for a event // mFitDrawOn=2, then keep incrementing mHitIdx to hold them until new page TGraphAsymmErrors* getGraph(int idx=-1); + TGraphAsymmErrors* getGraph(int det, int id); //measuring fit time void setMeasureTime(char* file) {mMeasureTime=file;} - TH1F* mTime; + TH1F* mTime = 0; //stage0 peak algo study TH2F* mTimeIntg[4]; //create makeTGraphAsymmErrors from given timebin and adc data - //with asymmetroc errors when adc is saturated + //with asymmetric errors when adc is saturated //Places graph at index 0 of internal TClonesArray and will be deleted in destructor TGraphAsymmErrors* makeTGraphAsymmErrors(int n, double* t, double* adc); TGraphAsymmErrors* makeTGraphAsymmErrors(TGraph* g); @@ -128,7 +131,8 @@ class StFcsWaveformFitMaker : public StMaker { float highest3(TGraphAsymmErrors* g, float *res); //! mEnergySelect=4 float gausFit (TGraphAsymmErrors* g, float *res, TF1*& f, float ped=0.0); //! mEnergySelect=10 float gausFitWithPed (TGraphAsymmErrors* g, float *res, TF1*& f); //! mEnergySelect=11 - float PulseFit(TGraphAsymmErrors* g, float* res, TF1*& f); //! mEnergySelect=12 + float PulseFit1(TGraphAsymmErrors* g, float* res, TF1*& f); //! mEnergySelect=12 + float PulseFit2(TGraphAsymmErrors* g, float* res, TF1*& f); //! mEnergySelect=13 //Pedestal Analysis methods //res[0] pedestal Value @@ -145,40 +149,107 @@ class StFcsWaveformFitMaker : public StMaker { // res[6] pedestal Sigma // res[7] pedestal Chi2/NDF float LedFit( TGraphAsymmErrors* g, float* res, TF1*& f);//! mEnergySelect=31 (TF1 is for LED pulse fit) - - //pulse shape functions - double pulseShape(double* x, double* p); - double multiPulseShape(double* x, double* p); - + //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); //!< Use to change the name of the file where test histograms will be saved void setFitDrawOn(int v=1) {mFitDrawOn=v;} //=1 to keep for a event, =2 for a page void setFitFilter(char* filter) {mFilter=filter; mFitDrawOn=2;} //Draw from David - StFcsPulseFit* davidFitter(){return mPulseFit;} + int centerTB()const{return mCenterTB;} + void setDavidFitter(StFcsPulseAna* v){if(mPulseFit==0){mPulseFit=v;}} + StFcsPulseAna* davidFitter(){return mPulseFit;} + StFcsPulseAna* InitFitter(Double_t ped=0); //! Sets up basic values needed by #StFcsPulseAna 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 printArray() const; + void drawFitter(Option_t* opt){ if(mPulseFit!=0){mPulseFit->Draw(opt);} } + void drawCh(UInt_t detid, UInt_t ch) const; protected: TClonesArray mChWaveData; //Contains all graph data void drawFit(TGraphAsymmErrors* g, TF1* func); - StFcsPulseFit* mPulseFit; + StFcsPulseAna* mPulseFit; + + /**@brief Variable to use when testing StFcsWaveformFitMaker algorithms + + Self contained analysis for testing various components/functions of #StFcsWaveformFitMaker + + - 0 = no testing + - 1 = test DEP algorithm + - 2 = test PeakAna vs. gausFit + - 3 = test PulseFit1 picking sum method + - 4 = test PulseFit1 all data with peaks + - 5 = test timing of gausFit() vs. PulseFit1() + - 6 = like test==3 but for PulseFit2() + - 7 = test PulseFit2() for overall quality doesn't include preshower + */ + int mTest = 0; + + TFile* mOutFile;//Root output file for testing + //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) + TH1F* mH1_NPeaksAkio = 0; //Number of peaks found by gausFit + TH1F* mH1_NPeaksFilteredAkio = 0; //Number of peaks found by gausFit for signals that had a triggered crossing + TH2F* mH2F_AdcTbAkio[6]; //Adc vs. Tb for different number of peaks Akio method + TH2F* mH2F_AdcTbMine[6]; //Adc vs. Tb for different number of peaks my method + TH2F* mH2F_AdcTbValidPeak[7];//Adc vs. Tb from my algorithm that had a peak at centerTb (Need an extra one for non-valid peaks) + TH2F* mH2F_SumFitvSumWin[6]; //Sum from Akio's Fit function vs. Sum from my found peak window for different number of peaks + TH2F* mH2F_APeakvMPeak[6]; //PeakLocations from Akio vs. Mine + TH1F* mH1F_PeakStart[6]; //PeakWindow Starting x-values + TH1F* mH1F_PeakEnd[6]; //PeakWindow Ending x-values + TH1F* mH1_PeakTiming = 0; //Timing for just peak finding. + + TH1F* mH1F_NPeaks[7]; //Number of peaks found by peak finder + TH1F* mH1F_NPeaksFiltered[7]; //Number of peaks for cases where a valid peak was found + TH2F* mH2_NPeakvsPeakXdiff = 0; //Number of peaks vs. Peak X diff + TH2F* mH2_NPeakvsPeakYratio = 0;//Number of peaks vs. Peak Y ratio + TH1F* mH1_VOverlap = 0; //Value of overlap + TH2F* mH2_NOvsNPeaks = 0; //NO (Number of overlaps) vs. Number of peaks + TH2F* mH2_VvsNOverlap = 0; //Compare value for that peak comparison vs. Overlap index + TH2F* mH2F_NOvsId[6]; //NO (number of overlaps) vs channel id for the 6 detector ids + TH1F* mH1F_Res0[7]; //Histogram of all res[0] regardless of method + TH1F* mH1F_Res0Zoom[7]; //Histogram of all res[0] with finer bining at low end + TH1F* mH1F_Sum8Res0[7]; //Histogram of "res[0]" using sum8 regardless of method called + TH1F* mH1F_Sum8Res0Zoom[7]; //same as above histogram with finer bining at low end + TH1F* mH1F_FitRes0[7]; //Histogram of "res[0]" from fit regardless of method called + TH1F* mH1F_FitRes0Zoom[7]; //same as above histogram with finer bining at low end + TH2F* mH2F_Sum8vFit[7]; //Histograms of Fit res[0] vs. sum8 res[0] + TH1F* mH1_TimeFitPulse = 0; //Histogram to time how long just the fitting takes in PulseFit() + + //For testing peak height vs. sigma + TH2F* mH2_HeightvsSigma = 0; //Histogram of all fitted peak heights vs. their sigma + TH2F* mH2_HeightvsSigmaTrig = 0; //Histogram of height of fitted peaks in triggered crossing vs. their sigma + TH1F* mH1_ChiNdf = 0; //Histogram of chi^2/ndf for all fits + TH2F* mH2_HeightvsChiNdf = 0; //Histogram of height vs. chi^2/ndf for all fits + TH2F* mH2_MeanvsChiNdf = 0; //Histogram of height vs. chi^2/ndf for all fits + TH2F* mH2_SigmavsChiNdf = 0; //Histogram of chi^2/ndf for all fits + + TH1F* mH1_PeakTimingGaus = 0; //Histogram to test timing of gausFit() + TH1F* mH1_PeakTimingPuls = 0; //Histogram to test timing of PulseFit() + TH2F* mH2_PeakTimingCompare = 0; //Histogram to test timing between gausFit() and PulseFit() + + void SetupDavidFitterMay2022(Double_t ped=0); //! This special function is used to set all the parameters for #StFcsPulseAna based on cosmic and Run 22 data. It is intended to be used only for Run 22 data + int PeakCompare(const PeakWindow& pwin1, const PeakWindow& pwin2 ); //Compare if two peaks overlap and return a bit vector of tests passed/failed for comparing pwin1 to pwin2. 0 means all tests passed and pwin1 does not overlap with pwin2 + int NPeaksPre2Post1(int& trigidx, Double_t& xmin, Double_t& xmax) const;//xmin and xmax will be the range of the pre-crossing -2 and post-crossing +1 peaks. trigidx is needed to pick up the triggered crossing in the new number of peaks private: StFcsDb* mDb=0; //! pointer to fcsDb + StFcsDbPulse* mDbPulse = 0; //! pointer to fcsPulse StFcsCollection* mFcsCollection=0; //! pointer to fcsCollection unsigned int mHitIdx = 0; //running index for the TClonesArray //Figures out error to set on the TGraphAsymmErrors for a given point and adc - void setTGraphAsymmErrors(TGraphAsymmErrors* gae, const int &i, const double &adc); - + //void setTGraphAsymmErrors(TGraphAsymmErrors* gae, const int &i, const double &adc); char *mMeasureTime=0; //! output file for measuring fitting time int mEnergySelect[3]; //! 0=MC (straight from dE), >0 see above @@ -207,10 +278,10 @@ class StFcsWaveformFitMaker : public StMaker { char* mFilename=0; char* mFilter=0; char mDetName[100]; - int mFitDrawOn=0; //! If set, it will also create a new TGraphAsymmErrors for each hit + int mFitDrawOn=0; //! If set to 1 it will create a new TGraphAsymmErrors for each hit, if set to 2 will create new TGraphAsmmErrors for each hit up to mMaxPage then reset virtual const Char_t *GetCVS() const {static const Char_t cvs[]="Tag " __DATE__ " " __TIME__ ; return cvs;} - ClassDef(StFcsWaveformFitMaker, 3) + ClassDef(StFcsWaveformFitMaker, 4) }; #endif // STROOT_STFCSWAVEFORMFITMAKER_STFCSWAVEFORMFITMAKER_H_