CMSSW Reference Manual

00001 /*
00002  * \file EcalEndcapDigisValidation.cc
00003  *
00004  * $Date: 2008/10/29 10:54:11 $
00005  * $Revision: 1.20 $
00006  * \author F. Cossutti
00007  *
00008 */
00009 
00010 #include <Validation/EcalDigis/interface/EcalEndcapDigisValidation.h>
00011 #include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialConditionRetriever.h"
00012 #include "DQMServices/Core/interface/DQMStore.h"
00013 
00014 using namespace cms;
00015 using namespace edm;
00016 using namespace std;
00017 
00018 EcalEndcapDigisValidation::EcalEndcapDigisValidation(const ParameterSet& ps):
00019   EEdigiCollection_(ps.getParameter<edm::InputTag>("EEdigiCollection"))
00020   {
00021  
00022   // verbosity switch
00023   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00024                                                                                                                                            
00025   dbe_ = 0;
00026                                                                                                                                           
00027   // get hold of back-end interface
00028   dbe_ = Service<DQMStore>().operator->();
00029                                                                                                                                           
00030   if ( dbe_ ) {
00031     if ( verbose_ ) {
00032       dbe_->setVerbose(1);
00033     } else {
00034       dbe_->setVerbose(0);
00035     }
00036   }
00037                                                                                                                                           
00038   if ( dbe_ ) {
00039     if ( verbose_ ) dbe_->showDirStructure();
00040   }
00041 
00042   gainConv_[1] = 1.;
00043   gainConv_[2] = 2.;
00044   gainConv_[3] = 12.;
00045   gainConv_[0] = 12.;  // saturated channels
00046   barrelADCtoGeV_ = 0.035;
00047   endcapADCtoGeV_ = 0.06;
00048  
00049   meEEDigiOccupancyzp_ = 0;
00050   meEEDigiOccupancyzm_ = 0;
00051  
00052   meEEDigiMultiplicityzp_ = 0;
00053   meEEDigiMultiplicityzm_ = 0;
00054 
00055   meEEDigiADCGlobal_ = 0;
00056 
00057   for (int i = 0; i < 10 ; i++ ) {
00058     meEEDigiADCAnalog_[i] = 0;
00059     meEEDigiADCgS_[i]  = 0;
00060     meEEDigiADCg1_[i]  = 0;
00061     meEEDigiADCg6_[i]  = 0;
00062     meEEDigiADCg12_[i] = 0;
00063     meEEDigiGain_[i] = 0;
00064   }
00065 
00066   meEEPedestal_ = 0;
00067                                  
00068   meEEMaximumgt100ADC_ = 0; 
00069                                  
00070   meEEMaximumgt20ADC_ = 0; 
00071 
00072   meEEnADCafterSwitch_ = 0;
00073  
00074   Char_t histo[200];
00075  
00076   
00077   if ( dbe_ ) {
00078     dbe_->setCurrentFolder("EcalDigisV/EcalDigiTask");
00079   
00080     sprintf (histo, "EcalDigiTask Endcap occupancy z+" ) ;
00081     meEEDigiOccupancyzp_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00082     
00083     sprintf (histo, "EcalDigiTask Endcap occupancy z-" ) ;
00084     meEEDigiOccupancyzm_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00085   
00086     sprintf (histo, "EcalDigiTask Endcap multiplicity z+" ) ;
00087     meEEDigiMultiplicityzp_ = dbe_->book1D(histo, histo, 100, 0., 7324.);
00088     
00089     sprintf (histo, "EcalDigiTask Endcap multiplicity z-" ) ;
00090     meEEDigiMultiplicityzm_ = dbe_->book1D(histo, histo, 100, 0., 7324.);
00091     
00092     sprintf (histo, "EcalDigiTask Endcap global pulse shape" ) ;
00093     meEEDigiADCGlobal_ = dbe_->bookProfile(histo, histo, 10, 0, 10, 10000, 0., 1000.) ;
00094 
00095     for (int i = 0; i < 10 ; i++ ) {
00096 
00097       sprintf (histo, "EcalDigiTask Endcap analog pulse %02d", i+1) ;
00098       meEEDigiADCAnalog_[i] = dbe_->book1D(histo, histo, 4000, 0., 400.);
00099 
00100       sprintf (histo, "EcalDigiTask Endcap ADC pulse %02d Gain 0 - Saturated", i+1) ;
00101       meEEDigiADCgS_[i] = dbe_->book1D(histo, histo, 4096, -0.5, 4095.5);
00102 
00103       sprintf (histo, "EcalDigiTask Endcap ADC pulse %02d Gain 1", i+1) ;
00104       meEEDigiADCg1_[i] = dbe_->book1D(histo, histo, 4096, -0.5, 4095.5);
00105 
00106       sprintf (histo, "EcalDigiTask Endcap ADC pulse %02d Gain 6", i+1) ;
00107       meEEDigiADCg6_[i] = dbe_->book1D(histo, histo, 4096, -0.5, 4095.5);
00108 
00109       sprintf (histo, "EcalDigiTask Endcap ADC pulse %02d Gain 12", i+1) ;
00110       meEEDigiADCg12_[i] = dbe_->book1D(histo, histo, 4096, -0.5, 4095.5);
00111 
00112       sprintf (histo, "EcalDigiTask Endcap gain pulse %02d", i+1) ;
00113       meEEDigiGain_[i] = dbe_->book1D(histo, histo, 4, 0, 4);
00114     }
00115     
00116     sprintf (histo, "EcalDigiTask Endcap pedestal for pre-sample" ) ;
00117     meEEPedestal_ = dbe_->book1D(histo, histo, 4096, -0.5, 4095.5) ;
00118 
00119     sprintf (histo, "EcalDigiTask Endcap maximum position gt 100 ADC" ) ;
00120     meEEMaximumgt100ADC_ = dbe_->book1D(histo, histo, 10, 0., 10.) ;
00121 
00122     sprintf (histo, "EcalDigiTask Endcap maximum position gt 20 ADC" ) ;
00123     meEEMaximumgt20ADC_ = dbe_->book1D(histo, histo, 10, 0., 10.) ;
00124 
00125     sprintf (histo, "EcalDigiTask Endcap ADC counts after gain switch" ) ;
00126     meEEnADCafterSwitch_ = dbe_->book1D(histo, histo, 10, 0., 10.) ;
00127 
00128   }
00129  
00130 }
00131 
00132 EcalEndcapDigisValidation::~EcalEndcapDigisValidation(){
00133  
00134 }
00135 
00136 void EcalEndcapDigisValidation::beginJob(const EventSetup& c){
00137 
00138   checkCalibrations(c);
00139 
00140 }
00141 
00142 void EcalEndcapDigisValidation::endJob(){
00143 
00144 }
00145 
00146 void EcalEndcapDigisValidation::analyze(const Event& e, const EventSetup& c){
00147 
00148   //LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
00149 
00150   Handle<EEDigiCollection> EcalDigiEE;
00151 
00152   e.getByLabel( EEdigiCollection_ , EcalDigiEE );
00153 
00154   // Return if no Endcap data available
00155   if( !EcalDigiEE.isValid() ) return;
00156 
00157   // ENDCAP
00158 
00159   // loop over Digis
00160 
00161   const EEDigiCollection * endcapDigi = EcalDigiEE.product () ;
00162 
00163   std::vector<double> eeAnalogSignal ;
00164   std::vector<double> eeADCCounts ;
00165   std::vector<double> eeADCGains ;
00166   eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
00167   eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
00168   eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
00169 
00170   int nDigiszp = 0;
00171   int nDigiszm = 0;
00172 
00173   for (unsigned int digis=0; digis<EcalDigiEE->size(); ++digis) {
00174     
00175     EEDataFrame eedf=(*endcapDigi)[digis];
00176     int nrSamples=eedf.size();
00177     
00178     EEDetId eeid = eedf.id () ;
00179     
00180     if (eeid.zside() > 0 ) {
00181       if (meEEDigiOccupancyzp_) meEEDigiOccupancyzp_->Fill( eeid.ix(), eeid.iy() );
00182       nDigiszp++;
00183     }
00184     else if (eeid.zside() < 0 ) {
00185       if (meEEDigiOccupancyzm_) meEEDigiOccupancyzm_->Fill( eeid.ix(), eeid.iy() );
00186       nDigiszm++;
00187     }
00188     
00189     double Emax = 0. ;
00190     int Pmax = 0 ;
00191     double pedestalPreSample = 0.;
00192     double pedestalPreSampleAnalog = 0.;
00193     int countsAfterGainSwitch = -1;
00194     double higherGain = 1.;
00195     int higherGainSample = 0;
00196     
00197     for (int sample = 0 ; sample < nrSamples; ++sample) {
00198       eeAnalogSignal[sample] = 0.;
00199       eeADCCounts[sample] = 0.;
00200       eeADCGains[sample] = 0.;
00201     }
00202     
00203     for (int sample = 0 ; sample < nrSamples; ++sample) {
00204       
00205       EcalMGPASample mySample = eedf[sample];
00206       
00207       eeADCCounts[sample] = (mySample.adc());
00208       eeADCGains[sample]  = (mySample.gainId()) ;
00209       eeAnalogSignal[sample] = (eeADCCounts[sample]*gainConv_[(int)eeADCGains[sample]]*endcapADCtoGeV_);
00210       
00211       if (Emax < eeAnalogSignal[sample] ) {
00212         Emax = eeAnalogSignal[sample] ;
00213         Pmax = sample ;
00214       }
00215       
00216       if ( sample < 3 ) {
00217         pedestalPreSample += eeADCCounts[sample] ;
00218         pedestalPreSampleAnalog += eeADCCounts[sample]*gainConv_[(int)eeADCGains[sample]]*endcapADCtoGeV_ ;
00219       }
00220 
00221       if (sample > 0 && ( ((eeADCGains[sample] > eeADCGains[sample-1]) && (eeADCGains[sample-1]!=0)) || (countsAfterGainSwitch<0 && eeADCGains[sample]==0)) ){  
00222         higherGain = eeADCGains[sample];
00223         higherGainSample = sample;
00224         countsAfterGainSwitch = 1;
00225       }
00226       
00227       if ( (higherGain > 1 && (higherGainSample != sample) && (eeADCGains[sample] == higherGain)) || (higherGain==3 && (higherGainSample != sample) && (eeADCGains[sample]==0)) || (higherGain==0 && (higherGainSample != sample) && ((eeADCGains[sample]==0) || (eeADCGains[sample]==3))) ) countsAfterGainSwitch++ ;
00228     }
00229     pedestalPreSample /= 3. ; 
00230     pedestalPreSampleAnalog /= 3. ; 
00231     
00232     LogDebug("DigiInfo") << "Endcap Digi for EEDetId = " << eeid.rawId() << " x,y " << eeid.ix() << " " << eeid.iy() ;
00233     for ( int i = 0; i < 10 ; i++ ) {
00234       LogDebug("DigiInfo") << "sample " << i << " ADC = " << eeADCCounts[i] << " gain = " << eeADCGains[i] << " Analog = " << eeAnalogSignal[i] ;
00235     }
00236     LogDebug("DigiInfo") << "Maximum energy = " << Emax << " in sample " << Pmax << " Pedestal from pre-sample = " << pedestalPreSampleAnalog;  
00237     if ( countsAfterGainSwitch > 0 ) LogDebug("DigiInfo") << "Counts after switch " << countsAfterGainSwitch; 
00238     
00239     if ( countsAfterGainSwitch > 0 && countsAfterGainSwitch < 5 ) {
00240       edm::LogWarning("DigiWarning") << "Wrong number of counts after gain switch before next switch! " << countsAfterGainSwitch ;
00241       for ( int i = 0; i < 10 ; i++ ) {
00242         edm::LogWarning("DigiWarning") << "sample " << i << " ADC = " << eeADCCounts[i] << " gain = " << eeADCGains[i] << " Analog = " << eeAnalogSignal[i];
00243       }
00244     }
00245     
00246     for ( int i = 0 ; i < 10 ; i++ ) {
00247       if (meEEDigiADCGlobal_ && (Emax-pedestalPreSampleAnalog*gainConv_[(int)eeADCGains[Pmax]]) > 100.*endcapADCtoGeV_) meEEDigiADCGlobal_->Fill( i , eeAnalogSignal[i] ) ;
00248       if (meEEDigiADCAnalog_[i]) meEEDigiADCAnalog_[i]->Fill( eeAnalogSignal[i] ) ;
00249       if ( eeADCGains[i] == 0 ) {
00250         if (meEEDigiADCgS_[i]) meEEDigiADCgS_[i]->Fill( eeADCCounts[i] ) ;
00251       }
00252       else if ( eeADCGains[i] == 3 ) {
00253         if (meEEDigiADCg1_[i]) meEEDigiADCg1_[i]->Fill( eeADCCounts[i] ) ;
00254       }
00255       else if ( eeADCGains[i] == 2 ) {
00256         if (meEEDigiADCg6_[i]) meEEDigiADCg6_[i]->Fill( eeADCCounts[i] ) ;
00257       }
00258       else if ( eeADCGains[i] == 1 ) {
00259         if (meEEDigiADCg12_[i]) meEEDigiADCg12_[i]->Fill( eeADCCounts[i] ) ;
00260       }
00261       if (meEEDigiGain_[i]) meEEDigiGain_[i]->Fill( eeADCGains[i] ) ;
00262     }
00263     
00264     if (meEEPedestal_) meEEPedestal_->Fill ( pedestalPreSample ) ;
00265     if (meEEMaximumgt20ADC_ && (Emax-pedestalPreSampleAnalog*gainConv_[(int)eeADCGains[Pmax]]) > 20.*endcapADCtoGeV_) meEEMaximumgt20ADC_->Fill( Pmax ) ;
00266     if (meEEMaximumgt100ADC_ && (Emax-pedestalPreSampleAnalog*gainConv_[(int)eeADCGains[Pmax]]) > 100.*endcapADCtoGeV_) meEEMaximumgt100ADC_->Fill( Pmax ) ;
00267     if (meEEnADCafterSwitch_) meEEnADCafterSwitch_->Fill(countsAfterGainSwitch);
00268     
00269   } 
00270   
00271   if ( meEEDigiMultiplicityzp_ ) meEEDigiMultiplicityzp_->Fill(nDigiszp);
00272   if ( meEEDigiMultiplicityzm_ ) meEEDigiMultiplicityzm_->Fill(nDigiszm);
00273   
00274 }
00275 
00276 void  EcalEndcapDigisValidation::checkCalibrations(const edm::EventSetup & eventSetup) 
00277 {
00278 
00279   // ADC -> GeV Scale
00280   edm::ESHandle<EcalADCToGeVConstant> pAgc;
00281   eventSetup.get<EcalADCToGeVConstantRcd>().get(pAgc);
00282   const EcalADCToGeVConstant* agc = pAgc.product();
00283   
00284   EcalMGPAGainRatio * defaultRatios = new EcalMGPAGainRatio();
00285 
00286   gainConv_[1] = 1.;
00287   gainConv_[2] = defaultRatios->gain12Over6() ;
00288   gainConv_[3] = gainConv_[2]*(defaultRatios->gain6Over1()) ;
00289   gainConv_[0] = gainConv_[2]*(defaultRatios->gain6Over1()) ;
00290 
00291   LogDebug("EcalDigi") << " Gains conversions: " << "\n" << " g0 = " << gainConv_[0] << "\n" << " g1 = " << gainConv_[1] << "\n" << " g2 = " << gainConv_[2] << "\n" << " g3 = " << gainConv_[3]; 
00292 
00293   delete defaultRatios;
00294 
00295   const double barrelADCtoGeV_  = agc->getEBValue();
00296   LogDebug("EcalDigi") << " Barrel GeV/ADC = " << barrelADCtoGeV_;
00297   const double endcapADCtoGeV_ = agc->getEEValue();
00298   LogDebug("EcalDigi") << " Endcap GeV/ADC = " << endcapADCtoGeV_;
00299 
00300 }
00301 
00302