pileuprejection/src/PileupRejector.cc

0001 #include "PileupRejector.h"
0002
0003 #include <mbd/MbdPmtContainer.h>
0004 #include <mbd/MbdPmtHit.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/getClass.h>
0007 #include <cmath>
0008 int PileupRejector::decodeEvent(PHCompositeNode *topNode)
0009 {
0010
0011   MbdPmtContainer *pmts_mbd = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0012
0013   float maxtime[2] = {-100, -100};
0014   float mintime[2] = {100, 100};
0015   float rmstime[2] = {0,0};
0016   float sumtime[2] = {0,0};
0017   float npmt[2] = {0,0};
0018
0019   m_chargesum = 0;
0020   m_prodsigma = 0;
0021   m_avgsigma = 0;
0022   m_maxsigma = 0;
0023   m_proddelta = 0;
0024   m_avgdelta = 0;
0025   m_maxdelta = 0;
0026
0027   if (pmts_mbd)
0028     {
0029       for (int i = 0 ; i < 128; i++)
0030     {
0031       MbdPmtHit *tmp_pmt = pmts_mbd->get_pmt(i);
0032
0033       float charge = tmp_pmt->get_q();
0034       float time = tmp_pmt->get_time();
0035
0036       if (fabs(time) < 25. && charge > 0.4)
0037         {
0038           m_chargesum += charge;
0039           int side = i/64;
0040           npmt[side]++;
0041           if (time > maxtime[side])
0042         maxtime[side] = time;
0043           if (time < mintime[side])
0044         mintime[side] = time;
0045
0046               rmstime[side] += time*time;
0047               sumtime[side] += time;
0048
0049         }
0050     }
0051       if (npmt[0] < 1 && npmt[1] < 1)
0052     {
0053       m_pileup = false;
0054       return 0;
0055     }
0056
0057       if (npmt[0] >= m_hitcut)
0058         {
0059           rmstime[0]/=(float)npmt[0];
0060           sumtime[0]/=(float)npmt[0];
0061           rmstime[0] = sqrt(rmstime[0] -sumtime[0]*sumtime[0]);
0062     }
0063       else
0064     {
0065       maxtime[0] = 0;
0066       mintime[0] = 0;
0067       rmstime[0] = 0;
0068       sumtime[0] = 0;
0069     }
0070       if (npmt[1] >= m_hitcut)
0071         {
0072           rmstime[1]/=(float)npmt[1];
0073           sumtime[1]/=(float)npmt[1];
0074           rmstime[1] = sqrt(rmstime[1] -sumtime[1]*sumtime[1]);
0075     }
0076       else
0077     {
0078       maxtime[1] = 0;
0079       mintime[1] = 0;
0080       rmstime[1] = 0;
0081       sumtime[1] = 0;
0082     }
0083
0084       m_prodsigma = rmstime[0]*rmstime[1];
0085       m_avgsigma = (rmstime[0] + rmstime[1])/2.;
0086       m_maxsigma = std::max(rmstime[0], rmstime[1]);
0087
0088       m_proddelta = (maxtime[1] - mintime[1])*(maxtime[0] - mintime[0]);
0089       m_avgdelta = ((maxtime[1] - mintime[1]) + (maxtime[0] - mintime[0]))/2.;
0090       m_maxdelta = std::max(maxtime[1] - mintime[1], maxtime[0] - mintime[0]);
0091
0092
0093     }
0094
0095   return 0;
0096 }
0097
0098 bool PileupRejector::isPileup()
0099 {
0100   m_pileup = false;
0101   if (m_cutstrength == PileupCutStrength::COMFORT)
0102     {
0103       if (m_prodsigma >= m_comfort_time_cut) m_pileup = true;
0104     }
0105   if (m_cutstrength == PileupCutStrength::STRICT)
0106     {
0107       if (m_avgsigma >= m_strict_time_cut) m_pileup = true;
0108     }
0109   if (m_cutstrength == PileupCutStrength::DRACONIAN)
0110     {
0111       if (m_avgsigma >= m_draconian_time_cut) m_pileup = true;
0112     }
0113   return m_pileup;
0114 }
0115 void PileupRejector::Print()
0116 {
0117
0118   std::cout << " PileupRejector: " << std::endl;
0119   std::cout << "    Is Pileup : " << (m_pileup ? "Yes" : "No" ) << std::endl;
0120
0121 }
0122