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