sepd/sepd_eventplanecalib/sEPD_TreeGen.cc

0001 #include "sEPD_TreeGen.h"
0002 #include "QVecDefs.h"
0003 #include "EventPlaneData.h"
0004
0005 // -- Calo
0006 #include <calobase/TowerInfo.h>
0007 #include <calobase/TowerInfoContainer.h>
0008
0009 // -- Vtx
0010 #include <globalvertex/GlobalVertex.h>
0011 #include <globalvertex/GlobalVertexMap.h>
0012
0013 // -- MB
0014 #include <calotrigger/MinimumBiasInfo.h>
0015
0016 #include <centrality/CentralityInfo.h>
0017
0018 // -- sEPD
0019 #include <epd/EpdGeom.h>
0020
0021 // -- event
0022 #include <ffaobjects/EventHeader.h>
0023
0024 // -- Fun4All
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026 #include <fun4all/Fun4AllServer.h>
0027
0028 // -- Nodes
0029 #include <phool/PHCompositeNode.h>
0030 #include <phool/PHNodeIterator.h>
0031 #include <phool/getClass.h>
0032
0033 // -- ROOT
0034 #include <TProfile.h>
0035 #include <TH2.h>
0036
0037 // -- c++
0038 #include <iomanip>
0039 #include <iostream>
0040
0041 //____________________________________________________________________________..
0042 sEPD_TreeGen::sEPD_TreeGen(const std::string &name)
0043   : SubsysReco(name)
0044 {
0045 }
0046
0047 //____________________________________________________________________________..
0048 int sEPD_TreeGen::Init(PHCompositeNode *topNode)
0049 {
0050   Fun4AllServer *se = Fun4AllServer::instance();
0051   if (Verbosity() > 0)
0052   {
0053     se->Print("NODETREE");
0054   }
0055   unsigned int bins_sepd_totalcharge{100};
0056   double sepd_totalcharge_low{0};
0057   double sepd_totalcharge_high{2e4};
0058
0059   unsigned int bins_centrality{80};
0060   double centrality_low{-0.5};
0061   double centrality_high{79.5};
0062
0063   hSEPD_Charge = new TProfile("hSEPD_Charge", "|z| < 10 cm and MB; Channel; Avg Charge", QVecShared::SEPD_CHANNELS, 0, QVecShared::SEPD_CHANNELS);
0064   hSEPD_Charge->Sumw2();
0065
0066   h2SEPD_totalcharge_centrality = new TH2F("h2SEPD_totalcharge_centrality",
0067                                            "|z| < 10 cm and MB; sEPD Total Charge; Centrality [%]",
0068                                            bins_sepd_totalcharge, sepd_totalcharge_low, sepd_totalcharge_high,
0069                                            bins_centrality, centrality_low, centrality_high);
0070
0071   se->registerHisto(hSEPD_Charge);
0072   se->registerHisto(h2SEPD_totalcharge_centrality);
0073
0074   PHNodeIterator node_itr(topNode);
0075   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(node_itr.findFirst("PHCompositeNode", "DST"));
0076
0077   if (!dstNode)
0078   {
0079     std::cout << PHWHERE << "DST node missing, cannot attach EventPlaneData." << std::endl;
0080     return Fun4AllReturnCodes::ABORTRUN;
0081   }
0082
0083   EventPlaneData *evtdata = findNode::getClass<EventPlaneData>(topNode, "EventPlaneData");
0084   if (!evtdata)
0085   {
0086     evtdata = new EventPlaneData();
0087     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(evtdata, "EventPlaneData", "PHObject");
0088     dstNode->addNode(newNode);
0089   }
0090
0091   return Fun4AllReturnCodes::EVENT_OK;
0092 }
0093
0094 //____________________________________________________________________________..
0095 int sEPD_TreeGen::process_event_check(PHCompositeNode *topNode)
0096 {
0097   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0098
0099   if (!vertexmap)
0100   {
0101     std::cout << PHWHERE << "GlobalVertexMap Node missing, doing nothing." << std::endl;
0102     return Fun4AllReturnCodes::ABORTRUN;
0103   }
0104
0105   if (vertexmap->empty())
0106   {
0107     if (Verbosity() > 1)
0108     {
0109       std::cout << PHWHERE << "GlobalVertexMap Empty, Skipping Event: " << m_data.event_id << std::endl;
0110     }
0111     return Fun4AllReturnCodes::ABORTEVENT;
0112   }
0113
0114   GlobalVertex *vtx = vertexmap->begin()->second;
0115   double zvtx = vtx->get_z();
0116
0117   MinimumBiasInfo *m_mb_info = findNode::getClass<MinimumBiasInfo>(topNode, "MinimumBiasInfo");
0118   if (!m_mb_info)
0119   {
0120     std::cout << PHWHERE << "MinimumBiasInfo Node missing, doing nothing." << std::endl;
0121     return Fun4AllReturnCodes::ABORTRUN;
0122   }
0123
0124   // skip event if not minimum bias
0125   if (!m_mb_info->isAuAuMinimumBias())
0126   {
0127     if (Verbosity() > 1)
0128     {
0129       std::cout << "Event: " << m_data.event_id << ", Not Min Bias, Skipping" << std::endl;
0130     }
0131     return Fun4AllReturnCodes::ABORTEVENT;
0132   }
0133
0134   // skip event if zvtx is too large
0135   if (std::abs(zvtx) >= m_cuts.m_zvtx_max)
0136   {
0137     if (Verbosity() > 1)
0138     {
0139       std::cout << "Event: " << m_data.event_id << ", Z: " << zvtx << " cm, Skipping" << std::endl;
0140     }
0141     return Fun4AllReturnCodes::ABORTEVENT;
0142   }
0143
0144   m_evtdata->set_event_zvertex(zvtx);
0145
0146   return Fun4AllReturnCodes::EVENT_OK;
0147 }
0148
0149 //____________________________________________________________________________..
0150 int sEPD_TreeGen::process_centrality(PHCompositeNode *topNode)
0151 {
0152   CentralityInfo *centInfo = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0153   if (!centInfo)
0154   {
0155     std::cout << PHWHERE << "CentralityInfo Node missing, doing nothing." << std::endl;
0156     return Fun4AllReturnCodes::ABORTRUN;
0157   }
0158
0159   double cent = centInfo->get_centile(CentralityInfo::PROP::mbd_NS) * 100;
0160
0161   // skip event if centrality is bad or too peripheral
0162   if (!std::isfinite(cent) || cent < 0 || cent >= m_cuts.m_cent_max)
0163   {
0164     if (Verbosity() > 1)
0165     {
0166       std::cout << "Event: " << m_data.event_id << ", Centrality: " << cent << ", Skipping" << std::endl;
0167     }
0168     return Fun4AllReturnCodes::ABORTEVENT;
0169   }
0170
0171   m_evtdata->set_event_centrality(cent);
0172   m_data.event_centrality = cent;
0173
0174   return Fun4AllReturnCodes::EVENT_OK;
0175 }
0176
0177 //____________________________________________________________________________..
0178 int sEPD_TreeGen::process_sEPD(PHCompositeNode *topNode)
0179 {
0180   TowerInfoContainer *towerinfosEPD = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_SEPD");
0181   if (!towerinfosEPD)
0182   {
0183     std::cout << PHWHERE << "TOWERINFO_CALIB_SEPD Node missing, doing nothing." << std::endl;
0184     return Fun4AllReturnCodes::ABORTRUN;
0185   }
0186
0187   EpdGeom *epdgeom = findNode::getClass<EpdGeom>(topNode, "TOWERGEOM_EPD");
0188   if (!epdgeom)
0189   {
0190     std::cout << PHWHERE << "TOWERGEOM_EPD Node missing, doing nothing." << std::endl;
0191     return Fun4AllReturnCodes::ABORTRUN;
0192   }
0193
0194   // sepd
0195   unsigned int sepd_channels = towerinfosEPD->size();
0196
0197   if(sepd_channels != QVecShared::SEPD_CHANNELS)
0198   {
0199     if (Verbosity() > 1)
0200     {
0201       std::cout << "Event: " << m_data.event_id << ", SEPD Channels = " << sepd_channels << " != " << QVecShared::SEPD_CHANNELS << std::endl;
0202     }
0203     return Fun4AllReturnCodes::ABORTEVENT;
0204   }
0205
0206   double sepd_totalcharge = 0;
0207
0208   for (unsigned int channel = 0; channel < sepd_channels; ++channel)
0209   {
0210     TowerInfo *tower = towerinfosEPD->get_tower_at_channel(channel);
0211
0212     if (!tower)
0213     {
0214       if (Verbosity() > 1)
0215       {
0216         std::cout << PHWHERE << "Null SEPD tower at channel " << channel << std::endl;
0217       }
0218       continue;
0219     }
0220
0221     double charge = tower->get_energy();
0222     bool isZS = tower->get_isZS();
0223
0224     // exclude ZS
0225     // exclude Nmips
0226     if (isZS || charge < m_cuts.m_sepd_charge_min)
0227     {
0228       continue;
0229     }
0230
0231     m_evtdata->set_sepd_charge(channel, charge);
0232
0233     sepd_totalcharge += charge;
0234
0235     hSEPD_Charge->Fill(channel, charge);
0236   }
0237
0238   m_evtdata->set_sepd_totalcharge(sepd_totalcharge);
0239   h2SEPD_totalcharge_centrality->Fill(sepd_totalcharge, m_data.event_centrality);
0240
0241   return Fun4AllReturnCodes::EVENT_OK;
0242 }
0243
0244 //____________________________________________________________________________..
0245 void sEPD_TreeGen::Print([[maybe_unused]] const std::string &what) const
0246 {
0247   // Only execute if Verbosity is high enough
0248   if (Verbosity() <= 2)
0249   {
0250     return;
0251   }
0252
0253   std::cout << "\n============================================================" << std::endl;
0254   std::cout << "sEPD_TreeGen::Print -> Event Data State" << std::endl;
0255
0256   if (!m_evtdata)
0257   {
0258     std::cout << " [WARNING] m_evtdata is null." << std::endl;
0259     return;
0260   }
0261
0262   // Verbosity > 2: Print basic scalars
0263   std::cout << "  Event ID:          " << m_evtdata->get_event_id() << std::endl;
0264   std::cout << "  Z-Vertex:          " << m_evtdata->get_event_zvertex() << " cm" << std::endl;
0265   std::cout << "  Centrality:        " << m_evtdata->get_event_centrality() << " %" << std::endl;
0266   std::cout << "  sEPD Total Charge: " << m_evtdata->get_sepd_totalcharge() << std::endl;
0267
0268   // Verbosity > 3: Print channel arrays
0269   if (Verbosity() > 3)
0270   {
0271     std::cout << "  Active Towers (Charge > 0):" << std::endl;
0272     for (int i = 0; i < QVecShared::SEPD_CHANNELS; ++i)
0273     {
0274       double charge = m_evtdata->get_sepd_charge(i);
0275       if (charge > 0)
0276       {
0277         std::cout << "    Channel: " << std::setw(3) << i
0278                   << " | Charge: " << std::fixed << std::setprecision(4) << charge << std::endl;
0279       }
0280     }
0281   }
0282   std::cout << "============================================================\n" << std::endl;
0283 }
0284
0285 //____________________________________________________________________________..
0286 int sEPD_TreeGen::process_event(PHCompositeNode *topNode)
0287 {
0288   EventHeader *eventInfo = findNode::getClass<EventHeader>(topNode, "EventHeader");
0289   if (!eventInfo)
0290   {
0291     std::cout << PHWHERE << "EventHeader Node missing, doing nothing." << std::endl;
0292     return Fun4AllReturnCodes::ABORTRUN;
0293   }
0294
0295   m_data.event_id = eventInfo->get_EvtSequence();
0296
0297   if (Verbosity() && m_event % PROGRESS_PRINT_INTERVAL == 0)
0298   {
0299     std::cout << "Progress: " << m_event << ", Global: " << m_data.event_id << std::endl;
0300   }
0301   ++m_event;
0302
0303   m_evtdata = findNode::getClass<EventPlaneData>(topNode, "EventPlaneData");
0304   if (!m_evtdata)
0305   {
0306     std::cout << PHWHERE << "EventPlaneData Node missing." << std::endl;
0307     return Fun4AllReturnCodes::ABORTRUN;
0308   }
0309
0310   m_evtdata->set_event_id(m_data.event_id);
0311
0312   int ret = process_event_check(topNode);
0313   if (ret)
0314   {
0315     return ret;
0316   }
0317
0318   ret = process_centrality(topNode);
0319   if (ret)
0320   {
0321     return ret;
0322   }
0323
0324   ret = process_sEPD(topNode);
0325   if (ret)
0326   {
0327     return ret;
0328   }
0329
0330   Print();
0331
0332   return Fun4AllReturnCodes::EVENT_OK;
0333 }
0334
0335 //____________________________________________________________________________..
0336 int sEPD_TreeGen::ResetEvent([[maybe_unused]] PHCompositeNode *topNode)
0337 {
0338   // Event
0339   m_data.event_id = -1;
0340   m_data.event_centrality = 9999;
0341
0342   return Fun4AllReturnCodes::EVENT_OK;
0343 }