// Gaudi
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiUtils/Aida2ROOT.h"
// Local
#include "TbTestMC.h"
/** @file TbTestMC.cpp
*
* Implementation of class : TbTestMC
* Author: Dan Saunders
*/
DECLARE_ALGORITHM_FACTORY(TbTestMC)
//=============================================================================
/// Standard constructor
//=============================================================================
TbTestMC::TbTestMC(const std::string& name, ISvcLocator* pSvcLocator)
: GaudiTupleAlg(name, pSvcLocator) {
declareProperty("InitAlignment", m_filename = "Alignment_raw.dat");
declareProperty("doMisAlign", m_misalign = false);
declareProperty("NTracks", m_nTracks = 100); // Per event. Homogenous in t.
declareProperty("RunDuration", m_EventLength = 10000); // Max clock time per event.
declareProperty("NNoiseClusters", m_nNoiseClusters = 20); // Per event, per plane. Homogenous in t.
declareProperty("HitTimeJitter", m_HitTimeJitter = 10); // Time jitter width (clocks).
declareProperty("ClusterPosnError", m_ClustPosnResolution = 0.01);
declareProperty("ChargeSharingWidth", m_ChargeSharingWidth = 0.01);
declareProperty("ThresholdCut", m_ThresholdCut = 50.0);
declareProperty("ClusterCharge", m_charge = 300.0);
declareProperty("ForceEfficiency", m_ForceEfficiency = true);
// ...A measure of charge sharing - see below.
m_ChipWidth = 14.08; // mm
m_PosnSpread = 4.0; // mm
m_ChargeSigma = 150;
m_Pitch = 0.055;
m_nEvents = 0;
m_nExportedHits = 0;
m_nExportedTracks = 0;
m_ExportHits = true;
m_ExportClusters = false;
m_ExportTracks = false;
}
//=============================================================================
/// Initialization
//=============================================================================
StatusCode TbTestMC::initialize() {
StatusCode sc = GaudiAlgorithm::initialize();
if (sc.isFailure()) return sc;
m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
setHistoTopDir("Tb/");
if (m_misalign) {
geomSvc()->readConditions(m_filename, m_modules);
geomSvc()->printAlignment(m_modules);
}
else m_modules = geomSvc()->modules();
m_nPlanes = m_modules.size();
m_Hits.resize(m_nPlanes);
// Initialize the random number generators.
sc = m_uniform.initialize(randSvc(), Rndm::Flat(0.0, 1.0));
if (!sc) {
error() << "Cannot initialize uniform random number generator." << endmsg;
return sc;
}
sc = m_gauss.initialize(randSvc(), Rndm::Gauss(0.0, 1.0));
if (!sc) {
error() << "Cannot initialize Gaussian random number generator." << endmsg;
return sc;
}
sc = m_landau.initialize(randSvc(), Rndm::Landau(m_charge, m_ChargeSigma));
if (!sc) {
error() << "Cannot initialize Landau random number generator." << endmsg;
return sc;
}
sc = m_gauss2.initialize(randSvc(), Rndm::Gauss(0.0, 1.0e-4));
if (!sc) {
error() << "Cannot initialize Gaussian random number generator." << endmsg;
return sc;
}
sc = m_uniform2.initialize(randSvc(), Rndm::Flat(0,m_ChipWidth ));
if (!sc) {
error() << "Cannot initialize uniform random number generator." << endmsg;
return sc;
}
return StatusCode::SUCCESS;
}
//=============================================================================
/// Main execution
//=============================================================================
StatusCode TbTestMC::execute() {
if (m_nEvents % 10 == 0)
std::cout << "Event number: " << m_nEvents << std::endl;
generate_tracks();
createClustFromTrack();
//make_noise();
//make_tracks();
sort_stuff();
add_to_TES();
m_nEvents++;
for (unsigned int i = 0; i < m_nPlanes; ++i) {
m_nExportedHits += m_Hits[i].size();
// Prepare for next event.
m_Hits[i].clear();
}
m_Clusters.clear();
m_Tracks.clear();
return StatusCode::SUCCESS;
}
//=============================================================================
// Finalize
//=============================================================================
StatusCode TbTestMC::finalize() {
info() << "------------------------------------------------------" << endmsg;
info() << " TbTestMC " << endmsg;
info() << "------------------------------------------------------" << endmsg;
info() << "Number of exported hits:\t"<<m_nExportedHits<<endmsg;
info() << "------------------------------------------------------" << endmsg;
// Finalise the base class.
return GaudiTupleAlg::finalize();
}
//=============================================================================
/// Sorting
//=============================================================================
void TbTestMC::add_to_TES() {
// Hits.
if (m_ExportHits) ExportHits();
// Clusters.
if (m_ExportClusters) {
LHCb::TbClusters* clusters = new LHCb::TbClusters();
for (std::vector<LHCb::TbCluster*>::iterator ic = m_Clusters.begin();
ic != m_Clusters.end(); ic++)
clusters->insert(*ic);
put(clusters, LHCb::TbClusterLocation::Default);
}
// Tracks.
if (m_ExportTracks) {
LHCb::TbTracks* tracks = new LHCb::TbTracks();
for (std::vector<LHCb::TbTrack*>::iterator it = m_Tracks.begin();
it != m_Tracks.end(); it++)
tracks->insert(*it);
put(tracks, LHCb::TbTrackLocation::Default);
}
}
//=============================================================================
/// Export hits to different TES locations.
//=============================================================================
void TbTestMC::ExportHits() {
std::vector<LHCb::TbHits*> all_hits;
for (unsigned int i = 0; i < m_nPlanes; i++) {
LHCb::TbHits* hits = new LHCb::TbHits();
put(hits, LHCb::TbHitLocation::Default + std::to_string(i));
for (std::vector<LHCb::TbHit*>::iterator ih = m_Hits[i].begin();
ih != m_Hits[i].end(); ih++) {
hits->insert(*ih);
}
}
}
//=============================================================================
/// Sorting - no sorting if not exporting.
//=============================================================================
void TbTestMC::sort_stuff() {
sort_by_chip();
if (m_ExportHits) {
for (unsigned int i = 0; i < m_nPlanes; ++i) {
hit_torder(i);
}
}
if (m_ExportClusters) cluster_torder();
if (m_ExportTracks) track_torder();
}
//=============================================================================
/// Sorting by chip
//=============================================================================
void TbTestMC::sort_by_chip() {
std::vector<LHCb::TbCluster*> old_Clusters(m_Clusters);
m_Clusters.clear();
std::vector<TbModule*>::iterator m;
for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
std::vector<LHCb::TbCluster*>::iterator ic;
for (ic = old_Clusters.begin(); ic != old_Clusters.end(); ++ic) {
if ((*ic)->plane() == (unsigned int)iplane) m_Clusters.push_back(*ic);
}
}
}
//=============================================================================
/// Track time ordering - honeycomb
//=============================================================================
void TbTestMC::track_torder() {
int N = m_Tracks.size();
float s_factor = 1.3;
LHCb::TbTrack* track1, *track2;
int gap = N / s_factor;
bool swapped = false;
// Start the swap loops.
while (gap > 1 || swapped) {
if (gap > 1) gap /= s_factor;
swapped = false; // Reset per swap loop.
// Do the swaps.
for (int i = 0; gap + i < N; ++i) {
track1 = m_Tracks[i];
track2 = m_Tracks[i + gap];
if (track1->time() > track2->time()) {
// Do the swap.
m_Tracks[i] = track2;
m_Tracks[i + gap] = track1;
swapped = true;
}
}
}
}
//=============================================================================
/// Hit time ordering - honeycomb
//=============================================================================
void TbTestMC::hit_torder(const unsigned int plane) {
int N = m_Hits[plane].size();
float s_factor = 1.3;
LHCb::TbHit* hit1, *hit2;
int gap = N / s_factor;
bool swapped = false;
// Start the swap loops.
while (gap > 1 || swapped) {
if (gap > 1) gap /= s_factor;
swapped = false; // Reset per swap loop.
// Do the swaps.
for (int i = 0; gap + i < N; ++i) {
hit1 = m_Hits[plane][i];
hit2 = m_Hits[plane][i + gap];
if (hit1->time() > hit2->time()) {
// Do the swap.
m_Hits[plane][i] = hit2;
m_Hits[plane][i + gap] = hit1;
swapped = true;
}
}
}
}
//=============================================================================
/// Cluster time ordering - honeycomb
//=============================================================================
void TbTestMC::cluster_torder() {
int N = m_Clusters.size();
float s_factor = 1.3;
LHCb::TbCluster* clust1, *clust2;
int gap = N / s_factor;
bool swapped = false;
// Start the swap loops.
while (gap > 1 || swapped) {
if (gap > 1) gap /= s_factor;
swapped = false; // Reset per swap loop.
// Do the swaps.
for (int i = 0; gap + i < N; ++i) {
clust1 = m_Clusters[i];
clust2 = m_Clusters[i + gap];
if (clust1->time() > clust2->time() &&
clust1->plane() == clust2->plane()) {
// Do the swap.
m_Clusters[i] = clust2;
m_Clusters[i + gap] = clust1;
swapped = true;
}
}
}
}
//=============================================================================
/// Noise cluster posn
//=============================================================================
Gaudi::XYZPoint TbTestMC::ClustGposn(int iplane) {
float x, y;
bool on = true;
while (on) {
x = 0.5 * m_ChipWidth + m_gauss() * m_PosnSpread;
y = 0.5 * m_ChipWidth + m_gauss() * m_PosnSpread;
if (x > 0.028 + (2 * 0.055) && x < 14.05 - (2 * 0.055) &&
y > 0.028 + (2 * 0.055) && y < 14.05 - (2 * 0.055))
on = false;
} // NO EDGES TAKEN INTO ACCOUNT YET.
Gaudi::XYZPoint posn(x, y, m_modules.at((unsigned int)iplane)->z());
return posn;
}
//=============================================================================
/// Noise cluster charge
//=============================================================================
int TbTestMC::ClustCharge() {
float charge;
bool on = true;
while (on) {
charge = m_landau();
if (charge > 3.0 && charge < 1000) on = false;
}
return (int)charge;
}
//=============================================================================
/// Tracks
//=============================================================================
void TbTestMC::make_tracks() {
for (int i = 0; i < m_nTracks; i++) {
// Make a track at a particular position. Clusters and hits saved auto.
int time = (int)(m_uniform() * m_EventLength);
LHCb::TbTrack* track = make_track(ClustGposn(0), time, i);
m_Tracks.push_back(track);
unsigned int nActualClusters = 0;
for (unsigned int i=0; i<track->clusters().size(); i++) {
if (track->clusters()[i]->size() != 0) nActualClusters++;
}
if (nActualClusters == m_nPlanes || nActualClusters == m_nPlanes-1)
m_nExportedTracks++;
}
}
//=============================================================================
/// Track
//=============================================================================
LHCb::TbTrack* TbTestMC::make_track(Gaudi::XYZPoint p, float t, int itrack) {
LHCb::TbTrack* track = new LHCb::TbTrack();
for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
Gaudi::XYZPoint ClustPosn(p.x(), p.y(), m_modules.at(iplane)->z());
LHCb::TbCluster* cluster = make_cluster(ClustPosn, t, iplane, itrack);
track->addToClusters(cluster);
m_Clusters.push_back(cluster);
}
m_trackFit->fit(track);
return track;
}
//=============================================================================
/// Noise
//=============================================================================
void TbTestMC::make_noise() {
for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
for (int i = 0; i < m_nNoiseClusters; i++) {
// Make a cluster at a particular position. Hits saved auto.
int time = (int)(m_uniform() * m_EventLength);
LHCb::TbCluster* cluster = make_cluster(ClustGposn(iplane), time, iplane, -1);
m_Clusters.push_back(cluster);
}
}
}
//=============================================================================
/// MC cluster maker.
//=============================================================================
LHCb::TbCluster* TbTestMC::make_cluster(Gaudi::XYZPoint pGlobal_perfect,
float t, int iplane, int track_tag) {
LHCb::TbCluster* clust = new LHCb::TbCluster();
clust->setTime(t);
double x = pGlobal_perfect.x();
double y = pGlobal_perfect.y();
x += m_gauss() * m_ClustPosnResolution;
y += m_gauss() * m_ClustPosnResolution;
Gaudi::XYZPoint pGlobal_actual(x, y, pGlobal_perfect.z());
Gaudi::XYZPoint pLocal = geomSvc()->globalToLocal(pGlobal_actual, iplane);
clust->setXloc(pLocal.x());
clust->setYloc(pLocal.y());
//std::cout<<clust->xloc()<<"\t"<<pGlobal_actual.x()<<"\t"<<clust->yloc()<<"\t"<<pGlobal_actual.y()<<std::endl;
clust->setX(pGlobal_actual.x());
clust->setY(pGlobal_actual.y());
clust->setZ(pGlobal_actual.z());
clust->setPlane(iplane);
setClusterHitsAndCharge(clust, track_tag); // Also sets Charge.
return clust;
}
//=============================================================================
/// Cluster hits
//=============================================================================
void TbTestMC::setClusterHitsAndCharge(LHCb::TbCluster* clust, int /*track_tag*/) {
// This version will model the charge cloud as a gaussian that is ceneted
// at the cluster position, and whose height is set by the cluster charge
// (albeit abit sneakily). Width is set above. The charge of surrounding
// pixels is then taken to be the height of the gaussian at that pixels
// position (if passing the threshold - kind of zero suppression);
// like poor mans integration.
double SigSq = m_ChargeSharingWidth*m_ChargeSharingWidth;
int clustChargeish = ClustCharge();
double N = clustChargeish*0.4;
int SeedHitX = (int)(clust->xloc()/m_Pitch);
int SeedHitY = (int)(clust->yloc()/m_Pitch);
int ClustCharge = 0;
int HalfAreaSide = 1; // 3x3
// Search over area around the seed.
//std::cout<<"New cluster at:\t"<<clust->xloc()<<"\t"<<clust->yloc()<<std::endl;
for (int ix=SeedHitX-HalfAreaSide; ix!=SeedHitX+HalfAreaSide+1; ix++){
for (int iy=SeedHitY-HalfAreaSide; iy!=SeedHitY+HalfAreaSide+1; iy++){
// Find distance to the pixel from the center of the gaussian.
double delx = (double(ix)+0.5)*m_Pitch - clust->xloc();
double dely = (double(iy)+0.5)*m_Pitch - clust->yloc();
double delrSq = delx*delx + dely*dely;
int pixToT = (int)N*exp(-delrSq/SigSq);
if ((float)pixToT > m_ThresholdCut) {
LHCb::TbHit * hit = new LHCb::TbHit();
hit->setDevice(clust->plane());
hit->setCol(ix);
hit->setRow(iy);
hit->setToT(pixToT);
hit->setTime(clust->time() + m_gauss() * m_HitTimeJitter);
clust->addToHits(hit);
m_Hits[clust->plane()].push_back(hit);
// hit->setTrackTag(track_tag);
ClustCharge += pixToT;
}
}
}
if (clust->size()==0 && m_ForceEfficiency){ // i.e. none of the hits passed the thresehold
LHCb::TbHit * hit = new LHCb::TbHit();
hit->setCol(SeedHitX);
hit->setRow(SeedHitY);
hit->setToT(clustChargeish);
hit->setTime(clust->time() + m_gauss() * m_HitTimeJitter);
// hit->setTrackTag(track_tag);
clust->addToHits(hit);
m_Hits[clust->plane()].push_back(hit);
ClustCharge += clustChargeish;
}
clust->setCharge(ClustCharge);
}
//=============================================================================
// Randomly generates tracks (C.Hombach)
//=============================================================================
void TbTestMC::generate_tracks()
{
double slope_xz;
double slope_yz;
double inter_x;
double inter_y;
for (int i = 0; i < m_nTracks; i++) {
//Generate slopes and intercepts equally distributed over the sensor
slope_xz = m_gauss2();
slope_yz = m_gauss2();
inter_x = m_uniform2();
inter_y = m_uniform2();
plot(slope_xz, "slopeXZ", "slopeXZ", -1.e-3,1.e-3,1000);
plot(slope_yz, "slopeYZ", "slopeYZ", -1.e-3,1.e-3,1000);
//Create track
LHCb::TbState state;// = new LHCb::TbState();
state.setTx(slope_xz);
state.setTy(slope_yz);
state.setX(inter_x);
state.setY(inter_y);
LHCb::TbTrack* track = new LHCb::TbTrack();
track->setFirstState(state);
m_Tracks.push_back(track);
}
}
//=============================================================================
//Calculates cluster from track intercepts with module planes
//=============================================================================
void TbTestMC::createClustFromTrack()
{
std::vector<LHCb::TbTrack*>::iterator it = m_Tracks.begin();
for( ; it != m_Tracks.end(); ++it)
{
int time = (int)(m_uniform() * m_EventLength);
std::vector<TbModule*>::iterator im = m_modules.begin();
for( ; im != m_modules.end(); ++im)
{
Gaudi::XYZPoint ip = getIntercept(im - m_modules.begin(), (*it));
Gaudi::XYZPoint pLocal = geomSvc()->globalToLocal(ip, im - m_modules.begin());
LHCb::TbCluster* clust = new LHCb::TbCluster();
clust->setTime(time);
clust->setXloc(pLocal.x());
clust->setYloc(pLocal.y());
clust->setX(ip.x());
clust->setY(ip.y());
clust->setZ(ip.z());
clust->setCharge(ClustCharge());
clust->setPlane(im - m_modules.begin());
setClusterHitsAndCharge(clust,-1);
m_Clusters.push_back(clust);
}
}
}
//============================================================================
// Calculates Intercept of track and a module (Needs to be moved somewhere else)
//============================================================================
Gaudi::XYZPoint TbTestMC::getIntercept(const unsigned int& plane, LHCb::TbTrack* track)
{
Gaudi::XYZPoint p1Local(0., 0., 0.), p2Local(0., 0., 1.);
Gaudi::XYZPoint p1Global =
geomSvc()->localToGlobal(p1Local, plane);
Gaudi::XYZPoint p2Global =
geomSvc()->localToGlobal(p2Local, plane);
Gaudi::XYZPoint normal(p2Global.x() - p1Global.x(),
p2Global.y() - p1Global.y(),
p2Global.z() - p1Global.z());
// std::cout << track->firstState().tx() << std::endl;
const double dir_r =
sqrt(track->firstState().tx() * track->firstState().tx() +
track->firstState().ty() * track->firstState().ty() + 1.);
const double dir_x = track->firstState().tx() / dir_r;
const double dir_y = track->firstState().ty() / dir_r;
const double dir_z = 1. / dir_r;
const double length =
((p1Global.x() - track->firstState().x()) * normal.x() +
(p1Global.y() - track->firstState().y()) * normal.y() +
(p1Global.z() - 0.) * normal.z()) /
(dir_x * normal.x() + dir_y * normal.y() + dir_z * normal.z());
const double x_inter = track->firstState().x() + dir_x * length;
const double y_inter = track->firstState().y() + dir_y * length;
const double z_inter = 0. + dir_z * length;
Gaudi::XYZPoint intersect_global(x_inter, y_inter, z_inter);
return intersect_global;
}
//=============================================================================
/// End
//=============================================================================
iaro