packages/TpcLineLaserDiode/DiodeReco.cc

0001 #include "DiodeReco.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/PHIODataNode.h>    // for PHIODataNode
0006 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0007 #include <phool/PHObject.h>        // for PHObject
0008 #include <phool/getClass.h>
0009
0010 #include <ffarawobjects/TpcDiodeContainerv1.h>
0011 #include <ffarawobjects/TpcDiodev1.h>
0012
0013 #include <Event/Event.h>
0014 #include <Event/caen_correction.h>
0015 #include <Event/packet.h>
0016
0017 #include <ffamodules/CDBInterface.h>
0018
0019 #include <TSystem.h>
0020
0021 #include <algorithm>
0022 #include <cstdint>                             // for uint16_t
0023 #include <cstdlib>                             // for exit, size_t
0024 #include <iostream>                             // for basic_ostream, operat...
0025
0026 DiodeReco::DiodeReco(const std::string &name)
0027   : SubsysReco(name)
0028   , m_DiodeContainerName("TPCDIODES")
0029 {
0030   for (int c = 0; c < 32; c++)
0031   {
0032     adc.clear();
0033   }
0034 }
0035
0036 int DiodeReco::InitRun(PHCompositeNode *topNode)
0037 {
0038    std::cout << "DiodeReco::InitRun(PHCompositeNode *topNode) Initializing" << std::endl;
0039
0040   m_cdb = CDBInterface::instance();
0041
0042   PHNodeIterator iter(topNode);
0043   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0044   if (!dstNode)
0045   {
0046     dstNode = new PHCompositeNode("DST");
0047     topNode->addNode(dstNode);
0048   }
0049
0050   PHNodeIterator iterDst(dstNode);
0051   PHCompositeNode *detNode = dynamic_cast<PHCompositeNode *>(iterDst.findFirst("PHCompositeNode", "TPC"));
0052   if (!detNode)
0053   {
0054     detNode = new PHCompositeNode("TPC");
0055     dstNode->addNode(detNode);
0056   }
0057
0058   TpcDiodeContainer *tpcdiodecont = findNode::getClass<TpcDiodeContainer>(detNode, m_DiodeContainerName);
0059   if (!tpcdiodecont)
0060   {
0061     tpcdiodecont = new TpcDiodeContainerv1();
0062     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(tpcdiodecont, m_DiodeContainerName, "PHObject");
0063     detNode->addNode(newNode);
0064   }
0065
0066   // char name[100];
0067   // char title[100];
0068
0069   // c_persistency_N = new TCanvas("c_persistency_N","c_persistency_N",800,600);
0070   // c_persistency_N->Divide(4,4);
0071   // c_persistency_S = new TCanvas("c_persistency_S","c_persistency_S",800,600);
0072   // c_persistency_S->Divide(4,4);
0073   // c_waveforms = new TCanvas("c_waveforms","c_waveforms",800,600);
0074
0075   // for(int c=0;c<32;c++){
0076   //   sprintf(name,"CAEN Persistency Waveform (Channel %d)",c);
0077   //   sprintf(title,"CAEN Persistency Waveform (Channel %d);time bins [ns];ADC",c);
0078   //   persistency[c] = new TH2F(name,title,1024,-0.5,1023.5,256,-250,4250);
0079   // }
0080   // waveforms = new TH2F("CAEN Waveforms","CAEN Waveforms;time bins [ns];channels",1024,0,1024,32,0,32);
0081   return Fun4AllReturnCodes::EVENT_OK;
0082 }
0083
0084 int DiodeReco::process_event(PHCompositeNode *topNode)
0085 {
0086   Event *_event = findNode::getClass<Event>(topNode, "PRDF");
0087   if (!_event)
0088   {
0089     std::cout << "RawdataUnpacker::Process_Event - Event not found" << std::endl;
0090     return Fun4AllReturnCodes::DISCARDEVENT;
0091   }
0092
0093   diodes = findNode::getClass<TpcDiodeContainer>(topNode, "TPCDIODES");
0094
0095   //_event->identify();
0096   caen_correction *cc{nullptr};
0097   Packet *p = _event->getPacket(2000);
0098   if (p)
0099   {
0100     if (!cc)
0101     {
0102       switch (p->iValue(0, "FREQUENCY"))
0103       {
0104       case 0:  // 5GS
0105     calibdir = m_cdb->getUrl("TPC_CAEN_CORRECTION_24974_5G");
0106         cc = new caen_correction(calibdir.c_str());
0107         break;
0108       case 1:  // 2.5GS
0109         calibdir = m_cdb->getUrl("TPC_CAEN_CORRECTION_24974_2G");
0110         cc = new caen_correction(calibdir.c_str());
0111         break;
0112       case 2:  // 1GS
0113         calibdir = m_cdb->getUrl("TPC_CAEN_CORRECTION_24974_1G");
0114         cc = new caen_correction(calibdir.c_str());
0115         break;
0116       default:
0117         std::cout << "Bad selection " << p->iValue(0, "FREQUENCY") << std::endl;
0118         gSystem->Exit(1);
0119         exit(1);
0120         break;
0121       }
0122     }
0123
0124     if (cc)
0125     {
0126       cc->init(p);
0127     }
0128
0129     TpcDiodev1 *newdiode = new TpcDiodev1();
0130
0131     for (int c = 0; c < 32; c++)
0132     {
0133       for (int t = 0; t < 1024; t++)
0134       {
0135         adc.push_back(cc->caen_corrected(t, c));
0136       }
0137
0138       PedestalCorrected(0, 200);
0139       uint16_t maxadc = MaxAdc(2, 200, 300);
0140       uint16_t maxbin = MaxBin(2);
0141       double integral = Integral(0, 1024);
0142       uint16_t nabovethreshold = NAboveThreshold(200, 100);
0143       double pulsewidth = PulseWidth(200, 100);
0144       const uint16_t samples = adc.size();
0145
0146       uint16_t packetid = 2000;
0147       uint16_t channel = c;
0148
0149       newdiode->set_packetid(packetid);
0150       newdiode->set_channel(channel);
0151       newdiode->set_maxadc(maxadc);
0152       newdiode->set_maxbin(maxbin);
0153       newdiode->set_integral(integral);
0154       newdiode->set_nabovethreshold(nabovethreshold);
0155       newdiode->set_pulsewidth(pulsewidth);
0156       newdiode->set_samples(samples);
0157
0158       for (uint16_t s = 0; s < samples; s++)
0159       {
0160         uint16_t adcval = adc[s];
0161         newdiode->set_adc(s, adcval);
0162       }
0163
0164       diodes->AddDiode(newdiode);
0165       adc.clear();
0166     }
0167
0168     // if(event==2){
0169     //    for(int c=0;c<32;c++)
0170     //      {
0171     //        for(int t=0;t<1024;t++){
0172     //      adc.push_back(cc->caen_corrected(t,c));
0173     //        }
0174
0175     //        PedestalCorrected(0,200);
0176
0177     //        if(nlaser>1)
0178     //      {
0179     //        std::cout << "More than one laser fired! Stopping code!" << std::endl;
0180     //        return -1;
0181     //      }
0182
0183     //        if(laser<0)
0184     //      {
0185     //        std::cout << "No laser fired in this event!" << std::endl;
0186     //      }
0187     //        for(int s=0;s<static_cast<int>(adc.size());s++){
0188     //      persistency[c]->Fill(s,adc[s]);
0189     //      waveforms->Fill(s,c,adc[s]);
0190     //        }
0191     //        if(c<16){
0192     //      c_persistency_N->cd(c+1);
0193     //      persistency[c]->Draw("colz");
0194     //        }
0195     //        else
0196     //      {
0197     //        c_persistency_S->cd(c-15);
0198     //        persistency[c]->Draw("colz");
0199     //      }
0200     //        adc.clear();
0201     //      }
0202     // }
0203     // p->dump();
0204     std::cout << diodes->get_Laser() << std::endl;
0205     // for(int c=0;c<32;c++){
0206     //  std::cout << diodes->get_diode(c)->get_maxadc() << std::endl;
0207     // }
0208     delete p;
0209   }
0210
0211   // c_waveforms->cd();
0212   // waveforms->Draw("LEGO");
0213
0214   return Fun4AllReturnCodes::EVENT_OK;
0215 }
0216
0217 int DiodeReco::End(PHCompositeNode *)
0218 {
0219   return Fun4AllReturnCodes::EVENT_OK;
0220 }
0221
0222 double DiodeReco::MaxAdc(int n, int low_bin, int high_bin)
0223 {
0224   double MaxSum = -99999;  // Maximum sum over n bins within the bin range
0225   int MaxMid = -1;         // Bin number of the middle bin used to calculate MaxSum
0226
0227   // Bracket the limits safely...
0228   int start = std::max(1 + n, low_bin + n);
0229   int end = std::min(int(adc.size()) - n, high_bin - n);
0230   for (int mid = start; mid < end; mid++)
0231   {
0232     double Sum = 0;
0233     for (int bin = mid - n; bin <= mid + n; bin++)
0234     {
0235       Sum += adc[bin];
0236     }
0237     if (Sum > MaxSum)
0238     {
0239       MaxSum = Sum;
0240       MaxMid = mid;
0241     }
0242   }
0243
0244   if (MaxMid < 0)
0245   {
0246     std::cout << "Error: Maximum ADC could not be found!" << std::endl;
0247   }
0248   return MaxSum / (2.0 * n + 1.0);
0249 }
0250
0251 int DiodeReco::MaxBin(int n)
0252 {
0253   double MaxSum = -99999;
0254   int MaxMid = -1;
0255   for (int mid = 1 + n; mid < static_cast<int>(adc.size()) - n; mid++)
0256   {
0257     double Sum = 0;
0258     for (int bin = mid - n; bin <= mid + n; bin++)
0259     {
0260       Sum += adc[bin];
0261     }
0262     if (Sum > MaxSum)
0263     {
0264       MaxSum = Sum;
0265       MaxMid = mid;
0266     }
0267   }
0268   if (MaxMid < 0)
0269   {
0270     std::cout << "Error: Maximum ADC could not be found!" << std::endl;
0271   }
0272   return MaxMid;
0273 }
0274
0275 double DiodeReco::Integral(int low_bin, int high_bin)
0276 {
0277   low_bin = std::max(low_bin, 0);
0278   high_bin = std::min(high_bin, static_cast<int>(adc.size()));
0279
0280   double SUM = 0;
0281   for (int i = low_bin; i < high_bin; i++)
0282   {
0283     SUM += adc[i];
0284   }
0285   return SUM;
0286 }
0287
0288 int DiodeReco::NAboveThreshold(double upper_thr, double lower_thr)
0289 {
0290   int nAbove = 0;
0291
0292   bool belowThreshold = true;
0293
0294   for (size_t i = 0; i < adc.size(); i++)
0295   {
0296     if (belowThreshold && adc[i] >= upper_thr)
0297     {
0298       nAbove++;
0299       belowThreshold = false;
0300     }
0301
0302     else if (!belowThreshold && adc[i] < lower_thr)
0303     {
0304       belowThreshold = true;
0305     }
0306   }
0307
0308   return nAbove;
0309 }
0310
0311 double DiodeReco::PulseWidth(double upper_thr, double lower_thr)
0312 {
0313   //  The results of this routine are ONLY valid
0314   //  if NAbove is one.
0315
0316   bool belowThreshold = true;
0317
0318   int left = 0;
0319   int right = 0;
0320
0321   for (int i = 0; i < static_cast<int>(adc.size()); i++)
0322   {
0323     if (belowThreshold && adc[i] >= upper_thr)
0324     {
0325       left = i;
0326       belowThreshold = false;
0327     }
0328
0329     else if (!belowThreshold && adc[i] < lower_thr)
0330     {
0331       right = i;
0332       belowThreshold = true;
0333     }
0334   }
0335
0336   return right - left;
0337 }
0338
0339 void DiodeReco::PedestalCorrected(int low_bin = -1, int high_bin = 9999)
0340 {
0341   low_bin = std::max(low_bin, 0);
0342   if (high_bin > static_cast<int>(adc.size()))
0343   {
0344     high_bin = adc.size();
0345   }
0346
0347   // Copy all voltages in the pedestal region & sort
0348   std::vector<double> sam;
0349   for (int i = low_bin; i < high_bin; i++)
0350   {
0351     sam.push_back(adc[i]);
0352   }
0353   sort(sam.begin(), sam.end());
0354
0355   // Assign the pedestal as the median of this distribution
0356   int n = sam.size();
0357   double PEDESTAL;
0358   if (n % 2 != 0)
0359   {
0360     PEDESTAL = sam[n / 2];
0361   }
0362   else
0363   {
0364     PEDESTAL = (sam[(n - 1) / 2] + sam[n / 2]) / 2.0;
0365   }
0366
0367   for (size_t i = 0; i < adc.size(); i++)
0368   {
0369     adc[i] = adc[i] - PEDESTAL;
0370   }
0371 }