Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    Castor
00004 // Class:      CastorTowerProducer
00005 // 
00006 
00013 //
00014 // Original Author:  Hans Van Haevermaet, Benoit Roland
00015 //         Created:  Wed Jul  9 14:00:40 CEST 2008
00016 // $Id: CastorTowerProducer.cc,v 1.12 2013/01/07 15:32:45 hvanhaev Exp $
00017 //
00018 //
00019 
00020 
00021 // system include 
00022 #include <memory>
00023 #include <vector>
00024 #include <iostream>
00025 #include <TMath.h>
00026 #include <TRandom3.h>
00027 
00028 // user include 
00029 #include "FWCore/Framework/interface/Frameworkfwd.h"
00030 #include "FWCore/Framework/interface/EDProducer.h"
00031 #include "FWCore/Framework/interface/Event.h"
00032 #include "FWCore/Framework/interface/MakerMacros.h"
00033 #include "FWCore/Framework/interface/ESHandle.h"
00034 #include "FWCore/Framework/interface/EventSetup.h"
00035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00036 #include "DataFormats/Math/interface/Point3D.h"
00037 
00038 // Castor Object include
00039 #include "DataFormats/HcalRecHit/interface/CastorRecHit.h"
00040 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
00041 #include "DataFormats/CastorReco/interface/CastorTower.h"
00042 
00043 // Channel quality
00044 #include "CondFormats/CastorObjects/interface/CastorChannelQuality.h"
00045 #include "CondFormats/CastorObjects/interface/CastorChannelStatus.h"
00046 #include "CondFormats/DataRecord/interface/CastorChannelQualityRcd.h"
00047 
00048 //
00049 // class declaration
00050 //
00051 
00052 class CastorTowerProducer : public edm::EDProducer {
00053    public:
00054       explicit CastorTowerProducer(const edm::ParameterSet&);
00055       ~CastorTowerProducer();
00056 
00057    private:
00058       virtual void beginJob() ;
00059       virtual void produce(edm::Event&, const edm::EventSetup&);
00060       virtual void endJob() ;
00061       virtual void ComputeTowerVariable(const edm::RefVector<edm::SortedCollection<CastorRecHit> >& usedRecHits, double&  Ehot, double& depth);
00062       
00063       // member data
00064       typedef math::XYZPointD Point;
00065       typedef ROOT::Math::RhoEtaPhiPoint TowerPoint;
00066       typedef ROOT::Math::RhoZPhiPoint CellPoint;
00067       typedef edm::SortedCollection<CastorRecHit> CastorRecHitCollection; 
00068       typedef std::vector<reco::CastorTower> CastorTowerCollection;
00069       typedef edm::RefVector<CastorRecHitCollection> CastorRecHitRefVector;
00070       std::string input_;
00071       double towercut_;
00072       double mintime_;
00073       double maxtime_;
00074 };
00075 
00076 //
00077 // constants, enums and typedefs
00078 //
00079 
00080 const double MYR2D = 180/M_PI;
00081 
00082 //
00083 // static data member definitions
00084 //
00085 
00086 //
00087 // constructor and destructor
00088 //
00089 
00090 CastorTowerProducer::CastorTowerProducer(const edm::ParameterSet& iConfig) :
00091   input_(iConfig.getParameter<std::string>("inputprocess")),
00092   towercut_(iConfig.getParameter<double>("towercut")),
00093   mintime_(iConfig.getParameter<double>("mintime")),
00094   maxtime_(iConfig.getParameter<double>("maxtime"))
00095 {
00096   //register your products
00097   produces<CastorTowerCollection>();
00098   //now do what ever other initialization is needed
00099 }
00100 
00101 
00102 CastorTowerProducer::~CastorTowerProducer()
00103 {
00104    // do anything here that needs to be done at desctruction time
00105    // (e.g. close files, deallocate resources etc.)
00106 }
00107 
00108 
00109 //
00110 // member functions
00111 //
00112 
00113 // ------------ method called to produce the data  ------------
00114 void CastorTowerProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
00115 
00116   using namespace edm;
00117   using namespace reco;
00118   using namespace TMath;
00119   
00120   // Produce CastorTowers from CastorCells
00121   
00122   edm::Handle<CastorRecHitCollection> InputRecHits;
00123   iEvent.getByLabel(input_,InputRecHits);
00124 
00125   std::auto_ptr<CastorTowerCollection> OutputTowers (new CastorTowerCollection);
00126    
00127   // get and check input size
00128   int nRecHits = InputRecHits->size();
00129 
00130   LogDebug("CastorTowerProducer")
00131     <<"2. entering CastorTowerProducer"<<std::endl;
00132 
00133   if (nRecHits==0)
00134     LogDebug("CastorTowerProducer") <<"Warning: You are trying to run the Tower algorithm with 0 input rechits.";
00135   
00136   // declare castor array
00137   // (0,x): Energies - (1,x): emEnergies - (2,x): hadEnergies - (3,x): phi position
00138   
00139   double poscastortowerarray[4][16]; 
00140   double negcastortowerarray[4][16];
00141 
00142   CastorRecHitRefVector poscastorusedrechits[16];
00143   CastorRecHitRefVector negcastorusedrechits[16];
00144 
00145   // set phi values and everything else to zero
00146   for (int j = 0; j < 16; j++) {
00147     poscastortowerarray[3][j] = -2.94524 + j*0.3927;
00148     poscastortowerarray[0][j] = 0.;
00149     poscastortowerarray[1][j] = 0.;
00150     poscastortowerarray[2][j] = 0.;
00151 
00152     negcastortowerarray[3][j] = -2.94524 + j*0.3927;
00153     negcastortowerarray[0][j] = 0.;
00154     negcastortowerarray[1][j] = 0.;
00155     negcastortowerarray[2][j] = 0.;
00156   }
00157   
00158   // retrieve the channel quality lists from database
00159   edm::ESHandle<CastorChannelQuality> p;
00160   iSetup.get<CastorChannelQualityRcd>().get(p);
00161   std::vector<DetId> channels = p->getAllChannels();
00162 
00163   // loop over rechits to build castortowerarray[4][16] and castorusedrechits[16] 
00164   for (unsigned int i = 0; i < InputRecHits->size(); i++) {
00165     
00166     edm::Ref<CastorRecHitCollection> rechit_p = edm::Ref<CastorRecHitCollection>(InputRecHits, i);
00167     
00168     HcalCastorDetId id = rechit_p->id();
00169     DetId genericID=(DetId)id;
00170     
00171     // first check if the rechit is in the BAD channel list
00172     bool bad = false;
00173     for (std::vector<DetId>::iterator channel = channels.begin();channel !=  channels.end();channel++) {        
00174         if (channel->rawId() == genericID.rawId()) {
00175                 // if the rechit is found in the list, set it bad
00176           bad = true; break;
00177         }
00178     }
00179     // if bad, continue the loop to the next rechit
00180     if (bad) continue;
00181     
00182     double Erechit = rechit_p->energy();
00183     int module = id.module();
00184     int sector = id.sector();
00185     double zrechit = 0;
00186     if (module < 3) zrechit = -14390 - 24.75 - 49.5*(module-1);
00187     if (module > 2) zrechit = -14390 - 99 - 49.5 - 99*(module-3); 
00188     double phirechit = -100;
00189     if (sector < 9) phirechit = 0.19635 + (sector-1)*0.3927;
00190     if (sector > 8) phirechit = -2.94524 + (sector - 9)*0.3927;
00191 
00192     // add time conditions for the rechit
00193     if (rechit_p->time() > mintime_ && rechit_p->time() < maxtime_) {
00194 
00195         // loop over the 16 towers possibilities
00196         for ( int j=0;j<16;j++) {
00197       
00198                 // phi matching condition
00199                 if (TMath::Abs(phirechit - poscastortowerarray[3][j]) < 0.0001) {
00200 
00201                         // condition over rechit z value
00202                         if (zrechit > 0.) {
00203                                 poscastortowerarray[0][j]+=Erechit;
00204                                 if (module < 3) {poscastortowerarray[1][j]+=Erechit;} else {poscastortowerarray[2][j]+=Erechit;}
00205                                 poscastorusedrechits[j].push_back(rechit_p);
00206                         } else {
00207                                 negcastortowerarray[0][j]+=Erechit;
00208                                 if (module < 3) {negcastortowerarray[1][j]+=Erechit;} else {negcastortowerarray[2][j]+=Erechit;}
00209                                 negcastorusedrechits[j].push_back(rechit_p);
00210                         } // end condition over rechit z value
00211                 } // end phi matching condition
00212         } // end loop over the 16 towers possibilities
00213     } // end time conditions
00214     
00215   } // end loop over rechits to build castortowerarray[4][16] and castorusedrechits[16]
00216   
00217   // make towers of the arrays
00218 
00219   double fem, Ehot, depth;
00220   double rhoTower = 88.5;
00221 
00222   // loop over the 16 towers possibilities
00223   for (int k=0;k<16;k++) {
00224     
00225     fem = 0;
00226     Ehot = 0;
00227     depth = 0;
00228 
00229     // select the positive towers with E > sqrt(Nusedrechits)*Ecut
00230     if (poscastortowerarray[0][k] > sqrt(poscastorusedrechits[k].size())*towercut_) {
00231       
00232       fem = poscastortowerarray[1][k]/poscastortowerarray[0][k];
00233       CastorRecHitRefVector usedRecHits = poscastorusedrechits[k];
00234       ComputeTowerVariable(usedRecHits,Ehot,depth);
00235 
00236       LogDebug("CastorTowerProducer")
00237         <<"tower "<<k+1<<": fem = "<<fem<<" ,depth = "<<depth<<" ,Ehot = "<<Ehot<<std::endl;
00238 
00239       TowerPoint temptowerposition(rhoTower,5.9,poscastortowerarray[3][k]);
00240       Point towerposition(temptowerposition);
00241 
00242       CastorTower newtower(poscastortowerarray[0][k],towerposition,poscastortowerarray[1][k],poscastortowerarray[2][k],fem,depth,Ehot,
00243                            poscastorusedrechits[k]);
00244       OutputTowers->push_back(newtower);
00245     } // end select the positive towers with E > Ecut
00246     
00247     // select the negative towers with E > sqrt(Nusedrechits)*Ecut
00248     if (negcastortowerarray[0][k] > sqrt(negcastorusedrechits[k].size())*towercut_) {
00249       
00250       fem = negcastortowerarray[1][k]/negcastortowerarray[0][k];
00251       CastorRecHitRefVector usedRecHits = negcastorusedrechits[k];
00252       ComputeTowerVariable(usedRecHits,Ehot,depth);
00253 
00254       LogDebug("CastorTowerProducer")
00255          <<"tower "<<k+1 << " energy = " << negcastortowerarray[0][k] << "EM = " << negcastortowerarray[1][k] << "HAD = " << negcastortowerarray[2][k] << "phi = " << negcastortowerarray[3][k] << ": fem = "<<fem<<" ,depth = "<<depth<<" ,Ehot = "<<Ehot<<std::endl;
00256 
00257       TowerPoint temptowerposition(rhoTower,-5.9,negcastortowerarray[3][k]);
00258       Point towerposition(temptowerposition);
00259 
00260       CastorTower newtower(negcastortowerarray[0][k],towerposition,negcastortowerarray[1][k],negcastortowerarray[2][k],fem,depth,Ehot,
00261                            negcastorusedrechits[k]);
00262       OutputTowers->push_back(newtower);
00263     } // end select the negative towers with E > Ecut
00264     
00265   } // end loop over the 16 towers possibilities
00266   
00267   iEvent.put(OutputTowers);
00268 } 
00269 
00270 
00271 // ------------ method called once each job just before starting event loop  ------------
00272 void CastorTowerProducer::beginJob() {
00273   LogDebug("CastorTowerProducer")
00274     <<"Starting CastorTowerProducer";
00275 }
00276 
00277 // ------------ method called once each job just after ending the event loop  ------------
00278 void CastorTowerProducer::endJob() {
00279   LogDebug("CastorTowerProducer")
00280     <<"Ending CastorTowerProducer";
00281 }
00282 
00283 void CastorTowerProducer::ComputeTowerVariable(const edm::RefVector<edm::SortedCollection<CastorRecHit> >& usedRecHits, double&  Ehot, double& depth) {
00284 
00285   using namespace reco;
00286 
00287   double Etot = 0;
00288 
00289   // loop over the cells used in the tower k
00290   for (CastorRecHitRefVector::iterator it = usedRecHits.begin(); it != usedRecHits.end(); it++) {
00291     edm::Ref<CastorRecHitCollection> rechit_p = *it;
00292 
00293     double Erechit = rechit_p->energy();
00294     HcalCastorDetId id = rechit_p->id();
00295     int module = id.module();
00296     double zrechit = 0;
00297     if (module < 3) zrechit = -14390 - 24.75 - 49.5*(module-1);
00298     if (module > 2) zrechit = -14390 - 99 - 49.5 - 99*(module-3); 
00299 
00300     if(Erechit > Ehot) Ehot = Erechit;
00301     depth+=Erechit*zrechit;
00302     Etot+=Erechit;
00303   }
00304 
00305   depth/=Etot;
00306   Ehot/=Etot;
00307 }
00308 
00309 //define this as a plug-in
00310 DEFINE_FWK_MODULE(CastorTowerProducer);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoLocalCalo/Castor/src/CastorTowerProducer.cc
