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 const&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 const&r, edm::EventSetup const & es){
00256   if (tsFromDB_==true)
00257     {
00258       delete paramTS; paramTS=0;
00259     }
00260   if (digiTimeFromDB_==true)
00261     {
00262       //DL delete HFDigiTimeParams; HFDigiTimeParams = 0;
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       for (i=digi->begin(); i!=digi->end(); i++) {
00322         HcalDetId cell = i->id();
00323         DetId detcell=(DetId)cell;
00324 
00325         if(tsFromDB_ || recoParamsFromDB_) {
00326           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00327           if(tsFromDB_) {
00328             firstSample_  = param_ts->firstSample();
00329             samplesToAdd_ = param_ts->samplesToAdd();
00330           }
00331           if(recoParamsFromDB_) {
00332              bool correctForTimeslew=param_ts->correctForTimeslew();
00333              bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
00334              float phaseNS=param_ts->correctionPhaseNS();
00335              useLeakCorrection_= param_ts->useLeakCorrection();
00336              correctTiming_ = param_ts->correctTiming();
00337              firstAuxTS_ = param_ts->firstAuxTS();
00338              int pileupCleaningID = param_ts->pileupCleaningID();
00339 
00340              /*      
00341              int sub     = cell.subdet();
00342              int depth   = cell.depth();
00343              int inteta  = cell.ieta();
00344              int intphi  = cell.iphi();
00345 
00346              std::cout << "HcalHitReconstructor::produce  cell:" 
00347                        << " sub, ieta, iphi, depth = " 
00348                        << sub << "  " << inteta << "  " << intphi 
00349                        << "  " << depth << std::endl
00350                        << "    first, toadd = " << firstSample_ << ", "
00351                        << samplesToAdd_ << std::endl
00352                        << "    correctForTimeslew " << correctForTimeslew
00353                        << std::endl
00354                        << "    correctForPhaseContainment " 
00355                        <<  correctForPhaseContainment << std::endl
00356                        << "    phaseNS " <<  phaseNS << std::endl
00357                        << "    useLeakCorrection  " << useLeakCorrection_ 
00358                        << std::endl 
00359                        << "    correctTiming " << correctTiming_ << std::endl
00360                        << "    firstAuxTS " << firstAuxTS_  << std::endl
00361                        << "    pileupCleaningID "  << pileupCleaningID
00362                        << std::endl;
00363              */
00364 
00365              reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
00366           }
00367         }
00368 
00369         int first = firstSample_;
00370         int toadd = samplesToAdd_;
00371 
00372         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00373         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00374         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00375         if (dropZSmarkedPassed_)
00376           if (i->zsMarkAndPass()) continue;
00377 
00378         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00379         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00380         const HcalQIEShape* shape = conditions->getHcalShape (channelCoder);
00381         HcalCoderDb coder (*channelCoder, *shape);
00382 
00383         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00384 
00385         // Set auxiliary flag
00386         int auxflag=0;
00387         int fTS = firstAuxTS_;
00388         if (fTS<0) fTS=0; // silly protection against time slice <0
00389         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00390           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
00391         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00392         auxflag+=((i->sample(fTS).capid())<<28);
00393         (rec->back()).setAux(auxflag);
00394 
00395         (rec->back()).setFlags(0);  // this sets all flag bits to 0
00396         // Set presample flag
00397         if (fTS>0)
00398           (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);
00399 
00400         if (hbheTimingShapedFlagSetter_!=0)
00401           hbheTimingShapedFlagSetter_->SetTimingShapedFlags(rec->back());
00402         if (setNoiseFlags_)
00403           hbheFlagSetter_->SetFlagsFromDigi(&(*topo),rec->back(),*i,coder,calibrations,first,toadd);
00404         if (setPulseShapeFlags_ == true)
00405           hbhePulseShapeFlagSetter_->SetPulseShapeFlags(rec->back(), *i, coder, calibrations);
00406         if (setSaturationFlags_)
00407           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00408         if (correctTiming_)
00409           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00410         if (setHSCPFlags_ && i->id().ietaAbs()<16)
00411           {
00412             double DigiEnergy=0;
00413             for(int j=0; j!=i->size(); DigiEnergy += i->sample(j++).nominal_fC());
00414             if(DigiEnergy > hbheHSCPFlagSetter_->EnergyThreshold())
00415               {
00416                 HBDigis.push_back(*i);
00417                 RecHitIndex.push_back(rec->size()-1);
00418               }
00419             
00420           } // if (set HSCPFlags_ && |ieta|<16)
00421       } // loop over HBHE digis
00422 
00423 
00424       if (setNoiseFlags_) hbheFlagSetter_->SetFlagsFromRecHits(&(*topo),*rec);
00425       if (setHSCPFlags_)  hbheHSCPFlagSetter_->hbheSetTimeFlagsFromDigi(rec.get(), HBDigis, RecHitIndex);
00426       // return result
00427       e.put(rec);
00428 
00429       //  HO ------------------------------------------------------------------
00430     } else if (subdet_==HcalOuter) {
00431       edm::Handle<HODigiCollection> digi;
00432       e.getByLabel(inputLabel_,digi);
00433       
00434       // create empty output
00435       std::auto_ptr<HORecHitCollection> rec(new HORecHitCollection);
00436       rec->reserve(digi->size());
00437       // run the algorithm
00438       HODigiCollection::const_iterator i;
00439 
00440       // Vote on majority TS0 CapId
00441       int favorite_capid = 0; 
00442       if (correctTiming_) {
00443         long capid_votes[4] = {0,0,0,0};
00444         for (i=digi->begin(); i!=digi->end(); i++) {
00445           capid_votes[(*i)[0].capid()]++;
00446         }
00447         for (int k = 0; k < 4; k++)
00448           if (capid_votes[k] > capid_votes[favorite_capid])
00449             favorite_capid = k;
00450       }
00451 
00452       for (i=digi->begin(); i!=digi->end(); i++) {
00453         HcalDetId cell = i->id();
00454         DetId detcell=(DetId)cell;
00455         // firstSample & samplesToAdd
00456         if(tsFromDB_ || recoParamsFromDB_) {
00457           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00458           if(tsFromDB_) {
00459             firstSample_  = param_ts->firstSample();
00460             samplesToAdd_ = param_ts->samplesToAdd();
00461           }
00462           if(recoParamsFromDB_) {
00463              bool correctForTimeslew=param_ts->correctForTimeslew();
00464              bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
00465              float phaseNS=param_ts->correctionPhaseNS();
00466              useLeakCorrection_= param_ts->useLeakCorrection();
00467              correctTiming_ = param_ts->correctTiming();
00468              firstAuxTS_ = param_ts->firstAuxTS();
00469              int pileupCleaningID = param_ts->pileupCleaningID();
00470              reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
00471           }
00472         }
00473 
00474         int first = firstSample_;
00475         int toadd = samplesToAdd_;
00476 
00477         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00478         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00479         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00480         if (dropZSmarkedPassed_)
00481           if (i->zsMarkAndPass()) continue;
00482 
00483         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00484         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00485         const HcalQIEShape* shape = conditions->getHcalShape (channelCoder);
00486         HcalCoderDb coder (*channelCoder, *shape);
00487 
00488         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00489 
00490         // Set auxiliary flag
00491         int auxflag=0;
00492         int fTS = firstAuxTS_;
00493         if (fTS<0) fTS=0; //silly protection against negative time slice values
00494         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00495           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
00496         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00497         auxflag+=((i->sample(fTS).capid())<<28);
00498         (rec->back()).setAux(auxflag);
00499 
00500         (rec->back()).setFlags(0);
00501         // Fill Presample ADC flag
00502         if (fTS>0)
00503           (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);
00504 
00505         if (setSaturationFlags_)
00506           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00507         if (correctTiming_)
00508           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00509       }
00510       // return result
00511       e.put(rec);    
00512 
00513       // HF -------------------------------------------------------------------
00514     } else if (subdet_==HcalForward) {
00515       edm::Handle<HFDigiCollection> digi;
00516       e.getByLabel(inputLabel_,digi);
00517 
00518 
00520       // create empty output
00521       std::auto_ptr<HFRecHitCollection> rec(new HFRecHitCollection);
00522       rec->reserve(digi->size());
00523       // run the algorithm
00524       HFDigiCollection::const_iterator i;
00525 
00526       // Vote on majority TS0 CapId
00527       int favorite_capid = 0; 
00528       if (correctTiming_) {
00529         long capid_votes[4] = {0,0,0,0};
00530         for (i=digi->begin(); i!=digi->end(); i++) {
00531           capid_votes[(*i)[0].capid()]++;
00532         }
00533         for (int k = 0; k < 4; k++)
00534           if (capid_votes[k] > capid_votes[favorite_capid])
00535             favorite_capid = k;
00536       }
00537 
00538       for (i=digi->begin(); i!=digi->end(); i++) {
00539         HcalDetId cell = i->id();
00540         DetId detcell=(DetId)cell;
00541 
00542         if(tsFromDB_ || recoParamsFromDB_) {
00543           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00544           if(tsFromDB_) {
00545             firstSample_  = param_ts->firstSample();
00546             samplesToAdd_ = param_ts->samplesToAdd();
00547           }
00548           if(recoParamsFromDB_) {
00549              bool correctForTimeslew=param_ts->correctForTimeslew();
00550              bool correctForPhaseContainment= param_ts->correctForPhaseContainment();
00551              float phaseNS=param_ts->correctionPhaseNS();
00552              useLeakCorrection_= param_ts->useLeakCorrection();
00553              correctTiming_ = param_ts->correctTiming();
00554              firstAuxTS_ = param_ts->firstAuxTS();
00555              int pileupCleaningID = param_ts->pileupCleaningID();
00556              reco_.setRecoParams(correctForTimeslew,correctForPhaseContainment,useLeakCorrection_,pileupCleaningID,phaseNS);
00557           }
00558         }
00559 
00560         int first = firstSample_;
00561         int toadd = samplesToAdd_;
00562 
00563         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00564         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00565         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00566         if (dropZSmarkedPassed_)
00567           if (i->zsMarkAndPass()) continue;
00568 
00569         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00570         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00571         const HcalQIEShape* shape = conditions->getHcalShape (channelCoder);
00572         HcalCoderDb coder (*channelCoder, *shape);
00573 
00574         // Set HFDigiTime flag values from digiTimeFromDB_
00575         if (digiTimeFromDB_==true && hfdigibit_!=0)
00576           {
00577             const HcalFlagHFDigiTimeParam* hfDTparam = HFDigiTimeParams->getValues(detcell.rawId());
00578             hfdigibit_->resetParamsFromDB(hfDTparam->HFdigiflagFirstSample(),
00579                                           hfDTparam->HFdigiflagSamplesToAdd(),
00580                                           hfDTparam->HFdigiflagExpectedPeak(),
00581                                           hfDTparam->HFdigiflagMinEThreshold(),
00582                                           hfDTparam->HFdigiflagCoefficients()
00583                                           );
00584           }
00585 
00586         //std::cout << "TOADDHF " << toadd << " " << first << " " << std::endl;
00587         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00588 
00589         // Set auxiliary flag
00590         int auxflag=0;
00591         int fTS = firstAuxTS_;
00592         if (fTS<0) fTS=0; // silly protection against negative time slice values
00593         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00594           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
00595         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00596         auxflag+=((i->sample(fTS).capid())<<28);
00597         (rec->back()).setAux(auxflag);
00598 
00599         // Clear flags
00600         (rec->back()).setFlags(0);
00601 
00602         // Fill Presample ADC flag
00603         if (fTS>0)
00604           (rec->back()).setFlagField((i->sample(fTS-1).adc()), HcalCaloFlagLabels::PresampleADC,7);
00605 
00606         // This calls the code for setting the HF noise bit determined from digi shape
00607         if (setNoiseFlags_) 
00608           hfdigibit_->hfSetFlagFromDigi(rec->back(),*i,coder,calibrations);
00609         if (setSaturationFlags_)
00610           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00611         if (setTimingTrustFlags_)
00612           HFTimingTrustFlagSetter_->setHFTimingTrustFlag(rec->back(),*i);
00613         if (correctTiming_)
00614           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00615       } // for (i=digi->begin(); i!=digi->end(); i++) -- loop on all HF digis
00616 
00617       // The following flags require the full set of rechits
00618       // These need to be set consecutively, so an energy check should be the first 
00619       // test performed on these hits (to minimize the loop time)
00620       if (setNoiseFlags_) 
00621         {
00622           // Step 1:  Set PET flag  (short fibers of |ieta|==29)
00623           // Neighbor/partner channels that are flagged by Pulse Shape algorithm (HFDigiTime)
00624           // won't be considered in these calculations
00625           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00626             {
00627               int depth=i->id().depth();
00628               int ieta=i->id().ieta();
00629               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00630               if (depth==2 || abs(ieta)==29 ) 
00631                 hfPET_->HFSetFlagFromPET(*i,*rec,myqual,mySeverity);
00632             }
00633 
00634           // Step 2:  Set S8S1 flag (short fibers or |ieta|==29)
00635           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00636             {
00637               int depth=i->id().depth();
00638               int ieta=i->id().ieta();
00639               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00640               if (depth==2 || abs(ieta)==29 ) 
00641                 hfS8S1_->HFSetFlagFromS9S1(*i,*rec,myqual,mySeverity);
00642             }
00643 
00644           // Set 3:  Set S9S1 flag (long fibers)
00645           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00646             {
00647               int depth=i->id().depth();
00648               int ieta=i->id().ieta();
00649               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00650               if (depth==1 && abs(ieta)!=29 ) 
00651                 hfS9S1_->HFSetFlagFromS9S1(*i,*rec,myqual, mySeverity);
00652             }
00653         }
00654 
00655       // return result
00656       e.put(rec);     
00657     } else if (subdet_==HcalOther && subdetOther_==HcalCalibration) {
00658       edm::Handle<HcalCalibDigiCollection> digi;
00659       e.getByLabel(inputLabel_,digi);
00660       
00661       // create empty output
00662       std::auto_ptr<HcalCalibRecHitCollection> rec(new HcalCalibRecHitCollection);
00663       rec->reserve(digi->size());
00664       // run the algorithm
00665       int first = firstSample_;
00666       int toadd = samplesToAdd_;
00667 
00668       HcalCalibDigiCollection::const_iterator i;
00669       for (i=digi->begin(); i!=digi->end(); i++) {
00670         HcalCalibDetId cell = i->id();
00671         //      HcalDetId cellh = i->id();
00672         DetId detcell=(DetId)cell;
00673         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00674         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00675         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00676         if (dropZSmarkedPassed_)
00677           if (i->zsMarkAndPass()) continue;
00678 
00679         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00680         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00681         const HcalQIEShape* shape = conditions->getHcalShape (channelCoder);
00682         HcalCoderDb coder (*channelCoder, *shape);
00683 
00684         // firstSample & samplesToAdd
00685         if(tsFromDB_) {
00686           const HcalRecoParam* param_ts = paramTS->getValues(detcell.rawId());
00687           first = param_ts->firstSample();    
00688           toadd = param_ts->samplesToAdd();    
00689         }
00690         rec->push_back(reco_.reconstruct(*i,first,toadd,coder,calibrations));
00691 
00692         /*
00693           // Flag setting not available for calibration rechits
00694         // Set auxiliary flag
00695         int auxflag=0;
00696         int fTS = firstAuxTS_;
00697         for (int xx=fTS; xx<fTS+4 && xx<i->size();++xx)
00698           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
00699         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00700         auxflag+=((i->sample(fTS).capid())<<28);
00701         (rec->back()).setAux(auxflag);
00702 
00703         (rec->back()).setFlags(0); // Not yet implemented for HcalCalibRecHit
00704         */
00705       }
00706       // return result
00707       e.put(rec);     
00708     }
00709   } 
00710   //DL  delete myqual;
00711 } // void HcalHitReconstructor::produce(...)
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc
