00001 #include "Validation/EcalClusters/interface/EgammaSuperClusters.h"
00002
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "FWCore/Utilities/interface/Exception.h"
00005 #include "FWCore/ServiceRegistry/interface/Service.h"
00006
00007 #include "DataFormats/Common/interface/Handle.h"
00008 #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
00009 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00010 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00011 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00012 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
00013
00014 #include "DataFormats/GeometryVector/interface/Pi.h"
00015 #include "DQMServices/Core/interface/DQMStore.h"
00016
00017 EgammaSuperClusters::EgammaSuperClusters( const edm::ParameterSet& ps )
00018 {
00019 outputFile_ = ps.getUntrackedParameter<std::string>("outputFile", "");
00020 CMSSW_Version_ = ps.getUntrackedParameter<std::string>("CMSSW_Version", "");
00021
00022 verboseDBE_ = ps.getUntrackedParameter<bool>("verboseDBE", false);
00023
00024 hist_min_Size_ = ps.getParameter<double>("hist_min_Size");
00025 hist_max_Size_ = ps.getParameter<double>("hist_max_Size");
00026 hist_bins_Size_ = ps.getParameter<int> ("hist_bins_Size");
00027
00028 hist_min_NumBC_ = ps.getParameter<double>("hist_min_NumBC");
00029 hist_max_NumBC_ = ps.getParameter<double>("hist_max_NumBC");
00030 hist_bins_NumBC_ = ps.getParameter<int> ("hist_bins_NumBC");
00031
00032 hist_min_ET_ = ps.getParameter<double>("hist_min_ET");
00033 hist_max_ET_ = ps.getParameter<double>("hist_max_ET");
00034 hist_bins_ET_ = ps.getParameter<int> ("hist_bins_ET");
00035
00036 hist_min_Eta_ = ps.getParameter<double>("hist_min_Eta");
00037 hist_max_Eta_ = ps.getParameter<double>("hist_max_Eta");
00038 hist_bins_Eta_ = ps.getParameter<int> ("hist_bins_Eta");
00039
00040 hist_min_Phi_ = ps.getParameter<double>("hist_min_Phi");
00041 hist_max_Phi_ = ps.getParameter<double>("hist_max_Phi");
00042 hist_bins_Phi_ = ps.getParameter<int> ("hist_bins_Phi");
00043
00044 hist_min_S1toS9_ = ps.getParameter<double>("hist_min_S1toS9");
00045 hist_max_S1toS9_ = ps.getParameter<double>("hist_max_S1toS9");
00046 hist_bins_S1toS9_ = ps.getParameter<int> ("hist_bins_S1toS9");
00047
00048 hist_min_S25toE_ = ps.getParameter<double>("hist_min_S25toE");
00049 hist_max_S25toE_ = ps.getParameter<double>("hist_max_S25toE");
00050 hist_bins_S25toE_ = ps.getParameter<int> ("hist_bins_S25toE");
00051
00052 hist_min_EToverTruth_ = ps.getParameter<double>("hist_min_EToverTruth");
00053 hist_max_EToverTruth_ = ps.getParameter<double>("hist_max_EToverTruth");
00054 hist_bins_EToverTruth_ = ps.getParameter<int> ("hist_bins_EToverTruth");
00055
00056 hist_min_deltaEta_ = ps.getParameter<double>("hist_min_deltaEta");
00057 hist_max_deltaEta_ = ps.getParameter<double>("hist_max_deltaEta");
00058 hist_bins_deltaEta_ = ps.getParameter<int> ("hist_bins_deltaEta");
00059
00060 MCTruthCollection_ = ps.getParameter<edm::InputTag>("MCTruthCollection");
00061 barrelSuperClusterCollection_ = ps.getParameter<edm::InputTag>("barrelSuperClusterCollection");
00062 endcapSuperClusterCollection_ = ps.getParameter<edm::InputTag>("endcapSuperClusterCollection");
00063 barrelRecHitCollection_ = ps.getParameter<edm::InputTag>("barrelRecHitCollection");
00064 endcapRecHitCollection_ = ps.getParameter<edm::InputTag>("endcapRecHitCollection");
00065 }
00066
00067 EgammaSuperClusters::~EgammaSuperClusters() {}
00068
00069 void EgammaSuperClusters::beginJob(edm::EventSetup const&)
00070 {
00071 dbe_ = edm::Service<DQMStore>().operator->();
00072
00073 if ( verboseDBE_ )
00074 {
00075 dbe_->setVerbose(1);
00076 dbe_->showDirStructure();
00077 }
00078 else
00079 dbe_->setVerbose(0);
00080
00081 dbe_->setCurrentFolder("EcalClustersV/CMSSW_"+CMSSW_Version_+"/EcalClusters/SuperClusters/");
00082
00083 hist_EB_SC_Size_
00084 = dbe_->book1D("hist_EB_SC_Size_","# Super Clusters in Barrel",
00085 hist_bins_Size_,hist_min_Size_,hist_max_Size_);
00086 hist_EE_SC_Size_
00087 = dbe_->book1D("hist_EE_SC_Size_","# Super Clusters in Endcap",
00088 hist_bins_Size_,hist_min_Size_,hist_max_Size_);
00089
00090 hist_EB_SC_NumBC_
00091 = dbe_->book1D("hist_EB_SC_NumBC_","# of Basic Clusters in Super Clusters in Barrel",
00092 hist_bins_NumBC_,hist_min_NumBC_,hist_max_NumBC_);
00093 hist_EE_SC_NumBC_
00094 = dbe_->book1D("hist_EE_SC_NumBC_","# of Basic Clusters in Super Clusters in Endcap",
00095 hist_bins_NumBC_,hist_min_NumBC_,hist_max_NumBC_);
00096
00097 hist_EB_SC_ET_
00098 = dbe_->book1D("hist_EB_SC_ET_","ET of Super Clusters in Barrel",
00099 hist_bins_ET_,hist_min_ET_,hist_max_ET_);
00100 hist_EE_SC_ET_
00101 = dbe_->book1D("hist_EE_SC_ET_","ET of Super Clusters in Endcap",
00102 hist_bins_ET_,hist_min_ET_,hist_max_ET_);
00103
00104 hist_EB_SC_Eta_
00105 = dbe_->book1D("hist_EB_SC_Eta_","Eta of Super Clusters in Barrel",
00106 hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_);
00107 hist_EE_SC_Eta_
00108 = dbe_->book1D("hist_EE_SC_Eta_","Eta of Super Clusters in Endcap",
00109 hist_bins_Eta_,hist_min_Eta_,hist_max_Eta_);
00110
00111 hist_EB_SC_Phi_
00112 = dbe_->book1D("hist_EB_SC_Phi_","Phi of Super Clusters in Barrel",
00113 hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_);
00114 hist_EE_SC_Phi_
00115 = dbe_->book1D("hist_EE_SC_Phi_","Phi of Super Clusters in Endcap",
00116 hist_bins_Phi_,hist_min_Phi_,hist_max_Phi_);
00117
00118 hist_EB_SC_S1toS9_
00119 = dbe_->book1D("hist_EB_SC_S1toS9_","S1/S9 of Super Clusters in Barrel",
00120 hist_bins_S1toS9_,hist_min_S1toS9_,hist_max_S1toS9_);
00121 hist_EE_SC_S1toS9_
00122 = dbe_->book1D("hist_EE_SC_S1toS9_","S1/S9 of Super Clusters in Endcap",
00123 hist_bins_S1toS9_,hist_min_S1toS9_,hist_max_S1toS9_);
00124
00125 hist_EB_SC_S25toE_
00126 = dbe_->book1D("hist_EB_SC_S25toE_","S25/E of Super Clusters in Barrel",
00127 hist_bins_S25toE_,hist_min_S25toE_,hist_max_S25toE_);
00128 hist_EE_SC_S25toE_
00129 = dbe_->book1D("hist_EE_SC_S25toE_","S25/E of Super Clusters in Endcap",
00130 hist_bins_S25toE_,hist_min_S25toE_,hist_max_S25toE_);
00131
00132 hist_EB_SC_EToverTruth_
00133 = dbe_->book1D("hist_EB_SC_EToverTruth_","ET/True ET of Super Clusters in Barrel",
00134 hist_bins_EToverTruth_,hist_min_EToverTruth_,hist_max_EToverTruth_);
00135 hist_EE_SC_EToverTruth_
00136 = dbe_->book1D("hist_EE_SC_EToverTruth_","ET/True ET of Super Clusters in Endcap",
00137 hist_bins_EToverTruth_,hist_min_EToverTruth_,hist_max_EToverTruth_);
00138
00139 hist_EB_SC_deltaEta_
00140 = dbe_->book1D("hist_EB_SC_deltaEta_","Eta-True Eta of Super Clusters in Barrel",
00141 hist_bins_deltaEta_,hist_min_deltaEta_,hist_max_deltaEta_);
00142 hist_EE_SC_deltaEta_
00143 = dbe_->book1D("hist_EE_SC_deltaEta_","Eta-True Eta of Super Clusters in Endcap",
00144 hist_bins_deltaEta_,hist_min_deltaEta_,hist_max_deltaEta_);
00145 }
00146
00147
00148 void EgammaSuperClusters::analyze( const edm::Event& evt, const edm::EventSetup& es )
00149 {
00150 edm::Handle<reco::SuperClusterCollection> pBarrelSuperClusters;
00151 evt.getByLabel(barrelSuperClusterCollection_, pBarrelSuperClusters);
00152 if (!pBarrelSuperClusters.isValid()) {
00153 edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label "
00154 << barrelSuperClusterCollection_.label();
00155 }
00156
00157 EcalClusterLazyTools lazyTool( evt, es, barrelRecHitCollection_, endcapRecHitCollection_ );
00158
00159 const reco::SuperClusterCollection* barrelSuperClusters = pBarrelSuperClusters.product();
00160 hist_EB_SC_Size_->Fill(barrelSuperClusters->size());
00161
00162 for(reco::SuperClusterCollection::const_iterator aClus = barrelSuperClusters->begin();
00163 aClus != barrelSuperClusters->end(); aClus++)
00164 {
00165 hist_EB_SC_NumBC_->Fill(aClus->clustersSize());
00166 hist_EB_SC_ET_->Fill(aClus->energy()/std::cosh(aClus->position().eta()));
00167 hist_EB_SC_Eta_->Fill(aClus->position().eta());
00168 hist_EB_SC_Phi_->Fill(aClus->position().phi());
00169
00170 const reco::BasicClusterRef seed = aClus->seed();
00171 hist_EB_SC_S1toS9_->Fill( lazyTool.eMax( *seed ) / lazyTool.e3x3( *seed ) );
00172 hist_EB_SC_S25toE_->Fill( lazyTool.e5x5( *seed ) / aClus->energy() );
00173 }
00174
00175 edm::Handle<reco::SuperClusterCollection> pEndcapSuperClusters;
00176 evt.getByLabel(endcapSuperClusterCollection_, pEndcapSuperClusters);
00177 if (!pEndcapSuperClusters.isValid()) {
00178 edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label "
00179 << endcapSuperClusterCollection_.label();
00180 }
00181
00182 const reco::SuperClusterCollection* endcapSuperClusters = pEndcapSuperClusters.product();
00183 hist_EE_SC_Size_->Fill(endcapSuperClusters->size());
00184
00185 for(reco::SuperClusterCollection::const_iterator aClus = endcapSuperClusters->begin();
00186 aClus != endcapSuperClusters->end(); aClus++)
00187 {
00188 hist_EE_SC_NumBC_->Fill(aClus->clustersSize());
00189 hist_EE_SC_ET_->Fill(aClus->energy()/std::cosh(aClus->position().eta()));
00190 hist_EE_SC_Eta_->Fill(aClus->position().eta());
00191 hist_EE_SC_Phi_->Fill(aClus->position().phi());
00192
00193 const reco::BasicClusterRef seed = aClus->seed();
00194 hist_EE_SC_S1toS9_->Fill( lazyTool.eMax( *seed ) / lazyTool.e3x3( *seed ) );
00195 hist_EE_SC_S25toE_->Fill( lazyTool.e5x5( *seed ) / aClus->energy() );
00196 }
00197
00198 edm::Handle<edm::HepMCProduct> pMCTruth ;
00199 evt.getByLabel(MCTruthCollection_, pMCTruth);
00200 if (!pMCTruth.isValid()) {
00201 edm::LogError("EgammaSuperClusters") << "Error! can't get collection with label "
00202 << MCTruthCollection_.label();
00203 }
00204
00205 const HepMC::GenEvent* genEvent = pMCTruth->GetEvent();
00206 for(HepMC::GenEvent::particle_const_iterator currentParticle = genEvent->particles_begin();
00207 currentParticle != genEvent->particles_end(); currentParticle++ )
00208 {
00209 if((*currentParticle)->status()==1)
00210 {
00211 HepMC::FourVector vtx = (*currentParticle)->production_vertex()->position();
00212 double phiTrue = (*currentParticle)->momentum().phi();
00213 double etaTrue = ecalEta((*currentParticle)->momentum().eta(), vtx.z()/10., vtx.perp()/10.);
00214 double etTrue = (*currentParticle)->momentum().e()/cosh(etaTrue);
00215
00216 if(std::fabs(etaTrue) < 1.479)
00217 {
00218 {
00219 double etaCurrent, etaFound = 0;
00220 double phiCurrent;
00221 double etCurrent, etFound = 0;
00222
00223 double closestParticleDistance = 999;
00224
00225 for(reco::SuperClusterCollection::const_iterator aClus = barrelSuperClusters->begin();
00226 aClus != barrelSuperClusters->end(); aClus++)
00227 {
00228 etaCurrent = aClus->position().eta();
00229 phiCurrent = aClus->position().phi();
00230 etCurrent = aClus->energy()/std::cosh(etaCurrent);
00231
00232 double deltaR = std::sqrt(std::pow(etaCurrent-etaTrue,2)+std::pow(phiCurrent-phiTrue,2));
00233
00234 if(deltaR < closestParticleDistance)
00235 {
00236 etFound = etCurrent;
00237 etaFound = etaCurrent;
00238 closestParticleDistance = deltaR;
00239 }
00240 }
00241
00242 if(closestParticleDistance < 0.3)
00243 {
00244 hist_EB_SC_EToverTruth_->Fill(etFound/etTrue);
00245 hist_EB_SC_deltaEta_->Fill(etaFound-etaTrue);
00246 }
00247 }
00248 }
00249 else
00250 {
00251 double etaCurrent, etaFound = 0;
00252 double phiCurrent;
00253 double etCurrent, etFound = 0;
00254
00255 double closestParticleDistance = 999;
00256
00257 for(reco::SuperClusterCollection::const_iterator aClus = endcapSuperClusters->begin();
00258 aClus != endcapSuperClusters->end(); aClus++)
00259 {
00260 etaCurrent = aClus->position().eta();
00261 phiCurrent = aClus->position().phi();
00262 etCurrent = aClus->energy()/std::cosh(etaCurrent);
00263
00264 double deltaR = std::sqrt(std::pow(etaCurrent-etaTrue,2)+std::pow(phiCurrent-phiTrue,2));
00265
00266 if(deltaR < closestParticleDistance)
00267 {
00268 etFound = etCurrent;
00269 etaFound = etaCurrent;
00270 closestParticleDistance = deltaR;
00271 }
00272 }
00273
00274 if(closestParticleDistance < 0.3)
00275 {
00276 hist_EE_SC_EToverTruth_->Fill(etFound/etTrue);
00277 hist_EE_SC_deltaEta_->Fill(etaFound-etaTrue);
00278 }
00279 }
00280 }
00281 }
00282 }
00283
00284 void EgammaSuperClusters::endJob()
00285 {
00286 if (outputFile_.size() != 0 && dbe_) dbe_->save(outputFile_);
00287 }
00288
00289 float EgammaSuperClusters::ecalEta(float EtaParticle , float Zvertex, float plane_Radius)
00290 {
00291 const float R_ECAL = 136.5;
00292 const float Z_Endcap = 328.0;
00293 const float etaBarrelEndcap = 1.479;
00294
00295 if(EtaParticle != 0.)
00296 {
00297 float Theta = 0.0 ;
00298 float ZEcal = (R_ECAL-plane_Radius)*sinh(EtaParticle)+Zvertex;
00299
00300 if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
00301 if(Theta<0.0) Theta = Theta+Geom::pi() ;
00302
00303 float ETA = - log(tan(0.5*Theta));
00304
00305 if( fabs(ETA) > etaBarrelEndcap )
00306 {
00307 float Zend = Z_Endcap ;
00308 if(EtaParticle<0.0 ) Zend = -Zend ;
00309 float Zlen = Zend - Zvertex ;
00310 float RR = Zlen/sinh(EtaParticle);
00311 Theta = atan((RR+plane_Radius)/Zend);
00312 if(Theta<0.0) Theta = Theta+Geom::pi() ;
00313 ETA = - log(tan(0.5*Theta));
00314 }
00315
00316 return ETA;
00317 }
00318 else
00319 {
00320 edm::LogWarning("") << "[EgammaSuperClusters::ecalEta] Warning: Eta equals to zero, not correcting" ;
00321 return EtaParticle;
00322 }
00323 }
1.5.4