tpc/TPCPedestalCalibration/TPCPedestalCalibration.cc

0001 #include "TPCPedestalCalibration.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/PHIODataNode.h>    // for PHIODataNode
0007 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0008 #include <phool/PHObject.h>        // for PHObject
0009 #include <phool/recoConsts.h>
0010 #include <phool/getClass.h>
0011 #include <phool/phool.h>
0012
0013 #include <Event/Event.h>
0014 #include <Event/EventTypes.h>
0015 #include <Event/packet.h>
0016
0017 #include <TFile.h>
0018 #include <TTree.h>
0019
0020 #include <memory>
0021 #include <cassert>
0022 #include <iostream>
0023
0024 /*************************************************************/
0025 /*                TPC Pedestal Calibration                   */
0026 /*               Thomas Marshall,Aditya Dash                 */
0027 /*        rosstom@ucla.edu,aditya55@physics.ucla.edu         */
0028 /*************************************************************/
0029
0030 TPCPedestalCalibration::TPCPedestalCalibration(const std::string &name)
0031  :SubsysReco("TPCPedestalCalibration")
0032  , m_fname(name)
0033  , m_writeToCDB(false)
0034 {
0035   // reserve memory for max ADC samples
0036   m_adcSamples.resize(1024, 0);
0037   for(int fee_no=0;fee_no<26;fee_no++)
0038   {
0039     for(int channel_no=0;channel_no<256;channel_no++)
0040     {
0041       m_aveADCFeeChannel[fee_no][channel_no]=0.0;
0042       m_stdADCFeeChannel[fee_no][channel_no]=0.0;
0043       m_countsADCFeeChannel[fee_no][channel_no]=0.0;
0044       m_aliveArrayFeeChannel[fee_no][channel_no]=1;
0045     }
0046   }
0047
0048 }
0049
0050 int TPCPedestalCalibration::InitRun(PHCompositeNode * /*unused*/)
0051 {
0052   m_cdbttree = new CDBTTree(m_fname);
0053
0054   return Fun4AllReturnCodes::EVENT_OK;
0055 }
0056
0057 //____________________________________________________________________________..
0058 int TPCPedestalCalibration::process_event(PHCompositeNode *topNode)
0059 {
0060   Event *_event = findNode::getClass<Event>(topNode, "PRDF");
0061   if (_event == nullptr)
0062   {
0063     std::cout << "TPCRawDataTree::Process_Event - Event not found" << std::endl;
0064     return -1;
0065   }
0066   if (_event->getEvtType() >= 8)  /// special events
0067   {
0068     return Fun4AllReturnCodes::DISCARDEVENT;
0069   }
0070
0071   for (int packet : m_packets)
0072   {
0073     if (Verbosity())
0074     {
0075       std::cout << __PRETTY_FUNCTION__ << " : decoding packet " << packet << std::endl;
0076     }
0077
0078     m_packet = packet;
0079
0080     std::unique_ptr<Packet> p (_event->getPacket(m_packet));
0081     if (!p)
0082     {
0083       if (Verbosity())
0084       {
0085         std::cout << __PRETTY_FUNCTION__ << " : missing packet " << packet << std::endl;
0086       }
0087
0088       continue;
0089     }
0090
0091     m_nWaveormInFrame = p->iValue(0, "NR_WF");
0092
0093     for (int wf = 0; wf < m_nWaveormInFrame; wf++)
0094     {
0095       if (m_firstBCO == true)
0096       {
0097         m_BCO = p->iValue(wf, "BCO");
0098         m_firstBCO = false;
0099       }
0100
0101       m_nSamples = p->iValue(wf, "SAMPLES");
0102       m_fee = p->iValue(wf, "FEE");
0103       m_Channel = p->iValue(wf, "CHANNEL");
0104
0105       if(m_nSamples==0)
0106       {
0107         m_aliveArrayFeeChannel[m_fee][m_Channel]=0;
0108         continue;
0109       }
0110
0111       assert(m_nSamples < (int) m_adcSamples.size());  // no need for movements in memory allocation
0112       for (int s = 0; s < m_nSamples; s++)
0113       {
0114         m_adcSamples[s] = p->iValue(wf, s);
0115       }
0116
0117       bool dead = false;
0118       for(int adc_sam_no=0;adc_sam_no<m_nSamples;adc_sam_no++)
0119       {
0120         if (m_adcSamples[adc_sam_no] == 0 || std::isnan(float(m_adcSamples[adc_sam_no])))
0121         {
0122           m_aliveArrayFeeChannel[m_fee][m_Channel]=0;
0123         }
0124       }
0125       if (dead) {continue;}
0126
0127       for(int adc_sam_no=0;adc_sam_no<m_nSamples;adc_sam_no++){
0128         m_aveADCFeeChannel[m_fee][m_Channel]+=m_adcSamples[adc_sam_no];
0129         m_stdADCFeeChannel[m_fee][m_Channel]+=pow(m_adcSamples[adc_sam_no],2);
0130         m_countsADCFeeChannel[m_fee][m_Channel]+=1.0;
0131       }
0132     }
0133   }  //   for (int packet : m_packets)
0134
0135   return Fun4AllReturnCodes::EVENT_OK;
0136 }
0137
0138 int TPCPedestalCalibration::EndRun(const int runnumber)
0139 {
0140   std::cout << "TPCPedestalCalibration::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0141
0142   for(int fee_no=0;fee_no<26;fee_no++)
0143   {
0144     for(int channel_no=0;channel_no<256;channel_no++)
0145     {
0146       if(m_countsADCFeeChannel[fee_no][channel_no] != 0.0)
0147       {
0148         float temp1 = m_aveADCFeeChannel[fee_no][channel_no]/m_countsADCFeeChannel[fee_no][channel_no];
0149         float temp2 = m_stdADCFeeChannel[fee_no][channel_no]/m_countsADCFeeChannel[fee_no][channel_no];
0150         m_aveADCFeeChannel[fee_no][channel_no] = temp1;
0151         m_stdADCFeeChannel[fee_no][channel_no] = temp2;
0152       }
0153       else
0154       {
0155         m_aveADCFeeChannel[fee_no][channel_no] = 0.0;
0156         m_stdADCFeeChannel[fee_no][channel_no] = 0.0;
0157         m_aliveArrayFeeChannel[fee_no][channel_no]=0;
0158       }
0159
0160       if(m_aveADCFeeChannel[fee_no][channel_no] > 200 || m_aveADCFeeChannel[fee_no][channel_no] < 10)
0161       {
0162         m_aliveArrayFeeChannel[fee_no][channel_no]=0;
0163       }
0164
0165       m_pedMean=m_aveADCFeeChannel[fee_no][channel_no];
0166       m_pedStd=sqrt(m_stdADCFeeChannel[fee_no][channel_no] - pow(m_aveADCFeeChannel[fee_no][channel_no],2));
0167       m_isAlive=m_aliveArrayFeeChannel[fee_no][channel_no];
0168       m_chan=channel_no;
0169       m_outFEE=fee_no;
0170       m_module=mod_arr[fee_no];
0171       m_slot=slot_arr[fee_no];
0172
0173       m_cdbttree->SetIntValue(fee_no*256 + channel_no,"isAlive",m_isAlive);
0174       m_cdbttree->SetFloatValue(fee_no*256 + channel_no,"pedMean",m_pedMean);
0175       m_cdbttree->SetFloatValue(fee_no*256 + channel_no,"pedStd",m_pedStd);
0176       m_cdbttree->SetIntValue(fee_no*256 + channel_no,"sector",m_sector);
0177       m_cdbttree->SetIntValue(fee_no*256 + channel_no,"fee",m_outFEE);
0178       m_cdbttree->SetIntValue(fee_no*256 + channel_no,"channel",m_chan);
0179       m_cdbttree->SetIntValue(fee_no*256 + channel_no,"module",m_module);
0180       m_cdbttree->SetIntValue(fee_no*256 + channel_no,"slot",m_slot);
0181     }
0182   }
0183
0184   return Fun4AllReturnCodes::EVENT_OK;
0185 }
0186
0187 void TPCPedestalCalibration::CDBInsert()
0188 {
0189    recoConsts *rc = recoConsts::instance();
0190
0191    CDBUtils *cdbInsert = new CDBUtils(rc->get_StringFlag("CDB_GLOBALTAG"));
0192    cdbInsert->createPayloadType("TPCPedestalCalibration");
0193    cdbInsert->insertPayload("TPCPedestalCalibration",m_fname,m_BCO); // uses first BCO value from first waveform
0194
0195    return;
0196 }
0197
0198 //____________________________________________________________________________..
0199 int TPCPedestalCalibration::End(PHCompositeNode * /*topNode*/)
0200 {
0201   m_cdbttree->Commit();
0202   if (Verbosity())
0203   {
0204     m_cdbttree->Print();
0205   }
0206   m_cdbttree->WriteCDBTTree();
0207   delete m_cdbttree;
0208
0209   if (m_writeToCDB)
0210   {
0211     std::cout << "Inserting file " << m_fname << " in CDB under username " << m_username << std::endl;
0212     CDBInsert();
0213   }
0214
0215   return Fun4AllReturnCodes::EVENT_OK;
0216 }