Reference Manual

Go to the documentation of this file.
00001 #include "HcalHitReconstructor.h"
00002 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00003 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00004 #include "DataFormats/Common/interface/EDCollection.h"
00005 #include "DataFormats/Common/interface/Handle.h"
00006 #include "FWCore/Framework/interface/Selector.h"
00007 #include "FWCore/Framework/interface/ESHandle.h"
00008 #include "FWCore/Framework/interface/EventSetup.h"
00009 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
00010 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
00011 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
00012 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
00013 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
00014 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
00015 #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
00016 
00017 #include <iostream>
00018 #include <fstream>
00019 
00020 
00021 /*  Hcal Hit reconstructor allows for CaloRecHits with status words */
00022 
00023 HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf):
00024   reco_(conf.getParameter<bool>("correctForTimeslew"),
00025         conf.getParameter<bool>("correctForPhaseContainment"),
00026         conf.getParameter<double>("correctionPhaseNS")),
00027   det_(DetId::Hcal),
00028   inputLabel_(conf.getParameter<edm::InputTag>("digiLabel")),
00029   correctTiming_(conf.getParameter<bool>("correctTiming")),
00030   setNoiseFlags_(conf.getParameter<bool>("setNoiseFlags")),
00031   setHSCPFlags_(conf.getParameter<bool>("setHSCPFlags")),
00032   setSaturationFlags_(conf.getParameter<bool>("setSaturationFlags")),
00033   setTimingTrustFlags_(conf.getParameter<bool>("setTimingTrustFlags")),
00034   setPulseShapeFlags_(conf.getParameter<bool>("setPulseShapeFlags")),
00035   dropZSmarkedPassed_(conf.getParameter<bool>("dropZSmarkedPassed")),
00036   firstAuxTS_(conf.getParameter<int>("firstAuxTS")),
00037   firstSample_(conf.getParameter<int>("firstSample")),
00038   samplesToAdd_(conf.getParameter<int>("samplesToAdd")),
00039   tsFromDB_(conf.getParameter<bool>("tsFromDB")),
00040   useLeakCorrection_( conf.getParameter<bool>("useLeakCorrection")) 
00041 {
00042 
00043   std::string subd=conf.getParameter<std::string>("Subdetector");
00044   //Set all FlagSetters to 0
00045   /* Important to do this!  Otherwise, if the setters are turned off,
00046      the "if (XSetter_) delete XSetter_;" commands can crash
00047   */
00048 
00049   recoParamsFromDB_ = conf.getParameter<bool>("recoParamsFromDB");
00050   //  recoParamsFromDB_ = false ; //  trun off for now.
00051 
00052   // std::cout<<"  HcalHitReconstructor   recoParamsFromDB_ "<<recoParamsFromDB_<<std::endl;
00053 
00054   hbheFlagSetter_             = 0;
00055   hbheHSCPFlagSetter_         = 0;
00056   hbhePulseShapeFlagSetter_   = 0;
00057   hbheTimingShapedFlagSetter_ = 0;
00058   hfdigibit_                  = 0;
00059 
00060   hfS9S1_                     = 0;
00061   hfS8S1_                     = 0;
00062   hfPET_                      = 0;
00063   saturationFlagSetter_       = 0;
00064   HFTimingTrustFlagSetter_    = 0;
00065   digiTimeFromDB_             = false; // only need for HF
00066   
00067   if (setSaturationFlags_)
00068     {
00069       const edm::ParameterSet& pssat      = conf.getParameter<edm::ParameterSet>("saturationParameters");
00070       saturationFlagSetter_ = new HcalADCSaturationFlag(pssat.getParameter<int>("maxADCvalue"));
00071     }
00072 
00073   if (!strcasecmp(subd.c_str(),"HBHE")) {
00074     subdet_=HcalBarrel;
00075     bool timingShapedCutsFlags = conf.getParameter<bool>("setTimingShapedCutsFlags");
00076     if (timingShapedCutsFlags)
00077       {
00078         const edm::ParameterSet& psTshaped = conf.getParameter<edm::ParameterSet>("timingshapedcutsParameters");
00079         hbheTimingShapedFlagSetter_ = new HBHETimingShapedFlagSetter(psTshaped.getParameter<std::vector<double> >("tfilterEnvelope"),
00080                                                                      psTshaped.getParameter<bool>("ignorelowest"),
00081                                                                      psTshaped.getParameter<bool>("ignorehighest"),
00082                                                                      psTshaped.getParameter<double>("win_offset"),
00083                                                                      psTshaped.getParameter<double>("win_gain"));
00084       }
00085       
00086     if (setNoiseFlags_)
00087       {
00088         const edm::ParameterSet& psdigi    =conf.getParameter<edm::ParameterSet>("flagParameters");
00089         hbheFlagSetter_=new HBHEStatusBitSetter(psdigi.getParameter<double>("nominalPedestal"),
00090                                                 psdigi.getParameter<double>("hitEnergyMinimum"),
00091                                                 psdigi.getParameter<int>("hitMultiplicityThreshold"),
00092                                                 psdigi.getParameter<std::vector<edm::ParameterSet> >("pulseShapeParameterSets")
00093          );
00094       } // if (setNoiseFlags_)
00095     if (setHSCPFlags_)
00096       {
00097         const edm::ParameterSet& psHSCP = conf.getParameter<edm::ParameterSet>("hscpParameters");
00098         hbheHSCPFlagSetter_ = new HBHETimeProfileStatusBitSetter(psHSCP.getParameter<double>("r1Min"),
00099                                                                  psHSCP.getParameter<double>("r1Max"),
00100                                                                  psHSCP.getParameter<double>("r2Min"),
00101                                                                  psHSCP.getParameter<double>("r2Max"),
00102                                                                  psHSCP.getParameter<double>("fracLeaderMin"),
00103                                                                  psHSCP.getParameter<double>("fracLeaderMax"),
00104                                                                  psHSCP.getParameter<double>("slopeMin"),
00105                                                                  psHSCP.getParameter<double>("slopeMax"),
00106                                                                  psHSCP.getParameter<double>("outerMin"),
00107                                                                  psHSCP.getParameter<double>("outerMax"),
00108                                                                  psHSCP.getParameter<double>("TimingEnergyThreshold"));
00109       } // if (setHSCPFlags_) 
00110     if (setPulseShapeFlags_)
00111       {
00112         const edm::ParameterSet &psPulseShape = conf.getParameter<edm::ParameterSet>("pulseShapeParameters");
00113         hbhePulseShapeFlagSetter_ = new HBHEPulseShapeFlagSetter(
00114                                                                  psPulseShape.getParameter<double>("MinimumChargeThreshold"),
00115                                                                  psPulseShape.getParameter<double>("TS4TS5ChargeThreshold"),
00116                                                                  psPulseShape.getParameter<unsigned int>("TrianglePeakTS"),
00117                                                                  psPulseShape.getParameter<std::vector<double> >("LinearThreshold"),
00118                                                                  psPulseShape.getParameter<std::vector<double> >("LinearCut"),
00119                                                                  psPulseShape.getParameter<std::vector<double> >("RMS8MaxThreshold"),
00120                                                                  psPulseShape.getParameter<std::vector<double> >("RMS8MaxCut"),
00121                                                                  psPulseShape.getParameter<std::vector<double> >("LeftSlopeThreshold"),
00122                                                                  psPulseShape.getParameter<std::vector<double> >("LeftSlopeCut"),
00123                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeThreshold"),
00124                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeCut"),
00125                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeSmallThreshold"),
00126                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeSmallCut"),
00127                                                                  psPulseShape.getParameter<std::vector<double> >("TS4TS5LowerThreshold"),
00128                                                                  psPulseShape.getParameter<std::vector<double> >("TS4TS5LowerCut"),
00129                                                                  psPulseShape.getParameter<std::vector<double> >("TS4TS5UpperThreshold"),
00130                                                                  psPulseShape.getParameter<std::vector<double> >("TS4TS5UpperCut"),
00131                                                                  psPulseShape.getParameter<bool>("UseDualFit"),
00132                          psPulseShape.getParameter<bool>("TriangleIgnoreSlow"));
00133       }  // if (setPulseShapeFlags_)
00134 
00135     produces<HBHERecHitCollection>();
00136   } else if (!strcasecmp(subd.c_str(),"HO")) {
00137     subdet_=HcalOuter;
00138     produces<HORecHitCollection>();
00139   } else if (!strcasecmp(subd.c_str(),"HF")) {
00140     subdet_=HcalForward;
00141     digiTimeFromDB_=conf.getParameter<bool>("digiTimeFromDB");
00142 
00143     if (setTimingTrustFlags_) {
00144       
00145       const edm::ParameterSet& pstrust      = conf.getParameter<edm::ParameterSet>("hfTimingTrustParameters");
00146       HFTimingTrustFlagSetter_=new HFTimingTrustFlag(pstrust.getParameter<int>("hfTimingTrustLevel1"),
00147                                                      pstrust.getParameter<int>("hfTimingTrustLevel2"));
00148     }
00149 
00150     if (setNoiseFlags_)
00151       {
00152         const edm::ParameterSet& psdigi    =conf.getParameter<edm::ParameterSet>("digistat");
00153         const edm::ParameterSet& psTimeWin =conf.getParameter<edm::ParameterSet>("HFInWindowStat");
00154         hfdigibit_=new HcalHFStatusBitFromDigis(psdigi,psTimeWin);
00155 
00156         const edm::ParameterSet& psS9S1   = conf.getParameter<edm::ParameterSet>("S9S1stat");
00157         hfS9S1_   = new HcalHF_S9S1algorithm(psS9S1.getParameter<std::vector<double> >("short_optimumSlope"),
00158                                              psS9S1.getParameter<std::vector<double> >("shortEnergyParams"),
00159                                              psS9S1.getParameter<std::vector<double> >("shortETParams"),
00160                                              psS9S1.getParameter<std::vector<double> >("long_optimumSlope"),
00161                                              psS9S1.getParameter<std::vector<double> >("longEnergyParams"),
00162                                              psS9S1.getParameter<std::vector<double> >("longETParams"),
00163                                              psS9S1.getParameter<int>("HcalAcceptSeverityLevel"),
00164                                              psS9S1.getParameter<bool>("isS8S1")
00165                                              );
00166 
00167         const edm::ParameterSet& psS8S1   = conf.getParameter<edm::ParameterSet>("S8S1stat");
00168         hfS8S1_   = new HcalHF_S9S1algorithm(psS8S1.getParameter<std::vector<double> >("short_optimumSlope"),
00169                                              psS8S1.getParameter<std::vector<double> >("shortEnergyParams"),
00170                                              psS8S1.getParameter<std::vector<double> >("shortETParams"),
00171                                              psS8S1.getParameter<std::vector<double> >("long_optimumSlope"),
00172                                              psS8S1.getParameter<std::vector<double> >("longEnergyParams"),
00173                                              psS8S1.getParameter<std::vector<double> >("longETParams"),
00174                                              psS8S1.getParameter<int>("HcalAcceptSeverityLevel"),
00175                                              psS8S1.getParameter<bool>("isS8S1")
00176                                              );
00177 
00178         const edm::ParameterSet& psPET    = conf.getParameter<edm::ParameterSet>("PETstat");
00179         hfPET_    = new HcalHF_PETalgorithm(psPET.getParameter<std::vector<double> >("short_R"),
00180                                             psPET.getParameter<std::vector<double> >("shortEnergyParams"),
00181                                             psPET.getParameter<std::vector<double> >("shortETParams"),
00182                                             psPET.getParameter<std::vector<double> >("long_R"),
00183                                             psPET.getParameter<std::vector<double> >("longEnergyParams"),
00184                                             psPET.getParameter<std::vector<double> >("longETParams"),
00185                                             psPET.getParameter<int>("HcalAcceptSeverityLevel"),
00186                                             psPET.getParameter<std::vector<double> >("short_R_29"),
00187                                             psPET.getParameter<std::vector<double> >("long_R_29")
00188                                             );
00189       }
00190     produces<HFRecHitCollection>();
00191   } else if (!strcasecmp(subd.c_str(),"ZDC")) {
00192     det_=DetId::Calo;
00193     subdet_=HcalZDCDetId::SubdetectorId;
00194     produces<ZDCRecHitCollection>();
00195   } else if (!strcasecmp(subd.c_str(),"CALIB")) {
00196     subdet_=HcalOther;
00197     subdetOther_=HcalCalibration;
00198     produces<HcalCalibRecHitCollection>();
00199   } else {
00200      std::cout << "HcalHitReconstructor is not associated with a specific subdetector!" << std::endl;
00201   }       
00202   
00203 }
00204 
00205 HcalHitReconstructor::~HcalHitReconstructor() {
00206   if (hbheFlagSetter_)        delete hbheFlagSetter_;
00207   if (hfdigibit_)             delete hfdigibit_;
00208   if (hbheHSCPFlagSetter_)    delete hbheHSCPFlagSetter_;
00209   if (hbhePulseShapeFlagSetter_) delete hbhePulseShapeFlagSetter_;
00210   if (hfS9S1_)                delete hfS9S1_;
00211   if (hfPET_)                 delete hfPET_;
00212 }
00213 
00214 void HcalHitReconstructor::beginRun(edm::Run&r, edm::EventSetup const & es){
00215 
00216   if ( tsFromDB_== true || recoParamsFromDB_ == true )
00217     {
00218       edm::ESHandle<HcalRecoParams> p;
00219       es.get<HcalRecoParamsRcd>().get(p);
00220       paramTS = new HcalRecoParams(*p.product());
00221 
00222       // std::cout<<" skdump in HcalHitReconstructor::beginRun   dupm RecoParams "<<std::endl;
00223       // std::ofstream skfile("skdumpRecoParamsNewFormat.txt");
00224       // HcalDbASCIIIO::dumpObject(skfile, (*paramTS) );
00225     }
00226 
00227   if (digiTimeFromDB_==true)
00228     {
00229       edm::ESHandle<HcalFlagHFDigiTimeParams> p;
00230       es.get<HcalFlagHFDigiTimeParamsRcd>().get(p);
00231       HFDigiTimeParams = new HcalFlagHFDigiTimeParams(*p.product());
00232     }
00233   reco_.beginRun(es);
00234 }
00235 
00236 void HcalHitReconstructor::endRun(edm::Run&r, edm::EventSetup const & es){
00237   if (tsFromDB_==true)
00238     {
00239       if (paramTS) delete paramTS;
00240     }
00241   if (digiTimeFromDB_==true)
00242     {
00243       if (HFDigiTimeParams) delete HFDigiTimeParams;
00244     }
00245   reco_.endRun();
00246 }
00247 
00248 void HcalHitReconstructor::produce(edm::Event& e, const edm::EventSetup& eventSetup)
00249 {
00250 
00251   // get conditions
00252   edm::ESHandle<HcalDbService> conditions;
00253   eventSetup.get<HcalDbRecord>().get(conditions);
00254   const HcalQIEShape* shape = conditions->getHcalShape (); // this one is generic
00255   // HACK related to HB- corrections
00256   if(e.isRealData()) reco_.setForData();    
00257   if(useLeakCorrection_) reco_.setLeakCorrection();
00258 
00259   edm::ESHandle<HcalChannelQuality> p;
00260   eventSetup.get<HcalChannelQualityRcd>().get(p);
00261   HcalChannelQuality* myqual = new HcalChannelQuality(*p.product());
00262 
00263   edm::ESHandle<HcalSeverityLevelComputer> mycomputer;
00264   eventSetup.get<HcalSeverityLevelComputerRcd>().get(mycomputer);
00265   const HcalSeverityLevelComputer* mySeverity = mycomputer.product();
00266   
00267   if (det_==DetId::Hcal) {
00268 
00269     // HBHE -------------------------------------------------------------------
00270     if (subdet_==HcalBarrel || subdet_==HcalEndcap) {
00271       edm::Handle<HBHEDigiCollection> digi;
00272       
00273       e.getByLabel(inputLabel_,digi);
00274       
00275       // create empty output
00276       std::auto_ptr<HBHERecHitCollection> rec(new HBHERecHitCollection);
00277       rec->reserve(digi->size());
00278       // run the algorithm
00279       if (setNoiseFlags_) hbheFlagSetter_->Clear();
00280       HBHEDigiCollection::const_iterator i;
00281       std::vector<HBHEDataFrame> HBDigis;
00282       std::vector<int> RecHitIndex;
00283 
00284       // Vote on majority TS0 CapId
00285       int favorite_capid = 0; 
00286       if (correctTiming_) {
00287         long capid_votes[4] = {0,0,0,0};
00288         for (i=digi->begin(); i!=digi->end(); i++) {
00289           capid_votes[(*i)[0].capid()]++;
00290         }
00291         for (int k = 0; k < 4; k++)
00292           if (capid_votes[k] > capid_votes[favorite_capid])
00293             favorite_capid = k;
00294       }
00295 
00296       int first = firstSample_;
00297       int toadd = samplesToAdd_;
00298 
00299       for (i=digi->begin(); i!=digi->end(); i++) {
00300         HcalDetId cell = i->id();
00301         DetId detcell=(DetId)cell;
00302 
00303         if(tsFromDB_ || recoParamsFromDB_) {
00304           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00305           firstSample_ = param_ts->firstSample();
00306           samplesToAdd_ = param_ts->samplesToAdd();
00307           if(recoParamsFromDB_) {
00308              bool correctForTimeslew=param_ts->correctForTimeslew();
00309              bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
00310              float phaseNS=param_ts->correctionPhaseNS();
00311              useLeakCorrection_= param_ts->useLeakCorrection();
00312              correctTiming_ = param_ts->correctTiming();
00313              firstAuxTS_ = param_ts->firstAuxTS();
00314              int pileupCleaningID = param_ts->pileupCleaningID();
00315 
00316              /*      
00317              int sub     = cell.subdet();
00318              int depth   = cell.depth();
00319              int inteta  = cell.ieta();
00320              int intphi  = cell.iphi();
00321 
00322              std::cout << "HcalHitReconstructor::produce  cell:" 
00323                        << " sub, ieta, iphi, depth = " 
00324                        << sub << "  " << inteta << "  " << intphi 
00325                        << "  " << depth << std::endl
00326                        << "    first, toadd = " << firstSample_ << ", "
00327                        << samplesToAdd_ << std::endl
00328                        << "    correctForTimeslew " << correctForTimeslew
00329                        << std::endl
00330                        << "    correctForPhaseContainment " 
00331                        <<  correctForPhaseContainment << std::endl
00332                        << "    phaseNS " <<  phaseNS << std::endl
00333                        << "    useLeakCorrection  " << useLeakCorrection_ 
00334                        << std::endl 
00335                        << "    correctTiming " << correctTiming_ << std::endl
00336                        << "    firstAuxTS " << firstAuxTS_  << std::endl
00337                        << "    pileupCleaningID "  << pileupCleaningID
00338                        << std::endl;
00339              */
00340 
00341              reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
00342           }
00343         }
00344 
00345         first = firstSample_;
00346         toadd = samplesToAdd_;
00347 
00348         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00349         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00350         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00351         if (dropZSmarkedPassed_)
00352           if (i->zsMarkAndPass()) continue;
00353 
00354         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00355         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00356         HcalCoderDb coder (*channelCoder, *shape);
00357 
00358         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00359 
00360         // Set auxiliary flag
00361         int auxflag=0;
00362         int fTS = firstAuxTS_;
00363         if (fTS<0) fTS=0; // silly protection against time slice <0
00364         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00365           auxflag+=(i->sample(xx).adc())<<(7*(xx-fTS)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00366         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00367         auxflag+=((i->sample(fTS).capid())<<28);
00368         (rec->back()).setAux(auxflag);
00369 
00370         (rec->back()).setFlags(0);  // this sets all flag bits to 0
00371         // Set presample flag
00372         if (fTS>0)
00373           (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);
00374 
00375         if (hbheTimingShapedFlagSetter_!=0)
00376           hbheTimingShapedFlagSetter_->SetTimingShapedFlags(rec->back());
00377         if (setNoiseFlags_)
00378           hbheFlagSetter_->SetFlagsFromDigi(rec->back(),*i,coder,calibrations,first,toadd);
00379         if (setPulseShapeFlags_ == true)
00380           hbhePulseShapeFlagSetter_->SetPulseShapeFlags(rec->back(), *i, coder, calibrations);
00381         if (setSaturationFlags_)
00382           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00383         if (correctTiming_)
00384           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00385         if (setHSCPFlags_ && i->id().ietaAbs()<16)
00386           {
00387             double DigiEnergy=0;
00388             for(int j=0; j!=i->size(); DigiEnergy += i->sample(j++).nominal_fC());
00389             if(DigiEnergy > hbheHSCPFlagSetter_->EnergyThreshold())
00390               {
00391                 HBDigis.push_back(*i);
00392                 RecHitIndex.push_back(rec->size()-1);
00393               }
00394             
00395           } // if (set HSCPFlags_ && |ieta|<16)
00396       } // loop over HBHE digis
00397 
00398 
00399       if (setNoiseFlags_) hbheFlagSetter_->SetFlagsFromRecHits(*rec);
00400       if (setHSCPFlags_)  hbheHSCPFlagSetter_->hbheSetTimeFlagsFromDigi(rec.get(), HBDigis, RecHitIndex);
00401       // return result
00402       e.put(rec);
00403 
00404       //  HO ------------------------------------------------------------------
00405     } else if (subdet_==HcalOuter) {
00406       edm::Handle<HODigiCollection> digi;
00407       e.getByLabel(inputLabel_,digi);
00408       
00409       // create empty output
00410       std::auto_ptr<HORecHitCollection> rec(new HORecHitCollection);
00411       rec->reserve(digi->size());
00412       // run the algorithm
00413       HODigiCollection::const_iterator i;
00414 
00415       // Vote on majority TS0 CapId
00416       int favorite_capid = 0; 
00417       if (correctTiming_) {
00418         long capid_votes[4] = {0,0,0,0};
00419         for (i=digi->begin(); i!=digi->end(); i++) {
00420           capid_votes[(*i)[0].capid()]++;
00421         }
00422         for (int k = 0; k < 4; k++)
00423           if (capid_votes[k] > capid_votes[favorite_capid])
00424             favorite_capid = k;
00425       }
00426 
00427       int first = firstSample_;
00428       int toadd = samplesToAdd_;
00429 
00430       for (i=digi->begin(); i!=digi->end(); i++) {
00431         HcalDetId cell = i->id();
00432         DetId detcell=(DetId)cell;
00433         // firstSample & samplesToAdd
00434         if(tsFromDB_ || recoParamsFromDB_) {
00435           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00436           firstSample_ = param_ts->firstSample();
00437           samplesToAdd_ = param_ts->samplesToAdd();
00438           if(recoParamsFromDB_) {
00439              bool correctForTimeslew=param_ts->correctForTimeslew();
00440              bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
00441              float phaseNS=param_ts->correctionPhaseNS();
00442              useLeakCorrection_= param_ts->useLeakCorrection();
00443              correctTiming_ = param_ts->correctTiming();
00444              firstAuxTS_ = param_ts->firstAuxTS();
00445              int pileupCleaningID = param_ts->pileupCleaningID();
00446              reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
00447           }
00448         }
00449 
00450         first = firstSample_;
00451         toadd = samplesToAdd_;
00452 
00453         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00454         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00455         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00456         if (dropZSmarkedPassed_)
00457           if (i->zsMarkAndPass()) continue;
00458 
00459         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00460         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00461         HcalCoderDb coder (*channelCoder, *shape);
00462 
00463         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00464 
00465         // Set auxiliary flag
00466         int auxflag=0;
00467         int fTS = firstAuxTS_;
00468         if (fTS<0) fTS=0; //silly protection against negative time slice values
00469         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00470           auxflag+=(i->sample(xx).adc())<<(7*(xx-fTS)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00471         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00472         auxflag+=((i->sample(fTS).capid())<<28);
00473         (rec->back()).setAux(auxflag);
00474 
00475         (rec->back()).setFlags(0);
00476         // Fill Presample ADC flag
00477         if (fTS>0)
00478           (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);
00479 
00480         if (setSaturationFlags_)
00481           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00482         if (correctTiming_)
00483           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00484       }
00485       // return result
00486       e.put(rec);    
00487 
00488       // HF -------------------------------------------------------------------
00489     } else if (subdet_==HcalForward) {
00490       edm::Handle<HFDigiCollection> digi;
00491       e.getByLabel(inputLabel_,digi);
00492 
00493 
00495       // create empty output
00496       std::auto_ptr<HFRecHitCollection> rec(new HFRecHitCollection);
00497       rec->reserve(digi->size());
00498       // run the algorithm
00499       HFDigiCollection::const_iterator i;
00500 
00501       // Vote on majority TS0 CapId
00502       int favorite_capid = 0; 
00503       if (correctTiming_) {
00504         long capid_votes[4] = {0,0,0,0};
00505         for (i=digi->begin(); i!=digi->end(); i++) {
00506           capid_votes[(*i)[0].capid()]++;
00507         }
00508         for (int k = 0; k < 4; k++)
00509           if (capid_votes[k] > capid_votes[favorite_capid])
00510             favorite_capid = k;
00511       }
00512 
00513       int first = firstSample_;
00514       int toadd = samplesToAdd_;
00515 
00516       for (i=digi->begin(); i!=digi->end(); i++) {
00517         HcalDetId cell = i->id();
00518         DetId detcell=(DetId)cell;
00519 
00520         if(tsFromDB_ || recoParamsFromDB_) {
00521           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00522           firstSample_ = param_ts->firstSample();
00523           samplesToAdd_ = param_ts->samplesToAdd();
00524           if(recoParamsFromDB_) {
00525              bool correctForTimeslew=param_ts->correctForTimeslew();
00526              bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
00527              float phaseNS=param_ts->correctionPhaseNS();
00528              useLeakCorrection_= param_ts->useLeakCorrection();
00529              correctTiming_ = param_ts->correctTiming();
00530              firstAuxTS_ = param_ts->firstAuxTS();
00531              int pileupCleaningID = param_ts->pileupCleaningID();
00532              reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
00533           }
00534         }
00535 
00536         first = firstSample_;
00537         toadd = samplesToAdd_;
00538 
00539         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00540         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00541         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00542         if (dropZSmarkedPassed_)
00543           if (i->zsMarkAndPass()) continue;
00544 
00545         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00546         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00547         HcalCoderDb coder (*channelCoder, *shape);
00548 
00549         // Set HFDigiTime flag values from digiTimeFromDB_
00550         if (digiTimeFromDB_==true && hfdigibit_!=0)
00551           {
00552             const HcalFlagHFDigiTimeParam* hfDTparam = HFDigiTimeParams->getValues(detcell.rawId());
00553             hfdigibit_->resetParamsFromDB(hfDTparam->HFdigiflagFirstSample(),
00554                                           hfDTparam->HFdigiflagSamplesToAdd(),
00555                                           hfDTparam->HFdigiflagExpectedPeak(),
00556                                           hfDTparam->HFdigiflagMinEThreshold(),
00557                                           hfDTparam->HFdigiflagCoefficients()
00558                                           );
00559           }
00560 
00561         //std::cout << "TOADDHF " << toadd << " " << first << " " << std::endl;
00562         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00563 
00564         // Set auxiliary flag
00565         int auxflag=0;
00566         int fTS = firstAuxTS_;
00567         if (fTS<0) fTS=0; // silly protection against negative time slice values
00568         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00569           auxflag+=(i->sample(xx).adc())<<(7*(xx-fTS)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00570         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00571         auxflag+=((i->sample(fTS).capid())<<28);
00572         (rec->back()).setAux(auxflag);
00573 
00574         // Clear flags
00575         (rec->back()).setFlags(0);
00576 
00577         // Fill Presample ADC flag
00578         if (fTS>0)
00579           (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);
00580 
00581         // This calls the code for setting the HF noise bit determined from digi shape
00582         if (setNoiseFlags_) 
00583           hfdigibit_->hfSetFlagFromDigi(rec->back(),*i,coder,calibrations);
00584         if (setSaturationFlags_)
00585           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00586         if (setTimingTrustFlags_)
00587           HFTimingTrustFlagSetter_->setHFTimingTrustFlag(rec->back(),*i);
00588         if (correctTiming_)
00589           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00590       } // for (i=digi->begin(); i!=digi->end(); i++) -- loop on all HF digis
00591 
00592       // The following flags require the full set of rechits
00593       // These need to be set consecutively, so an energy check should be the first 
00594       // test performed on these hits (to minimize the loop time)
00595       if (setNoiseFlags_) 
00596         {
00597           // Step 1:  Set PET flag  (short fibers of |ieta|==29)
00598           // Neighbor/partner channels that are flagged by Pulse Shape algorithm (HFDigiTime)
00599           // won't be considered in these calculations
00600           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00601             {
00602               int depth=i->id().depth();
00603               int ieta=i->id().ieta();
00604               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00605               if (depth==2 || abs(ieta)==29 ) 
00606                 hfPET_->HFSetFlagFromPET(*i,*rec,myqual,mySeverity);
00607             }
00608 
00609           // Step 2:  Set S8S1 flag (short fibers or |ieta|==29)
00610           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00611             {
00612               int depth=i->id().depth();
00613               int ieta=i->id().ieta();
00614               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00615               if (depth==2 || abs(ieta)==29 ) 
00616                 hfS8S1_->HFSetFlagFromS9S1(*i,*rec,myqual,mySeverity);
00617             }
00618 
00619           // Set 3:  Set S9S1 flag (long fibers)
00620           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00621             {
00622               int depth=i->id().depth();
00623               int ieta=i->id().ieta();
00624               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00625               if (depth==1 && abs(ieta)!=29 ) 
00626                 hfS9S1_->HFSetFlagFromS9S1(*i,*rec,myqual, mySeverity);
00627             }
00628         }
00629 
00630       // return result
00631       e.put(rec);     
00632     } else if (subdet_==HcalOther && subdetOther_==HcalCalibration) {
00633       edm::Handle<HcalCalibDigiCollection> digi;
00634       e.getByLabel(inputLabel_,digi);
00635       
00636       // create empty output
00637       std::auto_ptr<HcalCalibRecHitCollection> rec(new HcalCalibRecHitCollection);
00638       rec->reserve(digi->size());
00639       // run the algorithm
00640       int first = firstSample_;
00641       int toadd = samplesToAdd_;
00642 
00643       HcalCalibDigiCollection::const_iterator i;
00644       for (i=digi->begin(); i!=digi->end(); i++) {
00645         HcalCalibDetId cell = i->id();
00646         //      HcalDetId cellh = i->id();
00647         DetId detcell=(DetId)cell;
00648         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00649         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00650         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00651         if (dropZSmarkedPassed_)
00652           if (i->zsMarkAndPass()) continue;
00653 
00654         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00655         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00656         HcalCoderDb coder (*channelCoder, *shape);
00657 
00658         // firstSample & samplesToAdd
00659         if(tsFromDB_) {
00660           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00661           first = param_ts->firstSample();    
00662           toadd = param_ts->samplesToAdd();    
00663         }
00664         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00665 
00666         /*
00667           // Flag setting not available for calibration rechits
00668         // Set auxiliary flag
00669         int auxflag=0;
00670         int fTS = firstAuxTS_;
00671         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00672           auxflag+=(i->sample(xx).adc())<<(7*(xx-fTS)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00673         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00674         auxflag+=((i->sample(fTS).capid())<<28);
00675         (rec->back()).setAux(auxflag);
00676 
00677         (rec->back()).setFlags(0); // Not yet implemented for HcalCalibRecHit
00678         */
00679       }
00680       // return result
00681       e.put(rec);     
00682     }
00683   } 
00684   delete myqual;
00685 } // void HcalHitReconstructor::produce(...)
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc
