JetSpectrum/src/caloTreeGen.cc

0001 #include "caloTreeGen.h"
0002
0003 // Jet.
0004 #include <jetbase/Jetv1.h>
0005 #include <jetbase/Jetv2.h>
0006 #include <jetbase/JetContainer.h>
0007
0008 // Tower.
0009 #include <calobase/TowerInfo.h>
0010 #include <calobase/TowerInfov1.h>
0011 #include <calobase/TowerInfov2.h>
0012 #include <calobase/TowerInfov3.h>
0013 #include <calobase/TowerInfoContainer.h>
0014 #include <calobase/TowerInfoContainerv1.h>
0015 #include <calobase/TowerInfoContainerv2.h>
0016 #include <calobase/TowerInfoContainerv3.h>
0017 #include <calobase/TowerInfoContainerv4.h>
0018 #include <calobase/TowerInfoDefs.h>
0019
0020 // Vertex.
0021 #include <globalvertex/GlobalVertex.h>
0022 #include <globalvertex/GlobalVertexMap.h>
0023 #include <globalvertex/MbdVertex.h>
0024 #include <globalvertex/MbdVertexMap.h>
0025
0026 // GL1 Information.
0027 #include <ffarawobjects/Gl1Packet.h>
0028
0029 // Fun4All.
0030 #include <fun4all/Fun4AllReturnCodes.h>
0031 #include <fun4all/Fun4AllServer.h>
0032
0033 // Nodes.
0034 #include <phool/PHCompositeNode.h>
0035 #include <phool/getClass.h>
0036 #include <phool/phool.h>
0037 #include <phool/PHIODataNode.h>
0038 #include <phool/PHNode.h>
0039 #include <phool/PHNodeIterator.h>
0040 #include <phool/PHObject.h>
0041
0042 #include <pdbcalbase/PdbParameterMap.h>
0043 #include <phparameter/PHParameters.h>
0044
0045 // General.
0046 #include <cstdint>
0047 #include <iostream>
0048 #include <map>
0049 #include <utility>
0050
0051 class PHCompositeNode;
0052 class TowerInfoContainer;
0053 class caloTreeGen;
0054
0055 caloTreeGen::caloTreeGen(const std::string &name, const std::string &outfilename)
0056   :SubsysReco(name)
0057 {
0058   verbosity = 0;
0059   doRecoJet_radius = {0.4};
0060   foutname = outfilename;
0061   std::cout << "caloTreeGen::caloTreeGen(const std::string &name) Calling ctor" << std::endl;
0062 }
0063
0064 int caloTreeGen::Init(PHCompositeNode * /*topNode*/) {
0065   if (verbosity > 0) std::cout << "Processing initialization: CaloEmulatorTreeMaker::Init(PHCompositeNode *topNode)" << std::endl;
0066   file = new TFile( foutname.c_str(), "RECREATE");
0067   tree = new TTree("ttree","TTree for JES calibration");
0068
0069   // Vertex information.
0070   Initialize_z_vertex();
0071   // Jet information.
0072   for (int ir = 0; ir < static_cast<int>(doRecoJet_radius.size()); ir++) {
0073     Initialize_jet("unsubjet", doRecoJet_radius[ir]);
0074   }
0075   // Trigger information.
0076   Initialize_trigger_vector();
0077
0078   ievent = 0;
0079   return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 int caloTreeGen::process_event(PHCompositeNode *topNode) {
0083   if (verbosity >= 0) {
0084     if (ievent%100 == 0) std::cout << "Processing event " << ievent << std::endl;
0085   }
0086
0087   // Vertex information.
0088   Fill_z_vertex(topNode);
0089   // Jet information.
0090   for (int ir = 0; ir < static_cast<int>(doRecoJet_radius.size()); ir++) {
0091     Fill_recojet(topNode, "unsubjet", "unsubtracted", doRecoJet_radius[ir]);
0092   }
0093   // Trigger information.
0094   Fill_trigger_vector(topNode);
0095
0096   tree->Fill();
0097   ievent++;
0098   return Fun4AllReturnCodes::EVENT_OK;
0099 }
0100
0101 int caloTreeGen::ResetEvent(PHCompositeNode * /*topNode*/) {
0102   if (Verbosity() > 1) std::cout << "Resetting the tree branches" << std::endl;
0103
0104   Reset_jet();
0105
0106   return Fun4AllReturnCodes::EVENT_OK;
0107 }
0108
0109 int caloTreeGen::End(PHCompositeNode * /*topNode*/) {
0110   std::cout << "caloTreeGen::End(PHCompositeNode *topNode) Saving TTree" << std::endl;
0111   std::cout<<"Total events: "<<ievent<<std::endl;
0112   file->cd();
0113   tree->Write();
0114   file->Close();
0115   delete file;
0116   std::cout << "CaloEmulatorTreeMaker complete." << std::endl;
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 ////////// ********** Initialize functions ********** //////////
0121 void caloTreeGen::Initialize_z_vertex() {
0122   tree->Branch("z_vertex", &z_vertex, "z_vertex/F");
0123 }
0124
0125 void caloTreeGen::Initialize_jet(std::string jet_prefix , float radius) {
0126   // Prefix for the jet.
0127   std::string vector_prefix = jet_prefix + "0" + std::to_string((int)(radius * 10)) + "_";
0128   // Jet general information.
0129   std::string jet_e_vectorname = vector_prefix + "e"; jet_e_map[jet_e_vectorname] = std::vector<float>(); tree->Branch(jet_e_vectorname.c_str(), &jet_e_map[jet_e_vectorname]);
0130   std::string jet_et_vectorname = vector_prefix + "et"; jet_et_map[jet_et_vectorname] = std::vector<float>(); tree->Branch(jet_et_vectorname.c_str(), &jet_et_map[jet_et_vectorname]);
0131   std::string jet_pt_vectorname = vector_prefix + "pt"; jet_pt_map[jet_pt_vectorname] = std::vector<float>(); tree->Branch(jet_pt_vectorname.c_str(), &jet_pt_map[jet_pt_vectorname]);
0132   std::string jet_eta_vectorname = vector_prefix + "eta"; jet_eta_map[jet_eta_vectorname] = std::vector<float>(); tree->Branch(jet_eta_vectorname.c_str(), &jet_eta_map[jet_eta_vectorname]);
0133   std::string jet_phi_vectorname = vector_prefix + "phi"; jet_phi_map[jet_phi_vectorname] = std::vector<float>(); tree->Branch(jet_phi_vectorname.c_str(), &jet_phi_map[jet_phi_vectorname]);
0134   std::string jet_emcal_calo_e_vectorname = vector_prefix + "emcal_calo_e"; jet_calo_e_map[jet_emcal_calo_e_vectorname] = std::vector<float>(); tree->Branch(jet_emcal_calo_e_vectorname.c_str(), &jet_calo_e_map[jet_emcal_calo_e_vectorname]); // iHCaL tower information.
0135   std::string jet_ihcal_calo_e_vectorname = vector_prefix + "ihcal_calo_e"; jet_calo_e_map[jet_ihcal_calo_e_vectorname] = std::vector<float>(); tree->Branch(jet_ihcal_calo_e_vectorname.c_str(), &jet_calo_e_map[jet_ihcal_calo_e_vectorname]);
0136   std::string jet_ohcal_calo_e_vectorname = vector_prefix + "ohcal_calo_e"; jet_calo_e_map[jet_ohcal_calo_e_vectorname] = std::vector<float>(); tree->Branch(jet_ohcal_calo_e_vectorname.c_str(), &jet_calo_e_map[jet_ohcal_calo_e_vectorname]);
0137 }
0138
0139 void caloTreeGen::Initialize_trigger_vector() {
0140   tree->Branch("gl1_trigger_vector_live", gl1_trigger_vector_live, "gl1_trigger_vector_live[64]/I");
0141   tree->Branch("gl1_trigger_vector_scaled", gl1_trigger_vector_scaled, "gl1_trigger_vector_scaled[64]/I");
0142   tree->Branch("gl1_trigger_scalar_live", gl1_trigger_scalar_live, "gl1_trigger_scalar_live[64]/l");
0143   tree->Branch("gl1_trigger_scalar_scaled", gl1_trigger_scalar_scaled, "gl1_trigger_scalar_scaled[64]/l");
0144 }
0145
0146 ////////// ********** Fill functions ********** //////////
0147 void caloTreeGen::Fill_z_vertex(PHCompositeNode *topNode) {
0148   z_vertex = -999;
0149   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0150   if (vertexmap) {
0151     if (!vertexmap->empty()) {
0152       GlobalVertex *vtx = vertexmap->begin()->second;
0153       if (vtx) {
0154         z_vertex = vtx->get_z();
0155       }
0156     } else {
0157       std::cout << "GlobalVertexMap is empty" << std::endl;
0158     }
0159   } else {
0160     std::cout << "GlobalVertexMap is missing" << std::endl;
0161   }
0162 }
0163
0164 void caloTreeGen::Fill_recojet(PHCompositeNode *topNode, std::string jet_prefix, std::string node_prefix, float radius) {
0165   // Set vector name prefix for the jet map.
0166   std::string vector_prefix = jet_prefix + "0" + std::to_string((int)(radius * 10)) + "_";
0167   // Get the jet container.
0168   std::string nodename = "AntiKt_" + node_prefix + "_r0" + std::to_string((int)(radius * 10));
0169   if (node_prefix == "subtracted") {
0170     nodename = "AntiKt_Tower_r0" + std::to_string((int)(radius * 10)) + "_Sub1";
0171   }
0172   JetContainer* _jets = findNode::getClass<JetContainer>(topNode, nodename);
0173   if (!_jets) std::cout << "JetContainer for " << nodename << " is missing" << std::endl;
0174   // Get the towerinfo container.
0175   TowerInfoContainer *_emcal_towers = nullptr;
0176   TowerInfoContainer *_ihcal_towers = nullptr;
0177   TowerInfoContainer *_ohcal_towers = nullptr;
0178   if (node_prefix == "unsubtracted") {
0179     _emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0180     _ihcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0181     _ohcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0182   } else {
0183     _emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0184     _ihcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0185     _ohcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0186   }
0187   if (!_emcal_towers || !_ihcal_towers || !_ohcal_towers) std::cout << "TowerInfoContainer is missing" << std::endl;
0188   // Fill the jet information.
0189   for (auto _jet : *_jets) {
0190     if (_jet->get_pt() < 1) continue;
0191     jet_e_map[vector_prefix+"e"].push_back(_jet->get_e());
0192     jet_et_map[vector_prefix+"et"].push_back(_jet->get_et());
0193     jet_pt_map[vector_prefix+"pt"].push_back(_jet->get_pt());
0194     jet_eta_map[vector_prefix+"eta"].push_back(_jet->get_eta());
0195     jet_phi_map[vector_prefix+"phi"].push_back(_jet->get_phi());
0196
0197     float emcal_calo_e = 0;
0198     float ihcal_calo_e = 0;
0199     float ohcal_calo_e = 0;
0200     for (auto comp: _jet->get_comp_vec()) {
0201         unsigned int channel = comp.second;
0202       if (comp.first == 14 || comp.first == 29 || comp.first == 25 || comp.first == 28) {
0203         TowerInfo *_tower = _emcal_towers->get_tower_at_channel(channel);
0204           float jet_tower_e = _tower->get_energy();
0205         emcal_calo_e += jet_tower_e;
0206         } else if (comp.first == 15 ||  comp.first == 30 || comp.first == 26) {
0207         TowerInfo *_tower = _ihcal_towers->get_tower_at_channel(channel);
0208           float jet_tower_e = _tower->get_energy();
0209         ihcal_calo_e += jet_tower_e;
0210         } else if (comp.first == 16 || comp.first == 31 || comp.first == 27) {
0211         TowerInfo *_tower = _ohcal_towers->get_tower_at_channel(channel);
0212           float jet_tower_e = _tower->get_energy();
0213         ohcal_calo_e += jet_tower_e;
0214       }
0215       }
0216     jet_calo_e_map[vector_prefix+"emcal_calo_e"].push_back(emcal_calo_e);
0217     jet_calo_e_map[vector_prefix+"ihcal_calo_e"].push_back(ihcal_calo_e);
0218     jet_calo_e_map[vector_prefix+"ohcal_calo_e"].push_back(ohcal_calo_e);
0219   }
0220 }
0221
0222 void caloTreeGen::Fill_trigger_vector(PHCompositeNode *topNode) {
0223   Gl1Packet *_gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0224   if (_gl1PacketInfo) {
0225     uint64_t livetriggervec = _gl1PacketInfo->lValue(0, "TriggerVector");
0226     uint64_t scaledtriggervec = _gl1PacketInfo->lValue(0, "ScaledVector");
0227
0228     for (int iv = 0; iv < n_trigger_vec; ++iv) {
0229       gl1_trigger_vector_live[iv] = ((livetriggervec & 0x1U) == 0x1U);
0230       gl1_trigger_vector_scaled[iv] = ((scaledtriggervec & 0x1U) == 0x1U);
0231       livetriggervec = (livetriggervec >> 1U) & 0xffffffffU;
0232       scaledtriggervec = (scaledtriggervec >> 1U) & 0xffffffffU;
0233
0234         gl1_trigger_scalar_live[iv] = _gl1PacketInfo->lValue(iv, 1);
0235       gl1_trigger_scalar_scaled[iv] = _gl1PacketInfo->lValue(iv, 2);
0236       }
0237   } else {
0238     std::cout << "GL1Packet is missing" << std::endl;
0239   }
0240 }
0241
0242 ////////// ********** Reset functions ********** //////////
0243 void caloTreeGen::Reset_jet() {
0244   for (auto& entry : jet_e_map) entry.second.clear();
0245   for (auto& entry : jet_et_map) entry.second.clear();
0246   for (auto& entry : jet_pt_map) entry.second.clear();
0247   for (auto& entry : jet_eta_map) entry.second.clear();
0248   for (auto& entry : jet_phi_map) entry.second.clear();
0249   for (auto& entry : jet_calo_e_map) entry.second.clear();
0250 }