#include <TBPositionCalc.h>

Public Member Functions
CLHEP::Hep3Vector	CalculateCMSPos (std::vector< EBDetId > passedDetIds, int myCrystal, EcalRecHitCollection const *passedRecHitsMap)
CLHEP::Hep3Vector	CalculateTBPos (std::vector< EBDetId > passedDetIds, int myCrystal, EcalRecHitCollection const *passedRecHitsMap)
void	computeRotation (int myCrystal, CLHEP::HepRotation &CMStoTB)
	TBPositionCalc ()
	TBPositionCalc (std::map< std::string, double > providedParameters, const std::string &mapFile, const CaloSubdetectorGeometry *passedGeometry)
	~TBPositionCalc ()
Private Attributes
bool	param_LogWeighted_
Double32_t	param_T0_
Double32_t	param_W0_
Double32_t	param_X0_
const CaloSubdetectorGeometry *	theGeometry_
EcalTBCrystalMap *	theTestMap

Detailed Description

Definition at line 26 of file TBPositionCalc.h.

Constructor & Destructor Documentation

TBPositionCalc::TBPositionCalc	(	std::map< std::string, double >	providedParameters,
		const std::string &	mapFile,
		const CaloSubdetectorGeometry *	passedGeometry
	)

Definition at line 5 of file TBPositionCalc.cc.

References NULL.

{
  // barrel geometry initialization
  if(passedGeometry == NULL)
    throw(std::runtime_error("\n\n TBPositionCalc: wrong initialization.\n\n"));
  theGeometry_ = passedGeometry;
  
  // parameters initialization
  param_LogWeighted_ = providedParameters.find("LogWeighted")->second;
  param_X0_ = providedParameters.find("X0")->second;
  param_T0_ = providedParameters.find("T0")->second; 
  param_W0_ = providedParameters.find("W0")->second;

  theTestMap = new EcalTBCrystalMap(fullMapName);
}

TBPositionCalc::TBPositionCalc ( ) [inline]

Definition at line 32 of file TBPositionCalc.h.

{ };

TBPositionCalc::~TBPositionCalc ( )

Definition at line 21 of file TBPositionCalc.cc.

{
  if (theTestMap) delete theTestMap;
}

Member Function Documentation

CLHEP::Hep3Vector TBPositionCalc::CalculateCMSPos	(	std::vector< EBDetId >	passedDetIds,
		int	myCrystal,
		EcalRecHitCollection const *	passedRecHitsMap
	)

Definition at line 58 of file TBPositionCalc.cc.

