Newer
Older
TB_Chris / TbSimulation / src / .svn / text-base / TbTestMC.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3. #include "GaudiUtils/Aida2ROOT.h"
  4.  
  5. // Local
  6. #include "TbTestMC.h"
  7.  
  8. /** @file TbTestMC.cpp
  9. *
  10. * Implementation of class : TbTestMC
  11. * Author: Dan Saunders
  12. */
  13.  
  14. DECLARE_ALGORITHM_FACTORY(TbTestMC)
  15.  
  16. //=============================================================================
  17. /// Standard constructor
  18. //=============================================================================
  19. TbTestMC::TbTestMC(const std::string& name, ISvcLocator* pSvcLocator)
  20. : GaudiTupleAlg(name, pSvcLocator) {
  21. declareProperty("InitAlignment", m_filename = "Alignment_raw.dat");
  22. declareProperty("doMisAlign", m_misalign = false);
  23. declareProperty("NTracks", m_nTracks = 100); // Per event. Homogenous in t.
  24. declareProperty("RunDuration", m_EventLength = 10000); // Max clock time per event.
  25. declareProperty("NNoiseClusters", m_nNoiseClusters = 20); // Per event, per plane. Homogenous in t.
  26.  
  27.  
  28. declareProperty("HitTimeJitter", m_HitTimeJitter = 10); // Time jitter width (clocks).
  29. declareProperty("ClusterPosnError", m_ClustPosnResolution = 0.01);
  30.  
  31. declareProperty("ChargeSharingWidth", m_ChargeSharingWidth = 0.01);
  32. declareProperty("ThresholdCut", m_ThresholdCut = 50.0);
  33. declareProperty("ClusterCharge", m_charge = 300.0);
  34. declareProperty("ForceEfficiency", m_ForceEfficiency = true);
  35. // ...A measure of charge sharing - see below.
  36.  
  37. m_ChipWidth = 14.08; // mm
  38. m_PosnSpread = 4.0; // mm
  39.  
  40. m_ChargeSigma = 150;
  41. m_Pitch = 0.055;
  42. m_nEvents = 0;
  43. m_nExportedHits = 0;
  44. m_nExportedTracks = 0;
  45. m_ExportHits = true;
  46. m_ExportClusters = false;
  47. m_ExportTracks = false;
  48. }
  49.  
  50. //=============================================================================
  51. /// Initialization
  52. //=============================================================================
  53. StatusCode TbTestMC::initialize() {
  54.  
  55. StatusCode sc = GaudiAlgorithm::initialize();
  56. if (sc.isFailure()) return sc;
  57.  
  58. m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
  59. setHistoTopDir("Tb/");
  60.  
  61. if (m_misalign) {
  62. geomSvc()->readConditions(m_filename, m_modules);
  63. geomSvc()->printAlignment(m_modules);
  64. }
  65.  
  66. else m_modules = geomSvc()->modules();
  67. m_nPlanes = m_modules.size();
  68.  
  69. m_Hits.resize(m_nPlanes);
  70.  
  71. // Initialize the random number generators.
  72. sc = m_uniform.initialize(randSvc(), Rndm::Flat(0.0, 1.0));
  73. if (!sc) {
  74. error() << "Cannot initialize uniform random number generator." << endmsg;
  75. return sc;
  76. }
  77. sc = m_gauss.initialize(randSvc(), Rndm::Gauss(0.0, 1.0));
  78. if (!sc) {
  79. error() << "Cannot initialize Gaussian random number generator." << endmsg;
  80. return sc;
  81. }
  82. sc = m_landau.initialize(randSvc(), Rndm::Landau(m_charge, m_ChargeSigma));
  83. if (!sc) {
  84. error() << "Cannot initialize Landau random number generator." << endmsg;
  85. return sc;
  86. }
  87. sc = m_gauss2.initialize(randSvc(), Rndm::Gauss(0.0, 1.0e-4));
  88. if (!sc) {
  89. error() << "Cannot initialize Gaussian random number generator." << endmsg;
  90. return sc;
  91. }
  92. sc = m_uniform2.initialize(randSvc(), Rndm::Flat(0,m_ChipWidth ));
  93. if (!sc) {
  94. error() << "Cannot initialize uniform random number generator." << endmsg;
  95. return sc;
  96. }
  97. return StatusCode::SUCCESS;
  98. }
  99.  
  100. //=============================================================================
  101. /// Main execution
  102. //=============================================================================
  103. StatusCode TbTestMC::execute() {
  104. if (m_nEvents % 10 == 0)
  105. std::cout << "Event number: " << m_nEvents << std::endl;
  106.  
  107. generate_tracks();
  108. createClustFromTrack();
  109.  
  110. //make_noise();
  111. //make_tracks();
  112.  
  113. sort_stuff();
  114. add_to_TES();
  115.  
  116. m_nEvents++;
  117. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  118. m_nExportedHits += m_Hits[i].size();
  119. // Prepare for next event.
  120. m_Hits[i].clear();
  121. }
  122. m_Clusters.clear();
  123. m_Tracks.clear();
  124. return StatusCode::SUCCESS;
  125. }
  126.  
  127. //=============================================================================
  128. // Finalize
  129. //=============================================================================
  130. StatusCode TbTestMC::finalize() {
  131.  
  132. info() << "------------------------------------------------------" << endmsg;
  133. info() << " TbTestMC " << endmsg;
  134. info() << "------------------------------------------------------" << endmsg;
  135. info() << "Number of exported hits:\t"<<m_nExportedHits<<endmsg;
  136. info() << "------------------------------------------------------" << endmsg;
  137. // Finalise the base class.
  138. return GaudiTupleAlg::finalize();
  139.  
  140. }
  141.  
  142. //=============================================================================
  143. /// Sorting
  144. //=============================================================================
  145. void TbTestMC::add_to_TES() {
  146. // Hits.
  147. if (m_ExportHits) ExportHits();
  148.  
  149. // Clusters.
  150. if (m_ExportClusters) {
  151. LHCb::TbClusters* clusters = new LHCb::TbClusters();
  152. for (std::vector<LHCb::TbCluster*>::iterator ic = m_Clusters.begin();
  153. ic != m_Clusters.end(); ic++)
  154. clusters->insert(*ic);
  155.  
  156. put(clusters, LHCb::TbClusterLocation::Default);
  157. }
  158.  
  159. // Tracks.
  160. if (m_ExportTracks) {
  161. LHCb::TbTracks* tracks = new LHCb::TbTracks();
  162. for (std::vector<LHCb::TbTrack*>::iterator it = m_Tracks.begin();
  163. it != m_Tracks.end(); it++)
  164. tracks->insert(*it);
  165.  
  166. put(tracks, LHCb::TbTrackLocation::Default);
  167. }
  168. }
  169.  
  170. //=============================================================================
  171. /// Export hits to different TES locations.
  172. //=============================================================================
  173. void TbTestMC::ExportHits() {
  174. std::vector<LHCb::TbHits*> all_hits;
  175. for (unsigned int i = 0; i < m_nPlanes; i++) {
  176. LHCb::TbHits* hits = new LHCb::TbHits();
  177. put(hits, LHCb::TbHitLocation::Default + std::to_string(i));
  178. for (std::vector<LHCb::TbHit*>::iterator ih = m_Hits[i].begin();
  179. ih != m_Hits[i].end(); ih++) {
  180. hits->insert(*ih);
  181. }
  182. }
  183. }
  184.  
  185. //=============================================================================
  186. /// Sorting - no sorting if not exporting.
  187. //=============================================================================
  188. void TbTestMC::sort_stuff() {
  189. sort_by_chip();
  190. if (m_ExportHits) {
  191. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  192. hit_torder(i);
  193. }
  194. }
  195. if (m_ExportClusters) cluster_torder();
  196. if (m_ExportTracks) track_torder();
  197. }
  198.  
  199. //=============================================================================
  200. /// Sorting by chip
  201. //=============================================================================
  202. void TbTestMC::sort_by_chip() {
  203. std::vector<LHCb::TbCluster*> old_Clusters(m_Clusters);
  204. m_Clusters.clear();
  205. std::vector<TbModule*>::iterator m;
  206. for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
  207. std::vector<LHCb::TbCluster*>::iterator ic;
  208. for (ic = old_Clusters.begin(); ic != old_Clusters.end(); ++ic) {
  209. if ((*ic)->plane() == (unsigned int)iplane) m_Clusters.push_back(*ic);
  210. }
  211. }
  212.  
  213. }
  214.  
  215. //=============================================================================
  216. /// Track time ordering - honeycomb
  217. //=============================================================================
  218. void TbTestMC::track_torder() {
  219. int N = m_Tracks.size();
  220. float s_factor = 1.3;
  221. LHCb::TbTrack* track1, *track2;
  222. int gap = N / s_factor;
  223. bool swapped = false;
  224.  
  225. // Start the swap loops.
  226. while (gap > 1 || swapped) {
  227. if (gap > 1) gap /= s_factor;
  228.  
  229. swapped = false; // Reset per swap loop.
  230.  
  231. // Do the swaps.
  232. for (int i = 0; gap + i < N; ++i) {
  233. track1 = m_Tracks[i];
  234. track2 = m_Tracks[i + gap];
  235.  
  236. if (track1->time() > track2->time()) {
  237. // Do the swap.
  238. m_Tracks[i] = track2;
  239. m_Tracks[i + gap] = track1;
  240.  
  241. swapped = true;
  242. }
  243. }
  244. }
  245. }
  246.  
  247. //=============================================================================
  248. /// Hit time ordering - honeycomb
  249. //=============================================================================
  250. void TbTestMC::hit_torder(const unsigned int plane) {
  251. int N = m_Hits[plane].size();
  252. float s_factor = 1.3;
  253. LHCb::TbHit* hit1, *hit2;
  254. int gap = N / s_factor;
  255. bool swapped = false;
  256.  
  257. // Start the swap loops.
  258. while (gap > 1 || swapped) {
  259. if (gap > 1) gap /= s_factor;
  260.  
  261. swapped = false; // Reset per swap loop.
  262.  
  263. // Do the swaps.
  264. for (int i = 0; gap + i < N; ++i) {
  265. hit1 = m_Hits[plane][i];
  266. hit2 = m_Hits[plane][i + gap];
  267.  
  268. if (hit1->time() > hit2->time()) {
  269. // Do the swap.
  270. m_Hits[plane][i] = hit2;
  271. m_Hits[plane][i + gap] = hit1;
  272.  
  273. swapped = true;
  274. }
  275. }
  276. }
  277. }
  278.  
  279. //=============================================================================
  280. /// Cluster time ordering - honeycomb
  281. //=============================================================================
  282. void TbTestMC::cluster_torder() {
  283. int N = m_Clusters.size();
  284. float s_factor = 1.3;
  285. LHCb::TbCluster* clust1, *clust2;
  286. int gap = N / s_factor;
  287. bool swapped = false;
  288.  
  289. // Start the swap loops.
  290. while (gap > 1 || swapped) {
  291. if (gap > 1) gap /= s_factor;
  292.  
  293. swapped = false; // Reset per swap loop.
  294.  
  295. // Do the swaps.
  296. for (int i = 0; gap + i < N; ++i) {
  297. clust1 = m_Clusters[i];
  298. clust2 = m_Clusters[i + gap];
  299.  
  300. if (clust1->time() > clust2->time() &&
  301. clust1->plane() == clust2->plane()) {
  302.  
  303. // Do the swap.
  304. m_Clusters[i] = clust2;
  305. m_Clusters[i + gap] = clust1;
  306.  
  307. swapped = true;
  308. }
  309. }
  310. }
  311. }
  312.  
  313. //=============================================================================
  314. /// Noise cluster posn
  315. //=============================================================================
  316. Gaudi::XYZPoint TbTestMC::ClustGposn(int iplane) {
  317. float x, y;
  318. bool on = true;
  319. while (on) {
  320. x = 0.5 * m_ChipWidth + m_gauss() * m_PosnSpread;
  321. y = 0.5 * m_ChipWidth + m_gauss() * m_PosnSpread;
  322.  
  323. if (x > 0.028 + (2 * 0.055) && x < 14.05 - (2 * 0.055) &&
  324. y > 0.028 + (2 * 0.055) && y < 14.05 - (2 * 0.055))
  325. on = false;
  326. } // NO EDGES TAKEN INTO ACCOUNT YET.
  327.  
  328. Gaudi::XYZPoint posn(x, y, m_modules.at((unsigned int)iplane)->z());
  329. return posn;
  330. }
  331.  
  332. //=============================================================================
  333. /// Noise cluster charge
  334. //=============================================================================
  335. int TbTestMC::ClustCharge() {
  336. float charge;
  337. bool on = true;
  338. while (on) {
  339. charge = m_landau();
  340. if (charge > 3.0 && charge < 1000) on = false;
  341. }
  342. return (int)charge;
  343. }
  344.  
  345. //=============================================================================
  346. /// Tracks
  347. //=============================================================================
  348. void TbTestMC::make_tracks() {
  349. for (int i = 0; i < m_nTracks; i++) {
  350. // Make a track at a particular position. Clusters and hits saved auto.
  351. int time = (int)(m_uniform() * m_EventLength);
  352. LHCb::TbTrack* track = make_track(ClustGposn(0), time, i);
  353. m_Tracks.push_back(track);
  354.  
  355. unsigned int nActualClusters = 0;
  356. for (unsigned int i=0; i<track->clusters().size(); i++) {
  357. if (track->clusters()[i]->size() != 0) nActualClusters++;
  358. }
  359.  
  360. if (nActualClusters == m_nPlanes || nActualClusters == m_nPlanes-1)
  361. m_nExportedTracks++;
  362. }
  363. }
  364.  
  365. //=============================================================================
  366. /// Track
  367. //=============================================================================
  368. LHCb::TbTrack* TbTestMC::make_track(Gaudi::XYZPoint p, float t, int itrack) {
  369. LHCb::TbTrack* track = new LHCb::TbTrack();
  370. for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
  371. Gaudi::XYZPoint ClustPosn(p.x(), p.y(), m_modules.at(iplane)->z());
  372. LHCb::TbCluster* cluster = make_cluster(ClustPosn, t, iplane, itrack);
  373. track->addToClusters(cluster);
  374. m_Clusters.push_back(cluster);
  375. }
  376. m_trackFit->fit(track);
  377. return track;
  378. }
  379.  
  380. //=============================================================================
  381. /// Noise
  382. //=============================================================================
  383. void TbTestMC::make_noise() {
  384. for (unsigned int iplane = 0; iplane < m_nPlanes; iplane++) {
  385. for (int i = 0; i < m_nNoiseClusters; i++) {
  386. // Make a cluster at a particular position. Hits saved auto.
  387. int time = (int)(m_uniform() * m_EventLength);
  388. LHCb::TbCluster* cluster = make_cluster(ClustGposn(iplane), time, iplane, -1);
  389. m_Clusters.push_back(cluster);
  390. }
  391. }
  392. }
  393.  
  394. //=============================================================================
  395. /// MC cluster maker.
  396. //=============================================================================
  397. LHCb::TbCluster* TbTestMC::make_cluster(Gaudi::XYZPoint pGlobal_perfect,
  398. float t, int iplane, int track_tag) {
  399. LHCb::TbCluster* clust = new LHCb::TbCluster();
  400. clust->setTime(t);
  401.  
  402. double x = pGlobal_perfect.x();
  403. double y = pGlobal_perfect.y();
  404.  
  405. x += m_gauss() * m_ClustPosnResolution;
  406. y += m_gauss() * m_ClustPosnResolution;
  407.  
  408. Gaudi::XYZPoint pGlobal_actual(x, y, pGlobal_perfect.z());
  409.  
  410. Gaudi::XYZPoint pLocal = geomSvc()->globalToLocal(pGlobal_actual, iplane);
  411. clust->setXloc(pLocal.x());
  412. clust->setYloc(pLocal.y());
  413.  
  414. //std::cout<<clust->xloc()<<"\t"<<pGlobal_actual.x()<<"\t"<<clust->yloc()<<"\t"<<pGlobal_actual.y()<<std::endl;
  415.  
  416. clust->setX(pGlobal_actual.x());
  417. clust->setY(pGlobal_actual.y());
  418. clust->setZ(pGlobal_actual.z());
  419.  
  420. clust->setPlane(iplane);
  421. setClusterHitsAndCharge(clust, track_tag); // Also sets Charge.
  422.  
  423. return clust;
  424. }
  425.  
  426. //=============================================================================
  427. /// Cluster hits
  428. //=============================================================================
  429. void TbTestMC::setClusterHitsAndCharge(LHCb::TbCluster* clust, int /*track_tag*/) {
  430. // This version will model the charge cloud as a gaussian that is ceneted
  431. // at the cluster position, and whose height is set by the cluster charge
  432. // (albeit abit sneakily). Width is set above. The charge of surrounding
  433. // pixels is then taken to be the height of the gaussian at that pixels
  434. // position (if passing the threshold - kind of zero suppression);
  435. // like poor mans integration.
  436. double SigSq = m_ChargeSharingWidth*m_ChargeSharingWidth;
  437. int clustChargeish = ClustCharge();
  438. double N = clustChargeish*0.4;
  439. int SeedHitX = (int)(clust->xloc()/m_Pitch);
  440. int SeedHitY = (int)(clust->yloc()/m_Pitch);
  441.  
  442.  
  443. int ClustCharge = 0;
  444. int HalfAreaSide = 1; // 3x3
  445. // Search over area around the seed.
  446. //std::cout<<"New cluster at:\t"<<clust->xloc()<<"\t"<<clust->yloc()<<std::endl;
  447. for (int ix=SeedHitX-HalfAreaSide; ix!=SeedHitX+HalfAreaSide+1; ix++){
  448. for (int iy=SeedHitY-HalfAreaSide; iy!=SeedHitY+HalfAreaSide+1; iy++){
  449. // Find distance to the pixel from the center of the gaussian.
  450. double delx = (double(ix)+0.5)*m_Pitch - clust->xloc();
  451. double dely = (double(iy)+0.5)*m_Pitch - clust->yloc();
  452. double delrSq = delx*delx + dely*dely;
  453. int pixToT = (int)N*exp(-delrSq/SigSq);
  454.  
  455.  
  456. if ((float)pixToT > m_ThresholdCut) {
  457. LHCb::TbHit * hit = new LHCb::TbHit();
  458. hit->setDevice(clust->plane());
  459. hit->setCol(ix);
  460. hit->setRow(iy);
  461. hit->setToT(pixToT);
  462. hit->setTime(clust->time() + m_gauss() * m_HitTimeJitter);
  463. clust->addToHits(hit);
  464. m_Hits[clust->plane()].push_back(hit);
  465.  
  466. // hit->setTrackTag(track_tag);
  467. ClustCharge += pixToT;
  468. }
  469. }
  470. }
  471.  
  472. if (clust->size()==0 && m_ForceEfficiency){ // i.e. none of the hits passed the thresehold
  473. LHCb::TbHit * hit = new LHCb::TbHit();
  474. hit->setCol(SeedHitX);
  475. hit->setRow(SeedHitY);
  476. hit->setToT(clustChargeish);
  477. hit->setTime(clust->time() + m_gauss() * m_HitTimeJitter);
  478.  
  479. // hit->setTrackTag(track_tag);
  480.  
  481. clust->addToHits(hit);
  482. m_Hits[clust->plane()].push_back(hit);
  483. ClustCharge += clustChargeish;
  484. }
  485. clust->setCharge(ClustCharge);
  486. }
  487.  
  488. //=============================================================================
  489. // Randomly generates tracks (C.Hombach)
  490. //=============================================================================
  491.  
  492. void TbTestMC::generate_tracks()
  493. {
  494. double slope_xz;
  495. double slope_yz;
  496. double inter_x;
  497. double inter_y;
  498. for (int i = 0; i < m_nTracks; i++) {
  499. //Generate slopes and intercepts equally distributed over the sensor
  500. slope_xz = m_gauss2();
  501. slope_yz = m_gauss2();
  502. inter_x = m_uniform2();
  503. inter_y = m_uniform2();
  504. plot(slope_xz, "slopeXZ", "slopeXZ", -1.e-3,1.e-3,1000);
  505. plot(slope_yz, "slopeYZ", "slopeYZ", -1.e-3,1.e-3,1000);
  506. //Create track
  507. LHCb::TbState state;// = new LHCb::TbState();
  508. state.setTx(slope_xz);
  509. state.setTy(slope_yz);
  510. state.setX(inter_x);
  511. state.setY(inter_y);
  512. LHCb::TbTrack* track = new LHCb::TbTrack();
  513. track->setFirstState(state);
  514. m_Tracks.push_back(track);
  515. }
  516. }
  517.  
  518. //=============================================================================
  519. //Calculates cluster from track intercepts with module planes
  520. //=============================================================================
  521.  
  522. void TbTestMC::createClustFromTrack()
  523. {
  524. std::vector<LHCb::TbTrack*>::iterator it = m_Tracks.begin();
  525. for( ; it != m_Tracks.end(); ++it)
  526. {
  527. int time = (int)(m_uniform() * m_EventLength);
  528. std::vector<TbModule*>::iterator im = m_modules.begin();
  529. for( ; im != m_modules.end(); ++im)
  530. {
  531. Gaudi::XYZPoint ip = getIntercept(im - m_modules.begin(), (*it));
  532. Gaudi::XYZPoint pLocal = geomSvc()->globalToLocal(ip, im - m_modules.begin());
  533. LHCb::TbCluster* clust = new LHCb::TbCluster();
  534. clust->setTime(time);
  535. clust->setXloc(pLocal.x());
  536. clust->setYloc(pLocal.y());
  537. clust->setX(ip.x());
  538. clust->setY(ip.y());
  539. clust->setZ(ip.z());
  540. clust->setCharge(ClustCharge());
  541. clust->setPlane(im - m_modules.begin());
  542. setClusterHitsAndCharge(clust,-1);
  543. m_Clusters.push_back(clust);
  544. }
  545. }
  546. }
  547.  
  548. //============================================================================
  549. // Calculates Intercept of track and a module (Needs to be moved somewhere else)
  550. //============================================================================
  551.  
  552. Gaudi::XYZPoint TbTestMC::getIntercept(const unsigned int& plane, LHCb::TbTrack* track)
  553. {
  554. Gaudi::XYZPoint p1Local(0., 0., 0.), p2Local(0., 0., 1.);
  555. Gaudi::XYZPoint p1Global =
  556. geomSvc()->localToGlobal(p1Local, plane);
  557. Gaudi::XYZPoint p2Global =
  558. geomSvc()->localToGlobal(p2Local, plane);
  559. Gaudi::XYZPoint normal(p2Global.x() - p1Global.x(),
  560. p2Global.y() - p1Global.y(),
  561. p2Global.z() - p1Global.z());
  562. // std::cout << track->firstState().tx() << std::endl;
  563. const double dir_r =
  564. sqrt(track->firstState().tx() * track->firstState().tx() +
  565. track->firstState().ty() * track->firstState().ty() + 1.);
  566. const double dir_x = track->firstState().tx() / dir_r;
  567. const double dir_y = track->firstState().ty() / dir_r;
  568. const double dir_z = 1. / dir_r;
  569. const double length =
  570. ((p1Global.x() - track->firstState().x()) * normal.x() +
  571. (p1Global.y() - track->firstState().y()) * normal.y() +
  572. (p1Global.z() - 0.) * normal.z()) /
  573. (dir_x * normal.x() + dir_y * normal.y() + dir_z * normal.z());
  574. const double x_inter = track->firstState().x() + dir_x * length;
  575. const double y_inter = track->firstState().y() + dir_y * length;
  576. const double z_inter = 0. + dir_z * length;
  577. Gaudi::XYZPoint intersect_global(x_inter, y_inter, z_inter);
  578. return intersect_global;
  579. }
  580.  
  581.  
  582. //=============================================================================
  583. /// End
  584. //=============================================================================