packages/zdcinfo/ZdcReco.cc

0001 #include "ZdcReco.h"
0002 #include "Zdcinfo.h"
0003 #include "Zdcinfov1.h"
0004 #include "Zdcinfov2.h"
0005
0006 #include <calobase/TowerInfo.h>
0007 #include <calobase/TowerInfoContainer.h>
0008 #include <calobase/TowerInfoDefs.h>
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/SubsysReco.h>  // for SubsysReco
0012
0013 #include <cdbobjects/CDBTTree.h>  // for CDBTTree
0014 #include <ffamodules/CDBInterface.h>
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>
0018 #include <phool/PHNode.h>  // for PHNode
0019 #include <phool/PHNodeIterator.h>
0020 #include <phool/PHObject.h>  // for PHObject
0021 #include <phool/getClass.h>
0022 #include <phool/phool.h>  // for PHWHERE
0023 #include <phool/recoConsts.h>
0024
0025 #include <TSystem.h>
0026
0027 #include <array>  // for array
0028 #include <cfloat>
0029 #include <cmath>
0030 #include <cstdlib>  // for exit
0031 #include <iostream>
0032 #include <set>      // for _Rb_tree_const_iterator
0033 #include <utility>  // for pair
0034 #include <vector>   // for vector
0035
0036 ZdcReco::ZdcReco(const std::string &name)
0037   : SubsysReco(name)
0038 {
0039 }
0040
0041 ZdcReco::~ZdcReco()
0042 {
0043   delete cdbttree;
0044 }
0045
0046 int ZdcReco::InitRun(PHCompositeNode *topNode)
0047 {
0048   if (!m_overrideCalibName)
0049   {
0050     m_calibName = "data_driven_zdc_calib";
0051   }
0052   if (!m_overrideFieldName)
0053   {
0054     m_fieldname = "zdc_calib";
0055   }
0056   std::string calibdir = CDBInterface::instance()->getUrl(m_calibName);
0057   if (!calibdir.empty())
0058   {
0059     cdbttree = new CDBTTree(calibdir);
0060   }
0061   else
0062   {
0063     std::cout << "ZdcReco::::InitRun No calibration file for domain "
0064               << m_calibName << " found" << std::endl;
0065     exit(1);
0066   }
0067
0068   PHNodeIterator node_itr(topNode);
0069   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(
0070       node_itr.findFirst("PHCompositeNode", "DST"));
0071   if (!dstNode)
0072   {
0073     std::cout << "PHComposite node created: DST" << std::endl;
0074     dstNode = new PHCompositeNode("DST");
0075     topNode->addNode(dstNode);
0076   }
0077
0078   PHNodeIterator nodeItr(dstNode);
0079   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(
0080       nodeItr.findFirst("PHCompositeNode", m_Detector));
0081   if (!DetNode)
0082   {
0083     DetNode = new PHCompositeNode(m_Detector);
0084     dstNode->addNode(DetNode);
0085   }
0086
0087   m_zdcinfo = findNode::getClass<Zdcinfo>(topNode, "Zdcinfo");
0088   if (!m_zdcinfo)
0089   {
0090     m_zdcinfo = new Zdcinfov2();
0091     PHIODataNode<PHObject> *newNode =
0092         new PHIODataNode<PHObject>(m_zdcinfo, "Zdcinfo", "PHObject");
0093     DetNode->addNode(newNode);
0094   }
0095
0096   return Fun4AllReturnCodes::EVENT_OK;
0097 }
0098
0099 int ZdcReco::process_event(PHCompositeNode *topNode)
0100 {
0101   if (Verbosity() > 1)
0102   {
0103     std::cout << "ZdcReco::process_event -- entered" << std::endl;
0104   }
0105
0106   //---------------------------------
0107   // Get Objects off of the Node Tree
0108   //---------------------------------
0109
0110   ResetMe();
0111   _z_vertex = std::numeric_limits<float>::quiet_NaN();
0112   _radius_north = std::numeric_limits<float>::quiet_NaN();
0113   _radius_south = std::numeric_limits<float>::quiet_NaN();
0114   _sumSden = 0.;
0115   _sumNden = 0.;
0116   _sumSt = 0.;
0117   _sumNt = 0.;
0118   _sumS = 0.;
0119   _sumN = 0.;
0120   _nhor = 0;
0121   _nver = 0;
0122   _shor = 0;
0123   _sver = 0;
0124
0125   TowerInfoContainer *zdc_towerinfo =
0126       findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_ZDC");
0127   if (!zdc_towerinfo)
0128   {
0129     std::cout << PHWHERE << "::ERROR - cannot find TOWERS_ZDC" << std::endl;
0130     exit(-1);
0131   }
0132
0133   if (zdc_towerinfo)
0134   {
0135     if (Verbosity())
0136     {
0137       std::cout << "ZdcReco::process_event -  zdc_towerinfo" << std::endl;
0138     }
0139
0140     unsigned int ntowers = zdc_towerinfo->size();
0141
0142     if (ntowers != 52)
0143     {
0144       std::cout << "ZdcReco::process_event -  zdc size mismatch" << std::endl;
0145       exit(1);
0146     }
0147
0148     // get zdc-smd info
0149     for (unsigned int ch = 0; ch < ntowers; ch++)
0150     {
0151       TowerInfo *_tower = zdc_towerinfo->get_tower_at_channel(ch);
0152       float zdc_e = _tower->get_energy();
0153       float zdc_time = _tower->get_time_float();
0154       if (TowerInfoDefs::isZDC(ch))
0155       {
0156         vzdcadc.push_back(zdc_e);
0157         vzdctime.push_back(zdc_time);
0158       }
0159
0160       if (TowerInfoDefs::isSMD(ch))
0161       {
0162         vsmdadc.push_back(zdc_e);
0163         vsmdtime.push_back(zdc_time);
0164       }
0165     }
0166
0167     // check smd mapping
0168     int ssize = vsmdadc.size();
0169     if (ssize != 32)
0170     {
0171       std::cout << "smd channel mapping error" << std::endl;
0172       if (vsmdtime.size() != 32)
0173       {
0174         exit(1);
0175       }
0176     }
0177
0178     // check zdc mapping
0179     int zsize = vzdcadc.size();
0180     if (zsize != 16)
0181     {
0182       std::cout << "zdc channel mapping error" << std::endl;
0183       if (vzdctime.size() != 16)
0184       {
0185         exit(1);
0186       }
0187     }
0188
0189     // apply time cuts per smd side
0190     for (int j = 0; j < ssize; j++)
0191     {
0192       if (j < 16)
0193       {
0194         if (vsmdtime[j] > 9.0 && vsmdtime[j] < 14.0)
0195         {
0196           smd_adc[j] = vsmdadc[j];
0197           if (vsmdadc[j] > _smd_e)
0198           {
0199             if (j <= 7)
0200             {
0201               _nhor++;
0202             }
0203             if (j >= 8 && j <= 14)
0204             {
0205               _nver++;
0206             }
0207           }
0208         }
0209         else
0210         {
0211           smd_adc[j] = 0.0;
0212         }
0213       }
0214       else
0215       {
0216         if (vsmdtime[j] > 6.0 && vsmdtime[j] < 12.0)
0217         {
0218           smd_adc[j] = vsmdadc[j];
0219           if (vsmdadc[j] > _smd_e)
0220           {
0221             if (j >= 16 && j <= 23)
0222             {
0223               _shor++;
0224             }
0225             if (j >= 24 && j <= 30)
0226             {
0227               _sver++;
0228             }
0229           }
0230         }
0231         else
0232         {
0233           smd_adc[j] = 0.0;
0234         }
0235       }
0236     }
0237
0238     // get smd position
0239     CompSmdPos();
0240
0241     if (_nver > 1 && _nhor > 1)
0242     {
0243       smd_north_fired = true;
0244     }
0245     if (_sver > 1 && _shor > 1)
0246     {
0247       smd_south_fired = true;
0248     }
0249
0250     if (smd_north_fired)
0251     {
0252       _radius_north =
0253           std::sqrt(smd_pos[1] * smd_pos[1] + smd_pos[0] * smd_pos[0]);
0254     }
0255     if (smd_south_fired)
0256     {
0257       _radius_south =
0258           std::sqrt(smd_pos[3] * smd_pos[3] + smd_pos[2] * smd_pos[2]);
0259     }
0260
0261     // apply time cuts per zdc and get sums
0262     for (int i = 0; i < zsize; i++)
0263     {
0264       unsigned int key = TowerInfoDefs::encode_zdc(i);
0265       int arm = TowerInfoDefs::get_zdc_side(key);
0266
0267       if (vzdctime[i] > 5.0 && vzdctime[i] < 9.0)
0268       {
0269         if (arm == 0)
0270         {
0271           if (vzdcadc[6] > 50.)
0272           {
0273             _sumSt = vzdcadc[0] * cdbttree->GetFloatValue(0, m_fieldname) *
0274                          vzdctime[0] +
0275                      vzdcadc[2] * cdbttree->GetFloatValue(2, m_fieldname) *
0276                          vzdctime[2] +
0277                      vzdcadc[4] * cdbttree->GetFloatValue(4, m_fieldname) *
0278                          vzdctime[4];
0279
0280             _sumSden = vzdcadc[0] * cdbttree->GetFloatValue(0, m_fieldname) +
0281                        vzdcadc[2] * cdbttree->GetFloatValue(2, m_fieldname) +
0282                        vzdcadc[4] * cdbttree->GetFloatValue(4, m_fieldname);
0283           }
0284
0285           if (vzdcadc[0] > _zdc1_e && vzdcadc[2] > _zdc2_e)
0286           {
0287             _sumS = vzdcadc[0] * cdbttree->GetFloatValue(0, m_fieldname) +
0288                     vzdcadc[2] * cdbttree->GetFloatValue(2, m_fieldname) +
0289                     vzdcadc[4] * cdbttree->GetFloatValue(4, m_fieldname);
0290           }
0291         }
0292         else if (arm == 1)
0293         {
0294           if (vzdcadc[14] > 50.)
0295           {
0296             _sumNt = vzdcadc[8] * cdbttree->GetFloatValue(8, m_fieldname) *
0297                          vzdctime[8] +
0298                      vzdcadc[10] * cdbttree->GetFloatValue(10, m_fieldname) *
0299                          vzdctime[10] +
0300                      vzdcadc[12] * cdbttree->GetFloatValue(12, m_fieldname) *
0301                          vzdctime[12];
0302
0303             _sumNden = vzdcadc[8] * cdbttree->GetFloatValue(8, m_fieldname) +
0304                        vzdcadc[10] * cdbttree->GetFloatValue(10, m_fieldname) +
0305                        vzdcadc[12] * cdbttree->GetFloatValue(12, m_fieldname);
0306           }
0307
0308           if (vzdcadc[8] > _zdc1_e && vzdcadc[10] > _zdc2_e)
0309           {
0310             _sumN = vzdcadc[8] * cdbttree->GetFloatValue(8, m_fieldname) +
0311                     vzdcadc[10] * cdbttree->GetFloatValue(10, m_fieldname) +
0312                     vzdcadc[12] * cdbttree->GetFloatValue(12, m_fieldname);
0313           }
0314         }
0315       }
0316     }
0317
0318     // in ns
0319     double south_time_weighted = (_sumSt / _sumSden) * _t;
0320     double north_time_weighted = (_sumNt / _sumNden) * _t;
0321     double time_diff_ns_weighted = south_time_weighted - north_time_weighted;
0322     _z_vertex = time_diff_ns_weighted * 0.5 * _c;
0323
0324     m_zdcinfo->set_zdc_energy(0, _sumS);
0325     m_zdcinfo->set_zdc_energy(1, _sumN);
0326     m_zdcinfo->set_radius(0, _radius_south);
0327     m_zdcinfo->set_radius(1, _radius_north);
0328     m_zdcinfo->set_zvertex(_z_vertex);
0329   }
0330
0331   return Fun4AllReturnCodes::EVENT_OK;
0332 }
0333
0334 void ZdcReco::ResetMe()
0335 {
0336   vsmdadc.clear();
0337   vsmdtime.clear();
0338   vzdcadc.clear();
0339   vzdctime.clear();
0340   smd_north_fired = false;
0341   smd_south_fired = false;
0342 }
0343
0344 void ZdcReco::CompSmdPos()  // computing position with weighted averages
0345 {
0346   float w_ave[4];  // 0 -> north hor; 1 -> noth vert; 2 -> south hor; 3 -> south
0347                    // vert.
0348   float weights[4] = {0};
0349   memset(weights, 0, sizeof(weights));  // memset float works only for 0
0350   float w_sum[4];
0351   memset(w_sum, 0, sizeof(w_sum));
0352
0353   // these constants convert the SMD channel number into real dimensions (cm's)
0354   const float hor_scale = 2.0 * 11.0 / 10.5 * sin(M_PI / 4);  // from gsl_math.h
0355   const float ver_scale = 1.5 * 11.0 / 10.5;
0356   float hor_offset = (hor_scale * 8 / 2.0) * (7.0 / 8.0);
0357   float ver_offset = (ver_scale * 7 / 2.0) * (6.0 / 7.0);
0358
0359   for (int i = 0; i < 8; i++)
0360   {
0361     weights[0] += smd_adc[i];  // summing weights
0362     weights[2] += smd_adc[i + 16];
0363     w_sum[0] += (float) i * smd_adc[i];  // summing for the average
0364     w_sum[2] += ((float) i + 16.) * smd_adc[i + 16];
0365   }
0366   for (int i = 0; i < 7; i++)
0367   {
0368     weights[1] += smd_adc[i + 8];
0369     weights[3] += smd_adc[i + 24];
0370     w_sum[1] += ((float) i + 8.) * smd_adc[i + 8];
0371     w_sum[3] += ((float) i + 24.) * smd_adc[i + 24];
0372   }
0373
0374   if (weights[0] > 0.0)
0375   {
0376     w_ave[0] = w_sum[0] / weights[0];  // average = sum / sumn of weights...
0377     smd_pos[0] = hor_scale * w_ave[0] - hor_offset;
0378   }
0379   else
0380   {
0381     smd_pos[0] = 0;
0382   }
0383   if (weights[1] > 0.0)
0384   {
0385     w_ave[1] = w_sum[1] / weights[1];
0386     smd_pos[1] = ver_scale * (w_ave[1] - 8.0) - ver_offset;
0387   }
0388   else
0389   {
0390     smd_pos[1] = 0;
0391   }
0392
0393   if (weights[2] > 0.0)
0394   {
0395     w_ave[2] = w_sum[2] / weights[2];
0396     smd_pos[2] = hor_scale * (w_ave[2] - 16.0) - hor_offset;
0397   }
0398   else
0399   {
0400     smd_pos[2] = 0;
0401   }
0402
0403   if (weights[3] > 0.0)
0404   {
0405     w_ave[3] = w_sum[3] / weights[3];
0406     smd_pos[3] = ver_scale * (w_ave[3] - 24.0) - ver_offset;
0407   }
0408   else
0409   {
0410     smd_pos[3] = 0;
0411   }
0412 }