EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BEmcProfile.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BEmcProfile.cc
1 #include "BEmcProfile.h"
2 #include "BEmcCluster.h"
3 
4 #include <TFile.h>
5 #include <TH1.h> // for TH1F
6 #include <TMath.h>
7 #include <TROOT.h>
8 
9 #include <boost/format.hpp>
10 
11 #include <cmath> // for sqrt, log, pow, fabs
12 #include <cstdlib> // for abs
13 #include <iostream>
14 #include <memory> // for allocator_traits<>::value_type
15 
16 using namespace std;
17 
18 const int NP = 4; // Number of profiles in a bin
19 
20 BEmcProfile::BEmcProfile(const string& fname)
21  : bloaded(false)
22  , thresh(0.01)
23  , nen(0)
24  , nth(0)
25  , energy_array(0)
26  , theta_array(0)
27  , hmean(0)
28  , hsigma(0)
29  , m_Verbosity(0)
30 {
31  TFile* f = new TFile(fname.c_str());
32  if (f->IsZombie())
33  {
34  cout << "BEmcProfile: Error when opening profile data file " << fname << endl;
35  return;
36  }
37 
38  gROOT->cd();
39 
40  TH1F* hen = (TH1F*) f->Get("hen");
41  if (!hen)
42  {
43  cout << "BEmcProfile: Error when loading profile data: hen" << endl;
44  f->Close();
45  return;
46  }
47 
48  TH1F* hth = (TH1F*) f->Get("hth");
49  if (!hth)
50  {
51  cout << "BEmcProfile: Error when loading profile data: hth" << endl;
52  f->Close();
53  return;
54  }
55 
56  nen = hen->GetNbinsX();
57  nth = hth->GetNbinsX();
58 
59  energy_array = new float[nen];
60  for (int i = 0; i < nen; i++) energy_array[i] = hen->GetBinContent(i + 1);
61  theta_array = new float[nth];
62  for (int i = 0; i < nth; i++) theta_array[i] = hth->GetBinContent(i + 1);
63 
64  if (Verbosity())
65  {
66  cout << "BEmcProfile: Loading profile data from file " << fname << endl;
67 
68  cout << " " << nen << " bins in energy: ";
69  for (int i = 0; i < nen; i++)
70  {
71  cout << boost::format("%4.1f, ") % energy_array[i];
72  }
73  cout << endl;
74  cout << " " << nth << " bins in theta: ";
75  for (int i = 0; i < nth; i++)
76  {
77  cout << boost::format("%4.2f, ") % theta_array[i];
78  }
79  cout << endl;
80  }
81  hmean = new TH1F*[nen * nth * NP];
82  hsigma = new TH1F*[nen * nth * NP];
83  hr4 = new TH1F*[nen * nth];
84 
85  TH1F* hh;
86  int ii = 0;
87  int ii2 = 0;
88 
89  string hname;
90  for (int it = 0; it < nth; it++)
91  {
92  for (int ie = 0; ie < nen; ie++)
93  {
94  for (int ip = 0; ip < NP; ip++)
95  {
96  hname = boost::str(boost::format("hmean%d_en%d_t%d") % (ip + 1) % ie % it);
97  hh = (TH1F*) f->Get(hname.c_str());
98  if (!hh)
99  {
100  cout << "BEmcProfile: Could not load histogram " << hname
101  << ", Error when loading profile data for hmean it = "
102  << it << ", ie = " << ie << "ip = " << ip << endl;
103  f->Close();
104  return;
105  }
106  hmean[ii] = (TH1F*) hh->Clone();
107 
108  hname = boost::str(boost::format("hsigma%d_en%d_t%d") % (ip + 1) % ie % it);
109  hh = (TH1F*) f->Get(hname.c_str());
110  if (!hh)
111  {
112  cout << "BEmcProfile: Could not load histogram " << hname
113  << ", Error when loading profile data for hsigma it = "
114  << it << ", ie = " << ie << ", ip = " << ip << endl;
115  f->Close();
116  return;
117  }
118  hsigma[ii] = (TH1F*) hh->Clone();
119 
120  ii++;
121  }
122 
123  hname = boost::str(boost::format("hr4_en%d_t%d") % ie % it);
124  hh = (TH1F*) f->Get(hname.c_str());
125 
126  if (!hh)
127  {
128  cout << "BEmcProfile: Error when loading profile data for hr4 it = "
129  << it << ", ie = " << ie << endl;
130  f->Close();
131  return;
132  }
133  hr4[ii2] = (TH1F*) hh->Clone();
134  ii2++;
135 
136  }
137  }
138 
139  f->Close();
140  bloaded = true;
141 }
142 
143 BEmcProfile::~BEmcProfile()
144 {
145  if (bloaded)
146  {
147  delete[] energy_array;
148  delete[] theta_array;
149  for (int i = 0; i < nth * nen * NP; i++)
150  {
151  delete hmean[i];
152  delete hsigma[i];
153  }
154  delete[] hmean;
155  delete[] hsigma;
156 
157  for (int i = 0; i < nth * nen; i++)
158  {
159  delete hr4[i];
160  }
161  delete[] hr4;
162  }
163 }
164 
165 float BEmcProfile::GetProb(std::vector<EmcModule>* plist, int NX, float en, float theta, float phi)
166 {
167  float enoise = 0.01; // 10 MeV per tower
168 
169  if (!bloaded)
170  {
171  // cout << "Error in BEmcProfile::GetProb: profiles not loaded" << endl;
172  return -1;
173  }
174 
175  int nn = plist->size();
176  if (nn <= 0) return -1;
177 
178  // z coordinate below means x coordinate
179 
180  float ee;
181  int ich, iy, iz;
182 
183  int iy0 = -1, iz0 = -1;
184  float emax = 0;
185  for (int i = 0; i < nn; i++)
186  {
187  ee = (*plist)[i].amp;
188  if (ee > emax)
189  {
190  emax = ee;
191  ich = (*plist)[i].ich;
192  iy0 = ich / NX;
193  iz0 = ich % NX;
194  }
195  }
196  if (emax <= 0) return -1;
197 
198  float etot = 0;
199  float sz = 0;
200  float sy = 0;
201  for (int i = 0; i < nn; i++)
202  {
203  ee = (*plist)[i].amp;
204  ich = (*plist)[i].ich;
205  iy = ich / NX;
206  iz = ich % NX;
207  if (ee > thresh && abs(iz - iz0) <= 3 && abs(iy - iy0) <= 3)
208  {
209  etot += ee;
210  sz += iz * ee;
211  sy += iy * ee;
212  }
213  }
214  float zcg = sz / etot;
215  float ycg = sy / etot;
216  int iz0cg = int(zcg + 0.5);
217  int iy0cg = int(ycg + 0.5);
218  float ddz = fabs(zcg - iz0cg);
219  float ddy = fabs(ycg - iy0cg);
220 
221  int isz = 1;
222  if (zcg - iz0cg < 0) isz = -1;
223  int isy = 1;
224  if (ycg - iy0cg < 0) isy = -1;
225 
226  // 4 central towers: 43
227  // 12
228  // Tower 1 - central one
229  float e1, e2, e3, e4;
230  e1 = GetTowerEnergy(iy0cg, iz0cg, plist, NX);
231  e2 = GetTowerEnergy(iy0cg, iz0cg + isz, plist, NX);
232  e3 = GetTowerEnergy(iy0cg + isy, iz0cg + isz, plist, NX);
233  e4 = GetTowerEnergy(iy0cg + isy, iz0cg, plist, NX);
234  if (e1 < thresh) e1 = 0;
235  if (e2 < thresh) e2 = 0;
236  if (e3 < thresh) e3 = 0;
237  if (e4 < thresh) e4 = 0;
238 
239  float e1t = (e1 + e2 + e3 + e4) / etot;
240  float e2t = (e1 + e2 - e3 - e4) / etot;
241  float e3t = (e1 - e2 - e3 + e4) / etot;
242  float e4t = (e3) / etot;
243  // float rr = sqrt((0.5-ddz)*(0.5-ddz)+(0.5-ddy)*(0.5-ddy));
244 
245  // Predicted values
246  float ep[NP];
247  float err[NP];
248  for (int ip = 0; ip < NP; ip++)
249  {
250  PredictEnergy(ip, en, theta, phi, ddz, ddy, ep[ip], err[ip]);
251  if (ep[ip] < 0) return -1;
252  if (ip < 3)
253  err[ip] = sqrt(err[ip] * err[ip] + 4 * enoise * enoise / etot / etot);
254  else
255  err[ip] = sqrt(err[ip] * err[ip] + 1 * enoise * enoise / etot / etot);
256  }
257 
258  float chi2 = 0.;
259  chi2 += (ep[0] - e1t) * (ep[0] - e1t) / err[0] / err[0];
260  chi2 += (ep[1] - e2t) * (ep[1] - e2t) / err[1] / err[1];
261  chi2 += (ep[2] - e3t) * (ep[2] - e3t) / err[2] / err[2];
262  chi2 += (ep[3] - e4t) * (ep[3] - e4t) / err[3] / err[3];
263  int ndf = 4;
264 
265  float prob = TMath::Prob(chi2, ndf);
266 
267  return prob;
268 }
269 
270 /*
271 float BEmcProfile::GetProbTest(std::vector<EmcModule>* plist, int NX, float en, float theta, float& test_rr, float& test_et, float& test_ep, float& test_err)
272 // This is only for test purposes. Can be removed if not needed
273 {
274  int nn = plist->size();
275  if( nn <= 0 ) return -1;
276 
277  // z coordinate below means x coordinate
278 
279  float ee;
280  int ich, iy, iz;
281 
282  int iy0=-1, iz0=-1;
283  float emax=0;
284  for( int i=0; i<nn; i++ ) {
285  ee = (*plist)[i].amp;
286  if( ee>emax ) {
287  emax = ee;
288  ich = (*plist)[i].ich;
289  iy0 = ich/NX;
290  iz0 = ich%NX;
291  }
292  }
293  if( emax<=0 ) return -1;
294 
295  float etot=0;
296  float sz=0;
297  float sy=0;
298  for( int i=0; i<nn; i++ ) {
299  ee = (*plist)[i].amp;
300  ich = (*plist)[i].ich;
301  iy = ich/NX;
302  iz = ich%NX;
303  if( ee>thresh && abs(iz-iz0)<=3 && abs(iy-iy0)<=3 ) {
304  etot += ee;
305  sz += iz*ee;
306  sy += iy*ee;
307  }
308  }
309  float zcg = sz/etot;
310  float ycg = sy/etot;
311  int iz0cg = int(zcg+0.5);
312  int iy0cg = int(ycg+0.5);
313  float ddz = fabs(zcg-iz0cg);
314  float ddy = fabs(ycg-iy0cg);
315 
316  int isz = 1;
317  if( zcg-iz0cg < 0 ) isz = -1;
318  int isy = 1;
319  if( ycg-iy0cg < 0 ) isy = -1;
320 
321  // 4 central towers: 43
322  // 12
323  // Tower 1 - central one
324  float e1, e2, e3, e4;
325  e1 = GetTowerEnergy(iy0cg, iz0cg, plist, NX);
326  e2 = GetTowerEnergy(iy0cg, iz0cg+isz, plist, NX);
327  e3 = GetTowerEnergy(iy0cg+isy,iz0cg+isz, plist, NX);
328  e4 = GetTowerEnergy(iy0cg+isy,iz0cg, plist, NX);
329  if( e1<thresh ) e1 = 0;
330  if( e2<thresh ) e2 = 0;
331  if( e3<thresh ) e3 = 0;
332  if( e4<thresh ) e4 = 0;
333 
334  float e1t = (e1+e2+e3+e4)/etot;
335  float e2t = (e1+e2-e3-e4)/etot;
336  float e3t = (e1-e2-e3+e4)/etot;
337  float e4t = (e3)/etot;
338  float rr = sqrt((0.5-ddz)*(0.5-ddz)+(0.5-ddy)*(0.5-ddy));
339 
340  // Predicted values
341  float ep[NP];
342  float err[NP];
343  for( int ip=0; ip<NP; ip++ )
344  PredictEnergy(ip, en, theta, ddz, ddy, ep[ip], err[ip]);
345 
346  float chi2 = 0.;
347  chi2 += (ep[0]-e1t)*(ep[0]-e1t)/err[0]/err[0];
348  chi2 += (ep[1]-e2t)*(ep[1]-e2t)/err[1]/err[1];
349  chi2 += (ep[2]-e3t)*(ep[2]-e3t)/err[2]/err[2];
350  chi2 += (ep[3]-e4t)*(ep[3]-e4t)/err[3]/err[3];
351  int ndf = 4;
352 
353  float prob = TMath::Prob(chi2, ndf);
354 
355  test_rr = rr;
356  test_et = e1t;
357  test_ep = ep[0];
358  test_err = err[0];
359 
360  return prob;
361 }
362 */
363 
364 void BEmcProfile::PredictEnergy(int ip, float energy, float theta, float /*phi*/, float ddz, float ddy, float& ep, float& err)
365 // ip changes from 0 to NP-1, meaning the profile index 1,2,..,NP
366 {
367  ep = err = -1;
368 
369  if (!bloaded)
370  {
371  //cout << "Error in BEmcProfile::PredictEnergy: profiles not loaded" << endl;
372  return;
373  }
374 
375  if (ip < 0 || ip >= NP)
376  {
377  cout << "Error in BEmcProfile::PredictEnergy: profile index = "
378  << ip << " but should be from 0 to " << NP - 1 << endl;
379  return;
380  }
381 
382  if (energy <= 0)
383  {
384  cout << "Error in BEmcProfile::PredictEnergy: energy = "
385  << energy << " but should be >0" << endl;
386  return;
387  }
388  if (theta < 0)
389  {
390  cout << "Error in BEmcProfile::PredictEnergy: theta = "
391  << theta << " but should be >=0" << endl;
392  return;
393  }
394 
395  if (ddz < 0 || ddz > 0.5)
396  {
397  cout << "Error in BEmcProfile::PredictEnergy: ddz = "
398  << ddz << " but should be from 0 to 0.5" << endl;
399  }
400 
401  if (ddy < 0 || ddy > 0.5)
402  {
403  cout << "Error in BEmcProfile::PredictEnergy: ddy = "
404  << ddy << " but should be from 0 to 0.5" << endl;
405  }
406 
407  // Safety margin (slightly away from bin edge)
408  // if (ddz < 0.021) ddz = 0.021;
409  // if (ddz > 0.489) ddz = 0.489;
410  // if (ddy < 0.021) ddy = 0.021;
411  // if (ddy > 0.489) ddy = 0.489;
412 
413  // Energy bin
414  int ie2 = 0;
415  while (ie2 < nen && energy > energy_array[ie2]) ie2++;
416  if (ie2 == 0)
417  ie2 = 1;
418  else if (ie2 >= nen)
419  ie2 = nen - 1;
420  int ie1 = ie2 - 1;
421  // int ie1 = ie2-2; // For a test()
422 
423  // Theta bin
424  int it2 = 0;
425  while (it2 < nth && theta > theta_array[it2]) it2++;
426  if (it2 == 0)
427  it2 = 1;
428  else if (it2 >= nth)
429  it2 = nth - 1;
430  int it1 = it2 - 1;
431  // int it1 = it2-2; // For a test()
432 
433  // cout << "Energy bin= " << ie1 << " " << ie2 << " ("
434  // << energy << ") Theta bin= " << it1 << " " << it2
435  // << " (" << theta << ")" << endl;
436 
437  float rr = sqrt((0.5 - ddz) * (0.5 - ddz) + (0.5 - ddy) * (0.5 - ddy));
438 
439  float xx = rr;
440  if (ip == 1)
441  xx = ddy;
442  else if (ip == 2)
443  xx = ddz;
444 
445  float en1 = energy_array[ie1];
446  float en2 = energy_array[ie2];
447  float th1 = theta_array[it1];
448  float th2 = theta_array[it2];
449 
450  // 1st index - ie, second index - it
451  int ii11 = ip + ie1 * NP + it1 * nen * NP;
452  int ii21 = ip + ie2 * NP + it1 * nen * NP;
453  int ii12 = ip + ie1 * NP + it2 * nen * NP;
454  int ii22 = ip + ie2 * NP + it2 * nen * NP;
455 
456  int ibin = hmean[ii11]->FindBin(xx);
457 
458  // Log (1/sqrt) energy dependence of mean (sigma)
459  //
460  float pr11 = hmean[ii11]->GetBinContent(ibin);
461  float pr21 = hmean[ii21]->GetBinContent(ibin);
462  float prt1 = pr11 + (pr21 - pr11) / (log(en2) - log(en1)) * (log(energy) - log(en1));
463  if (prt1 < 0) prt1 = 0;
464 
465  float er11 = hsigma[ii11]->GetBinContent(ibin);
466  float er21 = hsigma[ii21]->GetBinContent(ibin);
467  float ert1 = er11 + (er21 - er11) / (1. / sqrt(en2) - 1. / sqrt(en1)) * (1. / sqrt(energy) - 1. / sqrt(en1));
468  if (ert1 < 0) ert1 = 0;
469 
470  float pr12 = hmean[ii12]->GetBinContent(ibin);
471  float pr22 = hmean[ii22]->GetBinContent(ibin);
472  float prt2 = pr12 + (pr22 - pr12) / (log(en2) - log(en1)) * (log(energy) - log(en1));
473  if (prt2 < 0) prt2 = 0;
474 
475  float er12 = hsigma[ii12]->GetBinContent(ibin);
476  float er22 = hsigma[ii22]->GetBinContent(ibin);
477  float ert2 = er12 + (er22 - er12) / (1. / sqrt(en2) - 1. / sqrt(en1)) * (1. / sqrt(energy) - 1. / sqrt(en1));
478  if (ert2 < 0) ert2 = 0;
479 
480  // Quadratic theta dependence of mean and sigma
481  //
482  float pr = prt1 + (prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
483  if (pr < 0) pr = 0;
484  float er = ert1 + (ert2 - ert1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
485  if (er < 0) er = 0;
486 
487  // Additional error due to binning in xx
488  //
489  int ibin1 = ibin;
490  if (ibin > 1) ibin1 = ibin - 1;
491  int ibin2 = ibin;
492  if (ibin < hmean[ii11]->GetNbinsX())
493  {
494  if (hmean[ii11]->GetBinContent(ibin + 1) > 0)
495  {
496  ibin2 = ibin + 1;
497  }
498  }
499  float dd = (hmean[ii11]->GetBinContent(ibin2) -
500  hmean[ii11]->GetBinContent(ibin1)) /
501  2.;
502  // if( fabs(dd)>er )
503  // {
504  // cout << "ie = " << ie1 << ", it = " << it1 << ", bin = "
505  // << ibin << ": " << er << " " << dd << endl;
506  // }
507  er = sqrt(er * er + dd * dd);
508 
509  ep = pr;
510  err = er;
511 }
512 
513 
514 float BEmcProfile::PredictEnergyR(float energy, float theta, float /*phi*/, float rr)
515 {
516 
517  if (!bloaded)
518  {
519  // cout << "Error in BEmcProfile::PredictEnergyR: profiles not loaded" << endl;
520  return 0;
521  }
522 
523  if (energy <= 0)
524  {
525  cout << "Error in BEmcProfile::PredictEnergyR: energy = " << energy
526  << " but should be >0" << endl;
527  return 0;
528  }
529  if (theta < 0)
530  {
531  cout << "Error in BEmcProfile::PredictEnergyR: theta = " << theta
532  << " but should be >=0" << endl;
533  return 0;
534  }
535 
536  if (rr < 1 )
537  {
538  cout << "Error in BEmcProfile::PredictEnergyR: rr = " << rr
539  << " but should be >1" << endl;
540  }
541 
542  // Energy bin
543  int ie2 = 0;
544  while (ie2 < nen && energy > energy_array[ie2]) ie2++;
545  if (ie2 == 0)
546  ie2 = 1;
547  else if (ie2 >= nen)
548  ie2 = nen - 1;
549  int ie1 = ie2 - 1;
550  // int ie1 = ie2-2; // For a test()
551 
552  // Theta bin
553  int it2 = 0;
554  while (it2 < nth && theta > theta_array[it2]) it2++;
555  if (it2 == 0)
556  it2 = 1;
557  else if (it2 >= nth)
558  it2 = nth - 1;
559  int it1 = it2 - 1;
560  // int it1 = it2-2; // For a test()
561 
562  // printf("Energy bin= %d %d (%f) Theta bin= %d %d (%f)\n",ie1,ie2,energy,it1,it2,theta);
563 
564  float en1 = energy_array[ie1];
565  float en2 = energy_array[ie2];
566  float th1 = theta_array[it1];
567  float th2 = theta_array[it2];
568 
569  // 1st index - ie, second index - it
570  int ii11 = ie1 + it1 * nen;
571  int ii21 = ie2 + it1 * nen;
572  int ii12 = ie1 + it2 * nen;
573  int ii22 = ie2 + it2 * nen;
574 
575  int ibin = hr4[ii11]->FindBin(rr);
576 
577  // Log (1/sqrt) energy dependence of mean (sigma)
578  //
579  float pr11 = hr4[ii11]->GetBinContent(ibin);
580  float pr21 = hr4[ii21]->GetBinContent(ibin);
581  float prt1 = pr11 + (pr21 - pr11) / (log(en2) - log(en1)) * (log(energy) - log(en1));
582  if (prt1 < 0) prt1 = 0;
583 
584  float pr12 = hr4[ii12]->GetBinContent(ibin);
585  float pr22 = hr4[ii22]->GetBinContent(ibin);
586  float prt2 = pr12 + (pr22 - pr12) / (log(en2) - log(en1)) * (log(energy) - log(en1));
587  if (prt2 < 0) prt2 = 0;
588 
589  // Quadratic theta dependence of mean and sigma
590  //
591  float pr = prt1 + (prt2 - prt1) / (pow(th2, 2) - pow(th1, 2)) * (pow(theta, 2) - pow(th1, 2));
592  if (pr < 0) pr = 0;
593 
594  return pr;
595 }
596 
597 
598 float BEmcProfile::GetTowerEnergy(int iy, int iz, std::vector<EmcModule>* plist, int NX)
599 {
600  int nn = plist->size();
601  if (nn <= 0) return 0;
602 
603  for (int i = 0; i < nn; i++)
604  {
605  int ich = (*plist)[i].ich;
606  int iyt = ich / NX;
607  int izt = ich % NX;
608  if (iy == iyt && iz == izt)
609  {
610  return (*plist)[i].amp;
611  }
612  }
613  return 0;
614 }