Newer
Older
TB_Chris / TbIO / src / .svn / text-base / TbTupleWriter.cpp.svn-base
  1. // Tb/TbEvent
  2. #include "Event/TbTrigger.h"
  3. #include "Event/TbHit.h"
  4. #include "Event/TbCluster.h"
  5. #include "Event/TbTrack.h"
  6.  
  7. // Local
  8. #include "TbTupleWriter.h"
  9.  
  10. DECLARE_ALGORITHM_FACTORY(TbTupleWriter)
  11.  
  12. //=============================================================================
  13. // Standard constructor
  14. //=============================================================================
  15. TbTupleWriter::TbTupleWriter(const std::string& name, ISvcLocator* pSvcLocator)
  16. : TbAlgorithm(name, pSvcLocator) {
  17.  
  18. declareProperty("WriteTriggers", m_writeTriggers = false);
  19. declareProperty("WriteHits", m_writeHits = false);
  20. declareProperty("WriteClusters", m_writeClusters = true);
  21. declareProperty("WriteTracks", m_writeTracks = false);
  22.  
  23. declareProperty("MaxClusterSize", m_maxclustersize = 200);
  24. declareProperty("WriteClusterHits", m_writeClusterHits = true);
  25.  
  26. declareProperty("TriggerLocation",
  27. m_triggerLocation = LHCb::TbTriggerLocation::Default);
  28. declareProperty("HitLocation", m_hitLocation = LHCb::TbHitLocation::Default);
  29. declareProperty("ClusterLocation",
  30. m_clusterLocation = LHCb::TbClusterLocation::Default);
  31. declareProperty("TrackLocation",
  32. m_trackLocation = LHCb::TbTrackLocation::Default);
  33.  
  34. /// Stripped nTuples for low rate external users,
  35. /// only output tracks that have associated triggers
  36. declareProperty("StrippedNTuple", m_strippedNTuple = false);
  37. declareProperty("MaxTriggers", m_maxTriggers = 20 );
  38. // Switch off output during finalize.
  39. setProperty("NTuplePrint", false);
  40. }
  41.  
  42. //=============================================================================
  43. // Destructor
  44. //=============================================================================
  45. TbTupleWriter::~TbTupleWriter() {}
  46.  
  47. //=============================================================================
  48. // Initialisation
  49. //=============================================================================
  50. StatusCode TbTupleWriter::initialize() {
  51.  
  52. // Initialise the base class.
  53. StatusCode sc = TbAlgorithm::initialize();
  54. if (sc.isFailure()) return sc;
  55.  
  56. if (m_writeTriggers) bookTriggers();
  57. if (m_writeHits) bookHits();
  58. if (m_writeTracks) bookTracks();
  59. if (m_writeClusters) bookClusters();
  60. m_evtNo = 0;
  61. return StatusCode::SUCCESS;
  62. }
  63.  
  64. //=============================================================================
  65. // Main execution
  66. //=============================================================================
  67. StatusCode TbTupleWriter::execute() {
  68.  
  69. if (m_writeTriggers) fillTriggers();
  70. if (m_writeHits) fillHits();
  71. if (m_writeTracks) fillTracks();
  72. if (m_writeClusters) fillClusters();
  73. ++m_evtNo;
  74. return StatusCode::SUCCESS;
  75. }
  76.  
  77. //=============================================================================
  78. // Book trigger tuple
  79. //=============================================================================
  80. void TbTupleWriter::bookTriggers() {
  81.  
  82. NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1");
  83. NTuplePtr nt(ntupleSvc(), "/NTUPLES/FILE1/TbTupleWriter/Trigger");
  84. // Check if already booked.
  85. if (nt) return;
  86. nt = ntupleSvc()->book("/NTUPLES/FILE1/TbTupleWriter/Trigger",
  87. CLID_ColumnWiseTuple, "nTuple of Triggers");
  88. nt->addItem("TgID", m_TgID);
  89. nt->addItem("TgTime", m_TgTime);
  90. nt->addItem("TgHTime", m_TgHTime);
  91. nt->addItem("TgEvt", m_TgEvt);
  92. nt->addItem("TgPlane", m_TgPlane);
  93. nt->addItem("TgCounter",m_TgCounter);
  94. }
  95.  
  96. //=============================================================================
  97. // Book hit tuple
  98. //=============================================================================
  99. void TbTupleWriter::bookHits() {
  100.  
  101. NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1");
  102. NTuplePtr nt(ntupleSvc(), "/NTUPLES/FILE1/TbTupleWriter/Hits");
  103. // Check if already booked.
  104. if (nt) return;
  105. nt = ntupleSvc()->book("/NTUPLES/FILE1/TbTupleWriter/Hits",
  106. CLID_ColumnWiseTuple, "nTuple of Hits");
  107. nt->addItem("hID", m_hID);
  108. nt->addItem("hCol", m_hCol);
  109. nt->addItem("hRow", m_hRow);
  110. nt->addItem("hTime", m_hTime);
  111. nt->addItem("hHTime", m_hHTime);
  112. nt->addItem("hToT", m_hToT);
  113. nt->addItem("hPlane", m_hPlane);
  114. }
  115.  
  116. //=============================================================================
  117. // Book cluster tuple
  118. //=============================================================================
  119. void TbTupleWriter::bookClusters() {
  120.  
  121. NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1");
  122. NTuplePtr nt(ntupleSvc(), "/NTUPLES/FILE1/TbTupleWriter/Clusters");
  123. // Check if already booked.
  124. if (nt) return;
  125. nt = ntupleSvc()->book("/NTUPLES/FILE1/TbTupleWriter/Clusters",
  126. CLID_ColumnWiseTuple, "nTuple of Clusters");
  127. nt->addItem("clID", m_clID);
  128. nt->addItem("clGx", m_clGx);
  129. nt->addItem("clGy", m_clGy);
  130. nt->addItem("clGz", m_clGz);
  131. nt->addItem("clLx", m_clLx);
  132. nt->addItem("clLy", m_clLy);
  133. nt->addItem("clTime", m_clTime);
  134. nt->addItem("clHTime", m_clHTime);
  135. nt->addItem("clSize", m_clSize);
  136. nt->addItem("clCharge", m_clCharge);
  137. nt->addItem("clIsTracked", m_clTracked);
  138. nt->addItem("clPlane", m_clPlane);
  139. nt->addItem("clEvtNo", m_clEvtNo);
  140. nt->addItem("clNHits", m_clN, 0, (int)m_maxclustersize);
  141. nt->addIndexedItem("hRow", m_clN, m_clhRow);
  142. nt->addIndexedItem("hCol", m_clN, m_clhCol);
  143. nt->addIndexedItem("sCol", m_clN, m_clsCol);
  144. nt->addIndexedItem("hHTime", m_clN, m_clhHTime);
  145. nt->addIndexedItem("hToT", m_clN, m_clhToT);
  146. }
  147.  
  148. //=============================================================================
  149. // Book track tuple
  150. //=============================================================================
  151. void TbTupleWriter::bookTracks() {
  152.  
  153. NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1");
  154. NTuplePtr nt(ntupleSvc(), "/NTUPLES/FILE1/TbTupleWriter/Tracks");
  155. // Check if already booked.
  156. if (nt) return;
  157. nt = ntupleSvc()->book("/NTUPLES/FILE1/TbTupleWriter/Tracks",
  158. CLID_ColumnWiseTuple, "nTuple of Tracks");
  159. nt->addItem("TkID", m_TkID);
  160. nt->addItem("TkTime", m_TkTime);
  161. nt->addItem("TkHTime", m_TkHTime);
  162. nt->addItem("TkNCl", m_TkNCl, 0, 10);
  163. nt->addItem("TkX", m_TkX0);
  164. nt->addItem("TkY", m_TkY0);
  165. nt->addItem("TkTx", m_TkTx);
  166. nt->addItem("TkTy", m_TkTy);
  167. nt->addItem("TkChi2PerNdof", m_TkChi2ndof);
  168. nt->addIndexedItem("TkClId", m_TkNCl, m_TkClId);
  169. nt->addItem("TkEvt",m_TkEvt);
  170. nt->addItem("TkNTg",m_TkNTg,0,(int)m_maxTriggers);
  171. nt->addIndexedItem("TkTgId",m_TkNTg,m_TkTgId);
  172. nt->addIndexedItem("TkTgPlane",m_TkNTg,m_TkTgPlane);
  173. nt->addIndexedItem("TkXResidual",m_TkNCl, m_TkXResidual);
  174. nt->addIndexedItem("TkYResidual",m_TkNCl, m_TkYResidual);
  175.  
  176. }
  177.  
  178. //=============================================================================
  179. // Fill trigger tuple
  180. //=============================================================================
  181. void TbTupleWriter::fillTriggers() {
  182.  
  183. const uint64_t evtOffset = (uint64_t)m_evtNo << 36;
  184. for (unsigned int i = 0 ; i < m_nPlanes; ++i) {
  185. const std::string location = m_triggerLocation + std::to_string(i);
  186. const LHCb::TbTriggers* triggers = getIfExists<LHCb::TbTriggers>(location);
  187. if (!triggers) continue;
  188. for (const LHCb::TbTrigger* trigger : *triggers) {
  189. // if (m_triggerChannel != 999 && m_triggerChannel != trigger->plane()) continue;
  190. const uint64_t offset = ((uint64_t)trigger->plane() << 32) + evtOffset;
  191. m_TgID = trigger->index() + offset;
  192. m_TgTime = trigger->time();
  193. m_TgHTime = trigger->htime();
  194. m_TgPlane = i;
  195. m_TgCounter = trigger->counter();
  196. ntupleSvc()->writeRecord("/NTUPLES/FILE1/TbTupleWriter/Trigger");
  197. m_TgEvt = m_evtNo;
  198. }
  199. }
  200. }
  201.  
  202. //=============================================================================
  203. // Fill hit tuple
  204. //=============================================================================
  205. void TbTupleWriter::fillHits() {
  206.  
  207. const uint64_t evtOffset = (uint64_t)m_evtNo << 36;
  208. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  209. const std::string location = m_hitLocation + std::to_string(i);
  210. const LHCb::TbHits* hits = getIfExists<LHCb::TbHits>(location);
  211. if (!hits) continue;
  212. const uint64_t offset = ((uint64_t)i << 32) + evtOffset;
  213. for (const LHCb::TbHit* hit : *hits) {
  214. m_hID = hit->index() + offset;
  215. m_hCol = hit->col();
  216. m_hRow = hit->row();
  217. m_hTime = hit->time();
  218. m_hHTime = hit->htime();
  219. m_hToT = hit->ToT();
  220. m_hPlane = i;
  221. ntupleSvc()->writeRecord("/NTUPLES/FILE1/TbTupleWriter/Hits");
  222. }
  223. }
  224. }
  225.  
  226. //=============================================================================
  227. // Fill cluster tuple
  228. //=============================================================================
  229. void TbTupleWriter::fillClusters() {
  230.  
  231. const uint64_t evtOffset = (uint64_t)m_evtNo << 36;
  232. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  233. const std::string location = m_clusterLocation + std::to_string(i);
  234. const LHCb::TbClusters* clusters = getIfExists<LHCb::TbClusters>(location);
  235. if (!clusters) continue;
  236. const uint64_t offset = ((uint64_t)i << 32) + evtOffset;
  237. for (const LHCb::TbCluster* cluster : *clusters) {
  238. const unsigned long long clid = cluster->index() + offset;
  239. m_clID = clid;
  240. m_clGx = cluster->x();
  241. m_clGy = cluster->y();
  242. m_clGz = cluster->z();
  243. m_clLx = cluster->xloc();
  244. m_clLy = cluster->yloc();
  245. m_clTime = cluster->time();
  246. m_clHTime = cluster->htime();
  247. m_clSize = cluster->hits().size();
  248. m_clCharge = cluster->charge();
  249. m_clTracked = cluster->associated();
  250. m_clPlane = cluster->plane();
  251. m_clEvtNo = m_evtNo;
  252. m_clN = m_clSize > m_maxclustersize ? m_maxclustersize : m_clSize;
  253. for (unsigned int j = 0; j < m_clN; ++j) {
  254. const auto& h = cluster->hits()[j];
  255. m_clhRow[j] = h->row();
  256. m_clhCol[j] = h->col();
  257. m_clsCol[j] = h->scol();
  258. m_clhToT[j] = h->ToT();
  259. m_clhHTime[j] = h->htime();
  260. }
  261. ntupleSvc()->writeRecord("/NTUPLES/FILE1/TbTupleWriter/Clusters");
  262. }
  263. }
  264. }
  265.  
  266. //=============================================================================
  267. // Fill track tuple
  268. //=============================================================================
  269. void TbTupleWriter::fillTracks() {
  270.  
  271. const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  272. if (!tracks) return;
  273. const uint64_t evtOffset = (uint64_t)m_evtNo << 36;
  274. for (const LHCb::TbTrack* track : *tracks) {
  275. if (m_strippedNTuple && track->triggers().empty()) continue;
  276. m_TkX0 = track->firstState().x();
  277. m_TkY0 = track->firstState().y();
  278. m_TkTx = track->firstState().tx();
  279. m_TkTy = track->firstState().ty();
  280. m_TkHTime = track->htime();
  281. m_TkTime = track->time();
  282. m_TkID = track->index() + evtOffset;
  283. m_TkNCl = track->clusters().size();
  284. m_TkChi2ndof = track->chi2PerNdof();
  285. m_TkEvt = m_evtNo;
  286. unsigned int i = 0;
  287. const SmartRefVector<LHCb::TbCluster>& clusters = track->clusters();
  288. for (auto it = clusters.cbegin(), end = clusters.cend(); it != end; ++it) {
  289. const uint64_t offset = ((uint64_t)(*it)->plane() << 32) + evtOffset;
  290. const unsigned long long clid = (*it)->index() + offset;
  291. m_TkClId[i] = clid;
  292. auto intercept = geomSvc()->intercept( track, (*it)->plane() );
  293. m_TkXResidual[(*it)->plane()] = (*it)->x() - intercept.x() ;
  294. m_TkYResidual[(*it)->plane()] = (*it)->y() - intercept.y() ;
  295. ++i;
  296. }
  297. const SmartRefVector<LHCb::TbCluster>& associatedClusters = track->associatedClusters();
  298. for( auto it = associatedClusters.cbegin(), end = associatedClusters.cend(); it != end; ++it ){
  299. const uint64_t offset = ((uint64_t)(*it)->plane() << 32) + evtOffset;
  300. const unsigned long long clid = (*it)->index() + offset;
  301. m_TkClId[i] = clid;
  302. auto intercept = geomSvc()->intercept( track, (*it)->plane() );
  303. m_TkXResidual[(*it)->plane()] = (*it)->x() - intercept.x() ;
  304. m_TkYResidual[(*it)->plane()] = (*it)->y() - intercept.y() ;
  305. ++i;
  306. }
  307. m_TkNTg = track->triggers().size();
  308. i = 0;
  309. const SmartRefVector<LHCb::TbTrigger>& triggers = track->triggers();
  310. for (auto it = triggers.cbegin(), end = triggers.cend(); it != end; ++it) {
  311. const uint64_t offset = ((uint64_t)(*it)->plane() << 32) + evtOffset;
  312. const unsigned long long clid = (*it)->index() + offset;
  313. if (i == 20) {
  314. warning() << "More than 20 triggers associated to track: skipping"
  315. << endmsg;
  316. break;
  317. }
  318. m_TkTgId[i] = clid;
  319. m_TkTgPlane[i] = (*it)->plane();
  320. ++i;
  321. }
  322. ntupleSvc()->writeRecord("/NTUPLES/FILE1/TbTupleWriter/Tracks");
  323. }
  324. }