packages/TpcLineLaserDiode/DiodeReco.cc

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