hbcup/src/UPCTrigStudy.cc

0001 #include "UPCTrigStudy.h"
0002
0003 /// Cluster/Calorimeter includes
0004 /*
0005 #include <calobase/RawCluster.h>
0006 #include <calobase/RawClusterContainer.h>
0007 #include <calobase/RawClusterUtility.h>
0008 #include <calobase/RawTower.h>
0009 #include <calobase/RawTowerContainer.h>
0010 #include <calobase/RawTowerGeom.h>
0011 #include <calobase/RawTowerGeomContainer.h>
0012 #include <calotrigger/CaloTriggerInfo.h>
0013 */
0014
0015 #include <calobase/TowerInfoContainer.h>
0016 #include <calobase/TowerInfo.h>
0017
0018 /// Tracking includes
0019 #include <globalvertex/GlobalVertex.h>
0020 #include <globalvertex/GlobalVertexMap.h>
0021 #include <globalvertex/SvtxVertex.h>
0022 #include <globalvertex/SvtxVertexMap.h>
0023
0024 ///
0025 #include <mbd/MbdOut.h>
0026 #include <ffarawobjects/Gl1Packet.h>
0027
0028 /// Fun4All includes
0029 #include <fun4all/Fun4AllHistoManager.h>
0030 #include <fun4all/Fun4AllReturnCodes.h>
0031 #include <g4main/PHG4Hit.h>
0032 #include <g4main/PHG4Particle.h>
0033 #include <g4main/PHG4TruthInfoContainer.h>
0034 #include <phool/PHCompositeNode.h>
0035 #include <phool/getClass.h>
0036 #include <ffaobjects/EventHeader.h>
0037
0038 /// ROOT includes
0039 #include <TFile.h>
0040 #include <TH1.h>
0041 #include <TH2.h>
0042 #include <TNtuple.h>
0043 #include <TTree.h>
0044 #include <Math/Vector4D.h>
0045 #include <Math/VectorUtil.h>
0046
0047
0048 /// C++ includes
0049 #include <cassert>
0050 #include <cmath>
0051 #include <sstream>
0052 #include <string>
0053
0054 /**
0055  * Constructor of module
0056  */
0057 UPCTrigStudy::UPCTrigStudy(const std::string &name, const std::string &filename)
0058   : SubsysReco(name)
0059   , m_outfilename(filename)
0060   , m_hm(nullptr)
0061   //, m_mincluspt(0.25)
0062   //, m_analyzeTracks(true)
0063   //, m_analyzeClusters(true)
0064   //, m_analyzeTruth(false)
0065 {
0066   /// Initialize variables and trees so we don't accidentally access
0067   /// memory that was never allocated
0068   resetVariables();
0069   initializeTrees();
0070 }
0071
0072 /**
0073  * Destructor of module
0074  */
0075 UPCTrigStudy::~UPCTrigStudy()
0076 {
0077   delete m_hm;
0078   delete _globaltree;
0079 }
0080
0081 /**
0082  * Initialize the module and prepare looping over events
0083  */
0084 int UPCTrigStudy::Init(PHCompositeNode * /*topNode*/)
0085 {
0086   if (Verbosity() > 5)
0087   {
0088     std::cout << "Beginning Init in UPCTrigStudy" << std::endl;
0089   }
0090
0091   m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
0092
0093   h_trig = new TH1F("h_trig", "trig", 16, 0, 16);
0094   h_ntracks = new TH1F("h_ntracks", "num tracks", 2000, 0, 2000);
0095
0096   h_zdce = new TH1F("h_zdce","ZDC Energy",820,-50,4050);
0097   h_zdcse = new TH1F("h_zdcse","ZDC.S Energy",500,0,250);
0098   h_zdcne = new TH1F("h_zdcne","ZDC.N Energy",500,0,250);
0099   h_zdctimecut = new TH1F("h_zdctimecut", "zdctimecut", 50, -17.5 , 32.5);
0100
0101   return 0;
0102 }
0103
0104 /**
0105  * Main workhorse function where each event is looped over and
0106  * data from each event is collected from the node tree for analysis
0107  */
0108 int UPCTrigStudy::GetNodes(PHCompositeNode *topNode)
0109 {
0110   /// EventHeader node
0111   evthdr = findNode::getClass<EventHeader>(topNode, "EventHeader");
0112
0113   // Get the raw gl1 data from event combined DST
0114   _gl1raw = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0115   if(!_gl1raw && f_evt<4) std::cout << PHWHERE << " GL1Packet node not found on node tree" << std::endl;
0116
0117   // MbdOut
0118   _mbdout = findNode::getClass<MbdOut>(topNode, "MbdOut");
0119   if(!_mbdout && f_evt<4) std::cout << PHWHERE << " MbdOut node not found on node tree" << std::endl;
0120
0121   _zdctowers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_ZDC");
0122   if(!_zdctowers && f_evt<4) std::cout << PHWHERE << " TOWERINFO_CALIB_ZDC node not found on node tree" << std::endl;
0123
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126
0127 /**
0128  * Main workhorse function where each event is looped over and
0129  * data from each event is collected from the node tree for analysis
0130  */
0131 int UPCTrigStudy::process_event(PHCompositeNode *topNode)
0132 {
0133   if (Verbosity() > 5)
0134   {
0135     std::cout << "Beginning process_event in UPCTrigStudy" << std::endl;
0136   }
0137
0138   resetVariables();
0139
0140   h_trig->Fill( 0 );  // event counter
0141   _nprocessed++;
0142
0143   /// Get all the data nodes
0144   int status = GetNodes(topNode);
0145   if ( status != Fun4AllReturnCodes::EVENT_OK )
0146   {
0147     return status;
0148   }
0149
0150   /// Get the run and eventnumber
0151   if ( evthdr )
0152   {
0153     f_run = evthdr->get_RunNumber();
0154     f_evt = evthdr->get_EvtSequence();
0155   }
0156
0157   // process triggers, skip if not a clock trigger
0158   status = process_triggers();
0159   if ( status != Fun4AllReturnCodes::EVENT_OK)
0160   {
0161     return status;
0162   }
0163
0164   /// Get MBD information
0165   f_mbdsn = _mbdout->get_npmt(0);
0166   f_mbdnn = _mbdout->get_npmt(1);
0167
0168   process_zdc();
0169
0170   /// Get calorimeter information
0171   /*
0172   if (m_analyzeClusters)
0173   {
0174     getEMCalClusters(topNode);
0175   }
0176   */
0177
0178   _globaltree->Fill();
0179
0180   return Fun4AllReturnCodes::EVENT_OK;
0181 }
0182
0183 /**
0184  * End the module and finish any data collection. Clean up any remaining
0185  * loose ends
0186  */
0187 int UPCTrigStudy::End(PHCompositeNode * /*topNode*/)
0188 {
0189   if (Verbosity() > 5)
0190   {
0191     std::cout << "Ending UPCTrigStudy analysis package" << std::endl;
0192   }
0193
0194   /// Change to the outfile
0195   m_outfile->cd();
0196
0197   _globaltree->Write();
0198
0199   /// If we analyzed the clusters, write them out
0200   /*
0201   if (m_analyzeClusters)
0202   {
0203     m_clustertree->Write();
0204   }
0205   */
0206
0207   /// Write and close the outfile
0208   m_outfile->Write();
0209   m_outfile->Close();
0210
0211   /// Clean up pointers and associated histos/trees in TFile
0212   delete m_outfile;
0213
0214   if (Verbosity() > 1)
0215   {
0216     std::cout << "Finished UPCTrigStudy analysis package" << std::endl;
0217   }
0218
0219   return 0;
0220 }
0221
0222 int UPCTrigStudy::process_triggers()
0223 {
0224   // Only use MBDNS triggered events
0225   if ( _gl1raw != nullptr )
0226   {
0227     const uint64_t CLOCK = 0x1;        // CLOCK trigger bit
0228     const uint64_t MBDTRIGS = 0x7c00;  // MBDNS trigger bits
0229     const uint64_t MBDWIDE = 0xc00;    // MBDNS wide bits (no vtx cut)
0230     //const uint64_t ZDCNS = 0x8;        // ZDCNS trigger bits
0231
0232     f_cross = _gl1raw->getBunchNumber();
0233     f_rtrig = _gl1raw->getTriggerVector();
0234     f_ltrig = _gl1raw->getLiveVector();
0235     f_strig = _gl1raw->getScaledVector();
0236
0237     if ( (f_strig&CLOCK) == 0 )
0238     {
0239       return Fun4AllReturnCodes::ABORTEVENT;
0240     }
0241
0242     if ( f_strig&MBDWIDE )
0243     {
0244       _mbdwide = 1;
0245     }
0246     else
0247     {
0248       _mbdwide = 0;
0249     }
0250
0251     if ( _nprocessed<10 )
0252     {
0253       std::cout << "Trig " << _nprocessed << std::endl;
0254       std::cout << std::hex;
0255       std::cout << "Raw\t" << f_rtrig << std::endl;
0256       std::cout << "Live\t" << f_ltrig << std::endl;
0257       std::cout << "Scaled\t" << f_strig << std::endl;
0258       std::cout << "AND\t" << std::hex << (f_strig&MBDTRIGS) << std::dec << std::endl;
0259       std::cout << "WIDE\t" << _mbdwide << std::endl;
0260       std::cout << std::dec;
0261
0262     }
0263   }
0264
0265   return Fun4AllReturnCodes::EVENT_OK;
0266 }
0267
0268
0269 void UPCTrigStudy::process_zdc()
0270 {
0271   int max_zdc_t = Getpeaktime( h_zdctimecut );
0272   int range = 1;
0273   float zdc_etot = 0.;
0274
0275   int size = _zdctowers->size(); //online towers should be the same!
0276   for (int ich = 0; ich < size; ich++)
0277   {
0278     TowerInfo* zdctower = _zdctowers->get_tower_at_channel(ich);
0279     float zdce = zdctower->get_energy();
0280     int zdct = _zdctowers->get_tower_at_channel(ich)->get_time();
0281     h_zdctimecut->Fill( zdct );
0282
0283     if ( (zdct  < (max_zdc_t - range)) ||  (zdct  > (max_zdc_t + range)) )
0284     {
0285       continue;
0286     }
0287
0288     //
0289     if (ich==0||ich==2||ich==4)
0290     {
0291       f_zdcse += zdce;
0292     }
0293     else if (ich == 8 || ich == 10 || ich == 12)
0294     {
0295       f_zdcne += zdce;
0296     }
0297
0298   }
0299
0300   h_zdcse->Fill( f_zdcse );
0301   h_zdcne->Fill( f_zdcne );
0302   zdc_etot = f_zdcse + f_zdcne;
0303   h_zdce->Fill( zdc_etot );
0304 }
0305
0306 int UPCTrigStudy::Getpeaktime(TH1 * h)
0307 {
0308   int getmaxtime, tcut = -1;
0309
0310   for(int bin = 1; bin < h->GetNbinsX()+1; bin++)
0311   {
0312     double c = h->GetBinContent(bin);
0313     double max = h->GetMaximum();
0314     int bincenter = h->GetBinCenter(bin);
0315     if(max == c)
0316     {
0317       getmaxtime = bincenter;
0318       if(getmaxtime != -1) tcut = getmaxtime;
0319     }
0320   }
0321
0322   return tcut;
0323 }
0324
0325
0326
0327 /**
0328  * This method gets clusters from the EMCal and stores them in a tree. It
0329  * also demonstrates how to get trigger emulator information. Clusters from
0330  * other containers can be obtained in a similar way (e.g. clusters from
0331  * the IHCal, etc.)
0332  */
0333 /*
0334    void UPCTrigStudy::getEMCalClusters(PHCompositeNode *topNode)
0335    {
0336 /// Get the raw cluster container
0337 /// Note: other cluster containers exist as well. Check out the node tree when
0338 /// you run a simulation
0339 RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0340
0341 if (!clusters)
0342 {
0343 std::cout << PHWHERE
0344 << "EMCal cluster node is missing, can't collect EMCal clusters"
0345 << std::endl;
0346 return;
0347 }
0348
0349 /// Get the global vertex to determine the appropriate pseudorapidity of the clusters
0350 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0351 if (!vertexmap)
0352 {
0353 std::cout << "UPCTrigStudy::getEmcalClusters - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0354 assert(vertexmap);  // force quit
0355
0356 return;
0357 }
0358
0359 if (vertexmap->empty())
0360 {
0361 std::cout << "UPCTrigStudy::getEmcalClusters - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0362 return;
0363 }
0364
0365 /// just take a bbc vertex
0366 GlobalVertex *vtx = nullptr;
0367 for(auto iter = vertexmap->begin(); iter!= vertexmap->end(); ++iter)
0368 {
0369 GlobalVertex* vertex = iter->second;
0370 if(vertex->find_vtxids(GlobalVertex::MBD) != vertex->end_vtxids())
0371 {
0372 vtx = vertex;
0373 }
0374 }
0375 if (vtx == nullptr)
0376 {
0377 return;
0378 }
0379
0380 /// Trigger emulator
0381 CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topNode, "CaloTriggerInfo");
0382
0383 /// Can obtain some trigger information if desired
0384 if (trigger)
0385 {
0386 m_E_4x4 = trigger->get_best_EMCal_4x4_E();
0387 }
0388 RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0389 RawClusterContainer::ConstIterator clusIter;
0390
0391 /// Loop over the EMCal clusters
0392 for (clusIter = begin_end.first; clusIter != begin_end.second; ++clusIter)
0393 {
0394 /// Get this cluster
0395 const RawCluster *cluster = clusIter->second;
0396
0397 /// Get cluster characteristics
0398 /// This helper class determines the photon characteristics
0399 /// depending on the vertex position
0400 /// This is important for e.g. eta determination and E_T determination
0401 CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0402 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
0403 m_clusenergy = E_vec_cluster.mag();
0404 m_cluseta = E_vec_cluster.pseudoRapidity();
0405 m_clustheta = E_vec_cluster.getTheta();
0406 m_cluspt = E_vec_cluster.perp();
0407 m_clusphi = E_vec_cluster.getPhi();
0408
0409 if (m_cluspt < m_mincluspt)
0410 {
0411   continue;
0412 }
0413
0414 m_cluspx = m_cluspt * cos(m_clusphi);
0415 m_cluspy = m_cluspt * sin(m_clusphi);
0416 m_cluspz = sqrt(m_clusenergy * m_clusenergy - m_cluspx * m_cluspx - m_cluspy * m_cluspy);
0417
0418 // fill the cluster tree with all emcal clusters
0419 m_clustertree->Fill();
0420 }
0421 }
0422 */
0423
0424 /**
0425  * This function puts all of the tree branch assignments in one place so as to not
0426  * clutter up the UPCTrigStudy::Init function.
0427  */
0428 void UPCTrigStudy::initializeTrees()
0429 {
0430   _globaltree = new TTree("gt", "Global Info");
0431   _globaltree->Branch("run", &f_run, "run/I");
0432   _globaltree->Branch("evt", &f_evt, "evt/I");
0433   _globaltree->Branch("rtrig",&f_rtrig,"rtrig/l");
0434   _globaltree->Branch("ltrig",&f_ltrig,"ltrig/l");
0435   _globaltree->Branch("strig",&f_strig,"strig/l");
0436   _globaltree->Branch("zdcse",&f_zdcse,"zdcse/F");
0437   _globaltree->Branch("zdcne",&f_zdcne,"zdcne/F");
0438   _globaltree->Branch("cross",&f_cross,"cross/S");
0439   _globaltree->Branch("mbdsn",&f_mbdsn,"mbdsn/S");
0440   _globaltree->Branch("mbdnn",&f_mbdnn,"mbdnn/S");
0441
0442   //_globaltree->Branch("ntrks", &m_ntracks, "ntrks/I");
0443
0444   /*
0445      m_clustertree = new TTree("clustertree", "A tree with emcal clusters");
0446      m_clustertree->Branch("m_clusenergy", &m_clusenergy, "m_clusenergy/D");
0447      m_clustertree->Branch("m_cluseta", &m_cluseta, "m_cluseta/D");
0448      m_clustertree->Branch("m_clustheta", &m_clustheta, "m_clustheta/D");
0449      m_clustertree->Branch("m_cluspt", &m_cluspt, "m_cluspt/D");
0450      m_clustertree->Branch("m_clusphi", &m_clusphi, "m_clusphi/D");
0451      m_clustertree->Branch("m_cluspx", &m_cluspx, "m_cluspx/D");
0452      m_clustertree->Branch("m_cluspy", &m_cluspy, "m_cluspy/D");
0453      m_clustertree->Branch("m_cluspz", &m_cluspz, "m_cluspz/D");
0454      m_clustertree->Branch("m_E_4x4", &m_E_4x4, "m_E_4x4/D");
0455      */
0456 }
0457
0458 /**
0459  * This function initializes all of the member variables in this class so that there
0460  * are no variables that might not be set before e.g. writing them to the output
0461  * trees.
0462  */
0463 void UPCTrigStudy::resetVariables()
0464 {
0465   f_evt = 0;
0466   f_rtrig = 0UL;
0467   f_ltrig = 0UL;
0468   f_strig = 0UL;
0469   f_zdcse = 0.;
0470   f_zdcne = 0.;
0471   f_mbdsn = 0;
0472   f_mbdnn = 0;
0473   f_ntracks = 0;
0474   _mbdwide = 0;
0475 }
0476