Bug fix for StFcsWaveformFitMaker and some minor additions

StRoot/StFcsWaveformFitMaker/PeakAna.cxx

@@ -210,7 +210,7 @@ bool PeakAna::ValidPeakIdx() const
 bool PeakAna::GoodWindow()
 {
   if( mComputedIndex<0 ){this->AnalyzeForPeak();}
-  //First check if start and end times are within our max timebin window of 0-1023
+  //First check if peak values are within the specified range
   if( mFoundPeak.mStartX<mXRangeMin || mFoundPeak.mStartX>mXRangeMax ){return false;}
   else if( mFoundPeak.mEndX<mXRangeMin || mFoundPeak.mEndX > mXRangeMax ){return false;}
   //Next check if values make physical sense
@@ -230,8 +230,8 @@ void PeakAna::GetXYMax(Double_t xmin, Double_t xmax)
   }
   if( mMaxY<0 && mMaxX<0 ){
     LOG_WARN << "Unable to find a maximum adc\nSetting to impossible values" << endm;
-    mMaxY = 5000;
-    mMaxX = 5000;
+    mMaxY = mYRangeMax;
+    mMaxX = mXRangeMax;
   }
 }
 
@@ -730,61 +730,53 @@ void PeakAna::GetPossiblePeaks()
 //Returns index that corresponds to the serach criteria
 Int_t PeakAna::SearchForPeak(const std::vector<PeakWindow> &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();
+  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:"<<PossiblePeaks.size();
+      LOG_DEBUG << "|Search Peak:"<<mSearch.mStartX;
+      LOG_DEBUG << "|Search Width:"<<mSearch.mEndX;
+      LOG_DEBUG << endm;
     }
-  else
-    {
-      if( GetDebug()>0 )
-	{
-	  LOG_DEBUG << "|Size of Possible peaks:"<<PossiblePeaks.size();
-	  LOG_DEBUG << "|Search Peak:"<<mSearch.mStartX;
-	  LOG_DEBUG << "|Search Width:"<<mSearch.mEndX;
+    Short_t peakindex=-1;
+    for( UShort_t ipeak=0; ipeak<PossiblePeaks.size(); ++ipeak ){
+      if( GetDebug()>1 ){
+	LOG_DEBUG << " - ";
+	LOG_DEBUG << "|Index:"<<ipeak;
+	PossiblePeaks.at(ipeak).Print("debug");
+	LOG_DEBUG << endm;
+      }
+      Double_t PeakLoc = PossiblePeaks.at(ipeak).mPeakX;
+      if( mSearch.mStartX-mSearch.mEndX<=PeakLoc && PeakLoc <= mSearch.mStartX+mSearch.mEndX ){
+	peakindex=ipeak;
+	if( GetDebug()>1 ){
+	  LOG_DEBUG << "   + ";
+	  LOG_DEBUG << "|TrueIndex:"<<peakindex;
 	  LOG_DEBUG << endm;
 	}
-      Short_t peakindex=-1;
-      for( UShort_t ipeak=0; ipeak<PossiblePeaks.size(); ++ipeak )
-	{
-	  if( GetDebug()>1 )
-	    {
-	      LOG_DEBUG << " - ";
-	      LOG_DEBUG << "|Index:"<<ipeak;
-	      PossiblePeaks.at(ipeak).Print("debug");
-	      LOG_DEBUG << endm;
-	    }
-	  Double_t PeakLoc = PossiblePeaks.at(ipeak).mPeakX;
-	  if( mSearch.mStartX-mSearch.mEndX<=PeakLoc && PeakLoc <= mSearch.mStartX+mSearch.mEndX )
-	    {
-	      peakindex=ipeak;
-	      if( GetDebug()>1 )
-		{
-		  LOG_DEBUG << "   + ";
-		  LOG_DEBUG << "|TrueIndex:"<<peakindex;
-		  LOG_DEBUG << endm;
-		}
-	    }
-	}
-      if( peakindex>=0 && mSearch.mStartX>0 && PossiblePeaks.at(peakindex).mP_Peak>0 )
-	{
-	  if(GetDebug()>1){PossiblePeaks.at(peakindex).Print("debug");LOG_DEBUG<<endm;}
-	  return peakindex;
-	}
-      else{ return PossiblePeaks.size();}
+      }
+    }
+    if( peakindex>=0 && mSearch.mStartX>0 && PossiblePeaks.at(peakindex).mP_Peak>0 ){
+      if(GetDebug()>1){PossiblePeaks.at(peakindex).Print("debug");LOG_DEBUG<<endm;}
+      return peakindex;
     }
+    else{ return PossiblePeaks.size();}
+  }
 }
 
 Int_t PeakAna::SearchForPeak(const std::vector<PeakWindow> &PossiblePeaks, const PeakWindow& search)
 {
-	mSearch = search;
-	return this->SearchForPeak(PossiblePeaks);
+  mSearch = search;
+  return this->SearchForPeak(PossiblePeaks);
 }
 
 Int_t PeakAna::SearchForPeak(const std::vector<PeakWindow> &PossiblePeaks, Double_t peak, Double_t width)
 {
-	SetSearchWindow(peak, width);
-	return this->SearchForPeak(PossiblePeaks);
+  SetSearchWindow(peak, width);
+  return this->SearchForPeak(PossiblePeaks);
 }
 
 void PeakAna::MergeByProbability(std::vector<PeakWindow>& newpeaks) const
