diff --git a/.gitignore b/.gitignore index df7ca6d..93d2033 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,10 @@ *.png *.root __init__.py +*.old +pippo* +*processDump.py *.pyc html/* +LSF* +*.list diff --git a/EcalTiming/data/crystalmap.txt b/EcalTiming/data/crystalmap.txt index 77c0241..ba02b12 100644 --- a/EcalTiming/data/crystalmap.txt +++ b/EcalTiming/data/crystalmap.txt @@ -1712,6 +1712,7 @@ -5 299 0 3076 24 4 -5 838863659 -5 300 0 3076 24 4 -5 838863660 -5 301 0 3201 25 1 -5 838863661 +-5 302 0 3201 25 1 -5 838863662 -5 303 0 3201 25 1 -5 838863663 -5 304 0 3201 25 1 -5 838863664 -5 305 0 3201 25 1 -5 838863665 @@ -57453,6 +57454,7 @@ 77 110 0 4287 33 63 77 838965870 77 111 0 4286 33 62 77 838965871 77 112 0 4286 33 62 77 838965872 +77 113 0 4286 33 62 77 838965873 77 114 0 4286 33 62 77 838965874 77 115 0 4286 33 62 77 838965875 77 116 0 4285 33 61 77 838965876 @@ -71331,6 +71333,11 @@ 50 98 1 6273 49 1 2 872438114 50 99 1 6273 49 1 1 872438115 50 100 1 6273 49 1 0 872438116 +51 1 1 6801 53 17 0 872438145 +51 2 1 6801 53 17 1 872438146 +51 3 1 6801 53 17 2 872438147 +51 4 1 6801 53 17 3 872438148 +51 5 1 6801 53 17 4 872438149 51 6 1 6800 53 16 5 872438150 51 7 1 6800 53 16 6 872438151 51 8 1 6800 53 16 7 872438152 @@ -71404,6 +71411,10 @@ 51 98 1 6145 48 1 2 872438242 51 99 1 6145 48 1 1 872438243 51 100 1 6145 48 1 0 872438244 +52 1 1 6801 53 17 0 872438273 +52 3 1 6801 53 17 2 872438275 +52 4 1 6801 53 17 3 872438276 +52 5 1 6801 53 17 4 872438277 52 6 1 6800 53 16 4 872438278 52 7 1 6800 53 16 6 872438279 52 8 1 6800 53 16 7 872438280 @@ -71477,6 +71488,11 @@ 52 98 1 6145 48 1 2 872438370 52 99 1 6145 48 1 1 872438371 52 100 1 6145 48 1 0 872438372 +53 1 1 6801 53 17 0 872438401 +53 2 1 6801 53 17 1 872438402 +53 3 1 6801 53 17 2 872438403 +53 4 1 6801 53 17 3 872438404 +53 5 1 6801 53 17 4 872438405 53 6 1 6800 53 16 4 872438406 53 7 1 6800 53 16 5 872438407 53 8 1 6800 53 16 6 872438408 @@ -71550,6 +71566,11 @@ 53 98 1 6145 48 1 2 872438498 53 99 1 6145 48 1 1 872438499 53 100 1 6145 48 1 0 872438500 +54 1 1 6801 53 17 0 872438529 +54 2 1 6801 53 17 0 872438530 +54 3 1 6801 53 17 1 872438531 +54 4 1 6801 53 17 2 872438532 +54 5 1 6801 53 17 3 872438533 54 6 1 6800 53 16 4 872438534 54 7 1 6800 53 16 5 872438535 54 8 1 6800 53 16 6 872438536 @@ -71623,6 +71644,11 @@ 54 98 1 6145 48 1 2 872438626 54 99 1 6145 48 1 1 872438627 54 100 1 6145 48 1 0 872438628 +55 1 1 6801 53 17 0 872438657 +55 2 1 6801 53 17 0 872438658 +55 3 1 6801 53 17 1 872438659 +55 4 1 6801 53 17 2 872438660 +55 5 1 6801 53 17 3 872438661 55 6 1 6800 53 16 4 872438662 55 7 1 6800 53 16 5 872438663 55 8 1 6800 53 16 6 872438664 @@ -71993,6 +72019,7 @@ 59 73 1 6158 48 14 26 872439241 59 74 1 6158 48 14 25 872439242 59 75 1 6158 48 14 24 872439243 +59 76 1 6157 48 13 23 872439244 59 77 1 6157 48 13 22 872439245 59 78 1 6157 48 13 21 872439246 59 79 1 6157 48 13 20 872439247 diff --git a/EcalTiming/interface/EcalCrystalTimingCalibration.h b/EcalTiming/interface/EcalCrystalTimingCalibration.h index e2c77fa..1a465f3 100644 --- a/EcalTiming/interface/EcalCrystalTimingCalibration.h +++ b/EcalTiming/interface/EcalCrystalTimingCalibration.h @@ -16,7 +16,6 @@ #define DS_RING 0x08 #define DS_CRYS 0x10 - class EcalCrystalTimingCalibration { public: diff --git a/EcalTiming/plugins/DummyRechitDigis.cc b/EcalTiming/plugins/DummyRechitDigis.cc new file mode 100644 index 0000000..f6033a5 --- /dev/null +++ b/EcalTiming/plugins/DummyRechitDigis.cc @@ -0,0 +1,234 @@ +// -*- C++ -*- +// +// Package: ECALlite/DummyRechitDigis +// Class: DummyRechitDigis +// +/**\class DummyRechitDigis DummyRechitDigis.cc ECALlite/DummyRechitDigis/plugins/DummyRechitDigis.cc + Description: [one line class summary] + Implementation: + [Notes on implementation] +*/ +// +// Original Author: Joshua Robert Hardenbrook +// Created: Mon, 01 Jun 2015 06:58:24 GMT +// +// + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/EgammaReco/interface/BasicCluster.h" +#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/CaloRecHit/interface/CaloID.h" +#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h" + +#include "CondFormats/EcalObjects/interface/EcalSampleMask.h" + +#include "EcalTiming/EcalTiming/plugins/EcalTimingCalibProducer.h" + +// +// class declaration +// + +class DummyRechitDigis : public edm::EDProducer { +public: + explicit DummyRechitDigis(const edm::ParameterSet&); + ~DummyRechitDigis(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + virtual void beginJob() override; + virtual void produce(edm::Event&, const edm::EventSetup&) override; + virtual void endJob() override; + + // ----------member data --------------------------- + edm::InputTag tag_barrelHitProducer_; + edm::InputTag tag_endcapHitProducer_; + const std::string barrelRecHitCollection_; + const std::string endcapRecHitCollection_; + edm::InputTag tag_barrelDigiProducer_; + edm::InputTag tag_endcapDigiProducer_; + const std::string barrelDigiCollection_; + const std::string endcapDigiCollection_; + const bool doDigi_; +}; + +DummyRechitDigis::DummyRechitDigis(const edm::ParameterSet& iConfig): + tag_barrelHitProducer_ (iConfig.getParameter< edm::InputTag > ("barrelHitProducer")), + tag_endcapHitProducer_ (iConfig.getParameter< edm::InputTag > ("endcapHitProducer")), + barrelRecHitCollection_ (iConfig.getUntrackedParameter("barrelRecHitCollection")), + endcapRecHitCollection_ (iConfig.getUntrackedParameter("endcapRecHitCollection")), + tag_barrelDigiProducer_ (iConfig.getParameter< edm::InputTag > ("barrelDigis")), + tag_endcapDigiProducer_ (iConfig.getParameter< edm::InputTag > ("endcapDigis")), + barrelDigiCollection_ (iConfig.getUntrackedParameter("barrelDigiCollection")), + endcapDigiCollection_ (iConfig.getUntrackedParameter("endcapDigiCollection")), + doDigi_ (iConfig.getUntrackedParameter("doDigi")) +{ + if(doDigi_) { + produces< EBDigiCollection >(barrelDigiCollection_); + produces< EEDigiCollection >(endcapDigiCollection_); + } + else { + produces< EcalRecHitCollection >(barrelRecHitCollection_); + produces< EcalRecHitCollection >(endcapRecHitCollection_); + } +} + +DummyRechitDigis::~DummyRechitDigis(){ } + +void DummyRechitDigis::produce(edm::Event& iEvent, const edm::EventSetup& iSetup){ + //std::cout << "\n-----------New Event ----------------\n " << std::endl; + using namespace edm; + // build an empty collection + // fake rechits + // handle to try to fill + edm::Handle barrelRecHitsHandle; + edm::Handle endcapRecHitsHandle; + // dummy collection to Put() + std::auto_ptr< EcalRecHitCollection > rechits_temp( new EcalRecHitCollection); + std::auto_ptr< EcalRecHitCollection > rechits_temp2( new EcalRecHitCollection); + + // EcalRecHit::EcalRecHit(const DetId& id, float energy, float time, uint32_t flags, uint32_t flagBits) + // add one rechit + EcalRecHit zero_rechit(0,0,0,0,0); + EcalRecHitCollection zero_collection; + zero_collection.push_back(zero_rechit); + *rechits_temp = zero_collection; + *rechits_temp2 = zero_collection; + + std::auto_ptr< EcalRecHitCollection > rechits_eb( new EcalRecHitCollection); + std::auto_ptr< EcalRecHitCollection > rechits_ee( new EcalRecHitCollection); + + // fake digis + // handle to try to fill + Handle digisEBHandle; + Handle digisEEHandle; + // dummy collection to Put() + std::auto_ptr outputEBDigiCollection( new EBDigiCollection ); + std::auto_ptr outputEEDigiCollection( new EEDigiCollection ); + + //Digi zero_digi; + EBDigiCollection ebfakecol; + EEDigiCollection eefakecol; + //ebfakecol.push_back(zerodigi); + //eefakecol.push_back(zerodigi); + + // fake empty collections + std::auto_ptr fakeEBDigiCollection( new EBDigiCollection ); + *fakeEBDigiCollection = ebfakecol; + std::auto_ptr fakeEEDigiCollection( new EEDigiCollection) ; + *fakeEEDigiCollection = eefakecol; + + if(!doDigi_) { + bool foundEBRechit = true; + bool foundEERechit = true; + std::cout<<"sono entrato nel !dodigi"<begin(); recHit_itr != barrelRecHitsHandle->end(); ++recHit_itr){ + std::cout<<"sono entrato nell'if eb"<size() > 0;} + // if you found the collection put it back into the event + iEvent.put( foundEBRechit ? rechits_eb : rechits_temp, barrelRecHitCollection_); + + // if you dont find the endcap rechits youre looking for, put in a fake one + try { + //std::cout<<"sono entrato nell'if ee dodigi"<size() > 0; + } + + iEvent.put( foundEERechit ? rechits_ee : rechits_temp2, endcapRecHitCollection_); + std::cout<<"Build fake digi!"<begin(); digiItr != digisEBHandle->end(); ++digiItr){ + //EBDetId id_crystal(digiItr->id()); + DetId id = digiItr->id(); + //std::cout<id()<size() > 0; + //std::cout<<"sono entrato nel catch eb"<begin(); digiItr != digisEEHandle->end(); ++digiItr){ + DetId id = digiItr->id(); + // std::cout<id()); + // std::cout<<"SUBDETID!! Endcap "<id()<id().subdetId() <<"\n"; + } + *outputEEDigiCollection = *(digisEEHandle.product()); + //std::cout<<"sono entrato nell'else ee"<size() > 0; + //std::cout<<"sono entrato nel catch ee"<("minRecHitEnergyStep")), _minRecHitEnergyNStep(iConfig.getParameter("minRecHitEnergyNStep")), - _energyThresholdOffsetEB(iConfig.getParameter("energyThresholdOffsetEB")), - _energyThresholdOffsetEE(iConfig.getParameter("energyThresholdOffsetEE")), + _energyThresholdOffsetEB(iConfig.getParameter("energyThresholdOffsetEB")), + _energyThresholdOffsetEE(iConfig.getParameter("energyThresholdOffsetEE")), + _chi2ThresholdOffsetEB(iConfig.getParameter("chi2ThresholdOffsetEB")),//added + _chi2ThresholdOffsetEE(iConfig.getParameter("chi2ThresholdOffsetEE")),//added _minEntries(iConfig.getParameter("minEntries")), _globalOffset(iConfig.getParameter("globalOffset")), _storeEvents(iConfig.getParameter("storeEvents")), @@ -101,9 +103,9 @@ bool EcalTimingCalibProducer::addRecHit(const EcalRecHit& recHit, EventTimeMap& //check if rechit is valid if(! recHit.checkFlags(_recHitFlags)) return false; - float energyThreshold = getEnergyThreshold(recHit.detid()) ; - if( recHit.energy() < (energyThreshold)) return false; // minRecHitEnergy in ADC for EB - //if(recHit.detid().subdetId() == EcalEndcap && recHit.energy() < 2 * (_minRecHitEnergy+_minRecHitEnergyStep*_iter)) return false; + float energyThreshold = getEnergyThreshold(recHit.detid()); + float chi2Threshold = getChi2Threshold(recHit.detid()); + if( (recHit.energy() < energyThreshold) || (recHit.chi2() > chi2Threshold)) return false; // add the EcalTimingEvent to the EcalCreateTimeCalibrations EcalTimingEvent timeEvent(recHit); @@ -171,6 +173,7 @@ bool EcalTimingCalibProducer::filter(edm::Event& iEvent, const edm::EventSetup& edm::ESHandle pG; iSetup.get().get(pG); EcalRingCalibrationTools::setCaloGeometry(&(*pG)); + endcapGeometry_ = pG->getSubdetectorGeometry(DetId::Ecal, EcalEndcap); barrelGeometry_ = pG->getSubdetectorGeometry(DetId::Ecal, EcalBarrel); @@ -199,12 +202,8 @@ bool EcalTimingCalibProducer::filter(edm::Event& iEvent, const edm::EventSetup& // recHit_itr is of type: edm::Handle::const_iterator for(auto recHit_itr = ebRecHitHandle->begin(); recHit_itr != ebRecHitHandle->end(); ++recHit_itr) { // add the recHit to the list of recHits used for calibration (with the relative information) + if(addRecHit(*recHit_itr, _eventTimeMap)) { - //EBDetId id(recHit.detid()); - // if(id.ieta() == -75 && id.iphi() == 119) { - // std::cout << "RawID\t" << id.rawId() << std::endl; - // return false; - // } timeEB.add(EcalTimingEvent(*recHit_itr), false); } } @@ -257,11 +256,11 @@ bool EcalTimingCalibProducer::filter(edm::Event& iEvent, const edm::EventSetup& EcalTimingEvent event = _isSplash ? correctGlobalOffset(it.second, splashDir, bunchCorr) : it.second; if(_makeEventPlots) plotRecHit(event); - _timeCalibMap[it.first].add(event,_storeEvents); + _timeCalibMap[it.first].add(event, _storeEvents); //Find the CCU(tower) that this crystal belongs to unsigned int elecID = getElecID(it.first); - _HWCalibrationMap[elecID].add(event,false); + _HWCalibrationMap[elecID].add(event, false); } @@ -282,8 +281,7 @@ void EcalTimingCalibProducer::endJob() // set the values in _calibConstants, _calibErrors, _offsetConstant #ifdef DEBUG - for (auto it : _HWCalibrationMap) - { + for (auto it : _HWCalibrationMap) { if (abs(it.second.mean()) > HW_UNIT * 1.5) std::cout << "HW: " << it.first << ' ' << it.second.mean() << std::endl; } #endif @@ -297,7 +295,6 @@ void EcalTimingCalibProducer::endJob() _timeCalibConstants.setValue(calibRecHit_itr->first.rawId(), correction); unsigned int ds = DS_NONE; - //TODO: This probably shouldn't be commented out. Move the stat check into the individual functions? if(calibRecHit_itr->second.num() > 50) { // check the asymmetry of the distribution: if asymmetric, dump the full set of events for further offline studies if(fabs(calibRecHit_itr->second.getSkewnessWithinNSigma(n_sigma, 10)) > _maxSkewnessForDump) { @@ -315,21 +312,19 @@ void EcalTimingCalibProducer::endJob() } unsigned int elecID = getElecID(calibRecHit_itr->first); - if( abs(_HWCalibrationMap[elecID].mean()) > HW_UNIT * 1.5) - { + if( abs(_HWCalibrationMap[elecID].mean()) > HW_UNIT * 1.5) { ds |= DS_CCU_OOT; } if((calibRecHit_itr->first.rawId() == EBCRYex) || - (calibRecHit_itr->first.rawId() == EECRYex)) ds |= DS_CRYS; + (calibRecHit_itr->first.rawId() == EECRYex)) ds |= DS_CRYS; int iRing = _ringTools.getRingIndexInSubdet(calibRecHit_itr->first); if(calibRecHit_itr->first.subdetId() == EcalBarrel && iRing == EBRING) ds |= DS_RING; else if(calibRecHit_itr->first.subdetId() == EcalEndcap && (iRing == EEmRING || iRing == EEpRING )) ds |= DS_RING; - if(ds != DS_NONE) - { + if(ds != DS_NONE) { int ix, iy, iz; if(calibRecHit_itr->first.subdetId() == EcalBarrel) { EBDetId id(calibRecHit_itr->first); @@ -344,7 +339,6 @@ void EcalTimingCalibProducer::endJob() } calibRecHit_itr->second.dumpToTree(dumpTree, ix, iy, iz, ds, elecID, iRing); } - // add filing Energy hists here } @@ -360,7 +354,6 @@ void EcalTimingCalibProducer::endJob() dumpCalibration(filename); sprintf(filename, "%s-corr.dat", _outputDumpFileName.substr(0, _outputDumpFileName.find(".root")).c_str()); //text file holding constants dumpCorrections(filename); - // save the xml } void EcalTimingCalibProducer::dumpCorrections(std::string filename) @@ -374,12 +367,12 @@ void EcalTimingCalibProducer::dumpCorrections(std::string filename) if(id_.subdetId() == EcalBarrel) { EBDetId id(id_); fout << id.ieta() << "\t" << id.iphi() << "\t" << 0 - << "\t" << calibRecHit_itr->second.getMeanWithinNSigma(2,10) << "\t" << calibRecHit_itr->second.stdDev() << "\t" << calibRecHit_itr->second.num() << "\t" << calibRecHit_itr->second.meanE() + << "\t" << calibRecHit_itr->second.getMeanWithinNSigma(2, 10) << "\t" << calibRecHit_itr->second.stdDev() << "\t" << calibRecHit_itr->second.num() << "\t" << calibRecHit_itr->second.meanE() << "\t" << id.rawId() << std::endl; } else { EEDetId id(id_); fout << id.ix() << "\t" << id.iy() << "\t" << id.zside() - << "\t" << calibRecHit_itr->second.getMeanWithinNSigma(2,10) << "\t" << calibRecHit_itr->second.stdDev() << "\t" << calibRecHit_itr->second.num() << "\t" << calibRecHit_itr->second.meanE() + << "\t" << calibRecHit_itr->second.getMeanWithinNSigma(2, 10) << "\t" << calibRecHit_itr->second.stdDev() << "\t" << calibRecHit_itr->second.num() << "\t" << calibRecHit_itr->second.meanE() << "\t" << id.rawId() << std::endl; } } @@ -410,7 +403,7 @@ void EcalTimingCalibProducer::dumpCalibration(std::string filename) for(unsigned int i = 0; i < _timeCalibConstants.endcapItems().size(); ++i) { EEDetId id(EEDetId::detIdFromDenseIndex(i)); // this is a stupid thing that I'm obliged to do due to the stupid structure of the ECAL container - fout << id.ix() << "\t" << id.iy() << "\t" << id.zside() << "\t" << _timeCalibConstants.endcapItems()[i] << "\t" << id.rawId() << std::endl; + fout << id.ix() << "\t" << id.iy() << "\t" << id.zside() << "\t" << _timeCalibConstants.endcapItems()[i] << "\t" << id.rawId() << std::endl; } fout.close(); } @@ -422,17 +415,20 @@ void EcalTimingCalibProducer::dumpCalibration(std::string filename) void EcalTimingCalibProducer::FillCalibrationCorrectionHists(EcalTimeCalibrationMap::const_iterator cal_itr) { - int ix,iy,iz; + int ix, iy, iz; int rawid = cal_itr->first.rawId(); if(cal_itr->first.subdetId() == EcalBarrel) { EBDetId id(cal_itr->first); // Fill Rechit Energy EneMapEB_->Fill(id.ieta(), id.iphi(), cal_itr->second.meanE()); // 2D energy map - TimeMapEB_->Fill(id.ieta(), id.iphi(), cal_itr->second.getMeanWithinNSigma(2,10)); // 2D time map - TimeErrorMapEB_->Fill(id.ieta(), id.iphi(), cal_itr->second.getMeanErrorWithinNSigma(2,10)); + TimeMapEB_->Fill(id.ieta(), id.iphi(), cal_itr->second.getMeanWithinNSigma(2, 10)); // 2D time map + TimeErrorMapEB_->Fill(id.ieta(), id.iphi(), cal_itr->second.getMeanErrorWithinNSigma(2, 10)); RechitEneEB_->Fill(cal_itr->second.meanE()); // 1D histogram - RechitTimeEB_->Fill(cal_itr->second.getMeanWithinNSigma(2,10)); // 1D histogram + RechitTimeEB_->Fill(cal_itr->second.getMeanWithinNSigma(2, 10)); // 1D histogram + + // Fill Occupancy histo + OccupancyEB_->Fill(id.ieta(), id.iphi(), cal_itr->second.getNumWithinNSigma(2, 10)); ix = id.ieta(); iy = id.iphi(); @@ -443,25 +439,30 @@ void EcalTimingCalibProducer::FillCalibrationCorrectionHists(EcalTimeCalibration EEDetId id(cal_itr->first); if(id.zside() < 0) { EneMapEEM_->Fill(id.ix(), id.iy(), cal_itr->second.meanE()); - TimeMapEEM_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanWithinNSigma(2,10)); - TimeErrorMapEEM_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanErrorWithinNSigma(2,10)); + TimeMapEEM_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanWithinNSigma(2, 10)); + TimeErrorMapEEM_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanErrorWithinNSigma(2, 10)); RechitEneEEM_->Fill(cal_itr->second.meanE()); - RechitTimeEEM_->Fill(cal_itr->second.getMeanWithinNSigma(2,10)); + RechitTimeEEM_->Fill(cal_itr->second.getMeanWithinNSigma(2, 10)); + + OccupancyEEM_->Fill(id.ix(), id.iy(), cal_itr->second.getNumWithinNSigma(2, 10)); + } else { EneMapEEP_->Fill(id.ix(), id.iy(), cal_itr->second.meanE()); - TimeMapEEP_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanWithinNSigma(2,10)); - TimeErrorMapEEP_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanErrorWithinNSigma(2,10)); + TimeMapEEP_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanWithinNSigma(2, 10)); + TimeErrorMapEEP_->Fill(id.ix(), id.iy(), cal_itr->second.getMeanErrorWithinNSigma(2, 10)); RechitEneEEP_->Fill(cal_itr->second.meanE()); - RechitTimeEEP_->Fill(cal_itr->second.getMeanWithinNSigma(2,10)); + RechitTimeEEP_->Fill(cal_itr->second.getMeanWithinNSigma(2, 10)); + + OccupancyEEP_->Fill(id.ix(), id.iy(), cal_itr->second.getNumWithinNSigma(2, 10)); + } ix = id.ix(); iy = id.iy(); iz = id.zside(); } - int iRing = _ringTools.getRingIndexInSubdet(cal_itr->first); cal_itr->second.dumpCalibToTree(timingTree, rawid, ix, iy, iz, getElecID(cal_itr->first), iRing); } @@ -470,6 +471,7 @@ void EcalTimingCalibProducer::FillHWCorrectionHists(EcalTimeCalibrationMap::cons { unsigned int elecID = getElecID(cal_itr->first); float time = _HWCalibrationMap[elecID].mean(); + if(cal_itr->first.subdetId() == EcalBarrel) { EBDetId id(cal_itr->first); // Fill Rechit Energy @@ -488,8 +490,8 @@ void EcalTimingCalibProducer::FillHWCorrectionHists(EcalTimeCalibrationMap::cons void EcalTimingCalibProducer::FillEnergyStabilityHists(EcalTimeCalibrationMap::const_iterator cal_itr, std::vector< std::pair > energyStability) { - int ix,iy,iz; - int rawid = cal_itr->first.rawId(); + int ix, iy, iz; + int rawid = cal_itr->first.rawId(); //choose which map to store in if(cal_itr->first.subdetId() == EcalBarrel) { EBDetId id(cal_itr->first); @@ -502,7 +504,6 @@ void EcalTimingCalibProducer::FillEnergyStabilityHists(EcalTimeCalibrationMap::c iy = id.iy(); iz = id.zside(); } - int iRing = _ringTools.getRingIndexInSubdet(cal_itr->first); // Add min_energy to the tree which gets filld inside the dump function @@ -514,10 +515,9 @@ void EcalTimingCalibProducer::FillEnergyStabilityHists(EcalTimeCalibrationMap::c if(energyStabilityTree->GetBranch("index") == NULL) energyStabilityTree->Branch("index", &index, "index/b"); energyStabilityTree->SetBranchAddress("index", &index); - for(auto it = energyStability.begin(); it!=energyStability.end(); it++) - { + for(auto it = energyStability.begin(); it != energyStability.end(); it++) { min_energy = it->first; - it->second->dumpCalibToTree(energyStabilityTree,rawid,ix,iy,iz,getElecID(cal_itr->first),iRing); + it->second->dumpCalibToTree(energyStabilityTree, rawid, ix, iy, iz, getElecID(cal_itr->first), iRing); index++; } @@ -540,6 +540,8 @@ void EcalTimingCalibProducer::initEventHists(TFileDirectory fdir) // void EcalTimingCalibProducer::initHists(TFileDirectory fdir) { + std::cout << "Initializing histos" << std::endl; + EneMapEB_ = fdir.make("EneMapEB", "RecHit Energy[GeV] EB profile map;i#eta; i#phi;E[GeV]", 171, -85, 86, 360, 1., 361.); TimeMapEB_ = fdir.make("TimeMapEB", "Mean Time [ns] EB profile map; i#eta; i#phi;Time[ns]", 171, -85, 86, 360, 1., 361.); @@ -566,11 +568,17 @@ void EcalTimingCalibProducer::initHists(TFileDirectory fdir) RechitEnergyTimeEB = fdir.make("RechitEnergyTimeEB", "RecHit Energy vs Time EB; Rechit Energy[GeV]; Time[ns]; Events", 200, 0.0, 100.0, 100, -15, 30); // HW Histograms - + HWTimeMapEEP_ = fdir.make("HWTimeMapEEP", "Mean HW Time[ns] profile map EE+;ix;iy; Time[ns]", 100, 1, 101, 100, 1, 101); HWTimeMapEEM_ = fdir.make("HWTimeMapEEM", "Mean HW Time[ns] profile map EE-;ix;iy; Time[ns]", 100, 1, 101, 100, 1, 101); HWTimeMapEB_ = fdir.make("HWTimeMapEB", "Mean HW Time[ns] EB profile map; i#eta; i#phi;Time[ns]", 171, -85, 86, 360, 1., 361.); + //Occupancy histograms + + OccupancyEB_ = fdir.make("OccupancyEB", "Occupancy EB; i#eta; i#phi; #Hits", 171, -85, 86, 361, 1., 361.); + OccupancyEEM_ = fdir.make("OccupancyEEM", "OccupancyEEM; iy; ix; #Hits", 100, 1, 101, 100, 1, 101); + OccupancyEEP_ = fdir.make("OccupancyEEP", "OccupancyEEP; iy; ix; #Hits", 100, 1, 101, 100, 1, 101); + } // diff --git a/EcalTiming/plugins/EcalTimingCalibProducer.h b/EcalTiming/plugins/EcalTimingCalibProducer.h index 596a3df..dd8f547 100644 --- a/EcalTiming/plugins/EcalTimingCalibProducer.h +++ b/EcalTiming/plugins/EcalTimingCalibProducer.h @@ -132,7 +132,7 @@ class EcalTimingCalibProducer : public edm::EDFilter EcalTimeCalibrationMap _timeCalibMap; ///< calibration map: contains the time shift for each crystal EventTimeMap _eventTimeMap; ///< container of recHits passing selection in the event (reset at each event) EcalHWCalibrationMap _HWCalibrationMap; //!< The keys for this map are EcalElectronicIds with xtalid = stripid = 1 - ///< calibration map for the CCU's (Hardware Constants). + ///< calibration map for the CCU's (Hardware Constants). // For finding averages for specific eta ring EcalCrystalTimingCalibration timeEEP; ///< global time calibration of EE+ @@ -149,9 +149,11 @@ class EcalTimingCalibProducer : public edm::EDFilter ~EcalTimingCalibProducer(); // default destructor - virtual void beginJob() override; + virtual void beginJob() override; virtual bool filter(edm::Event&, const edm::EventSetup&) override; - virtual void endJob() override; + virtual void endJob() override; + + private: // ----------member data --------------------------- /** @name Input Parameters @@ -168,11 +170,13 @@ class EcalTimingCalibProducer : public edm::EDFilter unsigned int _recHitMin; ///< require at least this many rec hits to count the event double _minRecHitEnergyStep; ///< to check step size to check energy stability double _minRecHitEnergyNStep; ///< number of steps to check energy stability - double _energyThresholdOffsetEB; ///< energy to add to the minimum energy thresholc - double _energyThresholdOffsetEE; ///< energy to add to the minimum energy thresholc + double _energyThresholdOffsetEB; ///< energy to add to the minimum energy thresholc + double _energyThresholdOffsetEE; ///< energy to add to the minimum energy thresholc + double _chi2ThresholdOffsetEB; //chi2 square thr for the Barrel --- Added + double _chi2ThresholdOffsetEE; //chi2 square thr for the Endcap --- Added unsigned int _minEntries; ///< require a minimum number of entries in a ring to do averages float _globalOffset; ///< time to subtract from every event - bool _storeEvents; + bool _storeEvents; bool _produceNewCalib; ///< true if you don't want to use the values in DB and what to extract new absolute calibrations, if false iteration does not work std::string _outputDumpFileName; ///< name of the output file for the calibration constants' dump float _maxSkewnessForDump; @@ -222,9 +226,21 @@ class EcalTimingCalibProducer : public edm::EDFilter float getEnergyThreshold(const DetId detid) { int iRing = _ringTools.getRingIndexInSubdet(detid); - return detid.subdetId() == EcalBarrel ? 13 * 0.04 + _energyThresholdOffsetEB : - 20 * (79.29 - 4.148 * iRing + 0.2442 * iRing * iRing ) / 1000 + _energyThresholdOffsetEE; + if (detid.subdetId() == EcalBarrel) { + return 13 * 0.04 + _energyThresholdOffsetEB; + } else { + return 20 * (79.29 - 4.148 * iRing + 0.2442 * iRing * iRing ) / 1000 + _energyThresholdOffsetEE; + } + } + float getChi2Threshold(const DetId detid) + { + if (detid.subdetId() == EcalBarrel) { + return _chi2ThresholdOffsetEB; + } else { + return _chi2ThresholdOffsetEE; + } } + std::map _CrysEnergyMap; edm::Service fileService_; @@ -279,11 +295,16 @@ class EcalTimingCalibProducer : public edm::EDFilter TH2F* RechitEnergyTimeEEM; TH2F* RechitEnergyTimeEEP; + //Occupancy plots + TH2D* OccupancyEB_; + TH2D* OccupancyEEM_; + TH2D* OccupancyEEP_; + EcalRingCalibrationTools _ringTools; const CaloSubdetectorGeometry * endcapGeometry_; const CaloSubdetectorGeometry * barrelGeometry_; - const EcalElectronicsMapping * elecMap_; + const EcalElectronicsMapping * elecMap_; unsigned int _iter; }; diff --git a/EcalTiming/python/PlotUtils.py b/EcalTiming/python/PlotUtils.py index 673960f..5dce24d 100644 --- a/EcalTiming/python/PlotUtils.py +++ b/EcalTiming/python/PlotUtils.py @@ -139,3 +139,15 @@ def drawMultipleSame(hists,labels,filename,colors=[], width = 500, height = 500, leg.Draw() canv.SaveAs(filename) + +def makePlotFolder(path): + import os + import errno + import shutil + try: + os.makedirs(path) + except OSError as exc: # Python >2.5 + if exc.errno == errno.EEXIST and os.path.isdir(path): + pass + else: raise + shutil.copy("plots/index.php", path) diff --git a/EcalTiming/python/ecalLocalRecoSequenceAlCaP0Stream_cff.py b/EcalTiming/python/ecalLocalRecoSequenceAlCaP0Stream_cff.py new file mode 100644 index 0000000..9d82aa1 --- /dev/null +++ b/EcalTiming/python/ecalLocalRecoSequenceAlCaP0Stream_cff.py @@ -0,0 +1,17 @@ +from RecoLocalCalo.Configuration.ecalLocalRecoSequence_cff import * + +ecalLocalRecoSequenceAlCaP0Stream = cms.Sequence (ecalMultiFitUncalibRecHit * + ecalRecHit) + +ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag("dummyHits","dummyBarrelDigis") +ecalMultiFitUncalibRecHit.EEdigiCollection = cms.InputTag("dummyHits","dummyEndcapDigis") + + +#ecalDetIdToBeRecovered = RecoLocalCalo.EcalRecProducers.ecalDetIdToBeRecovered_cfi.ecalDetIdToBeRecovered.clone() +ecalRecHit.killDeadChannels = cms.bool( False ) +ecalRecHit.recoverEBVFE = cms.bool( False ) +ecalRecHit.recoverEEVFE = cms.bool( False ) +ecalRecHit.recoverEBFE = cms.bool( False ) +ecalRecHit.recoverEEFE = cms.bool( False ) +ecalRecHit.recoverEEIsolatedChannels = cms.bool( False ) +ecalRecHit.recoverEBIsolatedChannels = cms.bool( False ) diff --git a/EcalTiming/python/loadOldCalib.py b/EcalTiming/python/loadOldCalib.py index 3541637..a22e3ac 100644 --- a/EcalTiming/python/loadOldCalib.py +++ b/EcalTiming/python/loadOldCalib.py @@ -1,12 +1,13 @@ import ROOT from EcalTiming.EcalTiming.txt2tree import txt2tree from collections import namedtuple +import math Crystal = namedtuple('Crystal',("ix", "iy", "iz", "elecID", "FED", "CCU", "iRing")) Calib = namedtuple('Calib',("time", "stddev", "num", "energy","rawid" )) -det_name = {0:"EB", -1:"EEM", 1:"EEP"} +det_name = {0:"EB", -1:"EEM", 1:"EEP", (0,-1):"EBM", (0,1):"EBP"} def getCalibFromTree(tree): calibMap = dict() @@ -136,6 +137,37 @@ def calcMeanBySD(map): return dict( [ (iz, sum[iz]/num[iz]) for iz in [0,-1,1] ] ) +def calcMeanByCCU(map): + global rawidMap + CCU_sum = dict() + CCU_num = dict() + for id in map: + key =(rawidMap[id].FED, rawidMap[id].CCU) + CCU_sum[key] = 0 + CCU_num[key] = 0 + + for id in map: + key =(rawidMap[id].FED, rawidMap[id].CCU) + CCU_sum[key] += map[id] + CCU_num[key] += 1 + + out_map = dict() + for id in map: + key =(rawidMap[id].FED, rawidMap[id].CCU) + out_map[id] = CCU_sum[key]/CCU_num[key] + + printed = set() + for id,time in out_map.items(): + if abs(time) > 1.: + c = rawidMap[id] + k = (c.FED, c.CCU) + if k not in printed: + print c.FED, c.CCU, c.ix, c.iy, c.iz, time + printed.add(k) + return out_map + + + def plot1d(map, dir, name, low, hi,xtitle="[ns]", ytitle="Events"): global rawidMap @@ -203,6 +235,66 @@ def plot2d(map, dir, name, low, hi,setLogz=False): h.Draw("colz") c.SaveAs(dir + '/' + h.GetName() + ".png") +def calcMeanByiRing(m,iRings): + global rawidMap + time_sum = dict(zip(iRings,[0]*len(iRings))) + time_sum2 = dict(zip(iRings,[0]*len(iRings))) + time_num = dict(zip(iRings,[0]*len(iRings))) + + for id,time in m.iteritems(): + try: + crys = rawidMap[id] + except KeyError: + print id, "not found" + continue + key = (crys.iz, crys.iRing) + if key in iRings: + time_sum[key] += time + time_sum2[key] += time*time + time_num[key] += 1 + + mean = dict() + stddev = dict() + for k in time_sum: + if time_num[k]: + mean[k] = time_sum[k]/time_num[k] + stddev[k] = math.sqrt(time_sum2[k]/time_num[k] - mean[k]**2) + return mean,stddev + +def calcMean(m): + global rawidMap + time_sum = dict(zip(det_name.keys(),[0]*len(det_name))) + time_sum2 = dict(zip(det_name.keys(),[0]*len(det_name))) + time_num = dict(zip(det_name.keys(),[0]*len(det_name))) + + for id,time in m.iteritems(): + try: + crys = rawidMap[id] + except KeyError: + print id, "not found" + continue + if crys.iz == 0: + side = math.copysign(1,crys.ix) + time_sum[(0,side)] += time + time_sum2[(0,side)] += time*time + time_num[(0,side)] += 1 + time_sum[0] += time + time_sum2[0] += time*time + time_num[0] += 1 + else: + time_sum[crys.iz] += time + time_sum2[crys.iz] += time*time + time_num[crys.iz] += 1 + + mean = dict() + stddev = dict() + for k in time_sum: + if time_num[k]: + mean[k] = time_sum[k]/time_num[k] + stddev[k] = math.sqrt(time_sum2[k]/time_num[k] - mean[k]**2) + return mean,stddev + + if __name__ == "__main__": print "hi" calib = getCalib(); diff --git a/EcalTiming/python/methDiff.py b/EcalTiming/python/methDiff.py index ff0b9bb..737ef08 100644 --- a/EcalTiming/python/methDiff.py +++ b/EcalTiming/python/methDiff.py @@ -1,9 +1,12 @@ import ROOT import EcalTiming.EcalTiming.loadOldCalib as cal -from EcalTiming.EcalTiming.PlotUtils import customROOTstyle +from EcalTiming.EcalTiming.PlotUtils import customROOTstyle, makePlotFolder import sys,os import shutil import errno +import numpy as np + +colors = [ROOT.kBlack, ROOT.kBlue, ROOT.kRed, ROOT.kGreen, ROOT.kCyan, ROOT.kMagenta, ROOT.kOrange - 3, ROOT.kYellow -2, ROOT.kGreen - 3, ROOT.kCyan - 3] def mapDiff(outdir, map1, map2, name1, name2, errs=[]): @@ -13,19 +16,213 @@ def mapDiff(outdir, map1, map2, name1, name2, errs=[]): if errs: resMap = cal.divideCalib(map1, errs) cal.plot1d(resMap,outdir, "resmap" + name2 + "_" + name1, -2, 2, xtitle="diff/timeError") - cal.plot1d(diffMap, outdir, "diff" + name2 + "_" + name1, -.4, .4) + cal.plot1d(diffMap, outdir, "diff" + name2 + "_" + name1, -2, 2) ROOT.gStyle.SetOptStat(0) - cal.plot2d(diffMap, outdir, "diffMap" + name2 + "_" + name1, -.4, .4) + cal.plot2d(diffMap, outdir, "diffMap" + name2 + "_" + name1, -1, 1) ringdiff = cal.plotiRing(diffMap, outdir, "iRing" + name2 + "_" + name1) return ringdiff +def plotAvePerMap(hist_name, maps, names): + + nmaps = len(maps) + if nmaps != len(names): + print "LENGTH MISMATCH" + return + + bfield_box = True + Bon_start = ["251244", "254790", "256630"] + Bon_end = ["254232", "254914", "256869"] + aves = dict() + #aves[0] = ("EB_ave_" + hist_name, "EB", np.zeros(nmaps)) + aves[(0,-1)] = ("EBM_ave_" + hist_name, "EB-", np.zeros(nmaps)) + aves[(0,1)] = ("EBP_ave_" + hist_name, "EB+", np.zeros(nmaps)) + aves[-1] = ("EEM_ave_" + hist_name, "EE-", np.zeros(nmaps)) + aves[1] = ("EEP_ave_" + hist_name, "EE+", np.zeros(nmaps)) + sds = dict() + #sds[0] = ("EB_stddev_" + hist_name, "EB", np.zeros(nmaps)) + sds[(0,-1)] = ("EBM_stddev_" + hist_name, "EB-", np.zeros(nmaps)) + sds[(0,1)] = ("EBP_stddev_" + hist_name, "EB+", np.zeros(nmaps)) + sds[-1] = ("EEM_stddev_" + hist_name, "EE-", np.zeros(nmaps)) + sds[1] = ("EEP_stddev_" + hist_name, "EE+", np.zeros(nmaps)) + + for i in range(nmaps): + ave,stddev = cal.calcMean(maps[i]) + for k in aves: + if k in ave: + aves[k][2][i] = ave[k] + sds[k][2][i] = stddev[k] + + + x1 = ROOT.gStyle.GetPadLeftMargin(); + x2 = 1 - ROOT.gStyle.GetPadRightMargin(); + y2 = 1 - ROOT.gStyle.GetPadTopMargin(); + y1 = y2 - .1 + leg = ROOT.TLegend(x1,y1,x2,y2) + leg.SetNColumns(4 + bfield_box) + c = ROOT.TCanvas("c","c",1600,1200) + h = ROOT.TH1F("test","",1000,-.1,1.1) + h.SetLineColor(0) + h.SetFillColor(19) + if bfield_box: + leg.AddEntry(h,"Bon","f") + x = np.linspace(0,1,nmaps) + h.GetXaxis().SetTitle("") + h.GetYaxis().SetTitle("Average Time[ns]") + binnames = ['251244', '254231', '256630','256869'] + x1 = [0]*len(Bon_start) + x2 = [0]*len(Bon_start) + for i in range(nmaps): + bin = h.FindBin(x[i]) + if names[i] in binnames: + h.GetXaxis().SetBinLabel(bin,names[i]) + if names[i] in Bon_start: + x1[Bon_start.index(names[i])] = x[i] + if names[i] in Bon_end: + x2[Bon_end.index(names[i])] = x[i] + + h.SetTitle("") + h.SetMinimum(-1) + h.SetMaximum(1) + h.Draw() + if bfield_box: + boxes = [ ROOT.TBox(x1[bi],-1,x2[bi],1) for bi in range(len(x1))] + [b.Draw() for b in boxes] + ci = 0 + mg_ave = ROOT.TMultiGraph() + for name,leg_name,ave in aves.values(): + g = ROOT.TGraph(nmaps,x,ave) + g.SetLineColor(colors[ci % len(colors)]) + g.SetLineWidth(2) + mg_ave.Add(g) + leg.AddEntry(g,leg_name,"lp") + ci += 1 + leg.Draw() + mg_ave.Draw("LP") + c.SaveAs(outdir + "/ave_" + hist_name + ".png") + + h.GetYaxis().SetTitle("Standard Deviation [ns]") + h.SetMinimum(0) + h.SetMaximum(2) + h.SetTitle("") + h.Draw() + if bfield_box: + boxes = [ ROOT.TBox(x1[bi],0,x2[bi],2) for bi in range(len(x1))] + [b.Draw() for b in boxes] + ci = 0 + mg_stddev = ROOT.TMultiGraph() + for name,__,stddev in sds.values(): + g = ROOT.TGraph(nmaps,x,stddev) + g.SetLineColor(colors[ci % len(colors)]) + g.SetLineWidth(2) + mg_stddev.Add(g) + ci += 1 + leg.Draw() + mg_stddev.Draw("LP") + c.SaveAs(outdir + "/stddev_" + hist_name + ".png") + +def drawBoxes(x1, x2, y1,y2): + print x1,x2, y1, y2 + boxes = [ ROOT.TBox(x1[bi],y1,x2[bi],y2) for bi in range(len(x1))] + [b.Draw() for b in boxes] + testbox = ROOT.TBox(.4, 0, .5, .1) + testbox.Draw() + +def plotRingAverPerMap(hist_name, maps, names): + + nmaps = len(maps) + if nmaps != len(names): + print "LENGTH MISMATCH" + return + + bfield_box = True + Bon_start = ["251244", "254790", "256630"] + Bon_end = ["254232", "254914", "256869"] + + iRingsEB = [ (0,1), (0,80), (0,-1), (0,-80)] + iRingsEEM= [(-1, 5), (-1,15), (-1,25), (-1,35)] + iRingsEEP= [(1, 5), (1,15), (1,25), (1,35)] + + iRings = iRingsEB + iRingsEEM + iRingsEEP + aves = dict() + sds = dict() + for key in iRings: + aves[key] = ( "(%d,%d)" % key, np.zeros(nmaps)) + sds[key] = ( "(%d,%d)" % key, np.zeros(nmaps)) + + for i in range(nmaps): + ave,stddev = cal.calcMeanByiRing(maps[i], iRings) + for k in aves: + if k in ave: + aves[k][1][i] = ave[k] + sds[k][1][i] = stddev[k] + + + h = ROOT.TH1F("test","",1000,-.1,1.1) + h.SetLineColor(0) + h.SetFillColor(19) + x = np.linspace(0,1,nmaps) + h.GetXaxis().SetTitle("") + h.GetYaxis().SetTitle("Average Time[ns]") + + binnames = ['251244', '254231', '256630','256869'] + box_x1 = [0]*len(Bon_start) + box_x2 = [0]*len(Bon_start) + for i in range(nmaps): + bin = h.FindBin(x[i]) + if names[i] in binnames: + h.GetXaxis().SetBinLabel(bin,names[i]) + if names[i] in Bon_start: + box_x1[Bon_start.index(names[i])] = x[i] + if names[i] in Bon_end: + box_x2[Bon_end.index(names[i])] = x[i] + + plot(h, x, aves, iRingsEB, "ave_EB_iRing" , -1, 1, box_x1, box_x2) + plot(h, x, aves, iRingsEEM, "ave_EEM_iRing", -1, 1, box_x1, box_x2) + plot(h, x, aves, iRingsEEP, "ave_EEP_iRing", -1, 1, box_x1, box_x2) + + h.GetYaxis().SetTitle("Standard Deviation [ns]") + plot(h, x, sds, iRingsEB , "stddev_EB_iRing", 0, 2, box_x1, box_x2) + plot(h, x, sds, iRingsEEM, "stddev_EEM_iRing", 0, 2, box_x1, box_x2) + plot(h, x, sds, iRingsEEP, "stddev_EEP_iRing", 0, 2, box_x1, box_x2) + +def plot(h, x, graphs, keys, name, min_y, max_y, box_x1, box_x2): + c = ROOT.TCanvas("c","c",1600,1200) + h.SetTitle("") + h.SetMinimum(min_y) + h.SetMaximum(max_y) + h.Draw() + + x1 = ROOT.gStyle.GetPadLeftMargin(); + x2 = 1 - ROOT.gStyle.GetPadRightMargin(); + y2 = 1 - ROOT.gStyle.GetPadTopMargin(); + y1 = y2 - .1 + + leg = ROOT.TLegend(x1,y1,x2,y2) + leg.SetNColumns(len(keys) + 1) + + #added boxes + leg.AddEntry(h,"Bon","f") + drawBoxes(box_x1,box_x2,min_y,max_y) + ci = 0 + mg_ave = ROOT.TMultiGraph() + for key in keys: + leg_name,ave = graphs[key] + g = ROOT.TGraph(len(x),x,ave) + g.SetLineColor(colors[ci % len(colors)]) + g.SetLineWidth(2) + mg_ave.Add(g) + leg.AddEntry(g,leg_name,"lp") + ci += 1 + leg.Draw() + mg_ave.Draw("LP") + c.SaveAs(outdir + "/" + name + ".png") + def listDiff(hist_name, maps1, maps2, names1, names2, errs=[]): mg = dict() mg[0] = ("EB_" + hist_name,ROOT.TMultiGraph()) mg[-1] = ("EEM_" + hist_name,ROOT.TMultiGraph()) mg[1] = ("EEP_" + hist_name,ROOT.TMultiGraph()) - colors = [ROOT.kBlack, ROOT.kBlue, ROOT.kRed, ROOT.kGreen, ROOT.kCyan, ROOT.kMagenta] x1 = ROOT.gStyle.GetPadLeftMargin(); x2 = 1 - ROOT.gStyle.GetPadRightMargin(); @@ -61,20 +258,12 @@ def listDiff(hist_name, maps1, maps2, names1, names2, errs=[]): dir, basename = os.path.split(file_pattern) dir = dir.split('/') - print dir[-2:] - outdir = '/'.join(dir[:-2]) + "/plots/" + "/" + sys.argv[2] + "/" + print dir[-1:] + outdir = '/'.join(dir[:-1]) + "/plots/" + "/" + sys.argv[2] + "/" outdir = os.path.normpath(outdir) - def mkdir_p(path): - try: - os.makedirs(path) - except OSError as exc: # Python >2.5 - if exc.errno == errno.EEXIST and os.path.isdir(path): - pass - else: raise - - mkdir_p(outdir) - shutil.copy("plots/index.php", outdir) + print outdir + makePlotFolder(outdir) subs = sys.argv[3:] files = [ file_pattern % s for s in subs ] @@ -83,6 +272,10 @@ def mkdir_p(path): print files print outdir + #do diff each step + plotRingAverPerMap("PerRuniRing", maps, subs) + plotAvePerMap("PerRun", maps, subs) #listDiff( sys.argv[2], maps[:-1], maps[1:], subs[:-1], subs[1:], errs=maps_err[:-1]) - listDiff(sys.argv[2], maps[:1]*(len(maps)-1), maps[1:], subs[:1]*(len(maps) -1), subs[1:], errs=maps_err[:-1]) + #do diff with respect to first + #listDiff(sys.argv[2] + "_base", maps[:1]*(len(maps)-1), maps[1:], subs[:1]*(len(maps) -1), subs[1:], errs=maps_err[:-1]) diff --git a/EcalTiming/python/plotMaps.py b/EcalTiming/python/plotMaps.py index e9bb4a8..f9cc62e 100644 --- a/EcalTiming/python/plotMaps.py +++ b/EcalTiming/python/plotMaps.py @@ -102,6 +102,8 @@ def addFitToPlot(h,outfilename): def plotMaps(tree, outdir, prefix=""): + low = -1 + hi = 1 # dictionaries to store histograms time = dict() time_rel2012 = dict() @@ -154,7 +156,7 @@ def plotMaps(tree, outdir, prefix=""): initHists(prefix,timeError, initMap, key, "timeError", "Time Error[ns]") initHists(prefix,timeError1d, inittime1d, key, "timeError1d", "Time Error Test [ns]", low=0, hi=.2) initHists(prefix,stdDev, initMap, key, "stdDev", "Std Dev [ns]") - initHists(prefix,time1d, inittime1d, key, "time1d", "time1d") + initHists(prefix,time1d, inittime1d, key, "time1d", "time1d", low=low, hi=hi) initHists(prefix,occupancy, initMap, key, "occupancy", "Occupancy") initHists(prefix,energy, initMap, key, "energy", "Energy [GeV]") @@ -211,7 +213,7 @@ def plotMaps(tree, outdir, prefix=""): c = ROOT.TCanvas("c","c",1600,1200) for key in CCU_maps: - CCU_maps[key].SetAxisRange(-2, 2, "Z") + CCU_maps[key].SetAxisRange(low,hi, "Z") CCU_maps[key].SetZTitle("[ns]") CCU_maps[key].Draw("colz") c.SaveAs(outdir + "/" + CCU_maps[key].GetName() + ".png") @@ -232,9 +234,9 @@ def plotMaps(tree, outdir, prefix=""): time[key].SetAxisRange(-5,5,"Z") time[key].Draw("colz") c.SaveAs(outdir + "/" + time[key].GetName() + ".5.png") - time[key].SetAxisRange(-2,2,"Z") + time[key].SetAxisRange(low,hi,"Z") time[key].Draw("colz") - c.SaveAs(outdir + "/" + time[key].GetName() + ".2.png") + c.SaveAs(outdir + "/" + time[key].GetName() + ".png") c = ROOT.TCanvas("c","c",1600,1200) for key in time_rel2012: @@ -340,9 +342,9 @@ def mkdir_p(path): else: raise mkdir_p(outdir) - shutil.copy("plots/index.php", outdir) + #shutil.copy("plots/index.php", outdir) file = ROOT.TFile.Open(filename) - tree = file.Get("filter/EcalSplashTiming/timingTree") + tree = file.Get("timing/EcalSplashTiming/timingTree") time = plotMaps(tree, outdir, prefix = prefix) diff --git a/EcalTiming/python/runInfo.py b/EcalTiming/python/runInfo.py new file mode 100644 index 0000000..07dcfc1 --- /dev/null +++ b/EcalTiming/python/runInfo.py @@ -0,0 +1,70 @@ +from collections import namedtuple +from FWCore.PythonUtilities.LumiList import LumiList + +jsonFolder="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions15/13TeV/" + +RunPeriod = namedtuple('RunPeriod',("name","bfield","firstRun","lastRun","json","bunchSpacing","GT")) +runPeriods = [] + +Run2015AB_json = jsonFolder + "Cert_246908-255031_13TeV_PromptReco_Collisions15_50ns_JSON_v2.txt" +Run2015AB_ZeroTesla_json = jsonFolder + "Cert_246908-252126_13TeV_PromptReco_Collisions15_ZeroTesla_50ns_JSON.txt" +Run2015CD_json = jsonFolder + "Cert_246908-256869_13TeV_PromptReco_Collisions15_25ns_JSON.txt" +Run2015CD_ZeroTesla_json = jsonFolder + "Cert_246908-256869_13TeV_PromptReco_Collisions15_ZeroTesla_25ns_JSON.txt" +Run254833_json = jsonFolder + "Cert_254833_13TeV_PromptReco_Collisions15_JSON.txt" + +#Run2015AB = LumiList(filename = Run2015AB_json ) +#Run2015AB_ZeroTesla = LumiList(filename = Run2015AB_ZeroTesla_json) +#Run2015CD = LumiList(filename = Run2015CD_json ) +#Run2015CD_ZeroTesla = LumiList(filename = Run2015CD_ZeroTesla_json) +#Run254833 = LumiList(filename = Run254833_json ) + +runPeriods.append(RunPeriod("Run2015A-v1", False, 247607, 250932, Run2015AB_ZeroTesla_json, 50, "GR_P_V56")) +runPeriods.append(RunPeriod("Run2015A-v1", True, 247607, 250932, Run2015AB_json, 50, "GR_P_V56")) +runPeriods.append(RunPeriod("Run2015B-v1", False, 250985, 253620, Run2015AB_ZeroTesla_json, 50, "74X_dataRun2_Prompt_v0")) +runPeriods.append(RunPeriod("Run2015B-v1", True, 250985, 253620, Run2015AB_json, 50, "74X_dataRun2_Prompt_v0")) +runPeriods.append(RunPeriod("Run2015C-v1", False, 253659, 256464, Run2015CD_ZeroTesla_json, 25, "74X_dataRun2_Prompt_v1")) +runPeriods.append(RunPeriod("Run2015C-v1", True, 253659, 256464, Run2015CD_json, 25, "74X_dataRun2_Prompt_v1")) +runPeriods.append(RunPeriod("Run2015D-v1", False, 256630, 999999, Run2015CD_ZeroTesla_json, 25, "74X_dataRun2_Prompt_v2")) +runPeriods.append(RunPeriod("Run2015D-v1", True, 256630, 999999, Run2015CD_json, 25, "74X_dataRun2_Prompt_v2")) + +runPeriods.append(RunPeriod("Run2015C-v1", True, 254833, 254833, Run254833_json, 50, "74X_dataRun2_Prompt_v1")) + +def LumiListForRunPeriod(rp): + ll = LumiList(filename = rp.json) + runs = [ run for run in map(int,ll.getRuns()) if run >= rp.firstRun and run <= rp.lastRun] + ll.selectRuns(runs) + return ll + +def getRuns(name=None,bfield=None,bunchSpacing=None): + ll = LumiList() + for rp in runPeriods: + if name is None or rp.name == name: + if bfield is None or rp.bfield == bfield: + if bunchSpacing is None or rp.bunchSpacing == bunchSpacing: + newll = LumiListForRunPeriod(rp) + ll += LumiListForRunPeriod(rp) + return ll.getRuns() + +def getRunPeriod(name,bfield,bunchSpacing): + for rp in runPeriods: + if rp.name == name: + if rp.bfield == bfield: + if rp.bunchSpacing == bunchSpacing: + return rp + +import sys + +if __name__ == "__main__": + name = sys.argv[1] + bfield = bool(int(sys.argv[2])) + bunchSpacing = int(sys.argv[3]) + rp = getRunPeriod(name,bfield,bunchSpacing) + ll = LumiListForRunPeriod(rp) + if not ll.getRuns(): + sys.exit(0) + options = "" + options += " -r=" +",".join(map(str,ll.getRuns())) + options += " --json=" + rp.json + options += " --tag=" + rp.GT + options += " --runperiod=" + rp.name + print options diff --git a/EcalTiming/python/storeTools_cff.py b/EcalTiming/python/storeTools_cff.py new file mode 100644 index 0000000..097280b --- /dev/null +++ b/EcalTiming/python/storeTools_cff.py @@ -0,0 +1,196 @@ +import FWCore.ParameterSet.Config as cms +import os,sys +import getopt +import commands +import re + +def natural_sort(l): + convert = lambda text: int(text) if text.isdigit() else text.lower() + alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ] + return sorted(l, key = alphanum_key) + +""" +lists the files available in castor +""" +def fillFromStore(dir,ffile=0,step=-1,generatePfn=True): + + localdataset=cms.untracked.vstring() + if(len(dir)==0) : return localdataset + + #check if it is a directory (check if it is castor, eos or local) + prefix='singlefile' + lsout=[dir] + if(dir.find('path=')>=0) : + print 'Using dbs to query %s'%(dir) + prefix='eoscms' + lsout=commands.getstatusoutput('dbs lsf --' + dir)[1].split() + + elif(dir.find('castor')>=0) : + prefix='rfio' + lscommand ='rfdir ' + dir + ' | awk \'{print $9}\'' + lsout = commands.getstatusoutput(lscommand)[1].split() + + elif(dir.find('/store/')==0): + prefix='eoscms' + lscommand = 'cmsLs ' + dir + ' | grep root | awk \'{print $5}\'' + lsouttmp = commands.getstatusoutput(lscommand)[1].split() + + #this is needed as cmsLs lists twice files staged from castor + #(only needed during transition to EOS, can remove now) + lsout=[] + nduplicate=0 + for l in lsouttmp : + if len(l)==0 : continue + if l in lsout : + nduplicate += 1 + continue + #filter out CMG trees and histograms + basename = os.path.basename(l) + if(basename.find('tree_')==0) : continue + if(basename.find('histogram')==0): continue + lsout.append(l) + print 'Discarded ' + str(nduplicate) + ' files duplicated in cmsLs output' + + elif(dir.find('.root')<0): + prefix='file' + lscommand='ls ' + dir + lsout = commands.getstatusoutput(lscommand)[1].split() + + #check for the files needed (first file, firstfile+step) + ifile=0 + for line in lsout : + if(type(line) is not str ) : continue + if(len(line)==0) : continue + if(line.find('root')<0) : continue + + if(ifile>=ffile): + if( (step<0) or (step>0 and ifile=0 or command_out[1].find("probably not closed")>=0 or command_out[1].find("Corrupted")>=0):return False + return True + + +""" +check store for duplicates +""" +def checkStoreForDuplicates(outdir): + ls_cms = "ls " + outdir + isEOS=False + isCastor=False + if(outdir.find('/store/')==0) : + isEOS=True + splitOnString=',' + ls_cms='cmsLs ' + outdir + ' | grep root | awk \'{print $5}\'' + elif(outdir.find('/store/')==0) : + isCastor = True + ls_cms = "rfdir " + outdir + + ' | grep root' + nOutFile = 0 + lsCmd_out = commands.getstatusoutput(ls_cms) + jobNumbers = [] + duplicatedJobs = [] + origFiles=[] + duplicatedFiles=[] + if lsCmd_out[0] == 0: + lsCmd_outLines = lsCmd_out[1].split("\n") + if len(lsCmd_outLines) != 0: + for fileLine in lsCmd_outLines: + if not "root" in fileLine: continue + fileName=fileLine + if(isCastor) : fileName = fileLine.split()[8] + + if(checkInputFile(fileName)==True): + jobNumber=-1 + try: + fileBaseName=os.path.basename(fileName) + jobNumber=int(fileBaseName.split("_")[1]) + except: + continue + + if jobNumber in jobNumbers: + if not jobNumber in duplicatedJobs: duplicatedJobs.append(jobNumber) + duplicatedFiles.append(fileName) + else : + jobNumbers.append(jobNumber) + origFiles.append(fileName) + nOutFile += 1 + else: + print(" #corrupted file found : " + fileName) + duplicatedFiles.append(fileName) + return natural_sort(duplicatedFiles) + + +""" +clean up for duplicats in the storage area +""" +def removeDuplicates(dir): + duplicatedFiles=checkStoreForDuplicates(dir) + print 'Removing ' + str(len(duplicatedFiles)) + ' duplicated files in ' + dir + isEOS=False + isCastor=False + if(dir.find('/store/')==0) : isEOS=True + if(dir.find('castor')>=0) : isCastor=True + for f in duplicatedFiles : + print f + if(isEOS) : commands.getstatusoutput('cmsRm ' + f) + elif(isCastor) : commands.getstatusoutput('rfrm ' +dir + '/' + f) + else : commands.getstatusoutput('rm ' +dir + '/' + f) + + +""" +wrapper to read the configuration from comand line +args: castor_directory output_file first_file step +""" +def configureSourceFromCommandLine() : + storeDir='' + outputFile='Events.root' + ffile=0 + step=-1 + try: + if(len(sys.argv)>2 ): + if(sys.argv[2].find('/')>=0 or sys.argv[2].find('.root')>0) : + storeDir=sys.argv[2] + if(len(sys.argv)>3 ): + if(sys.argv[3].find('.root')>0): + outputFile=sys.argv[3] + if(len(sys.argv)>4 ): + if(sys.argv[4].isdigit()) : ffile=int(sys.argv[4]) + if(len(sys.argv)>5 ): + if(sys.argv[5].isdigit()) : step=int(sys.argv[5]) + except: + print '[storeTools_cff] Could not configure from command line, will return default values' + + return outputFile, fillFromStore(storeDir,ffile,step) + + + +def addPrefixSuffixToFileList(Prefix, fileList, Suffix): + outList = [] + for s in fileList: + outList.append(Prefix+s+Suffix) + return natural_sort(outList) + diff --git a/EcalTiming/scripts/testCAF.sh b/EcalTiming/scripts/testCAF.sh index c86a070..37f7f5a 100755 --- a/EcalTiming/scripts/testCAF.sh +++ b/EcalTiming/scripts/testCAF.sh @@ -1,15 +1,25 @@ -RUNLIST="243479 243484 243506" -RUNLIST="248030" -RUNLIST=`cat runlist` -RUNLIST="251244 251251 251252 251521 251522 251548 251559 251560 251561 251562" +#RUNLIST="251244 251251 251252 251521 251522 251548 251559 251560 251561 251562" +#2015B +#RUNLIST="251643 251721 251883" +#RUNPERIOD=Run2015B-v1 +#2015C +#RUNLIST="254231 254232 254790 254879" #254852 +RUNPERIOD=Run2015C-v1 + +#RUNLIST="254292 254293 254294 254319 254332 254341 254342 254349 254458 254459 254608" +RUNLIST="254307 255981 256001 256002 256003 256167 256169 256171 256215 256237" #256245 256349 256350 256355 256406" STREAM=AlCaPhiSym NEVENTS=-1 QUEUE=2nd -DIR=/afs/cern.ch/work/p/phansen/public/EcalTiming/Validation_notCCUSubtracted/ -CONFIG=test/ecalTime_fromAlcaStream_cfg.py - -FROMRECO=NO +EOSPREFIX=root://eoscms//eos/cms/ +EOSDIR=/store/group/dpg_ecal/alca_ecalcalib/EcalTiming/ +#DIR=/afs/cern.ch/work/p/phansen/public/EcalTiming/RunII/ +CONFIG=$PWD/test/ecalTime_fromAlcaStream_cfg.py +NJOBS=1 +USER=$(whoami) + +STEP=RECOTIMEANALYSIS for i in "$@" do case $i in @@ -18,7 +28,7 @@ case $i in shift # past argument=value ;; -c=*|--cfg=*) - CONFIG="${i#*=}" + CONFIG=$PWD/"${i#*=}" shift # past argument=value ;; -d=*|--dir=*) @@ -46,45 +56,80 @@ case $i in SPLIT=YES shift # past argument with no value ;; - --fromreco) - FROMRECO=YES - shift # past argument with no value + --step=*) + STEP="${i#*=}" + shift # past argument=value + ;; + --json=*) + JSON="${i#*=}" + shift # past argument=value + ;; + --runperiod=*) + RUNPERIOD="${i#*=}" + shift # past argument=value + ;; + --tag=*) + GLOBALTAG="${i#*=}" + shift # past argument=value + ;; + --njobs=*) + NJOBS="${i#*=}" + shift # past argument=value ;; *) # unknown option echo option $i not defined + exit 1 ;; esac done +if [[ -z $GLOBALTAG ]] +then + echo Tag not set + exit 1 +fi + for RUN in ${RUNLIST} do - + if ! grep $RUN $JSON > /dev/null; + then + echo $RUN not in JSON $JSON + continue + fi echo "=== RUN = ${RUN}" - OUTDIR=$DIR/${STREAM}-${RUN}/ - mkdir -p ${OUTDIR} + OUTDIR=$EOSDIR/${STREAM}-${RUN}/ + /afs/cern.ch/project/eos/installation/0.3.84-aquamarine/bin/eos.select mkdir ${OUTDIR} + AFSDIR=/afs/cern.ch/work/p/phansen/public/EcalTiming/analysis/${STREAM}-${RUN}/ + mkdir -p ${AFSDIR} - if [ "$FROMRECO" == "NO" ] + if [[ $STEP == *"RECO"* ]] then - STEP=RECOTIMEANALYSIS - filelist=`das_client.py --query="file dataset=/MinimumBias/Commissioning2015-v1/RAW run=${RUN}" --limit=50 | sed '2 d'` - nfiles=`das_client.py --query="file dataset=/${STREAM}/Run2015B-v1/RAW run=${RUN} | count(file.name)" | sed '2 d'` + RECO="_RECO" + nfiles=`das_client.py --query="file dataset=/${STREAM}/${RUNPERIOD}/RAW run=${RUN} | count(file.name)" | sed '2 d'` nfiles=${nfiles:19} if ! [[ $nfiles =~ ^[0-9]+$ ]]; then echo "No Files found" - exit 1 + continue fi #echo Will run over $nfiles files - filelist=`das_client.py --query="file dataset=/${STREAM}/Run2015B-v1/RAW run=${RUN}" --limit=${nfiles} | sed '2 d'` + filelist=`das_client.py --query="file dataset=/${STREAM}/${RUNPERIOD}/RAW run=${RUN}" --limit=${nfiles} | sed '2 d'` # for file in ${filelist} # do # das_client.py --query="file=${file} | sum(file.nevents)" # done else - STEP=TIMEANALYSIS - filelist=`grep ${RUN} ~shervin/public/4peter/fileMap-sorted.dat | cut -d ' ' -f2` + #filelist=`grep ${RUN} ~shervin/public/4peter/fileMap-sorted.dat | cut -d ' ' -f2` + #filelist=(${OUTDIR}/ecalTiming_*_numEvent50000_RECO.root) + #filelist=${filelist[@]/#/file://} + filelist=${OUTDIR} + + fi + if [[ $STEP == *"TIME"* ]] + then + RECO="" fi if [ "$SPLIT" != "YES" ] @@ -93,23 +138,43 @@ do filelist=`echo ${filelist}| sed 's| |,|g;s|,$||'` fi - jsonFile=/afs/cern.ch/cms/CAF/CMSALCA/ALCA_ECALCALIB/json_ecalonly/251022-251562-Prompt-pfgEcal-noLaserProblem.json + #jsonFile=/afs/cern.ch/cms/CAF/CMSALCA/ALCA_ECALCALIB/json_ecalonly/251022-251562-Prompt-pfgEcal-noLaserProblem.json i=0 for file in $filelist do - name=${i}_${en}GeV - runcommand="cmsRun ${CONFIG} files=${filelist} output=${OUTDIR}/ecalTiming_${name}.root maxEvents=${NEVENTS} jsonFile=${jsonFile} minEnergyEB=1.5 minEnergyEE=2.5 step=${STEP}" + for job in `seq $NJOBS` + do + name=${i} + let n=$NEVENTS/$NJOBS + let skip=$n*$i + tmp_file=ecalTiming_${name}.root + if [ "$n" == "-1" ] + then + tmp_file_real=ecalTiming_${name}${RECO}.root + else + tmp_file_real=ecalTiming_${name}_numEvent${n}${RECO}.root + fi + runcommand="cmsRun ${CONFIG} files=${filelist} output=${tmp_file} maxEvents=${n} jsonFile=${JSON} minEnergyEB=1.5 minEnergyEE=2.5 step=${STEP} skipEvents=${skip} globaltag=${GLOBALTAG};" + #runcommand+="ls -l;" + #runcommand+="df -h;" + runcommand+="xrdcp ${tmp_file_real} ${EOSPREFIX}${OUTDIR}/;" + if [[ $STEP == *"TIME"* ]] + then + runcommand+="cp ecalTiming_* ${AFSDIR}/;" + fi if [ "$BATCH" == "YES" ] then - bsub -oo ${OUTDIR}/stdout-${name}-${NEVENTS}.log -eo ${OUTDIR}/stderr-${name}-${NEVENTS}.log -R "rusage[mem=4000]" -q ${QUEUE} "cd $PWD; eval \`scramv1 runtime -sh\`; + bsub -oo log/stdout-${STREAM}-${RUN}-${name}-${NEVENTS}.log -eo log/stderr-${STREAM}-${RUN}-${name}-${NEVENTS}.log -R "rusage[mem=4000]" -q ${QUEUE} "cd $PWD; eval \`scramv1 runtime -sh\`; cd - ${runcommand} " || exit 1 else - $runcommand + echo $runcommand + eval $runcommand fi let i=i+1 done + done done exit 0 diff --git a/EcalTiming/src/EcalCrystalTimingCalibration.cc b/EcalTiming/src/EcalCrystalTimingCalibration.cc index 708539e..cbaae06 100644 --- a/EcalTiming/src/EcalCrystalTimingCalibration.cc +++ b/EcalTiming/src/EcalCrystalTimingCalibration.cc @@ -224,3 +224,4 @@ void EcalCrystalTimingCalibration::dumpToTree(TTree *tree, int ix_, int iy_, int } } + diff --git a/EcalTiming/test/ecalTime_fromAlcaStreamHighPU_cfg.py b/EcalTiming/test/ecalTime_fromAlcaStreamHighPU_cfg.py new file mode 100644 index 0000000..fca5db1 --- /dev/null +++ b/EcalTiming/test/ecalTime_fromAlcaStreamHighPU_cfg.py @@ -0,0 +1,330 @@ +import FWCore.ParameterSet.Config as cms +import os, sys, imp, re +import FWCore.ParameterSet.VarParsing as VarParsing + +#sys.path(".") + +#new options to make everything easier for batch + +############################################################ +### SETUP OPTIONS + +options = VarParsing.VarParsing('standard') +options.register('jsonFile', + "", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "path and name of the json file") +options.register('step', + "RECOTIMEANALYSIS", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "Do reco, time analysis or both, RECO|TIMEANALYSIS|RECOTIMEANALYSIS") +options.register('offset', + 0.0, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "add this to each crystal time") +options.register('minEnergyEB', + 1.5, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "add this to minimum energy threshold") +options.register('minEnergyEE', + 3, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "add this to minimum energy threshold") +options.register('minChi2EB', + 50.0, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "add this to minimum chi2 threshold") +options.register('minChi2EE', + 50.0, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "add this to minimum chi2 threshold") +options.register('isSplash', + 0, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, + "0=false, 1=true" + ) +options.register('streamName', + 'AlCaPhiSym', + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "type of stream: AlCaPhiSym or AlCaP0") +options.register('loneBunch', + 1, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, + "0=No, 1=Yes" + ) + +#options.jsonFile="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions15/13TeV/Cert_246908-256869_13TeV_PromptReco_Collisions15_25ns_JSON.txt" +#options.jsonFile="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions15/13TeV/DCSOnly/json_DCSONLY.txt" +#options.jsonFile; +### setup any defaults you want +options.output="output/ecalTiming.root" +options.secondaryOutput="ntuple.root" + +if(options.streamName=="AlCaP0"): print "stream ",options.streamName#options.files = "/store/data/Commissioning2015/AlCaP0/RAW/v1/000/246/342/00000/048ECF48-F906-E511-95AC-02163E011909.root" +elif(options.streamName=="AlCaPhiSym"): print "stream ",options.streamName#options.files = "/store/data/Commissioning2015/AlCaPhiSym/RAW/v1/000/244/768/00000/A8219906-44FD-E411-8DA9-02163E0121C5.root" +else: + print "stream ",options.streamName," not foreseen" + exit + +#options.files = cms.untracked.vstring +#options.streamName = cms.untracked.vstring +options.maxEvents = -1# -1 means all events +### get and parse the command line arguments +options.parseArguments() +print options + +processname = options.step + +doReco = True +doAnalysis = True +if "RECO" not in processname: + doReco = False +if "TIME" not in processname: + doAnalysis = False + +process = cms.Process(processname) + +process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(1000) + +process.load('Configuration.EventContent.EventContent_cff') +process.load('SimGeneral.MixingModule.mixNoPU_cfi') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_38T_cff') + +if(options.isSplash==1): + ## Get Cosmic Reconstruction + process.load('Configuration/StandardSequences/ReconstructionCosmics_cff') + process.caloCosmics.remove(process.hbhereco) + process.caloCosmics.remove(process.hcalLocalRecoSequence) + process.caloCosmics.remove(process.hfreco) + process.caloCosmics.remove(process.horeco) + process.caloCosmics.remove(process.zdcreco) + process.caloCosmics.remove(process.ecalClusters) + process.caloCosmicOrSplashRECOSequence = cms.Sequence(process.caloCosmics )#+ process.egammaCosmics) +else: + process.load('Configuration/StandardSequences/Reconstruction_cff') + process.recoSequence = cms.Sequence(process.calolocalreco )#+ process.egammaCosmics) + +#process.load('PhiSym.EcalCalibAlgos.ecalPhiSymLocarecoWeights_cff') +#process.load('RecoLocalCalo.Configuration.ecalLocalRecoSequence_cff') +process.load('Configuration.StandardSequences.EndOfProcess_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff') +process.load('EcalTiming.EcalTiming.ecalLocalRecoSequenceAlCaStream_cff') +process.load('EcalTiming.EcalTiming.ecalLocalRecoSequenceAlCaP0Stream_cff') + +if(options.streamName=="AlCaP0"): + process.ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag("hltAlCaPi0EBRechitsToDigis","pi0EBDigis") + process.ecalMultiFitUncalibRecHit.EEdigiCollection = cms.InputTag("hltAlCaPi0EERechitsToDigis","pi0EEDigis") +else: + process.ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag("hltEcalPhiSymFilter","phiSymEcalDigisEB") + process.ecalMultiFitUncalibRecHit.EEdigiCollection = cms.InputTag("hltEcalPhiSymFilter","phiSymEcalDigisEE") + + +## Raw to Digi +process.load('Configuration/StandardSequences/RawToDigi_Data_cff') + +## HLT Filter Splash +import HLTrigger.HLTfilters.hltHighLevel_cfi +process.spashesHltFilter = HLTrigger.HLTfilters.hltHighLevel_cfi.hltHighLevel.clone( + throw = cms.bool(False), + HLTPaths = ['HLT_EG20*', 'HLT_SplashEcalSumET', 'HLT_Calibration','HLT_EcalCalibration','HLT_HcalCalibration','HLT_Random','HLT_Physics','HLT_HcalNZS','HLT_SplashEcalSumET','HLTriggerFinalPath' ] +) + + +## GlobalTag Conditions Related +from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, '74X_dataRun2_Prompt_v2', '') #run2_data', '') + +## Process Digi To Raw Step +process.digiStep = cms.Sequence(process.ecalDigis + process.ecalPreshowerDigis) + +## Process Reco + + + +### Print Out Some Messages +process.MessageLogger = cms.Service("MessageLogger", + cout = cms.untracked.PSet( + threshold = cms.untracked.string('WARNING') + ), + categories = cms.untracked.vstring('ecalTimeTree'), + destinations = cms.untracked.vstring('cout') +) + +# enable the TrigReport and TimeReport +process.options = cms.untracked.PSet( + wantSummary = cms.untracked.bool(True) +# SkipEvent = cms.untracked.vstring('ProductNotFound') +) + +SkipEvent = cms.untracked.vstring('ProductNotFound','EcalProblem') + +# dbs search --query "find file where dataset=/ExpressPhysics/BeamCommissioning09-Express-v2/FEVT and run=124020" | grep store | awk '{printf "\"%s\",\n", $1}' +# Input source +process.source = cms.Source("PoolSource", + secondaryFileNames = cms.untracked.vstring(), + fileNames = cms.untracked.vstring(options.files), +) + +if(len(options.jsonFile) > 0): + import FWCore.PythonUtilities.LumiList as LumiList + process.source.lumisToProcess = LumiList.LumiList(filename = options.jsonFile).getVLuminosityBlockRange() + + + +# Output definition +process.RECOoutput = cms.OutputModule("PoolOutputModule", +splitLevel = cms.untracked.int32(0), +eventAutoFlushCompressedSize = cms.untracked.int32(5242880), +outputCommands = cms.untracked.vstring('drop *',"keep *_ecalRecHitE*Selector_*_*"), +fileName = cms.untracked.string(options.output), +dataset = cms.untracked.PSet( + filterName = cms.untracked.string(''), + dataTier = cms.untracked.string('RECO') +) +) + + +## Histogram files +process.TFileService = cms.Service("TFileService", + fileName = cms.string(options.output), + closeFileFast = cms.untracked.bool(True) + ) + +### NumBer of events +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(options.maxEvents)) + +#DUMMY RECHIT +process.dummyHits = cms.EDProducer("DummyRechitDigis", + doDigi = cms.untracked.bool(True), + # rechits + barrelHitProducer = cms.InputTag('hltAlCaPi0EBUncalibrator','pi0EcalRecHitsEB' ,"HLT"), + endcapHitProducer = cms.InputTag('hltAlCaPi0EEUncalibrator','pi0EcalRecHitsEE' ,"HLT"), + barrelRecHitCollection = cms.untracked.string("dummyBarrelRechitsPi0"), + endcapRecHitCollection = cms.untracked.string("dummyEndcapRechitsPi0"), + # digis + barrelDigis = cms.InputTag('hltAlCaPi0EBRechitsToDigis','pi0EBDigis',"HLT"), + endcapDigis = cms.InputTag('hltAlCaPi0EERechitsToDigis','pi0EEDigis',"HLT"), #changed hltAlCaPi0EERechitsToDigis in LowPU....changed in the file -.- + barrelDigiCollection = cms.untracked.string("dummyBarrelDigisPi0"), + endcapDigiCollection = cms.untracked.string("dummyEndcapDigisPi0")) + +##ADDED +# TRIGGER RESULTS FILTER +process.triggerSelectionLoneBunch = cms.EDFilter( "TriggerResultsFilter", + triggerConditions = cms.vstring('L1_AlwaysTrue'), + hltResults = cms.InputTag( "TriggerResults", "", "HLT" ), + l1tResults = cms.InputTag( "hltGtDigis" ), + l1tIgnoreMask = cms.bool( False ), + l1techIgnorePrescales = cms.bool( False ), + daqPartitions = cms.uint32( 1 ), + throw = cms.bool( True ) + ) + +process.filter=cms.Sequence() +if(options.isSplash==1): + process.filter+=process.spashesHltFilter + process.reco_step = cms.Sequence(process.caloCosmicOrSplashRECOSequence) +else: + if(options.streamName=="AlCaP0"): + #from RecoLocalCalo.Configuration.ecalLocalRecoSequence_cff import * + ##process.reco_step = cms.Sequence(ecalMultiFitUncalibRecHit * + ## ecalRecHit) + #process.ecalMultiFitUncalibRecHit = RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi.ecalMultiFitUncalibRecHit.clone() + + #process.ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag('dummyHits','dummyBarrelDigis')#,'piZeroAnalysis') + #process.ecalMultiFitUncalibRecHit.EEdigiCollection = cms.InputTag('dummyHits','dummyEndcapDigis')#,'piZeroAnalysis') + #ecalRecHit.killDeadChannels = False + #ecalRecHit.recoverEBFE = False + if(options.loneBunch==1): + process.filter+=process.triggerSelectionLoneBunch + import RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi + process.ecalMultiFitUncalibRecHit = RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi.ecalMultiFitUncalibRecHit.clone() + process.ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag('dummyHits','dummyBarrelDigisPi0')#,'piZeroAnalysis') + process.ecalMultiFitUncalibRecHit.EEdigiCollection = cms.InputTag('dummyHits','dummyEndcapDigisPi0')#,'piZeroAnalysis') + + #UNCALIB to CALIB + from RecoLocalCalo.EcalRecProducers.ecalRecHit_cfi import * + process.ecalDetIdToBeRecovered = RecoLocalCalo.EcalRecProducers.ecalDetIdToBeRecovered_cfi.ecalDetIdToBeRecovered.clone() + process.ecalRecHit.killDeadChannels = cms.bool( False ) + process.ecalRecHit.recoverEBVFE = cms.bool( False ) + process.ecalRecHit.recoverEEVFE = cms.bool( False ) + process.ecalRecHit.recoverEBFE = cms.bool( False ) + process.ecalRecHit.recoverEEFE = cms.bool( False ) + process.ecalRecHit.recoverEEIsolatedChannels = cms.bool( False ) + process.ecalRecHit.recoverEBIsolatedChannels = cms.bool( False ) + + process.reco_step = cms.Sequence(process.dummyHits + * process.ecalMultiFitUncalibRecHit + * process.ecalRecHit) + else: + #process.reco_step = cms.Sequence(process.reconstruction_step_multiFit) + if(options.loneBunch==1): + process.filter+=process.triggerSelectionLoneBunch + process.reco_step = cms.Sequence(process.ecalLocalRecoSequenceAlCaStream) + + +### Process Full Path +if(options.isSplash==0): + process.digiStep = cms.Sequence() + + + +evtPlots = True if options.isSplash else False + + +#import Electronics mapping +process.load("Geometry.EcalCommonData.EcalOnly_cfi") +process.load("Geometry.EcalMapping.EcalMapping_cfi") +process.load("Geometry.EcalMapping.EcalMappingRecord_cfi") +#process.load("Geometry.CaloEventSetup.CaloGeometry_cff") + +#ESLooperProducer looper is imported here: +process.load('EcalTiming.EcalTiming.ecalTimingCalibProducer_cfi') +process.load('EcalTiming.EcalTiming.RecHitsSelector_cfi') + +#process.timing.outputDumpFile = process.TFileService.fileName #spostato per vedere se l'errore cambia o rimane o addirittura sparisce! +process.timing.recHitEBCollection = cms.InputTag("ecalRecHitEBSelector") +process.timing.recHitEECollection = cms.InputTag("ecalRecHitEESelector") +process.timing.isSplash= cms.bool(True if options.isSplash else False) +process.timing.makeEventPlots=evtPlots +process.timing.globalOffset = cms.double(options.offset) +process.timing.outputDumpFile = process.TFileService.fileName +process.timing.energyThresholdOffsetEB = cms.double(options.minEnergyEB) +process.timing.energyThresholdOffsetEE = cms.double(options.minEnergyEE) +process.timing.chi2ThresholdOffsetEB = cms.double(options.minChi2EB) +process.timing.chi2ThresholdOffsetEE = cms.double(options.minChi2EE) +process.timing.storeEvents = cms.bool(True) + + +process.analysis = cms.Sequence( process.timing ) +process.reco = cms.Sequence( (process.filter + + process.digiStep + + process.reco_step) + * (process.ecalRecHitEBSelector + process.ecalRecHitEESelector) + ) + + +process.seq = cms.Sequence() +if doReco: + process.seq += process.reco +if doAnalysis: + process.seq += process.analysis +else: + process.endp = cms.EndPath(process.RECOoutput) + +process.p = cms.Path(process.seq) + +processDumpFile = open('processDump.py', 'w') +print >> processDumpFile, process.dumpPython() diff --git a/EcalTiming/test/ecalTime_fromAlcaStream_cfg.py b/EcalTiming/test/ecalTime_fromAlcaStream_cfg.py index eec75d5..e5aca26 100644 --- a/EcalTiming/test/ecalTime_fromAlcaStream_cfg.py +++ b/EcalTiming/test/ecalTime_fromAlcaStream_cfg.py @@ -20,21 +20,36 @@ VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string, "Do reco, time analysis or both, RECO|TIMEANALYSIS|RECOTIMEANALYSIS") +options.register('skipEvents', + 0, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, + "Skip this many events") options.register('offset', 0.0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.float, "add this to each crystal time") options.register('minEnergyEB', - 0.0, + 1.5, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.float, "add this to minimum energy threshold") options.register('minEnergyEE', - 0.0, + 3.0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.float, "add this to minimum energy threshold") +options.register('minChi2EB', + 50.0, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "add this to minimum chi2 threshold") +options.register('minChi2EE', + 50.0, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float, + "add this to minimum chi2 threshold") options.register('isSplash', 0, VarParsing.VarParsing.multiplicity.singleton, @@ -46,35 +61,61 @@ VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string, "type of stream: AlCaPhiSym or AlCaP0") +options.register('loneBunch', + 1, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.int, + "0=No, 1=Yes" + ) +options.register('globaltag', + '', + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "Global tag to use, no default") +#options.jsonFile ="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions15/13TeV/Cert_246908-256869_13TeV_PromptReco_Collisions15_25ns_JSON.txt" +#options.jsonFile ="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions15/13TeV/Cert_254986-255031_13TeV_PromptReco_Collisions15_LOWPU_25ns_JSON_MuonPhys.txt" +#options.jsonFile + + ### setup any defaults you want options.output="output/ecalTiming.root" options.secondaryOutput="ntuple.root" -if(options.streamName=="AlCaP0"): options.files = "/store/data/Commissioning2015/AlCaP0/RAW/v1/000/246/342/00000/048ECF48-F906-E511-95AC-02163E011909.root" -elif(options.streamName=="AlCaPhiSym"): options.files = "/store/data/Commissioning2015/AlCaPhiSym/RAW/v1/000/244/768/00000/A8219906-44FD-E411-8DA9-02163E0121C5.root" +if(options.streamName=="AlCaP0"): print "stream ",options.streamName #options.files = "/store/data/Commissioning2015/AlCaP0/RAW/v1/000/246/342/00000/048ECF48-F906-E511-95AC-02163E011909.root" +elif(options.streamName=="AlCaPhiSym"): print "stream ",options.streamName #options.files = "/store/data/Commissioning2015/AlCaPhiSym/RAW/v1/000/244/768/00000/A8219906-44FD-E411-8DA9-02163E0121C5.root" else: print "stream ",options.streamName," not foreseen" exit + +#options.files = cms.untracked.vstring +#options.streamName = cms.untracked.vstring options.maxEvents = -1 # -1 means all events ### get and parse the command line arguments options.parseArguments() print options +# if the one file is a folder, grab all the files in it that are RECO +if len(options.files) == 1 and options.files[0][-1] == '/': + from EcalTiming.EcalTiming.storeTools_cff import fillFromStore + files = fillFromStore(options.files[0]) + options.files.pop() + options.files = [ f for f in files if "RECO" in f] + processname = options.step doReco = True doAnalysis = True if "RECO" not in processname: - doReco = False + doReco = False if "TIME" not in processname: - doAnalysis = False + doAnalysis = False process = cms.Process(processname) process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') process.load('FWCore.MessageService.MessageLogger_cfi') -process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(1000) +process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(5000) process.load('Configuration.EventContent.EventContent_cff') process.load('SimGeneral.MixingModule.mixNoPU_cfi') @@ -100,6 +141,7 @@ process.load('Configuration.StandardSequences.EndOfProcess_cff') process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff') process.load('EcalTiming.EcalTiming.ecalLocalRecoSequenceAlCaStream_cff') +process.load('EcalTiming.EcalTiming.ecalLocalRecoSequenceAlCaP0Stream_cff') if(options.streamName=="AlCaP0"): process.ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag("hltAlCaPi0EBRechitsToDigis","pi0EBDigis") @@ -122,7 +164,7 @@ ## GlobalTag Conditions Related from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, 'GR_P_V56', '') #run2_data', '') +process.GlobalTag = GlobalTag(process.GlobalTag, options.globaltag, '') ## Process Digi To Raw Step process.digiStep = cms.Sequence(process.ecalDigis + process.ecalPreshowerDigis) @@ -131,9 +173,6 @@ - - - ### Print Out Some Messages process.MessageLogger = cms.Service("MessageLogger", cout = cms.untracked.PSet( @@ -149,25 +188,31 @@ # SkipEvent = cms.untracked.vstring('ProductNotFound') ) +SkipEvent = cms.untracked.vstring('ProductNotFound','EcalProblem') + # dbs search --query "find file where dataset=/ExpressPhysics/BeamCommissioning09-Express-v2/FEVT and run=124020" | grep store | awk '{printf "\"%s\",\n", $1}' # Input source +print "source files:",options.files process.source = cms.Source("PoolSource", secondaryFileNames = cms.untracked.vstring(), fileNames = cms.untracked.vstring(options.files), + skipEvents = cms.untracked.uint32(options.skipEvents) ) if(len(options.jsonFile) > 0): - import FWCore.PythonUtilities.LumiList as LumiList - process.source.lumisToProcess = LumiList.LumiList(filename = options.jsonFile).getVLuminosityBlockRange() + import FWCore.PythonUtilities.LumiList as LumiList + process.source.lumisToProcess = LumiList.LumiList(filename = options.jsonFile).getVLuminosityBlockRange() +recofile = str(options.output) +recofile = recofile[:recofile.find(".root")] + "_RECO.root" # Output definition process.RECOoutput = cms.OutputModule("PoolOutputModule", splitLevel = cms.untracked.int32(0), eventAutoFlushCompressedSize = cms.untracked.int32(5242880), outputCommands = cms.untracked.vstring('drop *',"keep *_ecalRecHitE*Selector_*_*"), -fileName = cms.untracked.string(options.output), +fileName = cms.untracked.string(recofile), dataset = cms.untracked.PSet( filterName = cms.untracked.string(''), dataTier = cms.untracked.string('RECO') @@ -175,8 +220,9 @@ ) -## Histogram files -process.TFileService = cms.Service("TFileService", +if doAnalysis: + ## Histogram files + process.TFileService = cms.Service("TFileService", fileName = cms.string(options.output), closeFileFast = cms.untracked.bool(True) ) @@ -184,20 +230,76 @@ ### NumBer of events process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(options.maxEvents)) +#DUMMY RECHIT +process.dummyHits = cms.EDProducer("DummyRechitDigis", + doDigi = cms.untracked.bool(True), + # rechits + barrelHitProducer = cms.InputTag('hltAlCaPi0EBUncalibrator','pi0EcalRecHitsEB' ,"HLT"), + endcapHitProducer = cms.InputTag('hltAlCaPi0EEUncalibrator','pi0EcalRecHitsEE' ,"HLT"), + barrelRecHitCollection = cms.untracked.string("dummyBarrelRechitsPi0"), + endcapRecHitCollection = cms.untracked.string("dummyEndcapRechitsPi0"), + # digis + barrelDigis = cms.InputTag('hltAlCaPi0EBRechitsToDigis','pi0EBDigis',"HLT"), + endcapDigis = cms.InputTag('hltAlCaPi0EERechitsToDigisLowPU','pi0EEDigis',"HLT"), #changed hltAlCaPi0EERechitsToDigis in LowPU....changed in the file -.- + barrelDigiCollection = cms.untracked.string("dummyBarrelDigisPi0"), + endcapDigiCollection = cms.untracked.string("dummyEndcapDigisPi0")) + +##ADDED +# TRIGGER RESULTS FILTER +process.triggerSelectionLoneBunch = cms.EDFilter( "TriggerResultsFilter", + triggerConditions = cms.vstring('L1_AlwaysTrue'), + hltResults = cms.InputTag( "TriggerResults", "", "HLT" ), + l1tResults = cms.InputTag( "hltGtDigis" ), + l1tIgnoreMask = cms.bool( False ), + l1techIgnorePrescales = cms.bool( False ), + daqPartitions = cms.uint32( 1 ), + throw = cms.bool( True ) + ) process.filter=cms.Sequence() if(options.isSplash==1): process.filter+=process.spashesHltFilter process.reco_step = cms.Sequence(process.caloCosmicOrSplashRECOSequence) else: - #process.reco_step = cms.Sequence(process.reconstruction_step_multiFit) - process.reco_step = cms.Sequence(process.ecalLocalRecoSequenceAlCaStream) + if(options.streamName=="AlCaP0"): + #from RecoLocalCalo.Configuration.ecalLocalRecoSequence_cff import * + ##process.reco_step = cms.Sequence(ecalMultiFitUncalibRecHit * + ## ecalRecHit) + #process.ecalMultiFitUncalibRecHit = RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi.ecalMultiFitUncalibRecHit.clone() + + #process.ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag('dummyHits','dummyBarrelDigis')#,'piZeroAnalysis') + #process.ecalMultiFitUncalibRecHit.EEdigiCollection = cms.InputTag('dummyHits','dummyEndcapDigis')#,'piZeroAnalysis') + #ecalRecHit.killDeadChannels = False + #ecalRecHit.recoverEBFE = False + import RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi + process.ecalMultiFitUncalibRecHit = RecoLocalCalo.EcalRecProducers.ecalMultiFitUncalibRecHit_cfi.ecalMultiFitUncalibRecHit.clone() + process.ecalMultiFitUncalibRecHit.EBdigiCollection = cms.InputTag('dummyHits','dummyBarrelDigisPi0')#,'piZeroAnalysis') + process.ecalMultiFitUncalibRecHit.EEdigiCollection = cms.InputTag('dummyHits','dummyEndcapDigisPi0')#,'piZeroAnalysis') + + #UNCALIB to CALIB + from RecoLocalCalo.EcalRecProducers.ecalRecHit_cfi import * + process.ecalDetIdToBeRecovered = RecoLocalCalo.EcalRecProducers.ecalDetIdToBeRecovered_cfi.ecalDetIdToBeRecovered.clone() + process.ecalRecHit.killDeadChannels = cms.bool( False ) + process.ecalRecHit.recoverEBVFE = cms.bool( False ) + process.ecalRecHit.recoverEEVFE = cms.bool( False ) + process.ecalRecHit.recoverEBFE = cms.bool( False ) + process.ecalRecHit.recoverEEFE = cms.bool( False ) + process.ecalRecHit.recoverEEIsolatedChannels = cms.bool( False ) + process.ecalRecHit.recoverEBIsolatedChannels = cms.bool( False ) + + process.reco_step = cms.Sequence(process.dummyHits + * process.ecalMultiFitUncalibRecHit + * process.ecalRecHit) + else: + #process.reco_step = cms.Sequence(process.reconstruction_step_multiFit) + process.reco_step = cms.Sequence(process.ecalLocalRecoSequenceAlCaStream) + ### Process Full Path if(options.isSplash==0): process.digiStep = cms.Sequence() - - +if(options.loneBunch==1): + process.filter+=process.triggerSelectionLoneBunch evtPlots = True if options.isSplash else False @@ -209,22 +311,26 @@ #process.load("Geometry.CaloEventSetup.CaloGeometry_cff") #ESLooperProducer looper is imported here: -process.load('EcalTiming.EcalTiming.ecalTimingCalibProducer_cfi') process.load('EcalTiming.EcalTiming.RecHitsSelector_cfi') -process.timing.recHitEBCollection = cms.InputTag("ecalRecHitEBSelector") -process.timing.recHitEECollection = cms.InputTag("ecalRecHitEESelector") -process.timing.isSplash= cms.bool(True if options.isSplash else False) -process.timing.makeEventPlots=evtPlots -process.timing.globalOffset = cms.double(options.offset) -process.timing.outputDumpFile = process.TFileService.fileName -process.timing.energyThresholdOffsetEB = cms.double(options.minEnergyEB) -process.timing.energyThresholdOffsetEE = cms.double(options.minEnergyEE) -process.timing.storeEvents = cms.bool(True) +if doAnalysis: + process.load('EcalTiming.EcalTiming.ecalTimingCalibProducer_cfi') + process.timing.recHitEBCollection = cms.InputTag("ecalRecHitEBSelector") + process.timing.recHitEECollection = cms.InputTag("ecalRecHitEESelector") + process.timing.isSplash= cms.bool(True if options.isSplash else False) + process.timing.makeEventPlots=evtPlots + process.timing.globalOffset = cms.double(options.offset) + process.timing.outputDumpFile = process.TFileService.fileName + process.timing.energyThresholdOffsetEB = cms.double(options.minEnergyEB) + process.timing.energyThresholdOffsetEE = cms.double(options.minEnergyEE) + process.timing.chi2ThresholdOffsetEB = cms.double(options.minChi2EB) + process.timing.chi2ThresholdOffsetEE = cms.double(options.minChi2EE) + process.timing.storeEvents = cms.bool(True) + process.analysis = cms.Sequence( process.timing ) -process.analysis = cms.Sequence( process.timing ) -process.reco = cms.Sequence( (process.filter +if doReco: + process.reco = cms.Sequence( (process.filter + process.digiStep + process.reco_step) * (process.ecalRecHitEBSelector + process.ecalRecHitEESelector) @@ -233,11 +339,11 @@ process.seq = cms.Sequence() if doReco: - process.seq += process.reco + process.seq += process.reco if doAnalysis: - process.seq += process.analysis + process.seq += process.analysis else: - process.endp = cms.EndPath(process.RECOoutput) + process.endp = cms.EndPath(process.RECOoutput) process.p = cms.Path(process.seq) diff --git a/EcalTiming/test/launcher.sh b/EcalTiming/test/launcher.sh new file mode 100755 index 0000000..133548c --- /dev/null +++ b/EcalTiming/test/launcher.sh @@ -0,0 +1,109 @@ +#! /bin/bash + +############################################################################## +# USAGE # +# # +# ./launcher.sh /store/data/Run2015B/AlCaP0/RAW/v1/000/ 251 562 "AlCaP0" # +############################################################################## + +cd /afs/cern.ch/work/s/sgelli/public/CMSSW_7_4_12_patch4/src/ +export SCRAM_ARCH='slc6_amd64_gcc491' +eval `scramv1 runtime -sh` +source /afs/cern.ch/project/eos/installation/cms/etc/setup.sh +EOSCMD="/afs/cern.ch/project/eos/installation/cms/bin/eos.select" + +#FILES=$($EOSCMD ls /eos/cms/store/data/Run2015B/AlCaP0/RAW/v1/000/251/562/00000/) +#$EOSCMD ls /eos/cms/store/data/Run2015B/AlCaP0/RAW/v1/000/251/562/00000/ >test.txt + +#for files in $FILES +#do +# echo "Input " $files +#done + +#echo "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++" + +#while read line +#do +# echo "Input " $line +# +#done study only the run 251562 from PhiSym stream + + echo "++++++++++++++++++++++" + echo "+ Single Run Study +" + echo "++++++++++++++++++++++" + + runs=$($EOSCMD ls "/eos/cms/"$1"/"$2"/"$3"/00000/") + + echo "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++" + + #echo $TEST1 + + + #fileListSingle="'" + for run in $runs + do + fileList+=$"$1/$2/$3/00000/$run, " + #echo "Input test1 " $test1 + + done + fileList=${fileList%??} + #fileListSingle+="'" + echo $fileList + + #cmsRun /afs/cern.ch/work/s/sgelli/public/CMSSW_7_4_4_patch4/src/EcalTiming/EcalTiming/test/ecalTime_fromAlcaStream_cfg.py files="$fileList" + +fi + +echo $CMSSW_BASE +cd $CMSSW_BASE/src/EcalTiming/EcalTiming/test/ + + +if [ -z "$3" ]; then + cmsRun /afs/cern.ch/work/s/sgelli/public/CMSSW_7_4_12_patch4/src/EcalTiming/EcalTiming/test/ecalTime_fromAlcaStreamHighPU_cfg.py files="$fileList" output="/afs/cern.ch/work/s/sgelli/public/CMSSW_7_4_12_patch4/src/EcalTiming/EcalTiming/test/output/run"${firstFolder}${secondFolder}"_"${2}"_test_new.root" streamName=$2 jsonFile="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions15/13TeV/DCSOnly/json_DCSONLY.txt" +else + + cmsRun /afs/cern.ch/work/s/sgelli/public/CMSSW_7_4_12_patch4/src/EcalTiming/EcalTiming/test/ecalTime_fromAlcaStreamHighPU_cfg.py files="$fileList" output="/afs/cern.ch/work/s/sgelli/public/CMSSW_7_4_12_patch4/src/EcalTiming/EcalTiming/test/output/run"${2}${3}"_"${4}"_test_new.root" streamName=$4 jsonFile="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions15/13TeV/DCSOnly/json_DCSONLY.txt" + +fi + + +#cmsRun ecalTime_fromAlcaStream_cfg.py files='/store/data/Run2015A/AlCaP0/RAW/v1/000/246/908/00000/581F78FF-D909-E511-A356-02163E01231F.root, /store/data/Run2015A/AlCaP0/RAW/v1/000/246/908/00000/581F78FF-D909-E511-A356-02163E01231F.root' streamName=$4