analysis/EventDisplayMaker/displayEvents.cc

0001 #include "displayEvents.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <trackbase/InttDefs.h>              // for getLadderPhiId, getLad...
0006 #include <trackbase/MvtxDefs.h>              // for getChipId, getStaveId
0007 #include <trackbase/TpcDefs.h>               // for getSectorId, getSide
0008 #include <trackbase/TrkrDefs.h>              // for getLayer, getTrkrId
0009
0010 #include <phool/getClass.h>
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/recoConsts.h>
0013
0014 #include <TDatabasePDG.h>
0015
0016 #include <nlohmann/json.hpp>
0017
0018 using json = nlohmann::json;
0019
0020 //____________________________________________________________________________..
0021 displayEvents::displayEvents(const std::string &name)
0022  : SubsysReco(name)
0023 {
0024 }
0025
0026 //____________________________________________________________________________..
0027 displayEvents::~displayEvents()
0028 {
0029 }
0030
0031 //____________________________________________________________________________..
0032 int displayEvents::Init(PHCompositeNode *topNode)
0033 {
0034   load_nodes(topNode);
0035
0036   triggeranalyzer = new TriggerAnalyzer();
0037
0038   recoConsts *rc = recoConsts::instance();
0039   m_runNumber = rc->get_IntFlag("RUNNUMBER");
0040
0041   m_run_date = getDate();
0042
0043   return Fun4AllReturnCodes::EVENT_OK;
0044 }
0045
0046 //____________________________________________________________________________..
0047 int displayEvents::process_event(PHCompositeNode *topNode)
0048 {
0049   if (counter >= m_max_displays) return Fun4AllReturnCodes::EVENT_OK; //Made the max number of displays
0050
0051   load_nodes(topNode);
0052
0053   for (auto &vtx_iter : *m_vertexMap)
0054   {
0055     if (counter >= m_max_displays) break; //Made the max number of displays
0056
0057     all_tracks.clear();
0058
0059     m_vertex = vtx_iter.second;
0060
0061     //Get all tracks associated to this vertex
0062     for (SvtxVertex::TrackIter all_track_iter =  m_vertex->begin_tracks(); all_track_iter != m_vertex->end_tracks(); ++all_track_iter)
0063     {
0064       m_track = m_trackMap->get(*all_track_iter);
0065       all_tracks.push_back(m_track);
0066     }
0067
0068     if (loadCalos && (hasEMCal || hasIHCal || hasOHCal) && m_track->get_crossing() != 0) continue;
0069
0070     triggeranalyzer->decodeTriggers(topNode);
0071
0072     //Now fill the json file
0073     json data;
0074
0075     uint64_t m_bco = gl1packet->lValue(0, "BCO") + all_tracks[0]->get_crossing();
0076
0077     //Check if we have a list of BCOs we want a plot for, and this is one of them
0078     if (m_bco_list.size() != 0 && std::find(std::begin(m_bco_list), std::end(m_bco_list), m_bco) == std::end(m_bco_list)) continue;
0079
0080     std::string date_run_stamp = m_run_date + ", Run " + std::to_string(m_runNumber);
0081     std::string bco_stamp = "BCO: " + std::to_string(m_bco);
0082
0083     data[eventName]["evtid"] = 1;
0084     data[eventName]["runid"] = m_runNumber;
0085     data[eventName]["type"] = "Collision";
0086     data[eventName]["s_nn"] = 0;
0087     data[eventName]["B"] = 1.38;
0088     data[eventName]["pv"] = {m_vertex->get_x(), m_vertex->get_y(), m_vertex->get_z()};
0089
0090     if (m_approved)
0091     {
0092       data[eventName]["runstats"] = {"sPHENIX Tracking", date_run_stamp, m_decay_tag};
0093     }
0094     else
0095     {
0096       data[eventName]["runstats"] = {"sPHENIX Internal", date_run_stamp, bco_stamp, m_decay_tag};
0097     }
0098
0099     for (auto& hitsMetaSetup : hitsNameVector)
0100     {
0101       data[metadataName][hitsName][hitsMetaSetup]["type"] = "3D";
0102       data[metadataName][hitsName][hitsMetaSetup]["options"]["size"] = 1;
0103       data[metadataName][hitsName][hitsMetaSetup]["options"]["color"] = pidToColourMap[hitsMetaSetup];
0104     }
0105
0106     data[metadataName][trackName][allTrackName]["width"] = 10;
0107
0108     data[trackName]["B"] = 0.000014;
0109
0110     getJSONdata(all_tracks, data);
0111
0112     if (loadCalos)
0113     {
0114       data[metadataName][hitsName][emcalName] = caloMetadata(emcalName);
0115       data[metadataName][hitsName][ihcalName] = caloMetadata(ihcalName);
0116       data[metadataName][hitsName][ohcalName] = caloMetadata(ohcalName);
0117
0118       if (hasEMCal) getJSONcalo(data, emcalName, towersEM, geomEMCal);
0119       if (hasIHCal) getJSONcalo(data, ihcalName, towersIH, geoIHCal);
0120       if (hasOHCal) getJSONcalo(data, ohcalName, towersOH, geoOHCal);
0121     }
0122
0123     std::string m_output_file = m_evt_display_path + "/EvtDisplay_" + std::to_string(m_runNumber) + "_" + std::to_string(m_bco) + ".json";
0124
0125     json_output.open(m_output_file);
0126     json_output << data.dump(2);
0127     json_output.close();
0128
0129     ++counter;
0130   }
0131
0132   return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134
0135 //____________________________________________________________________________..
0136 int displayEvents::End(PHCompositeNode *topNode)
0137 {
0138   return Fun4AllReturnCodes::EVENT_OK;
0139 }
0140
0141 int displayEvents::load_nodes(PHCompositeNode *topNode)
0142 {
0143   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0144   if (!m_vertexMap)
0145   {
0146     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing SvtxVertexMap" << std::endl;
0147     return Fun4AllReturnCodes::ABORTRUN;
0148   }
0149
0150   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0151   if (!m_trackMap)
0152   {
0153     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing SvtxTrackMap" << std::endl;
0154     return Fun4AllReturnCodes::ABORTRUN;
0155   }
0156
0157   gl1packet = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0158   if (!gl1packet)
0159   {
0160     gl1packet = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0161     if (!gl1packet)
0162     {
0163       std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing Gl1Packet" << std::endl;
0164       return Fun4AllReturnCodes::ABORTRUN;
0165     }
0166   }
0167
0168   triggerruninfo = findNode::getClass<TriggerRunInfo>(topNode, "TriggerRunInfo");
0169   if (!triggerruninfo)
0170   {
0171     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing TriggerRunInfo" << std::endl;
0172     return Fun4AllReturnCodes::ABORTRUN;
0173   }
0174
0175   dst_clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0176   if (!dst_clustermap)
0177   {
0178     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing TRKR_CLUSTER" << std::endl;
0179     return Fun4AllReturnCodes::ABORTRUN;
0180   }
0181
0182   geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0183   if (!geometry)
0184   {
0185     std::cout << __PRETTY_FUNCTION__ << " Fatal Error : missing ActsGeometry" << std::endl;
0186     return Fun4AllReturnCodes::ABORTRUN;
0187   }
0188
0189   if (loadCalos)
0190   {
0191     if (!towersEM) towersEM = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0192     if (!towersEM) towersEM = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_CEMC");
0193     if (!geomEMCal) geomEMCal = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0194
0195     if (!towersIH) towersIH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0196     if (!towersIH) towersIH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALIN");
0197     if (!geoIHCal) geoIHCal = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0198
0199     if (!towersOH) towersOH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0200     if (!towersOH) towersOH = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALOUT");
0201     if (!geoOHCal) geoOHCal = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0202
0203     if (!towersEM || !geomEMCal) hasEMCal = false;
0204     if (!towersIH || !geoIHCal) hasIHCal = false;
0205     if (!towersOH || !geoOHCal) hasOHCal = false;
0206   }
0207
0208   return Fun4AllReturnCodes::EVENT_OK;
0209 }
0210
0211 std::string displayEvents::getParticleName(int ID)
0212 {
0213   return TDatabasePDG::Instance()->GetParticle(ID)->GetName();
0214 }
0215
0216 SvtxTrack *displayEvents::getTrack(unsigned int track_id, SvtxTrackMap *trackmap)
0217 {
0218   SvtxTrack *matched_track = nullptr;
0219
0220   for (auto &iter : *trackmap)
0221   {
0222     if (iter.first == track_id)
0223     {
0224       matched_track = iter.second;
0225     }
0226   }
0227
0228   return matched_track;
0229 }
0230
0231 void displayEvents::getJSONcalo(json &jsonData, std::string calo, TowerInfoContainer *towerCont, RawTowerGeomContainer *geom)
0232 {
0233   json calos;
0234
0235   for (size_t channel = 0; channel < towerCont->size(); channel++)
0236   {
0237     TowerInfo* tower = towerCont->get_tower_at_channel(channel);
0238     float raw_energy = tower->get_energy();
0239     float thresholdCut = caloThreshold[calo];
0240     if (raw_energy < thresholdCut)
0241     {
0242       continue;
0243     }
0244
0245     unsigned int towerkey = towerCont->encode_key(channel);
0246     int ieta = towerCont->getTowerEtaBin(towerkey);
0247     int iphi = towerCont->getTowerPhiBin(towerkey);
0248     calos["e"] = raw_energy;
0249     calos["phi"] = geom->get_phicenter(iphi);
0250     calos["eta"] = geom->get_etacenter(ieta);
0251     calos["event"] = 42;
0252
0253     jsonData[hitsName][calo] += calos;
0254   }
0255 }
0256
0257 void displayEvents::getJSONdata(std::vector<SvtxTrack*> tracks, json &jsonData)
0258 {
0259   json tracksJson, clustersJson;
0260
0261   int trackCounter = 0; //Needed for getting KFP objects
0262
0263   for (auto track : tracks)
0264   {
0265     float length = 0;
0266
0267     TrackSeed *silseed = track->get_silicon_seed();
0268     if (silseed)
0269     {
0270       for (auto cluster_iter = silseed->begin_cluster_keys(); cluster_iter != silseed->end_cluster_keys(); ++cluster_iter)
0271       {
0272         uint64_t clusterKey = *cluster_iter;
0273
0274         TrkrCluster *cluster = dst_clustermap->findCluster(clusterKey);
0275         auto global = geometry->getGlobalPosition(clusterKey, cluster);
0276
0277         clustersJson["x"] = (float) global.x();
0278         clustersJson["y"] = (float) global.y();
0279         clustersJson["z"] = (float) global.z();
0280         clustersJson["e"] = 0;
0281
0282         if (m_dont_show_track_clus == false) jsonData[hitsName][allTrackName] += clustersJson;
0283       }
0284     }
0285
0286     for (auto state_iter = track->begin_states();
0287          state_iter != track->end_states();
0288          ++state_iter)
0289     {
0290       SvtxTrackState* tstate = state_iter->second;
0291       if (tstate->get_pathlength() == 0) continue; //The first track state is an extrapolation so has no cluster
0292
0293       uint64_t clusterKey = tstate->get_cluskey();
0294
0295       TrkrCluster *cluster = dst_clustermap->findCluster(clusterKey);
0296       auto global = geometry->getGlobalPosition(clusterKey, cluster);
0297
0298       uint8_t id = TrkrDefs::getTrkrId(clusterKey);
0299
0300       if (id == TrkrDefs::mvtxId || id == TrkrDefs::inttId) continue; //Due to track state issue, use seeds to populate silicon clusters
0301
0302       clustersJson["x"] = id == TrkrDefs::tpcId ? tstate->get_x() : (float) global.x();
0303       clustersJson["y"] = id == TrkrDefs::tpcId ? tstate->get_y() : (float) global.y();
0304       clustersJson["z"] = id == TrkrDefs::tpcId ? tstate->get_z() : (float) global.z();
0305       clustersJson["e"] = 0;
0306
0307       if (m_dont_show_track_clus == false) jsonData[hitsName][allTrackName] += clustersJson;
0308
0309       length = std::max(tstate->get_pathlength(), length);
0310
0311     } //End of state loop
0312
0313     std::vector<float> trackPosition;
0314     std::vector<float> trackMomentum;
0315
0316     trackPosition = {track->get_x(), track->get_y(), track->get_z()};
0317     trackMomentum = {track->get_px(), track->get_py(), track->get_pz()};
0318
0319     tracksJson["pxyz"] = trackMomentum;
0320     tracksJson["xyz"] = trackPosition;
0321     tracksJson["trk_color"] = allColour;
0322     tracksJson["nh"] = 60;
0323     tracksJson["l"] = length;
0324     tracksJson["q"] = track->get_charge();
0325
0326     jsonData[trackName][allTrackName] += tracksJson;
0327
0328     ++trackCounter;
0329
0330   } //End of track loop
0331 }
0332
0333 //https://stackoverflow.com/questions/997946/how-to-get-current-time-and-date-in-c
0334 std::string displayEvents::getDate()
0335 {
0336     std::time_t t = std::time(0);   // get time now
0337     std::tm* now = std::localtime(&t);
0338
0339     std::stringstream date;
0340     date << (now->tm_year + 1900) << '-'
0341          << (now->tm_mon + 1) << '-'
0342          <<  now->tm_mday;
0343     return date.str();
0344 }
0345
0346 bool displayEvents::isInRange(float min, float value, float max)
0347 {
0348   return min <= value && value <= max;
0349 }
0350
0351 json displayEvents::caloMetadata(std::string type)
0352 {
0353   json metaData, options, minmax;
0354
0355   //Default for EMCal
0356   float rMin = 90;
0357   float rMax = 136.1;
0358   int color = emcalColor;
0359   float scaleMax = 2.71998;
0360
0361   if (type != emcalName)
0362   {
0363     rMin = type == ihcalName ? 147.27 : 183.3;
0364     rMax = type == ihcalName ? 175.0 : 348.634;
0365     color = type == ihcalName ? ihcalColor : ohcalColor;
0366     scaleMax = type == ihcalName ? 0.569264 : 2.57763;
0367   }
0368
0369   metaData["type"] = "PROJECTIVE";
0370
0371   options["rmin"] = rMin;
0372   options["rmax"] = rMax;
0373   options["deta"] = 0.025;
0374   options["dphi"] = 0.025;
0375   options["color"] = color;
0376   options["transparent"] = 0.6;
0377   minmax["min"] = 0.05;
0378   minmax["max"] = scaleMax;
0379   options["scaleminmax"] = minmax;
0380
0381   metaData["options"] = options;
0382
0383   return metaData;
0384 }