@@ -840,7 +832,6 @@ void PeakAna::MergeByChirality(std::vector<PeakWindow>& newpeaks) const
   //std::cout << "PeakAna::MergeByChirality - newpeaksize:" << newpeaks.size()  << std::endl;
 }
 
-
 short PeakAna::MergeLeftOrRight(UInt_t index) const
 {
   if( index>=mPeaks.size() ){return 0;}//invalid index
@@ -980,44 +971,10 @@ std::vector< std::pair<int,int> > PeakAna::MergeIndices(std::vector<short>& vec)
   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<double>(rx)*0.5; }
-  if( sy==0 ){ sy = static_cast<double>(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<nsizex; ++ix ){
-    for( int jy=0; jy<nsizey; ++jy ){
-      int xi = ix-rx;
-      int yi = jy-ry;
-      double product = 1.0;
-      if( sx!=0 ){ product *= exp( (-1.0*xi*xi)/(2.0*sx*sx) )/(sqrt(2.0*3.14159265358979323846)*sx); }
-      if( sy!=0 ){ product *= exp( (-1.0*yi*yi)/(2.0*sy*sy) )/(sqrt(2.0*3.14159265358979323846)*sy); }
-      GM_2D[ix*nsizey+jy] = product;
-      GM_sum += product;
-    }
-  }
-
-  if( kNorm ){
-    double invsum = 1.0/GM_sum;
-    for( int i=0; i<nsizex; ++i ){
-      for( int j=0; j<nsizey; ++j ){
-	GM_2D[i*nsizey+j] *= invsum;
-      }
-    }
-  }
-  return GM_2D;
-}
-
 void PeakAna::Draw(Option_t* opt)
 {
   AppendPad(opt);
+  //gPad->IncrementPaletteColor(1, mOption); (Doesn't work with ROOT v5) 
 }
 
 void PeakAna::Paint(Option_t* opt)
@@ -1278,3 +1235,38 @@ void PeakAna::SetAllPeakLineWidth(Width_t s_width, Width_t e_width)
   }
 }
 
