EMCTemplateFitting/src/Prototype4DSTReader.cc

0001 // $Id: Prototype4DSTReader.cc,v 1.11 2015/01/06 02:52:07 jinhuang Exp $
0002
0003 /*!
0004  * \file Prototype4DSTReader.cc
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision: 1.11 $
0008  * \date $Date: 2015/01/06 02:52:07 $
0009  */
0010
0011 #include "Prototype4DSTReader.h"
0012
0013 #include <calobase/RawTower.h>  // for RawTower
0014 #include <calobase/RawTowerContainer.h>
0015
0016 #include <pdbcalbase/PdbParameterMap.h>
0017
0018 #include <phparameter/PHParameters.h>
0019
0020 #include <fun4all/PHTFileServer.h>
0021 #include <fun4all/SubsysReco.h>  // for SubsysReco
0022
0023 #include <phool/getClass.h>
0024
0025 #include <TClass.h>        // for TClass
0026 #include <TClonesArray.h>  // for TClonesArray
0027 #include <TObject.h>       // for TObject
0028 #include <TString.h>       // for Form
0029 #include <TTree.h>
0030
0031 #include <cassert>
0032 #include <iostream>  // for operator<<, basic_ostream, endl
0033 #include <limits>
0034 #include <map>
0035 #include <utility>                       // for pair
0036
0037 class PHCompositeNode;
0038
0039 using namespace std;
0040
0041 Prototype4DSTReader::Prototype4DSTReader(const string &filename)
0042   : SubsysReco("Prototype4DSTReader")
0043   , nblocks(0)
0044   , _event(0)
0045   ,  //
0046   _out_file_name(filename)
0047   , /*_file(nullptr), */ _T(nullptr)
0048   ,  //
0049   _tower_zero_sup(-10000000)
0050 {
0051 }
0052
0053 Prototype4DSTReader::~Prototype4DSTReader()
0054 {
0055   cout << "Prototype4DSTReader::destructor - Clean ups" << endl;
0056
0057   if (_T)
0058   {
0059     _T->ResetBranchAddresses();
0060   }
0061
0062   _records.clear();
0063 }
0064
0065 int Prototype4DSTReader::Init(PHCompositeNode *)
0066 {
0067   const static int arr_size = 100;
0068
0069   if (_tower_postfix.size())
0070   {
0071     cout << "Prototype4DSTReader::Init - zero suppression for calorimeter "
0072             "towers = "
0073          << _tower_zero_sup << " GeV" << endl;
0074   }
0075   for (vector<string>::const_iterator it = _runinfo_list.begin();
0076        it != _runinfo_list.end(); ++it)
0077   {
0078     const string &nodenam = *it;
0079
0080     record rec;
0081     rec._cnt = 0;
0082     rec._name = nodenam;
0083     rec._arr = nullptr;
0084     rec._arr_ptr = nullptr;
0085     rec._dvalue = 0;
0086     rec._type = record::typ_runinfo;
0087
0088     _records.push_back(rec);
0089
0090     nblocks++;
0091   }
0092   for (vector<string>::const_iterator it = _eventinfo_list.begin();
0093        it != _eventinfo_list.end(); ++it)
0094   {
0095     const string &nodenam = *it;
0096
0097     record rec;
0098     rec._cnt = 0;
0099     rec._name = nodenam;
0100     rec._arr = nullptr;
0101     rec._arr_ptr = nullptr;
0102     rec._dvalue = 0;
0103     rec._type = record::typ_eventinfo;
0104
0105     _records.push_back(rec);
0106
0107     nblocks++;
0108   }
0109
0110   for (vector<string>::const_iterator it = _tower_postfix.begin();
0111        it != _tower_postfix.end(); ++it)
0112   {
0113     const char *class_name = RawTower_type::Class()->GetName();
0114
0115     const string &nodenam = *it;
0116
0117     string hname = Form("TOWER_%s", nodenam.c_str());
0118     //      _node_name.push_back(hname);
0119     cout << "Prototype4DSTReader::Init - saving raw tower info from node: "
0120          << hname << " - " << class_name << endl;
0121
0122     record rec;
0123     rec._cnt = 0;
0124     rec._name = hname;
0125     rec._arr = make_shared<TClonesArray>(class_name, arr_size);
0126     rec._arr_ptr = rec._arr.get();
0127     rec._dvalue = 0;
0128     rec._type = record::typ_tower;
0129
0130     _records.push_back(rec);
0131
0132     nblocks++;
0133   }
0134
0135   for (vector<string>::const_iterator it = _towertemp_postfix.begin();
0136        it != _towertemp_postfix.end(); ++it)
0137   {
0138     const string &nodenam = *it;
0139     string hname = Form("TOWER_TEMPERATURE_%s", nodenam.c_str());
0140
0141     cout << "Prototype4DSTReader::Init - saving average tower temperature info "
0142             "from node: "
0143          << hname << endl;
0144
0145     record rec;
0146     rec._cnt = 0;
0147     rec._name = hname;
0148     rec._arr = nullptr;
0149     rec._arr_ptr = nullptr;
0150     rec._dvalue = 0;
0151     rec._type = record::typ_towertemp;
0152
0153     _records.push_back(rec);
0154
0155     nblocks++;
0156   }
0157   cout << "Prototype4DSTReader::Init - requested " << nblocks << " nodes"
0158        << endl;
0159
0160   build_tree();
0161
0162   return 0;
0163 }
0164
0165 void Prototype4DSTReader::build_tree()
0166 {
0167   cout << "Prototype4DSTReader::build_tree - output to " << _out_file_name
0168        << endl;
0169
0170   static const int BUFFER_SIZE = 32000;
0171
0172   // open TFile
0173   PHTFileServer::get().open(_out_file_name, "RECREATE");
0174
0175   _T = new TTree("T", "Prototype4DSTReader");
0176
0177   nblocks = 0;
0178   for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0179   {
0180     record &rec = *it;
0181
0182     cout << "Prototype4DSTReader::build_tree - Add " << rec._name << endl;
0183
0184     if (rec._type == record::typ_runinfo)
0185     {
0186       const string name_cnt = rec._name;
0187       const string name_cnt_desc = name_cnt + "/D";
0188       _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
0189                  BUFFER_SIZE);
0190     }
0191     if (rec._type == record::typ_eventinfo)
0192     {
0193       const string name_cnt = rec._name;
0194       const string name_cnt_desc = name_cnt + "/D";
0195       _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
0196                  BUFFER_SIZE);
0197     }
0198     else if (rec._type == record::typ_tower)
0199     {
0200       const string name_cnt = "n_" + rec._name;
0201       const string name_cnt_desc = name_cnt + "/I";
0202       _T->Branch(name_cnt.c_str(), &(rec._cnt), name_cnt_desc.c_str(),
0203                  BUFFER_SIZE);
0204       _T->Branch(rec._name.c_str(), &(rec._arr_ptr), BUFFER_SIZE, 99);
0205     }
0206     else if (rec._type == record::typ_towertemp)
0207     {
0208       const string name_cnt = rec._name + "_AVG";
0209       const string name_cnt_desc = name_cnt + "/D";
0210       _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
0211                  BUFFER_SIZE);
0212     }
0213
0214     nblocks++;
0215   }
0216
0217   cout << "Prototype4DSTReader::build_tree - added " << nblocks << " nodes"
0218        << endl;
0219
0220   _T->SetAutoSave(16000);
0221 }
0222
0223 int Prototype4DSTReader::process_event(PHCompositeNode *topNode)
0224 {
0225   //  const double significand = _event / TMath::Power(10, (int)
0226   //  (log10(_event)));
0227   //
0228   //  if (fmod(significand, 1.0) == 0 && significand <= 10)
0229   //    cout << "Prototype4DSTReader::process_event - " << _event << endl;
0230   _event++;
0231
0232   for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0233   {
0234     record &rec = *it;
0235
0236     rec._cnt = 0;
0237
0238     if (rec._type == record::typ_hit)
0239     {
0240       assert(0);
0241     }  //      if (rec._type == record::typ_hit)
0242     else if (rec._type == record::typ_tower)
0243     {
0244       assert(rec._arr.get() == rec._arr_ptr);
0245       assert(rec._arr.get());
0246       rec._arr->Clear();
0247
0248       if (Verbosity() >= 2)
0249         cout << "Prototype4DSTReader::process_event - processing tower "
0250              << rec._name << endl;
0251
0252       RawTowerContainer *hits =
0253           findNode::getClass<RawTowerContainer>(topNode, rec._name);
0254       if (!hits)
0255       {
0256         if (_event < 2)
0257           cout << "Prototype4DSTReader::process_event - Error - can not find "
0258                   "node "
0259                << rec._name << endl;
0260       }
0261       else
0262       {
0263         RawTowerContainer::ConstRange hit_range = hits->getTowers();
0264
0265         if (Verbosity() >= 2)
0266           cout << "Prototype4DSTReader::process_event - processing "
0267                << rec._name << " and received " << hits->size() << " tower hits"
0268                << endl;
0269
0270         for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
0271              hit_iter != hit_range.second; hit_iter++)
0272         {
0273           RawTower *hit_raw = hit_iter->second;
0274
0275           RawTower_type *hit = dynamic_cast<RawTower_type *>(hit_raw);
0276           //                  RawTower * hit = hit_iter->second;
0277
0278           assert(hit);
0279
0280           if (hit->get_energy() < _tower_zero_sup)
0281           {
0282             if (Verbosity() >= 2)
0283               cout << "Prototype4DSTReader::process_event - suppress low "
0284                       "energy tower hit "
0285                    << rec._name
0286                    << " @ ("
0287                    //                            << hit->get_thetaMin()
0288                    //                            << ", " << hit->get_phiMin()
0289                    << "), Energy = " << hit->get_energy() << endl;
0290
0291             continue;
0292           }
0293
0294           new ((*(rec._arr.get()))[rec._cnt]) RawTower_type();
0295
0296           if (Verbosity() >= 2)
0297             cout << "Prototype4DSTReader::process_event - processing Tower hit "
0298                  << rec._name
0299                  << " @ ("
0300                  //                        << hit->get_thetaMin() << ", "
0301                  //                        << hit->get_phiMin()
0302                  << "), Energy = " << hit->get_energy() << " - "
0303                  << rec._arr.get()->At(rec._cnt)->ClassName() << endl;
0304
0305           //                  rec._arr->Print();
0306
0307           RawTower_type *new_hit =
0308               dynamic_cast<RawTower_type *>(rec._arr.get()->At(rec._cnt));
0309           assert(new_hit);
0310
0311           *new_hit = (*hit);
0312
0313           rec._cnt++;
0314         }
0315       }  // if (!hits)
0316     }    //      if (rec._type == record::typ_hit)
0317     else if (rec._type == record::typ_towertemp)
0318     {
0319       if (Verbosity() >= 2)
0320         cout << "Prototype4DSTReader::process_event - processing tower "
0321                 "temperature "
0322              << rec._name << endl;
0323
0324       RawTowerContainer *hits =
0325           findNode::getClass<RawTowerContainer>(topNode, rec._name);
0326       if (!hits)
0327       {
0328         if (_event < 2)
0329           cout << "Prototype4DSTReader::process_event - Error - can not find "
0330                   "node "
0331                << rec._name << endl;
0332       }
0333       else
0334       {
0335         RawTowerContainer::ConstRange hit_range = hits->getTowers();
0336
0337         if (Verbosity() >= 2)
0338           cout << "Prototype4DSTReader::process_event - processing "
0339                << rec._name << " and received " << hits->size() << " tower hits"
0340                << endl;
0341
0342         rec._cnt = 0;
0343
0344         for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
0345              hit_iter != hit_range.second; hit_iter++)
0346         {
0347           RawTower *hit_raw = hit_iter->second;
0348
0349           RawTowerT_type *hit = dynamic_cast<RawTowerT_type *>(hit_raw);
0350           //                  RawTower * hit = hit_iter->second;
0351
0352           assert(hit);
0353
0354           rec._dvalue +=
0355               hit->get_temperature_from_entry(hit->get_nr_entries() - 1);
0356
0357           ++rec._cnt;
0358         }
0359
0360         rec._dvalue /= rec._cnt;
0361       }  // if (!hits)
0362     }    //      if (rec._type == record::typ_hit)
0363     else if (rec._type == record::typ_runinfo)
0364     {
0365       PdbParameterMap *info =
0366           findNode::getClass<PdbParameterMap>(topNode, "RUN_INFO");
0367
0368       assert(info);
0369
0370       PHParameters run_info_copy("RunInfo");
0371       run_info_copy.FillFrom(info);
0372
0373       rec._dvalue = run_info_copy.get_double_param(rec._name);
0374
0375     }  //
0376     else if (rec._type == record::typ_eventinfo)
0377     {
0378       PdbParameterMap *info =
0379           findNode::getClass<PdbParameterMap>(topNode, "EVENT_INFO");
0380
0381       assert(info);
0382
0383       PHParameters event_info_copy("EVENT_INFO");
0384       event_info_copy.FillFrom(info);
0385
0386       rec._dvalue = numeric_limits<double>::quiet_NaN();
0387
0388       if (event_info_copy.exist_double_param(rec._name))
0389         rec._dvalue = event_info_copy.get_double_param(rec._name);
0390
0391     }  //      if (rec._type == record::typ_hit)
0392     else if (rec._type == record::typ_part)
0393     {
0394       assert(0);
0395     }  //      if (rec._type == record::typ_part)
0396     else if (rec._type == record::typ_vertex)
0397     {
0398       assert(0);
0399     }  //          if (_load_all_particle)
0400
0401   }  //  for (records_t::iterator it = _records.begin(); it != _records.end();
0402   //  ++it)
0403
0404   if (_T)
0405     _T->Fill();
0406
0407   return 0;
0408 }  //  for (records_t::iterator it = _records.begin(); it != _records.end();
0409 //  ++it)
0410
0411 int Prototype4DSTReader::End(PHCompositeNode * /*topNode*/)
0412 {
0413   cout << "Prototype4DSTReader::End - Clean ups" << endl;
0414
0415   if (_T)
0416   {
0417     PHTFileServer::get().cd(_out_file_name);
0418     _T->Write();
0419     _T->ResetBranchAddresses();
0420   }
0421
0422   _records.clear();
0423
0424   return 0;
0425 }