Reference Manual

Go to the documentation of this file.
00001 #include "RecoLocalCalo/HcalRecAlgos/interface/ZdcSimpleRecAlgo.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
00004 #include <algorithm> // for "max"
00005 #include <iostream>
00006 #include <math.h>
00007 
00008 
00009 
00010 static double MaximumFractionalError = 0.0005; // 0.05% error allowed from this source
00011 
00012 ZdcSimpleRecAlgo::ZdcSimpleRecAlgo(bool correctForTimeslew, bool correctForPulse, float phaseNS, int recoMethod) : 
00013   recoMethod_(recoMethod),
00014   correctForTimeslew_(correctForTimeslew),
00015   correctForPulse_(correctForPulse) {
00016 }
00017 
00018 ZdcSimpleRecAlgo::ZdcSimpleRecAlgo(int recoMethod) : 
00019   recoMethod_(recoMethod),
00020   correctForTimeslew_(false) {
00021 }
00022 void ZdcSimpleRecAlgo::initPulseCorr(int toadd) {
00023   if (correctForPulse_) {    
00024     pulseCorr_=std::auto_ptr<HcalPulseContainmentCorrection>(new HcalPulseContainmentCorrection(toadd,phaseNS_,MaximumFractionalError));
00025   }
00026 }
00027 static float timeshift_ns_zdc(float wpksamp);
00028 
00029 namespace ZdcSimpleRecAlgoImpl {
00030   template<class Digi, class RecHit>
00031   inline RecHit reco1(const Digi& digi, const HcalCoder& coder, const HcalCalibrations& calibs, 
00032                       const std::vector<unsigned int>& myNoiseTS, const std::vector<unsigned int>& mySignalTS, bool slewCorrect, const HcalPulseContainmentCorrection* corr, HcalTimeSlew::BiasSetting slewFlavor) {
00033     CaloSamples tool;
00034     coder.adc2fC(digi,tool);
00035     int ifirst = mySignalTS[0];
00036     int n = mySignalTS.size();
00037     double ampl=0; int maxI = -1; double maxA = -1e10; double ta=0;
00038     double fc_ampl=0;
00039     for (int i=ifirst; i<tool.size() && i<n+ifirst; i++) {
00040       int capid=digi[i].capid();
00041       ta = (tool[i]-calibs.pedestal(capid)); // pedestal subtraction
00042       fc_ampl+=ta; 
00043       ta*= calibs.respcorrgain(capid) ; // fC --> GeV
00044       ampl+=ta;
00045       if(ta>maxA){
00046         maxA=ta;
00047         maxI=i;
00048       }
00049     }
00050     
00051     double time=-9999;
00053     if(maxI==0 || maxI==(tool.size()-1)) {      
00054       LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reconstruct :" 
00055                                                << " Invalid max amplitude position, " 
00056                                                << " max Amplitude: "<< maxI
00057                                                << " first: "<<ifirst
00058                                                << " last: "<<(tool.size()-1)
00059                                                << std::endl;
00060     } else {
00061       maxA=fabs(maxA);
00062       int capid=digi[maxI-1].capid();
00063       double t0 = fabs((tool[maxI-1]-calibs.pedestal(capid))*calibs.respcorrgain(capid) );
00064       capid=digi[maxI+1].capid();
00065       double t2 = fabs((tool[maxI+1]-calibs.pedestal(capid))*calibs.respcorrgain(capid) );    
00066       double wpksamp = (t0 + maxA + t2);
00067       if (wpksamp!=0) wpksamp=(maxA + 2.0*t2) / wpksamp; 
00068       time = (maxI - digi.presamples())*25.0 + timeshift_ns_zdc(wpksamp);
00069       if (corr!=0) {
00070         // Apply phase-based amplitude correction:
00071         ampl *= corr->getCorrection(fc_ampl);
00072       }
00073     
00074       if (slewCorrect) time-=HcalTimeSlew::delay(std::max(1.0,fc_ampl),slewFlavor);
00075     }
00076 
00077     time=time-calibs.timecorr(); // time calibration
00078     return RecHit(digi.id(),ampl,time);    
00079   }
00080 }
00081 
00082 namespace ZdcSimpleRecAlgoImpl {
00083   template<class Digi, class RecHit>
00084   inline RecHit reco2(const Digi& digi, const HcalCoder& coder, const HcalCalibrations& calibs, 
00085                      const std::vector<unsigned int>& myNoiseTS, const std::vector<unsigned int>& mySignalTS, bool slewCorrect, const HcalPulseContainmentCorrection* corr, HcalTimeSlew::BiasSetting slewFlavor) {
00086     CaloSamples tool;
00087     coder.adc2fC(digi,tool);
00088     // Reads noiseTS and signalTS from database
00089     int ifirst = mySignalTS[0];
00090 //    int n = mySignalTS.size();
00091     double ampl=0; int maxI = -1; double maxA = -1e10; double ta=0;
00092     double Allnoise = 0; 
00093     int noiseslices = 0;
00094     int CurrentTS = 0;
00095     double noise = 0;
00096     double fc_ampl=0;
00097     
00098     for(unsigned int iv = 0; iv<myNoiseTS.size(); ++iv)
00099     {
00100       CurrentTS = myNoiseTS[iv];
00101       Allnoise += tool[CurrentTS];
00102       noiseslices++;
00103     }
00104     if(noiseslices != 0) {
00105       noise = (Allnoise)/double(noiseslices);
00106     } else {
00107       noise = 0;
00108     }
00109  // factor to multiply by noise to make 0 or 1 to handle negative noise situations
00110     double noisefactor=1.;
00111     for(unsigned int ivs = 0; ivs<mySignalTS.size(); ++ivs)
00112     {
00113       if(noise<0){
00114       // flag hit as having negative noise, and don't subtract anything, because
00115       // it will falsely increase the energy
00116          noisefactor=0.;
00117       } 
00118       CurrentTS = mySignalTS[ivs];
00119       int capid=digi[CurrentTS].capid();
00120       if(noise<0){
00121       // flag hit as having negative noise, and don't subtract anything, because
00122       // it will falsely increase the energy
00123          noisefactor=0.;
00124       } 
00125       ta = tool[CurrentTS]-noisefactor*noise;
00126       fc_ampl+=ta; 
00127       ta*= calibs.respcorrgain(capid) ; // fC --> GeV
00128       ampl+=ta;
00129       if(ta>maxA){
00130              maxA=ta;
00131              maxI=CurrentTS;
00132           }
00133     }  
00134     
00135 //    if(ta<0){
00136 //      // flag hits that have negative energy
00137 //    }
00138 
00139     double time=-9999;
00141     if(maxI==0 || maxI==(tool.size()-1)) {      
00142       LogDebug("HCAL Pulse") << "ZdcSimpleRecAlgo::reco2 :" 
00143                                                << " Invalid max amplitude position, " 
00144                                                << " max Amplitude: "<< maxI
00145                                                << " first: "<<ifirst
00146                                                << " last: "<<(tool.size()-1)
00147                                                << std::endl;
00148     } else {
00149       int capid=digi[maxI-1].capid();
00150       double Energy0 = ((tool[maxI-1])*calibs.respcorrgain(capid) );
00151 // if any of the energies used in the weight are negative, make them 0 instead
00152 // these are actually QIE values, not energy
00153       if(Energy0<0){Energy0=0.;}
00154       capid=digi[maxI].capid();
00155       double Energy1 = ((tool[maxI])*calibs.respcorrgain(capid) ) ;
00156       if(Energy1<0){Energy1=0.;}
00157       capid=digi[maxI+1].capid();
00158       double Energy2 = ((tool[maxI+1])*calibs.respcorrgain(capid) );
00159       if(Energy2<0){Energy2=0.;}
00160 //
00161       double TSWeightEnergy = ((maxI-1)*Energy0 + maxI*Energy1 + (maxI+1)*Energy2);
00162       double EnergySum=Energy0+Energy1+Energy2;
00163       double AvgTSPos=0.;
00164       if (EnergySum!=0) AvgTSPos=TSWeightEnergy/ EnergySum; 
00165 // If time is zero, set it to the "nonsensical" -99
00166 // Time should be between 75ns and 175ns (Timeslices 3-7)
00167       if(AvgTSPos==0){
00168          time=-99;
00169       } else {
00170          time = (AvgTSPos*25.0);
00171       }
00172       if (corr!=0) {
00173         // Apply phase-based amplitude correction:
00174                ampl *= corr->getCorrection(fc_ampl);
00175       }
00176     }
00177     return RecHit(digi.id(),ampl,time);    
00178   }
00179 }
00180 
00181 ZDCRecHit ZdcSimpleRecAlgo::reconstruct(const ZDCDataFrame& digi, const std::vector<unsigned int>& myNoiseTS, const std::vector<unsigned int>& mySignalTS,const HcalCoder& coder, const HcalCalibrations& calibs) const {
00182  
00183   if(recoMethod_ == 1)
00184    return ZdcSimpleRecAlgoImpl::reco1<ZDCDataFrame,ZDCRecHit>(digi,coder,calibs,
00185                                                               myNoiseTS,mySignalTS,false,
00186                                                               0,
00187                                                               HcalTimeSlew::Fast);
00188   if(recoMethod_ == 2)
00189    return ZdcSimpleRecAlgoImpl::reco2<ZDCDataFrame,ZDCRecHit>(digi,coder,calibs,
00190                                                               myNoiseTS,mySignalTS,false,
00191                                                               0,HcalTimeSlew::Fast);
00192 
00193      edm::LogError("ZDCSimpleRecAlgoImpl::reconstruct, recoMethod was not declared");
00194      throw cms::Exception("ZDCSimpleRecoAlgos::exiting process");
00195    
00196 }
00197 
00198 static const float wpksamp0_zdc = 0.500053;
00199 static const float scale_zdc    = 0.999683;
00200 static const int   num_bins_zdc = 100;
00201 
00202 static const float actual_ns_zdc[num_bins_zdc] = {
00203  0.00250, // 0.000-0.010
00204  0.08000, // 0.010-0.020
00205  0.16000, // 0.020-0.030
00206  0.23750, // 0.030-0.040
00207  0.31750, // 0.040-0.050
00208  0.39500, // 0.050-0.060
00209  0.47500, // 0.060-0.070
00210  0.55500, // 0.070-0.080
00211  0.63000, // 0.080-0.090
00212  0.70000, // 0.090-0.100
00213  0.77000, // 0.100-0.110
00214  0.84000, // 0.110-0.120
00215  0.91000, // 0.120-0.130
00216  0.98000, // 0.130-0.140
00217  1.05000, // 0.140-0.150
00218  1.12000, // 0.150-0.160
00219  1.19000, // 0.160-0.170
00220  1.26000, // 0.170-0.180
00221  1.33000, // 0.180-0.190
00222  1.40000, // 0.190-0.200
00223  1.47000, // 0.200-0.210
00224  1.54000, // 0.210-0.220
00225  1.61000, // 0.220-0.230
00226  1.68000, // 0.230-0.240
00227  1.75000, // 0.240-0.250
00228  1.82000, // 0.250-0.260
00229  1.89000, // 0.260-0.270
00230  1.96000, // 0.270-0.280
00231  2.03000, // 0.280-0.290
00232  2.10000, // 0.290-0.300
00233  2.17000, // 0.300-0.310
00234  2.24000, // 0.310-0.320
00235  2.31000, // 0.320-0.330
00236  2.38000, // 0.330-0.340
00237  2.45000, // 0.340-0.350
00238  2.52000, // 0.350-0.360
00239  2.59000, // 0.360-0.370
00240  2.68500, // 0.370-0.380
00241  2.79250, // 0.380-0.390
00242  2.90250, // 0.390-0.400
00243  3.01000, // 0.400-0.410
00244  3.11750, // 0.410-0.420
00245  3.22500, // 0.420-0.430
00246  3.33500, // 0.430-0.440
00247  3.44250, // 0.440-0.450
00248  3.55000, // 0.450-0.460
00249  3.73250, // 0.460-0.470
00250  4.02000, // 0.470-0.480
00251  4.30750, // 0.480-0.490
00252  4.59500, // 0.490-0.500
00253  6.97500, // 0.500-0.510
00254 10.98750, // 0.510-0.520
00255 13.03750, // 0.520-0.530
00256 14.39250, // 0.530-0.540
00257 15.39500, // 0.540-0.550
00258 16.18250, // 0.550-0.560
00259 16.85250, // 0.560-0.570
00260 17.42750, // 0.570-0.580
00261 17.91500, // 0.580-0.590
00262 18.36250, // 0.590-0.600
00263 18.76500, // 0.600-0.610
00264 19.11250, // 0.610-0.620
00265 19.46000, // 0.620-0.630
00266 19.76500, // 0.630-0.640
00267 20.03500, // 0.640-0.650
00268 20.30250, // 0.650-0.660
00269 20.57250, // 0.660-0.670
00270 20.79250, // 0.670-0.680
00271 21.00250, // 0.680-0.690
00272 21.21000, // 0.690-0.700
00273 21.42000, // 0.700-0.710
00274 21.62750, // 0.710-0.720
00275 21.79000, // 0.720-0.730
00276 21.95250, // 0.730-0.740
00277 22.11500, // 0.740-0.750
00278 22.27750, // 0.750-0.760
00279 22.44000, // 0.760-0.770
00280 22.60500, // 0.770-0.780
00281 22.73250, // 0.780-0.790
00282 22.86000, // 0.790-0.800
00283 22.98500, // 0.800-0.810
00284 23.11250, // 0.810-0.820
00285 23.23750, // 0.820-0.830
00286 23.36500, // 0.830-0.840
00287 23.49000, // 0.840-0.850
00288 23.61750, // 0.850-0.860
00289 23.71500, // 0.860-0.870
00290 23.81250, // 0.870-0.880
00291 23.91250, // 0.880-0.890
00292 24.01000, // 0.890-0.900
00293 24.10750, // 0.900-0.910
00294 24.20750, // 0.910-0.920
00295 24.30500, // 0.920-0.930
00296 24.40500, // 0.930-0.940
00297 24.50250, // 0.940-0.950
00298 24.60000, // 0.950-0.960
00299 24.68250, // 0.960-0.970
00300 24.76250, // 0.970-0.980
00301 24.84000, // 0.980-0.990
00302 24.92000  // 0.990-1.000
00303 };
00304 
00305 float timeshift_ns_zdc(float wpksamp) {
00306   float flx = (num_bins_zdc*(wpksamp - wpksamp0_zdc)/scale_zdc);
00307   int index = (int)flx;
00308   float yval;
00309   
00310   if      (index <    0)        return actual_ns_zdc[0];
00311   else if (index >= num_bins_zdc-1) return actual_ns_zdc[num_bins_zdc-1];
00312 
00313   // else interpolate:
00314   float y1       = actual_ns_zdc[index];
00315   float y2       = actual_ns_zdc[index+1];
00316 
00317   // float delta_x  = 1/(float)num_bins_zdc;
00318   // yval = y1 + (y2-y1)*(flx-(float)index)/delta_x;
00319 
00320   yval = y1 + (y2-y1)*(flx-(float)index);
00321   return yval;
00322 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoLocalCalo/HcalRecAlgos/src/ZdcSimpleRecAlgo.cc
