Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PROTOTYPE3_FEM.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PROTOTYPE3_FEM.cc
1 // $Id: $
2 
11 #include "PROTOTYPE3_FEM.h"
12 
13 #include <TCanvas.h>
14 #include <TF1.h>
15 #include <TGraph.h>
16 
17 #include <unistd.h> // for sleep
18 #include <cassert>
19 #include <cmath>
20 #include <cstdlib> // for exit
21 #include <iostream>
22 
23 using namespace std;
24 
25 int PROTOTYPE3_FEM::GetHBDCh(const std::string &caloname, int i_column,
26  int i_row)
27 {
28  if (caloname == "HCALIN")
29  {
30  return 64 + 8 * i_column + 2 * i_row;
31  }
32  else if (caloname == "HCALOUT")
33  {
34  return 112 + 8 * i_column + 2 * i_row;
35  }
36  else if (caloname == "EMCAL_PROTOTYPE2")
37  {
38  // EMcal cable mapping from John haggerty
39  assert(i_column >= 0);
40  assert(i_column < NCH_EMCAL_COLUMNS);
41  assert(i_row >= 0);
42  assert(i_row < NCH_EMCAL_ROWS);
43 
44  static int canmap[NCH_EMCAL_ROWS * NCH_EMCAL_COLUMNS] = {
45  // 1 ... 15
46  10 + 48, 11 + 48, 8 + 48, 9 + 48, 14 + 48, 15 + 48, 12 + 48, 13 + 48,
47  // 9 ... 16
48  2 + 48, 3 + 48, 0 + 48, 1 + 48, 6 + 48, 7 + 48, 4 + 48, 5 + 48,
49 
50  // 17 ... 24
51  10 + 32, 11 + 32, 8 + 32, 9 + 32, 14 + 32, 15 + 32, 12 + 32, 13 + 32,
52  // 25 ... 32
53  2 + 32, 3 + 32, 0 + 32, 1 + 32, 6 + 32, 7 + 32, 4 + 32, 5 + 32,
54 
55  // 33 ... 40
56  10 + 16, 11 + 16, 8 + 16, 9 + 16, 14 + 16, 15 + 16, 12 + 16, 13 + 16,
57  // 41 42 43 44 45 46 47 48
58  2 + 16, 3 + 16, 0 + 16, 1 + 16, 6 + 16, 7 + 16, 4 + 16, 5 + 16,
59 
60  // 49 50 51 52 53 54 55 56
61  10, 11, 8, 9, 14, 15, 12, 13,
62  // 57 58 59 60 61 62 63 64
63  2, 3, 0, 1, 6, 7, 4, 5};
64 
65  const int tower_index =
66  i_column + NCH_EMCAL_COLUMNS * (NCH_EMCAL_ROWS - 1 - i_row);
67 
68  assert(tower_index >= 0);
69  assert(tower_index < NCH_EMCAL_ROWS * NCH_EMCAL_COLUMNS);
70 
71  return canmap[tower_index];
72  }
73  else if (caloname == "EMCAL_HIGHETA")
74  {
75  // EMcal cable mapping from John haggerty
76  assert(i_column >= 0);
77  assert(i_column < NCH_EMCAL_COLUMNS);
78  assert(i_row >= 0);
79  assert(i_row < NCH_EMCAL_ROWS);
80 
81  static int canmap[] = {
82  // >
83  // https://docdb.sphenix.bnl.gov/0000/000034/001/T1044-2017a-2.xlsx
84  // Sean and John spot checked a number of these towers at
85  // BNL. There could be mistakes, but cosmics look
86  // reasonable.
87  // Front view (looking downstream, same as above) but in HBD
88  // channel numbers
89  // 3 , 2 , 19 , 18 , 35 , 34 , 51 , 50,
90  // 1 , 0 , 17 , 16 , 33 , 32 , 49 , 48,
91  // 7 , 6 , 23 , 22 , 39 , 38 , 55 , 54,
92  // 5 , 4 , 21 , 20 , 37 , 36 , 53 , 52,
93  // 11 , 10 , 27 , 26 , 43 , 42 , 59 , 58,
94  // 9 , 8 , 25 , 24 , 41 , 40 , 57 , 56,
95  // 15 , 14 , 31 , 30 , 47 , 46 , 63 , 62,
96  // 13 , 12 , 29 , 28 , 45 , 44 , 61 , 60,
97  3, 2, 19, 18, 35, 34, 51, 50, 1, 0, 17, 16, 33, 32, 49, 48, 7,
98  6, 23, 22, 39, 38, 55, 54, 5, 4, 21, 20, 37, 36, 53, 52, 11, 10,
99  27, 26, 43, 42, 59, 58, 9, 8, 25, 24, 41, 40, 57, 56, 15, 14, 31,
100  30, 47, 46, 63, 62, 13, 12, 29, 28, 45, 44, 61, 60, 0};
101 
102  const int tower_index =
103  i_column + NCH_EMCAL_COLUMNS * (NCH_EMCAL_ROWS - 1 - i_row);
104 
105  assert(tower_index >= 0);
106  assert(tower_index < NCH_EMCAL_ROWS * NCH_EMCAL_COLUMNS);
107 
108  return canmap[tower_index];
109  }
110 
111  std::cout << "PROTOTYPE3_FEM::GetHBDCh - invalid input caloname " << caloname
112  << " i_column " << i_column << " i_row " << i_row << std::endl;
113  exit(1);
114  return -9999;
115 }
116 
117 bool PROTOTYPE3_FEM::SampleFit_PowerLawExp( //
118  const std::vector<double> &samples, //
119  double &peak, //
120  double &peak_sample, //
121  double &pedstal, //
122  const int verbosity)
123 {
124  int peakPos = 0.;
125 
126  assert(samples.size() == NSAMPLES);
127 
128  TGraph gpulse(NSAMPLES);
129  for (int i = 0; i < NSAMPLES; i++)
130  {
131  (gpulse.GetX())[i] = i;
132 
133  (gpulse.GetY())[i] = samples[i];
134  }
135 
136  double pedestal = gpulse.GetY()[0]; //(double) PEDESTAL;
137  double peakval = pedestal;
138  const double risetime = 4;
139 
140  for (int iSample = 0; iSample < NSAMPLES; iSample++)
141  {
142  if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
143  {
144  peakval = gpulse.GetY()[iSample];
145  peakPos = iSample;
146  }
147  }
148  peakval -= pedestal;
149 
150  // fit function
151  TF1 fits("f_SignalShape_PowerLawExp", SignalShape_PowerLawExp, 0., 24., 6);
152 
153  double par[6] = {0};
154  par[0] = peakval; // /3.;
155  par[1] = peakPos - risetime;
156  if (par[1] < 0.)
157  par[1] = 0.;
158  par[2] = 4.;
159  par[3] = 1.5;
160  par[4] = pedestal;
161  par[5] = 0;
162  fits.SetParameters(par);
163  fits.SetParLimits(0, peakval * 0.9, peakval * 1.1);
164  fits.SetParLimits(1, 0, 24);
165  fits.SetParLimits(2, 2, 4.);
166  fits.SetParLimits(3, 1., 2.);
167  fits.SetParLimits(4, pedestal - abs(peakval), pedestal + abs(peakval));
168  // fits.SetParLimits(5, - abs(peakval), + abs(peakval));
169  fits.FixParameter(5, 0);
170 
171  // Saturation correction - Abhisek
172  for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
173  if ((gpulse.GetY())[ipoint] == 0 or
174  (gpulse.GetY())[ipoint] >=
175  4090) // drop point if touching max or low limit on ADCs
176  {
177  gpulse.RemovePoint(ipoint);
178  ipoint--;
179  }
180 
181  gpulse.Fit(&fits, "MQRN0", "goff", 0., (double) NSAMPLES);
182 
183  if (verbosity)
184  {
185  TCanvas *canvas = new TCanvas("PROTOTYPE3_FEM_SampleFit_PowerLawExp",
186  "PROTOTYPE3_FEM::SampleFit_PowerLawExp");
187  gpulse.DrawClone("ap*l");
188  fits.DrawClone("same");
189  fits.Print();
190  canvas->Update();
191  sleep(1);
192  }
193 
194  // peak = fits.GetParameter(0); // not exactly peak height
195  peak =
196  (fits.GetParameter(0) *
197  pow(fits.GetParameter(2) / fits.GetParameter(3), fits.GetParameter(2))) /
198  exp(fits.GetParameter(
199  2)); // exact peak height is (p0*Power(p2/p3,p2))/Power(E,p2)
200 
201  // peak_sample = fits.GetParameter(1); // signal start time
202  peak_sample = fits.GetParameter(1) +
203  fits.GetParameter(2) / fits.GetParameter(3); // signal peak time
204 
205  // peak integral = p0*Power(p3,-1 - p2)*Gamma(1 + p2). Note yet used in output
206 
207  pedstal = fits.GetParameter(4);
208 
209  return true;
210 }
211 
212 double PROTOTYPE3_FEM::SignalShape_PowerLawExp(double *x, double *par)
213 {
214  double pedestal =
215  par[4] + ((x[0] - 1.5 * par[1]) > 0) *
216  par[5]; // quick fix on exting tails on the signal function
217  if (x[0] < par[1])
218  return pedestal;
219  // double signal =
220  // (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
221  double signal =
222  par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
223  return pedestal + signal;
224 }