Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:57:59 

0001 #include "CalibCalorimetry/EcalPedestalOffsets/interface/TPedValues.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 
0004 #include "TAxis.h"
0005 #include "TF1.h"
0006 #include "TGraphErrors.h"
0007 #include <cassert>
0008 #include <cmath>
0009 #include <iostream>
0010 
0011 void reset(double vett[256]) {
0012   for (int i = 0; i < 256; ++i)
0013     vett[i] = 0.;
0014 }
0015 
0016 TPedValues::TPedValues(double RMSmax, int bestPedestal) : m_bestPedestal(bestPedestal), m_RMSmax(RMSmax) {
0017   LogDebug("EcalPedOffset") << "entering TPedValues ctor ...";
0018   for (int i = 0; i < 1700; ++i)
0019     endcapCrystalNumbers[i] = 0;
0020 }
0021 
0022 TPedValues::TPedValues(const TPedValues &orig) {
0023   LogDebug("EcalPedOffset") << "entering TPedValues copyctor ...";
0024   m_bestPedestal = orig.m_bestPedestal;
0025   m_RMSmax = orig.m_RMSmax;
0026 
0027   for (int gain = 0; gain < 3; ++gain)
0028     for (int crystal = 0; crystal < 1700; ++crystal)
0029       for (int DAC = 0; DAC < 256; ++DAC)
0030         m_entries[gain][crystal][DAC] = orig.m_entries[gain][crystal][DAC];
0031 
0032   for (int i = 0; i < 1700; ++i)
0033     endcapCrystalNumbers[i] = orig.endcapCrystalNumbers[i];
0034 }
0035 
0036 TPedValues::~TPedValues() {}
0037 
0038 void TPedValues::insert(const int gainId, const int crystal, const int DAC, const int pedestal, const int endcapIndex) {
0039   //  assert (gainId > 0) ;
0040   //  assert (gainId < 4) ;
0041   if (gainId <= 0 || gainId >= 4) {
0042     edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : gainId " << gainId << " does not exist, entry skipped";
0043     return;
0044   }
0045   //  assert (crystal > 0) ;
0046   //  assert (crystal <= 1700) ;
0047   if (crystal <= 0 || crystal > 1700) {
0048     edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : crystal " << crystal
0049                                      << " does not exist, entry skipped";
0050     return;
0051   }
0052   //  assert (DAC >= 0) ;
0053   //  assert (DAC < 256) ;
0054   if (DAC < 0 || DAC >= 256) {
0055     edm::LogWarning("EcalPedOffset") << "WARNING : TPedValues : DAC value " << DAC << " is out range, entry skipped";
0056     return;
0057   }
0058   m_entries[gainId - 1][crystal - 1][DAC].insert(pedestal);
0059   endcapCrystalNumbers[crystal - 1] = endcapIndex;
0060   return;
0061 }
0062 
0063 TPedResult TPedValues::terminate(const int &DACstart, const int &DACend) const {
0064   assert(DACstart >= 0);
0065   assert(DACend <= 256);
0066   //  checkEntries (DACstart, DACend) ;
0067 
0068   TPedResult bestDAC;
0069   //! loop over gains
0070   for (int gainId = 1; gainId < 4; ++gainId) {
0071     //! loop over crystals
0072     for (int crystal = 0; crystal < 1700; ++crystal) {
0073       //! find the DAC value with the average pedestal nearest to 200
0074       double delta = 1000;
0075       int dummyBestDAC = -1;
0076       bool hasDigis = false;
0077       //! loop over DAC values
0078       for (int DAC = DACstart; DAC < DACend; ++DAC) {
0079         double average = m_entries[gainId - 1][crystal][DAC].average();
0080         if (average == -1)
0081           continue;
0082         hasDigis = true;
0083         if (m_entries[gainId - 1][crystal][DAC].RMSSq() > m_RMSmax * m_RMSmax)
0084           continue;
0085         if (fabs(average - m_bestPedestal) < delta && average > 1) {
0086           delta = fabs(average - m_bestPedestal);
0087           dummyBestDAC = DAC;
0088         }
0089       }  //! loop over DAC values
0090 
0091       bestDAC.m_DACvalue[gainId - 1][crystal] = dummyBestDAC;
0092 
0093       if ((dummyBestDAC == (DACend - 1) || dummyBestDAC == -1) && hasDigis) {
0094         int gainHuman;
0095         if (gainId == 1)
0096           gainHuman = 12;
0097         else if (gainId == 2)
0098           gainHuman = 6;
0099         else if (gainId == 3)
0100           gainHuman = 1;
0101         else
0102           gainHuman = -1;
0103         edm::LogError("EcalPedOffset") << " TPedValues :  cannot find best DAC value for channel: "
0104                                        << endcapCrystalNumbers[crystal] << " gain: " << gainHuman;
0105       }
0106 
0107     }  // loop over crystals
0108   }    // loop over gains
0109   return bestDAC;
0110 }
0111 
0112 int TPedValues::checkEntries(const int &DACstart, const int &DACend) const {
0113   assert(DACstart >= 0);
0114   assert(DACend <= 256);
0115   int returnCode = 0;
0116   //! loop over gains
0117   for (int gainId = 1; gainId < 4; ++gainId) {
0118     //! loop over crystals
0119     for (int crystal = 0; crystal < 1700; ++crystal) {
0120       //! loop over DAC values
0121       for (int DAC = DACstart; DAC < DACend; ++DAC) {
0122         double average = m_entries[gainId - 1][crystal][DAC].average();
0123         if (average == -1) {
0124           ++returnCode;
0125           //! do something!
0126           /*
0127                             std::cerr << "[TPedValues][checkEntries] WARNING!"
0128                                       << "\tgainId " << gainId
0129                                       << "\tcrystal " << crystal+1
0130                                       << "\tDAC " << DAC
0131                                       << " : pedestal measurement missing"
0132                                       << std::endl ;
0133           */
0134         }
0135         /*
0136                       std::cout << "[pietro][RMS]: " <<
0137            m_entries[gainId-1][crystal][DAC].RMS ()     //FIXME
0138                                 << "\t" <<
0139            m_entries[gainId-1][crystal][DAC].RMSSq () //FIXME
0140                                 << "\t" << DAC //FIXME
0141                                 << "\t" << gainId //FIXME
0142                                 << "\t" << crystal << std::endl ; //FIXME
0143         */
0144       }  //! loop over DAC values
0145     }    // loop over crystals
0146   }      // loop over gains
0147   return returnCode;
0148 }
0149 
0150 //! create a plot of the DAC pedestal trend
0151 int TPedValues::makePlots(TFile *rootFile,
0152                           const std::string &dirName,
0153                           const double maxSlope,
0154                           const double minSlope,
0155                           const double maxChi2OverNDF) const {
0156   using namespace std;
0157   // prepare the ROOT file
0158   if (!rootFile->cd(dirName.c_str())) {
0159     rootFile->mkdir(dirName.c_str());
0160     rootFile->cd(dirName.c_str());
0161   }
0162 
0163   // loop over the crystals
0164   for (int xtl = 0; xtl < 1700; ++xtl) {
0165     // loop over the gains
0166     for (int gain = 0; gain < 3; ++gain) {
0167       vector<double> asseX;
0168       vector<double> sigmaX;
0169       vector<double> asseY;
0170       vector<double> sigmaY;
0171       asseX.reserve(256);
0172       sigmaX.reserve(256);
0173       asseY.reserve(256);
0174       sigmaY.reserve(256);
0175       // loop over DAC values
0176       for (int dac = 0; dac < 256; ++dac) {
0177         double average = m_entries[gain][xtl][dac].average();
0178         if (average > -1) {
0179           double rms = m_entries[gain][xtl][dac].RMS();
0180           asseX.push_back(dac);
0181           sigmaX.push_back(0);
0182           asseY.push_back(average);
0183           sigmaY.push_back(rms);
0184         }
0185       }  // loop over DAC values
0186       if (!asseX.empty()) {
0187         int lastBin = 0;
0188         while (lastBin < (int)asseX.size() - 1 && asseY[lastBin + 1] > 0 &&
0189                (asseY[lastBin + 1] - asseY[lastBin + 2]) != 0)
0190           lastBin++;
0191 
0192         int fitRangeEnd = (int)asseX[lastBin];
0193         int kinkPt = 64;
0194         if (fitRangeEnd < 66)
0195           kinkPt = fitRangeEnd - 4;
0196         TGraphErrors graph(asseX.size(), &(*asseX.begin()), &(*asseY.begin()), &(*sigmaX.begin()), &(*sigmaY.begin()));
0197         char funct[120];
0198         sprintf(funct, "(x<%d)*([0]*x+[1])+(x>=%d)*([2]*x+[3])", kinkPt, kinkPt);
0199         TF1 fitFunction("fitFunction", funct, asseX[0], fitRangeEnd);
0200         fitFunction.SetLineColor(2);
0201 
0202         char name[120];
0203         int gainHuman;
0204         if (gain == 0)
0205           gainHuman = 12;
0206         else if (gain == 1)
0207           gainHuman = 6;
0208         else if (gain == 2)
0209           gainHuman = 1;
0210         else
0211           gainHuman = -1;
0212         sprintf(name, "XTL%04d_GAIN%02d", endcapCrystalNumbers[xtl], gainHuman);
0213         graph.GetXaxis()->SetTitle("DAC value");
0214         graph.GetYaxis()->SetTitle("Average pedestal ADC");
0215         graph.Fit(&fitFunction, "RWQ");
0216         graph.Write(name);
0217 
0218         double slope1 = fitFunction.GetParameter(0);
0219         double slope2 = fitFunction.GetParameter(2);
0220 
0221         if (fitFunction.GetChisquare() / fitFunction.GetNDF() > maxChi2OverNDF ||
0222             fitFunction.GetChisquare() / fitFunction.GetNDF() < 0 || slope1 > 0 || slope2 > 0 ||
0223             ((slope1 < -29 || slope1 > -18) && slope1 < 0) || ((slope2 < -29 || slope2 > -18) && slope2 < 0)) {
0224           edm::LogError("EcalPedOffset") << "TPedValues : TGraph for channel:" << endcapCrystalNumbers[xtl]
0225                                          << " gain:" << gainHuman << " is not linear;"
0226                                          << "  slope of line1:" << fitFunction.GetParameter(0)
0227                                          << " slope of line2:" << fitFunction.GetParameter(2) << " reduced chi-squared:"
0228                                          << fitFunction.GetChisquare() / fitFunction.GetNDF();
0229         }
0230         // LogDebug("EcalPedOffset") << "TPedValues : TGraph for channel:" <<
0231         // xtl+1 << " gain:"
0232         //  << gainHuman << " has " << asseX.size() << " points...back is:" <<
0233         //  asseX.back()
0234         //  << " and front+1 is:" << asseX.front()+1;
0235         if ((asseX.back() - asseX.front() + 1) != asseX.size())
0236           edm::LogError("EcalPedOffset") << "TPedValues : Pedestal average not found "
0237                                          << "for all DAC values scanned in channel:" << endcapCrystalNumbers[xtl]
0238                                          << " gain:" << gainHuman;
0239       }
0240     }  // loop over the gains
0241   }    // (loop over the crystals)
0242 
0243   return 0;
0244 }
0245 
0246 // Look up the crystal number in the EE schema and return it
0247 int TPedValues::getCrystalNumber(int xtal) const { return endcapCrystalNumbers[xtal]; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!