00001
00002 #include <iostream>
00003 #include <vector>
00004 #include <memory>
00005
00006
00007 #include "FWCore/Framework/interface/Event.h"
00008 #include "FWCore/Framework/interface/EventSetup.h"
00009 #include "DataFormats/Common/interface/Handle.h"
00010 #include "FWCore/Framework/interface/ESHandle.h"
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012 #include "FWCore/Utilities/interface/Exception.h"
00013
00014
00015 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00016 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00018 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00019 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00020 #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
00021
00022
00023 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00024 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00025 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00026 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00027 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
00028 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
00029
00030
00031 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
00032 #include "RecoEcal/EgammaCoreTools/interface/ClusterShapeAlgo.h"
00033 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
00034 #include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"
00035
00036
00037 #include "RecoEcal/EgammaClusterProducers/interface/IslandClusterProducer.h"
00038
00039
00040 IslandClusterProducer::IslandClusterProducer(const edm::ParameterSet& ps)
00041 {
00042
00043 std::string verbosityString = ps.getParameter<std::string>("VerbosityLevel");
00044 if (verbosityString == "DEBUG") verbosity = IslandClusterAlgo::pDEBUG;
00045 else if (verbosityString == "WARNING") verbosity = IslandClusterAlgo::pWARNING;
00046 else if (verbosityString == "INFO") verbosity = IslandClusterAlgo::pINFO;
00047 else verbosity = IslandClusterAlgo::pERROR;
00048
00049
00050 barrelHitProducer_ = ps.getParameter<std::string>("barrelHitProducer");
00051 endcapHitProducer_ = ps.getParameter<std::string>("endcapHitProducer");
00052 barrelHitCollection_ = ps.getParameter<std::string>("barrelHitCollection");
00053 endcapHitCollection_ = ps.getParameter<std::string>("endcapHitCollection");
00054
00055
00056 barrelClusterCollection_ = ps.getParameter<std::string>("barrelClusterCollection");
00057 endcapClusterCollection_ = ps.getParameter<std::string>("endcapClusterCollection");
00058
00059
00060 double barrelSeedThreshold = ps.getParameter<double>("IslandBarrelSeedThr");
00061 double endcapSeedThreshold = ps.getParameter<double>("IslandEndcapSeedThr");
00062
00063
00064 std::map<std::string,double> providedParameters;
00065 providedParameters.insert(std::make_pair("LogWeighted",ps.getParameter<bool>("posCalc_logweight")));
00066 providedParameters.insert(std::make_pair("T0_barl",ps.getParameter<double>("posCalc_t0_barl")));
00067 providedParameters.insert(std::make_pair("T0_endc",ps.getParameter<double>("posCalc_t0_endc")));
00068 providedParameters.insert(std::make_pair("T0_endcPresh",ps.getParameter<double>("posCalc_t0_endcPresh")));
00069 providedParameters.insert(std::make_pair("W0",ps.getParameter<double>("posCalc_w0")));
00070 providedParameters.insert(std::make_pair("X0",ps.getParameter<double>("posCalc_x0")));
00071 posCalculator_ = PositionCalc(providedParameters);
00072 shapeAlgo_ = ClusterShapeAlgo(providedParameters);
00073
00074 clustershapecollectionEB_ = ps.getParameter<std::string>("clustershapecollectionEB");
00075 clustershapecollectionEE_ = ps.getParameter<std::string>("clustershapecollectionEE");
00076
00077
00078 barrelClusterShapeAssociation_ = ps.getParameter<std::string>("barrelShapeAssociation");
00079 endcapClusterShapeAssociation_ = ps.getParameter<std::string>("endcapShapeAssociation");
00080
00081
00082
00083 produces< reco::ClusterShapeCollection>(clustershapecollectionEE_);
00084 produces< reco::BasicClusterCollection >(endcapClusterCollection_);
00085 produces< reco::ClusterShapeCollection>(clustershapecollectionEB_);
00086 produces< reco::BasicClusterCollection >(barrelClusterCollection_);
00087 produces< reco::BasicClusterShapeAssociationCollection >(barrelClusterShapeAssociation_);
00088 produces< reco::BasicClusterShapeAssociationCollection >(endcapClusterShapeAssociation_);
00089
00090 island_p = new IslandClusterAlgo(barrelSeedThreshold, endcapSeedThreshold, posCalculator_,verbosity);
00091
00092 nEvt_ = 0;
00093 }
00094
00095
00096 IslandClusterProducer::~IslandClusterProducer()
00097 {
00098 delete island_p;
00099 }
00100
00101
00102 void IslandClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es)
00103 {
00104 clusterizeECALPart(evt, es, endcapHitProducer_, endcapHitCollection_, endcapClusterCollection_, endcapClusterShapeAssociation_, IslandClusterAlgo::endcap);
00105 clusterizeECALPart(evt, es, barrelHitProducer_, barrelHitCollection_, barrelClusterCollection_, barrelClusterShapeAssociation_, IslandClusterAlgo::barrel);
00106 nEvt_++;
00107 }
00108
00109
00110 const EcalRecHitCollection * IslandClusterProducer::getCollection(edm::Event& evt,
00111 const std::string& hitProducer_,
00112 const std::string& hitCollection_)
00113 {
00114 edm::Handle<EcalRecHitCollection> rhcHandle;
00115 evt.getByLabel(hitProducer_, hitCollection_, rhcHandle);
00116 if (!(rhcHandle.isValid()))
00117 {
00118 std::cout << "could not get a handle on the EcalRecHitCollection!" << std::endl;
00119 edm::LogError("IslandClusterProducerError") << "Error! can't get the product " << hitCollection_.c_str() ;
00120 return 0;
00121 } else {
00122 return rhcHandle.product();
00123 }
00124 }
00125
00126 void IslandClusterProducer::clusterizeECALPart(edm::Event &evt, const edm::EventSetup &es,
00127 const std::string& hitProducer,
00128 const std::string& hitCollection,
00129 const std::string& clusterCollection,
00130 const std::string& clusterShapeAssociation,
00131 const IslandClusterAlgo::EcalPart& ecalPart)
00132 {
00133
00134 const EcalRecHitCollection *hitCollection_p = getCollection(evt, hitProducer, hitCollection);
00135
00136
00137 edm::ESHandle<CaloGeometry> geoHandle;
00138 es.get<CaloGeometryRecord>().get(geoHandle);
00139
00140 const CaloSubdetectorGeometry *geometry_p;
00141 CaloSubdetectorTopology *topology_p;
00142
00143 std::string clustershapetag;
00144 if (ecalPart == IslandClusterAlgo::barrel)
00145 {
00146 geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
00147 topology_p = new EcalBarrelTopology(geoHandle);
00148 }
00149 else
00150 {
00151 geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
00152 topology_p = new EcalEndcapTopology(geoHandle);
00153 }
00154
00155 const CaloSubdetectorGeometry *geometryES_p;
00156 geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00157
00158
00159 reco::BasicClusterCollection clusters;
00160 clusters = island_p->makeClusters(hitCollection_p, geometry_p, topology_p, geometryES_p, ecalPart);
00161
00162
00163 std::vector <reco::ClusterShape> ClusVec;
00164 for (int erg=0;erg<int(clusters.size());++erg){
00165 reco::ClusterShape TestShape = shapeAlgo_.Calculate(clusters[erg],hitCollection_p,geometry_p,topology_p);
00166 ClusVec.push_back(TestShape);
00167 }
00168
00169
00170 std::auto_ptr< reco::ClusterShapeCollection> clustersshapes_p(new reco::ClusterShapeCollection);
00171 clustersshapes_p->assign(ClusVec.begin(), ClusVec.end());
00172 edm::OrphanHandle<reco::ClusterShapeCollection> clusHandle;
00173 if (ecalPart == IslandClusterAlgo::barrel)
00174 clusHandle= evt.put(clustersshapes_p, clustershapecollectionEB_);
00175 else
00176 clusHandle= evt.put(clustersshapes_p, clustershapecollectionEE_);
00177
00178
00179 std::auto_ptr< reco::BasicClusterCollection > clusters_p(new reco::BasicClusterCollection);
00180 clusters_p->assign(clusters.begin(), clusters.end());
00181 edm::OrphanHandle<reco::BasicClusterCollection> bccHandle;
00182 if (ecalPart == IslandClusterAlgo::barrel)
00183 bccHandle = evt.put(clusters_p, barrelClusterCollection_);
00184 else
00185 bccHandle = evt.put(clusters_p, endcapClusterCollection_);
00186
00187
00188
00189 std::auto_ptr<reco::BasicClusterShapeAssociationCollection> shapeAssocs_p(new reco::BasicClusterShapeAssociationCollection);
00190 for (unsigned int i = 0; i < clusHandle->size(); i++){
00191 shapeAssocs_p->insert(edm::Ref<reco::BasicClusterCollection>(bccHandle,i),edm::Ref<reco::ClusterShapeCollection>(clusHandle,i));
00192 }
00193 evt.put(shapeAssocs_p,clusterShapeAssociation);
00194
00195 delete topology_p;
00196 }