packages/Prototype2/PROTOTYPE2_FEM.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file PROTOTYPE2_FEM.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "PROTOTYPE2_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 PROTOTYPE2_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")
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
0074   std::cout << "PROTOTYPE2_FEM::GetHBDCh - invalid input caloname " << caloname
0075             << " i_column " << i_column << " i_row " << i_row << std::endl;
0076   exit(1);
0077   return -9999;
0078 }
0079
0080 bool PROTOTYPE2_FEM::SampleFit_PowerLawExp(  //
0081     const std::vector<double> &samples,      //
0082     double &peak,                            //
0083     double &peak_sample,                     //
0084     double &pedstal,                         //
0085     const int verbosity)
0086 {
0087   int peakPos = 0.;
0088
0089   assert(samples.size() == NSAMPLES);
0090
0091   TGraph gpulse(NSAMPLES);
0092   for (int i = 0; i < NSAMPLES; i++)
0093   {
0094     (gpulse.GetX())[i] = i;
0095
0096     (gpulse.GetY())[i] = samples[i];
0097   }
0098
0099   double pedestal = gpulse.GetY()[0];  //(double) PEDESTAL;
0100   double peakval = pedestal;
0101   const double risetime = 4;
0102
0103   for (int iSample = 0; iSample < NSAMPLES; iSample++)
0104   {
0105     if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
0106     {
0107       peakval = gpulse.GetY()[iSample];
0108       peakPos = iSample;
0109     }
0110   }
0111   peakval -= pedestal;
0112
0113   // fit function
0114   TF1 fits("f_SignalShape_PowerLawExp", SignalShape_PowerLawExp, 0., 24., 6);
0115
0116   double par[6] = {0};
0117   par[0] = peakval;  // /3.;
0118   par[1] = peakPos - risetime;
0119   if (par[1] < 0.)
0120     par[1] = 0.;
0121   par[2] = 4.;
0122   par[3] = 1.5;
0123   par[4] = pedestal;
0124   par[5] = 0;
0125   fits.SetParameters(par);
0126   fits.SetParLimits(0, peakval * 0.9, peakval * 1.1);
0127   fits.SetParLimits(1, 0, 24);
0128   fits.SetParLimits(2, 2, 4.);
0129   fits.SetParLimits(3, 1., 2.);
0130   fits.SetParLimits(4, pedestal - abs(peakval), pedestal + abs(peakval));
0131   //  fits.SetParLimits(5, - abs(peakval),  + abs(peakval));
0132   fits.FixParameter(5, 0);
0133
0134   // Saturation correction - Abhisek
0135   for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
0136     if ((gpulse.GetY())[ipoint] == 0)
0137     {
0138       gpulse.RemovePoint(ipoint);
0139       ipoint--;
0140     }
0141
0142   gpulse.Fit(&fits, "MQRN0", "goff", 0., (double) NSAMPLES);
0143
0144   if (verbosity)
0145   {
0146     TCanvas *canvas = new TCanvas("PROTOTYPE2_FEM_SampleFit_PowerLawExp",
0147                                   "PROTOTYPE2_FEM::SampleFit_PowerLawExp");
0148     gpulse.DrawClone("ap*l");
0149     fits.DrawClone("same");
0150     fits.Print();
0151     canvas->Update();
0152     sleep(1);
0153   }
0154
0155   //  peak = fits.GetParameter(0); // not exactly peak height
0156   peak =
0157       (fits.GetParameter(0) *
0158        pow(fits.GetParameter(2) / fits.GetParameter(3), fits.GetParameter(2))) /
0159       exp(fits.GetParameter(
0160           2));  // exact peak height is (p0*Power(p2/p3,p2))/Power(E,p2)
0161
0162   //  peak_sample = fits.GetParameter(1); // signal start time
0163   peak_sample = fits.GetParameter(1) +
0164                 fits.GetParameter(2) / fits.GetParameter(3);  // signal peak time
0165
0166   // peak integral = p0*Power(p3,-1 - p2)*Gamma(1 + p2). Note yet used in output
0167
0168   pedstal = fits.GetParameter(4);
0169
0170   return true;
0171 }
0172
0173 double PROTOTYPE2_FEM::SignalShape_PowerLawExp(double *x, double *par)
0174 {
0175   double pedestal =
0176       par[4] + ((x[0] - 1.5 * par[1]) > 0) *
0177                    par[5];  // quick fix on exting tails on the signal function
0178   if (x[0] < par[1])
0179     return pedestal;
0180   // double  signal =
0181   // (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0182   double signal =
0183       par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
0184   return pedestal + signal;
0185 }