Reference Manual

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