HF-Jet/HFMLTrigger/HFMLTriggerOccupancy.C

0001 #include "HFMLTriggerOccupancy.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <pdbcalbase/PdbParameterMap.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/PHTimeServer.h>
0007 #include <phool/PHTimer.h>
0008 #include <phool/getClass.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011
0012 #include <g4main/PHG4Hit.h>
0013 #include <g4main/PHG4Particle.h>
0014 #include <g4main/PHG4VtxPoint.h>
0015 #include <g4main/PHG4TruthInfoContainer.h>
0016
0017 #include <g4mvtx/PHG4CylinderGeom_MVTX.h>
0018
0019 #include <g4detectors/PHG4Cell.h>
0020 #include <g4detectors/PHG4CellContainer.h>
0021 #include <g4detectors/PHG4CylinderCellGeom.h>
0022 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0023 #include <g4detectors/PHG4CylinderGeom.h>
0024 #include <g4detectors/PHG4CylinderGeomContainer.h>
0025
0026 #include <trackbase_historic/SvtxCluster.h>
0027 #include <trackbase_historic/SvtxClusterMap.h>
0028 #include <trackbase_historic/SvtxHit.h>
0029 #include <trackbase_historic/SvtxHitMap.h>
0030 #include <trackbase_historic/SvtxTrack.h>
0031 #include <trackbase_historic/SvtxTrackMap.h>
0032 #include <trackbase_historic/SvtxVertex.h>
0033 #include <trackbase_historic/SvtxVertexMap.h>
0034
0035 #include <g4eval/SvtxEvalStack.h>
0036
0037 #include <HepMC/GenEvent.h>
0038 #include <HepMC/GenVertex.h>
0039 #include <phhepmc/PHHepMCGenEvent.h>
0040 #include <phhepmc/PHHepMCGenEventMap.h>
0041
0042 #include <TDatabasePDG.h>
0043 #include <TFile.h>
0044 #include <TH1D.h>
0045 #include <TH1F.h>
0046 #include <TH2F.h>
0047 #include <TH3F.h>
0048 #include <TLorentzVector.h>
0049 #include <TString.h>
0050 #include <TTree.h>
0051
0052 #include <rapidjson/document.h>
0053 #include <rapidjson/ostreamwrapper.h>
0054 #include <rapidjson/prettywriter.h>
0055
0056 #include <boost/format.hpp>
0057 #include <boost/property_tree/json_parser.hpp>
0058 #include <boost/property_tree/ptree.hpp>
0059 #include <boost/property_tree/xml_parser.hpp>
0060
0061 #include <cassert>
0062 #include <cmath>
0063 #include <cstddef>
0064 #include <iostream>
0065
0066 using namespace std;
0067
0068 HFMLTriggerOccupancy::HFMLTriggerOccupancy(std::string filename)
0069   : SubsysReco("HFMLTriggerOccupancy")
0070   , _ievent(0)
0071   , _f(nullptr)
0072   , _eta_min(-1)
0073   , _eta_max(1)
0074   , _embedding_id(1)
0075   , _nlayers_maps(3)
0076   , _svtxevalstack(nullptr)
0077   , m_hitMap(nullptr)
0078   , m_Geoms(nullptr)
0079   , m_truthInfo(nullptr)
0080   , m_Flags(nullptr)
0081   , m_hNorm(nullptr)
0082   , m_hNChEta(nullptr)
0083   , m_hVertexZ(nullptr)
0084   , m_hitStaveLayer(nullptr)
0085   , m_hitModuleHalfStave(nullptr)
0086   , m_hitChipModule(nullptr)
0087   , m_hitLayerMap(nullptr)
0088   , m_hitPixelPhiMap(nullptr)
0089   , m_hitPixelPhiMapHL(nullptr)
0090   , m_hitPixelZMap(nullptr)
0091   , m_Multiplicity(nullptr)
0092   , m_LayerMultiplicity(nullptr)
0093   , m_LayerChipMultiplicity(nullptr)
0094 {
0095   _foutname = filename;
0096 }
0097
0098 int HFMLTriggerOccupancy::Init(PHCompositeNode* topNode)
0099 {
0100   _ievent = 0;
0101
0102   _f = new TFile((_foutname + string(".root")).c_str(), "RECREATE");
0103
0104   //  m_jsonOut.open((_foutname + string(".json")).c_str(), fstream::out);
0105   //
0106   //  m_jsonOut << "{" << endl;
0107   //  m_jsonOut << "\"Events\" : [" << endl;
0108
0109   //  _h2 = new TH2D("h2", "", 100, 0, 100.0, 40, -2, +2);
0110   //  _h2_b = new TH2D("h2_b", "", 100, 0, 100.0, 40, -2, +2);
0111   //  _h2_c = new TH2D("h2_c", "", 100, 0, 100.0, 40, -2, +2);
0112   //  _h2all = new TH2D("h2all", "", 100, 0, 100.0, 40, -2, +2);
0113
0114   m_hNorm = new TH1F("hNormalization",  //
0115                      "Normalization;Items;Summed quantity", 10, .5, 10.5);
0116   int i = 1;
0117   m_hNorm->GetXaxis()->SetBinLabel(i++, "IntegratedLumi");
0118   m_hNorm->GetXaxis()->SetBinLabel(i++, "Event");
0119   m_hNorm->GetXaxis()->LabelsOption("v");
0120
0121   m_hNChEta = new TH1F("hNChEta",  //
0122                        "Charged particle #eta distribution;#eta;Count",
0123                        1000, -5, 5);
0124   m_hVertexZ = new TH1F("hVertexZ",  //
0125                         "Vertex z distribution;z [cm];Count",
0126                         1000, -200, 200);
0127
0128   m_hitLayerMap = new TH3F("hitLayerMap", "hitLayerMap", 600, -10, 10, 600, -10, 10, 10, -.5, 9.5);
0129   m_hitLayerMap->SetTitle("hitLayerMap;x [mm];y [mm];Half Layers");
0130
0131 //  m_hitPixelPhiMap = new TH3F("hitPixelPhiMap", "hitPixelPhiMap", 16000, -.5, 16000 - .5, 600, -M_PI, M_PI, 10, -.5, 9.5);
0132 //  m_hitPixelPhiMap->SetTitle("hitPixelPhiMap;PixelPhiIndex in layer;phi [rad];Half Layers Index");
0133 //  m_hitPixelPhiMapHL = new TH3F("hitPixelPhiMapHL", "hitPixelPhiMapHL", 16000, -.5, 16000 - .5, 600, -M_PI, M_PI, 10, -.5, 9.5);
0134 //  m_hitPixelPhiMapHL->SetTitle("hitPixelPhiMap;PixelPhiIndex in half layer;phi [rad];Half Layers Index");
0135 //  m_hitPixelZMap = new TH3F("hitPixelZMap", "hitPixelZMap", 16000, -.5, 16000 - .5, 600, 15, 15, 10, -.5, 9.5);
0136 //  m_hitPixelZMap->SetTitle("hitPixelZMap;hitPixelZMap;z [cm];Half Layers");
0137
0138   m_hitStaveLayer = new TH2F("hitStaveLayer", "hitStaveLayer", 100, -.5, 100 - .5, 10, -.5, 9.5);
0139   m_hitStaveLayer->SetTitle("hitStaveLayer;Stave index;Half Layers");
0140   m_hitModuleHalfStave = new TH2F("hitModuleHalfStave", "hitModuleHalfStave", 100, -.5, 100 - .5, 10, -.5, 9.5);
0141   m_hitModuleHalfStave->SetTitle("hitModuleHalfStave;Module index;Half Stave");
0142   m_hitChipModule = new TH2F("hitChipModule", "hitChipModule", 100, -.5, 100 - .5, 10, -.5, 9.5);
0143   m_hitChipModule->SetTitle("hitChipModule;Chip;Module");
0144
0145   m_Multiplicity = new TH1F("Multiplicity", "Multiplicity", 10000, -.5, 10000 - .5);
0146   m_Multiplicity->SetTitle("Multiplicity;Chip Multiplicity;Count");
0147
0148   m_LayerMultiplicity = new TH2F("LayerMultiplicity", "LayerMultiplicity", 3, -.5, 2.5, 1000, -.5, 1000 - .5);
0149   m_LayerMultiplicity->SetTitle("LayerMultiplicity;Layer;Chip Multiplicity");
0150
0151   m_LayerChipMultiplicity = new TH3F("LayerChipMultiplicity", "LayerChipMultiplicity", 3, -.5, 2.5, kNCHip, -.5, kNCHip - .5, 1000, -.5, 1000 - .5);
0152   m_LayerChipMultiplicity->SetTitle("LayerChipMultiplicity;Layer;Chip");
0153
0154   return Fun4AllReturnCodes::EVENT_OK;
0155 }
0156
0157 int HFMLTriggerOccupancy::InitRun(PHCompositeNode* topNode)
0158 {
0159   m_hitMap = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
0160   if (!m_hitMap)
0161   {
0162     std::cout << PHWHERE << "ERROR: Can't find node SvtxHitMap" << std::endl;
0163     return Fun4AllReturnCodes::ABORTRUN;
0164   }
0165   m_Geoms =
0166       findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0167   if (!m_Geoms)
0168   {
0169     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_MVTX" << std::endl;
0170     return Fun4AllReturnCodes::ABORTRUN;
0171   }
0172
0173   m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0174   if (!m_truthInfo)
0175   {
0176     std::cout << PHWHERE << "ERROR: Can't find node G4TruthInfo" << std::endl;
0177     return Fun4AllReturnCodes::ABORTRUN;
0178   }
0179
0180   m_Flags = findNode::getClass<PdbParameterMap>(topNode, "HFMLTrigger_HepMCTriggerFlags");
0181   if (!m_Flags)
0182   {
0183     cout << "HFMLTriggerOccupancy::InitRun - WARNING - missing HFMLTrigger_HepMCTriggerFlags" << endl;
0184   }
0185
0186   return Fun4AllReturnCodes::EVENT_OK;
0187 }
0188
0189 int HFMLTriggerOccupancy::process_event(PHCompositeNode* topNode)
0190 {
0191   if (!_svtxevalstack)
0192   {
0193     _svtxevalstack = new SvtxEvalStack(topNode);
0194     _svtxevalstack->set_strict(0);
0195     _svtxevalstack->set_verbosity(Verbosity() + 1);
0196   }
0197   else
0198   {
0199     _svtxevalstack->next_event(topNode);
0200   }
0201
0202   //  SvtxVertexEval* vertexeval = _svtxevalstack->get_vertex_eval();
0203   //  SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
0204   //  SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
0205   SvtxHitEval* hiteval = _svtxevalstack->get_hit_eval();
0206   //    SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
0207
0208   //  PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0209   //  if (!geneventmap)
0210   //  {
0211   //    std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
0212   //    return Fun4AllReturnCodes::ABORTRUN;
0213   //  }
0214   //
0215   //  PHHepMCGenEvent* genevt = geneventmap->get(_embedding_id);
0216   //  if (!genevt)
0217   //  {
0218   //    std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
0219   //    std::cout << ". Print PHHepMCGenEventMap:";
0220   //    geneventmap->identify();
0221   //    return Fun4AllReturnCodes::ABORTRUN;
0222   //  }
0223   //
0224   //  HepMC::GenEvent* theEvent = genevt->getEvent();
0225   //  assert(theEvent);
0226   //  if (Verbosity())
0227   //  {
0228   //    cout << "HFMLTriggerOccupancy::process_event - process HepMC::GenEvent with signal_process_id = "
0229   //         << theEvent->signal_process_id();
0230   //    if (theEvent->signal_process_vertex())
0231   //    {
0232   //      cout << " and signal_process_vertex : ";
0233   //      theEvent->signal_process_vertex()->print();
0234   //    }
0235   //    cout << "  - Event record:" << endl;
0236   //    theEvent->print();
0237   //  }
0238
0239   assert(m_hNorm);
0240   m_hNorm->Fill("Event", 1.);
0241
0242   assert(m_truthInfo);
0243   assert(m_hNChEta);
0244   assert(m_hVertexZ);
0245
0246   const PHG4VtxPoint* vtx =
0247   m_truthInfo->GetPrimaryVtx(m_truthInfo->GetPrimaryVertexIndex());
0248   if (vtx)
0249   {
0250     m_hVertexZ -> Fill(vtx->get_z());
0251   }
0252
0253   const auto primary_range =
0254       m_truthInfo->GetPrimaryParticleRange();
0255   for (auto particle_iter =
0256            primary_range.first;
0257        particle_iter != primary_range.second;
0258        ++particle_iter)
0259   {
0260     const PHG4Particle *p = particle_iter->second;
0261     assert(p);
0262     TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(
0263         p->get_pid());
0264     assert(pdg_p);
0265     if (fabs(pdg_p->Charge()) > 0)
0266     {
0267       TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
0268       if (pvec.Perp2() > 0)
0269       {
0270         assert(m_hNChEta);
0271         m_hNChEta->Fill(pvec.PseudoRapidity());
0272       }
0273     }
0274   }  //          if (_load_all_particle) else
0275
0276
0277   // chip counting
0278   vector<vector<vector<int> > >
0279       multiplicity_vec(_nlayers_maps);
0280
0281   //
0282   for (unsigned int layer = 0; layer < _nlayers_maps; ++layer)
0283   {
0284     PHG4CylinderGeom_MVTX* geom = dynamic_cast<PHG4CylinderGeom_MVTX*>(m_Geoms->GetLayerGeom(layer));
0285     assert(geom);
0286
0287     const int n_stave = geom->get_N_staves();
0288
0289     assert(multiplicity_vec[layer].empty());
0290     multiplicity_vec[layer].resize(n_stave, vector<int>(kNCHip, 0));
0291   }
0292
0293   set<unsigned int> mapsHits;  //internal consistency check for later stages of truth tracks
0294   assert(m_hitMap);
0295   for (SvtxHitMap::Iter iter = m_hitMap->begin();
0296        iter != m_hitMap->end();
0297        ++iter)
0298   {
0299     SvtxHit* hit = iter->second;
0300     assert(hit);
0301     unsigned int layer = hit->get_layer();
0302     if (layer < _nlayers_maps)
0303     {
0304       unsigned int hitID = hit->get_id();
0305       mapsHits.insert(hitID);
0306
0307       PHG4Cell* cell = hiteval->get_cell(hit);
0308       assert(cell);
0309       PHG4CylinderGeom_MVTX* geom = dynamic_cast<PHG4CylinderGeom_MVTX*>(m_Geoms->GetLayerGeom(layer));
0310       assert(geom);
0311
0312       TVector3 local_coords = geom->get_local_coords_from_pixel(cell->get_pixel_index());
0313       TVector3 world_coords = geom->get_world_from_local_coords(cell->get_stave_index(), cell->get_half_stave_index(), cell->get_module_index(), cell->get_chip_index(), local_coords);
0314
0315       unsigned int pixel_x(cell->get_pixel_index() % geom->get_NX());
0316       unsigned int pixel_z(cell->get_pixel_index() / geom->get_NX());
0317       unsigned int halflayer = (int) pixel_x >= geom->get_NX() / 2 ? 0 : 1;
0318
0319       assert((int) pixel_x < geom->get_NX());
0320       assert((int) pixel_z < geom->get_NZ());
0321
0322       unsigned int halfLayerIndex(layer * 2 + halflayer);
0323 //      unsigned int pixelPhiIndex(
0324 //          cell->get_stave_index() * geom->get_NX() + pixel_x);
0325 //      unsigned int pixelPhiIndexHL(
0326 //          cell->get_stave_index() * geom->get_NX() / 2 + pixel_x % (geom->get_NX() / 2));
0327 //      unsigned int pixelZIndex(cell->get_chip_index() * geom->get_NZ() + pixel_z);
0328
0329       //      //      ptree hitTree;
0330       //      rapidjson::Value hitTree(rapidjson::kObjectType);
0331       //
0332       //      //      ptree hitIDTree;
0333       //      rapidjson::Value hitIDTree(rapidjson::kObjectType);
0334       //      hitIDTree.AddMember("HitSequenceInEvent", hitID, alloc);
0335       //
0336       //      hitIDTree.AddMember("PixelHalfLayerIndex", halfLayerIndex, alloc);
0337       //      hitIDTree.AddMember("PixelPhiIndexInLayer", pixelPhiIndex, alloc);
0338       //      hitIDTree.AddMember("PixelPhiIndexInHalfLayer", pixelPhiIndexHL, alloc);
0339       //      hitIDTree.AddMember("PixelZIndex", pixelZIndex, alloc);
0340       //
0341       //      hitIDTree.AddMember("Layer", layer, alloc);
0342       //      hitIDTree.AddMember("HalfLayer", halflayer, alloc);
0343       //      hitIDTree.AddMember("Stave", cell->get_stave_index(), alloc);
0344       //      //      hitIDTree.put("HalfStave", cell->get_half_stave_index());
0345       //      //      hitIDTree.put("Module", cell->get_module_index());
0346       //      hitIDTree.AddMember("Chip", cell->get_chip_index(), alloc);
0347       //      hitIDTree.AddMember("Pixel", cell->get_pixel_index(), alloc);
0348       //      hitTree.AddMember("ID", hitIDTree, alloc);
0349       //
0350       //      hitTree.AddMember("Coordinate",
0351       //                        loadCoordinate(world_coords.x(),
0352       //                                       world_coords.y(),
0353       //                                       world_coords.z()),
0354       //                        alloc);
0355
0356       if (Verbosity() >= 2)
0357         cout << "HFMLTriggerOccupancy::process_event - MVTX hit "
0358              << "layer " << layer << "\t"
0359              << "Stave " << cell->get_stave_index() << "\t"
0360              << "Chip " << cell->get_chip_index() << "\t"
0361              << "Pixel " << cell->get_pixel_index() << "\t"
0362              << "Coordinate"
0363              << "(" << world_coords.x() << "," << world_coords.y() << "," << world_coords.z() << ")" << endl;
0364
0365       //      rawHitsTree.add_child("MVTXHit", hitTree);
0366       //      rawHitsTree.PushBack(hitTree, alloc);
0367
0368       m_hitStaveLayer->Fill(cell->get_stave_index(), halfLayerIndex);
0369       m_hitModuleHalfStave->Fill(cell->get_module_index(), cell->get_half_stave_index());
0370       m_hitChipModule->Fill(cell->get_chip_index(), cell->get_module_index());
0371
0372       m_hitLayerMap->Fill(world_coords.x(), world_coords.y(), halfLayerIndex);
0373 //      m_hitPixelPhiMap->Fill(pixelPhiIndex, atan2(world_coords.y(), world_coords.x()), halfLayerIndex);
0374 //      m_hitPixelPhiMapHL->Fill(pixelPhiIndexHL, atan2(world_coords.y(), world_coords.x()), halfLayerIndex);
0375 //      m_hitPixelZMap->Fill(pixelZIndex, world_coords.z(), halfLayerIndex);
0376
0377       assert(layer < multiplicity_vec.size());
0378       auto& multiplicity_layer = multiplicity_vec[layer];
0379       assert(static_cast<size_t>(cell->get_stave_index()) < multiplicity_layer.size());
0380       auto& multiplicity_stave = multiplicity_layer[cell->get_stave_index()];
0381       assert(static_cast<size_t>(cell->get_chip_index()) < multiplicity_stave.size());
0382       multiplicity_stave[cell->get_chip_index()]++;
0383
0384     }  //    if (layer < _nlayers_maps)
0385
0386   }  //   for (SvtxHitMap::Iter iter = m_hitMap->begin();
0387      //  rawHitTree.AddMember("MVTXHits", rawHitsTree, alloc);
0388
0389   //  m_Multiplicity = new TH1F("Multiplicity", "Multiplicity", 10000, -.5, 10000 - .5);
0390   //  m_Multiplicity->SetTitle("Multiplicity;Chip Multiplicity;Count");
0391   //
0392   //  m_LayerMultiplicity = new TH2F("LayerMultiplicity", "LayerMultiplicity", 3, -.5, 2.5, 1000, -.5, 1000 - .5);
0393   //  m_LayerMultiplicity->SetTitle("LayerMultiplicity;Layer;Chip Multiplicity");
0394   //
0395   //  m_LayerChipMultiplicity = new TH3F("LayerChipMultiplicity", "LayerChipMultiplicity", 3, -.5, 2.5, 9, -.5, 8.5, 1000, -.5, 1000 - .5);
0396   //  m_LayerChipMultiplicity->SetTitle("LayerChipMultiplicity;Layer;Chip");
0397
0398   int layer = 0;
0399   for (auto& multiplicity_layer : multiplicity_vec)
0400   {
0401     int stave = 0;
0402     for (auto& multiplicity_stave : multiplicity_layer)
0403     {
0404       int chip = 0;
0405       for (auto& multiplicity_chip : multiplicity_stave)
0406       {
0407         m_Multiplicity->Fill(multiplicity_chip);
0408         m_LayerMultiplicity->Fill(layer, multiplicity_chip);
0409         m_LayerChipMultiplicity->Fill(layer, chip, multiplicity_chip);
0410
0411         chip++;
0412       }
0413
0414       if (Verbosity() >= 2)
0415         cout << "HFMLTriggerOccupancy::process_event - fill layer " << layer << " stave. " << stave << "Accumulated chips = "
0416              << m_Multiplicity->GetSum()
0417              << endl;
0418
0419       stave++;
0420     }
0421     layer++;
0422   }
0423
0424   ++_ievent;
0425
0426   return Fun4AllReturnCodes::EVENT_OK;
0427 }
0428
0429 int HFMLTriggerOccupancy::End(PHCompositeNode* topNode)
0430 {
0431   if (_f)
0432   {
0433     _f->cd();
0434     _f->Write();
0435     //_f->Close();
0436
0437     //    m_hitLayerMap->Write();
0438   }
0439
0440   //  if (m_jsonOut.is_open())
0441   //  {
0442   //    m_jsonOut << "]" << endl;
0443   //    m_jsonOut << "}" << endl;
0444   //
0445   //    m_jsonOut.close();
0446   //  }
0447
0448   cout << "HFMLTriggerOccupancy::End - output to " << _foutname << ".*" << endl;
0449
0450   return Fun4AllReturnCodes::EVENT_OK;
0451 }