PreshowerHitMaker.cc

Go to the documentation of this file.

#include "FastSimulation/CaloHitMakers/interface/PreshowerHitMaker.h"
#
#include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
#include "FastSimulation/Utilities/interface/LandauFluctuationGenerator.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"

#include "Math/GenVector/Plane3D.h"

#include <cmath>

typedef ROOT::Math::Plane3D Plane3D;
typedef ROOT::Math::Transform3DPJ::Point Point;

// LandauFluctuationGenerator PreshowerHitMaker::theGenerator=LandauFluctuationGenerator();

PreshowerHitMaker::PreshowerHitMaker(CaloGeometryHelper* calo,
                                     const XYZPoint& layer1entrance,
                                     const XYZVector& layer1dir,
                                     const XYZPoint& layer2entrance,
                                     const XYZVector& layer2dir,
                                     const LandauFluctuationGenerator* aGenerator,
                                     const RandomEngineAndDistribution* engine)
    : CaloHitMaker(calo, DetId::Ecal, EcalPreshower, 2),
      psLayer1Entrance_(layer1entrance),
      psLayer1Dir_(layer1dir),
      psLayer2Entrance_(layer2entrance),
      psLayer2Dir_(layer2dir),
      totalLayer1_(0.),
      totalLayer2_(0.),
      theGenerator(aGenerator),
      random(engine) {
  double dummyt;
  anglecorrection1_ = 0.;
  anglecorrection2_ = 0.;
  // Check if the entrance points are really on the wafers
  // Layer 1
  layer1valid_ = (layer1entrance.Mag2() > 0.);
  if (layer1valid_) {
    int z = (psLayer1Entrance_.z() > 0) ? 1 : -1;
    Plane3D plan1(0., 0., 1., -z * myCalorimeter->preshowerZPosition(1));

    psLayer1Entrance_ = intersect(plan1, layer1entrance, layer1entrance + layer1dir, dummyt, false);

    XYZVector zaxis(0, 0, 1);
    XYZVector planeVec1 = (zaxis.Cross(layer1dir)).Unit();
    locToGlobal1_ = Transform3D(Point(0, 0, 0),
                                Point(0, 0, 1),
                                Point(1, 0, 0),
                                (Point)psLayer1Entrance_,
                                (Point)(psLayer1Entrance_ + layer1dir),
                                (Point)(psLayer1Entrance_ + planeVec1));

    anglecorrection1_ = fabs(zaxis.Dot(layer1dir));
    if (anglecorrection1_ != 0.)
      anglecorrection1_ = 1. / anglecorrection1_;
    //      std::cout << " Layer 1 entrance " << psLayer1Entrance_ << std::endl;
    //      std::cout << " Layer 1 corr " << anglecorrection1_ << std::endl;
  }

  // Layer 2
  layer2valid_ = (layer2entrance.Mag2() > 0.);
  if (layer2valid_) {
    int z = (psLayer2Entrance_.z() > 0) ? 1 : -1;
    Plane3D plan2(0., 0., 1., -z * myCalorimeter->preshowerZPosition(2));

    psLayer2Entrance_ = intersect(plan2, layer2entrance, layer2entrance + layer2dir, dummyt, false);

    XYZVector zaxis(0, 0, 1);
    XYZVector planeVec2 = (zaxis.Cross(layer2dir)).Unit();
    locToGlobal2_ = Transform3D(Point(0, 0, 0),
                                Point(0, 0, 1),
                                Point(1, 0, 0),
                                (Point)psLayer2Entrance_,
                                (Point)(psLayer2Entrance_ + layer2dir),
                                (Point)(psLayer2Entrance_ + planeVec2));

    anglecorrection2_ = fabs(zaxis.Dot(layer2dir));
    if (anglecorrection2_ != 0.)
      anglecorrection2_ = 1. / anglecorrection2_;
    //      std::cout << " Layer 2 entrance " << psLayer2Entrance_ << std::endl;
    //      std::cout << " Layer 2 corr " << anglecorrection2_ << std::endl;
  }
  //  theGenerator=LandauFluctuationGenerator();
}

bool PreshowerHitMaker::addHit(double r, double phi, unsigned layer) {
  if ((layer == 1 && !layer1valid_) || ((layer == 2 && !layer2valid_)))
    return false;

  r *= moliereRadius;
  XYZPoint point(r * std::cos(phi), r * std::sin(phi), 0.);
  point = (layer == 1) ? locToGlobal1_((Point)point) : locToGlobal2_((Point)point);
  //  std::cout << "  Point " << point  << std::endl;
  int z = (point.z() > 0) ? 1 : -1;
  point = XYZPoint(point.x(), point.y(), z * myCalorimeter->preshowerZPosition(layer));
  //  std::cout << "r " << r << "  Point after " << point  << std::endl;
  //  std::cout << " Layer " << layer << " " << point << std::endl;
  DetId strip = myCalorimeter->getEcalPreshowerGeometry()->getClosestCellInPlane(
      GlobalPoint(point.x(), point.y(), point.z()), layer);

  float meanspot = (layer == 1) ? mip1_ : mip2_;
  float spote = meanspot + 0.000021 * theGenerator->landau(random);
  spote *= ((layer == 1) ? anglecorrection1_ : anglecorrection2_);

  if (!strip.null()) {
    //calculate time of flight
    double tof =
        (myCalorimeter->getEcalPreshowerGeometry()->getGeometry(strip)->getPosition().mag()) / 29.98;  //speed of light
    CaloHitID current_id(strip.rawId(), tof, 0);                                                       //no track yet
    std::map<CaloHitID, float>::iterator cellitr;
    cellitr = hitMap_.find(current_id);
    if (cellitr == hitMap_.end()) {
      hitMap_.insert(std::pair<CaloHitID, float>(current_id, spote));
    } else {
      cellitr->second += spote;
    }
    //      std::cout << " found " << stripNumber << " " << spote <<std::endl;
    if (layer == 1) {
      totalLayer1_ += spote;
    } else if (layer == 2) {
      totalLayer2_ += spote;
    }
    return true;
  }
  //  std::cout << "  Could not find a cell " << point << std::endl;
  return false;
}