CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/CalibCalorimetry/EcalPedestalOffsets/src/TPedValues.cc

Go to the documentation of this file.
00001 #include "CalibCalorimetry/EcalPedestalOffsets/interface/TPedValues.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 
00004 #include <math.h>
00005 #include <iostream>
00006 #include <cassert>
00007 #include "TGraphErrors.h"
00008 #include "TAxis.h"
00009 #include "TF1.h"
00010 
00011 void reset (double vett[256]) 
00012 {
00013    for (int i=0 ; i<256; ++i) vett[i] = 0. ; 
00014 }
00015 
00016 
00017 TPedValues::TPedValues (double RMSmax, int bestPedestal) :
00018   m_bestPedestal (bestPedestal) ,
00019   m_RMSmax (RMSmax) 
00020 {
00021   LogDebug ("EcalPedOffset") << "entering TPedValues ctor ..." ;
00022   for(int i=0; i<1700;++i)
00023     endcapCrystalNumbers[i] = 0;
00024 }
00025 
00026 
00027 TPedValues::TPedValues (const TPedValues & orig) 
00028 {
00029   LogDebug ("EcalPedOffset") << "entering TPedValues copyctor ..." ;
00030   m_bestPedestal = orig.m_bestPedestal ;
00031   m_RMSmax = orig.m_RMSmax ;
00032 
00033   for (int gain = 0 ; gain < 3 ; ++gain)
00034     for (int crystal = 0 ; crystal < 1700 ; ++crystal)
00035         for (int DAC = 0 ; DAC < 256 ; ++DAC)
00036           m_entries[gain][crystal][DAC] = orig.m_entries[gain][crystal][DAC] ;
00037 
00038   for(int i=0; i<1700;++i)
00039     endcapCrystalNumbers[i] = orig.endcapCrystalNumbers[i];
00040 }
00041 
00042 
00043 
00044 TPedValues::~TPedValues () {}
00045   
00046 
00047 void TPedValues::insert (const int gainId, 
00048                          const int crystal, 
00049                          const int DAC, 
00050                          const int pedestal,
00051                          const int endcapIndex) 
00052 {
00053 //  assert (gainId > 0) ;
00054 //  assert (gainId < 4) ;
00055   if (gainId <= 0 || gainId >= 4)
00056     {
00057       edm::LogWarning ("EcalPedOffset") << "WARNING : TPedValues : gainId " << gainId
00058                                         << " does not exist, entry skipped" ;
00059       return ;    
00060     }
00061 //  assert (crystal > 0) ;
00062 //  assert (crystal <= 1700) ;
00063   if (crystal <= 0 || crystal > 1700)
00064     {
00065       edm::LogWarning ("EcalPedOffset") << "WARNING : TPedValues : crystal " << crystal
00066                                         << " does not exist, entry skipped" ;
00067       return ;    
00068     }
00069 //  assert (DAC >= 0) ; 
00070 //  assert (DAC < 256) ;
00071   if (DAC < 0 || DAC >= 256)
00072     {
00073       edm::LogWarning ("EcalPedOffset") << "WARNING : TPedValues : DAC value " << DAC
00074                                         << " is out range, entry skipped" ;
00075       return ;    
00076     }
00077   m_entries[gainId-1][crystal-1][DAC].insert (pedestal) ;
00078   endcapCrystalNumbers[crystal-1] = endcapIndex;
00079   return ;
00080 }
00081     
00082 
00083 TPedResult TPedValues::terminate (const int & DACstart, const int & DACend) const
00084 {
00085   assert (DACstart >= 0) ;
00086   assert (DACend <= 256) ;
00087 //  checkEntries (DACstart, DACend) ;
00088 
00089   TPedResult bestDAC ;
00091   for (int gainId = 1 ; gainId < 4 ; ++gainId)
00092     {
00094       for (int crystal = 0 ; crystal < 1700 ; ++crystal)
00095         {
00097           double delta = 1000 ;
00098           int dummyBestDAC = -1 ;
00099           bool hasDigis = false;
00101           for (int DAC = DACstart ; DAC < DACend ; ++DAC)
00102             {
00103               double average = m_entries[gainId-1][crystal][DAC].average () ;
00104               if (average == -1) continue ;
00105               hasDigis = true;
00106               if (m_entries[gainId-1][crystal][DAC].RMSSq () > m_RMSmax * m_RMSmax) continue ;
00107               if (fabs (average - m_bestPedestal) < delta &&   average>1 ) 
00108                 {
00109                   delta = fabs (average - m_bestPedestal) ;
00110                   dummyBestDAC = DAC ;
00111                 }
00112             } 
00113 
00114           bestDAC.m_DACvalue[gainId-1][crystal] = dummyBestDAC ;
00115           
00116           if ((dummyBestDAC == (DACend-1) || dummyBestDAC == -1) && hasDigis)
00117           {
00118             int gainHuman;
00119             if      (gainId ==1) gainHuman =12;
00120             else if (gainId ==2) gainHuman =6;
00121             else if (gainId ==3) gainHuman =1;
00122             else                 gainHuman =-1;
00123             edm::LogError("EcalPedOffset")
00124               << " TPedValues :  cannot find best DAC value for channel: "
00125               << endcapCrystalNumbers[crystal]
00126               << " gain: " << gainHuman;
00127           }
00128           
00129         } // loop over crystals
00130     } // loop over gains
00131   return bestDAC ;
00132 }
00133     
00134 
00135 int TPedValues::checkEntries (const int & DACstart, const int & DACend) const
00136 {
00137   assert (DACstart >= 0) ;
00138   assert (DACend <= 256) ;
00139   int returnCode = 0 ;
00141   for (int gainId = 1 ; gainId < 4 ; ++gainId)
00142     {
00144       for (int crystal = 0 ; crystal < 1700 ; ++crystal)
00145         {
00147           for (int DAC = DACstart ; DAC < DACend ; ++DAC)
00148             {
00149               double average = m_entries[gainId-1][crystal][DAC].average () ;
00150               if (average == -1) 
00151                 {
00152                   ++returnCode ;
00154 /*
00155                   std::cerr << "[TPedValues][checkEntries] WARNING!"
00156                             << "\tgainId " << gainId
00157                             << "\tcrystal " << crystal+1
00158                             << "\tDAC " << DAC
00159                             << " : pedestal measurement missing" 
00160                             << std::endl ;
00161 */
00162                 }
00163 /*
00164               std::cout << "[pietro][RMS]: " << m_entries[gainId-1][crystal][DAC].RMS ()     //FIXME
00165                         << "\t" << m_entries[gainId-1][crystal][DAC].RMSSq () //FIXME
00166                         << "\t" << DAC //FIXME
00167                         << "\t" << gainId //FIXME
00168                         << "\t" << crystal << std::endl ; //FIXME
00169 */              
00170             } 
00171         } // loop over crystals
00172     } // loop over gains
00173   return returnCode ;
00174 }
00175 
00176 
00178 int TPedValues::makePlots (TFile * rootFile, const std::string & dirName,
00179      const double maxSlope, const double minSlope, const double maxChi2OverNDF) const 
00180 {
00181   using namespace std;
00182   // prepare the ROOT file
00183   if (!rootFile->cd (dirName.c_str ())) 
00184     {
00185       rootFile->mkdir (dirName.c_str ()) ;
00186       rootFile->cd (dirName.c_str ()) ;
00187     }
00188   
00189   // loop over the crystals
00190   for (int xtl=0 ; xtl<1700 ; ++xtl)
00191   {
00192     // loop over the gains
00193     for (int gain=0 ; gain<3 ; ++gain)
00194       {
00195         vector<double> asseX;
00196         vector<double> sigmaX;
00197         vector<double> asseY;
00198         vector<double> sigmaY;
00199         asseX.reserve(256);
00200         sigmaX.reserve(256);
00201         asseY.reserve(256);
00202         sigmaY.reserve(256);
00203         // loop over DAC values
00204         for (int dac=0 ; dac<256 ; ++dac)
00205           {
00206             double average = m_entries[gain][xtl][dac].average();
00207             if(average > -1)
00208             {
00209               double rms = m_entries[gain][xtl][dac].RMS();
00210               asseX.push_back(dac);
00211               sigmaX.push_back(0);
00212               asseY.push_back(average);
00213               sigmaY.push_back(rms);
00214             }
00215           } // loop over DAC values
00216         if(asseX.size() > 0)
00217         {
00218           int lastBin = 0;
00219           while(lastBin<(int)asseX.size()-1 && asseY[lastBin+1]>0
00220               && (asseY[lastBin+1]-asseY[lastBin+2])!=0)
00221             lastBin++;
00222           
00223           int fitRangeEnd = (int)asseX[lastBin];
00224           int kinkPt = 64;
00225           if(fitRangeEnd < 66)
00226             kinkPt = fitRangeEnd-4;
00227           TGraphErrors graph(asseX.size(),&(*asseX.begin()),&(*asseY.begin()),
00228               &(*sigmaX.begin()),&(*sigmaY.begin()));
00229           char funct[120];
00230           sprintf(funct,"(x<%d)*([0]*x+[1])+(x>=%d)*([2]*x+[3])",kinkPt,kinkPt);
00231           TF1 fitFunction("fitFunction",funct,asseX[0],fitRangeEnd);
00232           fitFunction.SetLineColor(2);
00233           
00234           char name[120] ;
00235           int gainHuman;
00236           if      (gain ==0) gainHuman =12;
00237           else if (gain ==1) gainHuman =6;
00238           else if (gain ==2) gainHuman =1;
00239           else               gainHuman =-1;
00240           sprintf (name,"XTL%04d_GAIN%02d",endcapCrystalNumbers[xtl],gainHuman) ;
00241           graph.GetXaxis()->SetTitle("DAC value");
00242           graph.GetYaxis()->SetTitle("Average pedestal ADC");
00243           graph.Fit("fitFunction","RWQ");
00244           graph.Write (name);
00245           
00246           double slope1 = fitFunction.GetParameter(0);
00247           double slope2 = fitFunction.GetParameter(2);
00248 
00249           if(fitFunction.GetChisquare()/fitFunction.GetNDF()>maxChi2OverNDF ||
00250               fitFunction.GetChisquare()/fitFunction.GetNDF()<0 ||
00251               slope1>0 || slope2>0 ||
00252               ((slope1<-29 || slope1>-18) && slope1<0) || 
00253               ((slope2<-29 || slope2>-18) && slope2<0))
00254           {
00255             edm::LogError("EcalPedOffset") << "TPedValues : TGraph for channel:" << 
00256               endcapCrystalNumbers[xtl] << 
00257               " gain:" << gainHuman << " is not linear;" << "  slope of line1:" << 
00258               fitFunction.GetParameter(0) << " slope of line2:" << 
00259               fitFunction.GetParameter(2) << " reduced chi-squared:" << 
00260               fitFunction.GetChisquare()/fitFunction.GetNDF();
00261           }
00262           //LogDebug("EcalPedOffset") << "TPedValues : TGraph for channel:" << xtl+1 << " gain:"
00263           //  << gainHuman << " has " << asseX.size() << " points...back is:" << asseX.back() 
00264           //  << " and front+1 is:" << asseX.front()+1;
00265           if((asseX.back()-asseX.front()+1)!=asseX.size())
00266             edm::LogError("EcalPedOffset") << "TPedValues : Pedestal average not found " <<
00267               "for all DAC values scanned in channel:" << endcapCrystalNumbers[xtl]
00268               << " gain:" << gainHuman;
00269         }
00270       } // loop over the gains
00271   }     // (loop over the crystals)
00272   
00273   return 0 ;
00274 }
00275 
00276 // Look up the crystal number in the EE schema and return it
00277 int TPedValues::getCrystalNumber(int xtal) const
00278 {
00279   return endcapCrystalNumbers[xtal];
00280 }