packages/eventplaneinfo/EventPlaneReco.cc

0001 #include "EventPlaneReco.h"
0002
0003 #include "Eventplaneinfo.h"
0004 #include "EventplaneinfoMap.h"
0005 #include "EventplaneinfoMapv1.h"
0006 #include "Eventplaneinfov1.h"
0007
0008 #include <calobase/TowerInfo.h>
0009 #include <calobase/TowerInfoContainer.h>
0010 #include <calobase/TowerInfoDefs.h>
0011
0012 #include <epd/EpdGeom.h>
0013
0014 #include <mbd/MbdGeom.h>
0015 #include <mbd/MbdOut.h>
0016 #include <mbd/MbdPmtContainer.h>
0017 #include <mbd/MbdPmtHit.h>
0018
0019 #include <globalvertex/GlobalVertex.h>
0020 #include <globalvertex/GlobalVertexMap.h>
0021 #include <globalvertex/MbdVertex.h>
0022 #include <globalvertex/MbdVertexMap.h>
0023
0024 #include <cdbobjects/CDBHistos.h>
0025 #include <cdbobjects/CDBTTree.h>
0026 #include <ffamodules/CDBInterface.h>
0027
0028 #include <fun4all/Fun4AllReturnCodes.h>
0029 #include <fun4all/SubsysReco.h> // for SubsysReco
0030
0031 #include <phool/PHCompositeNode.h>
0032 #include <phool/PHIODataNode.h>
0033 #include <phool/PHNode.h> // for PHNode
0034 #include <phool/PHNodeIterator.h>
0035 #include <phool/PHObject.h> // for PHObject
0036 #include <phool/getClass.h>
0037 #include <phool/phool.h> // for PHWHERE
0038 #include <phool/recoConsts.h>
0039
0040 #include <TProfile2D.h>
0041
0042 #include <array> // for array
0043 #include <cfloat>
0044 #include <cmath>
0045 #include <cstdlib> // for exit
0046 #include <format>
0047 #include <iostream>
0048 #include <set>     // for _Rb_tree_const_iterator
0049 #include <utility> // for pair
0050 #include <vector>  // for vector
0051
0052 EventPlaneReco::EventPlaneReco(const std::string &name) : SubsysReco(name) {
0053
0054   south_q.resize(m_MaxOrder);
0055   north_q.resize(m_MaxOrder);
0056   northsouth_q.resize(m_MaxOrder);
0057
0058   south_q_subtract.resize(m_MaxOrder);
0059   north_q_subtract.resize(m_MaxOrder);
0060   northsouth_q_subtract.resize(m_MaxOrder);
0061
0062   shift_north.resize(m_MaxOrder);
0063   shift_south.resize(m_MaxOrder);
0064   shift_northsouth.resize(m_MaxOrder);
0065   tmp_south_psi.resize(m_MaxOrder);
0066   tmp_north_psi.resize(m_MaxOrder);
0067   tmp_northsouth_psi.resize(m_MaxOrder);
0068
0069   for (auto &vec : south_q) {
0070     vec.resize(2);
0071   }
0072
0073   for (auto &vec : north_q) {
0074     vec.resize(2);
0075   }
0076
0077   for (auto &vec : northsouth_q) {
0078     vec.resize(2);
0079   }
0080
0081   for (auto &vec : south_q_subtract) {
0082     vec.resize(2);
0083   }
0084
0085   for (auto &vec : north_q_subtract) {
0086     vec.resize(2);
0087   }
0088
0089   for (auto &vec : northsouth_q_subtract) {
0090     vec.resize(2);
0091   }
0092
0093   ring_q_north.resize(nRings);
0094   ring_q_south.resize(nRings);
0095
0096   for (auto &rq : ring_q_north) {
0097     rq.resize(m_MaxOrder, std::vector<double>(2, 0.0));
0098   }
0099   for (auto &rq : ring_q_south) {
0100     rq.resize(m_MaxOrder, std::vector<double>(2, 0.0));
0101   }
0102
0103   all_ring_Qvecs_north.assign(
0104       nRings, std::vector<std::pair<double, double>>(m_MaxOrder, {0.0, 0.0}));
0105
0106   all_ring_Qvecs_south.assign(
0107       nRings, std::vector<std::pair<double, double>>(m_MaxOrder, {0.0, 0.0}));
0108 }
0109
0110 int EventPlaneReco::InitRun(PHCompositeNode *topNode) {
0111
0112   FileName = "EVENTPLANE_CORRECTION";
0113   if (_isSim) {
0114     FileName = "EVENTPLANE_CORRECTION_SIM";
0115   }
0116
0117   std::string calibdir = CDBInterface::instance()->getUrl(FileName);
0118
0119   if (calibdir.empty()) {
0120     std::cout << PHWHERE << "No Eventplane calibration file for domain "
0121               << FileName << " found" << std::endl;
0122     std::cout << PHWHERE
0123               << "Will only produce raw Q vectors and event plane angles "
0124               << std::endl;
0125   }
0126
0127   CDBHistos *cdbhistosIn = new CDBHistos(calibdir);
0128   cdbhistosIn->LoadCalibrations();
0129
0130   // Get phiweights
0131   h_phi_weight_south_input =
0132       dynamic_cast<TH1F *>(cdbhistosIn->getHisto("h_phi_weight_south", false));
0133   h_phi_weight_north_input =
0134       dynamic_cast<TH1F *>(cdbhistosIn->getHisto("h_phi_weight_north", false));
0135
0136   // Get recentering histograms
0137   for (unsigned int order = 0; order < m_MaxOrder; order++) {
0138     tprof_mean_cos_south_epd_input[order] =
0139         dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0140             std::format("tprof_mean_cos_south_epd_order_{}", order), false));
0141     tprof_mean_sin_south_epd_input[order] =
0142         dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0143             std::format("tprof_mean_sin_south_epd_order_{}", order), false));
0144     tprof_mean_cos_north_epd_input[order] =
0145         dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0146             std::format("tprof_mean_cos_north_epd_order_{}", order), false));
0147     tprof_mean_sin_north_epd_input[order] =
0148         dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0149             std::format("tprof_mean_sin_north_epd_order_{}", order), false));
0150     tprof_mean_cos_northsouth_epd_input[order] =
0151         dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0152             std::format("tprof_mean_cos_northsouth_epd_order_{}", order),
0153             false));
0154     tprof_mean_sin_northsouth_epd_input[order] =
0155         dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0156             std::format("tprof_mean_sin_northsouth_epd_order_{}", order),
0157             false));
0158   }
0159
0160   // Get shifting histograms
0161   for (unsigned int order = 0; order < m_MaxOrder; order++) {
0162     for (int p = 0; p < _imax; p++) {
0163       tprof_cos_north_epd_shift_input[order][p] =
0164           dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0165               std::format("tprof_cos_north_epd_shift_order_{}_{}", order, p),
0166               false));
0167       tprof_sin_north_epd_shift_input[order][p] =
0168           dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0169               std::format("tprof_sin_north_epd_shift_order_{}_{}", order, p),
0170               false));
0171       tprof_cos_south_epd_shift_input[order][p] =
0172           dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0173               std::format("tprof_cos_south_epd_shift_order_{}_{}", order, p),
0174               false));
0175       tprof_sin_south_epd_shift_input[order][p] =
0176           dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0177               std::format("tprof_sin_south_epd_shift_order_{}_{}", order, p),
0178               false));
0179       tprof_cos_northsouth_epd_shift_input[order][p] =
0180           dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0181               std::format("tprof_cos_northsouth_epd_shift_order_{}_{}", order,
0182                           p),
0183               false));
0184       tprof_sin_northsouth_epd_shift_input[order][p] =
0185           dynamic_cast<TProfile2D *>(cdbhistosIn->getHisto(
0186               std::format("tprof_sin_northsouth_epd_shift_order_{}_{}", order,
0187                           p),
0188               false));
0189     }
0190   }
0191
0192   if (Verbosity() > 1) {
0193     cdbhistosIn->Print();
0194   }
0195
0196   return CreateNodes(topNode);
0197 }
0198
0199 int EventPlaneReco::process_event(PHCompositeNode *topNode) {
0200   if (Verbosity() > 1) {
0201     std::cout << "EventPlaneReco::process_event -- entered" << std::endl;
0202   }
0203
0204   //---------------------------------
0205   // Get Objects off of the Node Tree
0206   //---------------------------------
0207
0208   if (_isSim) {
0209     // Use GlobalVertexMap for simulation
0210     GlobalVertexMap *vertexmap =
0211         findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0212     if (!vertexmap) {
0213       std::cout << PHWHERE << "::ERROR - cannot find GlobalVertexMap"
0214                 << std::endl;
0215       exit(-1);
0216     }
0217
0218     if (!vertexmap->empty()) {
0219       GlobalVertex *vtx = vertexmap->begin()->second;
0220       if (vtx) {
0221         _mbdvtx = vtx->get_z();
0222       }
0223     }
0224   } else {
0225     // Use MbdVertexMap for data
0226     MbdVertexMap *mbdvtxmap =
0227         findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0228     if (!mbdvtxmap) {
0229       std::cout << PHWHERE << "::ERROR - cannot find MbdVertexMap" << std::endl;
0230       exit(-1);
0231     }
0232
0233     MbdVertex *mvertex = nullptr;
0234     if (mbdvtxmap) {
0235       for (MbdVertexMap::ConstIter mbditer = mbdvtxmap->begin();
0236            mbditer != mbdvtxmap->end(); ++mbditer) {
0237         mvertex = mbditer->second;
0238       }
0239       if (mvertex) {
0240         _mbdvtx = mvertex->get_z();
0241       }
0242     }
0243   }
0244
0245   EventplaneinfoMap *epmap =
0246       findNode::getClass<EventplaneinfoMap>(topNode, "EventplaneinfoMap");
0247   if (!epmap) {
0248     std::cout << PHWHERE << "::ERROR - cannot find EventplaneinfoMap"
0249               << std::endl;
0250     exit(-1);
0251   }
0252
0253   if (_sepdEpReco) {
0254
0255     TowerInfoContainer *epd_towerinfo =
0256         findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_SEPD");
0257     if (!epd_towerinfo) {
0258       epd_towerinfo = findNode::getClass<TowerInfoContainer>(
0259           topNode, "TOWERINFO_CALIB_EPD");
0260       if (!epd_towerinfo) {
0261         std::cout << PHWHERE
0262                   << "::ERROR - cannot find sEPD Calibrated TowerInfoContainer"
0263                   << std::endl;
0264         exit(-1);
0265       }
0266     }
0267
0268     EpdGeom *_epdgeom = findNode::getClass<EpdGeom>(topNode, "TOWERGEOM_EPD");
0269     if (!_epdgeom) {
0270       std::cout << PHWHERE << "::ERROR - cannot find TOWERGEOM_EPD"
0271                 << std::endl;
0272       exit(-1);
0273     }
0274
0275     ResetMe();
0276
0277     if ((std::fabs(_mbdvtx) < _mbd_vertex_cut)) {
0278
0279       unsigned int ntowers = epd_towerinfo->size();
0280       for (unsigned int ch = 0; ch < ntowers; ch++) {
0281         TowerInfo *_tower = epd_towerinfo->get_tower_at_channel(ch);
0282         float epd_e = _tower->get_energy();
0283         bool isZS = _tower->get_isZS();
0284         if (!isZS) // exclude ZS
0285         {
0286           unsigned int key = TowerInfoDefs::encode_epd(ch);
0287           int arm = TowerInfoDefs::get_epd_arm(key);
0288           if (arm == 0) {
0289             _ssum += epd_e;
0290           } else if (arm == 1) {
0291             _nsum += epd_e;
0292           }
0293         }
0294       }
0295
0296       if (_ssum > _epd_charge_min && _nsum > _epd_charge_min &&
0297           _ssum < _epd_charge_max && _nsum < _epd_charge_max) {
0298         _do_ep = true;
0299       }
0300
0301       if (_do_ep) {
0302
0303         // Apply phi weights in builiding ring Q-vectors
0304         for (unsigned int ch = 0; ch < ntowers; ch++) {
0305           TowerInfo *_tower = epd_towerinfo->get_tower_at_channel(ch);
0306           float epd_e = _tower->get_energy();
0307           bool isZS = _tower->get_isZS();
0308           if (!isZS) // exclude ZS
0309           {
0310             if (epd_e < 0.2) // expecting Nmips
0311             {
0312               continue;
0313             }
0314             unsigned int key = TowerInfoDefs::encode_epd(ch);
0315             float tile_phi = _epdgeom->get_phi(key);
0316             int arm = TowerInfoDefs::get_epd_arm(key);
0317             int rbin = TowerInfoDefs::get_epd_rbin(key);
0318             int phibin = TowerInfoDefs::get_epd_phibin(key);
0319
0320             float truncated_e =
0321                 (epd_e < _epd_e) ? epd_e : _epd_e; // set cutoff at _epd_e
0322
0323             float TileWeight = truncated_e; // default
0324
0325             if (h_phi_weight_south_input && h_phi_weight_north_input) {
0326               if (arm == 0) {
0327                 TileWeight =
0328                     truncated_e * h_phi_weight_south_input->GetBinContent(
0329                                       phibin + 1); // scale by 1/<N(φ_{i})>
0330               } else if (arm == 1) {
0331                 TileWeight =
0332                     truncated_e * h_phi_weight_north_input->GetBinContent(
0333                                       phibin + 1); // scale by 1/<N(φ_{i})>
0334               }
0335             }
0336
0337             for (unsigned int order = 0; order < m_MaxOrder; ++order) {
0338               double Cosine = cos(tile_phi * (double)(order + 1));
0339               double Sine = sin(tile_phi * (double)(order + 1));
0340
0341               // Arm-specific Q-vectors
0342               if (arm == 0) {
0343                 south_q[order][0] += truncated_e * Cosine;
0344                 south_q[order][1] += truncated_e * Sine;
0345                 ring_q_south[rbin][order][0] += TileWeight * Cosine;
0346                 ring_q_south[rbin][order][1] += TileWeight * Sine;
0347
0348               } else if (arm == 1) {
0349                 north_q[order][0] += truncated_e * Cosine;
0350                 north_q[order][1] += truncated_e * Sine;
0351                 ring_q_north[rbin][order][0] += TileWeight * Cosine;
0352                 ring_q_north[rbin][order][1] += TileWeight * Sine;
0353               }
0354
0355               // Combined Q-vectors
0356               northsouth_q[order][0] += truncated_e * Cosine;
0357               northsouth_q[order][1] += truncated_e * Sine;
0358             }
0359           }
0360         }
0361
0362         _totalcharge = _nsum + _ssum;
0363
0364         // Get recentering histograms and do recentering
0365         // Recentering: subtract Qn,x and Qn,y values averaged over all events
0366         for (unsigned int order = 0; order < m_MaxOrder; order++) {
0367           if (tprof_mean_cos_south_epd_input[order]) // check if recentering
0368                                                      // histograms exist
0369           {
0370             // south
0371             TAxis *south_xaxis =
0372                 tprof_mean_cos_south_epd_input[order]->GetXaxis();
0373             TAxis *south_yaxis =
0374                 tprof_mean_cos_south_epd_input[order]->GetYaxis();
0375             int xbin_south = south_xaxis->FindBin(_ssum);
0376             int ybin_south = south_yaxis->FindBin(_mbdvtx);
0377
0378             double event_ave_cos_south =
0379                 tprof_mean_cos_south_epd_input[order]->GetBinContent(
0380                     xbin_south, ybin_south);
0381             double event_ave_sin_south =
0382                 tprof_mean_sin_south_epd_input[order]->GetBinContent(
0383                     xbin_south, ybin_south);
0384             south_q_subtract[order][0] = _ssum * event_ave_cos_south;
0385             south_q_subtract[order][1] = _ssum * event_ave_sin_south;
0386             south_q[order][0] -= south_q_subtract[order][0];
0387             south_q[order][1] -= south_q_subtract[order][1];
0388
0389             // north
0390             TAxis *north_xaxis =
0391                 tprof_mean_cos_north_epd_input[order]->GetXaxis();
0392             TAxis *north_yaxis =
0393                 tprof_mean_cos_north_epd_input[order]->GetYaxis();
0394             int xbin_north = north_xaxis->FindBin(_nsum);
0395             int ybin_north = north_yaxis->FindBin(_mbdvtx);
0396
0397             double event_ave_cos_north =
0398                 tprof_mean_cos_north_epd_input[order]->GetBinContent(
0399                     xbin_north, ybin_north);
0400             double event_ave_sin_north =
0401                 tprof_mean_sin_north_epd_input[order]->GetBinContent(
0402                     xbin_north, ybin_north);
0403             north_q_subtract[order][0] = _nsum * event_ave_cos_north;
0404             north_q_subtract[order][1] = _nsum * event_ave_sin_north;
0405             north_q[order][0] -= north_q_subtract[order][0];
0406             north_q[order][1] -= north_q_subtract[order][1];
0407
0408             // northsouth
0409             TAxis *northsouth_xaxis =
0410                 tprof_mean_cos_northsouth_epd_input[order]->GetXaxis();
0411             TAxis *northsouth_yaxis =
0412                 tprof_mean_cos_northsouth_epd_input[order]->GetYaxis();
0413             int xbin_northsouth = northsouth_xaxis->FindBin(_totalcharge);
0414             int ybin_northsouth = northsouth_yaxis->FindBin(_mbdvtx);
0415
0416             double event_ave_cos_northsouth =
0417                 tprof_mean_cos_northsouth_epd_input[order]->GetBinContent(
0418                     xbin_northsouth, ybin_northsouth);
0419             double event_ave_sin_northsouth =
0420                 tprof_mean_sin_northsouth_epd_input[order]->GetBinContent(
0421                     xbin_northsouth, ybin_northsouth);
0422             northsouth_q_subtract[order][0] =
0423                 _totalcharge * event_ave_cos_northsouth;
0424             northsouth_q_subtract[order][1] =
0425                 _totalcharge * event_ave_sin_northsouth;
0426             northsouth_q[order][0] -= northsouth_q_subtract[order][0];
0427             northsouth_q[order][1] -= northsouth_q_subtract[order][1];
0428           }
0429         }
0430
0431         // Get recentered psi_n
0432         Eventplaneinfo *epinfo = new Eventplaneinfov1();
0433         for (unsigned int order = 0; order < m_MaxOrder; order++) {
0434           double n = order + 1.0;
0435           if (tprof_mean_cos_south_epd_input[order]) // if present, Qs are
0436                                                      // recentered
0437           {
0438             tmp_south_psi[order] =
0439                 epinfo->GetPsi(south_q[order][0], south_q[order][1], n);
0440             tmp_north_psi[order] =
0441                 epinfo->GetPsi(north_q[order][0], north_q[order][1], n);
0442             tmp_northsouth_psi[order] = epinfo->GetPsi(
0443                 northsouth_q[order][0], northsouth_q[order][1], n);
0444           } else {
0445             tmp_south_psi[order] = NAN;
0446             tmp_north_psi[order] = NAN;
0447             tmp_northsouth_psi[order] = NAN;
0448           }
0449         }
0450
0451         // Get shifting histograms and calculate shift
0452         for (unsigned int order = 0; order < m_MaxOrder; order++) {
0453           for (int p = 0; p < _imax; p++) {
0454             if (tprof_cos_south_epd_shift_input[order][p]) // check if shifting
0455                                                            // histograms exist
0456             {
0457               double terms = p + 1.0;
0458               double n = order + 1.0;
0459               double tmp = n * terms;
0460               double prefactor = 2.0 / terms;
0461
0462               // south
0463               TAxis *south_xaxis =
0464                   tprof_cos_south_epd_shift_input[order][p]->GetXaxis();
0465               TAxis *south_yaxis =
0466                   tprof_cos_south_epd_shift_input[order][p]->GetYaxis();
0467               int xbin_south = south_xaxis->FindBin(_ssum);
0468               int ybin_south = south_yaxis->FindBin(_mbdvtx);
0469
0470               // north
0471               TAxis *north_xaxis =
0472                   tprof_cos_north_epd_shift_input[order][p]->GetXaxis();
0473               TAxis *north_yaxis =
0474                   tprof_cos_north_epd_shift_input[order][p]->GetYaxis();
0475               int xbin_north = north_xaxis->FindBin(_nsum);
0476               int ybin_north = north_yaxis->FindBin(_mbdvtx);
0477
0478               //                   // northsouth
0479               TAxis *northsouth_xaxis =
0480                   tprof_cos_northsouth_epd_shift_input[order][p]->GetXaxis();
0481               TAxis *northsouth_yaxis =
0482                   tprof_cos_northsouth_epd_shift_input[order][p]->GetYaxis();
0483               int xbin_northsouth = northsouth_xaxis->FindBin(_totalcharge);
0484               int ybin_northsouth = northsouth_yaxis->FindBin(_mbdvtx);
0485
0486               //                    Equation (6) of arxiv:nucl-ex/9805001
0487               //                    i = terms; n = order; i*n = tmp
0488               //                    (2 / i ) * <cos(i*n*psi_n)> * sin(i*n*psi_n)
0489               //                    - <sin(i*n*psi_n)> * cos(i*n*psi_n)
0490
0491               // north
0492               shift_north[order] +=
0493                   prefactor *
0494                   (tprof_cos_north_epd_shift_input[order][p]->GetBinContent(
0495                        xbin_north, ybin_north) *
0496                        sin(tmp * tmp_north_psi[order]) -
0497                    tprof_sin_north_epd_shift_input[order][p]->GetBinContent(
0498                        xbin_north, ybin_north) *
0499                        cos(tmp * tmp_north_psi[order]));
0500
0501               // south
0502               shift_south[order] +=
0503                   prefactor *
0504                   (tprof_cos_south_epd_shift_input[order][p]->GetBinContent(
0505                        xbin_south, ybin_south) *
0506                        sin(tmp * tmp_south_psi[order]) -
0507                    tprof_sin_south_epd_shift_input[order][p]->GetBinContent(
0508                        xbin_south, ybin_south) *
0509                        cos(tmp * tmp_south_psi[order]));
0510
0511               //                     // northsouth
0512               shift_northsouth[order] +=
0513                   prefactor *
0514                   (tprof_cos_northsouth_epd_shift_input[order][p]
0515                            ->GetBinContent(xbin_northsouth, ybin_northsouth) *
0516                        sin(tmp * tmp_northsouth_psi[order]) -
0517                    tprof_sin_northsouth_epd_shift_input[order][p]
0518                            ->GetBinContent(xbin_northsouth, ybin_northsouth) *
0519                        cos(tmp * tmp_northsouth_psi[order]));
0520             }
0521           }
0522         }
0523
0524         // n * deltapsi_n = (2 / i ) * <cos(i*n*psi_n)> * sin(i*n*psi_n) -
0525         // <sin(i*n*psi_n)> * cos(i*n*psi_n) Divide out n
0526         for (unsigned int order = 0; order < m_MaxOrder; order++) {
0527           double n = order + 1.0;
0528           shift_north[order] /= n;
0529           shift_south[order] /= n;
0530           shift_northsouth[order] /= n;
0531         }
0532
0533         // Now add shift to psi_n to flatten it
0534         for (unsigned int order = 0; order < m_MaxOrder; order++) {
0535           if (tprof_cos_north_epd_shift_input[0][0]) {
0536             tmp_south_psi[order] += shift_south[order];
0537             tmp_north_psi[order] += shift_north[order];
0538             tmp_northsouth_psi[order] += shift_northsouth[order];
0539           }
0540         }
0541
0542         // Now enforce the range
0543         for (unsigned int order = 0; order < m_MaxOrder; order++) {
0544           if (tprof_cos_north_epd_shift_input[0][0]) {
0545             double range = M_PI / (double)(order + 1);
0546             if (tmp_south_psi[order] < -1.0 * range) {
0547               tmp_south_psi[order] += 2.0 * range;
0548             }
0549             if (tmp_south_psi[order] > range) {
0550               tmp_south_psi[order] -= 2.0 * range;
0551             }
0552             if (tmp_north_psi[order] < -1.0 * range) {
0553               tmp_north_psi[order] += 2.0 * range;
0554             }
0555             if (tmp_north_psi[order] > range) {
0556               tmp_north_psi[order] -= 2.0 * range;
0557             }
0558             if (tmp_northsouth_psi[order] < -1.0 * range) {
0559               tmp_northsouth_psi[order] += 2.0 * range;
0560             }
0561             if (tmp_northsouth_psi[order] > range) {
0562               tmp_northsouth_psi[order] -= 2.0 * range;
0563             }
0564           }
0565         }
0566
0567         for (unsigned int order = 0; order < m_MaxOrder; order++) {
0568           south_Qvec.emplace_back(south_q[order][0], south_q[order][1]);
0569           north_Qvec.emplace_back(north_q[order][0], north_q[order][1]);
0570           northsouth_Qvec.emplace_back(northsouth_q[order][0],
0571                                        northsouth_q[order][1]);
0572         }
0573
0574         for (int rbin = 0; rbin < nRings; ++rbin) {
0575           for (unsigned int order = 0; order < m_MaxOrder; ++order) {
0576             all_ring_Qvecs_north[rbin][order] = std::make_pair(
0577                 ring_q_north[rbin][order][0], ring_q_north[rbin][order][1]);
0578             all_ring_Qvecs_south[rbin][order] = std::make_pair(
0579                 ring_q_south[rbin][order][0], ring_q_south[rbin][order][1]);
0580           }
0581         }
0582
0583         if (epd_towerinfo) {
0584           Eventplaneinfo *sepds = new Eventplaneinfov1();
0585           sepds->set_qvector(south_Qvec);
0586           sepds->set_shifted_psi(tmp_south_psi);
0587           epmap->insert(sepds, EventplaneinfoMap::sEPDS);
0588
0589           Eventplaneinfo *sepdn = new Eventplaneinfov1();
0590           sepdn->set_qvector(north_Qvec);
0591           sepdn->set_shifted_psi(tmp_north_psi);
0592           epmap->insert(sepdn, EventplaneinfoMap::sEPDN);
0593
0594           Eventplaneinfo *sepdns = new Eventplaneinfov1();
0595           sepdns->set_qvector(northsouth_Qvec);
0596           sepdns->set_shifted_psi(tmp_northsouth_psi);
0597           epmap->insert(sepdns, EventplaneinfoMap::sEPDNS);
0598
0599           Eventplaneinfo *epring_south = new Eventplaneinfov1();
0600           epring_south->set_ring_qvector(all_ring_Qvecs_south);
0601           epmap->insert(epring_south, EventplaneinfoMap::sEPDRING_SOUTH);
0602
0603           Eventplaneinfo *epring_north = new Eventplaneinfov1();
0604           epring_north->set_ring_qvector(all_ring_Qvecs_north);
0605           epmap->insert(epring_north, EventplaneinfoMap::sEPDRING_NORTH);
0606
0607           if (Verbosity() > 1) {
0608             sepds->identify();
0609             sepdn->identify();
0610             sepdns->identify();
0611             epring_south->identify();
0612             epring_north->identify();
0613           }
0614         }
0615       }
0616     }
0617   }
0618
0619   if (_mbdEpReco) {
0620     ResetMe();
0621
0622     MbdPmtContainer *mbdpmts =
0623         findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0624     if (!mbdpmts) {
0625       std::cout << PHWHERE << "::ERROR - cannot find MbdPmtContainer"
0626                 << std::endl;
0627       exit(-1);
0628     }
0629
0630     MbdGeom *mbdgeom = findNode::getClass<MbdGeom>(topNode, "MbdGeom");
0631     if (!mbdgeom) {
0632       std::cout << PHWHERE << "::ERROR - cannot find MbdGeom" << std::endl;
0633       exit(-1);
0634     }
0635
0636     if (mbdpmts) {
0637       if (Verbosity()) {
0638         std::cout << "EventPlaneReco::process_event -  mbdpmts" << std::endl;
0639       }
0640
0641       for (int ipmt = 0; ipmt < mbdpmts->get_npmt(); ipmt++) {
0642         float mbd_q = mbdpmts->get_pmt(ipmt)->get_q();
0643         _mbdQ += mbd_q;
0644       }
0645
0646       for (int ipmt = 0; ipmt < mbdpmts->get_npmt(); ipmt++) {
0647         float mbd_q = mbdpmts->get_pmt(ipmt)->get_q();
0648         float phi = mbdgeom->get_phi(ipmt);
0649         int arm = mbdgeom->get_arm(ipmt);
0650
0651         if (_mbdQ < _mbd_e) {
0652           continue;
0653         }
0654
0655         if (arm == 0) {
0656           for (unsigned int order = 0; order < m_MaxOrder; order++) {
0657             double Cosine = cos(phi * (double)(order + 1));
0658             double Sine = sin(phi * (double)(order + 1));
0659             south_q[order][0] += mbd_q * Cosine; // south Qn,x
0660             south_q[order][1] += mbd_q * Sine;   // south Qn,y
0661           }
0662         } else if (arm == 1) {
0663           for (unsigned int order = 0; order < m_MaxOrder; order++) {
0664             double Cosine = cos(phi * (double)(order + 1));
0665             double Sine = sin(phi * (double)(order + 1));
0666             north_q[order][0] += mbd_q * Cosine; // north Qn,x
0667             north_q[order][1] += mbd_q * Sine;   // north Qn,y
0668           }
0669         }
0670       }
0671     }
0672
0673     for (unsigned int order = 0; order < m_MaxOrder; order++) {
0674       south_Qvec.emplace_back(south_q[order][0], south_q[order][1]);
0675       north_Qvec.emplace_back(north_q[order][0], north_q[order][1]);
0676     }
0677
0678     if (mbdpmts) {
0679       Eventplaneinfo *mbds = new Eventplaneinfov1();
0680       mbds->set_qvector(south_Qvec);
0681       epmap->insert(mbds, EventplaneinfoMap::MBDS);
0682
0683       Eventplaneinfo *mbdn = new Eventplaneinfov1();
0684       mbdn->set_qvector(north_Qvec);
0685       epmap->insert(mbdn, EventplaneinfoMap::MBDN);
0686
0687       if (Verbosity() > 1) {
0688         mbds->identify();
0689         mbdn->identify();
0690       }
0691     }
0692
0693     ResetMe();
0694   }
0695
0696   if (Verbosity()) {
0697     epmap->identify();
0698   }
0699
0700   return Fun4AllReturnCodes::EVENT_OK;
0701 }
0702
0703 int EventPlaneReco::CreateNodes(PHCompositeNode *topNode) {
0704   PHNodeIterator iter(topNode);
0705
0706   PHCompositeNode *dstNode =
0707       dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0708   if (!dstNode) {
0709     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0710     return Fun4AllReturnCodes::ABORTRUN;
0711   }
0712
0713   PHCompositeNode *globalNode = dynamic_cast<PHCompositeNode *>(
0714       iter.findFirst("PHCompositeNode", "GLOBAL"));
0715   if (!globalNode) {
0716     globalNode = new PHCompositeNode("GLOBAL");
0717     dstNode->addNode(globalNode);
0718   }
0719
0720   EventplaneinfoMap *eps =
0721       findNode::getClass<EventplaneinfoMap>(topNode, "EventplaneinfoMap");
0722   if (!eps) {
0723     eps = new EventplaneinfoMapv1();
0724     PHIODataNode<PHObject> *EpMapNode =
0725         new PHIODataNode<PHObject>(eps, "EventplaneinfoMap", "PHObject");
0726     globalNode->addNode(EpMapNode);
0727   }
0728   return Fun4AllReturnCodes::EVENT_OK;
0729 }
0730
0731 void EventPlaneReco::ResetMe() {
0732   for (auto &vec : south_q) {
0733     std::fill(vec.begin(), vec.end(), 0.);
0734   }
0735
0736   for (auto &vec : north_q) {
0737     std::fill(vec.begin(), vec.end(), 0.);
0738   }
0739
0740   for (auto &vec : northsouth_q) {
0741     std::fill(vec.begin(), vec.end(), 0.);
0742   }
0743
0744   for (auto &order_vec : ring_q_north) {
0745     for (auto &xy_vec : order_vec) {
0746       std::fill(xy_vec.begin(), xy_vec.end(), 0.0);
0747     }
0748   }
0749
0750   for (auto &order_vec : ring_q_south) {
0751     for (auto &xy_vec : order_vec) {
0752       std::fill(xy_vec.begin(), xy_vec.end(), 0.0);
0753     }
0754   }
0755
0756   south_Qvec.clear();
0757   north_Qvec.clear();
0758   northsouth_Qvec.clear();
0759
0760   for (auto &ring : all_ring_Qvecs_north) {
0761     for (auto &q : ring) {
0762       q = {0.0, 0.0};
0763     }
0764   }
0765
0766   for (auto &ring : all_ring_Qvecs_south) {
0767     for (auto &q : ring) {
0768       q = {0.0, 0.0};
0769     }
0770   }
0771
0772   for (auto &vec : south_q_subtract) {
0773     std::fill(vec.begin(), vec.end(), 0.);
0774   }
0775
0776   for (auto &vec : north_q_subtract) {
0777     std::fill(vec.begin(), vec.end(), 0.);
0778   }
0779
0780   for (auto &vec : northsouth_q_subtract) {
0781     std::fill(vec.begin(), vec.end(), 0.);
0782   }
0783
0784   std::fill(shift_north.begin(), shift_north.end(), 0.);
0785   std::fill(shift_south.begin(), shift_south.end(), 0.);
0786   std::fill(shift_northsouth.begin(), shift_northsouth.end(), 0.);
0787
0788   std::fill(tmp_south_psi.begin(), tmp_south_psi.end(), NAN);
0789   std::fill(tmp_north_psi.begin(), tmp_north_psi.end(), NAN);
0790   std::fill(tmp_northsouth_psi.begin(), tmp_northsouth_psi.end(), NAN);
0791
0792   _nsum = 0.;
0793   _ssum = 0.;
0794   _do_ep = false;
0795   _mbdQ = 0.;
0796   _totalcharge = 0.;
0797 }
0798
0799 int EventPlaneReco::End(PHCompositeNode * /*topNode*/) {
0800
0801   std::cout << " EventPlaneReco::End() " << std::endl;
0802   return Fun4AllReturnCodes::EVENT_OK;
0803 }