packages/Prototype2/Prototype2DSTReader.cc

0001 // $Id: Prototype2DSTReader.cc,v 1.11 2015/01/06 02:52:07 jinhuang Exp $
0002
0003 /*!
0004  * \file Prototype2DSTReader.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 "Prototype2DSTReader.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 <map>
0034 #include <sstream>
0035 #include <utility>                       // for pair
0036
0037 class PHCompositeNode;
0038
0039 using namespace std;
0040
0041 Prototype2DSTReader::Prototype2DSTReader(const string &filename)
0042   : SubsysReco("Prototype2DSTReader")
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 Prototype2DSTReader::~Prototype2DSTReader()
0054 {
0055   cout << "Prototype2DSTReader::destructor - Clean ups" << endl;
0056
0057   if (_T)
0058   {
0059     _T->ResetBranchAddresses();
0060   }
0061
0062   _records.clear();
0063 }
0064
0065 int Prototype2DSTReader::Init(PHCompositeNode *)
0066 {
0067   const static int arr_size = 100;
0068
0069   if (_tower_postfix.size())
0070   {
0071     cout << "Prototype2DSTReader::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
0093   for (vector<string>::const_iterator it = _tower_postfix.begin();
0094        it != _tower_postfix.end(); ++it)
0095   {
0096     const char *class_name = RawTower_type::Class()->GetName();
0097
0098     const string &nodenam = *it;
0099
0100     string hname = Form("TOWER_%s", nodenam.c_str());
0101     //      _node_name.push_back(hname);
0102     cout << "Prototype2DSTReader::Init - saving raw tower info from node: "
0103          << hname << " - " << class_name << endl;
0104
0105     record rec;
0106     rec._cnt = 0;
0107     rec._name = hname;
0108     rec._arr = make_shared<TClonesArray>(class_name, arr_size);
0109     rec._arr_ptr = rec._arr.get();
0110     rec._dvalue = 0;
0111     rec._type = record::typ_tower;
0112
0113     _records.push_back(rec);
0114
0115     nblocks++;
0116   }
0117
0118   for (vector<string>::const_iterator it = _towertemp_postfix.begin();
0119        it != _towertemp_postfix.end(); ++it)
0120   {
0121     const string &nodenam = *it;
0122     string hname = Form("TOWER_TEMPERATURE_%s", nodenam.c_str());
0123
0124     cout << "Prototype2DSTReader::Init - saving average tower temperature info "
0125             "from node: "
0126          << hname << endl;
0127
0128     record rec;
0129     rec._cnt = 0;
0130     rec._name = hname;
0131     rec._arr = nullptr;
0132     rec._arr_ptr = nullptr;
0133     rec._dvalue = 0;
0134     rec._type = record::typ_towertemp;
0135
0136     _records.push_back(rec);
0137
0138     nblocks++;
0139   }
0140   cout << "Prototype2DSTReader::Init - requested " << nblocks << " nodes"
0141        << endl;
0142
0143   build_tree();
0144
0145   return 0;
0146 }
0147
0148 void Prototype2DSTReader::build_tree()
0149 {
0150   cout << "Prototype2DSTReader::build_tree - output to " << _out_file_name
0151        << endl;
0152
0153   static const int BUFFER_SIZE = 32000;
0154
0155   // open TFile
0156   PHTFileServer::get().open(_out_file_name, "RECREATE");
0157
0158   _T = new TTree("T", "Prototype2DSTReader");
0159
0160   nblocks = 0;
0161   for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0162   {
0163     record &rec = *it;
0164
0165     cout << "Prototype2DSTReader::build_tree - Add " << rec._name << endl;
0166
0167     if (rec._type == record::typ_runinfo)
0168     {
0169       const string name_cnt = rec._name;
0170       const string name_cnt_desc = name_cnt + "/D";
0171       _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
0172                  BUFFER_SIZE);
0173     }
0174     else if (rec._type == record::typ_tower)
0175     {
0176       const string name_cnt = "n_" + rec._name;
0177       const string name_cnt_desc = name_cnt + "/I";
0178       _T->Branch(name_cnt.c_str(), &(rec._cnt), name_cnt_desc.c_str(),
0179                  BUFFER_SIZE);
0180       _T->Branch(rec._name.c_str(), &(rec._arr_ptr), BUFFER_SIZE, 99);
0181     }
0182     else if (rec._type == record::typ_towertemp)
0183     {
0184       const string name_cnt = rec._name + "_AVG";
0185       const string name_cnt_desc = name_cnt + "/D";
0186       _T->Branch(name_cnt.c_str(), &(rec._dvalue), name_cnt_desc.c_str(),
0187                  BUFFER_SIZE);
0188     }
0189
0190     nblocks++;
0191   }
0192
0193   cout << "Prototype2DSTReader::build_tree - added " << nblocks << " nodes"
0194        << endl;
0195
0196   _T->SetAutoSave(16000);
0197 }
0198
0199 int Prototype2DSTReader::process_event(PHCompositeNode *topNode)
0200 {
0201   //  const double significand = _event / TMath::Power(10, (int)
0202   //  (log10(_event)));
0203   //
0204   //  if (fmod(significand, 1.0) == 0 && significand <= 10)
0205   //    cout << "Prototype2DSTReader::process_event - " << _event << endl;
0206   _event++;
0207
0208   for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0209   {
0210     record &rec = *it;
0211
0212     rec._cnt = 0;
0213
0214     if (rec._type == record::typ_hit)
0215     {
0216       assert(0);
0217     }  //      if (rec._type == record::typ_hit)
0218     else if (rec._type == record::typ_tower)
0219     {
0220       assert(rec._arr.get() == rec._arr_ptr);
0221       assert(rec._arr.get());
0222       rec._arr->Clear();
0223
0224       if (Verbosity() >= 2)
0225         cout << "Prototype2DSTReader::process_event - processing tower "
0226              << rec._name << endl;
0227
0228       RawTowerContainer *hits =
0229           findNode::getClass<RawTowerContainer>(topNode, rec._name);
0230       if (!hits)
0231       {
0232         if (_event < 2)
0233           cout << "Prototype2DSTReader::process_event - Error - can not find "
0234                   "node "
0235                << rec._name << endl;
0236       }
0237       else
0238       {
0239         RawTowerContainer::ConstRange hit_range = hits->getTowers();
0240
0241         if (Verbosity() >= 2)
0242           cout << "Prototype2DSTReader::process_event - processing "
0243                << rec._name << " and received " << hits->size() << " tower hits"
0244                << endl;
0245
0246         for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
0247              hit_iter != hit_range.second; hit_iter++)
0248         {
0249           RawTower *hit_raw = hit_iter->second;
0250
0251           RawTower_type *hit = dynamic_cast<RawTower_type *>(hit_raw);
0252           //                  RawTower * hit = hit_iter->second;
0253
0254           assert(hit);
0255
0256           if (hit->get_energy() < _tower_zero_sup)
0257           {
0258             if (Verbosity() >= 2)
0259               cout << "Prototype2DSTReader::process_event - suppress low "
0260                       "energy tower hit "
0261                    << rec._name
0262                    << " @ ("
0263                    //                            << hit->get_thetaMin()
0264                    //                            << ", " << hit->get_phiMin()
0265                    << "), Energy = " << hit->get_energy() << endl;
0266
0267             continue;
0268           }
0269
0270           new ((*(rec._arr.get()))[rec._cnt]) RawTower_type();
0271
0272           if (Verbosity() >= 2)
0273             cout << "Prototype2DSTReader::process_event - processing Tower hit "
0274                  << rec._name
0275                  << " @ ("
0276                  //                        << hit->get_thetaMin() << ", "
0277                  //                        << hit->get_phiMin()
0278                  << "), Energy = " << hit->get_energy() << " - "
0279                  << rec._arr.get()->At(rec._cnt)->ClassName() << endl;
0280
0281           //                  rec._arr->Print();
0282
0283           RawTower_type *new_hit =
0284               dynamic_cast<RawTower_type *>(rec._arr.get()->At(rec._cnt));
0285           assert(new_hit);
0286
0287           *new_hit = (*hit);
0288
0289           rec._cnt++;
0290         }
0291       }  // if (!hits)
0292     }    //      if (rec._type == record::typ_hit)
0293     else if (rec._type == record::typ_towertemp)
0294     {
0295       if (Verbosity() >= 2)
0296         cout << "Prototype2DSTReader::process_event - processing tower "
0297                 "temperature "
0298              << rec._name << endl;
0299
0300       RawTowerContainer *hits =
0301           findNode::getClass<RawTowerContainer>(topNode, rec._name);
0302       if (!hits)
0303       {
0304         if (_event < 2)
0305           cout << "Prototype2DSTReader::process_event - Error - can not find "
0306                   "node "
0307                << rec._name << endl;
0308       }
0309       else
0310       {
0311         RawTowerContainer::ConstRange hit_range = hits->getTowers();
0312
0313         if (Verbosity() >= 2)
0314           cout << "Prototype2DSTReader::process_event - processing "
0315                << rec._name << " and received " << hits->size() << " tower hits"
0316                << endl;
0317
0318         rec._cnt = 0;
0319
0320         for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
0321              hit_iter != hit_range.second; hit_iter++)
0322         {
0323           RawTower *hit_raw = hit_iter->second;
0324
0325           RawTowerT_type *hit = dynamic_cast<RawTowerT_type *>(hit_raw);
0326           //                  RawTower * hit = hit_iter->second;
0327
0328           assert(hit);
0329
0330           rec._dvalue +=
0331               hit->get_temperature_from_entry(hit->get_nr_entries() - 1);
0332
0333           ++rec._cnt;
0334         }
0335
0336         rec._dvalue /= rec._cnt;
0337       }  // if (!hits)
0338     }    //      if (rec._type == record::typ_hit)
0339     else if (rec._type == record::typ_runinfo)
0340     {
0341       PdbParameterMap *info =
0342           findNode::getClass<PdbParameterMap>(topNode, "RUN_INFO");
0343
0344       assert(info);
0345
0346       PHParameters run_info_copy("RunInfo");
0347       run_info_copy.FillFrom(info);
0348
0349       rec._dvalue = run_info_copy.get_double_param(rec._name);
0350
0351     }  //      if (rec._type == record::typ_hit)
0352     else if (rec._type == record::typ_part)
0353     {
0354       assert(0);
0355     }  //      if (rec._type == record::typ_part)
0356     else if (rec._type == record::typ_vertex)
0357     {
0358       assert(0);
0359     }  //          if (_load_all_particle)
0360
0361   }  //  for (records_t::iterator it = _records.begin(); it != _records.end();
0362   //  ++it)
0363
0364   if (_T)
0365     _T->Fill();
0366
0367   return 0;
0368 }  //  for (records_t::iterator it = _records.begin(); it != _records.end();
0369 //  ++it)
0370
0371 int Prototype2DSTReader::End(PHCompositeNode * /*topNode*/)
0372 {
0373   cout << "Prototype2DSTReader::End - Clean ups" << endl;
0374
0375   if (_T)
0376   {
0377     PHTFileServer::get().cd(_out_file_name);
0378     _T->Write();
0379     _T->ResetBranchAddresses();
0380   }
0381
0382   _records.clear();
0383
0384   return 0;
0385 }