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