EMCal/ShowerCalib_PositionDependent/Proto3ShowerCalib.C

0001 #include "Proto3ShowerCalib.h"
0002 #include "TemperatureCorrection.h"
0003
0004 #include <prototype3/RawTower_Temperature.h>
0005 #include <prototype3/RawTower_Prototype3.h>
0006 #include <calobase/RawTowerContainer.h>
0007 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0008 #include <pdbcalbase/PdbParameterMap.h>
0009 #include <phparameter/PHParameters.h>
0010 #include <ffaobjects/EventHeader.h>
0011
0012 #include <fun4all/SubsysReco.h>
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/PHTFileServer.h>
0015 #include <phool/PHCompositeNode.h>
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 #include <phool/getClass.h>
0018 #include <fun4all/Fun4AllHistoManager.h>
0019
0020 #include <phool/PHCompositeNode.h>
0021
0022 #include <g4main/PHG4TruthInfoContainer.h>
0023 #include <g4main/PHG4Particle.h>
0024 #include <g4main/PHG4VtxPoint.h>
0025
0026 #include <TNtuple.h>
0027 #include <TFile.h>
0028 #include <TH1F.h>
0029 #include <TH2F.h>
0030 #include <TVector3.h>
0031 #include <TLorentzVector.h>
0032
0033 #include <exception>
0034 #include <stdexcept>
0035 #include <iostream>
0036 #include <sstream>
0037 #include <vector>
0038 #include <set>
0039 #include <algorithm>
0040 #include <cassert>
0041 #include <cmath>
0042
0043 using namespace std;
0044
0045 ClassImp(Proto3ShowerCalib::Eval_Cluster);
0046 ClassImp(Proto3ShowerCalib::Eval_Run);
0047
0048 Proto3ShowerCalib::Proto3ShowerCalib(const std::string &filename) :
0049   SubsysReco("Proto3ShowerCalib"), _is_sim(false), _use_hodoscope_calibs(false),_filename(filename), _ievent(
0050         0)
0051 {
0052
0053   verbosity = 1;
0054
0055   _eval_run.reset();
0056   _eval_3x3_raw.reset();
0057   _eval_5x5_raw.reset();
0058   _eval_3x3_prod.reset();
0059   _eval_5x5_prod.reset();
0060   _eval_3x3_temp.reset();
0061   _eval_5x5_temp.reset();
0062   _eval_3x3_recalib.reset();
0063   _eval_5x5_recalib.reset();
0064
0065   for (int col = 0; col < n_size; ++col)
0066     for (int row = 0; row < n_size; ++row)
0067       {
0068         _recalib_const[make_pair(col, row)] = 0;
0069       }
0070
0071 }
0072
0073 Proto3ShowerCalib::~Proto3ShowerCalib()
0074 {
0075 }
0076
0077 Fun4AllHistoManager *
0078 Proto3ShowerCalib::get_HistoManager()
0079 {
0080
0081   Fun4AllServer *se = Fun4AllServer::instance();
0082   Fun4AllHistoManager *hm = se->getHistoManager("Proto3ShowerCalib_HISTOS");
0083
0084   if (not hm)
0085     {
0086       cout
0087           << "Proto3ShowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto3ShowerCalib_HISTOS"
0088           << endl;
0089       hm = new Fun4AllHistoManager("Proto3ShowerCalib_HISTOS");
0090       se->registerHistoManager(hm);
0091     }
0092
0093   assert(hm);
0094
0095   return hm;
0096 }
0097
0098 int
0099 Proto3ShowerCalib::InitRun(PHCompositeNode *topNode)
0100 {
0101   if (verbosity)
0102     cout << "Proto3ShowerCalib::InitRun" << endl;
0103
0104   _ievent = 0;
0105   cout<<"IS THIS SIMULATION: "<<_is_sim<<endl;
0106   PHNodeIterator iter(topNode);
0107   PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
0108       "PHCompositeNode", "DST"));
0109   if (!dstNode)
0110     {
0111       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0112       throw runtime_error(
0113           "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0114     }
0115   for(int i=0; i<nfmodbins; i++)
0116     fmodbins[i] = 0. + i * 2. / (float)(nfmodbins-1);
0117
0118   return Fun4AllReturnCodes::EVENT_OK;
0119 }
0120
0121 int
0122 Proto3ShowerCalib::End(PHCompositeNode *topNode)
0123 {
0124   cout << "Proto3ShowerCalib::End - write to " << _filename << endl;
0125   PHTFileServer::get().cd(_filename);
0126
0127   Fun4AllHistoManager *hm = get_HistoManager();
0128   assert(hm);
0129   for (unsigned int i = 0; i < hm->nHistos(); i++)
0130     hm->getHisto(i)->Write();
0131
0132 //  if (_T_EMCalTrk)
0133 //    _T_EMCalTrk->Write();
0134
0135   fdata.close();
0136
0137   return Fun4AllReturnCodes::EVENT_OK;
0138 }
0139
0140 int
0141 Proto3ShowerCalib::Init(PHCompositeNode *topNode)
0142 {
0143
0144   _ievent = 0;
0145
0146   cout << "Proto3ShowerCalib::get_HistoManager - Making PHTFileServer "
0147       << _filename << endl;
0148   PHTFileServer::get().open(_filename, "RECREATE");
0149
0150   fdata.open(_filename + ".dat", std::fstream::out);
0151
0152   Fun4AllHistoManager *hm = get_HistoManager();
0153   assert(hm);
0154
0155   TH2F * hCheck_Cherenkov = new TH2F("hCheck_Cherenkov", "hCheck_Cherenkov",
0156       1000, -2000, 2000, 5, .5, 5.5);
0157   hCheck_Cherenkov->GetYaxis()->SetBinLabel(1, "C1");
0158   hCheck_Cherenkov->GetYaxis()->SetBinLabel(2, "C2 in");
0159   hCheck_Cherenkov->GetYaxis()->SetBinLabel(3, "C2 out");
0160   hCheck_Cherenkov->GetYaxis()->SetBinLabel(4, "C2 sum");
0161   hm->registerHisto(hCheck_Cherenkov);
0162
0163   TH1F * hNormalization = new TH1F("hNormalization", "hNormalization", 10, .5,
0164       10.5);
0165   hCheck_Cherenkov->GetXaxis()->SetBinLabel(1, "ALL");
0166   hCheck_Cherenkov->GetXaxis()->SetBinLabel(2, "C2-e");
0167   hCheck_Cherenkov->GetXaxis()->SetBinLabel(3, "trigger_veto_pass");
0168   hCheck_Cherenkov->GetXaxis()->SetBinLabel(4, "valid_hodo_h");
0169   hCheck_Cherenkov->GetXaxis()->SetBinLabel(5, "valid_hodo_v");
0170   hCheck_Cherenkov->GetXaxis()->SetBinLabel(6, "good_e");
0171   hCheck_Cherenkov->GetXaxis()->SetBinLabel(7, "good_data");
0172   hm->registerHisto(hNormalization);
0173
0174   hm->registerHisto(new TH1F("hCheck_Veto", "hCheck_Veto", 1000, -500, 500));
0175   hm->registerHisto(
0176       new TH1F("hCheck_Hodo_H", "hCheck_Hodo_H", 1000, -500, 500));
0177   hm->registerHisto(
0178       new TH1F("hCheck_Hodo_V", "hCheck_Hodo_V", 1000, -500, 500));
0179
0180   hm->registerHisto(new TH1F("hBeam_Mom", "hBeam_Mom", 1200, -120, 120));
0181
0182   hm->registerHisto(new TH2F("hEoP", "hEoP", 1000, 0, 1.5, 120, .5, 120.5));
0183
0184   hm->registerHisto(new TH1F("hTemperature", "hTemperature", 500, 0, 50));
0185
0186   // help index files with TChain
0187   TTree * T = new TTree("T", "T");
0188   assert(T);
0189   hm->registerHisto(T);
0190
0191   T->Branch("info", &_eval_run);
0192   T->Branch("clus_3x3_raw", &_eval_3x3_raw);
0193   T->Branch("clus_5x5_raw", &_eval_5x5_raw);
0194   T->Branch("clus_3x3_prod", &_eval_3x3_prod);
0195   T->Branch("clus_5x5_prod", &_eval_5x5_prod);
0196   T->Branch("clus_3x3_temp", &_eval_3x3_temp);
0197   T->Branch("clus_5x5_temp", &_eval_5x5_temp);
0198   T->Branch("clus_3x3_recalib", &_eval_3x3_recalib);
0199   T->Branch("clus_5x5_recalib", &_eval_5x5_recalib);
0200
0201
0202
0203
0204   //load in hodoscope calibration values
0205   ifstream hodoinfile;
0206   hodoinfile.open("/sphenix/user/jdosbo/Prototype3/EMCal/JoeShowerCalib/hodoscope_calibs_thirdjointenergyscan.txt",ifstream::in);
0207
0208   int horiz=0;
0209   if(hodoinfile.is_open()){
0210     while (!hodoinfile.eof()){
0211       if(horiz==8)
0212     break;
0213       float hodo0, hodo1, hodo2, hodo3, hodo4, hodo5, hodo6, hodo7;
0214       hodoinfile>>hodo0>>hodo1>>hodo2>>hodo3>>hodo4>>hodo5>>hodo6>>hodo7;
0215
0216
0217       hodoscope_recalibrations[horiz][0]=hodo0;
0218       hodoscope_recalibrations[horiz][1]=hodo1;
0219       hodoscope_recalibrations[horiz][2]=hodo2;
0220       hodoscope_recalibrations[horiz][3]=hodo3;
0221       hodoscope_recalibrations[horiz][4]=hodo4;
0222       hodoscope_recalibrations[horiz][5]=hodo5;
0223       hodoscope_recalibrations[horiz][6]=hodo6;
0224       hodoscope_recalibrations[horiz][7]=hodo7;
0225
0226
0227       horiz++;
0228     }
0229   }
0230   if(!hodoinfile.is_open())
0231     cout<<"HODO CALIBRATION FILE DIDNT OPEN"<<endl;
0232
0233
0234   if(_is_sim){
0235
0236     ifstream infile;
0237     infile.open("/sphenix/user/jdosbo/Prototype3/EMCal/JoeShowerCalib/3x3clus_posdep_recals_fromsim.txt");
0238
0239     int eta = 0;
0240     if(infile.is_open()){
0241       while(!infile.eof()){
0242
0243     if(eta==nfmodbins-1)
0244       break;
0245
0246     float vals[nfmodbins-1];
0247     infile>>vals[0]>>vals[1]>>vals[2]>>vals[3]>>vals[4]>>vals[5]>>vals[6]>>vals[7]
0248           >>vals[8]>>vals[9]>>vals[10]>>vals[11]>>vals[12]>>vals[13]>>vals[14]>>vals[15];
0249
0250     for(int i=0; i<nfmodbins-1; i++)
0251       clus_3x3_sim_recalibs[eta][i] = vals[i];
0252
0253     eta++;
0254
0255       }
0256     }
0257     else
0258       cout<<"NO SIMULATION POSITION DEPENDENT RECAL OPEN"<<endl;
0259
0260     ifstream infile1;
0261     infile1.open("/sphenix/user/jdosbo/Prototype3/EMCal/JoeShowerCalib/5x5clus_posdep_recals_fromsim.txt");
0262
0263     eta = 0;
0264     if(infile1.is_open()){
0265       while(!infile1.eof()){
0266
0267     if(eta==nfmodbins-1)
0268       break;
0269
0270     float vals[nfmodbins-1];
0271     infile1>>vals[0]>>vals[1]>>vals[2]>>vals[3]>>vals[4]>>vals[5]>>vals[6]>>vals[7]
0272            >>vals[8]>>vals[9]>>vals[10]>>vals[11]>>vals[12]>>vals[13]>>vals[14]>>vals[15];
0273
0274     for(int i=0; i<nfmodbins-1; i++)
0275       clus_5x5_sim_recalibs[eta][i] = vals[i];
0276
0277     eta++;
0278
0279       }
0280     }
0281     else
0282       cout<<"NO SIMULATION POSITION DEPENDENT RECAL 5X5 OPEN"<<endl;
0283
0284   }
0285
0286
0287
0288
0289   return Fun4AllReturnCodes::EVENT_OK;
0290 }
0291
0292 int
0293 Proto3ShowerCalib::process_event(PHCompositeNode *topNode)
0294 {
0295
0296   if (verbosity > 2)
0297     cout << "Proto3ShowerCalib::process_event() entered" << endl;
0298
0299   // init eval objects
0300   _eval_run.reset();
0301   _eval_3x3_raw.reset();
0302   _eval_5x5_raw.reset();
0303   _eval_3x3_prod.reset();
0304   _eval_5x5_prod.reset();
0305   _eval_3x3_temp.reset();
0306   _eval_5x5_temp.reset();
0307   _eval_3x3_recalib.reset();
0308   _eval_5x5_recalib.reset();
0309
0310   Fun4AllHistoManager *hm = get_HistoManager();
0311   assert(hm);
0312
0313   if (not _is_sim)
0314     {
0315       PdbParameterMap *info = findNode::getClass<PdbParameterMap>(topNode,
0316           "RUN_INFO");
0317
0318       assert(info);
0319
0320       PHParameters run_info_copy("RunInfo");
0321       run_info_copy.FillFrom(info);
0322
0323       _eval_run.beam_mom = run_info_copy.get_double_param("beam_MTNRG_GeV");
0324
0325       TH1F * hBeam_Mom = dynamic_cast<TH1F *>(hm->getHisto("hBeam_Mom"));
0326       assert(hBeam_Mom);
0327
0328       hBeam_Mom->Fill(_eval_run.beam_mom);
0329
0330       _eval_run.beam_2CH_mm = run_info_copy.get_double_param("beam_2CH_mm");
0331       _eval_run.beam_2CV_mm = run_info_copy.get_double_param("beam_2CV_mm");
0332
0333     }
0334
0335   EventHeader* eventheader = findNode::getClass<EventHeader>(topNode,
0336       "EventHeader");
0337   if (not _is_sim)
0338     {
0339       assert(eventheader);
0340
0341       _eval_run.run = eventheader->get_RunNumber();
0342       if (verbosity > 4)
0343         cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
0344
0345       _eval_run.event = eventheader->get_EvtSequence();
0346     }
0347
0348   if (_is_sim)
0349     {
0350
0351       PHG4TruthInfoContainer* truthInfoList = findNode::getClass<
0352           PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0353
0354       assert(truthInfoList);
0355
0356       _eval_run.run = -1;
0357
0358       const PHG4Particle * p =
0359           truthInfoList->GetPrimaryParticleRange().first->second;
0360       assert(p);
0361
0362       const PHG4VtxPoint * v = truthInfoList->GetVtx(p->get_vtx_id());
0363       assert(v);
0364
0365       _eval_run.beam_mom = sqrt(
0366           p->get_px() * p->get_px() + p->get_py() * p->get_py()
0367               + p->get_pz() * p->get_pz());
0368       _eval_run.truth_y = v->get_y();
0369       _eval_run.truth_z = v->get_z();
0370
0371     }
0372
0373   // normalization
0374   TH1F * hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
0375   assert(hNormalization);
0376
0377   hNormalization->Fill("ALL", 1);
0378
0379   RawTowerContainer* TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
0380       topNode, _is_sim ? "TOWER_RAW_LG_CEMC" : "TOWER_RAW_CEMC");
0381   assert(TOWER_RAW_CEMC);
0382   RawTowerContainer* TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
0383       topNode, _is_sim ? "TOWER_CALIB_LG_CEMC" : "TOWER_CALIB_CEMC");
0384   assert(TOWER_CALIB_CEMC);
0385
0386   // other nodes
0387   RawTowerContainer* TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
0388       RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
0389   if (not _is_sim)
0390     {
0391       assert(TOWER_CALIB_TRIGGER_VETO);
0392     }
0393
0394   RawTowerContainer* TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
0395       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
0396   if (not _is_sim)
0397     {
0398       assert(TOWER_CALIB_HODO_HORIZONTAL);
0399     }
0400   RawTowerContainer* TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
0401       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
0402   if (not _is_sim)
0403     {
0404       assert(TOWER_CALIB_HODO_VERTICAL);
0405     }
0406
0407   RawTowerContainer* TOWER_TEMPERATURE_EMCAL = findNode::getClass<
0408       RawTowerContainer>(topNode, "TOWER_TEMPERATURE_EMCAL");
0409   if (not _is_sim)
0410     {
0411       assert(TOWER_TEMPERATURE_EMCAL);
0412     }
0413
0414   RawTowerContainer* TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
0415       topNode, "TOWER_CALIB_C1");
0416   if (not _is_sim)
0417     {
0418       assert(TOWER_CALIB_C1);
0419     }
0420   RawTowerContainer* TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
0421       topNode, "TOWER_CALIB_C2");
0422   if (not _is_sim)
0423     {
0424       assert(TOWER_CALIB_C2);
0425     }
0426
0427
0428   RawTowerGeomContainer_Cylinderv1 *towergeom = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode,"TOWERGEOM_CEMC");
0429
0430   if( _is_sim)
0431     {
0432       assert(towergeom);
0433     }
0434
0435   // Cherenkov
0436   bool cherekov_e = false;
0437   if (not _is_sim)
0438     {
0439       RawTower * t_c2_in = NULL;
0440       RawTower * t_c2_out = NULL;
0441
0442       assert(eventheader);
0443       if (eventheader->get_RunNumber() >= 2105)
0444         {
0445           t_c2_in = TOWER_CALIB_C2->getTower(10);
0446           t_c2_out = TOWER_CALIB_C2->getTower(11);
0447         }
0448       else
0449         {
0450           t_c2_in = TOWER_CALIB_C2->getTower(0);
0451           t_c2_out = TOWER_CALIB_C2->getTower(1);
0452         }
0453       assert(t_c2_in);
0454       assert(t_c2_out);
0455
0456       const double c2_in = t_c2_in->get_energy();
0457       const double c2_out = t_c2_out->get_energy();
0458       const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
0459
0460       _eval_run.C2_sum = c2_in + c2_out;
0461       cherekov_e = (_eval_run.C2_sum)
0462           > (abs(_eval_run.beam_mom) >= 10 ? 100 : 240);
0463       hNormalization->Fill("C2-e", cherekov_e);
0464
0465       TH2F * hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
0466           "hCheck_Cherenkov"));
0467       assert(hCheck_Cherenkov);
0468       hCheck_Cherenkov->Fill(c1, "C1", 1);
0469       hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
0470       hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
0471       hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
0472     }
0473
0474   // veto
0475   TH1F * hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
0476   assert(hCheck_Veto);
0477   bool trigger_veto_pass = true;
0478   if (not _is_sim)
0479     {
0480       auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
0481       for (auto it = range.first; it != range.second; ++it)
0482         {
0483           RawTower* tower = it->second;
0484           assert(tower);
0485
0486           hCheck_Veto->Fill(tower->get_energy());
0487
0488           if (abs(tower->get_energy()) > 15)
0489             trigger_veto_pass = false;
0490         }
0491     }
0492   hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
0493   _eval_run.trigger_veto_pass = trigger_veto_pass;
0494
0495   // hodoscope
0496   TH1F * hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
0497   assert(hCheck_Hodo_H);
0498   int hodo_h_count = 0;
0499   if (not _is_sim)
0500     {
0501       auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
0502       for (auto it = range.first; it != range.second; ++it)
0503         {
0504           RawTower* tower = it->second;
0505           assert(tower);
0506
0507           hCheck_Hodo_H->Fill(tower->get_energy());
0508
0509           if (abs(tower->get_energy()) > 30)
0510             {
0511               hodo_h_count++;
0512               _eval_run.hodo_h = tower->get_id();
0513             }
0514         }
0515     }
0516   const bool valid_hodo_h = hodo_h_count == 1;
0517   hNormalization->Fill("valid_hodo_h", valid_hodo_h);
0518   _eval_run.valid_hodo_h = valid_hodo_h;
0519
0520   TH1F * hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
0521   assert(hCheck_Hodo_V);
0522   int hodo_v_count = 0;
0523   if (not _is_sim)
0524     {
0525       auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
0526       for (auto it = range.first; it != range.second; ++it)
0527         {
0528           RawTower* tower = it->second;
0529           assert(tower);
0530
0531           hCheck_Hodo_V->Fill(tower->get_energy());
0532
0533           if (abs(tower->get_energy()) > 30)
0534             {
0535               hodo_v_count++;
0536               _eval_run.hodo_v = tower->get_id();
0537             }
0538         }
0539     }
0540   const bool valid_hodo_v = hodo_v_count == 1;
0541   _eval_run.valid_hodo_v = valid_hodo_v;
0542   hNormalization->Fill("valid_hodo_v", valid_hodo_v);
0543
0544   const bool good_e = (valid_hodo_v and valid_hodo_h and cherekov_e
0545       and trigger_veto_pass) and (not _is_sim);
0546   hNormalization->Fill("good_e", good_e);
0547
0548   // simple clustering
0549   pair<int, int> max_3x3 = find_max(TOWER_CALIB_CEMC, 3);
0550   pair<int, int> max_5x5 = find_max(TOWER_CALIB_CEMC, 5);
0551
0552   _eval_3x3_raw.max_col = max_3x3.first;
0553   _eval_3x3_raw.max_row = max_3x3.second;
0554   _eval_3x3_prod.max_col = max_3x3.first;
0555   _eval_3x3_prod.max_row = max_3x3.second;
0556   _eval_3x3_temp.max_col = max_3x3.first;
0557   _eval_3x3_temp.max_row = max_3x3.second;
0558   _eval_3x3_recalib.max_col = max_3x3.first;
0559   _eval_3x3_recalib.max_row = max_3x3.second;
0560
0561   _eval_5x5_raw.max_col = max_5x5.first;
0562   _eval_5x5_raw.max_row = max_5x5.second;
0563   _eval_5x5_prod.max_col = max_5x5.first;
0564   _eval_5x5_prod.max_row = max_5x5.second;
0565   _eval_5x5_temp.max_col = max_5x5.first;
0566   _eval_5x5_temp.max_row = max_5x5.second;
0567   _eval_5x5_recalib.max_col = max_5x5.first;
0568   _eval_5x5_recalib.max_row = max_5x5.second;
0569
0570   // tower
0571   bool good_temp = true;
0572   double sum_energy_calib = 0;
0573   double sum_energy_T = 0;
0574   TH1F * hTemperature = dynamic_cast<TH1F *>(hm->getHisto("hTemperature"));
0575   assert(hTemperature);
0576
0577   stringstream sdata;
0578
0579
0580   //code for position dependence in simulation
0581   std::vector<float> toweretas3x3;
0582   std::vector<float> towerphis3x3;
0583   std::vector<float> towerenergies3x3;
0584   std::vector<float> toweretas5x5;
0585   std::vector<float> towerphis5x5;
0586   std::vector<float> towerenergies5x5;
0587
0588
0589
0590   if (good_e)
0591     sdata << abs(_eval_run.beam_mom) << "\t";
0592
0593   // tower temperature and recalibration
0594     {
0595       auto range = TOWER_CALIB_CEMC->getTowers();
0596       for (auto it = range.first; it != range.second; ++it)
0597         {
0598
0599           RawTowerDefs::keytype key = it->first;
0600           RawTower* tower = it->second;
0601           assert(tower);
0602
0603           const int col = tower->get_bineta();
0604           const int row = tower->get_binphi();
0605       const int hodorow = _eval_run.hodo_h;
0606       const int hodocol = _eval_run.hodo_v;
0607
0608           if (col < 0 or col >= 8)
0609             continue;
0610
0611           if (row < 0 or row >= 8)
0612             continue;
0613
0614       if((hodorow<0 or hodorow>7) and !_is_sim)
0615         continue;
0616       if((hodocol<0 or hodocol>7) and !_is_sim)
0617         continue;
0618
0619           const double energy_calib = tower->get_energy();
0620           sum_energy_calib += energy_calib;
0621
0622           RawTower* tower_raw = TOWER_RAW_CEMC->getTower(key);
0623           assert(tower_raw);
0624
0625           double energy_T = energy_calib;
0626
0627 //          if (not _is_sim)
0628 //            {
0629 //              RawTower_Temperature * temp_t =
0630 //                  dynamic_cast<RawTower_Temperature *>(TOWER_TEMPERATURE_EMCAL->getTower(
0631 //                      tower->get_row(), tower->get_column())); // note swap of col/row in temperature storage
0632 //              assert(temp_t);
0633 //
0634 //              const double T = temp_t->get_temperature_from_time(
0635 //                  eventheader->get_TimeStamp());
0636 //              hTemperature->Fill(T);
0637 //
0638 //              if (T < 25 or T > 35)
0639 //                good_temp = false;
0640 //
0641 //              energy_T = TemperatureCorrection::Apply(energy_calib, T);
0642 //            }
0643
0644           // recalibration
0645           assert(
0646               _recalib_const.find(make_pair(col, row)) != _recalib_const.end());
0647       double energy_recalib=0;
0648
0649       if(_use_hodoscope_calibs)
0650         energy_recalib = energy_T
0651           *hodoscope_recalibrations[hodorow][7-hodocol];
0652
0653
0654
0655
0656
0657           // energy sums
0658           sum_energy_T += energy_T;
0659
0660           // calibration file
0661 //          sdata << tower->get_energy() << "\t";
0662           // calibration file - only output 5x5 towers
0663           if (col >= max_5x5.first - 1 and col <= max_5x5.first + 1
0664               and row >= max_5x5.second - 1 and row <= max_5x5.second + 1)
0665             {
0666               sdata << tower->get_energy() << "\t";
0667             }
0668           else
0669             {
0670               sdata << 0 << "\t";
0671             }
0672
0673
0674
0675           // cluster 3x3
0676           if (col >= max_3x3.first - 1 and col <= max_3x3.first + 1)
0677             if (row >= max_3x3.second - 1 and row <= max_3x3.second + 1)
0678               {
0679                 // in cluster
0680
0681                 _eval_3x3_raw.average_col += abs(tower_raw->get_energy()) * col;
0682                 _eval_3x3_raw.average_row += abs(tower_raw->get_energy()) * row;
0683                 _eval_3x3_raw.sum_E += abs(tower_raw->get_energy());
0684
0685                 _eval_3x3_prod.average_col += energy_calib * col;
0686                 _eval_3x3_prod.average_row += energy_calib * row;
0687                 _eval_3x3_prod.sum_E += energy_calib;
0688
0689                 _eval_3x3_temp.average_col += energy_T * col;
0690                 _eval_3x3_temp.average_row += energy_T * row;
0691                 _eval_3x3_temp.sum_E += energy_T;
0692
0693                 _eval_3x3_recalib.average_col += energy_recalib * col;
0694                 _eval_3x3_recalib.average_row += energy_recalib * row;
0695                 _eval_3x3_recalib.sum_E += energy_recalib;
0696
0697
0698         toweretas3x3.push_back(col);
0699         towerphis3x3.push_back(row);
0700         towerenergies3x3.push_back(energy_calib);
0701
0702
0703
0704               }
0705
0706           // cluster 5x5
0707           if (col >= max_5x5.first - 2 and col <= max_5x5.first + 2)
0708             if (row >= max_5x5.second - 2 and row <= max_5x5.second + 2)
0709               {
0710                 // in cluster
0711
0712                 _eval_5x5_raw.average_col += abs(tower_raw->get_energy()) * col;
0713                 _eval_5x5_raw.average_row += abs(tower_raw->get_energy()) * row;
0714                 _eval_5x5_raw.sum_E += abs(tower_raw->get_energy());
0715
0716                 _eval_5x5_prod.average_col += energy_calib * col;
0717                 _eval_5x5_prod.average_row += energy_calib * row;
0718                 _eval_5x5_prod.sum_E += energy_calib;
0719
0720                 _eval_5x5_temp.average_col += energy_T * col;
0721                 _eval_5x5_temp.average_row += energy_T * row;
0722                 _eval_5x5_temp.sum_E += energy_T;
0723
0724                 _eval_5x5_recalib.average_col += energy_recalib * col;
0725                 _eval_5x5_recalib.average_row += energy_recalib * row;
0726                 _eval_5x5_recalib.sum_E += energy_recalib;
0727
0728         toweretas5x5.push_back(col);
0729         towerphis5x5.push_back(row);
0730         towerenergies5x5.push_back(energy_calib);
0731
0732               }
0733         }
0734     }
0735
0736   _eval_3x3_raw.reweight_clus_pol();
0737   _eval_5x5_raw.reweight_clus_pol();
0738   _eval_3x3_prod.reweight_clus_pol();
0739   _eval_5x5_prod.reweight_clus_pol();
0740   _eval_3x3_temp.reweight_clus_pol();
0741   _eval_5x5_temp.reweight_clus_pol();
0742   _eval_3x3_recalib.reweight_clus_pol();
0743   _eval_5x5_recalib.reweight_clus_pol();
0744
0745
0746
0747
0748   //enter simulation position dependent weighting code
0749
0750
0751   //do 3x3 clusters first
0752   int ntowers = toweretas3x3.size();
0753   //assert(ntowers >= 1);
0754   const int nphibin = towergeom->get_phibins();
0755   if(ntowers>0){
0756
0757     //loop over the towers to determine the energy
0758     //weighted eta and phi position of the cluster
0759
0760     float etamult = 0;
0761     float etasum = 0;
0762     float phimult = 0;
0763     float phisum = 0;
0764
0765     for (int j = 0; j < ntowers; j++)
0766       {
0767     float energymult = towerenergies3x3.at(j) * toweretas3x3.at(j);
0768     etamult += energymult;
0769     etasum += towerenergies3x3.at(j);
0770
0771     int phibin = towerphis3x3.at(j);
0772
0773     if (phibin - towerphis3x3.at(0) < -nphibin / 2)
0774       phibin += nphibin;
0775     else if (phibin - towerphis3x3.at(0) > +nphibin / 2)
0776       phibin -= nphibin;
0777     assert(abs(phibin - towerphis3x3.at(0)) <= nphibin / 2);
0778
0779     energymult = towerenergies3x3.at(j) * phibin;
0780     phimult += energymult;
0781     phisum += towerenergies3x3.at(j);
0782       }
0783
0784
0785     float avgphi = phimult / phisum;
0786     float avgeta = etamult / etasum;
0787
0788     if (avgphi<0) avgphi += nphibin;
0789
0790     //this determines the position of the cluster in the 2x2 block
0791     float fmodphi3x3 = fmod(avgphi, 2.);
0792     float fmodeta3x3 = fmod(avgeta, 2.);
0793
0794     _eval_3x3_prod.setfmods(fmodeta3x3, fmodphi3x3);
0795
0796
0797     int fmodetabin = -99;
0798     int fmodphibin = -99;
0799     for(int k=0; k<nfmodbins-1; k++)
0800       if(fmodeta3x3 >= fmodbins[k] && fmodeta3x3 < fmodbins[k+1])
0801     fmodetabin = k;
0802     for(int k=0; k<nfmodbins-1; k++)
0803       if(fmodphi3x3 >= fmodbins[k] && fmodphi3x3 < fmodbins[k+1])
0804     fmodphibin = k;
0805
0806     if(fmodetabin>-1 && fmodphibin>-1 && !_use_hodoscope_calibs)
0807       _eval_3x3_recalib.sum_E = _eval_3x3_prod.getsumE()/clus_3x3_sim_recalibs[fmodetabin][fmodphibin];
0808     _eval_3x3_recalib.setfmods(fmodeta3x3,fmodphi3x3);
0809
0810
0811   }
0812
0813   //NOW DO THE 5X5 CLUSTERS
0814
0815   ntowers = toweretas5x5.size();
0816   //assert(ntowers >= 1);
0817   if(ntowers>0){
0818     float etamult = 0;
0819     float etasum = 0;
0820     float phimult = 0;
0821     float phisum = 0;
0822     for (int j = 0; j < ntowers; j++)
0823       {
0824     float energymult = towerenergies5x5.at(j) * toweretas5x5.at(j);
0825     etamult += energymult;
0826     etasum += towerenergies5x5.at(j);
0827
0828     int phibin = towerphis5x5.at(j);
0829
0830     if (phibin - towerphis5x5.at(0) < -nphibin / 2)
0831       phibin += nphibin;
0832     else if (phibin - towerphis5x5.at(0) > +nphibin / 2)
0833       phibin -= nphibin;
0834     assert(abs(phibin - towerphis5x5.at(0)) <= nphibin / 2);
0835
0836     energymult = towerenergies5x5.at(j) * phibin;
0837     phimult += energymult;
0838     phisum += towerenergies5x5.at(j);
0839       }
0840
0841     float avgphi = phimult / phisum;
0842     float avgeta = etamult / etasum;
0843
0844     if (avgphi < 0) avgphi += nphibin;
0845
0846     float fmodphi5x5 = fmod(avgphi , 2.);
0847     float fmodeta5x5 = fmod(avgeta , 2.);
0848
0849     _eval_5x5_prod.setfmods(fmodeta5x5, fmodphi5x5);
0850
0851
0852     int fmodetabin = -99;
0853     int fmodphibin = -99;
0854     for(int k=0; k<nfmodbins-1; k++)
0855       if(fmodeta5x5 >= fmodbins[k] && fmodeta5x5 < fmodbins[k+1])
0856     fmodetabin = k;
0857     for(int k=0; k<nfmodbins-1; k++)
0858       if(fmodphi5x5 >= fmodbins[k] && fmodphi5x5 < fmodbins[k+1])
0859     fmodphibin = k;
0860
0861     if(fmodetabin > -1 && fmodphibin > -1 && !_use_hodoscope_calibs)
0862       _eval_5x5_recalib.sum_E = _eval_5x5_prod.getsumE()/clus_5x5_sim_recalibs[fmodetabin][fmodphibin];
0863
0864     _eval_5x5_recalib.setfmods(fmodeta5x5,fmodphi5x5);
0865   }
0866
0867
0868
0869
0870
0871
0872
0873
0874
0875
0876
0877
0878
0879   const double EoP = sum_energy_T / abs(_eval_run.beam_mom);
0880   hNormalization->Fill("good_temp", good_temp);
0881
0882 //  bool good_data = good_e and good_temp;
0883   bool good_data = good_e;
0884   hNormalization->Fill("good_data", good_data);
0885
0886   _eval_run.good_temp = good_temp;
0887   _eval_run.good_e = good_e;
0888   _eval_run.good_data = good_data;
0889   _eval_run.sum_energy_T = sum_energy_T;
0890   _eval_run.EoP = EoP;
0891
0892   // E/p
0893   if (good_data)
0894     {
0895       if (verbosity >= 3)
0896         cout << __PRETTY_FUNCTION__ << " sum_energy_calib = "
0897             << sum_energy_calib << " sum_energy_T = " << sum_energy_T
0898             << " _eval_run.beam_mom = " << _eval_run.beam_mom << endl;
0899
0900       TH2F * hEoP = dynamic_cast<TH2F *>(hm->getHisto("hEoP"));
0901       assert(hEoP);
0902
0903       hEoP->Fill(EoP, abs(_eval_run.beam_mom));
0904     }
0905
0906   // calibration file
0907   if (good_data and abs(_eval_run.beam_mom) >= 4
0908       and abs(_eval_run.beam_mom) <= 8
0909       and abs(_eval_3x3_raw.average_col - round(_eval_3x3_raw.average_col))
0910           < 0.1
0911       and abs(_eval_3x3_raw.average_row - round(_eval_3x3_raw.average_row))
0912           < 0.1)
0913     {
0914       assert(fdata.is_open());
0915
0916       fdata << sdata.str();
0917
0918       fdata << endl;
0919     }
0920
0921   TTree * T = dynamic_cast<TTree *>(hm->getHisto("T"));
0922   assert(T);
0923   T->Fill();
0924
0925   return Fun4AllReturnCodes::EVENT_OK;
0926 }
0927
0928 pair<int, int>
0929 Proto3ShowerCalib::find_max(RawTowerContainer* towers, int cluster_size)
0930 {
0931   const int clus_edge_size = (cluster_size - 1) / 2;
0932   assert(clus_edge_size >= 0);
0933
0934   pair<int, int> max(-1000, -1000);
0935   double max_e = 0;
0936
0937   for (int col = 0; col < n_size; ++col)
0938     for (int row = 0; row < n_size; ++row)
0939       {
0940         double energy = 0;
0941
0942         for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
0943             ++dcol)
0944           for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
0945               ++drow)
0946             {
0947               if (dcol < 0 or drow < 0)
0948                 continue;
0949
0950               RawTower * t = towers->getTower(dcol, drow);
0951               if (t)
0952                 energy += t->get_energy();
0953             }
0954
0955         if (energy > max_e)
0956           {
0957             max = make_pair(col, row);
0958             max_e = energy;
0959           }
0960       }
0961
0962   return max;
0963 }
0964
0965 int
0966 Proto3ShowerCalib::LoadRecalibMap(const std::string & file)
0967 {
0968   if (verbosity)
0969     {
0970       cout << __PRETTY_FUNCTION__ << " - input recalibration constant from "
0971           << file << endl;
0972     }
0973
0974   ifstream fcalib(file);
0975
0976   assert(fcalib.is_open());
0977
0978   string line;
0979
0980   while (not fcalib.eof())
0981     {
0982       getline(fcalib, line);
0983
0984       if (verbosity > 10)
0985         {
0986           cout << __PRETTY_FUNCTION__ << " get line " << line << endl;
0987         }
0988       istringstream sline(line);
0989
0990       int col = -1;
0991       int row = -1;
0992       double calib = 0;
0993
0994       sline >> col >> row >> calib;
0995
0996       if (not sline.fail())
0997         {
0998           if (verbosity)
0999             {
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