Update to StFcsWaveformFitMaker to use new algorithm

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];
-}
+}

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

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 = "<<name<<endm;
   mFcsDb = new StFcsDb("fcsDb");
   AddData(mFcsDb,".const");
+  mFcsDbPulse = new StFcsDbPulse("fcsPulse");
+  AddData(mFcsDbPulse,".const");
 }; 
 
 StFcsDbMaker::~StFcsDbMaker() {
-    //delete mFcsDb;   //already deleted by chain because AddData in constructor
+  //delete mFcsDb; //already deleted by chain because AddData in constructor
+  //delete mFcsDbPulse; //already deleted by chain because AddData in constructor
 }
 
 int StFcsDbMaker::Init(){
   mFcsDb->Init();
+  mFcsDbPulse->Init();
   return kStOK;
   return StMaker::Init();
 }

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

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(adc<mAdcSaturation) { HighY=Yerr; }
+  else                   { HighY=YerrSat; }
+  gae->SetPointError(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<npulse; i++){
+	ret += pulseShape(x, &p[2+i*3]);
+    }    
+    return ret;
+}
+
+TF1* StFcsDbPulse::createPulse(double xlow, double xhigh, int npars)
+{
+  if( npars<5 ){LOG_WARN << "StFcsDbPulse::createPulse - npars must be greater than or equal to 5 paramaters" <<endm; return 0;}
+  else if( xlow>=xhigh ){LOG_ERROR << "StFcsDbPulse::createPulse - Invalid range" <<endm; return 0;}
+    else{return new TF1("F1_fcsDbPulse",this,&StFcsDbPulse::multiPulseShape,xlow,xhigh,npars);}//More unique name in future??
+}
+
+int StFcsDbPulse::GenericPadPos(int value, int Nvals, int PadNums )
+{
+  if( value<=0 ){return ceil( static_cast<double>(value+(Nvals*PadNums))/static_cast<double>(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; i<gae->GetN(); ++i )
+    {
+      Double_t X; Double_t Y;
+      gae->GetPoint(i,X,Y);
+      if( X<xmin || X>xmax ){continue;}
+      if( Y>MaxY ){ MaxY=Y; index=i;}
+      if( Y<MinY ){ MinY=Y; }
+    }
+  if( index<0 ){std::cout << "Unable to find a maximum ADC value" << std::endl; return index;}
+  Ymin=MinY;
+  Ymax=MaxY;
+  return index;
+}
+
+void StFcsDbPulse::Print(Option_t* opt)
+{
+  std::cout << "Constants"
+	    << "|sqrtpi:"  << sqrtpi()
+	    << "|sqrt2pi:" << sqrt2pi()
+	    << std::endl;
+  std::cout << "TimebinInfo"
+	    << "|TBPerRC:" << TBPerRC()
+	    << "|nsecPerTB:"<< nsecPerTB()
+	    << "|BeamLengthSig:"<< BeamLengthSig()
+	    << std::endl;
+  std::cout << "PulseInfo"
+	    << "|GSigma:" << GSigma()
+	    << "|A1:"     << A1()
+	    << "|A2:"     << A2()
+	    << "|Xoff1:"  << Xoff1()
+	    << "|Xoff2:"  << Xoff2()
+	    << "|Tau1:"   << Tau1()
+	    << "|Tau2:"   << Tau2()
+	    << "|P1:"     << P1()
+	    << "|P2:"     << P2()
+	    << std::endl;
+}
+