CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/DQM/HcalMonitorTasks/src/HcalDetDiagNoiseMonitor.cc

Go to the documentation of this file.
00001 #include "DQM/HcalMonitorTasks/interface/HcalDetDiagNoiseMonitor.h"
00002 #include "DataFormats/FEDRawData/interface/FEDRawData.h"
00003 #include "DataFormats/FEDRawData/interface/FEDRawDataCollection.h"
00004 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
00005 #include "DataFormats/HcalDigi/interface/HcalCalibrationEventTypes.h"
00006 #include "EventFilter/HcalRawToDigi/interface/HcalDCCHeader.h"
00007 #include "FWCore/Common/interface/TriggerNames.h"
00008 
00009 #include "TFile.h"
00010 #include "TTree.h"
00011 #include <TVector2.h>
00012 #include <TVector3.h>
00013 
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015 
00016 #include "DataFormats/Math/interface/deltaR.h"
00017 #include "DataFormats/Common/interface/TriggerResults.h"
00018 
00019 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
00020 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
00021 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
00022 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
00023 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
00024 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
00025 #include "DataFormats/L1Trigger/interface/L1EtMissParticle.h"
00026 #include "DataFormats/L1Trigger/interface/L1EtMissParticleFwd.h"
00027 #include "DataFormats/L1GlobalCaloTrigger/interface/L1GctHFRingEtSums.h"
00028 #include "DataFormats/L1GlobalCaloTrigger/interface/L1GctHFBitCounts.h"
00029 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTExtendedCand.h"
00030 #include "DataFormats/Candidate/interface/Candidate.h"
00031 #include "L1Trigger/RegionalCaloTrigger/interface/L1RCTProducer.h" 
00032 #include "DataFormats/L1CaloTrigger/interface/L1CaloCollections.h"
00033 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
00034 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00035 #include "HLTrigger/HLTanalyzers/interface/JetUtil.h"
00036 #include "DataFormats/METReco/interface/CaloMET.h"
00037 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00038 // #include "DataFormats/L1Trigger/interface/L1ParticleMap.h"
00039 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00040 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
00041 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
00042 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapFwd.h"
00043 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
00044 //#include "DataFormats/L1GlobalTrigger/interface/L1GtLogicParser.h"
00045 #include "DataFormats/TrackReco/interface/Track.h"
00046 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00047 #include "DataFormats/TrackReco/interface/TrackBase.h"
00048 #include "RecoMET/METAlgorithms/interface/HcalNoiseRBXArray.h"
00049 #include "DataFormats/VertexReco/interface/Vertex.h"
00050 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00051 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00052 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00053 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00054 #include "DQMServices/Core/interface/DQMStore.h"
00055 #include "DQMServices/Core/interface/MonitorElement.h"
00056 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
00057 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
00058 
00059 // this is to retrieve HCAL LogicalMap
00060 #include "CalibCalorimetry/HcalAlgos/interface/HcalLogicalMapGenerator.h"
00061 #include "CondFormats/HcalObjects/interface/HcalLogicalMap.h"
00062 // to retrive trigger information (local runs only)
00063 #include "TBDataFormats/HcalTBObjects/interface/HcalTBTriggerData.h"
00064 
00065 
00066 #include <math.h>
00067 
00068 using namespace reco;
00069 
00071 static const float adc2fC[128]={-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5, 10.5,11.5,12.5,
00072                    13.5,15.,17.,19.,21.,23.,25.,27.,29.5,32.5,35.5,38.5,42.,46.,50.,54.5,59.5,
00073                    64.5,59.5,64.5,69.5,74.5,79.5,84.5,89.5,94.5,99.5,104.5,109.5,114.5,119.5,
00074                    124.5,129.5,137.,147.,157.,167.,177.,187.,197.,209.5,224.5,239.5,254.5,272.,
00075                    292.,312.,334.5,359.5,384.5,359.5,384.5,409.5,434.5,459.5,484.5,509.5,534.5,
00076                    559.5,584.5,609.5,634.5,659.5,684.5,709.5,747.,797.,847.,897.,947.,997.,
00077                    1047.,1109.5,1184.5,1259.5,1334.5,1422.,1522.,1622.,1734.5,1859.5,1984.5,
00078                    1859.5,1984.5,2109.5,2234.5,2359.5,2484.5,2609.5,2734.5,2859.5,2984.5,
00079                    3109.5,3234.5,3359.5,3484.5,3609.5,3797.,4047.,4297.,4547.,4797.,5047.,
00080                    5297.,5609.5,5984.5,6359.5,6734.5,7172.,7672.,8172.,8734.5,9359.5,9984.5};
00082 static std::string subdets[11]={"HBM","HBP","HEM","HEP","HO1M","HO0","HO1P","HO2M","HO2P","HFM","HFP"};
00083 static std::string HB_RBX[36]={
00084 "HBM01","HBM02","HBM03","HBM04","HBM05","HBM06","HBM07","HBM08","HBM09","HBM10","HBM11","HBM12","HBM13","HBM14","HBM15","HBM16","HBM17","HBM18",
00085 "HBP01","HBP02","HBP03","HBP04","HBP05","HBP06","HBP07","HBP08","HBP09","HBP10","HBP11","HBP12","HBP13","HBP14","HBP15","HBP16","HBP17","HBP18"};
00086 static std::string HE_RBX[36]={
00087 "HEM01","HEM02","HEM03","HEM04","HEM05","HEM06","HEM07","HEM08","HEM09","HEM10","HEM11","HEM12","HEM13","HEM14","HEM15","HEM16","HEM17","HEM18",
00088 "HEP01","HEP02","HEP03","HEP04","HEP05","HEP06","HEP07","HEP08","HEP09","HEP10","HEP11","HEP12","HEP13","HEP14","HEP15","HEP16","HEP17","HEP18"};
00089 static std::string HO_RBX[36]={
00090 "HO2M02","HO2M04","HO2M06","HO2M08","HO2M10","HO2M12","HO1M02","HO1M04","HO1M06","HO1M08","HO1M10","HO1M12",
00091 "HO001","HO002","HO003","HO004","HO005","HO006","HO007","HO008","HO009","HO010","HO011","HO012",
00092 "HO1P02","HO1P04","HO1P06","HO1P08","HO1P10","HO1P12","HO2P02","HO2P04","HO2P06","HO2P08","HO2P10","HO2P12",
00093 };
00094 
00095 
00096 class HcalDetDiagNoiseRMData{
00097 public:
00098   HcalDetDiagNoiseRMData(){
00099      reset();
00100      reset_LS();
00101   };
00102   void reset(){
00103     n_th_hi=n_th_300=n_pix_1=n_pix_8=pix=n_pix=0;
00104   }  
00105   void reset_LS(){
00106     n_th_hi_LS=n_th_300_LS=0;
00107   }
00108   int    n_th_hi;
00109   int    n_th_300;
00110   int    n_pix_1;
00111   int    n_pix_8;
00112   int    pix;
00113   int    n_pix; 
00114   int    n_th_hi_LS;
00115   int    n_th_300_LS;
00116 };
00117 class HcalDetDiagNoiseRMEvent{
00118 public:
00119   HcalDetDiagNoiseRMEvent(){
00120     reset();
00121   }
00122   void reset(){
00123    n_pix_hi=n_pix_lo=energy=n_zero=0;
00124   }
00125   int n_pix_hi;
00126   int n_pix_lo;
00127   int n_zero;
00128   float energy;
00129 };
00130 
00131 class HcalDetDiagNoiseRMSummary{
00132 public:
00133   HcalDetDiagNoiseRMSummary(){ 
00134      reset();
00135   }
00136   void reset(){
00137      for(int i=0;i<HcalFrontEndId::maxRmIndex;i++) rm[i].reset(); 
00138   }
00139   void reset_LS(){
00140      for(int i=0;i<HcalFrontEndId::maxRmIndex;i++) rm[i].reset_LS(); 
00141   }
00142   int GetRMindex(const std::string& rbx,int rm){
00143       if(rbx.substr(0,3)=="HO0"){
00144          int sect=atoi(rbx.substr(3,2).c_str());
00145          if(sect>12) return -1;
00146          if(rm==1 && (sect==2  || sect==3 || sect==6 || sect==7 || sect==10 || sect==11)) return -1;
00147          if(rm==4 && (sect==12 || sect==1 || sect==4 || sect==5 || sect==8  || sect==9 )) return -1;
00148       }
00149       if(rbx.substr(0,3)=="HO1" || rbx.substr(0,3)=="HO2"){ 
00150          int sect=atoi(rbx.substr(4,2).c_str());
00151          if(sect>12) return -1;
00152          if(sect==1 || sect==3 || sect==5 || sect==7 || sect==9 || sect==11) return -1;
00153       }
00154       HcalFrontEndId id(rbx,rm,1,1,1,1,1);
00155       if(id.rawId()==0) return -1;
00156       return id.rmIndex(); 
00157   }
00158   HcalDetDiagNoiseRMData rm[HcalFrontEndId::maxRmIndex];
00159 };
00160 
00161 HcalDetDiagNoiseMonitor::HcalDetDiagNoiseMonitor(const edm::ParameterSet& ps) 
00162 {
00163   ievt_=0;
00164   run_number=-1;
00165   NoisyEvents=0;
00166   LocalRun=false; 
00167   dataset_seq_number=1;
00168   FirstOrbit=FirstOrbitLS=0xFFFFFFFF;
00169   LastOrbit=LastOrbitLS=0;
00170 
00171   Online_                = ps.getUntrackedParameter<bool>("online",false);
00172   mergeRuns_             = ps.getUntrackedParameter<bool>("mergeRuns",false);
00173   enableCleanup_         = ps.getUntrackedParameter<bool>("enableCleanup",false);
00174   debug_                 = ps.getUntrackedParameter<int>("debug",0);
00175   prefixME_              = ps.getUntrackedParameter<std::string>("subSystemFolder","Hcal/");
00176   if (prefixME_.substr(prefixME_.size()-1,prefixME_.size())!="/")
00177     prefixME_.append("/");
00178   subdir_                = ps.getUntrackedParameter<std::string>("TaskFolder","DetDiagNoiseMonitor_Hcal"); 
00179   if (subdir_.size()>0 && subdir_.substr(subdir_.size()-1,subdir_.size())!="/")
00180     subdir_.append("/");
00181   subdir_=prefixME_+subdir_;
00182   AllowedCalibTypes_     = ps.getUntrackedParameter<std::vector<int> > ("AllowedCalibTypes");
00183   skipOutOfOrderLS_      = ps.getUntrackedParameter<bool>("skipOutOfOrderLS",false);
00184   NLumiBlocks_           = ps.getUntrackedParameter<int>("NLumiBlocks",4000);
00185   makeDiagnostics_       = ps.getUntrackedParameter<bool>("makeDiagnostics",false);
00186 
00187   UseDB            = ps.getUntrackedParameter<bool>  ("UseDB"  , false);
00188   OutputFilePath   = ps.getUntrackedParameter<std::string>("OutputFilePath", "");
00189   HPDthresholdHi   = ps.getUntrackedParameter<double>("NoiseThresholdHPDhi",49.0);
00190   HPDthresholdLo   = ps.getUntrackedParameter<double>("NoiseThresholdHPDlo",10.0);
00191   SpikeThreshold   = ps.getUntrackedParameter<double>("NoiseSpikeThreshold",0.5);
00192   Overwrite        = ps.getUntrackedParameter<bool>  ("Overwrite",true);
00193 
00194   rawDataLabel_ =  ps.getUntrackedParameter<edm::InputTag>("RawDataLabel",edm::InputTag("source",""));
00195   digiLabel_     = ps.getUntrackedParameter<edm::InputTag>("digiLabel",edm::InputTag("hcalDigis"));
00196   L1ADataLabel_  = ps.getUntrackedParameter<edm::InputTag>("gtLabel");
00197   
00198   RMSummary = 0;
00199 
00200 }
00201 
00202 void HcalDetDiagNoiseMonitor::cleanup(){
00203   if(dbe_){
00204     dbe_->setCurrentFolder(subdir_);
00205     dbe_->removeContents();
00206     dbe_ = 0;
00207   }
00208 } 
00209 void HcalDetDiagNoiseMonitor::reset(){}
00210 
00211 
00212 void HcalDetDiagNoiseMonitor::beginRun(const edm::Run& run, const edm::EventSetup& c)
00213 {
00214   if (debug_>1) std::cout <<"HcalDetDiagNoiseMonitor::beginRun"<<std::endl;
00215   HcalBaseDQMonitor::beginRun(run,c);
00216 
00217   if (tevt_==0) this->setup(); // set up histograms if they have not been created before
00218   if (mergeRuns_==false)
00219     this->reset();
00220 
00221   return;
00222 
00223 } 
00224 
00225 void HcalDetDiagNoiseMonitor::setup(){
00226   // Call base class setup
00227   HcalBaseDQMonitor::setup();
00228   if (!dbe_) return;
00229   RMSummary = new HcalDetDiagNoiseRMSummary();
00230 
00231   std::string name;
00232   if(dbe_!=NULL){    
00233      dbe_->setCurrentFolder(subdir_);   
00234      meEVT_ = dbe_->bookInt("HcalNoiseMonitor Event Number");
00235      dbe_->setCurrentFolder(subdir_+"Common Plots");
00236      
00237      name="RBX Pixel multiplicity";     PixelMult        = dbe_->book1D(name,name,73,0,73);
00238      name="HPD energy";                 HPDEnergy        = dbe_->book1D(name,name,200,0,2500);
00239      name="RBX energy";                 RBXEnergy        = dbe_->book1D(name,name,200,0,3500);
00240      name="Number of zero TS per RBX";  NZeroes          = dbe_->book1D(name,name,100,0,100);
00241      name="Trigger BX Tbit11";          TriggerBx11      = dbe_->book1D(name,name,4000,0,4000);
00242      name="Trigger BX Tbit12";          TriggerBx12      = dbe_->book1D(name,name,4000,0,4000);
00243 
00244      dbe_->setCurrentFolder(subdir_+"HBHE Plots");
00245      name="HBP HPD Noise Rate Pixel above 50fC"; HBP_Rate50    = dbe_->book1D(name,name,73,0,73);
00246      name="HBM HPD Noise Rate Pixel above 50fC"; HBM_Rate50    = dbe_->book1D(name,name,73,0,73);
00247      name="HEP HPD Noise Rate Pixel above 50fC"; HEP_Rate50    = dbe_->book1D(name,name,73,0,73);
00248      name="HEM HPD Noise Rate Pixel above 50fC"; HEM_Rate50    = dbe_->book1D(name,name,73,0,73);
00249      name="HBP HPD Noise Rate HPD above 300fC";  HBP_Rate300   = dbe_->book1D(name,name,73,0,73);
00250      name="HBM HPD Noise Rate HPD above 300fC";  HBM_Rate300   = dbe_->book1D(name,name,73,0,73);
00251      name="HEP HPD Noise Rate HPD above 300fC";  HEP_Rate300   = dbe_->book1D(name,name,73,0,73);
00252      name="HEM HPD Noise Rate HPD above 300fC";  HEM_Rate300   = dbe_->book1D(name,name,73,0,73);
00253 
00254      dbe_->setCurrentFolder(subdir_+"HO Plots");
00255      name="HO0  HPD Noise Rate Pixel above 50fC"; HO0_Rate50   = dbe_->book1D(name,name,49,0,49);
00256      name="HO1P HPD Noise Rate Pixel above 50fC"; HO1P_Rate50   = dbe_->book1D(name,name,48,0,48);
00257      name="HO1M HPD Noise Rate Pixel above 50fC"; HO1M_Rate50   = dbe_->book1D(name,name,48,0,48);
00258      name="HO0 HPD Noise Rate HPD above 300fC";   HO0_Rate300  = dbe_->book1D(name,name,48,0,48);
00259      name="HO1P HPD Noise Rate HPD abGetRMindexove 300fC";  HO1P_Rate300 = dbe_->book1D(name,name,48,0,48);
00260      name="HO1M HPD Noise Rate HPD above 300fC";  HO1M_Rate300 = dbe_->book1D(name,name,48,0,48);
00261       
00262 
00263      dbe_->setCurrentFolder(subdir_+"Noise Spike Plots");
00264 
00265      name="HB RM Spike Map";          HB_RBXmapSpikeCnt= dbe_->book2D(name,name,4,0.5,4.5,36,0.5,36.5);
00266      name="HE RM Spike Map";          HE_RBXmapSpikeCnt= dbe_->book2D(name,name,4,0.5,4.5,36,0.5,36.5);
00267      name="HO RM Spike Map";          HO_RBXmapSpikeCnt= dbe_->book2D(name,name,4,0.5,4.5,36,0.5,36.5);
00268 
00269      std::string title="RM";
00270      HB_RBXmapSpikeCnt->setAxisTitle(title);
00271      HE_RBXmapSpikeCnt->setAxisTitle(title);
00272      HO_RBXmapSpikeCnt->setAxisTitle(title);
00273  
00274      for(int i=0;i<36;i++){
00275         HB_RBXmapSpikeCnt->setBinLabel(i+1,HB_RBX[i],2);
00276         HE_RBXmapSpikeCnt->setBinLabel(i+1,HE_RBX[i],2);
00277         HO_RBXmapSpikeCnt->setBinLabel(i+1,HO_RBX[i],2);
00278      }
00279   } 
00280 
00281   gen =new HcalLogicalMapGenerator();
00282   lmap =new HcalLogicalMap(gen->createMap());
00283 
00284   return;
00285 } 
00286 
00287 void HcalDetDiagNoiseMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup){
00288   if (!IsAllowedCalibType()) return;
00289   if (LumiInOrder(iEvent.luminosityBlock())==false) return;
00290   HcalBaseDQMonitor::analyze(iEvent, iSetup);
00291   bool isNoiseEvent=false;  
00292   if(!dbe_) return;
00293   int orbit=-1111;
00294   int bx=-1111;
00295 
00296   // for local runs 
00297   edm::Handle<HcalTBTriggerData> trigger_data;
00298   iEvent.getByType(trigger_data);
00299   if(trigger_data.isValid()){
00300       if(trigger_data->triggerWord()>1000) isNoiseEvent=true;
00301       LocalRun=true;
00302   }
00303 
00304   // We do not want to look at Abort Gap events
00305   edm::Handle<FEDRawDataCollection> rawdata;
00306   iEvent.getByLabel(rawDataLabel_,rawdata);
00307   //checking FEDs for calibration information
00308   for(int i=FEDNumbering::MINHCALFEDID;i<=FEDNumbering::MAXHCALFEDID; i++) {
00309       const FEDRawData& fedData = rawdata->FEDData(i) ;
00310       if ( fedData.size() < 24 ) continue ;
00311       orbit= ((const HcalDCCHeader*)(fedData.data()))->getOrbitNumber();
00312       bx=((const HcalDCCHeader*)(fedData.data()))->getBunchId();
00313       if(((const HcalDCCHeader*)(fedData.data()))->getCalibType()!=hc_Null) return;
00314   }
00315 
00316   // Check GCT trigger bits
00317   edm::Handle< L1GlobalTriggerReadoutRecord > gtRecord;
00318   iEvent.getByLabel(L1ADataLabel_, gtRecord);
00319   if(gtRecord.isValid()){
00320     const TechnicalTriggerWord tWord = gtRecord->technicalTriggerWord();
00321     if(tWord.at(11) || tWord.at(12)) isNoiseEvent=true;
00322     if(tWord.at(11)){ TriggerBx11->Fill(bx);}
00323     if(tWord.at(12)){ TriggerBx12->Fill(bx);}
00324   }
00325  
00326   if(!isNoiseEvent) return;
00327   if(ievt_==0){ FirstOrbit=orbit; FirstOrbitLS=orbit; newLS=true;}
00328   if(LastOrbit <orbit) LastOrbit=orbit; 
00329   if(FirstOrbit>orbit) FirstOrbit=orbit;
00330   if(LastOrbitLS <orbit) LastOrbitLS=orbit; 
00331   if(FirstOrbitLS>orbit) FirstOrbitLS=orbit;
00332   if(newLS){ 
00333      FirstOrbitLS=orbit; 
00334      newLS=false;
00335   }
00336 
00337   if(!LocalRun){
00338      double TIME=(double)(LastOrbit-FirstOrbit)/11223.0;
00339      if(TIME>1800.0){
00340         UpdateHistos();
00341         SaveRates();
00342         RMSummary->reset();
00343         FirstOrbit=orbit; 
00344      }
00345   }
00346 
00347   meEVT_->Fill(++ievt_);
00348 
00349   run_number=iEvent.id().run();
00350 
00351   HcalDetDiagNoiseRMEvent RMs[HcalFrontEndId::maxRmIndex];
00352    
00353    edm::Handle<HBHEDigiCollection> hbhe; 
00354    iEvent.getByLabel(digiLabel_,hbhe);
00355    for(HBHEDigiCollection::const_iterator digi=hbhe->begin();digi!=hbhe->end();digi++){
00356      double max=-100,sum,energy=0; int n_zero=0;
00357      for(int i=0;i<digi->size();i++){
00358        sum=adc2fC[digi->sample(i).adc()&0xff]; 
00359        if(max<adc2fC[digi->sample(i).adc()&0xff]) max=adc2fC[digi->sample(i).adc()&0xff];
00360        if(adc2fC[digi->sample(i).adc()&0xff]==0) n_zero++;
00361      }
00362      HcalFrontEndId lmap_entry=lmap->getHcalFrontEndId(digi->id());
00363      int index=lmap_entry.rmIndex(); if(index>=HcalFrontEndId::maxRmIndex) continue;
00364      RMs[index].n_zero++;
00365      if(max>HPDthresholdLo){
00366        for(int i=0;i<digi->size();i++) energy+=adc2fC[digi->sample(i).adc()&0xff]-2.5;
00367        RMs[index].n_pix_lo++;
00368        if(max>HPDthresholdHi){ RMs[index].n_pix_hi++; isNoiseEvent=true;}
00369        RMs[index].energy+=energy;
00370      }
00371    }
00372 
00373    edm::Handle<HODigiCollection> ho; 
00374    iEvent.getByLabel(digiLabel_,ho);
00375    for(HODigiCollection::const_iterator digi=ho->begin();digi!=ho->end();digi++){
00376      double max=-100,energy=0; int Eta=digi->id().ieta(); int Phi=digi->id().iphi(); int n_zero=0;
00377      for(int i=0;i<digi->size()-1;i++){
00378        if(max<adc2fC[digi->sample(i).adc()&0xff]) max=adc2fC[digi->sample(i).adc()&0xff];
00379        if(adc2fC[digi->sample(i).adc()&0xff]==0) n_zero++;
00380      }
00381      if((Eta>=11 && Eta<=15 && Phi>=59 && Phi<=70) || (Eta>=5 && Eta<=10 && Phi>=47 && Phi<=58)){
00382        continue; // ignory SiPMs
00383      }else{
00384        HcalFrontEndId lmap_entry=lmap->getHcalFrontEndId(digi->id());
00385        int index=lmap_entry.rmIndex(); if(index>=HcalFrontEndId::maxRmIndex) continue;
00386        RMs[index].n_zero++;
00387        if(max>HPDthresholdLo){
00388          for(int i=0;i<digi->size();i++) energy+=adc2fC[digi->sample(i).adc()&0xff]-2.5;
00389          RMs[index].n_pix_lo++;
00390          if(max>HPDthresholdHi){ RMs[index].n_pix_hi++; isNoiseEvent=true;}
00391          RMs[index].energy+=energy;
00392        }
00393      }                    
00394    }   
00395 
00396    NoisyEvents++;
00397       
00398    // RMs loop
00399    for(int i=0;i<HcalFrontEndId::maxRmIndex;i++){
00400       if(RMs[i].n_pix_hi>0){
00401          HPDEnergy->Fill(RMs[i].energy);
00402          RMSummary->rm[i].n_th_hi++;
00403          RMSummary->rm[i].n_th_hi_LS++;
00404          if(RMs[i].energy>300) RMSummary->rm[i].n_th_300++;
00405          if(RMs[i].energy>300) RMSummary->rm[i].n_th_300_LS++;
00406          if(RMs[i].n_pix_lo>1) RMSummary->rm[i].n_pix_1++;
00407          if(RMs[i].n_pix_lo>8) RMSummary->rm[i].n_pix_8++;
00408          RMSummary->rm[i].pix+=RMs[i].n_pix_lo;
00409          RMSummary->rm[i].n_pix++;
00410       }
00411    }
00412 
00413    // RBX loop
00414    for(int sd=0;sd<7;sd++) for(int sect=1;sect<=18;sect++){
00415      std::stringstream tempss;
00416      tempss << std::setw(2) << std::setfill('0') << sect;
00417      std::string rbx= subdets[sd]+tempss.str();
00418      
00419      double rbx_energy=0;int pix_mult=0; int n_zero=0; bool isValidRBX=false;
00420      for(int rm=1;rm<=4;rm++){
00421        int index=RMSummary->GetRMindex(rbx,rm);
00422        if(index>0 && index<HcalFrontEndId::maxRmIndex){
00423          rbx_energy+=RMs[index].energy;
00424          pix_mult+=RMs[index].n_pix_lo; 
00425          n_zero+=RMs[index].n_zero; 
00426          isValidRBX=true;
00427        }
00428      }
00429      if(isValidRBX){
00430        PixelMult->Fill(pix_mult);
00431        RBXEnergy->Fill(rbx_energy);
00432        NZeroes->Fill(n_zero);
00433      }
00434    }
00435 
00436    if((ievt_%100)==0 && debug_>0)
00437      std::cout <<ievt_<<"\t"<<NoisyEvents<<std::endl;
00438    return;
00439 }
00440 
00441 void HcalDetDiagNoiseMonitor::endLuminosityBlock(const edm::LuminosityBlock& lumiSeg,const edm::EventSetup& c){
00442 int first_rbx=0,last_rbx=0;  
00443   //double TIME=(double)(LastOrbitLS-FirstOrbitLS)/11223.0;
00444   double TIME=23.0;
00445   newLS=true;
00446   if(TIME==0) return;
00447  
00448   for(int sd=0;sd<9;sd++){
00449           if(sd==0){ first_rbx=0;  last_rbx=18;} //HBM
00450           if(sd==1){ first_rbx=18; last_rbx=36;} //HBP
00451           if(sd==0 || sd==1){  // update HB plots
00452             for(int rbx=first_rbx;rbx<last_rbx;rbx++)for(int rm=1;rm<=4;rm++){
00453                int index=RMSummary->GetRMindex(HB_RBX[rbx],rm);
00454                if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00455                double val=RMSummary->rm[index].n_th_hi_LS/TIME;
00456                if(val>SpikeThreshold){
00457                   HB_RBXmapSpikeCnt->Fill(rm,rbx+1,1);
00458                   //printf("%s %i %f (%f)\n",HO_RBX[rbx].c_str(),rm,RMSummary->rm[index].n_th_hi_LS/TIME,TIME);
00459                }
00460                RMSummary->rm[index].reset_LS();
00461             }
00462           }
00463           if(sd==2){ first_rbx=0;  last_rbx=18;} //HEM
00464           if(sd==3){ first_rbx=18; last_rbx=36;} //HEP
00465           if(sd==2 || sd==3){  // update HB plots
00466             for(int rbx=first_rbx;rbx<last_rbx;rbx++)for(int rm=1;rm<=4;rm++){
00467               int index=RMSummary->GetRMindex(HE_RBX[rbx],rm);
00468               if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00469               double val=RMSummary->rm[index].n_th_hi_LS/TIME;
00470               if(val>SpikeThreshold){
00471                   HE_RBXmapSpikeCnt->Fill(rm,rbx+1,1);
00472                   //printf("%s %i %f (%f)\n",HO_RBX[rbx].c_str(),rm,RMSummary->rm[index].n_th_hi_LS/TIME,TIME);
00473               }
00474               RMSummary->rm[index].reset_LS();
00475             }
00476           }
00477           if(sd==4){ first_rbx=6;  last_rbx=12;}  //HO1M
00478           if(sd==5){ first_rbx=12;  last_rbx=24;} //HO0
00479           if(sd==6){ first_rbx=24;  last_rbx=30;} //HO1P
00480           if(sd>3){
00481             for(int rbx=first_rbx;rbx<last_rbx;rbx++)for(int rm=1;rm<=4;rm++){
00482               int index=RMSummary->GetRMindex(HO_RBX[rbx],rm);
00483               if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00484               double val=RMSummary->rm[index].n_th_hi_LS/TIME;
00485               if(val>SpikeThreshold){
00486                   HO_RBXmapSpikeCnt->Fill(rm,rbx+1,1);
00487                   //printf("%s %i %f (%f)\n",HO_RBX[rbx].c_str(),rm,RMSummary->rm[index].n_th_hi_LS/TIME,TIME);
00488               }
00489               RMSummary->rm[index].reset_LS();
00490            }
00491          }
00492     } //sd=0;sd<9
00493 }
00494 
00495 void HcalDetDiagNoiseMonitor::UpdateHistos(){ 
00496    int first_rbx=0; 
00497    double TIME=(double)(LastOrbit-FirstOrbit)/11223.0;
00498    if(TIME==0) return;
00499    for(int sd=0;sd<9;sd++){
00500       if(sd==0){ first_rbx=0; } //HBM
00501       if(sd==1){ first_rbx=18;} //HBP
00502       if(sd==0 || sd==1){  // update HB plots
00503           for(int rbx=0;rbx<18;rbx++)for(int rm=1;rm<=4;rm++){
00504              int index=RMSummary->GetRMindex(HB_RBX[rbx+first_rbx],rm);
00505              if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00506              if(sd==0){
00507                  HBM_Rate50->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_hi/TIME);
00508                  HBM_Rate300->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_300/TIME);
00509              }
00510              if(sd==1){
00511                  HBP_Rate50->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_hi/TIME);
00512                  HBP_Rate300->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_300/TIME);
00513              }
00514           }     
00515       }
00516       if(sd==2){ first_rbx=0;} //HEM
00517       if(sd==3){ first_rbx=18;} //HEP
00518       if(sd==2 || sd==3){  // update HB plots
00519           for(int rbx=0;rbx<18;rbx++)for(int rm=1;rm<=4;rm++){
00520              int index=RMSummary->GetRMindex(HE_RBX[rbx+first_rbx],rm);
00521              if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00522              if(sd==2){
00523                  HEM_Rate50->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_hi/TIME);
00524                  HEM_Rate300->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_300/TIME);
00525              }
00526              if(sd==3){
00527                  HEP_Rate50->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_hi/TIME);
00528                  HEP_Rate300->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_300/TIME);
00529              }
00530           }     
00531       }
00532       int n=0;
00533       if(sd==4){ first_rbx=6; n=6;}  //HO1M
00534       if(sd==5){ first_rbx=12;n=12;} //HO0
00535       if(sd==6){ first_rbx=24;n=6;} //HO1P
00536       if(sd>3){ // update HO plots
00537           for(int rbx=0;rbx<n;rbx++)for(int rm=1;rm<=4;rm++){
00538              int index=RMSummary->GetRMindex(HO_RBX[rbx+first_rbx],rm);
00539              if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00540              if(sd==4){
00541                  HO1M_Rate50->setBinContent(rbx*4*2+rm,RMSummary->rm[index].n_th_hi/TIME);
00542                  HO1M_Rate300->setBinContent(rbx*4*2+rm,RMSummary->rm[index].n_th_300/TIME);
00543              }
00544              if(sd==5){
00545                  HO0_Rate50->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_hi/TIME);
00546                  HO0_Rate300->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_300/TIME);
00547              } 
00548              if(sd==5){
00549                  HO1P_Rate50->setBinContent(rbx*4*2+rm,RMSummary->rm[index].n_th_hi/TIME);
00550                  HO1P_Rate300->setBinContent(rbx*4+rm,RMSummary->rm[index].n_th_300/TIME);
00551              } 
00552          }
00553       }
00554    } //sd=0;sd<9
00555 } 
00556 
00557 void HcalDetDiagNoiseMonitor::SaveRates(){
00558 char   RBX[20];
00559 int    RM;
00560 float VAL1,VAL2,VAL3,VAL4,VAL5;
00561 char str[500]; 
00562     double TIME=(double)(LastOrbit-FirstOrbit)/11223.0;
00563     if(TIME==0) return;
00564     if(OutputFilePath.size()>0){
00565        if(!Overwrite){
00566           sprintf(str,"%sHcalDetDiagNoiseData_run%06i_%i.root",OutputFilePath.c_str(),run_number,dataset_seq_number);
00567        }else{
00568           sprintf(str,"%sHcalDetDiagNoiseData.root",OutputFilePath.c_str());
00569        }
00570        TFile *theFile = new TFile(str, "RECREATE");
00571        if(!theFile->IsOpen()) return;
00572        theFile->cd();
00573        sprintf(str,"%d",run_number); TObjString run(str);    run.Write("run number");
00574        sprintf(str,"%d",ievt_);      TObjString events(str); events.Write("Total events processed");
00575        sprintf(str,"%d",dataset_seq_number);      TObjString dsnum(str);  dsnum.Write("Dataset number");
00576        Long_t t; t=time(0); strftime(str,30,"%F %T",localtime(&t)); TObjString tm(str);  tm.Write("Dataset creation time");
00577 
00578        TTree *tree   =new TTree("HCAL Noise data","HCAL Noise data");
00579        if(tree==0)   return;
00580        tree->Branch("RBX",            &RBX,     "RBX/C");
00581        tree->Branch("rm",             &RM,      "rm/I");
00582        tree->Branch("RATE_50",        &VAL1,    "RATE_50");
00583        tree->Branch("RATE_300",       &VAL2,    "RATE_300");
00584        tree->Branch("RATE_PIX1",      &VAL3,    "RATE_PIX1");
00585        tree->Branch("RATE_PIX8",      &VAL4,    "RATE_PIX8");
00586        tree->Branch("RATE_PIXMEAN",   &VAL5,    "RATE_PIXMEAN");
00587        for(int rbx=0;rbx<36;rbx++) for(int rm=1;rm<=4;rm++){
00588            int index=RMSummary->GetRMindex(HB_RBX[rbx],rm);
00589            if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00590                sprintf(RBX,"%s",HB_RBX[rbx].c_str());
00591                RM=rm;
00592                VAL1=RMSummary->rm[index].n_th_hi/TIME;
00593                VAL2=RMSummary->rm[index].n_th_300/TIME;
00594                VAL3=RMSummary->rm[index].n_pix_1/TIME;
00595                VAL4=RMSummary->rm[index].n_pix_8/TIME;
00596                if(RMSummary->rm[index].n_pix>0)VAL5=RMSummary->rm[index].pix/RMSummary->rm[index].n_pix; else VAL5=0;
00597                tree->Fill();
00598        }
00599        for(int rbx=0;rbx<36;rbx++) for(int rm=1;rm<=4;rm++){
00600            int index=RMSummary->GetRMindex(HE_RBX[rbx],rm);
00601            if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00602                sprintf(RBX,"%s",HE_RBX[rbx].c_str());
00603                RM=rm;
00604                VAL1=RMSummary->rm[index].n_th_hi/TIME;
00605                VAL2=RMSummary->rm[index].n_th_300/TIME;
00606                VAL3=RMSummary->rm[index].n_pix_1/TIME;
00607                VAL4=RMSummary->rm[index].n_pix_8/TIME;
00608                if(RMSummary->rm[index].n_pix>0)VAL5=RMSummary->rm[index].pix/RMSummary->rm[index].n_pix; else VAL5=0;
00609                tree->Fill();
00610        }
00611        for(int rbx=0;rbx<36;rbx++) for(int rm=1;rm<=4;rm++){
00612            int index=RMSummary->GetRMindex(HO_RBX[rbx],rm);
00613            if(index<0 || index>=HcalFrontEndId::maxRmIndex) continue;
00614                sprintf(RBX,"%s",HO_RBX[rbx].c_str());
00615                RM=rm;
00616                VAL1=RMSummary->rm[index].n_th_hi/TIME;
00617                VAL2=RMSummary->rm[index].n_th_300/TIME;
00618                VAL3=RMSummary->rm[index].n_pix_1/TIME;
00619                VAL4=RMSummary->rm[index].n_pix_8/TIME;
00620                if(RMSummary->rm[index].n_pix>0)VAL5=RMSummary->rm[index].pix/RMSummary->rm[index].n_pix; else VAL5=0;
00621                tree->Fill();
00622        }
00623        theFile->Write();
00624        theFile->Close();
00625        dataset_seq_number++;
00626 
00627    }
00628 }
00629 
00630 
00631 void HcalDetDiagNoiseMonitor::done(){}
00632  
00633 HcalDetDiagNoiseMonitor::~HcalDetDiagNoiseMonitor(){if(LocalRun) UpdateHistos(); SaveRates(); }
00634 
00635 DEFINE_FWK_MODULE (HcalDetDiagNoiseMonitor);