CaloHitMakers/src/PreshowerHitMaker.cc

0001 #include "FastSimulation/CaloHitMakers/interface/PreshowerHitMaker.h"
0002 #
0003 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0004 #include "FastSimulation/Utilities/interface/LandauFluctuationGenerator.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0007 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0008
0009 #include "Math/GenVector/Plane3D.h"
0010
0011 #include <cmath>
0012
0013 typedef ROOT::Math::Plane3D Plane3D;
0014 typedef ROOT::Math::Transform3DPJ::Point Point;
0015
0016 // LandauFluctuationGenerator PreshowerHitMaker::theGenerator=LandauFluctuationGenerator();
0017
0018 PreshowerHitMaker::PreshowerHitMaker(CaloGeometryHelper* calo,
0019                                      const XYZPoint& layer1entrance,
0020                                      const XYZVector& layer1dir,
0021                                      const XYZPoint& layer2entrance,
0022                                      const XYZVector& layer2dir,
0023                                      const LandauFluctuationGenerator* aGenerator,
0024                                      const RandomEngineAndDistribution* engine)
0025     : CaloHitMaker(calo, DetId::Ecal, EcalPreshower, 2),
0026       psLayer1Entrance_(layer1entrance),
0027       psLayer1Dir_(layer1dir),
0028       psLayer2Entrance_(layer2entrance),
0029       psLayer2Dir_(layer2dir),
0030       totalLayer1_(0.),
0031       totalLayer2_(0.),
0032       theGenerator(aGenerator),
0033       random(engine) {
0034   double dummyt;
0035   anglecorrection1_ = 0.;
0036   anglecorrection2_ = 0.;
0037   // Check if the entrance points are really on the wafers
0038   // Layer 1
0039   layer1valid_ = (layer1entrance.Mag2() > 0.);
0040   if (layer1valid_) {
0041     int z = (psLayer1Entrance_.z() > 0) ? 1 : -1;
0042     Plane3D plan1(0., 0., 1., -z * myCalorimeter->preshowerZPosition(1));
0043
0044     psLayer1Entrance_ = intersect(plan1, layer1entrance, layer1entrance + layer1dir, dummyt, false);
0045
0046     XYZVector zaxis(0, 0, 1);
0047     XYZVector planeVec1 = (zaxis.Cross(layer1dir)).Unit();
0048     locToGlobal1_ = Transform3D(Point(0, 0, 0),
0049                                 Point(0, 0, 1),
0050                                 Point(1, 0, 0),
0051                                 (Point)psLayer1Entrance_,
0052                                 (Point)(psLayer1Entrance_ + layer1dir),
0053                                 (Point)(psLayer1Entrance_ + planeVec1));
0054
0055     anglecorrection1_ = fabs(zaxis.Dot(layer1dir));
0056     if (anglecorrection1_ != 0.)
0057       anglecorrection1_ = 1. / anglecorrection1_;
0058     //      std::cout << " Layer 1 entrance " << psLayer1Entrance_ << std::endl;
0059     //      std::cout << " Layer 1 corr " << anglecorrection1_ << std::endl;
0060   }
0061
0062   // Layer 2
0063   layer2valid_ = (layer2entrance.Mag2() > 0.);
0064   if (layer2valid_) {
0065     int z = (psLayer2Entrance_.z() > 0) ? 1 : -1;
0066     Plane3D plan2(0., 0., 1., -z * myCalorimeter->preshowerZPosition(2));
0067
0068     psLayer2Entrance_ = intersect(plan2, layer2entrance, layer2entrance + layer2dir, dummyt, false);
0069
0070     XYZVector zaxis(0, 0, 1);
0071     XYZVector planeVec2 = (zaxis.Cross(layer2dir)).Unit();
0072     locToGlobal2_ = Transform3D(Point(0, 0, 0),
0073                                 Point(0, 0, 1),
0074                                 Point(1, 0, 0),
0075                                 (Point)psLayer2Entrance_,
0076                                 (Point)(psLayer2Entrance_ + layer2dir),
0077                                 (Point)(psLayer2Entrance_ + planeVec2));
0078
0079     anglecorrection2_ = fabs(zaxis.Dot(layer2dir));
0080     if (anglecorrection2_ != 0.)
0081       anglecorrection2_ = 1. / anglecorrection2_;
0082     //      std::cout << " Layer 2 entrance " << psLayer2Entrance_ << std::endl;
0083     //      std::cout << " Layer 2 corr " << anglecorrection2_ << std::endl;
0084   }
0085   //  theGenerator=LandauFluctuationGenerator();
0086 }
0087
0088 bool PreshowerHitMaker::addHit(double r, double phi, unsigned layer) {
0089   if ((layer == 1 && !layer1valid_) || ((layer == 2 && !layer2valid_)))
0090     return false;
0091
0092   r *= moliereRadius;
0093   XYZPoint point(r * std::cos(phi), r * std::sin(phi), 0.);
0094   point = (layer == 1) ? locToGlobal1_((Point)point) : locToGlobal2_((Point)point);
0095   //  std::cout << "  Point " << point  << std::endl;
0096   int z = (point.z() > 0) ? 1 : -1;
0097   point = XYZPoint(point.x(), point.y(), z * myCalorimeter->preshowerZPosition(layer));
0098   //  std::cout << "r " << r << "  Point after " << point  << std::endl;
0099   //  std::cout << " Layer " << layer << " " << point << std::endl;
0100   DetId strip = myCalorimeter->getEcalPreshowerGeometry()->getClosestCellInPlane(
0101       GlobalPoint(point.x(), point.y(), point.z()), layer);
0102
0103   float meanspot = (layer == 1) ? mip1_ : mip2_;
0104   float spote = meanspot + 0.000021 * theGenerator->landau(random);
0105   spote *= ((layer == 1) ? anglecorrection1_ : anglecorrection2_);
0106
0107   if (!strip.null()) {
0108     //calculate time of flight
0109     double tof =
0110         (myCalorimeter->getEcalPreshowerGeometry()->getGeometry(strip)->getPosition().mag()) / 29.98;  //speed of light
0111     CaloHitID current_id(strip.rawId(), tof, 0);                                                       //no track yet
0112     std::map<CaloHitID, float>::iterator cellitr;
0113     cellitr = hitMap_.find(current_id);
0114     if (cellitr == hitMap_.end()) {
0115       hitMap_.insert(std::pair<CaloHitID, float>(current_id, spote));
0116     } else {
0117       cellitr->second += spote;
0118     }
0119     //      std::cout << " found " << stripNumber << " " << spote <<std::endl;
0120     if (layer == 1) {
0121       totalLayer1_ += spote;
0122     } else if (layer == 2) {
0123       totalLayer2_ += spote;
0124     }
0125     return true;
0126   }
0127   //  std::cout << "  Could not find a cell " << point << std::endl;
0128   return false;
0129 }