EcalTiming: plugins/EcalTimingCalibProducer.h Source File

00001 #ifndef ecaltimingcalibproducer_h
00002 #define ecaltimingcalibproducer_h
00003 /** \class EcalTimingCalibProducer EcalTimingCalibProducer.h EcalTiming/EcalTiming/plugins/EcalTimingCalibProducer.h
00004     \brief Plugin that derives the calibration constants
00005 
00006  This plugin runs over the events, selects the recHits according to the criteria defined in addRecHit
00007 
00008  \todo Exit condition based on convergence
00009 
00010 */
00011 
00012 
00013 /**
00014    Module description:
00015  - digi to calibrated recHit reconstruction
00016  - selection of events with reasonable activity: minimum number of recHits
00017  - select recHits based on:
00018    - recoFlag -> only good recHits (exclude OOT pileup contribution with MultiFit)
00019    - minimum energy, ring based threshold
00020  - save all the time events (recHits) passing the selection
00021  - discard those not within 2 sigma of the stdDev distribution for the single channel (then excluding OOT spurious events)
00022  - exclude events with large time error (> 3ns) or null time error (==0 ns)
00023  - verify that the distribution is symmetric
00024  - check stability vs energy of the single channel within uncertainty
00025  - save full dump of time events in TTree for further checks if:
00026    - distribution not symmetric
00027    - time calibration not stable vs energy
00028  - calculate global time for EB, EE
00029 
00030 
00031  */
00032 //#define DEBUG
00033 #define RAWIDCRY 838904321
00034 
00035 /// The entire set of events are saved for one channel in EB and one channel in EE in order to debug
00036 #define EBCRYex 838861346
00037 #define EECRYex 872422180
00038 
00039 /// For one ring in EB, one in EE+ and one in EE- the entire set of events are saved to debug
00040 #define EBRING 1
00041 #define EEmRING 20
00042 #define EEpRING 20
00043 
00044 #define SPEEDOFLIGHT 30.0 // (cm/ns)
00045 
00046 // system include files
00047 #include <memory>
00048 #include <iostream>
00049 #include <fstream>
00050 
00051 // Make Histograms the way!!
00052 #include "FWCore/ServiceRegistry/interface/Service.h"
00053 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00054 
00055 // user include files
00056 #include "FWCore/Framework/interface/Frameworkfwd.h"
00057 #include "FWCore/Framework/interface/EDConsumerBase.h"
00058 
00059 #include "FWCore/Framework/interface/EDProducer.h"
00060 #include "FWCore/Framework/interface/Event.h"
00061 //#include "FWCore/Framework/interface/Handle.h"
00062 #include "FWCore/Framework/interface/MakerMacros.h"
00063 #include "FWCore/Utilities/interface/InputTag.h"
00064 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00065 
00066 #include "DataFormats/Common/interface/Handle.h"
00067 #include "FWCore/Framework/interface/LooperFactory.h"
00068 #include "FWCore/Framework/interface/ESProducerLooper.h"
00069 #include "FWCore/Framework/interface/ESHandle.h"
00070 #include "FWCore/Framework/interface/ESProducts.h"
00071 #include "FWCore/Framework/interface/Event.h"
00072 // input collections
00073 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00074 //RingTools
00075 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
00076 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00077 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00078 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00079 
00080 // record to be produced:
00081 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
00082 #include "CondFormats/DataRecord/interface/EcalTimeCalibErrorsRcd.h"
00083 #include "CondFormats/DataRecord/interface/EcalTimeOffsetConstantRcd.h"
00084 #include "CondTools/Ecal/interface/EcalTimeCalibConstantsXMLTranslator.h"
00085 #include "CondTools/Ecal/interface/EcalTimeCalibErrorsXMLTranslator.h"
00086 #include "CondTools/Ecal/interface/EcalTimeOffsetXMLTranslator.h"
00087 #include "CondTools/Ecal/interface/EcalCondHeader.h"
00088 
00089 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
00090 #include "CondFormats/EcalObjects/interface/EcalTimeCalibErrors.h"
00091 #include "CondFormats/EcalObjects/interface/EcalTimeOffsetConstant.h"
00092 
00093 #include "EcalTiming/EcalTiming/interface/EcalTimingEvent.h"
00094 #include "EcalTiming/EcalTiming/interface/EcalCrystalTimingCalibration.h"
00095 
00096 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00097 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00098 
00099 #include "TProfile.h"
00100 #include "TGraphErrors.h"
00101 #include "TGraph.h"
00102 #include "TH1F.h"
00103 #include "TH2F.h"
00104 #include "TFile.h"
00105 #include "TProfile2D.h"
00106 
00107 #include <fstream>
00108 #include <string>
00109 #include <vector>
00110 #include <iostream>
00111 #include <map>
00112 #include <vector>
00113 #include <algorithm>
00114 #include <functional>
00115 #include <set>
00116 #include <assert.h>
00117 #include <time.h>
00118 
00119 #include <TMath.h>
00120 #include <Math/VectorUtil.h>
00121 //#include <boost/tokenizer.hpp>
00122 
00123 #include "EcalTiming/EcalTiming/interface/EcalTimeCalibrationMapFwd.h"
00124 
00125 class EcalTimingCalibProducer : public edm::ESProducerLooper
00126 {
00127 
00128 private:
00129     EcalTimeCalibrationMap _timeCalibMap; ///< calibration map: contains the time shift for each crystal
00130     EventTimeMap _eventTimeMap;           ///< container of recHits passing selection in the event (reset at each event)
00131 
00132     // For finding averages for specific eta ring
00133     EcalCrystalTimingCalibration timeEEP; ///< global time calibration of EE+
00134     EcalCrystalTimingCalibration timeEEM; ///< global time calibration of EE-
00135     EcalCrystalTimingCalibration timeEB; ///< global time calibration of EB
00136     EcalCrystalTimingCalibration timeEBRing; ///< global time calibration of one EB ring
00137     EcalCrystalTimingCalibration timeEEmRing; ///< global time calibration of one EE- ring
00138     EcalCrystalTimingCalibration timeEEpRing; ///< global time calibration of one EE+ ring
00139     EcalCrystalTimingCalibration timeEBCRYex; ///< global time calibration of one EB channel
00140     EcalCrystalTimingCalibration timeEECRYex; ///< global time calibration of one EE channel
00141 
00142     float nearEndcapTime;
00143     float farEndcapTime;
00144 
00145 public:
00146     EcalTimingCalibProducer(const edm::ParameterSet&); // default constructor
00147     ~EcalTimingCalibProducer();                        // default destructor
00148 
00149 
00150     virtual void beginOfJob(const edm::EventSetup&);
00151     virtual void startingNewLoop(unsigned int ) ;
00152     virtual Status duringLoop(const edm::Event&, const edm::EventSetup&) ;
00153     virtual Status endOfLoop(const edm::EventSetup&, unsigned int);
00154     virtual void endOfJob();
00155 private:
00156     // ----------member data ---------------------------
00157     /** @name Input Parameters
00158      * Parameters defined in the config file _cfi,py
00159      */
00160     ///@{
00161 
00162     unsigned int _maxLoop; ///< maximum number of loops for intercalibration
00163     bool _isSplash; ///< flag to activate for splash analysis
00164     bool _makeEventPlots; ///< flag for making plots for each event (meant for splashes)
00165     edm::InputTag _ecalRecHitsEBTAG; ///< input collection
00166     edm::InputTag _ecalRecHitsEETAG;///< input collection
00167     std::vector<int> _recHitFlags; ///< vector containing list of valid rec hit flags for calibration
00168     unsigned int _recHitMin; ///< require at least this many rec hits to count the event
00169     double       _minRecHitEnergy; ///< minimum energy for the recHit to be considered for timing studies
00170     double _minRecHitEnergyStep; ///< at each iteration the code increases the minimum energy required for the calibration
00171     unsigned int _minEntries; ///< require a minimum number of entries in a ring to do averages
00172     float        _globalOffset;    ///< time to subtract from every event
00173     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
00174     std::string _outputDumpFileName; ///< name of the output file for the calibration constants' dump
00175     float _noiseRMSThreshold;
00176     float _noiseTimeThreshold;
00177     float _maxSkewnessForDump;
00178 /// @}
00179 
00180     void dumpCalibration(std::string filename);
00181     void dumpCorrections(std::string filename);
00182 
00183 // plotting
00184 ///fill histograms with the measured shifts (that will become -corrections for the next step)
00185     void FillCalibrationCorrectionHists(EcalTimeCalibrationMap::const_iterator cal_itr);
00186     void initHists(TFileDirectory dir);
00187     void initEventHists(TFileDirectory dir);
00188     void initTree(TFileDirectory dir);
00189 
00190     // Create calibration container objects -> to be used in beginOfJob
00191     void createConstants(const edm::EventSetup& iSetup)
00192     {
00193 
00194         // time calib constants
00195         edm::ESHandle<EcalTimeCalibConstants> ecalTimeCalibConstantsHandle;
00196         iSetup.get<EcalTimeCalibConstantsRcd>().get( ecalTimeCalibConstantsHandle);
00197         _timeCalibConstants = *ecalTimeCalibConstantsHandle;
00198 #ifdef DEBUG
00199         // why they are 0 when setWhatProduced is uncommented? -> I've not found a way to read the default DB at iter 0
00200         for(auto t_itr = _timeCalibConstants.begin(); t_itr != _timeCalibConstants.end() && t_itr - _timeCalibConstants.begin() < 10; ++t_itr) {
00201             std::cout << "t_itr = " << *t_itr << std::endl;
00202         }
00203 #endif
00204         produceCalibConstants(iSetup.get<EcalTimeCalibConstantsRcd>());
00205 
00206         // time calib errors
00207         edm::ESHandle<EcalTimeCalibErrors> ecalTimeCalibErrorsHandle;
00208         //iSetup.get<EcalTimeCalibErrorsRcd>().get( ecalTimeCalibErrorsHandle);
00209         //_timeCalibErrors = *ecalTimeCalibErrorsHandle;
00210         //produceCalibErrors(iSetup.get<EcalTimeCalibErrorsRcd>());
00211 
00212         // time offset constant
00213         edm::ESHandle<EcalTimeOffsetConstant> ecalTimeOffsetConstantHandle;
00214         iSetup.get<EcalTimeOffsetConstantRcd>().get( ecalTimeOffsetConstantHandle);
00215         _timeOffsetConstant = *ecalTimeOffsetConstantHandle;
00216         produceOffsetConstant(iSetup.get<EcalTimeOffsetConstantRcd>());
00217 
00218     }
00219 
00220     EcalTimeCalibConstants _timeCalibConstants; ///< container of calibrations updated iter by iter
00221     boost::shared_ptr<EcalTimeCalibConstants> _calibConstants; // corrections used during this iter
00222     inline boost::shared_ptr<EcalTimeCalibConstants> produceCalibConstants(const EcalTimeCalibConstantsRcd& iRecord)
00223     {
00224         // new shared_ptr initialized with previous iter: _timeCalibConstants
00225         _calibConstants = boost::shared_ptr<EcalTimeCalibConstants>( new EcalTimeCalibConstants(_timeCalibConstants) );
00226         return  _calibConstants;
00227     }
00228 
00229     EcalTimeCalibErrors    _timeCalibErrors;
00230     boost::shared_ptr<EcalTimeCalibErrors> _calibErrors;
00231     inline boost::shared_ptr<EcalTimeCalibErrors>& produceCalibErrors(const EcalTimeCalibErrorsRcd& iRecord)
00232     {
00233         _calibErrors = boost::shared_ptr<EcalTimeCalibErrors>( new EcalTimeCalibErrors(_timeCalibErrors) );
00234         return _calibErrors;
00235     }
00236 
00237     EcalTimeOffsetConstant _timeOffsetConstant;
00238     boost::shared_ptr<EcalTimeOffsetConstant> _offsetConstant;
00239     inline boost::shared_ptr<EcalTimeOffsetConstant> produceOffsetConstant(const EcalTimeOffsetConstantRcd& iRecord)
00240     {
00241         _offsetConstant = boost::shared_ptr<EcalTimeOffsetConstant>( new EcalTimeOffsetConstant(_timeOffsetConstant) );
00242         return _offsetConstant;
00243     }
00244 
00245 
00246 
00247     /**
00248       \brief  If recHit passes the selection it is added to the list of recHits to be used for calibration
00249 
00250       The recHit is used (accepted) if:
00251         - recHit flag in the list of _recHitFlags defined in the config file
00252         - recHit energy in EB > _minRecHitEnergy defined in the config file
00253         - recHit energy in EE is x2 the threshold of EB
00254         - at each iteration the recHit threshold is raised by _minRecHitEnergyStep
00255         If the recHit is used, the time information is added to _eventTimeMap
00256     */
00257     bool addRecHit(const EcalRecHit& recHit, EventTimeMap& eventTimeMap_);
00258 
00259 
00260     /// Adds the recHit to the per Event histograms
00261     void plotRecHit(const EcalTimingEvent& recHit);
00262     ///
00263     /// Returns an EcalTimingEvent with a new time, which has been adjusted
00264     /// so that the upstream endcap is 0.
00265     ///
00266     EcalTimingEvent correctGlobalOffset(const EcalTimingEvent& ev, int splashDir, float bunchCorr);
00267 
00268     std::map<DetId, float>  _CrysEnergyMap;
00269 
00270     edm::Service<TFileService> fileService_;
00271     TFileDirectory histDir_;
00272     // Tree
00273     TTree *timeEBCRYexTree, *timeEECRYexTree, *timeEEpRingTree, *timeEEmRingTree, *timeEBRingTree;
00274     TTree *_unstableEnergyTree, *_highSkewnessTree;
00275 
00276     // Mean Histograms
00277     TProfile2D* EneMapEEP_; /// Using TProfile2D so we don't paint empty bins.
00278     TProfile2D* EneMapEEM_;
00279     TProfile2D* TimeMapEEP_;
00280     TProfile2D* TimeMapEEM_;
00281 
00282     TProfile2D* EneMapEB_;
00283     TProfile2D* TimeMapEB_;
00284 
00285     // Error Histograms
00286     TProfile2D* TimeErrorMapEEP_;
00287     TProfile2D* TimeErrorMapEEM_;
00288 
00289     TProfile2D* TimeErrorMapEB_;
00290 
00291     // Event Based Plots
00292     TProfile2D * Event_EneMapEEP_;
00293     TProfile2D * Event_EneMapEEM_;
00294     TProfile2D * Event_EneMapEB_;
00295 
00296     TProfile2D* Event_TimeMapEEP_;
00297     TProfile2D* Event_TimeMapEEM_;
00298     TProfile2D* Event_TimeMapEB_;
00299 
00300     TProfile2D* Event_TimeMapEEP_OOT;
00301     TProfile2D* Event_TimeMapEEM_OOT;
00302     TProfile2D* Event_TimeMapEB_OOT;
00303 
00304     TH1F* RechitEneEB_;
00305     TH1F* RechitTimeEB_;
00306     TH1F* RechitEneEEM_;
00307     TH1F* RechitTimeEEM_;
00308     TH1F* RechitEneEEP_;
00309     TH1F* RechitTimeEEP_;
00310 
00311     TH2F* RechitEnergyTimeEB;
00312     TH2F* RechitEnergyTimeEEM;
00313     TH2F* RechitEnergyTimeEEP;
00314 
00315     EcalRingCalibrationTools _ringTools;
00316     const CaloSubdetectorGeometry * endcapGeometry_;
00317     const CaloSubdetectorGeometry * barrelGeometry_;
00318 
00319     std::vector<int> _noisyXtals;
00320     std::vector<TH1F*> _noisyXtalsHists;
00321     unsigned int _iter;
00322 };
00323 
00324 #endif