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