added timing PhiSymProducer, implemented as a new class

EcalCalibAlgos/bin/BuildFile.xml

@@ -19,8 +19,4 @@
       name="TimeEvolution"
       file="TimeEvolution.cc">
   </bin>
-  <bin
-      name="fuckingEndcaps"
-      file="fuckingEndcaps.cc">
-  </bin>
 </environment>

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():

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<PhiSymInfo>")
+handlePhiSymRecHitsEB  = Handle ("std::vector<PhiSymRecHit>")
+handlePhiSymRecHitsEE  = Handle ("std::vector<PhiSymRecHit>")
+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()
+

EcalCalibAlgos/plugins/PhiSymProducer.cc

@@ -79,6 +79,9 @@ class PhiSymProducer : public edm::one::EDProducer<edm::EndLuminosityBlockProduc
     int            lumisToSum_;
     int            statusThreshold_;
     int            nLumis_;
+    std::vector<int> 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<int>("lumisToSum")),
     statusThreshold_(pSet.getParameter<int>("statusThreshold")),
     nLumis_(0),
+    recHitFlags_(pSet.getParameter<std::vector<int> >("recHitFlags")),
+    storeTimes_(pSet.getUntrackedParameter<bool>("storeTimes")),
     makeSpectraTreeEB_(pSet.getUntrackedParameter<bool>("makeSpectraTreeEB")),
     makeSpectraTreeEE_(pSet.getUntrackedParameter<bool>("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_)

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)
 )

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

EcalCalibDataFormats/interface/PhiSymRecHit.h

@@ -16,13 +16,16 @@
 #include <vector>
 
 #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
 
 class PhiSymRecHit
 {
 public:
+    typedef std::vector<int16_t> 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_;

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 <iostream>
+#include <iomanip>
+#include <vector>
+
+//---define the number of allowed mis-calibrated values for etSum_ (+ the central value)
+#define MAX_TIME 32.767
+
+typedef std::vector<int16_t> 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
+

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; i<N_MISCALIB_VALUES; ++i)
         etSum_[i] = 0;
+    time_collection.Reset();
 }
 
 //**********operators*********************************************************************
@@ -61,6 +63,8 @@ PhiSymRecHit& PhiSymRecHit::operator+=(const PhiSymRecHit& rhs)
     for(short i=0; i<N_MISCALIB_VALUES; ++i)
         etSum_[i] += rhs.GetSumEt(i);
 
+    for(auto t : rhs.time_collection.GetTimes())
+       this->time_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;