+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<double>(rx)*0.5; }
+  if( sy==0 ){ sy = static_cast<double>(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<nsizex; ++ix ){
+    for( int jy=0; jy<nsizey; ++jy ){
+      int xi = ix-rx;
+      int yi = jy-ry;
+      double product = 1.0;
+      if( sx!=0 ){ product *= exp( (-1.0*xi*xi)/(2.0*sx*sx) )/(sqrt(2.0*3.14159265358979323846)*sx); }
+      if( sy!=0 ){ product *= exp( (-1.0*yi*yi)/(2.0*sy*sy) )/(sqrt(2.0*3.14159265358979323846)*sy); }
+      GM_2D[ix*nsizey+jy] = product;
+      GM_sum += product;
+    }
+  }
+
+  if( kNorm ){
+    double invsum = 1.0/GM_sum;
+    for( int i=0; i<nsizex; ++i ){
+      for( int j=0; j<nsizey; ++j ){
+	GM_2D[i*nsizey+j] *= invsum;
+      }
+    }
+  }
+  return GM_2D;
+}
+

StRoot/StFcsWaveformFitMaker/PeakAna.h

@@ -532,3 +532,4 @@ class PeakAna : public TObject
 };
 
 #endif
+

StRoot/StFcsWaveformFitMaker/PeakAnaPainter.h

@@ -1,5 +1,8 @@
 /*
 Author: David Kapukchyan
+@[October 20, 2022]
+> Added doxygen style comments
+
 @[September 29, 2022]
 > Got rid of the virtual painter
 

StRoot/StFcsWaveformFitMaker/PeakWindow.cxx

@@ -334,6 +334,7 @@ TLine* PeakWindow::GetStartLine(Double_t ymin, Double_t ymax)
   }
   return mStartLine;
 }
+
 Color_t PeakWindow::GetStartLineColor() const
 { if( mStartLine!=0 ){return mStartLine->GetLineColor(); } return 0; }
 Style_t PeakWindow::GetStartLineStyle() const
@@ -357,6 +358,7 @@ TMarker* PeakWindow::GetPeakMarker()
   }
   return mPeakMarker;
 }
+
 Color_t PeakWindow::GetPeakMarkerColor() const
 { if( mPeakMarker!=0 ){return mPeakMarker->GetMarkerColor(); } return 0; }
 Style_t PeakWindow::GetPeakMarkerStyle() const
@@ -388,6 +390,7 @@ TLine* PeakWindow::GetEndLine(Double_t ymin, Double_t ymax)
   }
   return mEndLine;
 }
+
 Color_t PeakWindow::GetEndLineColor() const
 { if( mEndLine!=0 ){ return mEndLine->GetLineColor(); } return 0; }
 Style_t PeakWindow::GetEndLineStyle() const

StRoot/StFcsWaveformFitMaker/PeakWindow.h

@@ -1,5 +1,8 @@
 /*
 Author: David Kapukchyan
+@[October 20, 2022]
+> Added doxygen style comments
+
 @[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.
 

StRoot/StFcsWaveformFitMaker/StFcsPulseAna.cxx

@@ -620,17 +620,22 @@ Double_t StFcsPulseAna::PulseFit(Int_t Start, Int_t End)
   else{ mFitSum = true; return mSumAdcFit; }//If fit failed return 0 don't do a histogram sum
 }
 
-void StFcsPulseAna::SetFitPars(TF1* func)
+void StFcsPulseAna::SetFitPars(TF1*& func, int start, int end)
 {
-  if( mDbPulse==0 ){
-    if( func!=0 ){
-      func->FixParameter(0,0);
-      func->FixParameter(1,0);
-    }
-    return;
+  if( func!=0 ){
+    func->FixParameter(0,0);
+    func->FixParameter(1,0);
   }
+  if( mDbPulse==0 ){ return; }
   if( FoundPeakIndex()<0 ){ this->AnalyzeForPeak(); }
   int npeaks = NPeaks();
+  if( start>=npeaks ){ return; }
+  if( end>=npeaks ){ return; }
+  if( start<0 ){ start = 0; }
+  if( end<0 ){ end = npeaks-1; }
+  if( start>end ){ return; }
+
+  npeaks = end-start+1;  //Number of peaks to be used in function is end-start+1 since including end
   if( func==0 ){
     func = mDbPulse->createPulse(mXRangeMin,mXRangeMax,2+npeaks*3);
   }
@@ -640,25 +645,26 @@ void StFcsPulseAna::SetFitPars(TF1* func)
   para[1] = Baseline();
   func->SetParName(0,"NPeaks");
   func->SetParName(1,"Ped");
-  for( Int_t i=0; i<npeaks; ++i ){
+  for( int i=0; i<npeaks; ++i ){
     int j = 1+i*3+1;
-    char name[10];
-    sprintf(name,"P%d_A",i);
+    char name[40];
+    sprintf(name,"P%d_A",start+i);
     func->SetParName(j,name);
-    para[j++] = GetPeak(i).mPeakY;
-    sprintf(name,"P%d_M",i);
+    para[j++] = GetPeak(start+i).mPeakY;
+    sprintf(name,"P%d_M",start+i);
     func->SetParName(j,name);
-    para[j++] = GetPeak(i).mPeakX;
-    sprintf(name,"P%d_S",i);
+    para[j++] = GetPeak(start+i).mPeakX;
+    sprintf(name,"P%d_S",start+i);
     func->SetParName(j,name);
     para[j] = mDbPulse->GSigma();
   }
   func->SetParameters(para);
+
   func->FixParameter(0,npeaks);
   func->FixParameter(1,Baseline());
   for(int i=0; i<npeaks; i++){
     int j=1+i*3+1;
-    func->SetParLimits(j++,0.0,50000.0);       //limit peak not to go negative
+    func->SetParLimits(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
   }

StRoot/StFcsWaveformFitMaker/StFcsPulseAna.h

@@ -1,5 +1,8 @@
 //Author:David Kapukchyan
 /*
+@[December 17, 2022]
+> SetFitPars can now set a function parameters for a range of in npeaks as opposed to all peaks
+
 @[June 25, 2022]
  > Added Refinements to drawing and copying. Also made some changes to how the tunnel probability is used to merge the peaks.
 
@@ -8,7 +11,7 @@
 */
 
 /*!
-  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.
+  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 preshower.
 */
 
 #ifndef STFCSPULSEANA_H
