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