packages/Prototype3/Prototype3DSTReader.cc

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