CaloHitMakers/test/testEcalHitMaker.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    testCaloCaloGeometryTools
0004 // Class:      testCaloCaloGeometryTools
0005 //
0006 /**\class testCaloCaloGeometryTools testEcalHitMaker.cc test/testEcalHitMaker/src/testEcalHitMaker.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014
0015 // system include files
0016 #include <iomanip>
0017 #include <memory>
0018
0019 // user include files
0020 #include "FWCore/Framework/interface/stream/EDAnalyzer.h"
0021
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026
0027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0028 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0029 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031
0032 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0033 #include "DataFormats/Math/interface/Vector3D.h"
0034 #include "DataFormats/Math/interface/LorentzVector.h"
0035
0036 #include "FastSimulation/ShowerDevelopment/interface/EMECALShowerParametrization.h"
0037 #include "FastSimulation/ShowerDevelopment/interface/EMShower.h"
0038 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0039 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
0040 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
0041 #include "CommonTools/BaseParticlePropagator/interface/RawParticle.h"
0042 #include "FastSimulation/Event/interface/FSimTrack.h"
0043 #include "FastSimulation/Event/interface/FSimEvent.h"
0044 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0045 #include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
0046
0047 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0048
0049 //
0050 // class decleration
0051 //
0052
0053 typedef math::XYZVector XYZVector;
0054 typedef math::XYZVector XYZPoint;
0055
0056 class testEcalHitMaker : public edm::stream::EDAnalyzer<> {
0057 public:
0058   explicit testEcalHitMaker(const edm::ParameterSet&);
0059   ~testEcalHitMaker();
0060   virtual void beginRun(edm::Run const&, edm::EventSetup const&);
0061
0062   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0063
0064 private:
0065   // ----------member data ---------------------------
0066   void testBorderCrossing();
0067   int pass_;
0068
0069   const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tokCaloTopo_;
0070   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokCaloGeom_;
0071   const edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord> tokPdt_;
0072
0073   CaloGeometryHelper* myGeometry;
0074   GammaFunctionGenerator* aGammaGenerator;
0075   FSimEvent* mySimEvent;
0076 };
0077
0078 //
0079 // constants, enums and typedefs
0080 //
0081
0082 //
0083 // static data member definitions
0084 //
0085
0086 //
0087 // constructors and destructor
0088 //
0089 testEcalHitMaker::testEcalHitMaker(const edm::ParameterSet& iConfig)
0090     : tokCaloTopo_(esConsumes<edm::Transition::BeginRun>()),
0091       tokCaloGeom_(esConsumes<edm::Transition::BeginRun>()),
0092       tokPdt_(esConsumes<edm::Transition::BeginRun>()) {
0093   aGammaGenerator = new GammaFunctionGenerator();
0094
0095   mySimEvent = new FSimEvent(iConfig.getParameter<edm::ParameterSet>("TestParticleFilter"));
0096
0097   myGeometry = new CaloGeometryHelper(iConfig.getParameter<edm::ParameterSet>("Calorimetry"));
0098 }
0099
0100 void testEcalHitMaker::beginRun(edm::Run const& run, edm::EventSetup const& iSetup) {
0101   const edm::ESHandle<CaloTopology>& theCaloTopology = iSetup.getHandle(tokCaloTopo_);
0102   const edm::ESHandle<CaloGeometry>& pG = iSetup.getHandle(tokCaloGeom_);
0103
0104   // Setup the tools
0105   double bField000 = 4.;
0106   myGeometry->setupGeometry(*(pG.product()));
0107   myGeometry->setupTopology(*(theCaloTopology.product()));
0108   myGeometry->initialize(bField000);
0109
0110   // init Particle data table (from Pythia)
0111   const edm::ESHandle<HepPDT::ParticleDataTable>& pdt = iSetup.getHandle(tokPdt_);
0112   mySimEvent->initializePdt(pdt.product());
0113   std::cout << " done with beginRun " << std::endl;
0114 }
0115
0116 testEcalHitMaker::~testEcalHitMaker() {
0117   // do anything here that needs to be done at desctruction time
0118   // (e.g. close files, deallocate resources etc.)
0119 }
0120
0121 // ------------ method called to produce the data  ------------
0122 void testEcalHitMaker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0123   using namespace edm;
0124
0125   RandomEngineAndDistribution random(iEvent.streamID());
0126
0127   math::XYZTLorentzVectorD theMomentum(10., 0., 5., sqrt(125));
0128
0129   // no need actually define it at the ECAL entrance: the fill of FSimEvent will do the
0130   // propagation
0131   math::XYZVectorD thePositionatEcalEntrance(129., 0., 60);
0132
0133   std::vector<SimTrack> mySimTracks;
0134   SimTrack myTrack(11, theMomentum, 0, -1, thePositionatEcalEntrance, theMomentum);
0135   mySimTracks.push_back(myTrack);
0136   std::vector<SimVertex> mySimVertices;
0137   mySimVertices.push_back(SimVertex(thePositionatEcalEntrance, 0.));
0138
0139   mySimEvent->fill(mySimTracks, mySimVertices);
0140
0141   //   RawParticle myPart(11,theMomentum);
0142   //   myPart.setVertex(thePositionatEcalEntrance.X(),
0143   //            thePositionatEcalEntrance.Y(),
0144   //            thePositionatEcalEntrance.Z(),
0145   //            0.);
0146   //   FSimTrack mySimTrack(&myPart,-1,-1,11,mySimEvent);
0147   //   mySimTrack.setTkPosition(thePositionatEcalEntrance);
0148   //   mySimTrack.setTkMomentum(theMomentum);
0149   //   mySimTrack.setEcal(myPart,1);
0150   //   // put dummy quantities
0151   //   mySimTrack.setLayer1(RawParticle(),0);
0152   //   mySimTrack.setLayer2(RawParticle(),0);
0153
0154   FSimTrack& mySimTrack(mySimEvent->track(0));
0155   std::cout << mySimTrack << std::endl;
0156   std::cout << " done " << std::endl;
0157   RawParticle myPart = mySimTrack.ecalEntrance();
0158   std::vector<const RawParticle*> thePart;
0159   thePart.push_back(&myPart);
0160   // no preshower
0161   XYZPoint ecalentrance = myPart.vertex().Vect();
0162
0163   std::cout << " on ECAL/HCAL etc " << mySimTrack.onEcal() << " ";
0164   std::cout << mySimTrack.onHcal() << " ";
0165   std::cout << mySimTrack.onLayer1() << " ";
0166   std::cout << mySimTrack.onLayer2() << " " << std::endl;
0167
0168   // ask me if you want details - makes the simulation faster
0169   std::vector<double> coreParams;
0170   coreParams.push_back(100);
0171   coreParams.push_back(0.1);
0172   std::vector<double> tailParams;
0173   tailParams.push_back(1);
0174   tailParams.push_back(0.1);
0175   tailParams.push_back(100);
0176   tailParams.push_back(1);
0177
0178   // define the calorimeter properties
0179   EMECALShowerParametrization showerparam(myGeometry->ecalProperties(mySimTrack.onEcal()),
0180                                           myGeometry->hcalProperties(mySimTrack.onHcal()),
0181                                           myGeometry->layer1Properties(mySimTrack.onLayer1()),
0182                                           myGeometry->layer2Properties(mySimTrack.onLayer2()),
0183                                           coreParams,
0184                                           tailParams);
0185
0186   //define the shower parameters
0187   EMShower theShower(&random, aGammaGenerator, &showerparam, &thePart);
0188
0189   // you might want to replace this with something elese
0190   DetId pivot(myGeometry->getClosestCell(ecalentrance, true, mySimTrack.onEcal() == 1));
0191
0192   // define the 7x7 grid
0193   EcalHitMaker myGrid(myGeometry, ecalentrance, pivot, mySimTrack.onEcal(), 7, 0, &random);
0194   myGrid.setCrackPadSurvivalProbability(0.9);  // current parameters  in the Fast Sim
0195   myGrid.setRadiusFactor(1.096);               // current parameters
0196
0197   // define the track parameters
0198   myGrid.setTrackParameters(myPart.Vect().Unit(), 0., mySimTrack);
0199
0200   HcalHitMaker myHcalHitMaker(myGrid, (unsigned)0);
0201
0202   theShower.setGrid(&myGrid);
0203   theShower.setHcal(&myHcalHitMaker);
0204   theShower.compute();
0205
0206   // print the result
0207   std::map<CaloHitID, float>::const_iterator mapitr;
0208   std::map<CaloHitID, float>::const_iterator endmapitr = myGrid.getHits().end();
0209   for (mapitr = myGrid.getHits().begin(); mapitr != endmapitr; ++mapitr) {
0210     if (mapitr->second != 0)
0211       std::cout << "DetId " << EBDetId(mapitr->first.unitID()) << " " << std::setw(8) << std::setprecision(4)
0212                 << mapitr->second << std::endl;
0213   }
0214 }
0215
0216 //define this as a plug-in
0217 //DEFINE_SEAL_MODULE();
0218 //DEFINE_ANOTHER_FWK_MODULE(testEcalHitMaker);
0219 DEFINE_FWK_MODULE(testEcalHitMaker);