diff --git a/EcalCalibAlgos/bin/BuildFile.xml b/EcalCalibAlgos/bin/BuildFile.xml index 943ac81..cb93905 100644 --- a/EcalCalibAlgos/bin/BuildFile.xml +++ b/EcalCalibAlgos/bin/BuildFile.xml @@ -19,8 +19,4 @@ name="TimeEvolution" file="TimeEvolution.cc"> - - diff --git a/EcalCalibAlgos/macros/drawSimpleChecks.py b/EcalCalibAlgos/macros/drawSimpleChecks.py index 5ad2615..7743895 100644 --- a/EcalCalibAlgos/macros/drawSimpleChecks.py +++ b/EcalCalibAlgos/macros/drawSimpleChecks.py @@ -34,6 +34,9 @@ histos["EB_EtMeanMap"]=ROOT.TH2F("EB_EtMeanMap","EB_EtMeanMap",360,0.5,360.5,171,-85.5,85.5) histos["EB_LCSumMap"]=ROOT.TH2F("EB_LCSumMap","EB_LCSumMap",360,0.5,360.5,171,-85.5,85.5) histos["EB_LCMap"]=ROOT.TH2F("EB_LCMap","EB_LCMap",360,0.5,360.5,171,-85.5,85.5) +histos["EB_TimeMean"]=ROOT.TH2F("EB_TimeMean","EB_TimeMean",360,0.5,360.5,171,-85.5,85.5) +histos["EB_TimeSum"]=ROOT.TH2F("EB_TimeSum","EB_TimeSum",360,0.5,360.5,171,-85.5,85.5) +histos["EB_TimeN"]=ROOT.TH2F("EB_TimeN","EB_TimeN",360,0.5,360.5,171,-85.5,85.5) histos["EEM_OccupancyMap"]=ROOT.TH2F("EEM_OccupancyMap","EEM_OccupancyMap",100,0.5,100.5,100,0.5,100.5) histos["EEM_EtMap"]=ROOT.TH2F("EEM_EtMap","EEM_EtMap",100,0.5,100.5,100,0.5,100.5) @@ -62,6 +65,8 @@ histos["EB_OccupancyMap"].Fill(myId.iphi(),myId.ieta(),hit.GetNhits()) histos["EB_EtMap"].Fill(myId.iphi(),myId.ieta(),hit.GetSumEt(0)) histos["EB_LCSumMap"].Fill(myId.iphi(),myId.ieta(),hit.GetLCSum()) + histos["EB_TimeSum"].Fill(myId.iphi(),myId.ieta(),hit.time_collection.GetTimeSum()) + histos["EB_TimeN"].Fill(myId.iphi(),myId.ieta(),hit.time_collection.GetTimeN()) for hit in phiSymRecHitsEE: myId=ROOT.EEDetId(hit.GetRawId()) if (myId.zside()<0): @@ -75,9 +80,12 @@ for iEta in range(1, 172): iBin = histos["EB_LCMap"].GetBin(iPhi, iEta) nHits = histos["EB_OccupancyMap"].GetBinContent(iBin) + nTimeHits = histos["EB_TimeN"].GetBinContent(iBin) if nHits > 0: histos["EB_EtMeanMap"].SetBinContent(iPhi, iEta, histos["EB_EtMap"].GetBinContent(iBin)/nHits) histos["EB_LCMap"].SetBinContent(iPhi, iEta, histos["EB_LCSumMap"].GetBinContent(iBin)/nHits) + if nTimeHits > 0: + histos["EB_TimeMean"].SetBinContent(iPhi, iEta, histos["EB_TimeSum"].GetBinContent(iBin)/nTimeHits) outFile=ROOT.TFile("phiSymStreamCheck.root","RECREATE") for histo in histos.keys(): diff --git a/EcalCalibAlgos/macros/timingAnalysis.py b/EcalCalibAlgos/macros/timingAnalysis.py new file mode 100644 index 0000000..c4755f2 --- /dev/null +++ b/EcalCalibAlgos/macros/timingAnalysis.py @@ -0,0 +1,98 @@ +# import ROOT in batch mode +import sys +oldargv = sys.argv[:] +sys.argv = [ '-b-' ] +import ROOT +ROOT.gROOT.SetBatch(True) +sys.argv = oldargv + +# load FWLite C++ librarie +ROOT.gSystem.Load("libFWCoreFWLite.so"); +ROOT.gSystem.Load("libDataFormatsFWLite.so"); +ROOT.gSystem.Load("libDataFormatsEcalDetId.so"); +ROOT.gSystem.Load("libPhiSymEcalCalibDataFormats.so"); +ROOT.AutoLibraryLoader.enable() + +# load FWlite python libraries +from DataFormats.FWLite import Handle, Events, Lumis + +# open file (you can use 'edmFileUtil -d /store/whatever.root' to get the physical file name) +#lumis = Lumis("file:phisym.root") +#lumis = Lumis("root://xrootd-cms.infn.it//store/user/spigazzi/AlCaPhiSym/crab_PHISYM-CMSSW_741-weights-GR_P_V56-Run2015B_v1/150714_150558/0000/phisym_weights_1lumis_13.root") +lumis = Lumis("file:phisym_weights_1lumis_withtime_numEvent500000.root") + +handlePhiSymInfo = Handle ("std::vector") +handlePhiSymRecHitsEB = Handle ("std::vector") +handlePhiSymRecHitsEE = Handle ("std::vector") +labelPhiSymInfo = ("PhiSymProducer") +labelPhiSymRecHitsEB = ("PhiSymProducer","EB") +labelPhiSymRecHitsEE = ("PhiSymProducer","EE") + +histos={} + +histos["EB_TimeMean"]=ROOT.TH2F("EB_TimeMean","EB_TimeMean",360,0.5,360.5,171,-85.5,85.5) +histos["EB_TimeSum"]=ROOT.TH2F("EB_TimeSum","EB_TimeSum",360,0.5,360.5,171,-85.5,85.5) +histos["EB_TimeN"]=ROOT.TH2F("EB_TimeN","EB_TimeN",360,0.5,360.5,171,-85.5,85.5) + +histos["EEM_TimeMean"]=ROOT.TH2F("EEM_TimeMean","EEM_TimeMean",100,0.5,100.5,100,0.5,100.5) +histos["EEM_TimeSum"]=ROOT.TH2F("EEM_TimeSum","EEM_TimeSum",100,0.5,100.5,100,0.5,100.5) +histos["EEM_TimeN"]=ROOT.TH2F("EEM_TimeN","EEM_TimeN",100,0.5,100.5,100,0.5,100.5) + +histos["EEP_TimeMean"]=ROOT.TH2F("EEP_TimeMean","EEP_TimeMean",100,0.5,100.5,100,0.5,100.5) +histos["EEP_TimeSum"]=ROOT.TH2F("EEP_TimeSum","EEP_TimeSum",100,0.5,100.5,100,0.5,100.5) +histos["EEP_TimeN"]=ROOT.TH2F("EEP_TimeN","EEP_TimeN",100,0.5,100.5,100,0.5,100.5) + +for i,lumi in enumerate(lumis): + print "====>" + lumi.getByLabel (labelPhiSymInfo,handlePhiSymInfo) + lumi.getByLabel (labelPhiSymRecHitsEB,handlePhiSymRecHitsEB) + lumi.getByLabel (labelPhiSymRecHitsEE,handlePhiSymRecHitsEE) + phiSymInfo = handlePhiSymInfo.product() + phiSymRecHitsEB = handlePhiSymRecHitsEB.product() + phiSymRecHitsEE = handlePhiSymRecHitsEE.product() + print "Run "+str(phiSymInfo.back().getStartLumi().run())+" Lumi "+str(phiSymInfo.back().getStartLumi().luminosityBlock()) + print "NEvents in this LS "+str(phiSymInfo.back().GetNEvents()) + print "TotHits EB "+str(phiSymInfo.back().GetTotHitsEB())+" Avg occ EB "+str(float(phiSymInfo.back().GetTotHitsEB())/phiSymInfo.back().GetNEvents()) + print "TotHits EE "+str(phiSymInfo.back().GetTotHitsEE())+" Avg occ EE "+str(float(phiSymInfo.back().GetTotHitsEE())/phiSymInfo.back().GetNEvents()) + + print "EB PhiSymRecHits "+str(phiSymRecHitsEB.size()) + print "EE PhiSymRecHits "+str(phiSymRecHitsEE.size()) + + for hit in phiSymRecHitsEB: + myId=ROOT.EBDetId(hit.GetRawId()) + histos["EB_TimeSum"].Fill(myId.iphi(),myId.ieta(),hit.GetTimeSum()) + histos["EB_TimeN"].Fill(myId.iphi(),myId.ieta(),hit.GetTimeN()) + for hit in phiSymRecHitsEE: + myId=ROOT.EEDetId(hit.GetRawId()) + if (myId.zside()<0): + histos["EEM_TimeSum"].Fill(myId.ix(),myId.iy(),hit.GetTimeSum()) + histos["EEM_TimeN"].Fill(myId.ix(),myId.iy(),hit.GetTimeN()) + else: + histos["EEP_TimeSum"].Fill(myId.ix(),myId.iy(),hit.GetTimeSum()) + histos["EEP_TimeN"].Fill(myId.ix(),myId.iy(),hit.GetTimeN()) + +for iPhi in range(1, 361): + for iEta in range(1, 172): + iBin = histos["EB_TimeN"].GetBin(iPhi, iEta) + nTimeHits = histos["EB_TimeN"].GetBinContent(iBin) + if nTimeHits > 0: + histos["EB_TimeMean"].SetBinContent(iPhi, iEta, histos["EB_TimeSum"].GetBinContent(iBin)/nTimeHits) + +for ix in range(1, 101): + for iy in range(1, 101): + iBin = histos["EEP_TimeN"].GetBin(ix, iy) + nTimeHits = histos["EEP_TimeN"].GetBinContent(iBin) + if nTimeHits > 0: + histos["EEP_TimeMean"].SetBinContent(ix, iy, histos["EEP_TimeSum"].GetBinContent(iBin)/nTimeHits) + + iBin = histos["EEM_TimeN"].GetBin(ix, iy) + nTimeHits = histos["EEM_TimeN"].GetBinContent(iBin) + if nTimeHits > 0: + histos["EEM_TimeMean"].SetBinContent(ix, iy, histos["EEM_TimeSum"].GetBinContent(iBin)/nTimeHits) + +outFile=ROOT.TFile("phiSymTimeAnalysis.root","RECREATE") +for histo in histos.keys(): + histos[histo].Write() +outFile.Write() +outFile.Close() + diff --git a/EcalCalibAlgos/plugins/PhiSymProducer.cc b/EcalCalibAlgos/plugins/PhiSymProducer.cc index a6e9c51..bbe4215 100644 --- a/EcalCalibAlgos/plugins/PhiSymProducer.cc +++ b/EcalCalibAlgos/plugins/PhiSymProducer.cc @@ -79,6 +79,9 @@ class PhiSymProducer : public edm::one::EDProducer recHitFlags_; + bool storeTimes_; + //---geometry EcalRingCalibrationTools calibRing_; static const short kNRingsEB = EcalRingCalibrationTools::N_RING_BARREL; @@ -117,6 +120,8 @@ PhiSymProducer::PhiSymProducer(const edm::ParameterSet& pSet): lumisToSum_(pSet.getParameter("lumisToSum")), statusThreshold_(pSet.getParameter("statusThreshold")), nLumis_(0), + recHitFlags_(pSet.getParameter >("recHitFlags")), + storeTimes_(pSet.getUntrackedParameter("storeTimes")), makeSpectraTreeEB_(pSet.getUntrackedParameter("makeSpectraTreeEB")), makeSpectraTreeEE_(pSet.getUntrackedParameter("makeSpectraTreeEE")) { @@ -324,6 +329,9 @@ void PhiSymProducer::produce(edm::Event& event, const edm::EventSetup& setup) recHitCollEB_->at(ebHit.denseIndex()).AddHit(etValues, laser.product()->getLaserCorrection(recHit.id(), evtTimeStamp)); + if(storeTimes_ && etValues[0] && recHit.checkFlags(recHitFlags_)) + recHitCollEB_->at(ebHit.denseIndex()).time_collection.AddTime(recHit.time()); + //---fill the plain tree if(makeSpectraTreeEB_) { @@ -381,6 +389,9 @@ void PhiSymProducer::produce(edm::Event& event, const edm::EventSetup& setup) //---update the rechHit sumEt recHitCollEE_->at(eeHit.denseIndex()).AddHit(etValues, laser.product()->getLaserCorrection(recHit.id(), evtTimeStamp)); + + if(storeTimes_ && recHit.checkFlags(recHitFlags_)) + recHitCollEE_->at(eeHit.denseIndex()).time_collection.AddTime(recHit.time()); //---fill the plain tree if(makeSpectraTreeEE_) diff --git a/EcalCalibAlgos/python/PhiSymProducer_cfi.py b/EcalCalibAlgos/python/PhiSymProducer_cfi.py index f410352..9248f4e 100644 --- a/EcalCalibAlgos/python/PhiSymProducer_cfi.py +++ b/EcalCalibAlgos/python/PhiSymProducer_cfi.py @@ -16,6 +16,8 @@ misCalibRangeEE = cms.vdouble(0.90, 1.10), lumisToSum = cms.int32(1), statusThreshold = cms.int32(0), + recHitFlags = cms.vint32([0]), # only recHits with these flags are accepted for calibration + storeTimes = cms.untracked.bool(True), makeSpectraTreeEB = cms.untracked.bool(False), makeSpectraTreeEE = cms.untracked.bool(False) ) diff --git a/EcalCalibAlgos/test/PhiSymProducer_cfg.py b/EcalCalibAlgos/test/PhiSymProducer_cfg.py index d85f221..9cff920 100644 --- a/EcalCalibAlgos/test/PhiSymProducer_cfg.py +++ b/EcalCalibAlgos/test/PhiSymProducer_cfg.py @@ -4,7 +4,7 @@ # parse commad line options options = VarParsing('analysis') -options.maxEvents = -1 +options.maxEvents = 50000 options.outputFile = 'phisym_weights_1lumis.root' options.parseArguments() @@ -29,7 +29,7 @@ # import of standard configurations process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(50000) + input = cms.untracked.int32(options.maxEvents) ) # skip bad events diff --git a/EcalCalibDataFormats/interface/PhiSymRecHit.h b/EcalCalibDataFormats/interface/PhiSymRecHit.h index 2fe64a8..a07af70 100644 --- a/EcalCalibDataFormats/interface/PhiSymRecHit.h +++ b/EcalCalibDataFormats/interface/PhiSymRecHit.h @@ -16,6 +16,7 @@ #include #include "DataFormats/DetId/interface/DetId.h" +#include "PhiSym/EcalCalibDataFormats/interface/PhiSymTimeCollection.h" //---define the number of allowed mis-calibrated values for etSum_ (+ the central value) #define N_MISCALIB_VALUES 11 @@ -23,6 +24,8 @@ class PhiSymRecHit { public: + typedef std::vector timevec; + //---ctors--- PhiSymRecHit(); PhiSymRecHit(uint32_t id, float* etValues=NULL); @@ -40,12 +43,18 @@ class PhiSymRecHit //---utils--- void AddHit(float* etValues, float laserCorr=0); + int16_t CompressTime(float t) const; + float UncompressTime(int16_t t) const; void Reset(); //---operators--- PhiSymRecHit& operator+=(const PhiSymRecHit& rhs); friend std::ostream& operator<<(std::ostream& out, const PhiSymRecHit& obj); + //---Public Members--- + + PhiSymTimeCollection time_collection; + private: uint32_t id_; diff --git a/EcalCalibDataFormats/interface/PhiSymTimeCollection.h b/EcalCalibDataFormats/interface/PhiSymTimeCollection.h new file mode 100644 index 0000000..b17fb8e --- /dev/null +++ b/EcalCalibDataFormats/interface/PhiSymTimeCollection.h @@ -0,0 +1,61 @@ +#ifndef DATAFORMAT_PHISYMTIMECOLLECTION_H +#define DATAFORMAT_PHISYMTIMECOLLECTION_H + +/** \class PhiSymTimeCollection + * + * Dataformat dedicated to Phi Symmetry ecal time calibration + * + */ + +#include +#include +#include + +//---define the number of allowed mis-calibrated values for etSum_ (+ the central value) +#define MAX_TIME 32.767 + +typedef std::vector TimeCollection; + +class PhiSymTimeCollection +{ +public: + //---ctors--- + PhiSymTimeCollection(); + + //---dtor--- + ~PhiSymTimeCollection(); + + //---getters--- + inline const TimeCollection GetTimes() const {return times_;}; + float GetTimeSum() const; + float GetTimeSum2() const; + size_t GetTimeN() const; + float GetTimeMean() const; + float GetTimeStdDev() const; + //---getters for time in range--- + const TimeCollection GetTimes(float low, float hi) const; + float GetTimeMean(float low, float hi) const; + float GetTimeSum(float low, float hi) const; + float GetTimeSum2(float low, float hi) const; + size_t GetTimeN(float low, float hi) const; + float GetTimeStdDev(float low, float hi) const; + + //---utils--- + void AddTime(float t); + int16_t CompressTime(float t) const; + float UncompressTime(int16_t t) const; + void Reset(); + + //---operators--- + PhiSymTimeCollection& operator+=(const PhiSymTimeCollection& rhs); + friend std::ostream& operator<<(std::ostream& out, const PhiSymTimeCollection& obj); + +private: + + TimeCollection times_; + +}; + + +#endif + diff --git a/EcalCalibDataFormats/src/PhiSymRecHit.cc b/EcalCalibDataFormats/src/PhiSymRecHit.cc index 62967f5..a57806e 100644 --- a/EcalCalibDataFormats/src/PhiSymRecHit.cc +++ b/EcalCalibDataFormats/src/PhiSymRecHit.cc @@ -1,4 +1,5 @@ #include "PhiSym/EcalCalibDataFormats/interface/PhiSymRecHit.h" +#include "TMath.h" //**********constructors****************************************************************** PhiSymRecHit::PhiSymRecHit(): @@ -46,6 +47,7 @@ void PhiSymRecHit::Reset() lc2Sum_ = 0; for(short i=0; itime_collection.AddTime(t); return *this; } @@ -75,6 +79,7 @@ std::ostream& operator<<(std::ostream& out, const PhiSymRecHit& obj) out << std::setw(20) << "sumEt2: " << obj.GetSumEt2() << std::endl; out << std::setw(20) << "laser-corr sum: " << obj.GetLCSum() << std::endl; out << std::setw(20) << "laser-corr^2 sum: " << obj.GetLC2Sum() << std::endl; + out << obj.time_collection; out << std::endl; return out; diff --git a/EcalCalibDataFormats/src/PhiSymTimeCollection.cc b/EcalCalibDataFormats/src/PhiSymTimeCollection.cc new file mode 100644 index 0000000..ac9711c --- /dev/null +++ b/EcalCalibDataFormats/src/PhiSymTimeCollection.cc @@ -0,0 +1,129 @@ +#include "PhiSym/EcalCalibDataFormats/interface/PhiSymTimeCollection.h" +#include "TMath.h" + +//**********constructors****************************************************************** +PhiSymTimeCollection::PhiSymTimeCollection() +{ +} + +//**********destructor******************************************************************** +PhiSymTimeCollection::~PhiSymTimeCollection() +{} + +//**********getters*********************************************************************** +float PhiSymTimeCollection::GetTimeSum() const +{ + return GetTimeSum(0.0f, -1.0f); +} + +float PhiSymTimeCollection::GetTimeSum2() const +{ + return GetTimeSum2(0.0f, -1.0f); +} + +size_t PhiSymTimeCollection::GetTimeN() const +{ + return GetTimeN(0.0f, -1.0f); +} + + +const TimeCollection PhiSymTimeCollection::GetTimes(float low, float hi) const +{ + TimeCollection ret; + bool range = true; + if(hi < low) + range = false; + auto low_comp = CompressTime(low); + auto hi_comp = CompressTime(hi); + for(auto t : times_) + { + if(!range || (t >= low_comp && t <= hi_comp)) + ret.push_back(t); + } + return ret; +} + +float PhiSymTimeCollection::GetTimeSum(float low, float hi) const +{ + float time = 0; + for( auto t : GetTimes(low,hi)) + { + time += UncompressTime(t); + } + return time; +} +float PhiSymTimeCollection::GetTimeSum2(float low, float hi) const +{ + float time2 = 0; + for( auto t : GetTimes(low,hi)) + { + time2 += TMath::Power(UncompressTime(t),2); + } + return time2; +} +size_t PhiSymTimeCollection::GetTimeN(float low, float hi) const +{ + return GetTimes(low,hi).size(); +} +//**********utils************************************************************************* + +void PhiSymTimeCollection::AddTime(float t) +{ + //Range of 16bit int + if(t >= - MAX_TIME && t <= MAX_TIME) + times_.push_back(CompressTime(t)); +} + +int16_t PhiSymTimeCollection::CompressTime(float t) const +{ + return (int16_t)TMath::Nint(t*1000); +} +float PhiSymTimeCollection::UncompressTime(int16_t t) const +{ + return float(t)/1000.0f; +} + +void PhiSymTimeCollection::Reset() +{ + times_.clear(); +} + +float PhiSymTimeCollection::GetTimeMean() const +{ + return GetTimeMean(0.0f, -1.0f); +} + +float PhiSymTimeCollection::GetTimeMean(float low, float hi) const +{ + return GetTimeSum(low,hi)/GetTimeN(low,hi); +} + +float PhiSymTimeCollection::GetTimeStdDev() const +{ + return GetTimeStdDev(0.0f, -1.0f); +} + +float PhiSymTimeCollection::GetTimeStdDev(float low, float hi) const +{ + float mean = GetTimeMean(low,hi); + return TMath::Sqrt( GetTimeSum2(low,hi)/GetTimeN(low,hi) - mean*mean); +} + +//**********operators********************************************************************* + +PhiSymTimeCollection& PhiSymTimeCollection::operator+=(const PhiSymTimeCollection& rhs) +{ + for(auto t : rhs.GetTimes()) + this->times_.push_back(t); + return *this; +} + +std::ostream& operator<<(std::ostream& out, const PhiSymTimeCollection& obj) +{ + //---dump all the informations + out << std::setw(20) << "number of times:" << obj.GetTimeN() << std::endl; + out << std::endl; + + return out; +} + diff --git a/EcalCalibDataFormats/src/classes.h b/EcalCalibDataFormats/src/classes.h index 8059bff..2c25979 100644 --- a/EcalCalibDataFormats/src/classes.h +++ b/EcalCalibDataFormats/src/classes.h @@ -1,6 +1,7 @@ #include "DataFormats/Common/interface/Wrapper.h" #include "PhiSym/EcalCalibDataFormats/interface/PhiSymRecHit.h" +#include "PhiSym/EcalCalibDataFormats/interface/PhiSymTimeCollection.h" #include "PhiSym/EcalCalibDataFormats/interface/PhiSymInfo.h" #include "PhiSym/EcalCalibDataFormats/interface/CalibrationFile.h" #include "PhiSym/EcalCalibDataFormats/interface/PhiSymFile.h" @@ -25,6 +26,9 @@ namespace CrystalsEBTree dummy42; CrystalsEETree dummy43; CalibrationFile dummy44; + + PhiSymTimeCollection dummy51; + edm::Wrapper dummy52; }; } diff --git a/EcalCalibDataFormats/src/classes_def.xml b/EcalCalibDataFormats/src/classes_def.xml index ee9de8d..a2cd9cb 100644 --- a/EcalCalibDataFormats/src/classes_def.xml +++ b/EcalCalibDataFormats/src/classes_def.xml @@ -13,4 +13,7 @@ + + +