References funct::cos(), edm::SortedCollection< T, SORT >::find(), funct::log(), M_PI, max(), PV3DBase< T, PVType, FrameType >::phi(), CosmicsPD_Skims::radius, funct::sin(), EBDetId::SMCRYSTALMODE, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::theta(), CommonMethods::weight(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                                                                                            {
  
  // Calculate the total energy
  double thisEne = 0;
  double eTot = 0;
  EBDetId myId;
  std::vector<EBDetId>::iterator myIt;
  for (myIt = passedDetIds.begin(); myIt !=  passedDetIds.end(); myIt++) {
    myId = (*myIt);
    EcalRecHitCollection::const_iterator itt = passedRecHitsMap->find(myId);
    thisEne = itt->energy();
    eTot += thisEne;
  }

  // Calculate shower depth
  float depth = 0.;
  if(eTot<=0.) {
    edm::LogError("NegativeClusterEnergy") << "cluster with negative energy: " << eTot << ", setting depth to 0.";
  } else {
    depth = param_X0_ * (param_T0_ + log(eTot));
  }

  // Get position of the central crystal from shower depth 
  EBDetId maxId_ = EBDetId(1, myCrystal, EBDetId::SMCRYSTALMODE);
  const CaloCellGeometry* center_cell = theGeometry_ -> getGeometry(maxId_);
  GlobalPoint center_pos = 
    (dynamic_cast<const TruncatedPyramid*>(center_cell))->getPosition(depth);

  // Loop over the hits collection
  double weight        = 0;
  double total_weight  = 0;
  double cluster_theta = 0;
  double cluster_phi   = 0;
  std::vector<EBDetId>::iterator myIt2;
  for (myIt2 = passedDetIds.begin(); myIt2 != passedDetIds.end(); myIt2++) {

    // getting weights
    myId = (*myIt2);
    EcalRecHitCollection::const_iterator itj = passedRecHitsMap->find(myId);
    double ener = itj->energy();

    if (param_LogWeighted_) {
      if(eTot<=0.) { weight = 0.; } 
      else { weight = std::max(0., param_W0_ + log( fabs(ener)/eTot) ); }
    } else {
      weight = ener/eTot;
    }
    total_weight += weight;

    // weighted position of this detId
    const CaloCellGeometry* jth_cell = theGeometry_->getGeometry(myId);
    GlobalPoint jth_pos = dynamic_cast<const TruncatedPyramid*>(jth_cell)->getPosition(depth);    
    cluster_theta += weight*jth_pos.theta();
    cluster_phi   += weight*jth_pos.phi();
  }
  
  // normalizing
  cluster_theta /= total_weight;
  cluster_phi /= total_weight;
  if (cluster_phi > M_PI) { cluster_phi -= 2.*M_PI; }
  if (cluster_phi < -M_PI){ cluster_phi += 2.*M_PI; }

  // position in the cms frame
  double radius = sqrt(center_pos.x()*center_pos.x() + center_pos.y()*center_pos.y() + center_pos.z()*center_pos.z());
  double xpos = radius*cos(cluster_phi)*sin(cluster_theta);
  double ypos = radius*sin(cluster_phi)*sin(cluster_theta); 
  double zpos = radius*cos(cluster_theta);

  return CLHEP::Hep3Vector(xpos, ypos, zpos);
}

CLHEP::Hep3Vector TBPositionCalc::CalculateTBPos	(	std::vector< EBDetId >	passedDetIds,
		int	myCrystal,
		EcalRecHitCollection const *	passedRecHitsMap
	)

Definition at line 26 of file TBPositionCalc.cc.

References edm::SortedCollection< T, SORT >::end(), edm::SortedCollection< T, SORT >::find(), and NULL.

                                                                                                                                           {
  
  // throw an error if the cluster was not initialized properly  
  if(passedRecHitsMap == NULL)
    throw(std::runtime_error("\n\n TBPositionCalc::CalculateTBPos called uninitialized.\n\n"));
  
  // check DetIds are nonzero
  std::vector<EBDetId> validDetIds;
  std::vector<EBDetId>::iterator iter;
  for (iter = passedDetIds.begin(); iter != passedDetIds.end(); iter++) {
    if (((*iter) != DetId(0)) 
        && (passedRecHitsMap->find(*iter) != passedRecHitsMap->end()))
      validDetIds.push_back(*iter);
  }
  passedDetIds.clear();
  passedDetIds = validDetIds;
  
  // computing the position in the cms frame
  CLHEP::Hep3Vector cmsPos = CalculateCMSPos(passedDetIds, myCrystal, passedRecHitsMap);

  // computing the rotation matrix (from CMS to TB)
  CLHEP::HepRotation *CMStoTB = new CLHEP::HepRotation();
  computeRotation(myCrystal, (*CMStoTB));

  // moving to testbeam frame
  CLHEP::Hep3Vector tbPos = (*CMStoTB)*cmsPos;
  delete CMStoTB;

  return tbPos;
}

void TBPositionCalc::computeRotation	(	int	myCrystal,
		CLHEP::HepRotation &	CMStoTB
	)

Definition at line 130 of file TBPositionCalc.cc.

References funct::exp().

                                                                            {
  
  // taking eta/phi of the crystal

  double myEta   = 999.;
  double myPhi   = 999.;
  double myTheta = 999.;
  theTestMap->findCrystalAngles(MyCrystal, myEta, myPhi);
  myTheta = 2.0*atan(exp(-myEta));

  // matrix
  CLHEP::HepRotation * fromCMStoTB = new CLHEP::HepRotation();
  double angle1 = 90.*deg - myPhi;
  CLHEP::HepRotationZ * r1 = new CLHEP::HepRotationZ(angle1);
  double angle2 = myTheta;
  CLHEP::HepRotationX * r2 = new CLHEP::HepRotationX(angle2);
  double angle3 = 90.*deg;
  CLHEP::HepRotationZ * r3 = new CLHEP::HepRotationZ(angle3);
  (*fromCMStoTB) *= (*r3);
  (*fromCMStoTB) *= (*r2);
  (*fromCMStoTB) *= (*r1);
  
  CMStoTB = (*fromCMStoTB);

  delete fromCMStoTB;
  delete r1;
  delete r2;
  delete r3;
}