Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:20:03 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!