TPedValues.cc

Go to the documentation of this file.

#include "CalibCalorimetry/EcalPedestalOffsets/interface/TPedValues.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "TAxis.h"
#include "TF1.h"
#include "TGraphErrors.h"
#include <cassert>
#include <cmath>
#include <iostream>

void reset(double vett[256]) {
  for (int i = 0; i < 256; ++i)
    vett[i] = 0.;
}

TPedValues::TPedValues(double RMSmax, int bestPedestal) : m_bestPedestal(bestPedestal), m_RMSmax(RMSmax) {
  LogDebug("EcalPedOffset") << "entering TPedValues ctor ...";
  for (int i = 0; i < 1700; ++i)
    endcapCrystalNumbers[i] = 0;
}

TPedValues::TPedValues(const TPedValues &orig) {
  LogDebug("EcalPedOffset") << "entering TPedValues copyctor ...";
  m_bestPedestal = orig.m_bestPedestal;
  m_RMSmax = orig.m_RMSmax;

  for (int gain = 0; gain < 3; ++gain)
    for (int crystal = 0; crystal < 1700; ++crystal)
      for (int DAC = 0; DAC < 256; ++DAC)
        m_entries[gain][crystal][DAC] = orig.m_entries[gain][crystal][DAC];

  for (int i = 0; i < 1700; ++i)
    endcapCrystalNumbers[i] = orig.endcapCrystalNumbers[i];
}

TPedValues::~TPedValues() {}

void TPedValues::insert(const int gainId, const int crystal, const int DAC, const int pedestal, const int endcapIndex) {
  //  assert (gainId > 0) ;
  //  assert (gainId < 4) ;
  if (gainId <= 0 || gainId >= 4) {
    edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : gainId " << gainId << " does not exist, entry skipped";
    return;
  }
  //  assert (crystal > 0) ;
  //  assert (crystal <= 1700) ;
  if (crystal <= 0 || crystal > 1700) {
    edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : crystal " << crystal
                                     << " does not exist, entry skipped";
    return;
  }
  //  assert (DAC >= 0) ;
  //  assert (DAC < 256) ;
  if (DAC < 0 || DAC >= 256) {
    edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : DAC value " << DAC << " is out range, entry skipped";
    return;
  }
  m_entries[gainId - 1][crystal - 1][DAC].insert(pedestal);
  endcapCrystalNumbers[crystal - 1] = endcapIndex;
  return;
}

TPedResult TPedValues::terminate(const int &DACstart, const int &DACend) const {
  assert(DACstart >= 0);
  assert(DACend <= 256);
  //  checkEntries (DACstart, DACend) ;

  TPedResult bestDAC;
  for (int gainId = 1; gainId < 4; ++gainId) {
    for (int crystal = 0; crystal < 1700; ++crystal) {
      double delta = 1000;
      int dummyBestDAC = -1;
      bool hasDigis = false;
      for (int DAC = DACstart; DAC < DACend; ++DAC) {
        double average = m_entries[gainId - 1][crystal][DAC].average();
        if (average == -1)
          continue;
        hasDigis = true;
        if (m_entries[gainId - 1][crystal][DAC].RMSSq() > m_RMSmax * m_RMSmax)
          continue;
        if (fabs(average - m_bestPedestal) < delta && average > 1) {
          delta = fabs(average - m_bestPedestal);
          dummyBestDAC = DAC;
        }
      }  

      bestDAC.m_DACvalue[gainId - 1][crystal] = dummyBestDAC;

      if ((dummyBestDAC == (DACend - 1) || dummyBestDAC == -1) && hasDigis) {
        int gainHuman;
        if (gainId == 1)
          gainHuman = 12;
        else if (gainId == 2)
          gainHuman = 6;
        else if (gainId == 3)
          gainHuman = 1;
        else
          gainHuman = -1;
        edm::LogError("EcalPedOffset") << " TPedValues :  cannot find best DAC value for channel: "
                                       << endcapCrystalNumbers[crystal] << " gain: " << gainHuman;
      }

    }  // loop over crystals
  }  // loop over gains
  return bestDAC;
}

int TPedValues::checkEntries(const int &DACstart, const int &DACend) const {
  assert(DACstart >= 0);
  assert(DACend <= 256);
  int returnCode = 0;
  for (int gainId = 1; gainId < 4; ++gainId) {
    for (int crystal = 0; crystal < 1700; ++crystal) {
      for (int DAC = DACstart; DAC < DACend; ++DAC) {
        double average = m_entries[gainId - 1][crystal][DAC].average();
        if (average == -1) {
          ++returnCode;
          /*
                            std::cerr << "[TPedValues][checkEntries] WARNING!"
                                      << "\tgainId " << gainId
                                      << "\tcrystal " << crystal+1
                                      << "\tDAC " << DAC
                                      << " : pedestal measurement missing"
                                      << std::endl ;
          */
        }
        /*
                      std::cout << "[pietro][RMS]: " <<
           m_entries[gainId-1][crystal][DAC].RMS ()     //FIXME
                                << "\t" <<
           m_entries[gainId-1][crystal][DAC].RMSSq () //FIXME
                                << "\t" << DAC //FIXME
                                << "\t" << gainId //FIXME
                                << "\t" << crystal << std::endl ; //FIXME
        */
      }  
    }  // loop over crystals
  }  // loop over gains
  return returnCode;
}