@@ -71,7 +74,7 @@ class StFcsPulseAna : public PeakAna
   Double_t MBFit(Int_t Start=0,Int_t End=0);        //!< Fit a Maxwell-Boltzmann distribution to #mFoundPeak and return the integral minus the baseline
   Double_t PulseFit(Int_t Start=0, Int_t End=0);    //!< Fit the pulse shape defined in #StFcsDbPulse::multiPulseShape() to all peaks and return the integral of the found peak minus the baseline
 
-  void SetFitPars(TF1* func);                       //!< Set the parameters of an external TF1 function that has the form of #StFcsDbPulse::multiPulseShape()
+  void SetFitPars(TF1*& func, int start=-1, int end=-1);   //!< Set the parameters of an external TF1 function that has the form of #StFcsDbPulse::multiPulseShape(), optionally only set fit paramaters for peaks from index start up to and including end
 
   static void FillAdc(TGraphAsymmErrors* g, unsigned short& counter, int Start, unsigned short* adcdata);//!< Needed for SumDep0, et al. It basically fills in tb vs. adc sequentially so all timebins from a given start value have an adc.
   static int SumDep0(TGraphAsymmErrors* gdata, int Start, int ped=0);     //!< Test of sum method in DEP board
@@ -80,7 +83,7 @@ class StFcsPulseAna : public PeakAna
   void SetFitSignal(TF1* func){mF1_SignalFit=func;}      //!< @param func sets #mF1_SignalFit
   void SetBaselineFit(TF1* func){mF1_BaselineFit=func;}  //!< @param func sets #mF1_BaselineFit
 
-  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 SetZS(){SetBaseline(0,0.39);} //!< 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