Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:12:21

0001 // $Id: $
0002 
0003 /*!
0004  * \file PROTOTYPE4_FEM.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010 
0011 #include "PROTOTYPE4_FEM.h"
0012 
0013 #include <TAttMarker.h>  // for kFullCircle
0014 #include <TCanvas.h>
0015 #include <TF1.h>
0016 #include <TGraph.h>
0017 
0018 #include <cassert>
0019 #include <cmath>
0020 #include <cstdlib>  // for exit
0021 #include <iostream>
0022 #include <limits>
0023 #include <memory>  // for allocator_traits<>::value_type
0024 #include <string>
0025 
0026 
0027 using namespace std;
0028 
0029 int PROTOTYPE4_FEM::GetChannelNumber(const std::string &caloname, int i_column,
0030                                      int i_row)
0031 {
0032   if (caloname == "HCALIN")
0033   {
0034     assert(i_row >= 0);
0035     assert(i_row < NCH_IHCAL_ROWS);
0036     assert(i_column >= 0);
0037     assert(i_column < NCH_IHCAL_COLUMNS);
0038 
0039     /*
0040     static const int hbdchanIHC_col_0[4] = { 4, 3, 2, 1}; // i_row = 0, 1, 2, 3
0041     static const int hbdchanIHC_col_1[4] = { 8, 7, 6, 5};
0042     static const int hbdchanIHC_col_2[4] = {12,11,10, 9};
0043     static const int hbdchanIHC_col_3[4] = {16,15,14,13};
0044 
0045     if(i_column == 0) return hbdchanIHC_col_0[i_row] + 64 - 1;
0046     if(i_column == 1) return hbdchanIHC_col_1[i_row] + 64 - 1;
0047     if(i_column == 2) return hbdchanIHC_col_2[i_row] + 64 - 1;
0048     if(i_column == 3) return hbdchanIHC_col_3[i_row] + 64 - 1;
0049     */
0050 
0051     static const int hbdchanIHC_col_0[4] = {67, 66, 65,
0052                                             64};  // i_row = 0, 1, 2, 3
0053     static const int hbdchanIHC_col_1[4] = {71, 70, 69, 68};
0054     static const int hbdchanIHC_col_2[4] = {75, 74, 73, 72};
0055     static const int hbdchanIHC_col_3[4] = {79, 78, 77, 76};
0056 
0057     if (i_column == 0)
0058       return hbdchanIHC_col_0[i_row];
0059     if (i_column == 1)
0060       return hbdchanIHC_col_1[i_row];
0061     if (i_column == 2)
0062       return hbdchanIHC_col_2[i_row];
0063     if (i_column == 3)
0064       return hbdchanIHC_col_3[i_row];
0065   }
0066   else if (caloname == "HCALOUT")
0067   {
0068     assert(i_row >= 0);
0069     assert(i_row < NCH_OHCAL_ROWS);
0070     assert(i_column >= 0);
0071     assert(i_column < NCH_OHCAL_COLUMNS);
0072 
0073     static const int hbdchanOHC_col_0[4] = {116, 118, 112,
0074                                             114};  // i_row = 0, 1, 2, 3
0075     static const int hbdchanOHC_col_1[4] = {124, 126, 120, 122};
0076     static const int hbdchanOHC_col_2[4] = {132, 134, 128, 130};
0077     static const int hbdchanOHC_col_3[4] = {140, 142, 136, 138};
0078 
0079     if (i_column == 0)
0080       return hbdchanOHC_col_0[i_row];
0081     if (i_column == 1)
0082       return hbdchanOHC_col_1[i_row];
0083     if (i_column == 2)
0084       return hbdchanOHC_col_2[i_row];
0085     if (i_column == 3)
0086       return hbdchanOHC_col_3[i_row];
0087   }
0088   else if (caloname == "EMCAL")
0089   {
0090     assert(i_row >= 0);
0091     assert(i_row < NCH_EMCAL_ROWS);
0092     assert(i_column >= 0);
0093     assert(i_column < NCH_EMCAL_COLUMNS);
0094 
0095     //    > Anthony Hodges
0096     //    > PhD. Student, Georgia State University
0097     //    > Nuclear and High Energy Physics
0098     //    > ahodges21@student.gsu.edu
0099 
0100     // mapping taken from John Haggerty's emcalall.C found here:
0101     // /gpfs/mnt/gpfs02/sphenix/data/data01/caladc/wd/wd409/macros
0102     // This map here takes in a channel number and gives you the corresponding
0103     // canvas position Presumably we want the opposite, put in canvas position,
0104     // output channel number So we'll work backwards
0105     //    Float_t canmap[64] = {6, 7, 14, 15, 4, 5, 12, 13, 2, 3, 10, 11, 0, 1,
0106     //    8, 9, 22, 23, 30, 31, 20, 21, 28, 29,
0107     //                          18, 19, 26, 27, 16, 17, 24, 25, 38, 39, 46, 47,
0108     //                          36, 37, 44, 45, 34, 35, 42, 43, 32, 33, 40, 41,
0109     //                          54, 55, 62, 63, 52, 53, 60, 61, 50, 51, 58, 59,
0110     //                          48, 49, 56, 57};
0111     //    static int hbdchanEMC[8][8];  //I'm gonna fill this boy with the above
0112     //    channels for (int chan = 0; chan < 64; chan++)
0113     //    {
0114     //      hbdchanEMC[(int) floor(canmap[chan] / 8)][((int) canmap[chan]) % 8]
0115     //      = chan;
0116     //    }
0117 
0118     // Revision from Martin with static reverse:
0119     //    This is now lining up towers from 0....63, and tells you
0120     //
0121     //    for tower i, what is the actual ADC index I have to go to?
0122     //
0123     //    tower[0] -> chvector[0] = 12  is then the adc channel nr.
0124     //
0125     //
0126     //
0127     //
0128     //
0129     //
0130     //     static const int  chvector[]=  {12 , 13 , 8 , 9 , 4 ,5 ,0 ,1 ,14 ,15
0131     //    ,10 ,11 ,6 ,7 ,2 ,3 ,28 ,29 ,24 ,25 ,20 ,21 ,
0132     //                                      16 ,17 ,30 ,31 ,26 ,27 ,22 ,23 ,18
0133     //                                      ,19
0134     //    ,44 ,45 ,40 ,41 ,36 ,37 ,32 ,33 ,46 ,47 ,42 ,43 ,38 ,
0135     //                                      39 ,34 ,35 ,60 ,61 ,56 ,57 ,52 ,53
0136     //                                      ,48
0137     //    ,49 ,62 ,63 ,58 ,59 ,54 ,55 ,50 ,51};
0138     const static int canmap[64] = {
0139         12, 13, 8, 9, 4, 5, 0, 1, 14, 15, 10, 11, 6, 7, 2, 3,
0140         28, 29, 24, 25, 20, 21, 16, 17, 30, 31, 26, 27, 22, 23, 18, 19,
0141         44, 45, 40, 41, 36, 37, 32, 33, 46, 47, 42, 43, 38, 39, 34, 35,
0142         60, 61, 56, 57, 52, 53, 48, 49, 62, 63, 58, 59, 54, 55, 50, 51};
0143 
0144     const int linear_channel_ID =
0145         (NCH_EMCAL_ROWS - i_row - 1) * NCH_EMCAL_COLUMNS + i_column;
0146 
0147     assert(linear_channel_ID >= 0);
0148     assert(linear_channel_ID < 64);
0149 
0150     return canmap[linear_channel_ID];
0151   }
0152 
0153   std::cout << "PROTOTYPE4_FEM::GetHBDCh - invalid input caloname " << caloname
0154             << " i_column " << i_column << " i_row " << i_row << std::endl;
0155   exit(1);
0156   return -9999;
0157 }
0158 
0159 bool PROTOTYPE4_FEM::SampleFit_PowerLawExp(  //
0160     const std::vector<double> &samples,      //
0161     double &peak,                            //
0162     double &peak_sample,                     //
0163     double &pedestal,                        //
0164     const int verbosity)
0165 {
0166   int peakPos = 0.;
0167 
0168   assert(samples.size() == NSAMPLES);
0169 
0170   TGraph gpulse(NSAMPLES);
0171   for (int i = 0; i < NSAMPLES; i++)
0172   {
0173     (gpulse.GetX())[i] = i;
0174 
0175     (gpulse.GetY())[i] = samples[i];
0176   }
0177 
0178   pedestal = gpulse.GetY()[0];  //(double) PEDESTAL;
0179   double peakval = pedestal;
0180   const double risetime = 4;
0181 
0182   for (int iSample = 0; iSample < NSAMPLES; iSample++)
0183   {
0184     if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
0185     {
0186       peakval = gpulse.GetY()[iSample];
0187       peakPos = iSample;
0188     }
0189   }
0190   peakval -= pedestal;
0191 
0192   // fit function
0193   TF1 fits("f_SignalShape_PowerLawExp", SignalShape_PowerLawExp, 0., NSAMPLES,
0194            5);
0195 
0196   double par[6] = {0};
0197   par[0] = peakval;  // /3.;
0198   par[1] = peakPos - risetime;
0199   if (par[1] < 0.)
0200     par[1] = 0.;
0201   par[2] = 4.;
0202   par[3] = 1.5;
0203   par[4] = pedestal;
0204   par[5] = 0;
0205   fits.SetParameters(par);
0206   fits.SetParNames("Amplitude", "Sample Start", "Power", "Decay", "pedestal",
0207                    "Baseline shift");
0208   fits.SetParLimits(0, peakval * 0.5, peakval * 10);
0209   fits.SetParLimits(1, 0, NSAMPLES);
0210   fits.SetParLimits(2, 0, 10.);
0211   fits.SetParLimits(3, 0, 10);
0212   fits.SetParLimits(4, pedestal - abs(peakval), pedestal + abs(peakval));
0213   //  fits.SetParLimits(5, - abs(peakval),  + abs(peakval));
0214   //  fits.FixParameter(5, 0);
0215 
0216   // Saturation correction - Abhisek
0217   for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
0218     if ((gpulse.GetY())[ipoint] <= 10 or
0219         (gpulse.GetY())[ipoint] >=
0220             ((1 << 14) - 10))  // drop point if touching max or low limit on ADCs
0221     {
0222       gpulse.RemovePoint(ipoint);
0223       ipoint--;
0224     }
0225 
0226   if (verbosity <= 1)
0227     gpulse.Fit(&fits, "MQRN0W", "goff", 0., (double) NSAMPLES);
0228   else
0229     gpulse.Fit(&fits, "MRN0VW", "goff", 0., (double) NSAMPLES);
0230 
0231   if (verbosity)
0232   {
0233     static int id = 0;
0234     ++id;
0235 
0236     string c_name(string("PROTOTYPE4_FEM_SampleFit_PowerLawExp_") +
0237                   to_string(id));
0238 
0239     TCanvas *canvas = new TCanvas(c_name.c_str(), c_name.c_str());
0240     canvas->Update();
0241 
0242     TGraph *g_plot = static_cast<TGraph *>(gpulse.DrawClone("ap*l"));
0243     g_plot->SetTitle("ADC data and fit;Sample number;ADC value");
0244     fits.DrawClone("same");
0245     fits.Print();
0246     canvas->Update();
0247     //    sleep(1);
0248   }
0249 
0250   //  peak = fits.GetParameter(0); // not exactly peak height
0251   peak =
0252       (fits.GetParameter(0) *
0253        pow(fits.GetParameter(2) / fits.GetParameter(3), fits.GetParameter(2))) /
0254       exp(fits.GetParameter(
0255           2));  // exact peak height is (p0*Power(p2/p3,p2))/Power(E,p2)
0256 
0257   //  peak_sample = fits.GetParameter(1); // signal start time
0258   peak_sample = fits.GetParameter(1) +
0259                 fits.GetParameter(2) / fits.GetParameter(3);  // signal peak time
0260 
0261   // peak integral = p0*Power(p3,-1 - p2)*Gamma(1 + p2). Note yet used in output
0262 
0263   pedestal = fits.GetParameter(4);
0264 
0265   return true;
0266 }
0267 
0268 double PROTOTYPE4_FEM::SignalShape_PowerLawExp(double *x, double *par)
0269 {
0270   double pedestal = par[4];
0271   //                        + ((x[0] - 1.5 * par[1]) > 0) * par[5];  // quick
0272   //                        fix on exting tails on the signal function
0273   if (x[0] < par[1])
0274     return pedestal;
0275   // double  signal =
0276   // (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0277   double signal =
0278       par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
0279   return pedestal + signal;
0280 }
0281 
0282 bool PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(  //
0283     const std::vector<double> &samples,            //
0284     double &peak,                                  //
0285     double &peak_sample,                           //
0286     double &pedestal,                              //
0287     std::map<int, double> &parameters_io, const int verbosity)
0288 {
0289   static const int n_parameter = 7;
0290 
0291   // inital guesses
0292   int peakPos = 0.;
0293 
0294   assert(samples.size() == NSAMPLES);
0295 
0296   TGraph gpulse(NSAMPLES);
0297   for (int i = 0; i < NSAMPLES; i++)
0298   {
0299     (gpulse.GetX())[i] = i;
0300 
0301     (gpulse.GetY())[i] = samples[i];
0302   }
0303 
0304   // Saturation correction - Abhisek
0305   for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
0306     if ((gpulse.GetY())[ipoint] <= 10 or
0307         (gpulse.GetY())[ipoint] >=
0308             ((1 << 14) - 10)  // drop point if touching max or low limit on ADCs
0309         or (not isnormal((gpulse.GetY())[ipoint])))
0310     {
0311       gpulse.RemovePoint(ipoint);
0312       ipoint--;
0313     }
0314 
0315   pedestal = gpulse.GetY()[0];  //(double) PEDESTAL;
0316   double peakval = pedestal;
0317   const double risetime = 2;
0318 
0319   for (int iSample = 0; iSample < NSAMPLES - risetime * 3; iSample++)
0320   {
0321     if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
0322     {
0323       peakval = gpulse.GetY()[iSample];
0324       peakPos = iSample;
0325     }
0326   }
0327   peakval -= pedestal;
0328 
0329   if (verbosity)
0330   {
0331     cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - "
0332          << "pedestal = " << pedestal << ", "
0333          << "peakval = " << peakval << ", "
0334          << "peakPos = " << peakPos << endl;
0335   }
0336 
0337   // build default value
0338   struct default_values_t
0339   {
0340     default_values_t(double default_value, double min_value, double max_value)
0341       : def(default_value)
0342       , min(min_value)
0343       , max(max_value)
0344     {
0345     }
0346     double def;
0347     double min;
0348     double max;
0349   };
0350 
0351   vector<default_values_t> default_values(
0352       n_parameter, default_values_t(numeric_limits<double>::signaling_NaN(),
0353                                     numeric_limits<double>::signaling_NaN(),
0354                                     numeric_limits<double>::signaling_NaN()));
0355 
0356   default_values[0] =
0357       default_values_t(peakval * .7, peakval * -1.5, peakval * 1.5);
0358   default_values[1] = default_values_t(
0359       peakPos - risetime, peakPos - 3 * risetime, peakPos + risetime);
0360   default_values[2] = default_values_t(2., 1, 5.);
0361   default_values[3] = default_values_t(5, risetime * .5, risetime * 4);
0362   default_values[4] = default_values_t(pedestal, pedestal - abs(peakval),
0363                                        pedestal + abs(peakval));
0364   default_values[5] = default_values_t(.3, 0, 1);
0365   default_values[6] = default_values_t(5, risetime * .5, risetime * 4);
0366 
0367   // fit function
0368   TF1 fits("f_SignalShape_PowerLawDoubleExp", SignalShape_PowerLawDoubleExp, 0.,
0369            NSAMPLES, n_parameter);
0370   fits.SetParNames("Amplitude", "Sample Start", "Power", "Peak Time 1",
0371                    "Pedestal", "Amplitude ratio", "Peak Time 2");
0372 
0373   for (int i = 0; i < n_parameter; ++i)
0374   {
0375     if (parameters_io.find(i) == parameters_io.end())
0376     {
0377       fits.SetParameter(i, default_values[i].def);
0378       fits.SetParLimits(i, default_values[i].min, default_values[i].max);
0379 
0380       if (verbosity)
0381       {
0382         cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - parameter [" << i
0383              << "]: "
0384              << "default value = " << default_values[i].def
0385              << ", min value = " << default_values[i].min
0386              << ", max value = " << default_values[i].max << endl;
0387       }
0388     }
0389     else
0390     {
0391       fits.SetParLimits(i, parameters_io[i], parameters_io[i]);
0392       fits.SetParameter(i, parameters_io[i]);
0393 
0394       if (verbosity)
0395       {
0396         cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - parameter [" << i
0397              << "]: fixed to " << parameters_io[i] << endl;
0398       }
0399     }
0400   }
0401 
0402   if (verbosity <= 1)
0403     gpulse.Fit(&fits, "QRN0W", "goff", 0., (double) NSAMPLES);
0404   else
0405     gpulse.Fit(&fits, "RN0VW", "goff", 0., (double) NSAMPLES);
0406 
0407   if (verbosity)
0408   {
0409     static int id = 0;
0410     ++id;
0411 
0412     string c_name(string("PROTOTYPE4_FEM_SampleFit_PowerLawDoubleExp_") +
0413                   to_string(id));
0414 
0415     TCanvas *canvas = new TCanvas(c_name.c_str(), c_name.c_str());
0416     canvas->Update();
0417 
0418     TGraph *g_plot = static_cast<TGraph *>(gpulse.DrawClone("ap*l"));
0419     g_plot->SetTitle("ADC data and fit;Sample number;ADC value");
0420 
0421     fits.SetLineColor(kMagenta);
0422     fits.DrawClone("same");
0423     fits.Print();
0424 
0425     TF1 f1("f_SignalShape_PowerLawExp1", SignalShape_PowerLawExp, 0., NSAMPLES,
0426            5);
0427     f1.SetParameters(fits.GetParameter(0) * (1 - fits.GetParameter(5)) /
0428                          pow(fits.GetParameter(3), fits.GetParameter(2)) *
0429                          exp(fits.GetParameter(2)),
0430                      fits.GetParameter(1), fits.GetParameter(2),
0431                      fits.GetParameter(2) / fits.GetParameter(3),
0432                      fits.GetParameter(4));
0433     f1.SetLineColor(kBlue);
0434     f1.DrawClone("same");
0435 
0436     TF1 f2("f_SignalShape_PowerLawExp2", SignalShape_PowerLawExp, 0., NSAMPLES,
0437            5);
0438     f2.SetParameters(fits.GetParameter(0) * fits.GetParameter(5) /
0439                          pow(fits.GetParameter(6), fits.GetParameter(2)) *
0440                          exp(fits.GetParameter(2)),
0441                      fits.GetParameter(1), fits.GetParameter(2),
0442                      fits.GetParameter(2) / fits.GetParameter(6),
0443                      fits.GetParameter(4));
0444     f2.SetLineColor(kRed);
0445     f2.DrawClone("same");
0446 
0447     canvas->Update();
0448   }
0449 
0450   // store results
0451   pedestal = fits.GetParameter(4);
0452 
0453   const double peakpos1 = fits.GetParameter(3);
0454   const double peakpos2 = fits.GetParameter(6);
0455   double max_peakpos =
0456       fits.GetParameter(1) + (peakpos1 > peakpos2 ? peakpos1 : peakpos2);
0457   if (max_peakpos > NSAMPLES - 1)
0458     max_peakpos = NSAMPLES - 1;
0459 
0460   if (fits.GetParameter(0) > 0)
0461     peak_sample = fits.GetMaximumX(fits.GetParameter(1), max_peakpos);
0462   else
0463     peak_sample = fits.GetMinimumX(fits.GetParameter(1), max_peakpos);
0464 
0465   peak = fits.Eval(peak_sample) - pedestal;
0466 
0467   if (verbosity)
0468   {
0469     TGraph g_max(NSAMPLES);
0470 
0471     g_max.GetX()[0] = peak_sample;
0472     g_max.GetY()[0] = peak + pedestal;
0473 
0474     g_max.SetMarkerStyle(kFullCircle);
0475     g_max.SetMarkerSize(2);
0476     g_max.SetMarkerColor(kRed);
0477 
0478     g_max.DrawClone("p");
0479   }
0480 
0481   for (int i = 0; i < n_parameter; ++i)
0482   {
0483     parameters_io[i] = fits.GetParameter(i);
0484   }
0485 
0486   if (verbosity)
0487   {
0488     cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - "
0489          << "peak_sample = " << peak_sample << ", "
0490          << "max_peakpos = " << max_peakpos << ", "
0491          << "fits.GetParameter(1) = " << fits.GetParameter(1) << ", "
0492          << "peak = " << peak << ", "
0493          << "pedestal = " << pedestal << endl;
0494   }
0495 
0496   return true;
0497 }
0498 
0499 double PROTOTYPE4_FEM::SignalShape_PowerLawDoubleExp(double *x, double *par)
0500 {
0501   double pedestal = par[4];
0502   //                        + ((x[0] - 1.5 * par[1]) > 0) * par[5];  // quick
0503   //                        fix on exting tails on the signal function
0504   if (x[0] < par[1])
0505     return pedestal;
0506   // double  signal =
0507   // (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0508   //  peak / pow(fits.GetParameter(2) / fits.GetParameter(3),
0509   //  fits.GetParameter(2)) * exp(fits.GetParameter(2)) = fits.GetParameter(0);
0510   //  // exact peak height is (p0*Power(p2/p3,p2))/Power(E,p2)
0511   //  fits.GetParameter(2) / peak_shift =  fits.GetParameter(3);  // signal peak
0512   //  time
0513 
0514   double signal =                     //
0515       par[0]                          //
0516       * pow((x[0] - par[1]), par[2])  //
0517       * (((1. - par[5]) / pow(par[3], par[2]) * exp(par[2])) *
0518              exp(-(x[0] - par[1]) * (par[2] / par[3]))  //
0519          + (par[5] / pow(par[6], par[2]) * exp(par[2])) *
0520                exp(-(x[0] - par[1]) * (par[2] / par[6]))  //
0521         );
0522   return pedestal + signal;
0523 }
0524 
0525 bool PROTOTYPE4_FEM::SampleFit_PeakSample(  //
0526     const std::vector<double> &samples,     //
0527     double &peak,                           //
0528     double &peak_sample,                    //
0529     double &pedestal,                       //
0530     const int verbosity)
0531 {
0532   int peakPos = 0.;
0533 
0534   assert(samples.size() == NSAMPLES);
0535 
0536   const static int N_pedestal = 3;
0537   pedestal = 0;
0538   for (int iSample = 0; iSample < N_pedestal; iSample++)
0539   {
0540     pedestal += samples[iSample];
0541   }
0542   pedestal /= N_pedestal;
0543 
0544   double peakval = pedestal;
0545   for (int iSample = N_pedestal; iSample < NSAMPLES; iSample++)
0546   {
0547     if (abs(samples[iSample] - pedestal) > abs(peakval - pedestal))
0548     {
0549       peakval = samples[iSample];
0550       peakPos = iSample;
0551     }
0552   }
0553 
0554   peak = peakval - pedestal;
0555   peak_sample = peakPos;
0556 
0557   if (verbosity)
0558   {
0559     cout << "PROTOTYPE4_FEM::SampleFit_PeakSample - "
0560          << "pedestal = " << pedestal << ", "
0561          << "peakval = " << peakval << ", "
0562          << "peakPos = " << peakPos << endl;
0563   }
0564 
0565   return true;
0566 }