int TPedValues::makePlots(TFile *rootFile,
                          const std::string &dirName,
                          const double maxSlope,
                          const double minSlope,
                          const double maxChi2OverNDF) const {
  using namespace std;
  // prepare the ROOT file
  if (!rootFile->cd(dirName.c_str())) {
    rootFile->mkdir(dirName.c_str());
    rootFile->cd(dirName.c_str());
  }

  // loop over the crystals
  for (int xtl = 0; xtl < 1700; ++xtl) {
    // loop over the gains
    for (int gain = 0; gain < 3; ++gain) {
      vector<double> asseX;
      vector<double> sigmaX;
      vector<double> asseY;
      vector<double> sigmaY;
      asseX.reserve(256);
      sigmaX.reserve(256);
      asseY.reserve(256);
      sigmaY.reserve(256);
      // loop over DAC values
      for (int dac = 0; dac < 256; ++dac) {
        double average = m_entries[gain][xtl][dac].average();
        if (average > -1) {
          double rms = m_entries[gain][xtl][dac].RMS();
          asseX.push_back(dac);
          sigmaX.push_back(0);
          asseY.push_back(average);
          sigmaY.push_back(rms);
        }
      }  // loop over DAC values
      if (!asseX.empty()) {
        int lastBin = 0;
        while (lastBin < (int)asseX.size() - 1 && asseY[lastBin + 1] > 0 &&
               (asseY[lastBin + 1] - asseY[lastBin + 2]) != 0)
          lastBin++;

        int fitRangeEnd = (int)asseX[lastBin];
        int kinkPt = 64;
        if (fitRangeEnd < 66)
          kinkPt = fitRangeEnd - 4;
        TGraphErrors graph(asseX.size(), &(*asseX.begin()), &(*asseY.begin()), &(*sigmaX.begin()), &(*sigmaY.begin()));
        char funct[120];
        sprintf(funct, "(x<%d)*([0]*x+[1])+(x>=%d)*([2]*x+[3])", kinkPt, kinkPt);
        TF1 fitFunction("fitFunction", funct, asseX[0], fitRangeEnd);
        fitFunction.SetLineColor(2);

        char name[120];
        int gainHuman;
        if (gain == 0)
          gainHuman = 12;
        else if (gain == 1)
          gainHuman = 6;
        else if (gain == 2)
          gainHuman = 1;
        else
          gainHuman = -1;
        sprintf(name, "XTL%04d_GAIN%02d", endcapCrystalNumbers[xtl], gainHuman);
        graph.GetXaxis()->SetTitle("DAC value");
        graph.GetYaxis()->SetTitle("Average pedestal ADC");
        graph.Fit(&fitFunction, "RWQ");
        graph.Write(name);

        double slope1 = fitFunction.GetParameter(0);
        double slope2 = fitFunction.GetParameter(2);

        if (fitFunction.GetChisquare() / fitFunction.GetNDF() > maxChi2OverNDF ||
            fitFunction.GetChisquare() / fitFunction.GetNDF() < 0 || slope1 > 0 || slope2 > 0 ||
            ((slope1 < -29 || slope1 > -18) && slope1 < 0) || ((slope2 < -29 || slope2 > -18) && slope2 < 0)) {
          edm::LogError("EcalPedOffset") << "TPedValues : TGraph for channel:" << endcapCrystalNumbers[xtl]
                                         << " gain:" << gainHuman << " is not linear;"
                                         << "  slope of line1:" << fitFunction.GetParameter(0)
                                         << " slope of line2:" << fitFunction.GetParameter(2) << " reduced chi-squared:"
                                         << fitFunction.GetChisquare() / fitFunction.GetNDF();
        }
        // LogDebug("EcalPedOffset") << "TPedValues : TGraph for channel:" <<
        // xtl+1 << " gain:"
        //  << gainHuman << " has " << asseX.size() << " points...back is:" <<
        //  asseX.back()
        //  << " and front+1 is:" << asseX.front()+1;
        if ((asseX.back() - asseX.front() + 1) != asseX.size())
          edm::LogError("EcalPedOffset") << "TPedValues : Pedestal average not found "
                                         << "for all DAC values scanned in channel:" << endcapCrystalNumbers[xtl]
                                         << " gain:" << gainHuman;
      }
    }  // loop over the gains
  }  // (loop over the crystals)

  return 0;
}

// Look up the crystal number in the EE schema and return it
int TPedValues::getCrystalNumber(int xtal) const { return endcapCrystalNumbers[xtal]; }