Reference Manual

Go to the documentation of this file.
00001 // system include files
00002 #include <vector>
00003 #include <memory>
00004 
00005 // user include files
00006 #include "FWCore/Framework/interface/Frameworkfwd.h"
00007 #include "FWCore/Framework/interface/EDAnalyzer.h"
00008 
00009 #include "FWCore/Framework/interface/Event.h"
00010 #include "FWCore/Framework/interface/EventSetup.h"
00011 #include "DataFormats/Common/interface/Handle.h"
00012 #include "FWCore/Framework/interface/ESHandle.h"
00013 #include "FWCore/Framework/interface/MakerMacros.h"
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015 #include "FWCore/Utilities/interface/Exception.h"
00016 
00017 // Reconstruction Classes
00018 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00019 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00020 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00021 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00022 #include "DataFormats/EgammaReco/interface/PreshowerClusterFwd.h"
00023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00024 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00025 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00026 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00027 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
00028 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
00029 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
00030 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00031 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00032 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00033 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00034 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00035 #include "DataFormats/EcalDetId/interface/ESDetId.h"
00036 #include "CondFormats/DataRecord/interface/ESGainRcd.h"
00037 #include "CondFormats/DataRecord/interface/ESMIPToGeVConstantRcd.h"
00038 #include "CondFormats/DataRecord/interface/ESEEIntercalibConstantsRcd.h"
00039 #include "CondFormats/DataRecord/interface/ESChannelStatusRcd.h"
00040 #include <fstream>
00041 
00042 #include "RecoEcal/EgammaClusterProducers/interface/PreshowerClusterProducer.h"
00043 
00044 
00045 using namespace edm;
00046 using namespace std;
00047 
00048 
00049 PreshowerClusterProducer::PreshowerClusterProducer(const edm::ParameterSet& ps) {
00050 
00051   // use configuration file to setup input/output collection names 
00052   preshHitProducer_   = ps.getParameter<edm::InputTag>("preshRecHitProducer");
00053 
00054   // Name of a SuperClusterCollection to make associations:
00055   endcapSClusterProducer_   = ps.getParameter<edm::InputTag>("endcapSClusterProducer");
00056 
00057   // Output collections:
00058   preshClusterCollectionX_ = ps.getParameter<std::string>("preshClusterCollectionX");
00059   preshClusterCollectionY_ = ps.getParameter<std::string>("preshClusterCollectionY");
00060   preshNclust_             = ps.getParameter<int>("preshNclust");
00061 
00062   etThresh_ =  ps.getParameter<double>("etThresh");
00063 
00064   assocSClusterCollection_ = ps.getParameter<std::string>("assocSClusterCollection");
00065 
00066   produces< reco::PreshowerClusterCollection >(preshClusterCollectionX_);
00067   produces< reco::PreshowerClusterCollection >(preshClusterCollectionY_);
00068   produces< reco::SuperClusterCollection >(assocSClusterCollection_);
00069 
00070   float preshStripECut = ps.getParameter<double>("preshStripEnergyCut");
00071     int preshSeededNst = ps.getParameter<int>("preshSeededNstrip");
00072   preshClustECut = ps.getParameter<double>("preshClusterEnergyCut");
00073 
00074   // The debug level
00075   std::string debugString = ps.getParameter<std::string>("debugLevel");
00076   if      (debugString == "DEBUG")   debugL = PreshowerClusterAlgo::pDEBUG;
00077   else if (debugString == "INFO")    debugL = PreshowerClusterAlgo::pINFO;
00078   else                               debugL = PreshowerClusterAlgo::pERROR;
00079 
00080   presh_algo = new PreshowerClusterAlgo(preshStripECut,preshClustECut,preshSeededNst,debugL);
00081 
00082   nEvt_ = 0;  
00083 
00084 }
00085 
00086 
00087 PreshowerClusterProducer::~PreshowerClusterProducer() {
00088    delete presh_algo;
00089 }
00090 
00091 
00092 void PreshowerClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
00093 
00094   edm::Handle< EcalRecHitCollection >   pRecHits;
00095   edm::Handle< reco::SuperClusterCollection > pSuperClusters;
00096 
00097   // get the ECAL geometry:
00098   edm::ESHandle<CaloGeometry> geoHandle;
00099   es.get<CaloGeometryRecord>().get(geoHandle);
00100 
00101   // retrieve ES-EE intercalibration constants and channel status
00102   set(es);
00103   const ESChannelStatus *channelStatus = esChannelStatus_.product();
00104 
00105   const CaloSubdetectorGeometry *geometry = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00106   const CaloSubdetectorGeometry *& geometry_p = geometry;
00107 
00108    // create auto_ptr to a PreshowerClusterCollection
00109    std::auto_ptr< reco::PreshowerClusterCollection > clusters_p1(new reco::PreshowerClusterCollection);
00110    std::auto_ptr< reco::PreshowerClusterCollection > clusters_p2(new reco::PreshowerClusterCollection);
00111    // create new collection of corrected super clusters
00112    std::auto_ptr< reco::SuperClusterCollection > superclusters_p(new reco::SuperClusterCollection);
00113 
00114    CaloSubdetectorTopology * topology_p=0;
00115    if (geometry)
00116      topology_p  = new EcalPreshowerTopology(geoHandle);
00117 
00118    // fetch the product (pSuperClusters)
00119    evt.getByLabel(endcapSClusterProducer_, pSuperClusters);   
00120    const reco::SuperClusterCollection* SClusts = pSuperClusters.product();
00121    if ( debugL <= PreshowerClusterAlgo::pINFO ) std::cout <<"### Total # Endcap Superclusters: " << SClusts->size() << std::endl;
00122    
00123    // fetch the product (RecHits)
00124    evt.getByLabel( preshHitProducer_, pRecHits);
00125    // pointer to the object in the product
00126    const EcalRecHitCollection* rechits = pRecHits.product(); // EcalRecHitCollection hit_collection = *rhcHandle;
00127    if ( debugL == PreshowerClusterAlgo::pDEBUG ) std::cout << "PreshowerClusterProducerInfo: ### Total # of preshower RecHits: " 
00128                                                            << rechits->size() << std::endl;
00129 
00130   // make the map of rechits:
00131    std::map<DetId, EcalRecHit> rechits_map;
00132    EcalRecHitCollection::const_iterator it;
00133    for (it = rechits->begin(); it != rechits->end(); it++) {
00134      // remove bad ES rechits
00135      if (it->recoFlag()==14 || (it->recoFlag()<=10 && it->recoFlag()>=5)) continue;
00136      //Make the map of DetID, EcalRecHit pairs
00137      rechits_map.insert(std::make_pair(it->id(), *it));   
00138    }
00139    // The set of used DetID's for a given event:
00140    std::set<DetId> used_strips;
00141    used_strips.clear();
00142    
00143   if ( debugL <= PreshowerClusterAlgo::pINFO ) std::cout << "PreshowerClusterProducerInfo: ### rechits_map of size " 
00144                                                          << rechits_map.size() <<" was created!" << std::endl;   
00145   
00146   reco::PreshowerClusterCollection clusters1, clusters2;   // output collection of corrected PCs
00147   reco::SuperClusterCollection new_SC; // output collection of corrected SCs
00148 
00149   if ( debugL == PreshowerClusterAlgo::pDEBUG ) std::cout << " Making a cycle over Superclusters ..." << std::endl; 
00150   //make cycle over super clusters
00151   reco::SuperClusterCollection::const_iterator it_super;
00152   int isc = 0;
00153   for ( it_super=SClusts->begin();  it_super!=SClusts->end(); it_super++ ) {     
00154        float e1=0;
00155        float e2=0;
00156        float deltaE=0;
00157        reco::CaloClusterPtrVector new_BC; 
00158        ++isc;
00159 
00160        if ( debugL <= PreshowerClusterAlgo::pINFO ) std::cout << " superE = " << it_super->energy() << " superETA = " << it_super->eta() 
00161                                                        << " superPHI = " << it_super->phi() << std::endl;
00162        if ( debugL == PreshowerClusterAlgo::pINFO ) std::cout << " This SC contains " << it_super->clustersSize() << " BCs" << std::endl;
00163 
00164        int nBC = 0;
00165        int condP1 = 1; // 0: dead channel; 1: active channel
00166        int condP2 = 1;
00167        reco::CaloCluster_iterator bc_iter = it_super->clustersBegin();
00168        for ( ; bc_iter !=it_super->clustersEnd(); ++bc_iter ) {  
00169          if (geometry) {
00170            
00171            // Get strip position at intersection point of the line EE - Vertex:
00172            double X = (*bc_iter)->x();
00173            double Y = (*bc_iter)->y();
00174            double Z = (*bc_iter)->z();        
00175            const GlobalPoint point(X,Y,Z);    
00176            
00177            DetId tmp1 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 1);
00178            DetId tmp2 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 2);
00179            ESDetId strip1 = (tmp1 == DetId(0)) ? ESDetId(0) : ESDetId(tmp1);
00180            ESDetId strip2 = (tmp2 == DetId(0)) ? ESDetId(0) : ESDetId(tmp2);     
00181            
00182            if (nBC == 0) {
00183              if (strip1 != ESDetId(0) && strip2 != ESDetId(0)) {
00184                ESChannelStatusMap::const_iterator status_p1 = channelStatus->getMap().find(strip1);
00185                ESChannelStatusMap::const_iterator status_p2 = channelStatus->getMap().find(strip2);
00186                if (status_p1->getStatusCode() == 1) condP1 = 0;
00187                if (status_p2->getStatusCode() == 1) condP2 = 0;
00188              } else if (strip1 == ESDetId(0))
00189                condP1 = 0;
00190              else if (strip2 == ESDetId(0))
00191                condP2 = 0;
00192            }
00193            
00194            if ( debugL <= PreshowerClusterAlgo::pINFO ) {
00195              if ( strip1 != ESDetId(0) && strip2 != ESDetId(0) ) {
00196                std::cout << " Intersected preshower strips are: " << std::endl;
00197                std::cout << strip1 << std::endl;
00198                std::cout << strip2 << std::endl;
00199              } else if ( strip1 == ESDetId(0) )
00200                std::cout << " No intersected strip in plane 1 " << std::endl;
00201              else if ( strip2 == ESDetId(0) )
00202                std::cout << " No intersected strip in plane 2 " << std::endl;
00203            }        
00204            
00205            // Get a vector of ES clusters (found by the PreshSeeded algorithm) associated with a given EE basic cluster.           
00206            for (int i=0; i<preshNclust_; i++) {
00207              reco::PreshowerCluster cl1 = presh_algo->makeOneCluster(strip1,&used_strips,&rechits_map,geometry_p,topology_p);   
00208              cl1.setBCRef(*bc_iter);
00209              if (cl1.energy() > preshClustECut) {
00210                clusters1.push_back(cl1);
00211                e1 += cl1.energy();       
00212              }
00213              reco::PreshowerCluster cl2 = presh_algo->makeOneCluster(strip2,&used_strips,&rechits_map,geometry_p,topology_p); 
00214              cl2.setBCRef(*bc_iter);
00215              
00216              if ( cl2.energy() > preshClustECut) {
00217                clusters2.push_back(cl2);
00218                e2 += cl2.energy();
00219              }                                 
00220            } // end of cycle over ES clusters
00221          }
00222 
00223          new_BC.push_back(*bc_iter);
00224          nBC++;
00225        }  // end of cycle over BCs
00226        
00227        if ( debugL <= PreshowerClusterAlgo::pINFO ) std::cout << " For SC #" << isc-1 << ", containing " << it_super->clustersSize() 
00228                                                               << " basic clusters, PreshowerClusterAlgo made " 
00229                                                               << clusters1.size() << " in X plane and " << clusters2.size() 
00230                                                               << " in Y plane " << " preshower clusters " << std::endl;
00231        
00232        // update energy of the SuperCluster    
00233        if(e1+e2 > 1.0e-10) {
00234          // GeV to #MIPs
00235          e1 = e1 / mip_;
00236          e2 = e2 / mip_;
00237 
00238          if (condP1 == 1 && condP2 == 1) deltaE = gamma0_*(e1 + alpha0_*e2);       
00239          else if (condP1 == 1 && condP2 == 0) deltaE = gamma1_*(e1 + alpha1_*e2);       
00240          else if (condP1 == 0 && condP2 == 1) deltaE = gamma2_*(e1 + alpha2_*e2);       
00241          else if (condP1 == 0 && condP2 == 0) deltaE = gamma3_*(e1 + alpha3_*e2);       
00242        }
00243        
00244        //corrected Energy
00245        float E = it_super->energy() + deltaE;
00246        
00247        if (debugL == PreshowerClusterAlgo::pDEBUG) std::cout << " Creating corrected SC " << std::endl;
00248        reco::SuperCluster sc(E, it_super->position(), it_super->seed(), new_BC, deltaE);
00249        if (condP1 == 1 && condP2 == 1) sc.setPreshowerPlanesStatus(0);
00250        else if (condP1 == 1 && condP2 == 0) sc.setPreshowerPlanesStatus(1);
00251        else if (condP1 == 0 && condP2 == 1) sc.setPreshowerPlanesStatus(2);
00252        else if (condP1 == 0 && condP2 == 0) sc.setPreshowerPlanesStatus(3);
00253        
00254        if (sc.energy()*sin(sc.position().theta())>etThresh_)
00255          new_SC.push_back(sc);
00256        if (debugL <= PreshowerClusterAlgo::pINFO) std::cout << " SuperClusters energies: new E = " << sc.energy() 
00257                                         << " vs. old E =" << it_super->energy() << std::endl;
00258 
00259    } // end of cycle over SCs
00260 
00261    // copy the preshower clusters into collections and put in the Event:
00262    clusters_p1->assign(clusters1.begin(), clusters1.end());
00263    clusters_p2->assign(clusters2.begin(), clusters2.end());
00264    // put collection of preshower clusters to the event
00265    evt.put( clusters_p1, preshClusterCollectionX_ );
00266    evt.put( clusters_p2, preshClusterCollectionY_ );
00267    if ( debugL <= PreshowerClusterAlgo::pINFO ) std::cout << "Preshower clusters added to the event" << std::endl;
00268 
00269    // put collection of corrected super clusters to the event
00270    superclusters_p->assign(new_SC.begin(), new_SC.end());
00271    evt.put(superclusters_p, assocSClusterCollection_);
00272    if ( debugL <= PreshowerClusterAlgo::pINFO ) std::cout << "Corrected SClusters added to the event" << std::endl;
00273 
00274    if (topology_p)
00275      delete topology_p;
00276 
00277    nEvt_++;
00278 
00279 }
00280 
00281 void PreshowerClusterProducer::set(const edm::EventSetup& es) {
00282 
00283   es.get<ESGainRcd>().get(esgain_);
00284   const ESGain *gain = esgain_.product();
00285 
00286   double ESGain = gain->getESGain();
00287 
00288   es.get<ESMIPToGeVConstantRcd>().get(esMIPToGeV_);
00289   const ESMIPToGeVConstant *mipToGeV = esMIPToGeV_.product();
00290 
00291   mip_ = (ESGain == 1) ? mipToGeV->getESValueLow() : mipToGeV->getESValueHigh(); 
00292 
00293   es.get<ESChannelStatusRcd>().get(esChannelStatus_);
00294 
00295   es.get<ESEEIntercalibConstantsRcd>().get(esEEInterCalib_);
00296   const ESEEIntercalibConstants *esEEInterCalib = esEEInterCalib_.product();
00297 
00298   // both planes work
00299   gamma0_ = (ESGain == 1) ? esEEInterCalib->getGammaLow0() : esEEInterCalib->getGammaHigh0();
00300   alpha0_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow0() : esEEInterCalib->getAlphaHigh0();
00301 
00302   // only first plane works 
00303   gamma1_ = (ESGain == 1) ? esEEInterCalib->getGammaLow1() : esEEInterCalib->getGammaHigh1();
00304   alpha1_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow1() : esEEInterCalib->getAlphaHigh1();
00305 
00306   // only second plane works
00307   gamma2_ = (ESGain == 1) ? esEEInterCalib->getGammaLow2() : esEEInterCalib->getGammaHigh2();
00308   alpha2_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow2() : esEEInterCalib->getAlphaHigh2();
00309 
00310   // both planes do not work
00311   gamma3_ = (ESGain == 1) ? esEEInterCalib->getGammaLow3() : esEEInterCalib->getGammaHigh3();
00312   alpha3_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow3() : esEEInterCalib->getAlphaHigh3();
00313 
00314 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoEcal/EgammaClusterProducers/src/PreshowerClusterProducer.cc
