Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Proto3ShowerCalib.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Proto3ShowerCalib.C
1 #include "Proto3ShowerCalib.h"
2 #include "TemperatureCorrection.h"
3 
4 #include <prototype3/RawTower_Temperature.h>
5 #include <prototype3/RawTower_Prototype3.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
8 #include <pdbcalbase/PdbParameterMap.h>
9 #include <phparameter/PHParameters.h>
10 #include <ffaobjects/EventHeader.h>
11 
12 #include <fun4all/SubsysReco.h>
13 #include <fun4all/Fun4AllServer.h>
14 #include <fun4all/PHTFileServer.h>
15 #include <phool/PHCompositeNode.h>
16 #include <fun4all/Fun4AllReturnCodes.h>
17 #include <phool/getClass.h>
18 #include <fun4all/Fun4AllHistoManager.h>
19 
20 #include <phool/PHCompositeNode.h>
21 
22 #include <g4main/PHG4TruthInfoContainer.h>
23 #include <g4main/PHG4Particle.h>
24 #include <g4main/PHG4VtxPoint.h>
25 
26 #include <TNtuple.h>
27 #include <TFile.h>
28 #include <TH1F.h>
29 #include <TH2F.h>
30 #include <TVector3.h>
31 #include <TLorentzVector.h>
32 
33 #include <exception>
34 #include <stdexcept>
35 #include <iostream>
36 #include <sstream>
37 #include <vector>
38 #include <set>
39 #include <algorithm>
40 #include <cassert>
41 #include <cmath>
42 
43 using namespace std;
44 
45 ClassImp(Proto3ShowerCalib::Eval_Cluster);
46 ClassImp(Proto3ShowerCalib::Eval_Run);
47 
48 Proto3ShowerCalib::Proto3ShowerCalib(const std::string &filename) :
49  SubsysReco("Proto3ShowerCalib"), _is_sim(false), _use_hodoscope_calibs(false),_filename(filename), _ievent(
50  0)
51 {
52 
53  verbosity = 1;
54 
55  _eval_run.reset();
56  _eval_3x3_raw.reset();
57  _eval_5x5_raw.reset();
58  _eval_3x3_prod.reset();
59  _eval_5x5_prod.reset();
60  _eval_3x3_temp.reset();
61  _eval_5x5_temp.reset();
62  _eval_3x3_recalib.reset();
63  _eval_5x5_recalib.reset();
64 
65  for (int col = 0; col < n_size; ++col)
66  for (int row = 0; row < n_size; ++row)
67  {
68  _recalib_const[make_pair(col, row)] = 0;
69  }
70 
71 }
72 
73 Proto3ShowerCalib::~Proto3ShowerCalib()
74 {
75 }
76 
77 Fun4AllHistoManager *
78 Proto3ShowerCalib::get_HistoManager()
79 {
80 
81  Fun4AllServer *se = Fun4AllServer::instance();
82  Fun4AllHistoManager *hm = se->getHistoManager("Proto3ShowerCalib_HISTOS");
83 
84  if (not hm)
85  {
86  cout
87  << "Proto3ShowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto3ShowerCalib_HISTOS"
88  << endl;
89  hm = new Fun4AllHistoManager("Proto3ShowerCalib_HISTOS");
90  se->registerHistoManager(hm);
91  }
92 
93  assert(hm);
94 
95  return hm;
96 }
97 
98 int
99 Proto3ShowerCalib::InitRun(PHCompositeNode *topNode)
100 {
101  if (verbosity)
102  cout << "Proto3ShowerCalib::InitRun" << endl;
103 
104  _ievent = 0;
105  cout<<"IS THIS SIMULATION: "<<_is_sim<<endl;
106  PHNodeIterator iter(topNode);
107  PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
108  "PHCompositeNode", "DST"));
109  if (!dstNode)
110  {
111  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
112  throw runtime_error(
113  "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
114  }
115  for(int i=0; i<nfmodbins; i++)
116  fmodbins[i] = 0. + i * 2. / (float)(nfmodbins-1);
117 
118  return Fun4AllReturnCodes::EVENT_OK;
119 }
120 
121 int
122 Proto3ShowerCalib::End(PHCompositeNode *topNode)
123 {
124  cout << "Proto3ShowerCalib::End - write to " << _filename << endl;
125  PHTFileServer::get().cd(_filename);
126 
127  Fun4AllHistoManager *hm = get_HistoManager();
128  assert(hm);
129  for (unsigned int i = 0; i < hm->nHistos(); i++)
130  hm->getHisto(i)->Write();
131 
132 // if (_T_EMCalTrk)
133 // _T_EMCalTrk->Write();
134 
135  fdata.close();
136 
137  return Fun4AllReturnCodes::EVENT_OK;
138 }
139 
140 int
141 Proto3ShowerCalib::Init(PHCompositeNode *topNode)
142 {
143 
144  _ievent = 0;
145 
146  cout << "Proto3ShowerCalib::get_HistoManager - Making PHTFileServer "
147  << _filename << endl;
148  PHTFileServer::get().open(_filename, "RECREATE");
149 
150  fdata.open(_filename + ".dat", std::fstream::out);
151 
152  Fun4AllHistoManager *hm = get_HistoManager();
153  assert(hm);
154 
155  TH2F * hCheck_Cherenkov = new TH2F("hCheck_Cherenkov", "hCheck_Cherenkov",
156  1000, -2000, 2000, 5, .5, 5.5);
157  hCheck_Cherenkov->GetYaxis()->SetBinLabel(1, "C1");
158  hCheck_Cherenkov->GetYaxis()->SetBinLabel(2, "C2 in");
159  hCheck_Cherenkov->GetYaxis()->SetBinLabel(3, "C2 out");
160  hCheck_Cherenkov->GetYaxis()->SetBinLabel(4, "C2 sum");
161  hm->registerHisto(hCheck_Cherenkov);
162 
163  TH1F * hNormalization = new TH1F("hNormalization", "hNormalization", 10, .5,
164  10.5);
165  hCheck_Cherenkov->GetXaxis()->SetBinLabel(1, "ALL");
166  hCheck_Cherenkov->GetXaxis()->SetBinLabel(2, "C2-e");
167  hCheck_Cherenkov->GetXaxis()->SetBinLabel(3, "trigger_veto_pass");
168  hCheck_Cherenkov->GetXaxis()->SetBinLabel(4, "valid_hodo_h");
169  hCheck_Cherenkov->GetXaxis()->SetBinLabel(5, "valid_hodo_v");
170  hCheck_Cherenkov->GetXaxis()->SetBinLabel(6, "good_e");
171  hCheck_Cherenkov->GetXaxis()->SetBinLabel(7, "good_data");
172  hm->registerHisto(hNormalization);
173 
174  hm->registerHisto(new TH1F("hCheck_Veto", "hCheck_Veto", 1000, -500, 500));
175  hm->registerHisto(
176  new TH1F("hCheck_Hodo_H", "hCheck_Hodo_H", 1000, -500, 500));
177  hm->registerHisto(
178  new TH1F("hCheck_Hodo_V", "hCheck_Hodo_V", 1000, -500, 500));
179 
180  hm->registerHisto(new TH1F("hBeam_Mom", "hBeam_Mom", 1200, -120, 120));
181 
182  hm->registerHisto(new TH2F("hEoP", "hEoP", 1000, 0, 1.5, 120, .5, 120.5));
183 
184  hm->registerHisto(new TH1F("hTemperature", "hTemperature", 500, 0, 50));
185 
186  // help index files with TChain
187  TTree * T = new TTree("T", "T");
188  assert(T);
189  hm->registerHisto(T);
190 
191  T->Branch("info", &_eval_run);
192  T->Branch("clus_3x3_raw", &_eval_3x3_raw);
193  T->Branch("clus_5x5_raw", &_eval_5x5_raw);
194  T->Branch("clus_3x3_prod", &_eval_3x3_prod);
195  T->Branch("clus_5x5_prod", &_eval_5x5_prod);
196  T->Branch("clus_3x3_temp", &_eval_3x3_temp);
197  T->Branch("clus_5x5_temp", &_eval_5x5_temp);
198  T->Branch("clus_3x3_recalib", &_eval_3x3_recalib);
199  T->Branch("clus_5x5_recalib", &_eval_5x5_recalib);
200 
201 
202 
203 
204  //load in hodoscope calibration values
205  ifstream hodoinfile;
206  hodoinfile.open("/sphenix/user/jdosbo/Prototype3/EMCal/JoeShowerCalib/hodoscope_calibs_thirdjointenergyscan.txt",ifstream::in);
207 
208  int horiz=0;
209  if(hodoinfile.is_open()){
210  while (!hodoinfile.eof()){
211  if(horiz==8)
212  break;
213  float hodo0, hodo1, hodo2, hodo3, hodo4, hodo5, hodo6, hodo7;
214  hodoinfile>>hodo0>>hodo1>>hodo2>>hodo3>>hodo4>>hodo5>>hodo6>>hodo7;
215 
216 
217  hodoscope_recalibrations[horiz][0]=hodo0;
218  hodoscope_recalibrations[horiz][1]=hodo1;
219  hodoscope_recalibrations[horiz][2]=hodo2;
220  hodoscope_recalibrations[horiz][3]=hodo3;
221  hodoscope_recalibrations[horiz][4]=hodo4;
222  hodoscope_recalibrations[horiz][5]=hodo5;
223  hodoscope_recalibrations[horiz][6]=hodo6;
224  hodoscope_recalibrations[horiz][7]=hodo7;
225 
226 
227  horiz++;
228  }
229  }
230  if(!hodoinfile.is_open())
231  cout<<"HODO CALIBRATION FILE DIDNT OPEN"<<endl;
232 
233 
234  if(_is_sim){
235 
236  ifstream infile;
237  infile.open("/sphenix/user/jdosbo/Prototype3/EMCal/JoeShowerCalib/3x3clus_posdep_recals_fromsim.txt");
238 
239  int eta = 0;
240  if(infile.is_open()){
241  while(!infile.eof()){
242 
243  if(eta==nfmodbins-1)
244  break;
245 
246  float vals[nfmodbins-1];
247  infile>>vals[0]>>vals[1]>>vals[2]>>vals[3]>>vals[4]>>vals[5]>>vals[6]>>vals[7]
248  >>vals[8]>>vals[9]>>vals[10]>>vals[11]>>vals[12]>>vals[13]>>vals[14]>>vals[15];
249 
250  for(int i=0; i<nfmodbins-1; i++)
251  clus_3x3_sim_recalibs[eta][i] = vals[i];
252 
253  eta++;
254 
255  }
256  }
257  else
258  cout<<"NO SIMULATION POSITION DEPENDENT RECAL OPEN"<<endl;
259 
260  ifstream infile1;
261  infile1.open("/sphenix/user/jdosbo/Prototype3/EMCal/JoeShowerCalib/5x5clus_posdep_recals_fromsim.txt");
262 
263  eta = 0;
264  if(infile1.is_open()){
265  while(!infile1.eof()){
266 
267  if(eta==nfmodbins-1)
268  break;
269 
270  float vals[nfmodbins-1];
271  infile1>>vals[0]>>vals[1]>>vals[2]>>vals[3]>>vals[4]>>vals[5]>>vals[6]>>vals[7]
272  >>vals[8]>>vals[9]>>vals[10]>>vals[11]>>vals[12]>>vals[13]>>vals[14]>>vals[15];
273 
274  for(int i=0; i<nfmodbins-1; i++)
275  clus_5x5_sim_recalibs[eta][i] = vals[i];
276 
277  eta++;
278 
279  }
280  }
281  else
282  cout<<"NO SIMULATION POSITION DEPENDENT RECAL 5X5 OPEN"<<endl;
283 
284  }
285 
286 
287 
288 
289  return Fun4AllReturnCodes::EVENT_OK;
290 }
291 
292 int
293 Proto3ShowerCalib::process_event(PHCompositeNode *topNode)
294 {
295 
296  if (verbosity > 2)
297  cout << "Proto3ShowerCalib::process_event() entered" << endl;
298 
299  // init eval objects
300  _eval_run.reset();
301  _eval_3x3_raw.reset();
302  _eval_5x5_raw.reset();
303  _eval_3x3_prod.reset();
304  _eval_5x5_prod.reset();
305  _eval_3x3_temp.reset();
306  _eval_5x5_temp.reset();
307  _eval_3x3_recalib.reset();
308  _eval_5x5_recalib.reset();
309 
310  Fun4AllHistoManager *hm = get_HistoManager();
311  assert(hm);
312 
313  if (not _is_sim)
314  {
315  PdbParameterMap *info = findNode::getClass<PdbParameterMap>(topNode,
316  "RUN_INFO");
317 
318  assert(info);
319 
320  PHParameters run_info_copy("RunInfo");
321  run_info_copy.FillFrom(info);
322 
323  _eval_run.beam_mom = run_info_copy.get_double_param("beam_MTNRG_GeV");
324 
325  TH1F * hBeam_Mom = dynamic_cast<TH1F *>(hm->getHisto("hBeam_Mom"));
326  assert(hBeam_Mom);
327 
328  hBeam_Mom->Fill(_eval_run.beam_mom);
329 
330  _eval_run.beam_2CH_mm = run_info_copy.get_double_param("beam_2CH_mm");
331  _eval_run.beam_2CV_mm = run_info_copy.get_double_param("beam_2CV_mm");
332 
333  }
334 
335  EventHeader* eventheader = findNode::getClass<EventHeader>(topNode,
336  "EventHeader");
337  if (not _is_sim)
338  {
339  assert(eventheader);
340 
341  _eval_run.run = eventheader->get_RunNumber();
342  if (verbosity > 4)
343  cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
344 
345  _eval_run.event = eventheader->get_EvtSequence();
346  }
347 
348  if (_is_sim)
349  {
350 
351  PHG4TruthInfoContainer* truthInfoList = findNode::getClass<
352  PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
353 
354  assert(truthInfoList);
355 
356  _eval_run.run = -1;
357 
358  const PHG4Particle * p =
359  truthInfoList->GetPrimaryParticleRange().first->second;
360  assert(p);
361 
362  const PHG4VtxPoint * v = truthInfoList->GetVtx(p->get_vtx_id());
363  assert(v);
364 
365  _eval_run.beam_mom = sqrt(
366  p->get_px() * p->get_px() + p->get_py() * p->get_py()
367  + p->get_pz() * p->get_pz());
368  _eval_run.truth_y = v->get_y();
369  _eval_run.truth_z = v->get_z();
370 
371  }
372 
373  // normalization
374  TH1F * hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
375  assert(hNormalization);
376 
377  hNormalization->Fill("ALL", 1);
378 
379  RawTowerContainer* TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
380  topNode, _is_sim ? "TOWER_RAW_LG_CEMC" : "TOWER_RAW_CEMC");
381  assert(TOWER_RAW_CEMC);
382  RawTowerContainer* TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
383  topNode, _is_sim ? "TOWER_CALIB_LG_CEMC" : "TOWER_CALIB_CEMC");
384  assert(TOWER_CALIB_CEMC);
385 
386  // other nodes
387  RawTowerContainer* TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
388  RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
389  if (not _is_sim)
390  {
391  assert(TOWER_CALIB_TRIGGER_VETO);
392  }
393 
394  RawTowerContainer* TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
395  RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
396  if (not _is_sim)
397  {
398  assert(TOWER_CALIB_HODO_HORIZONTAL);
399  }
400  RawTowerContainer* TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
401  RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
402  if (not _is_sim)
403  {
404  assert(TOWER_CALIB_HODO_VERTICAL);
405  }
406 
407  RawTowerContainer* TOWER_TEMPERATURE_EMCAL = findNode::getClass<
408  RawTowerContainer>(topNode, "TOWER_TEMPERATURE_EMCAL");
409  if (not _is_sim)
410  {
411  assert(TOWER_TEMPERATURE_EMCAL);
412  }
413 
414  RawTowerContainer* TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
415  topNode, "TOWER_CALIB_C1");
416  if (not _is_sim)
417  {
418  assert(TOWER_CALIB_C1);
419  }
420  RawTowerContainer* TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
421  topNode, "TOWER_CALIB_C2");
422  if (not _is_sim)
423  {
424  assert(TOWER_CALIB_C2);
425  }
426 
427 
428  RawTowerGeomContainer_Cylinderv1 *towergeom = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode,"TOWERGEOM_CEMC");
429 
430  if( _is_sim)
431  {
432  assert(towergeom);
433  }
434 
435  // Cherenkov
436  bool cherekov_e = false;
437  if (not _is_sim)
438  {
439  RawTower * t_c2_in = NULL;
440  RawTower * t_c2_out = NULL;
441 
442  assert(eventheader);
443  if (eventheader->get_RunNumber() >= 2105)
444  {
445  t_c2_in = TOWER_CALIB_C2->getTower(10);
446  t_c2_out = TOWER_CALIB_C2->getTower(11);
447  }
448  else
449  {
450  t_c2_in = TOWER_CALIB_C2->getTower(0);
451  t_c2_out = TOWER_CALIB_C2->getTower(1);
452  }
453  assert(t_c2_in);
454  assert(t_c2_out);
455 
456  const double c2_in = t_c2_in->get_energy();
457  const double c2_out = t_c2_out->get_energy();
458  const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
459 
460  _eval_run.C2_sum = c2_in + c2_out;
461  cherekov_e = (_eval_run.C2_sum)
462  > (abs(_eval_run.beam_mom) >= 10 ? 100 : 240);
463  hNormalization->Fill("C2-e", cherekov_e);
464 
465  TH2F * hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
466  "hCheck_Cherenkov"));
467  assert(hCheck_Cherenkov);
468  hCheck_Cherenkov->Fill(c1, "C1", 1);
469  hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
470  hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
471  hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
472  }
473 
474  // veto
475  TH1F * hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
476  assert(hCheck_Veto);
477  bool trigger_veto_pass = true;
478  if (not _is_sim)
479  {
480  auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
481  for (auto it = range.first; it != range.second; ++it)
482  {
483  RawTower* tower = it->second;
484  assert(tower);
485 
486  hCheck_Veto->Fill(tower->get_energy());
487 
488  if (abs(tower->get_energy()) > 15)
489  trigger_veto_pass = false;
490  }
491  }
492  hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
493  _eval_run.trigger_veto_pass = trigger_veto_pass;
494 
495  // hodoscope
496  TH1F * hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
497  assert(hCheck_Hodo_H);
498  int hodo_h_count = 0;
499  if (not _is_sim)
500  {
501  auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
502  for (auto it = range.first; it != range.second; ++it)
503  {
504  RawTower* tower = it->second;
505  assert(tower);
506 
507  hCheck_Hodo_H->Fill(tower->get_energy());
508 
509  if (abs(tower->get_energy()) > 30)
510  {
511  hodo_h_count++;
512  _eval_run.hodo_h = tower->get_id();
513  }
514  }
515  }
516  const bool valid_hodo_h = hodo_h_count == 1;
517  hNormalization->Fill("valid_hodo_h", valid_hodo_h);
518  _eval_run.valid_hodo_h = valid_hodo_h;
519 
520  TH1F * hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
521  assert(hCheck_Hodo_V);
522  int hodo_v_count = 0;
523  if (not _is_sim)
524  {
525  auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
526  for (auto it = range.first; it != range.second; ++it)
527  {
528  RawTower* tower = it->second;
529  assert(tower);
530 
531  hCheck_Hodo_V->Fill(tower->get_energy());
532 
533  if (abs(tower->get_energy()) > 30)
534  {
535  hodo_v_count++;
536  _eval_run.hodo_v = tower->get_id();
537  }
538  }
539  }
540  const bool valid_hodo_v = hodo_v_count == 1;
541  _eval_run.valid_hodo_v = valid_hodo_v;
542  hNormalization->Fill("valid_hodo_v", valid_hodo_v);
543 
544  const bool good_e = (valid_hodo_v and valid_hodo_h and cherekov_e
545  and trigger_veto_pass) and (not _is_sim);
546  hNormalization->Fill("good_e", good_e);
547 
548  // simple clustering
549  pair<int, int> max_3x3 = find_max(TOWER_CALIB_CEMC, 3);
550  pair<int, int> max_5x5 = find_max(TOWER_CALIB_CEMC, 5);
551 
552  _eval_3x3_raw.max_col = max_3x3.first;
553  _eval_3x3_raw.max_row = max_3x3.second;
554  _eval_3x3_prod.max_col = max_3x3.first;
555  _eval_3x3_prod.max_row = max_3x3.second;
556  _eval_3x3_temp.max_col = max_3x3.first;
557  _eval_3x3_temp.max_row = max_3x3.second;
558  _eval_3x3_recalib.max_col = max_3x3.first;
559  _eval_3x3_recalib.max_row = max_3x3.second;
560 
561  _eval_5x5_raw.max_col = max_5x5.first;
562  _eval_5x5_raw.max_row = max_5x5.second;
563  _eval_5x5_prod.max_col = max_5x5.first;
564  _eval_5x5_prod.max_row = max_5x5.second;
565  _eval_5x5_temp.max_col = max_5x5.first;
566  _eval_5x5_temp.max_row = max_5x5.second;
567  _eval_5x5_recalib.max_col = max_5x5.first;
568  _eval_5x5_recalib.max_row = max_5x5.second;
569 
570  // tower
571  bool good_temp = true;
572  double sum_energy_calib = 0;
573  double sum_energy_T = 0;
574  TH1F * hTemperature = dynamic_cast<TH1F *>(hm->getHisto("hTemperature"));
575  assert(hTemperature);
576 
577  stringstream sdata;
578 
579 
580  //code for position dependence in simulation
581  std::vector<float> toweretas3x3;
582  std::vector<float> towerphis3x3;
583  std::vector<float> towerenergies3x3;
584  std::vector<float> toweretas5x5;
585  std::vector<float> towerphis5x5;
586  std::vector<float> towerenergies5x5;
587 
588 
589 
590  if (good_e)
591  sdata << abs(_eval_run.beam_mom) << "\t";
592 
593  // tower temperature and recalibration
594  {
595  auto range = TOWER_CALIB_CEMC->getTowers();
596  for (auto it = range.first; it != range.second; ++it)
597  {
598 
599  RawTowerDefs::keytype key = it->first;
600  RawTower* tower = it->second;
601  assert(tower);
602 
603  const int col = tower->get_bineta();
604  const int row = tower->get_binphi();
605  const int hodorow = _eval_run.hodo_h;
606  const int hodocol = _eval_run.hodo_v;
607 
608  if (col < 0 or col >= 8)
609  continue;
610 
611  if (row < 0 or row >= 8)
612  continue;
613 
614  if((hodorow<0 or hodorow>7) and !_is_sim)
615  continue;
616  if((hodocol<0 or hodocol>7) and !_is_sim)
617  continue;
618 
619  const double energy_calib = tower->get_energy();
620  sum_energy_calib += energy_calib;
621 
622  RawTower* tower_raw = TOWER_RAW_CEMC->getTower(key);
623  assert(tower_raw);
624 
625  double energy_T = energy_calib;
626 
627 // if (not _is_sim)
628 // {
629 // RawTower_Temperature * temp_t =
630 // dynamic_cast<RawTower_Temperature *>(TOWER_TEMPERATURE_EMCAL->getTower(
631 // tower->get_row(), tower->get_column())); // note swap of col/row in temperature storage
632 // assert(temp_t);
633 //
634 // const double T = temp_t->get_temperature_from_time(
635 // eventheader->get_TimeStamp());
636 // hTemperature->Fill(T);
637 //
638 // if (T < 25 or T > 35)
639 // good_temp = false;
640 //
641 // energy_T = TemperatureCorrection::Apply(energy_calib, T);
642 // }
643 
644  // recalibration
645  assert(
646  _recalib_const.find(make_pair(col, row)) != _recalib_const.end());
647  double energy_recalib=0;
648 
649  if(_use_hodoscope_calibs)
650  energy_recalib = energy_T
651  *hodoscope_recalibrations[hodorow][7-hodocol];
652 
653 
654 
655 
656 
657  // energy sums
658  sum_energy_T += energy_T;
659 
660  // calibration file
661 // sdata << tower->get_energy() << "\t";
662  // calibration file - only output 5x5 towers
663  if (col >= max_5x5.first - 1 and col <= max_5x5.first + 1
664  and row >= max_5x5.second - 1 and row <= max_5x5.second + 1)
665  {
666  sdata << tower->get_energy() << "\t";
667  }
668  else
669  {
670  sdata << 0 << "\t";
671  }
672 
673 
674 
675  // cluster 3x3
676  if (col >= max_3x3.first - 1 and col <= max_3x3.first + 1)
677  if (row >= max_3x3.second - 1 and row <= max_3x3.second + 1)
678  {
679  // in cluster
680 
681  _eval_3x3_raw.average_col += abs(tower_raw->get_energy()) * col;
682  _eval_3x3_raw.average_row += abs(tower_raw->get_energy()) * row;
683  _eval_3x3_raw.sum_E += abs(tower_raw->get_energy());
684 
685  _eval_3x3_prod.average_col += energy_calib * col;
686  _eval_3x3_prod.average_row += energy_calib * row;
687  _eval_3x3_prod.sum_E += energy_calib;
688 
689  _eval_3x3_temp.average_col += energy_T * col;
690  _eval_3x3_temp.average_row += energy_T * row;
691  _eval_3x3_temp.sum_E += energy_T;
692 
693  _eval_3x3_recalib.average_col += energy_recalib * col;
694  _eval_3x3_recalib.average_row += energy_recalib * row;
695  _eval_3x3_recalib.sum_E += energy_recalib;
696 
697 
698  toweretas3x3.push_back(col);
699  towerphis3x3.push_back(row);
700  towerenergies3x3.push_back(energy_calib);
701 
702 
703 
704  }
705 
706  // cluster 5x5
707  if (col >= max_5x5.first - 2 and col <= max_5x5.first + 2)
708  if (row >= max_5x5.second - 2 and row <= max_5x5.second + 2)
709  {
710  // in cluster
711 
712  _eval_5x5_raw.average_col += abs(tower_raw->get_energy()) * col;
713  _eval_5x5_raw.average_row += abs(tower_raw->get_energy()) * row;
714  _eval_5x5_raw.sum_E += abs(tower_raw->get_energy());
715 
716  _eval_5x5_prod.average_col += energy_calib * col;
717  _eval_5x5_prod.average_row += energy_calib * row;
718  _eval_5x5_prod.sum_E += energy_calib;
719 
720  _eval_5x5_temp.average_col += energy_T * col;
721  _eval_5x5_temp.average_row += energy_T * row;
722  _eval_5x5_temp.sum_E += energy_T;
723 
724  _eval_5x5_recalib.average_col += energy_recalib * col;
725  _eval_5x5_recalib.average_row += energy_recalib * row;
726  _eval_5x5_recalib.sum_E += energy_recalib;
727 
728  toweretas5x5.push_back(col);
729  towerphis5x5.push_back(row);
730  towerenergies5x5.push_back(energy_calib);
731 
732  }
733  }
734  }
735 
736  _eval_3x3_raw.reweight_clus_pol();
737  _eval_5x5_raw.reweight_clus_pol();
738  _eval_3x3_prod.reweight_clus_pol();
739  _eval_5x5_prod.reweight_clus_pol();
740  _eval_3x3_temp.reweight_clus_pol();
741  _eval_5x5_temp.reweight_clus_pol();
742  _eval_3x3_recalib.reweight_clus_pol();
743  _eval_5x5_recalib.reweight_clus_pol();
744 
745 
746 
747 
748  //enter simulation position dependent weighting code
749 
750 
751  //do 3x3 clusters first
752  int ntowers = toweretas3x3.size();
753  //assert(ntowers >= 1);
754  const int nphibin = towergeom->get_phibins();
755  if(ntowers>0){
756 
757  //loop over the towers to determine the energy
758  //weighted eta and phi position of the cluster
759 
760  float etamult = 0;
761  float etasum = 0;
762  float phimult = 0;
763  float phisum = 0;
764 
765  for (int j = 0; j < ntowers; j++)
766  {
767  float energymult = towerenergies3x3.at(j) * toweretas3x3.at(j);
768  etamult += energymult;
769  etasum += towerenergies3x3.at(j);
770 
771  int phibin = towerphis3x3.at(j);
772 
773  if (phibin - towerphis3x3.at(0) < -nphibin / 2)
774  phibin += nphibin;
775  else if (phibin - towerphis3x3.at(0) > +nphibin / 2)
776  phibin -= nphibin;
777  assert(abs(phibin - towerphis3x3.at(0)) <= nphibin / 2);
778 
779  energymult = towerenergies3x3.at(j) * phibin;
780  phimult += energymult;
781  phisum += towerenergies3x3.at(j);
782  }
783 
784 
785  float avgphi = phimult / phisum;
786  float avgeta = etamult / etasum;
787 
788  if (avgphi<0) avgphi += nphibin;
789 
790  //this determines the position of the cluster in the 2x2 block
791  float fmodphi3x3 = fmod(avgphi, 2.);
792  float fmodeta3x3 = fmod(avgeta, 2.);
793 
794  _eval_3x3_prod.setfmods(fmodeta3x3, fmodphi3x3);
795 
796 
797  int fmodetabin = -99;
798  int fmodphibin = -99;
799  for(int k=0; k<nfmodbins-1; k++)
800  if(fmodeta3x3 >= fmodbins[k] && fmodeta3x3 < fmodbins[k+1])
801  fmodetabin = k;
802  for(int k=0; k<nfmodbins-1; k++)
803  if(fmodphi3x3 >= fmodbins[k] && fmodphi3x3 < fmodbins[k+1])
804  fmodphibin = k;
805 
806  if(fmodetabin>-1 && fmodphibin>-1 && !_use_hodoscope_calibs)
807  _eval_3x3_recalib.sum_E = _eval_3x3_prod.getsumE()/clus_3x3_sim_recalibs[fmodetabin][fmodphibin];
808  _eval_3x3_recalib.setfmods(fmodeta3x3,fmodphi3x3);
809 
810 
811  }
812 
813  //NOW DO THE 5X5 CLUSTERS
814 
815  ntowers = toweretas5x5.size();
816  //assert(ntowers >= 1);
817  if(ntowers>0){
818  float etamult = 0;
819  float etasum = 0;
820  float phimult = 0;
821  float phisum = 0;
822  for (int j = 0; j < ntowers; j++)
823  {
824  float energymult = towerenergies5x5.at(j) * toweretas5x5.at(j);
825  etamult += energymult;
826  etasum += towerenergies5x5.at(j);
827 
828  int phibin = towerphis5x5.at(j);
829 
830  if (phibin - towerphis5x5.at(0) < -nphibin / 2)
831  phibin += nphibin;
832  else if (phibin - towerphis5x5.at(0) > +nphibin / 2)
833  phibin -= nphibin;
834  assert(abs(phibin - towerphis5x5.at(0)) <= nphibin / 2);
835 
836  energymult = towerenergies5x5.at(j) * phibin;
837  phimult += energymult;
838  phisum += towerenergies5x5.at(j);
839  }
840 
841  float avgphi = phimult / phisum;
842  float avgeta = etamult / etasum;
843 
844  if (avgphi < 0) avgphi += nphibin;
845 
846  float fmodphi5x5 = fmod(avgphi , 2.);
847  float fmodeta5x5 = fmod(avgeta , 2.);
848 
849  _eval_5x5_prod.setfmods(fmodeta5x5, fmodphi5x5);
850 
851 
852  int fmodetabin = -99;
853  int fmodphibin = -99;
854  for(int k=0; k<nfmodbins-1; k++)
855  if(fmodeta5x5 >= fmodbins[k] && fmodeta5x5 < fmodbins[k+1])
856  fmodetabin = k;
857  for(int k=0; k<nfmodbins-1; k++)
858  if(fmodphi5x5 >= fmodbins[k] && fmodphi5x5 < fmodbins[k+1])
859  fmodphibin = k;
860 
861  if(fmodetabin > -1 && fmodphibin > -1 && !_use_hodoscope_calibs)
862  _eval_5x5_recalib.sum_E = _eval_5x5_prod.getsumE()/clus_5x5_sim_recalibs[fmodetabin][fmodphibin];
863 
864  _eval_5x5_recalib.setfmods(fmodeta5x5,fmodphi5x5);
865  }
866 
867 
868 
869 
870 
871 
872 
873 
874 
875 
876 
877 
878 
879  const double EoP = sum_energy_T / abs(_eval_run.beam_mom);
880  hNormalization->Fill("good_temp", good_temp);
881 
882 // bool good_data = good_e and good_temp;
883  bool good_data = good_e;
884  hNormalization->Fill("good_data", good_data);
885 
886  _eval_run.good_temp = good_temp;
887  _eval_run.good_e = good_e;
888  _eval_run.good_data = good_data;
889  _eval_run.sum_energy_T = sum_energy_T;
890  _eval_run.EoP = EoP;
891 
892  // E/p
893  if (good_data)
894  {
895  if (verbosity >= 3)
896  cout << __PRETTY_FUNCTION__ << " sum_energy_calib = "
897  << sum_energy_calib << " sum_energy_T = " << sum_energy_T
898  << " _eval_run.beam_mom = " << _eval_run.beam_mom << endl;
899 
900  TH2F * hEoP = dynamic_cast<TH2F *>(hm->getHisto("hEoP"));
901  assert(hEoP);
902 
903  hEoP->Fill(EoP, abs(_eval_run.beam_mom));
904  }
905 
906  // calibration file
907  if (good_data and abs(_eval_run.beam_mom) >= 4
908  and abs(_eval_run.beam_mom) <= 8
909  and abs(_eval_3x3_raw.average_col - round(_eval_3x3_raw.average_col))
910  < 0.1
911  and abs(_eval_3x3_raw.average_row - round(_eval_3x3_raw.average_row))
912  < 0.1)
913  {
914  assert(fdata.is_open());
915 
916  fdata << sdata.str();
917 
918  fdata << endl;
919  }
920 
921  TTree * T = dynamic_cast<TTree *>(hm->getHisto("T"));
922  assert(T);
923  T->Fill();
924 
925  return Fun4AllReturnCodes::EVENT_OK;
926 }
927 
928 pair<int, int>
929 Proto3ShowerCalib::find_max(RawTowerContainer* towers, int cluster_size)
930 {
931  const int clus_edge_size = (cluster_size - 1) / 2;
932  assert(clus_edge_size >= 0);
933 
934  pair<int, int> max(-1000, -1000);
935  double max_e = 0;
936 
937  for (int col = 0; col < n_size; ++col)
938  for (int row = 0; row < n_size; ++row)
939  {
940  double energy = 0;
941 
942  for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
943  ++dcol)
944  for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
945  ++drow)
946  {
947  if (dcol < 0 or drow < 0)
948  continue;
949 
950  RawTower * t = towers->getTower(dcol, drow);
951  if (t)
952  energy += t->get_energy();
953  }
954 
955  if (energy > max_e)
956  {
957  max = make_pair(col, row);
958  max_e = energy;
959  }
960  }
961 
962  return max;
963 }
964 
965 int
966 Proto3ShowerCalib::LoadRecalibMap(const std::string & file)
967 {
968  if (verbosity)
969  {
970  cout << __PRETTY_FUNCTION__ << " - input recalibration constant from "
971  << file << endl;
972  }
973 
974  ifstream fcalib(file);
975 
976  assert(fcalib.is_open());
977 
978  string line;
979 
980  while (not fcalib.eof())
981  {
982  getline(fcalib, line);
983 
984  if (verbosity > 10)
985  {
986  cout << __PRETTY_FUNCTION__ << " get line " << line << endl;
987  }
988  istringstream sline(line);
989 
990  int col = -1;
991  int row = -1;
992  double calib = 0;
993 
994  sline >> col >> row >> calib;
995 
996  if (not sline.fail())
997  {
998  if (verbosity)
999  {
1000  cout << __PRETTY_FUNCTION__ << " - recalibration constant "
1001  << col << "," << row << " = " << calib << endl;
1002  }
1003 
1004  _recalib_const[make_pair(col, row)] = calib;
1005  }
1006 
1007  }
1008 
1009  return _recalib_const.size();
1010 }
1011