Reference Manual

Go to the documentation of this file.
00001 //using namespace std;
00002 #include "HcalHitReconstructor.h"
00003 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00004 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00005 #include "DataFormats/Common/interface/EDCollection.h"
00006 #include "DataFormats/Common/interface/Handle.h"
00007 #include "FWCore/Framework/interface/Selector.h"
00008 #include "FWCore/Framework/interface/ESHandle.h"
00009 #include "FWCore/Framework/interface/EventSetup.h"
00010 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
00011 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
00012 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
00013 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
00014 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
00015 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
00016 
00017 #include <iostream>
00018 
00019 /*  Hcal Hit reconstructor allows for CaloRecHits with status words */
00020 
00021 HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf):
00022   reco_(conf.getParameter<int>("firstSample"),
00023         conf.getParameter<int>("samplesToAdd"),
00024         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>("firstSample")+conf.getParameter<int>("firstAuxOffset"))
00037 {
00038   std::string subd=conf.getParameter<std::string>("Subdetector");
00039   //Set all FlagSetters to 0
00040   /* Important to do this!  Otherwise, if the setters are turned off,
00041      the "if (XSetter_) delete XSetter_;" commands can crash
00042   */
00043   hbheFlagSetter_             = 0;
00044   hbheHSCPFlagSetter_         = 0;
00045   hbhePulseShapeFlagSetter_   = 0;
00046   hbheTimingShapedFlagSetter_ = 0;
00047   hfdigibit_                  = 0;
00048 
00049   hfS9S1_                     = 0;
00050   hfS8S1_                     = 0;
00051   hfPET_                      = 0;
00052   saturationFlagSetter_       = 0;
00053   HFTimingTrustFlagSetter_    = 0;
00054   
00055   if (setSaturationFlags_)
00056     {
00057       const edm::ParameterSet& pssat      = conf.getParameter<edm::ParameterSet>("saturationParameters");
00058       saturationFlagSetter_ = new HcalADCSaturationFlag(pssat.getParameter<int>("maxADCvalue"));
00059     }
00060 
00061   if (!strcasecmp(subd.c_str(),"HBHE")) {
00062     subdet_=HcalBarrel;
00063     bool timingShapedCutsFlags = conf.getParameter<bool>("setTimingShapedCutsFlags");
00064     if (timingShapedCutsFlags)
00065       {
00066         const edm::ParameterSet& psTshaped = conf.getParameter<edm::ParameterSet>("timingshapedcutsParameters");
00067         hbheTimingShapedFlagSetter_ = new HBHETimingShapedFlagSetter(psTshaped.getParameter<std::vector<double> >("tfilterEnvelope"),
00068                                                                      psTshaped.getParameter<bool>("ignorelowest"),
00069                                                                      psTshaped.getParameter<bool>("ignorehighest"),
00070                                                                      psTshaped.getParameter<double>("win_offset"),
00071                                                                      psTshaped.getParameter<double>("win_gain"));
00072       }
00073       
00074     if (setNoiseFlags_)
00075       {
00076         const edm::ParameterSet& psdigi    =conf.getParameter<edm::ParameterSet>("flagParameters");
00077         hbheFlagSetter_=new HBHEStatusBitSetter(psdigi.getParameter<double>("nominalPedestal"),
00078                                                 psdigi.getParameter<double>("hitEnergyMinimum"),
00079                                                 psdigi.getParameter<int>("hitMultiplicityThreshold"),
00080                                                 psdigi.getParameter<std::vector<edm::ParameterSet> >("pulseShapeParameterSets"),
00081                                                 conf.getParameter<int>("firstSample"),
00082                                                 conf.getParameter<int>("samplesToAdd"));
00083       } // if (setNoiseFlags_)
00084     if (setHSCPFlags_)
00085       {
00086         const edm::ParameterSet& psHSCP = conf.getParameter<edm::ParameterSet>("hscpParameters");
00087         hbheHSCPFlagSetter_ = new HBHETimeProfileStatusBitSetter(psHSCP.getParameter<double>("r1Min"),
00088                                                                  psHSCP.getParameter<double>("r1Max"),
00089                                                                  psHSCP.getParameter<double>("r2Min"),
00090                                                                  psHSCP.getParameter<double>("r2Max"),
00091                                                                  psHSCP.getParameter<double>("fracLeaderMin"),
00092                                                                  psHSCP.getParameter<double>("fracLeaderMax"),
00093                                                                  psHSCP.getParameter<double>("slopeMin"),
00094                                                                  psHSCP.getParameter<double>("slopeMax"),
00095                                                                  psHSCP.getParameter<double>("outerMin"),
00096                                                                  psHSCP.getParameter<double>("outerMax"),
00097                                                                  psHSCP.getParameter<double>("TimingEnergyThreshold"));
00098       } // if (setHSCPFlags_) 
00099     if (setPulseShapeFlags_)
00100       {
00101         const edm::ParameterSet &psPulseShape = conf.getParameter<edm::ParameterSet>("pulseShapeParameters");
00102         hbhePulseShapeFlagSetter_ = new HBHEPulseShapeFlagSetter(
00103                                                                  psPulseShape.getParameter<double>("MinimumChargeThreshold"),
00104                                                                  psPulseShape.getParameter<unsigned int>("TrianglePeakTS"),
00105                                                                  psPulseShape.getParameter<std::vector<double> >("LinearThreshold"),
00106                                                                  psPulseShape.getParameter<std::vector<double> >("LinearCut"),
00107                                                                  psPulseShape.getParameter<std::vector<double> >("RMS8MaxThreshold"),
00108                                                                  psPulseShape.getParameter<std::vector<double> >("RMS8MaxCut"),
00109                                                                  psPulseShape.getParameter<std::vector<double> >("LeftSlopeThreshold"),
00110                                                                  psPulseShape.getParameter<std::vector<double> >("LeftSlopeCut"),
00111                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeThreshold"),
00112                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeCut"),
00113                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeSmallThreshold"),
00114                                                                  psPulseShape.getParameter<std::vector<double> >("RightSlopeSmallCut"),
00115                                                                  psPulseShape.getParameter<bool>("UseDualFit"),
00116                          psPulseShape.getParameter<bool>("TriangleIgnoreSlow"));
00117       }  // if (setPulseShapeFlags_)
00118 
00119     produces<HBHERecHitCollection>();
00120   } else if (!strcasecmp(subd.c_str(),"HO")) {
00121     subdet_=HcalOuter;
00122     produces<HORecHitCollection>();
00123   } else if (!strcasecmp(subd.c_str(),"HF")) {
00124     subdet_=HcalForward;
00125 
00126     if (setTimingTrustFlags_) {
00127       
00128       const edm::ParameterSet& pstrust      = conf.getParameter<edm::ParameterSet>("hfTimingTrustParameters");
00129       HFTimingTrustFlagSetter_=new HFTimingTrustFlag(pstrust.getParameter<int>("hfTimingTrustLevel1"),
00130                                                      pstrust.getParameter<int>("hfTimingTrustLevel2"));
00131     }
00132 
00133     if (setNoiseFlags_)
00134       {
00135         const edm::ParameterSet& psdigi    =conf.getParameter<edm::ParameterSet>("digistat");
00136         const edm::ParameterSet& psTimeWin =conf.getParameter<edm::ParameterSet>("HFInWindowStat");
00137         hfdigibit_=new HcalHFStatusBitFromDigis(conf.getParameter<int>("firstSample"),
00138                                                 conf.getParameter<int>("samplesToAdd"),
00139                                                 psdigi, psTimeWin);
00140 
00141         const edm::ParameterSet& psS9S1   = conf.getParameter<edm::ParameterSet>("S9S1stat");
00142         hfS9S1_   = new HcalHF_S9S1algorithm(psS9S1.getParameter<std::vector<double> >("short_optimumSlope"),
00143                                              psS9S1.getParameter<std::vector<double> >("shortEnergyParams"),
00144                                              psS9S1.getParameter<std::vector<double> >("shortETParams"),
00145                                              psS9S1.getParameter<std::vector<double> >("long_optimumSlope"),
00146                                              psS9S1.getParameter<std::vector<double> >("longEnergyParams"),
00147                                              psS9S1.getParameter<std::vector<double> >("longETParams"),
00148                                              psS9S1.getParameter<int>("flagsToSkip"),
00149                                              psS9S1.getParameter<bool>("isS8S1")
00150                                              );
00151 
00152         const edm::ParameterSet& psS8S1   = conf.getParameter<edm::ParameterSet>("S8S1stat");
00153         hfS8S1_   = new HcalHF_S9S1algorithm(psS8S1.getParameter<std::vector<double> >("short_optimumSlope"),
00154                                              psS8S1.getParameter<std::vector<double> >("shortEnergyParams"),
00155                                              psS8S1.getParameter<std::vector<double> >("shortETParams"),
00156                                              psS8S1.getParameter<std::vector<double> >("long_optimumSlope"),
00157                                              psS8S1.getParameter<std::vector<double> >("longEnergyParams"),
00158                                              psS8S1.getParameter<std::vector<double> >("longETParams"),
00159                                              psS8S1.getParameter<int>("flagsToSkip"),
00160                                              psS8S1.getParameter<bool>("isS8S1")
00161                                              );
00162 
00163         const edm::ParameterSet& psPET    = conf.getParameter<edm::ParameterSet>("PETstat");
00164         hfPET_    = new HcalHF_PETalgorithm(psPET.getParameter<std::vector<double> >("short_R"),
00165                                             psPET.getParameter<std::vector<double> >("shortEnergyParams"),
00166                                             psPET.getParameter<std::vector<double> >("shortETParams"),
00167                                             psPET.getParameter<std::vector<double> >("long_R"),
00168                                             psPET.getParameter<std::vector<double> >("longEnergyParams"),
00169                                             psPET.getParameter<std::vector<double> >("longETParams"),
00170                                             psPET.getParameter<int>("flagsToSkip"),
00171                                             psPET.getParameter<std::vector<double> >("short_R_29"),
00172                                             psPET.getParameter<std::vector<double> >("long_R_29")
00173                                             );
00174       }
00175     produces<HFRecHitCollection>();
00176   } else if (!strcasecmp(subd.c_str(),"ZDC")) {
00177     det_=DetId::Calo;
00178     subdet_=HcalZDCDetId::SubdetectorId;
00179     produces<ZDCRecHitCollection>();
00180   } else if (!strcasecmp(subd.c_str(),"CALIB")) {
00181     subdet_=HcalOther;
00182     subdetOther_=HcalCalibration;
00183     produces<HcalCalibRecHitCollection>();
00184   } else {
00185     std::cout << "HcalHitReconstructor is not associated with a specific subdetector!" << std::endl;
00186   }       
00187   
00188 }
00189 
00190 HcalHitReconstructor::~HcalHitReconstructor() {
00191   if (hbheFlagSetter_)        delete hbheFlagSetter_;
00192   if (hfdigibit_)             delete hfdigibit_;
00193   if (hbheHSCPFlagSetter_)    delete hbheHSCPFlagSetter_;
00194   if (hbhePulseShapeFlagSetter_) delete hbhePulseShapeFlagSetter_;
00195   if (hfS9S1_)                delete hfS9S1_;
00196   if (hfPET_)                 delete hfPET_;
00197 }
00198 
00199 void HcalHitReconstructor::produce(edm::Event& e, const edm::EventSetup& eventSetup)
00200 {
00201   //bool isData=e.isRealData(); // some flags should only be applied to real data
00202 
00203   // get conditions
00204   edm::ESHandle<HcalDbService> conditions;
00205   eventSetup.get<HcalDbRecord>().get(conditions);
00206   const HcalQIEShape* shape = conditions->getHcalShape (); // this one is generic
00207   
00208   edm::ESHandle<HcalChannelQuality> p;
00209   eventSetup.get<HcalChannelQualityRcd>().get(p);
00210   HcalChannelQuality* myqual = new HcalChannelQuality(*p.product());
00211 
00212   edm::ESHandle<HcalSeverityLevelComputer> mycomputer;
00213   eventSetup.get<HcalSeverityLevelComputerRcd>().get(mycomputer);
00214   const HcalSeverityLevelComputer* mySeverity = mycomputer.product();
00215   
00216   if (det_==DetId::Hcal) {
00217     if (subdet_==HcalBarrel || subdet_==HcalEndcap) {
00218       edm::Handle<HBHEDigiCollection> digi;
00219       
00220       e.getByLabel(inputLabel_,digi);
00221       
00222       // create empty output
00223       std::auto_ptr<HBHERecHitCollection> rec(new HBHERecHitCollection);
00224       rec->reserve(digi->size());
00225       // run the algorithm
00226       if (setNoiseFlags_) hbheFlagSetter_->Clear();
00227       HBHEDigiCollection::const_iterator i;
00228       std::vector<HBHEDataFrame> HBDigis;
00229       std::vector<int> RecHitIndex;
00230 
00231       // Vote on majority TS0 CapId
00232       int favorite_capid = 0; 
00233       if (correctTiming_) {
00234         long capid_votes[4] = {0,0,0,0};
00235         for (i=digi->begin(); i!=digi->end(); i++) {
00236           capid_votes[(*i)[0].capid()]++;
00237         }
00238         for (int k = 0; k < 4; k++)
00239           if (capid_votes[k] > capid_votes[favorite_capid])
00240             favorite_capid = k;
00241       }
00242 
00243       for (i=digi->begin(); i!=digi->end(); i++) {
00244         HcalDetId cell = i->id();
00245         DetId detcell=(DetId)cell;
00246 
00247         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00248         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00249         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00250         if (dropZSmarkedPassed_)
00251           if (i->zsMarkAndPass()) continue;
00252 
00253         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00254         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00255         HcalCoderDb coder (*channelCoder, *shape);
00256         rec->push_back(reco_.reconstruct(*i,coder,calibrations));
00257 
00258         // Set auxiliary flag
00259         int auxflag=0;
00260         for (int xx=firstauxTS_;xx<firstauxTS_+4 && xx<i->size();++xx)
00261           auxflag+=(i->sample(xx).adc())<<(7*(xx-firstauxTS_)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00262         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00263         auxflag+=((i->sample(firstauxTS_).capid())<<28);
00264         (rec->back()).setAux(auxflag);
00265 
00266         (rec->back()).setFlags(0);  // this sets all flag bits to 0
00267         if (hbheTimingShapedFlagSetter_!=0)
00268           hbheTimingShapedFlagSetter_->SetTimingShapedFlags(rec->back());
00269         if (setNoiseFlags_)
00270           hbheFlagSetter_->SetFlagsFromDigi(rec->back(),*i,coder,calibrations);
00271         if (setPulseShapeFlags_ == true)
00272           hbhePulseShapeFlagSetter_->SetPulseShapeFlags(rec->back(), *i, coder, calibrations);
00273         if (setSaturationFlags_)
00274           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00275         if (correctTiming_)
00276           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00277         if (setHSCPFlags_ && i->id().ietaAbs()<16)
00278           {
00279             double DigiEnergy=0;
00280             for(int j=0; j!=i->size(); DigiEnergy += i->sample(j++).nominal_fC());
00281             if(DigiEnergy > hbheHSCPFlagSetter_->EnergyThreshold())
00282               {
00283                 HBDigis.push_back(*i);
00284                 RecHitIndex.push_back(rec->size()-1);
00285               }
00286             
00287           } // if (set HSCPFlags_ && |ieta|<16)
00288       } // loop over HBHE digis
00289 
00290 
00291       if (setNoiseFlags_) hbheFlagSetter_->SetFlagsFromRecHits(*rec);
00292       if (setHSCPFlags_)  hbheHSCPFlagSetter_->hbheSetTimeFlagsFromDigi(rec.get(), HBDigis, RecHitIndex);
00293       // return result
00294       e.put(rec);
00295     } else if (subdet_==HcalOuter) {
00296       edm::Handle<HODigiCollection> digi;
00297       e.getByLabel(inputLabel_,digi);
00298       
00299       // create empty output
00300       std::auto_ptr<HORecHitCollection> rec(new HORecHitCollection);
00301       rec->reserve(digi->size());
00302       // run the algorithm
00303       HODigiCollection::const_iterator i;
00304 
00305       // Vote on majority TS0 CapId
00306       int favorite_capid = 0; 
00307       if (correctTiming_) {
00308         long capid_votes[4] = {0,0,0,0};
00309         for (i=digi->begin(); i!=digi->end(); i++) {
00310           capid_votes[(*i)[0].capid()]++;
00311         }
00312         for (int k = 0; k < 4; k++)
00313           if (capid_votes[k] > capid_votes[favorite_capid])
00314             favorite_capid = k;
00315       }
00316 
00317       for (i=digi->begin(); i!=digi->end(); i++) {
00318         HcalDetId cell = i->id();
00319         DetId detcell=(DetId)cell;
00320         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00321         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00322         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00323         if (dropZSmarkedPassed_)
00324           if (i->zsMarkAndPass()) continue;
00325 
00326         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00327         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00328         HcalCoderDb coder (*channelCoder, *shape);
00329         rec->push_back(reco_.reconstruct(*i,coder,calibrations));
00330 
00331         // Set auxiliary flag
00332         int auxflag=0;
00333         for (int xx=firstauxTS_;xx<firstauxTS_+4 && xx<i->size();++xx)
00334           auxflag+=(i->sample(xx).adc())<<(7*(xx-firstauxTS_)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00335         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00336         auxflag+=((i->sample(firstauxTS_).capid())<<28);
00337         (rec->back()).setAux(auxflag);
00338 
00339         (rec->back()).setFlags(0);
00340         if (setSaturationFlags_)
00341           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00342         if (correctTiming_)
00343           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00344       }
00345       // return result
00346       e.put(rec);    
00347     } else if (subdet_==HcalForward) {
00348       edm::Handle<HFDigiCollection> digi;
00349       e.getByLabel(inputLabel_,digi);
00350 
00351       // ugly hack only for purposes of 3.11 HF treatment
00352       if (e.isRealData() && e.run() <= 153943)
00353         {
00354           reco_.resetTimeSamples(3,4);
00355           if (hfdigibit_) hfdigibit_->resetTimeSamples(3,4);
00356           firstauxTS_=3; // hard-code starting position of aux word
00357         }
00358       else
00359         {
00360           reco_.resetTimeSamples(4,2);
00361           if (hfdigibit_) hfdigibit_->resetTimeSamples(3,3); // flag uses 3 TS, even if reco uses 2 TS
00362           firstauxTS_=3; // hard-code 
00363         }
00364 
00366       // create empty output
00367       std::auto_ptr<HFRecHitCollection> rec(new HFRecHitCollection);
00368       rec->reserve(digi->size());
00369       // run the algorithm
00370       HFDigiCollection::const_iterator i;
00371 
00372       // Vote on majority TS0 CapId
00373       int favorite_capid = 0; 
00374       if (correctTiming_) {
00375         long capid_votes[4] = {0,0,0,0};
00376         for (i=digi->begin(); i!=digi->end(); i++) {
00377           capid_votes[(*i)[0].capid()]++;
00378         }
00379         for (int k = 0; k < 4; k++)
00380           if (capid_votes[k] > capid_votes[favorite_capid])
00381             favorite_capid = k;
00382       }
00383 
00384       for (i=digi->begin(); i!=digi->end(); i++) {
00385         HcalDetId cell = i->id();
00386         DetId detcell=(DetId)cell;
00387         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00388         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00389         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00390         if (dropZSmarkedPassed_)
00391           if (i->zsMarkAndPass()) continue;
00392 
00393         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00394         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00395         HcalCoderDb coder (*channelCoder, *shape);
00396         rec->push_back(reco_.reconstruct(*i,coder,calibrations));
00397 
00398         // Set auxiliary flag
00399         int auxflag=0;
00400         for (int xx=firstauxTS_;xx<firstauxTS_+4 && xx<i->size();++xx)
00401           auxflag+=(i->sample(xx).adc())<<(7*(xx-firstauxTS_)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00402         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00403         auxflag+=((i->sample(firstauxTS_).capid())<<28);
00404         (rec->back()).setAux(auxflag);
00405 
00406         // Clear flags
00407         (rec->back()).setFlags(0);
00408         // This calls the code for setting the HF noise bit determined from digi shape
00409         if (setNoiseFlags_) hfdigibit_->hfSetFlagFromDigi(rec->back(),*i,coder,calibrations);
00410         if (setSaturationFlags_)
00411           saturationFlagSetter_->setSaturationFlag(rec->back(),*i);
00412         if (setTimingTrustFlags_)
00413           HFTimingTrustFlagSetter_->setHFTimingTrustFlag(rec->back(),*i);
00414         if (correctTiming_)
00415           HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid);
00416       } // for (i=digi->begin(); i!=digi->end(); i++) -- loop on all HF digis
00417 
00418       // The following flags require the full set of rechits
00419       // These need to be set consecutively, so an energy check should be the first 
00420       // test performed on these hits (to minimize the loop time)
00421       if (setNoiseFlags_) 
00422         {
00423           // Step 1:  Set PET flag  (short fibers of |ieta|==29)
00424           // Neighbor/partner channels that are flagged by Pulse Shape algorithm (HFDigiTime)
00425           // won't be considered in these calculations
00426           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00427             {
00428               int depth=i->id().depth();
00429               int ieta=i->id().ieta();
00430               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00431               if (depth==2 || abs(ieta)==29 ) 
00432                 hfPET_->HFSetFlagFromPET(*i,*rec,myqual,mySeverity);
00433             }
00434 
00435           // Step 2:  Set S8S1 flag (short fibers or |ieta|==29)
00436           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00437             {
00438               int depth=i->id().depth();
00439               int ieta=i->id().ieta();
00440               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00441               if (depth==2 || abs(ieta)==29 ) 
00442                 hfS8S1_->HFSetFlagFromS9S1(*i,*rec,myqual,mySeverity);
00443             }
00444 
00445           // Set 3:  Set S9S1 flag (long fibers)
00446           for (HFRecHitCollection::iterator i = rec->begin();i!=rec->end();++i)
00447             {
00448               int depth=i->id().depth();
00449               int ieta=i->id().ieta();
00450               // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
00451               if (depth==1 && abs(ieta)!=29 ) 
00452                 hfS9S1_->HFSetFlagFromS9S1(*i,*rec,myqual, mySeverity);
00453             }
00454         }
00455 
00456       // return result
00457       e.put(rec);     
00458     } else if (subdet_==HcalOther && subdetOther_==HcalCalibration) {
00459       edm::Handle<HcalCalibDigiCollection> digi;
00460       e.getByLabel(inputLabel_,digi);
00461       
00462       // create empty output
00463       std::auto_ptr<HcalCalibRecHitCollection> rec(new HcalCalibRecHitCollection);
00464       rec->reserve(digi->size());
00465       // run the algorithm
00466       HcalCalibDigiCollection::const_iterator i;
00467       for (i=digi->begin(); i!=digi->end(); i++) {
00468         HcalCalibDetId cell = i->id();
00469         DetId detcell=(DetId)cell;
00470         // check on cells to be ignored and dropped: (rof,20.Feb.09)
00471         const HcalChannelStatus* mydigistatus=myqual->getValues(detcell.rawId());
00472         if (mySeverity->dropChannel(mydigistatus->getValue() ) ) continue;
00473         if (dropZSmarkedPassed_)
00474           if (i->zsMarkAndPass()) continue;
00475 
00476         const HcalCalibrations& calibrations=conditions->getHcalCalibrations(cell);
00477         const HcalQIECoder* channelCoder = conditions->getHcalCoder (cell);
00478         HcalCoderDb coder (*channelCoder, *shape);
00479         rec->push_back(reco_.reconstruct(*i,coder,calibrations));
00480 
00481         /*
00482           // Flag setting not available for calibration rechits
00483         // Set auxiliary flag
00484         int auxflag=0;
00485         for (int xx=firstauxTS_;xx<firstauxTS_+4 && xx<i->size();++xx)
00486           auxflag+=(i->sample(xx).adc())<<(7*(xx-firstauxTS_)); // store the time slices in the first 28 bits of aux, a set of 4 7-bit adc values
00487         // bits 28 and 29 are reserved for capid of the first time slice saved in aux
00488         auxflag+=((i->sample(firstauxTS_).capid())<<28);
00489         (rec->back()).setAux(auxflag);
00490 
00491         (rec->back()).setFlags(0); // Not yet implemented for HcalCalibRecHit
00492         */
00493       }
00494       // return result
00495       e.put(rec);     
00496     }
00497   } 
00498 
00499   delete myqual;
00500 } // void HcalHitReconstructor::produce(...)
/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc
