packages/Prototype3/PROTOTYPE3_FEM.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file PROTOTYPE3_FEM.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "PROTOTYPE3_FEM.h"
0012
0013 #include <TCanvas.h>
0014 #include <TF1.h>
0015 #include <TGraph.h>
0016
0017 #include <unistd.h>  // for sleep
0018 #include <cassert>
0019 #include <cmath>
0020 #include <cstdlib>  // for exit
0021 #include <iostream>
0022
0023 using namespace std;
0024
0025 int PROTOTYPE3_FEM::GetHBDCh(const std::string &caloname, int i_column,
0026                              int i_row)
0027 {
0028   if (caloname == "HCALIN")
0029   {
0030     return 64 + 8 * i_column + 2 * i_row;
0031   }
0032   else if (caloname == "HCALOUT")
0033   {
0034     return 112 + 8 * i_column + 2 * i_row;
0035   }
0036   else if (caloname == "EMCAL_PROTOTYPE2")
0037   {
0038     // EMcal cable mapping from John haggerty
0039     assert(i_column >= 0);
0040     assert(i_column < NCH_EMCAL_COLUMNS);
0041     assert(i_row >= 0);
0042     assert(i_row < NCH_EMCAL_ROWS);
0043
0044     static int canmap[NCH_EMCAL_ROWS * NCH_EMCAL_COLUMNS] = {
0045         // 1 ... 15
0046         10 + 48, 11 + 48, 8 + 48, 9 + 48, 14 + 48, 15 + 48, 12 + 48, 13 + 48,
0047         // 9 ... 16
0048         2 + 48, 3 + 48, 0 + 48, 1 + 48, 6 + 48, 7 + 48, 4 + 48, 5 + 48,
0049
0050         // 17 ... 24
0051         10 + 32, 11 + 32, 8 + 32, 9 + 32, 14 + 32, 15 + 32, 12 + 32, 13 + 32,
0052         // 25 ... 32
0053         2 + 32, 3 + 32, 0 + 32, 1 + 32, 6 + 32, 7 + 32, 4 + 32, 5 + 32,
0054
0055         // 33 ... 40
0056         10 + 16, 11 + 16, 8 + 16, 9 + 16, 14 + 16, 15 + 16, 12 + 16, 13 + 16,
0057         // 41 42 43 44 45 46 47 48
0058         2 + 16, 3 + 16, 0 + 16, 1 + 16, 6 + 16, 7 + 16, 4 + 16, 5 + 16,
0059
0060         // 49 50 51 52 53 54 55 56
0061         10, 11, 8, 9, 14, 15, 12, 13,
0062         // 57 58 59 60 61 62 63 64
0063         2, 3, 0, 1, 6, 7, 4, 5};
0064
0065     const int tower_index =
0066         i_column + NCH_EMCAL_COLUMNS * (NCH_EMCAL_ROWS - 1 - i_row);
0067
0068     assert(tower_index >= 0);
0069     assert(tower_index < NCH_EMCAL_ROWS * NCH_EMCAL_COLUMNS);
0070
0071     return canmap[tower_index];
0072   }
0073   else if (caloname == "EMCAL_HIGHETA")
0074   {
0075     // EMcal cable mapping from John haggerty
0076     assert(i_column >= 0);
0077     assert(i_column < NCH_EMCAL_COLUMNS);
0078     assert(i_row >= 0);
0079     assert(i_row < NCH_EMCAL_ROWS);
0080
0081     static int canmap[] = {
0082         //            >
0083         //            https://docdb.sphenix.bnl.gov/0000/000034/001/T1044-2017a-2.xlsx
0084         //            Sean and John spot checked a number of these towers at
0085         //            BNL.  There could be mistakes, but cosmics look
0086         //            reasonable.
0087         //            Front view (looking downstream, same as above) but in HBD
0088         //            channel numbers
0089         //            3  , 2  , 19 , 18 , 35 , 34 , 51 , 50,
0090         //            1  , 0  , 17 , 16 , 33 , 32 , 49 , 48,
0091         //            7  , 6  , 23 , 22 , 39 , 38 , 55 , 54,
0092         //            5  , 4  , 21 , 20 , 37 , 36 , 53 , 52,
0093         //            11 , 10 , 27 , 26 , 43 , 42 , 59 , 58,
0094         //            9  , 8  , 25 , 24 , 41 , 40 , 57 , 56,
0095         //            15 , 14 , 31 , 30 , 47 , 46 , 63 , 62,
0096         //            13 , 12 , 29 , 28 , 45 , 44 , 61 , 60,
0097         3, 2, 19, 18, 35, 34, 51, 50, 1, 0, 17, 16, 33, 32, 49, 48, 7,
0098         6, 23, 22, 39, 38, 55, 54, 5, 4, 21, 20, 37, 36, 53, 52, 11, 10,
0099         27, 26, 43, 42, 59, 58, 9, 8, 25, 24, 41, 40, 57, 56, 15, 14, 31,
0100         30, 47, 46, 63, 62, 13, 12, 29, 28, 45, 44, 61, 60, 0};
0101
0102     const int tower_index =
0103         i_column + NCH_EMCAL_COLUMNS * (NCH_EMCAL_ROWS - 1 - i_row);
0104
0105     assert(tower_index >= 0);
0106     assert(tower_index < NCH_EMCAL_ROWS * NCH_EMCAL_COLUMNS);
0107
0108     return canmap[tower_index];
0109   }
0110
0111   std::cout << "PROTOTYPE3_FEM::GetHBDCh - invalid input caloname " << caloname
0112             << " i_column " << i_column << " i_row " << i_row << std::endl;
0113   exit(1);
0114   return -9999;
0115 }
0116
0117 bool PROTOTYPE3_FEM::SampleFit_PowerLawExp(  //
0118     const std::vector<double> &samples,      //
0119     double &peak,                            //
0120     double &peak_sample,                     //
0121     double &pedstal,                         //
0122     const int verbosity)
0123 {
0124   int peakPos = 0.;
0125
0126   assert(samples.size() == NSAMPLES);
0127
0128   TGraph gpulse(NSAMPLES);
0129   for (int i = 0; i < NSAMPLES; i++)
0130   {
0131     (gpulse.GetX())[i] = i;
0132
0133     (gpulse.GetY())[i] = samples[i];
0134   }
0135
0136   double pedestal = gpulse.GetY()[0];  //(double) PEDESTAL;
0137   double peakval = pedestal;
0138   const double risetime = 4;
0139
0140   for (int iSample = 0; iSample < NSAMPLES; iSample++)
0141   {
0142     if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
0143     {
0144       peakval = gpulse.GetY()[iSample];
0145       peakPos = iSample;
0146     }
0147   }
0148   peakval -= pedestal;
0149
0150   // fit function
0151   TF1 fits("f_SignalShape_PowerLawExp", SignalShape_PowerLawExp, 0., 24., 6);
0152
0153   double par[6] = {0};
0154   par[0] = peakval;  // /3.;
0155   par[1] = peakPos - risetime;
0156   if (par[1] < 0.)
0157     par[1] = 0.;
0158   par[2] = 4.;
0159   par[3] = 1.5;
0160   par[4] = pedestal;
0161   par[5] = 0;
0162   fits.SetParameters(par);
0163   fits.SetParLimits(0, peakval * 0.9, peakval * 1.1);
0164   fits.SetParLimits(1, 0, 24);
0165   fits.SetParLimits(2, 2, 4.);
0166   fits.SetParLimits(3, 1., 2.);
0167   fits.SetParLimits(4, pedestal - abs(peakval), pedestal + abs(peakval));
0168   //  fits.SetParLimits(5, - abs(peakval),  + abs(peakval));
0169   fits.FixParameter(5, 0);
0170
0171   // Saturation correction - Abhisek
0172   for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
0173     if ((gpulse.GetY())[ipoint] == 0 or
0174         (gpulse.GetY())[ipoint] >=
0175             4090)  // drop point if touching max or low limit on ADCs
0176     {
0177       gpulse.RemovePoint(ipoint);
0178       ipoint--;
0179     }
0180
0181   gpulse.Fit(&fits, "MQRN0", "goff", 0., (double) NSAMPLES);
0182
0183   if (verbosity)
0184   {
0185     TCanvas *canvas = new TCanvas("PROTOTYPE3_FEM_SampleFit_PowerLawExp",
0186                                   "PROTOTYPE3_FEM::SampleFit_PowerLawExp");
0187     gpulse.DrawClone("ap*l");
0188     fits.DrawClone("same");
0189     fits.Print();
0190     canvas->Update();
0191     sleep(1);
0192   }
0193
0194   //  peak = fits.GetParameter(0); // not exactly peak height
0195   peak =
0196       (fits.GetParameter(0) *
0197        pow(fits.GetParameter(2) / fits.GetParameter(3), fits.GetParameter(2))) /
0198       exp(fits.GetParameter(
0199           2));  // exact peak height is (p0*Power(p2/p3,p2))/Power(E,p2)
0200
0201   //  peak_sample = fits.GetParameter(1); // signal start time
0202   peak_sample = fits.GetParameter(1) +
0203                 fits.GetParameter(2) / fits.GetParameter(3);  // signal peak time
0204
0205   // peak integral = p0*Power(p3,-1 - p2)*Gamma(1 + p2). Note yet used in output
0206
0207   pedstal = fits.GetParameter(4);
0208
0209   return true;
0210 }
0211
0212 double PROTOTYPE3_FEM::SignalShape_PowerLawExp(double *x, double *par)
0213 {
0214   double pedestal =
0215       par[4] + ((x[0] - 1.5 * par[1]) > 0) *
0216                    par[5];  // quick fix on exting tails on the signal function
0217   if (x[0] < par[1])
0218     return pedestal;
0219   // double  signal =
0220   // (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0221   double signal =
0222       par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
0223   return pedestal + signal;
0224 }