Event-Display/src/caloTreeGen.cc

0001 #include "caloTreeGen.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <phool/PHCompositeNode.h>
0006
0007 //Fun4All
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009 #include <fun4all/Fun4AllServer.h>
0010 #include <fun4all/Fun4AllHistoManager.h>
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>
0014 #include <ffaobjects/EventHeader.h>
0015
0016 //ROOT stuff
0017 #include <TH1F.h>
0018 #include <TH2F.h>
0019 #include <TH3F.h>
0020 #include <TFile.h>
0021 #include <TLorentzVector.h>
0022 #include <TTree.h>
0023
0024 //For clusters and geometry
0025 #include <calobase/RawCluster.h>
0026 #include <calobase/RawClusterContainer.h>
0027 #include <calobase/RawClusterUtility.h>
0028 #include <calobase/RawTowerGeomContainer.h>
0029
0030 //Tower stuff
0031 #include <calobase/TowerInfoContainer.h>
0032 #include <calobase/TowerInfo.h>
0033 #include <calobase/TowerInfoDefs.h>
0034
0035 //thingy
0036 #include <CLHEP/Geometry/Point3D.h>
0037
0038
0039 //____________________________________________________________________________..
0040 caloTreeGen::caloTreeGen(const std::string &name,
0041                          const std::string &eventFile,
0042                          const int eventNumber,
0043                          const int runid,
0044                          const float tow_cemc_min,
0045                          const float tow_hcalin_min,
0046                          const float tow_hcalout_min):
0047   SubsysReco(name)
0048   ,T(nullptr)
0049   ,Outfile(name)
0050   ,eventFile(eventFile)
0051   ,totalCaloE(0)
0052   ,iEvent(0)
0053   ,eventNumber(eventNumber)
0054   ,runid(runid)
0055   ,tow_cemc_min(tow_cemc_min)
0056   ,tow_hcalin_min(tow_hcalin_min)
0057   ,tow_hcalout_min(tow_hcalout_min) {
0058
0059   std::cout << "caloTreeGen::caloTreeGen(const std::string &name) Calling ctor" << std::endl;
0060 }
0061
0062 //____________________________________________________________________________..
0063 caloTreeGen::~caloTreeGen() {
0064   std::cout << "caloTreeGen::~caloTreeGen() Calling dtor" << std::endl;
0065 }
0066
0067 //____________________________________________________________________________..
0068 int caloTreeGen::Init(PHCompositeNode *topNode) {
0069
0070   eventOutput.open(eventFile.c_str());
0071
0072   out = new TFile(Outfile.c_str(),"RECREATE");
0073
0074   T = new TTree("T","T");
0075
0076   //emc
0077   T -> Branch("emcTowE",&m_emcTowE);
0078   T -> Branch("emcTowEta",&m_emcTowEta);
0079   T -> Branch("emcTowPhi",&m_emcTowPhi);
0080   T -> Branch("emcTiming",&m_emcTiming);
0081   T -> Branch("emciEta",&m_emciEta);
0082   T -> Branch("emciPhi",&m_emciPhi);
0083
0084   T -> Branch("ohcTowE",&m_ohcTowE);
0085   T -> Branch("ohcTowEta",&m_ohciTowEta);
0086   T -> Branch("ohcTowPhi",&m_ohciTowPhi);
0087
0088   T -> Branch("ihcTowE",&m_ihcTowE);
0089   T -> Branch("ihcTowEta",&m_ihciTowEta);
0090   T -> Branch("ihcTowPhi",&m_ihciTowPhi);
0091
0092   T -> Branch("totalCaloE",&totalCaloE);
0093
0094   //so that the histos actually get written out
0095   Fun4AllServer *se = Fun4AllServer::instance();
0096   se -> Print("NODETREE");
0097
0098   eventOutput << "{\n \
0099       \"EVENT\": {\n \
0100       \"runid\": " << runid << ",\n \
0101       \"evtid\": " << eventNumber << ",\n \
0102        \"time\": 0,\n \
0103        \"type\": \"Collision\",\n \
0104        \"s_nn\": 200,\n \
0105           \"B\": 3.0,\n \
0106          \"pv\": [0,0,0]\n \
0107       },\n \
0108       \"META\": {\n \
0109         \"HITS\": {\n \
0110           \"INNERTRACKER\": {\n \
0111             \"type\": \"3D\",\n \
0112             \"options\": {\n \
0113               \"size\": 5,\n \
0114               \"color\": 16777215\n \
0115             }\n \
0116           },\n \
0117           \"CEMC\": {\n \
0118             \"type\": \"PROJECTIVE\",\n \
0119               \"options\": {\n \
0120                 \"rmin\": 90,\n \
0121                 \"rmax\": 116.1,\n \
0122                 \"deta\": 0.025,\n \
0123                 \"dphi\": 0.025,\n \
0124                 \"color\": 16766464,\n \
0125                 \"transparent\": 0.6,\n \
0126                 \"scaleminmax\": false\n \
0127               }\n \
0128           },\n \
0129           \"HCALIN\": {\n \
0130             \"type\": \"PROJECTIVE\",\n \
0131               \"options\": {\n \
0132                 \"rmin\": 117.27,\n \
0133                 \"rmax\": 165.0,\n \
0134                 \"deta\": 0.025,\n \
0135                 \"dphi\": 0.025,\n \
0136                 \"color\": 4290445312,\n \
0137                 \"transparent\": 0.6,\n \
0138                 \"scaleminmax\": false\n \
0139               }\n \
0140           },\n \
0141           \"HCALOUT\": {\n \
0142             \"type\": \"PROJECTIVE\",\n \
0143               \"options\": {\n \
0144                 \"rmin\": 183.3,\n \
0145                 \"rmax\": 274.317,\n \
0146                 \"deta\": 0.025,\n \
0147                 \"dphi\": 0.025,\n \
0148                 \"color\": 24773,\n \
0149                 \"transparent\": 0.6,\n \
0150                 \"scaleminmax\": false\n \
0151               }\n \
0152           }\n \
0153         }\n \
0154       },\n \
0155   \"HITS\": {\n";
0156
0157   std::cout << "caloTreeGen::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0158   return Fun4AllReturnCodes::EVENT_OK;
0159 }
0160
0161 //____________________________________________________________________________..
0162 int caloTreeGen::InitRun(PHCompositeNode *topNode) {
0163   std::cout << "caloTreeGen::InitRun(PHCompositeNode *topNode) Initializing for Run " << runid << std::endl;
0164   return Fun4AllReturnCodes::EVENT_OK;
0165 }
0166
0167 //____________________________________________________________________________..
0168 int caloTreeGen::process_event(PHCompositeNode *topNode) {
0169
0170   // skip events until event number is reached
0171   if(iEvent != eventNumber) {
0172     ++iEvent;
0173     return 0;
0174   }
0175   else std::cout << "Processing Event: " << iEvent << std::endl;
0176
0177   //tower information
0178   TowerInfoContainer *emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_CEMC");
0179   // TowerInfoContainer *emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_CEMC_RETOWER");
0180
0181   if(!emcTowerContainer) {
0182     std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERS_CEMC"  << std::endl;
0183     return Fun4AllReturnCodes::ABORTEVENT;
0184   }
0185
0186
0187   //Tower geometry node for location information
0188   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0189   // RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0190
0191   if (!towergeom) {
0192     std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERGEOM_CEMC"  << std::endl;
0193     return Fun4AllReturnCodes::ABORTEVENT;
0194   }
0195
0196   float calib = 1.;
0197
0198   //grab all the towers and fill their energies.
0199   unsigned int tower_range = emcTowerContainer->size();
0200   totalCaloE = 0;
0201
0202   bool first_entry = true;
0203
0204   for(unsigned int iter = 0; iter < tower_range; iter++) {
0205     unsigned int towerkey = emcTowerContainer->encode_key(iter);
0206     unsigned int ieta = getCaloTowerEtaBin(towerkey);
0207     unsigned int iphi = getCaloTowerPhiBin(towerkey);
0208     m_emciEta.push_back(ieta);
0209     m_emciPhi.push_back(iphi);
0210
0211     double energy = emcTowerContainer -> get_tower_at_channel(iter) -> get_energy()/calib;
0212     totalCaloE += energy;
0213     double time = emcTowerContainer -> get_tower_at_channel(iter) -> get_time();
0214
0215     double eta = towergeom -> get_etacenter(ieta);
0216     double phi = towergeom -> get_phicenter(iphi);
0217     m_emcTowEta.push_back(eta);
0218     m_emcTowPhi.push_back(phi);
0219     m_emcTowE.push_back(energy);
0220     m_emcTiming.push_back(time);
0221
0222     if(energy < tow_cemc_min) continue;
0223
0224     if(first_entry) {
0225         eventOutput << "\"CEMC\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << energy << "}\n";
0226         first_entry = false;
0227     }
0228     else eventOutput << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << energy << "}\n";
0229   }
0230   if(tower_range == 0) eventOutput << "\"CEMC\": [{}],\n";
0231   else          eventOutput << "],\n";
0232
0233   first_entry = true;
0234
0235   //tower information
0236   TowerInfoContainer *ohcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_HCALOUT");
0237   if(!ohcTowerContainer) {
0238     std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERINFO_CALIB_OHCAL"  << std::endl;
0239     return Fun4AllReturnCodes::ABORTEVENT;
0240   }
0241
0242   //Tower geometry node for location information
0243   towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0244
0245   if (!towergeom) {
0246     std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERGEOM_HCALOUT"  << std::endl;
0247     return Fun4AllReturnCodes::ABORTEVENT;
0248   }
0249
0250   tower_range = ohcTowerContainer->size();
0251
0252   for(unsigned int iter = 0; iter < tower_range; iter++) {
0253     unsigned int towerkey = ohcTowerContainer->encode_key(iter);
0254     unsigned int ieta = getCaloTowerEtaBin(towerkey);
0255     unsigned int iphi = getCaloTowerPhiBin(towerkey);
0256     m_ohciTowEta.push_back(ieta);
0257     m_ohciTowPhi.push_back(iphi);
0258
0259     // double eta = eta_map[ieta];
0260     // double phi = phi_map[iphi];
0261     double eta = towergeom -> get_etacenter(ieta);
0262     double phi = towergeom -> get_phicenter(iphi);
0263
0264     double energy = ohcTowerContainer -> get_tower_at_channel(iter) -> get_energy();
0265     m_ohcTowE.push_back(energy);
0266
0267     if(energy < tow_hcalout_min) continue;
0268
0269     if(first_entry) {
0270       eventOutput << "\"HCALOUT\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << energy << "}\n";
0271       first_entry = false;
0272     }
0273     else eventOutput << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << energy << "}\n";
0274   }
0275
0276   if(tower_range == 0) eventOutput << "\"HCALOUT\": [{}],\n";
0277   else          eventOutput << "],\n";
0278
0279   first_entry = true;
0280
0281   TowerInfoContainer *ihcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_HCALIN");
0282   if(!ihcTowerContainer) {
0283     std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERINFO_CALIB_HCALIN"  << std::endl;
0284     return Fun4AllReturnCodes::ABORTEVENT;
0285   }
0286
0287   //Tower geometry node for location information
0288   towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0289
0290   if (!towergeom) {
0291     std::cout << PHWHERE << "caloTreeGen::process_event Could not find node TOWERGEOM_HCALIN"  << std::endl;
0292     return Fun4AllReturnCodes::ABORTEVENT;
0293   }
0294
0295   tower_range = ihcTowerContainer->size();
0296
0297   for(unsigned int iter = 0; iter < tower_range; iter++) {
0298     unsigned int towerkey = ihcTowerContainer->encode_key(iter);
0299     unsigned int ieta = getCaloTowerEtaBin(towerkey);
0300     unsigned int iphi = getCaloTowerPhiBin(towerkey);
0301     m_ihciTowEta.push_back(ieta);
0302     m_ihciTowPhi.push_back(iphi);
0303
0304     // double eta = eta_map[ieta];
0305     // double phi = phi_map[iphi];
0306     double eta = towergeom -> get_etacenter(ieta);
0307     double phi = towergeom -> get_phicenter(iphi);
0308
0309     double energy = ihcTowerContainer -> get_tower_at_channel(iter) -> get_energy();
0310     m_ihcTowE.push_back(energy);
0311
0312     if(energy < tow_hcalin_min) continue;
0313
0314     if(first_entry) {
0315       eventOutput << "\"HCALIN\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << energy << "}\n";
0316       first_entry = false;
0317     }
0318     else eventOutput << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << energy << "}\n";
0319   }
0320
0321   if(tower_range == 0) eventOutput << "\"HCALIN\": [{}]\n";
0322   else          eventOutput << "]\n";
0323
0324   T -> Fill();
0325
0326   return Fun4AllReturnCodes::EVENT_OK;
0327 }
0328
0329 //____________________________________________________________________________..
0330 int caloTreeGen::ResetEvent(PHCompositeNode *topNode)
0331 {
0332   m_emcTowPhi.clear();
0333   m_emcTowEta.clear();
0334   m_emcTowE.clear();
0335   m_emcTiming.clear();
0336   m_emciEta.clear();
0337   m_emciPhi.clear();
0338
0339   m_ihcTowE.clear();
0340   m_ihciTowEta.clear();
0341   m_ihciTowPhi.clear();
0342
0343   m_ohcTowE.clear();
0344   m_ohciTowEta.clear();
0345   m_ohciTowPhi.clear();
0346
0347   return Fun4AllReturnCodes::EVENT_OK;
0348 }
0349
0350 //____________________________________________________________________________..
0351 int caloTreeGen::EndRun(const int runnumber)
0352 {
0353   std::cout << "caloTreeGen::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0354   return Fun4AllReturnCodes::EVENT_OK;
0355 }
0356
0357 //____________________________________________________________________________..
0358 int caloTreeGen::End(PHCompositeNode *topNode)
0359 {
0360   std::cout << "caloTreeGen::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0361
0362   out -> cd();
0363   T -> Write();
0364   out -> Close();
0365   delete out;
0366
0367   eventOutput << "},\n \
0368   \"TRACKS\": {\n \
0369   \"INNERTRACKER\": []\n \
0370   }}";
0371   eventOutput.close();
0372   //hm -> dumpHistos(Outfile.c_str(), "UPDATE");
0373   return Fun4AllReturnCodes::EVENT_OK;
0374 }
0375
0376 //____________________________________________________________________________..
0377 int caloTreeGen::Reset(PHCompositeNode *topNode)
0378 {
0379   std::cout << "caloTreeGen::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0380   return Fun4AllReturnCodes::EVENT_OK;
0381 }
0382
0383 //____________________________________________________________________________..
0384 void caloTreeGen::Print(const std::string &what) const
0385 {
0386   std::cout << "caloTreeGen::Print(const std::string &what) const Printing info for " << what << std::endl;
0387 }
0388 // convert from calorimeter key to phi bin
0389 //____________________________________________________________________________..
0390 unsigned int caloTreeGen::getCaloTowerPhiBin(const unsigned int key)
0391 {
0392   unsigned int etabin = key >> 16U;
0393   unsigned int phibin = key - (etabin << 16U);
0394   return phibin;
0395 }
0396
0397 //____________________________________________________________________________..
0398 // convert from calorimeter key to eta bin
0399 unsigned int caloTreeGen::getCaloTowerEtaBin(const unsigned int key)
0400 {
0401   unsigned int etabin = key >> 16U;
0402   return etabin;
0403 }
0404 //____________________________________________________________________________..
0405 float caloTreeGen::getMaxTowerE(RawCluster *cluster, TowerInfoContainer *towerContainer)
0406 {
0407   RawCluster::TowerConstRange towers = cluster -> get_towers();
0408   RawCluster::TowerConstIterator toweriter;
0409
0410   float maxEnergy = 0;
0411   for(toweriter = towers.first; toweriter != towers.second; toweriter++)
0412   {
0413     float towE = toweriter -> second;
0414
0415     if( towE > maxEnergy)  maxEnergy = towE;
0416   }
0417   return maxEnergy;
0418 }
0419 //____________________________________________________________________________..
0420 std::vector<int> caloTreeGen::returnClusterTowEta(RawCluster *cluster, TowerInfoContainer *towerContainer)
0421 {
0422   RawCluster::TowerConstRange towers = cluster -> get_towers();
0423   RawCluster::TowerConstIterator toweriter;
0424
0425   std::vector<int> towerIDsEta;
0426   for(toweriter = towers.first; toweriter != towers.second; toweriter++) towerIDsEta.push_back(RawTowerDefs::decode_index1(toweriter -> first));
0427
0428   return towerIDsEta;
0429 }
0430 //____________________________________________________________________________..
0431 std::vector<int> caloTreeGen::returnClusterTowPhi(RawCluster *cluster, TowerInfoContainer *towerContainer)
0432 {
0433   RawCluster::TowerConstRange towers = cluster -> get_towers();
0434   RawCluster::TowerConstIterator toweriter;
0435
0436   std::vector<int> towerIDsPhi;
0437   for(toweriter = towers.first; toweriter != towers.second; toweriter++) towerIDsPhi.push_back(RawTowerDefs::decode_index2(toweriter -> first));
0438   return towerIDsPhi;
0439 }
0440 //____________________________________________________________________________..
0441 std::vector<float> caloTreeGen::returnClusterTowE(RawCluster *cluster, TowerInfoContainer *towerContainer)
0442 {
0443   RawCluster::TowerConstRange towers = cluster -> get_towers();
0444   RawCluster::TowerConstIterator toweriter;
0445
0446   std::vector<float> towerE;
0447   for(toweriter = towers.first; toweriter != towers.second; toweriter++) towerE.push_back(toweriter -> second);
0448
0449   return towerE;
0450 }