Castor/src/CastorTowerProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Castor
0004 // Class:      CastorTowerProducer
0005 //
0006
0007 /**\class CastorTowerProducer CastorTowerProducer.cc RecoLocalCalo/Castor/src/CastorTowerProducer.cc
0008
0009  Description: CastorTower Reconstruction Producer. Produce CastorTowers from CastorCells.
0010  Implementation:
0011 */
0012
0013 //
0014 // Original Author:  Hans Van Haevermaet, Benoit Roland
0015 //         Created:  Wed Jul  9 14:00:40 CEST 2008
0016 //
0017 //
0018
0019 // system include
0020 #include <memory>
0021 #include <vector>
0022 #include <iostream>
0023 #include <TMath.h>
0024 #include <TRandom3.h>
0025
0026 // user include
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/stream/EDProducer.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/Framework/interface/ESHandle.h"
0032 #include "FWCore/Framework/interface/EventSetup.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "FWCore/Utilities/interface/ESGetToken.h"
0035 #include "DataFormats/Math/interface/Point3D.h"
0036
0037 // Castor Object include
0038 #include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
0039 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
0040 #include "DataFormats/CastorReco/interface/CastorTower.h"
0041
0042 // Channel quality
0043 #include "CondFormats/CastorObjects/interface/CastorChannelQuality.h"
0044 #include "CondFormats/CastorObjects/interface/CastorChannelStatus.h"
0045 #include "CondFormats/DataRecord/interface/CastorChannelQualityRcd.h"
0046
0047 //
0048 // class declaration
0049 //
0050
0051 class CastorTowerProducer : public edm::stream::EDProducer<> {
0052 public:
0053   explicit CastorTowerProducer(const edm::ParameterSet&);
0054   ~CastorTowerProducer() override;
0055
0056 private:
0057   void produce(edm::Event&, const edm::EventSetup&) override;
0058   virtual void ComputeTowerVariable(const edm::RefVector<edm::SortedCollection<CastorRecHit> >& usedRecHits,
0059                                     double& Ehot,
0060                                     double& depth);
0061
0062   // member data
0063   typedef math::XYZPointD Point;
0064   typedef ROOT::Math::RhoEtaPhiPoint TowerPoint;
0065   typedef ROOT::Math::RhoZPhiPoint CellPoint;
0066   typedef edm::SortedCollection<CastorRecHit> CastorRecHitCollection;
0067   typedef std::vector<reco::CastorTower> CastorTowerCollection;
0068   typedef edm::RefVector<CastorRecHitCollection> CastorRecHitRefVector;
0069   edm::EDGetTokenT<CastorRecHitCollection> tok_input_;
0070   edm::ESGetToken<CastorChannelQuality, CastorChannelQualityRcd> tok_channelQuality_;
0071   double towercut_;
0072   double mintime_;
0073   double maxtime_;
0074 };
0075
0076 //
0077 // constants, enums and typedefs
0078 //
0079
0080 //
0081 // static data member definitions
0082 //
0083
0084 //
0085 // constructor and destructor
0086 //
0087
0088 CastorTowerProducer::CastorTowerProducer(const edm::ParameterSet& iConfig)
0089     : towercut_(iConfig.getParameter<double>("towercut")),
0090       mintime_(iConfig.getParameter<double>("mintime")),
0091       maxtime_(iConfig.getParameter<double>("maxtime")) {
0092   tok_input_ = consumes<CastorRecHitCollection>(iConfig.getParameter<std::string>("inputprocess"));
0093   tok_channelQuality_ = esConsumes<CastorChannelQuality, CastorChannelQualityRcd>();
0094   //register your products
0095   produces<CastorTowerCollection>();
0096   //now do what ever other initialization is needed
0097 }
0098
0099 CastorTowerProducer::~CastorTowerProducer() {
0100   // do anything here that needs to be done at desctruction time
0101   // (e.g. close files, deallocate resources etc.)
0102 }
0103
0104 //
0105 // member functions
0106 //
0107
0108 // ------------ method called to produce the data  ------------
0109 void CastorTowerProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0110   using namespace edm;
0111   using namespace reco;
0112   using namespace TMath;
0113
0114   // Produce CastorTowers from CastorCells
0115
0116   edm::Handle<CastorRecHitCollection> InputRecHits;
0117   iEvent.getByToken(tok_input_, InputRecHits);
0118
0119   auto OutputTowers = std::make_unique<CastorTowerCollection>();
0120
0121   // get and check input size
0122   int nRecHits = InputRecHits->size();
0123
0124   LogDebug("CastorTowerProducer") << "2. entering CastorTowerProducer" << std::endl;
0125
0126   if (nRecHits == 0)
0127     LogDebug("CastorTowerProducer") << "Warning: You are trying to run the Tower algorithm with 0 input rechits.";
0128
0129   // declare castor array
0130   // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position
0131
0132   double poscastortowerarray[4][16];
0133   double negcastortowerarray[4][16];
0134
0135   CastorRecHitRefVector poscastorusedrechits[16];
0136   CastorRecHitRefVector negcastorusedrechits[16];
0137
0138   // set phi values and everything else to zero
0139   for (int j = 0; j < 16; j++) {
0140     poscastortowerarray[3][j] = -2.94524 + j * 0.3927;
0141     poscastortowerarray[0][j] = 0.;
0142     poscastortowerarray[1][j] = 0.;
0143     poscastortowerarray[2][j] = 0.;
0144
0145     negcastortowerarray[3][j] = -2.94524 + j * 0.3927;
0146     negcastortowerarray[0][j] = 0.;
0147     negcastortowerarray[1][j] = 0.;
0148     negcastortowerarray[2][j] = 0.;
0149   }
0150
0151   // retrieve the channel quality lists from database
0152   edm::ESHandle<CastorChannelQuality> p = iSetup.getHandle(tok_channelQuality_);
0153   std::vector<DetId> channels = p->getAllChannels();
0154
0155   // loop over rechits to build castortowerarray[4][16] and castorusedrechits[16]
0156   for (unsigned int i = 0; i < InputRecHits->size(); i++) {
0157     edm::Ref<CastorRecHitCollection> rechit_p = edm::Ref<CastorRecHitCollection>(InputRecHits, i);
0158
0159     HcalCastorDetId id = rechit_p->id();
0160     DetId genericID = (DetId)id;
0161
0162     // first check if the rechit is in the BAD channel list
0163     bool bad = false;
0164     for (std::vector<DetId>::iterator channel = channels.begin(); channel != channels.end(); channel++) {
0165       if (channel->rawId() == genericID.rawId()) {
0166         // if the rechit is found in the list, set it bad
0167         bad = true;
0168         break;
0169       }
0170     }
0171     // if bad, continue the loop to the next rechit
0172     if (bad)
0173       continue;
0174
0175     double Erechit = rechit_p->energy();
0176     int module = id.module();
0177     int sector = id.sector();
0178     double zrechit = 0;
0179     if (module < 3)
0180       zrechit = -14390 - 24.75 - 49.5 * (module - 1);
0181     if (module > 2)
0182       zrechit = -14390 - 99 - 49.5 - 99 * (module - 3);
0183     double phirechit = -100;
0184     if (sector < 9)
0185       phirechit = 0.19635 + (sector - 1) * 0.3927;
0186     if (sector > 8)
0187       phirechit = -2.94524 + (sector - 9) * 0.3927;
0188
0189     // add time conditions for the rechit
0190     if (rechit_p->time() > mintime_ && rechit_p->time() < maxtime_) {
0191       // loop over the 16 towers possibilities
0192       for (int j = 0; j < 16; j++) {
0193         // phi matching condition
0194         if (TMath::Abs(phirechit - poscastortowerarray[3][j]) < 0.0001) {
0195           // condition over rechit z value
0196           if (zrechit > 0.) {
0197             poscastortowerarray[0][j] += Erechit;
0198             if (module < 3) {
0199               poscastortowerarray[1][j] += Erechit;
0200             } else {
0201               poscastortowerarray[2][j] += Erechit;
0202             }
0203             poscastorusedrechits[j].push_back(rechit_p);
0204           } else {
0205             negcastortowerarray[0][j] += Erechit;
0206             if (module < 3) {
0207               negcastortowerarray[1][j] += Erechit;
0208             } else {
0209               negcastortowerarray[2][j] += Erechit;
0210             }
0211             negcastorusedrechits[j].push_back(rechit_p);
0212           }  // end condition over rechit z value
0213         }  // end phi matching condition
0214       }  // end loop over the 16 towers possibilities
0215     }  // end time conditions
0216
0217   }  // end loop over rechits to build castortowerarray[4][16] and castorusedrechits[16]
0218
0219   // make towers of the arrays
0220
0221   double fem, Ehot, depth;
0222   double rhoTower = 88.5;
0223
0224   // loop over the 16 towers possibilities
0225   for (int k = 0; k < 16; k++) {
0226     Ehot = 0;
0227     depth = 0;
0228
0229     // select the positive towers with E > sqrt(Nusedrechits)*Ecut
0230     if (poscastortowerarray[0][k] > sqrt(poscastorusedrechits[k].size()) * towercut_) {
0231       fem = poscastortowerarray[1][k] / poscastortowerarray[0][k];
0232       CastorRecHitRefVector usedRecHits = poscastorusedrechits[k];
0233       ComputeTowerVariable(usedRecHits, Ehot, depth);
0234
0235       LogDebug("CastorTowerProducer") << "tower " << k + 1 << ": fem = " << fem << " ,depth = " << depth
0236                                       << " ,Ehot = " << Ehot << std::endl;
0237
0238       TowerPoint temptowerposition(rhoTower, 5.9, poscastortowerarray[3][k]);
0239       Point towerposition(temptowerposition);
0240
0241       CastorTower newtower(poscastortowerarray[0][k],
0242                            towerposition,
0243                            poscastortowerarray[1][k],
0244                            poscastortowerarray[2][k],
0245                            fem,
0246                            depth,
0247                            Ehot,
0248                            poscastorusedrechits[k]);
0249       OutputTowers->push_back(newtower);
0250     }  // end select the positive towers with E > Ecut
0251
0252     // select the negative towers with E > sqrt(Nusedrechits)*Ecut
0253     if (negcastortowerarray[0][k] > sqrt(negcastorusedrechits[k].size()) * towercut_) {
0254       fem = negcastortowerarray[1][k] / negcastortowerarray[0][k];
0255       CastorRecHitRefVector usedRecHits = negcastorusedrechits[k];
0256       ComputeTowerVariable(usedRecHits, Ehot, depth);
0257
0258       LogDebug("CastorTowerProducer") << "tower " << k + 1 << " energy = " << negcastortowerarray[0][k]
0259                                       << "EM = " << negcastortowerarray[1][k] << "HAD = " << negcastortowerarray[2][k]
0260                                       << "phi = " << negcastortowerarray[3][k] << ": fem = " << fem
0261                                       << " ,depth = " << depth << " ,Ehot = " << Ehot << std::endl;
0262
0263       TowerPoint temptowerposition(rhoTower, -5.9, negcastortowerarray[3][k]);
0264       Point towerposition(temptowerposition);
0265
0266       CastorTower newtower(negcastortowerarray[0][k],
0267                            towerposition,
0268                            negcastortowerarray[1][k],
0269                            negcastortowerarray[2][k],
0270                            fem,
0271                            depth,
0272                            Ehot,
0273                            negcastorusedrechits[k]);
0274       OutputTowers->push_back(newtower);
0275     }  // end select the negative towers with E > Ecut
0276
0277   }  // end loop over the 16 towers possibilities
0278
0279   iEvent.put(std::move(OutputTowers));
0280 }
0281
0282 void CastorTowerProducer::ComputeTowerVariable(const edm::RefVector<edm::SortedCollection<CastorRecHit> >& usedRecHits,
0283                                                double& Ehot,
0284                                                double& depth) {
0285   using namespace reco;
0286
0287   double Etot = 0;
0288
0289   // loop over the cells used in the tower k
0290   for (CastorRecHitRefVector::iterator it = usedRecHits.begin(); it != usedRecHits.end(); it++) {
0291     edm::Ref<CastorRecHitCollection> rechit_p = *it;
0292
0293     double Erechit = rechit_p->energy();
0294     HcalCastorDetId id = rechit_p->id();
0295     int module = id.module();
0296     double zrechit = 0;
0297     if (module < 3)
0298       zrechit = -14390 - 24.75 - 49.5 * (module - 1);
0299     if (module > 2)
0300       zrechit = -14390 - 99 - 49.5 - 99 * (module - 3);
0301
0302     if (Erechit > Ehot)
0303       Ehot = Erechit;
0304     depth += Erechit * zrechit;
0305     Etot += Erechit;
0306   }
0307
0308   depth /= Etot;
0309   Ehot /= Etot;
0310 }
0311
0312 //define this as a plug-in
0313 DEFINE_FWK_MODULE(CastorTowerProducer);