Belle II Software  release-05-02-19
TRGGRLMatchModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Chunhua LI, Yun-Tsung Lai, Junhao Yin *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // include
12 #include <trg/grl/modules/trggrl/TRGGRLMatchModule.h>
13 #include <trg/grl/dataobjects/TRGGRLMATCH.h>
14 #include <trg/grl/dataobjects/TRGGRLMATCHKLM.h>
15 #include <trg/grl/dataobjects/TRGGRLPHOTON.h>
16 #include <trg/cdc/dataobjects/CDCTriggerTrack.h>
17 #include <trg/ecl/dataobjects/TRGECLCluster.h>
18 #include <trg/klm/dataobjects/KLMTriggerTrack.h>
19 #include <trg/cdc/dataobjects/CDCTriggerSegmentHit.h>
20 #include <trg/grl/dataobjects/TRGGRLInfo.h>
21 #include <trg/grl/dataobjects/TRGGRLShortTrack.h>
22 
23 // framework - DataStore
24 #include <framework/datastore/StoreArray.h>
25 #include <framework/datastore/StoreObjPtr.h>
26 
27 //framework aux
28 #include <framework/logging/Logger.h>
29 #include <framework/core/ModuleParamList.templateDetails.h>
30 
31 #include <stdlib.h>
32 
33 
34 using namespace Belle2;
35 //-----------------------------------------------------------------
36 // Register the Module
37 //-----------------------------------------------------------------
39 
40 //-----------------------------------------------------------------
41 // Implementation
42 //-----------------------------------------------------------------
43 
45 {
46  // Set module properties
47  setDescription("match CDC trigger tracks and ECL trigger clusters");
48  setPropertyFlags(c_ParallelProcessingCertified);
49  addParam("SimulationMode", m_simulationMode, "TRGGRL simulation switch", 1);
50  addParam("FastSimulationMode", m_fastSimulationMode, "TRGGRL fast simulation mode", m_fastSimulationMode);
51  addParam("FirmwareSimulationMode", m_firmwareSimulationMode, "TRGGRL firmware simulation mode", m_firmwareSimulationMode);
52 
53  addParam("DrMatch", m_dr_threshold, "the threshold of dr between track and cluster if they are matched successfully", 25.);
54  addParam("DzMatch", m_dz_threshold, "the threshold of dz between track and cluster if they are matched successfully", 30.);
55  addParam("DphidMatch", m_dphi_d_threshold, "the threshold of dphi_d between track and cluster if they are matched successfully", 2);
56  addParam("Ephoton", m_e_threshold, "the threshold of cluster energy as a photon", 1.0);
57  addParam("KLMMatch", m_dphi_klm_threshold,
58  "the threshold of dphi (in degree) between track and KLM sector if they are matched successfully", 65.0);
59  addParam("2DtrackCollection", m_2d_tracklist, "the 2d track list used in the match", std::string("TRGCDC2DFinderTracks"));
60  addParam("3DtrackCollection", m_3d_tracklist, "the 3d track list used in the match", std::string("TRGCDCNeuroTracks"));
61  addParam("TRGECLClusterCollection", m_clusterlist, "the cluster list used in the match", std::string("TRGECLClusters"));
62  addParam("KLMTriggerTrack", m_klmtracklist, "the KLM track list used in the match", std::string("TRGKLMTracks"));
63  addParam("2DmatchCollection", m_2dmatch_tracklist, "the 2d tracklist with associated cluster", std::string("TRG2DMatchTracks"));
64  addParam("PhimatchCollection", m_phimatch_tracklist, "the 2d tracklist with associated cluster", std::string("TRGPhiMatchTracks"));
65  addParam("3DmatchCollection", m_3dmatch_tracklist, "the 3d NN tracklist with associated cluster", std::string("TRG3DMatchTracks"));
66  addParam("KLMmatchCollection", m_klmmatch_tracklist, "the 2d tracklist with associated KLM track",
67  std::string("TRGKLMMatchTracks"));
68  addParam("GRLphotonCollection", m_grlphotonlist, "the isolated cluster list", std::string("TRGGRLPhotons"));
69  addParam("hitCollectionName", m_hitCollectionName, "Name of the input StoreArray of CDCTriggerSegmentHits.", std::string(""));
70  addParam("TrgGrlInformation", m_TrgGrlInformationName,
71  "Name of the StoreArray holding the information of tracks and clusters from cdc ecl klm.",
72  std::string("TRGGRLObjects"));
73  addParam("grlstCollectionName", m_grlstCollectionName, "Name of the output StoreArray of TRGGRLShortTrack.",
74  std::string("TRGGRLShortTracks"));
75 
76 
77 }
78 
80 {
81 }
82 
84 {
85  B2DEBUG(100, "TRGGRLMatchModule processing");
88  track2Dlist.isRequired();
89  track3Dlist.isRequired();
91  clusterslist.isRequired();
92  clusterslist.registerRelationTo(track2Dlist);
93  clusterslist.registerRelationTo(track3Dlist);
95  klmtracklist.isRequired();
96  klmtracklist.registerRelationTo(track2Dlist);
97 
99  tslist.isRequired();
100 
101 
102 // track2Dlist.registerRelationTo(clusterslist);
103 // track3Dlist.registerRelationTo(clusterslist);
104 
105 // modified by ytlai 2017/12/15: registerPersistent should be replaced by registerInDataStore
106 // StoreArray<TRGGRLMATCH>::registerPersistent(m_2dmatch_tracklist);
107 // StoreArray<TRGGRLMATCH>::registerPersistent(m_phimatch_tracklist);
108 // StoreArray<TRGGRLMATCH>::register(m_2dmatch_tracklist);
109 // StoreArray<TRGGRLMATCH>::registerPersistent(m_3dmatch_tracklist);
110 // StoreArray<TRGGRLMATCHKLM>::registerPersistent(m_klmmatch_tracklist);
111 
112 // StoreArray<TRGGRLMATCH> track2Dmatch(m_2dmatch_tracklist);
113  StoreArray<TRGGRLMATCH> track2Dmatch;
114  track2Dmatch.registerInDataStore(m_2dmatch_tracklist);
115  track2Dmatch.registerRelationTo(track2Dlist);
116  track2Dmatch.registerRelationTo(clusterslist);
117 
118 // StoreArray<TRGGRLMATCH> trackphimatch(m_phimatch_tracklist);
119  StoreArray<TRGGRLMATCH> trackphimatch;
120  trackphimatch.registerInDataStore(m_phimatch_tracklist);
121  trackphimatch.registerRelationTo(track2Dlist);
122  trackphimatch.registerRelationTo(clusterslist);
123 
124 // StoreArray<TRGGRLMATCH> track3Dmatch(m_3dmatch_tracklist);
125  StoreArray<TRGGRLMATCH> track3Dmatch;
126  track3Dmatch.registerInDataStore(m_3dmatch_tracklist);
127  track3Dmatch.registerRelationTo(clusterslist);
128  track3Dmatch.registerRelationTo(track3Dlist);
129 
130 // StoreArray<TRGGRLMATCHKLM> trackKLMmatch(m_klmmatch_tracklist);
131  StoreArray<TRGGRLMATCHKLM> trackKLMmatch;
132  trackKLMmatch.registerInDataStore(m_klmmatch_tracklist);
133  trackKLMmatch.registerRelationTo(track2Dlist);
134  trackKLMmatch.registerRelationTo(klmtracklist);
135 
136  StoreArray<TRGGRLPHOTON> grlphoton;
137  grlphoton.registerInDataStore(m_grlphotonlist);
138  grlphoton.registerRelationTo(clusterslist);
139 
141  grlst.registerInDataStore(m_grlstCollectionName);
142 
143  m_TRGGRLInfo.registerInDataStore(m_TrgGrlInformationName);
144 
145 //-- Fill the patterns for short tracking
146 
148 
149  for (int p = 0; p < 137; p++) {
150  int x0 = patterns_base2[p][0];
151  int x1 = patterns_base2[p][1];
152  int x2 = 0;
153  int x3 = patterns_base2[p][2];
154  int x4 = patterns_base2[p][3];
155  int d = x2 - x0;
156  x1 += d;
157  x2 += d;
158  x3 += d;
159  x4 += d;
160  patterns_base0.push_back({x1, x2, x3, x4});
161  }
162 
163 
164 }
165 
167 {
168 }
169 
171 {
172 
185  trgInfo.create();
186 
187 //initialize the phi map
188 
189  track_phimap.clear();
190  track_phimap_i.clear();
191 
192  for (int i = 0; i < 36; i++) {
193  track_phimap.push_back(false);
194  track_phimap_i.push_back(false);
195  }
196 
197 //do 2d track match with cluster
198  for (int i = 0; i < track2Dlist.getEntries(); i++) {
199 
200  double dr_tmp = 99999.;
201  int dphi_d_tmp = 100;
202  int dphi_klm_tmp = 100;
203  int cluster_ind = -1;
204  int cluster_ind_phi = -1;
205  int klmtrack_ind_phi = -1;
206 
207  // do 2d track match with KLM track
208  for (int j = 0; j < klmtracklist.getEntries(); j++) {
209  double dphi_klm = 99999.9;
210  sectormatching_klm(track2Dlist[i], klmtracklist[j], dphi_klm);
211  if (dphi_klm_tmp > dphi_klm) {
212  dphi_klm_tmp = dphi_klm;
213  klmtrack_ind_phi = j;
214  }
215  }
216 
217  for (int j = 0; j < clusterlist.getEntries(); j++) {
218  // skip the end-cap cluster
219  double _cluster_x = clusterlist[j]->getPositionX();
220  double _cluster_y = clusterlist[j]->getPositionY();
221  double _cluster_z = clusterlist[j]->getPositionZ();
222  double _cluster_theta = atan(_cluster_z / (sqrt(_cluster_x * _cluster_x + _cluster_y * _cluster_y)));
223  _cluster_theta = 0.5 * M_PI - _cluster_theta;
224  if (_cluster_theta < M_PI * 35.0 / 180.0 || _cluster_theta > M_PI * 126.0 / 180.0) continue;
225 
226  double ds_ct[2] = {99999., 99999.};
227  calculationdistance(track2Dlist[i], clusterlist[j], ds_ct, 0);
228  int dphi_d = 0;
229  calculationphiangle(track2Dlist[i], clusterlist[j], dphi_d, track_phimap, track_phimap_i);
230 
231  if (dr_tmp > ds_ct[0]) {
232  dr_tmp = ds_ct[0];
233  cluster_ind = j;
234  }
235  if (dphi_d_tmp > dphi_d) {
236  dphi_d_tmp = dphi_d;
237  cluster_ind_phi = j;
238  }
239 
240  }
241 
242  if (dr_tmp < m_dr_threshold && cluster_ind != -1) {
243  TRGGRLMATCH* mat2d = track2Dmatch.appendNew();
244  mat2d->setDeltaR(dr_tmp);
245  mat2d->addRelationTo(track2Dlist[i]);
246  mat2d->addRelationTo(clusterlist[cluster_ind]);
247  // track2Dlist[i]->addRelationTo(clusterlist[cluster_ind]);
248  clusterlist[cluster_ind]->addRelationTo(track2Dlist[i]);
249  }
250  if (dphi_d_tmp < m_dphi_d_threshold && cluster_ind_phi != -1) {
251  TRGGRLMATCH* matphi = trackphimatch.appendNew();
252  matphi->set_dphi_d(dphi_d_tmp);
253  matphi->addRelationTo(track2Dlist[i]);
254  matphi->addRelationTo(clusterlist[cluster_ind_phi]);
255  matphi->set_e(clusterlist[cluster_ind_phi]->getEnergyDep());
256  // track2Dlist[i]->addRelationTo(clusterlist[cluster_ind]);
257  clusterlist[cluster_ind_phi]->addRelationTo(track2Dlist[i]);
258  }
259 
260  if (dphi_klm_tmp < m_dphi_klm_threshold * M_PI * 0.5 / 180.0 && klmtrack_ind_phi != -1) {
261  TRGGRLMATCHKLM* matklm = trackKLMmatch.appendNew();
262  matklm->set_dphi(dphi_klm_tmp);
263  matklm->addRelationTo(track2Dlist[i]);
264  matklm->addRelationTo(klmtracklist[klmtrack_ind_phi]);
265  klmtracklist[klmtrack_ind_phi]->addRelationTo(track2Dlist[i]);
266  }
267 
268  }
269 
270 
271 //do 3d track match with cluster
272  for (int i = 0; i < track3Dlist.getEntries(); i++) {
273 
274  double dr_tmp = 99999.;
275  double dz_tmp = 99999.;
276  int cluster_ind = -1;
277  for (int j = 0; j < clusterlist.getEntries(); j++) {
278  // skip the end-cap cluster
279  double _cluster_x = clusterlist[j]->getPositionX();
280  double _cluster_y = clusterlist[j]->getPositionY();
281  double _cluster_z = clusterlist[j]->getPositionZ();
282  double _cluster_theta = atan(_cluster_z / (sqrt(_cluster_x * _cluster_x + _cluster_y * _cluster_y)));
283  _cluster_theta = 0.5 * M_PI - _cluster_theta;
284  if (_cluster_theta < M_PI * 35.0 / 180.0 || _cluster_theta > M_PI * 126.0 / 180.0) continue;
285 
286  double ds_ct[2] = {99999., 99999.};
287  calculationdistance(track3Dlist[i], clusterlist[j], ds_ct, 1);
288  if (dr_tmp > ds_ct[0]) {
289  dr_tmp = ds_ct[0];
290  dz_tmp = ds_ct[1];
291  cluster_ind = j;
292  }
293  }
294  if (dr_tmp < m_dr_threshold && dz_tmp < m_dz_threshold && cluster_ind != -1) {
295  TRGGRLMATCH* mat3d = track3Dmatch.appendNew();
296  mat3d->setDeltaR(dr_tmp);
297  mat3d->setDeltaZ(dz_tmp);
298  mat3d->addRelationTo(track3Dlist[i]);
299  mat3d->addRelationTo(clusterlist[cluster_ind]);
300  // if(mat3d->getRelatedTo<CDCTriggerTrack>())std::cout<<"get match-track3D" <<std::endl;
301  // track3Dlist[i]->addRelationTo(clusterlist[cluster_ind]);
302  // if(track3Dlist[i]->getRelatedTo<TRGECLCluster>())std::cout<<"get trk-cluster" <<std::endl;
303  clusterlist[cluster_ind]->addRelationTo(track3Dlist[i]);
304  //if(clusterlist[cluster_ind]->getRelatedTo<CDCTriggerTrack>())std::cout<<"get cluster-trk" <<std::endl;
305  //if(track3Dlist[i]->getRelatedFrom<TRGECLCluster>())std::cout<<"from trk-cluster" <<std::endl;
306  }
307  }
308 
309 //pick up isolated clusters as photons with energy thrshold
310  for (int j = 0; j < clusterlist.getEntries(); j++) {
311  if (photon_cluster(clusterlist[j], track_phimap, m_e_threshold)) {
312  TRGGRLPHOTON* photon = grlphoton.appendNew();
313  photon->set_e(clusterlist[j]->getEnergyDep());
314  photon->addRelationTo(clusterlist[j]);
315  }
316  }
317 
318 // Short tracking
319  std::vector<bool> map_veto(64, 0);
320  make_veto_map(track2Dlist, map_veto);
321  short_tracking(tslist, map_veto, track_phimap_i, patterns_base0, patterns_base2, grlst, trgInfo);
322 
323 }
324 
326 {
327 }
328 
330 {
331 }
332 
333 void TRGGRLMatchModule::calculationdistance(CDCTriggerTrack* _track, TRGECLCluster* _cluster, double* ds, int _match3D)
334 {
335 
336 //double _pt = _track->getTransverseMomentum(1.5);
337  double _r = 1.0 / _track->getOmega() ;
338  double _phi = _track->getPhi0() ;
339 
340  //-- cluster/TRGECL information
341  double _cluster_x = _cluster->getPositionX();
342  double _cluster_y = _cluster->getPositionY();
343  double _cluster_z = _cluster->getPositionZ();
344  double _R = sqrt(_cluster_x * _cluster_x + _cluster_y * _cluster_y);
345 //double _D = sqrt(_cluster_x * _cluster_x + _cluster_y * _cluster_y + _cluster_z * _cluster_z);
346 //double _re_scaled_p = _pt * _D / _R;
347 
348  //-- calculation
349  if (_R > abs(2 * _r)) {
350  ds[0] = 99999.;
351  } else {
352  double theta0 = _phi - asin(_R / (2 * _r));
353 
354  double ex_x0 = _R * cos(theta0), ex_y0 = _R * sin(theta0);
355  ds[0] = sqrt((ex_x0 - _cluster_x) * (ex_x0 - _cluster_x) + (ex_y0 - _cluster_y) * (ex_y0 - _cluster_y));
356  }
357  //z information
358  if (_match3D == 1) {
359  double _z0 = _track->getZ0();
360  double _slope = _track->getCotTheta();
361  double _ex_z = _z0 + _slope * 2 * _r * asin(_R / (2 * _r));
362  ds[1] = fabs(_cluster_z - _ex_z);
363 
364  }
365 
366 }
367 
369  std::vector<bool>& phimap, std::vector<bool>& phimap_i)
370 {
371 
372  //-- 2D track information
373  double _r = 1.0 / _track->getOmega() ;
374  double _phi = _track->getPhi0() ;
375 
376  //-- 2D phi angle calculation
377  double phi_p = acos(126.0 / (2 * fabs(_r))); // adjustment angle between 0 to 0.5*M_PI
378  int charge = 0;
379  if (_r > 0) {charge = 1;}
380  else if (_r < 0) {charge = -1;}
381  else {charge = 0;}
382 
383  double phi_CDC = 0.0;
384  if (charge == 1) {
385  phi_CDC = _phi + phi_p - 0.5 * M_PI;
386  } else if (charge == -1) {
387  phi_CDC = _phi - phi_p + 0.5 * M_PI;
388  } else {
389  phi_CDC = _phi;
390  }
391 
392  if (phi_CDC > 2 * M_PI) {phi_CDC = phi_CDC - 2 * M_PI;}
393  else if (phi_CDC < 0) {phi_CDC = phi_CDC + 2 * M_PI;}
394  if (_phi > 2 * M_PI) {_phi = _phi - 2 * M_PI;}
395  else if (_phi < 0) {_phi = _phi + 2 * M_PI;}
396 
397  //-- cluster/TRGECL information
398  double _cluster_x = _cluster->getPositionX();
399  double _cluster_y = _cluster->getPositionY();
400 
401  // -- ECL phi angle
402  double phi_ECL = 0.0;
403  if (_cluster_x >= 0 && _cluster_y >= 0) {phi_ECL = atan(_cluster_y / _cluster_x);}
404  else if (_cluster_x < 0 && _cluster_y >= 0) {phi_ECL = atan(_cluster_y / _cluster_x) + M_PI;}
405  else if (_cluster_x < 0 && _cluster_y < 0) {phi_ECL = atan(_cluster_y / _cluster_x) + M_PI;}
406  else if (_cluster_x >= 0 && _cluster_y < 0) {phi_ECL = atan(_cluster_y / _cluster_x) + 2 * M_PI;}
407 
408  int phi_ECL_d = 0, phi_CDC_d = 0, phi_i_d = 0;
409  // digitization on both angle
410  for (int i = 0; i < 36; i++) {
411  if (phi_ECL > i * M_PI / 18 && phi_ECL < (i + 1)*M_PI / 18) {phi_ECL_d = i;}
412  if (_phi > i * M_PI / 18 && _phi < (i + 1)*M_PI / 18) {phi_i_d = i;}
413  if (phi_CDC > i * M_PI / 18 && phi_CDC < (i + 1)*M_PI / 18) {phi_CDC_d = i;}
414  }
415 
416  phimap[phi_CDC_d] = true;
417  phimap_i[phi_i_d] = true;
418 
419  if (abs(phi_ECL_d - phi_CDC_d) == 0 || abs(phi_ECL_d - phi_CDC_d) == 36) {dphi_d = 0;}
420  else if (abs(phi_ECL_d - phi_CDC_d) == 1 || abs(phi_ECL_d - phi_CDC_d) == 35) {dphi_d = 1;}
421  else if (abs(phi_ECL_d - phi_CDC_d) == 2 || abs(phi_ECL_d - phi_CDC_d) == 34) {dphi_d = 2;}
422  else if (abs(phi_ECL_d - phi_CDC_d) == 3 || abs(phi_ECL_d - phi_CDC_d) == 33) {dphi_d = 3;}
423  else if (abs(phi_ECL_d - phi_CDC_d) == 4 || abs(phi_ECL_d - phi_CDC_d) == 32) {dphi_d = 4;}
424  else if (abs(phi_ECL_d - phi_CDC_d) == 5 || abs(phi_ECL_d - phi_CDC_d) == 31) {dphi_d = 5;}
425  else if (abs(phi_ECL_d - phi_CDC_d) == 6 || abs(phi_ECL_d - phi_CDC_d) == 30) {dphi_d = 6;}
426  else if (abs(phi_ECL_d - phi_CDC_d) == 7 || abs(phi_ECL_d - phi_CDC_d) == 29) {dphi_d = 7;}
427  else if (abs(phi_ECL_d - phi_CDC_d) == 8 || abs(phi_ECL_d - phi_CDC_d) == 28) {dphi_d = 8;}
428  else if (abs(phi_ECL_d - phi_CDC_d) == 9 || abs(phi_ECL_d - phi_CDC_d) == 27) {dphi_d = 9;}
429  else if (abs(phi_ECL_d - phi_CDC_d) == 10 || abs(phi_ECL_d - phi_CDC_d) == 26) {dphi_d = 10;}
430  else if (abs(phi_ECL_d - phi_CDC_d) == 11 || abs(phi_ECL_d - phi_CDC_d) == 25) {dphi_d = 11;}
431  else if (abs(phi_ECL_d - phi_CDC_d) == 12 || abs(phi_ECL_d - phi_CDC_d) == 24) {dphi_d = 12;}
432  else if (abs(phi_ECL_d - phi_CDC_d) == 13 || abs(phi_ECL_d - phi_CDC_d) == 23) {dphi_d = 13;}
433  else if (abs(phi_ECL_d - phi_CDC_d) == 14 || abs(phi_ECL_d - phi_CDC_d) == 22) {dphi_d = 14;}
434  else if (abs(phi_ECL_d - phi_CDC_d) == 15 || abs(phi_ECL_d - phi_CDC_d) == 21) {dphi_d = 15;}
435  else if (abs(phi_ECL_d - phi_CDC_d) == 16 || abs(phi_ECL_d - phi_CDC_d) == 20) {dphi_d = 16;}
436  else if (abs(phi_ECL_d - phi_CDC_d) == 17 || abs(phi_ECL_d - phi_CDC_d) == 19) {dphi_d = 17;}
437  else if (abs(phi_ECL_d - phi_CDC_d) == 18) {dphi_d = 18;}
438 
439 }
440 
442 {
443 
444  //-- 2D track information
445  double _r = 1.0 / _track->getOmega() ;
446  double _phi = _track->getPhi0() ;
447 
448  //-- 2D phi angle calculation (extrapolating up to superconducting coil)
449  double phi_p = acos(176.0 / (2 * fabs(_r))); // adjustment angle between 0 to 0.5*M_PI
450  int charge = 0;
451  if (_r > 0) {charge = 1;}
452  else if (_r < 0) {charge = -1;}
453  else {charge = 0;}
454 
455  double phi_CDC = 0.0;
456  if (charge == 1) {
457  phi_CDC = _phi + phi_p - 0.5 * M_PI;
458  } else if (charge == -1) {
459  phi_CDC = _phi - phi_p + 0.5 * M_PI;
460  } else {
461  phi_CDC = _phi;
462  }
463 
464  if (phi_CDC > 2 * M_PI) {phi_CDC = phi_CDC - 2 * M_PI;}
465  else if (phi_CDC < 0) {phi_CDC = phi_CDC + 2 * M_PI;}
466 
467  // KLM track's sector central phi
468  int _sector = _klmtrack->getSector();
469  double _sector_central = 0.25 * M_PI * _sector;
470 
471  if (fabs(phi_CDC - _sector_central) < M_PI) { dphi = fabs(phi_CDC - _sector_central); }
472  else { dphi = 2 * M_PI - fabs(phi_CDC - _sector_central); }
473 
474 }
475 
476 bool TRGGRLMatchModule::photon_cluster(TRGECLCluster* _cluster, std::vector<bool> phimap, double e_threshold)
477 {
478 
479  //-- cluster/TRGECL information
480  double _cluster_x = _cluster->getPositionX();
481  double _cluster_y = _cluster->getPositionY();
482  double _cluster_z = _cluster->getPositionZ();
483  double _cluster_theta = atan(_cluster_z / (sqrt(_cluster_x * _cluster_x + _cluster_y * _cluster_y)));
484  _cluster_theta = 0.5 * M_PI - _cluster_theta;
485  bool barrel = true;
486  if (_cluster_theta < M_PI * 35.0 / 180.0 || _cluster_theta > M_PI * 126.0 / 180.0) {barrel = false;}
487  double _cluster_e = _cluster->getEnergyDep();
488 
489  // -- ECL phi angle
490  double phi_ECL = 0.0;
491  if (_cluster_x >= 0 && _cluster_y >= 0) {phi_ECL = atan(_cluster_y / _cluster_x);}
492  else if (_cluster_x < 0 && _cluster_y >= 0) {phi_ECL = atan(_cluster_y / _cluster_x) + M_PI;}
493  else if (_cluster_x < 0 && _cluster_y < 0) {phi_ECL = atan(_cluster_y / _cluster_x) + M_PI;}
494  else if (_cluster_x >= 0 && _cluster_y < 0) {phi_ECL = atan(_cluster_y / _cluster_x) + 2 * M_PI;}
495 
496  int phi_ECL_d = 0;
497  // digitization on both angle
498  for (int i = 0; i < 36; i++) {
499  if (phi_ECL > i * M_PI / 18 && phi_ECL < (i + 1)*M_PI / 18) {phi_ECL_d = i;}
500  }
501 
502  int index = phi_ECL_d, index_p = phi_ECL_d + 1, index_m = phi_ECL_d - 1;
503  if (index_p > 35) {index_p = index_p - 36;}
504  if (index_m < 0) {index_m = index_m + 36;}
505 
506  if (!phimap[index] && !phimap[index_p] && !phimap[index_m] && _cluster_e >= e_threshold && barrel) {return true;}
507  else if (!barrel) {return true;}
508  else {return false;}
509 
510 }
511 
513 {
514  if (x > 63) x -= 64;
515  if (x < 0) x += 64;
516  return x;
517 }
518 
520 {
521  if (x > 35) x -= 36;
522  if (x < 0) x += 36;
523  return x;
524 }
525 
526 void TRGGRLMatchModule::fill_pattern_base2(std::vector< std::vector<int> >& patt)
527 {
528  patt.push_back({ 0, 0, 0, 0});
529  patt.push_back({ 0, -1, 0, 0});
530  patt.push_back({ 0, -1, 1, 0});
531  patt.push_back({ 0, -1, -1, 0});
532  patt.push_back({ 0, -2, 0, 0});
533  patt.push_back({ 0, -2, 1, 0});
534  patt.push_back({ 0, -2, 2, 0});
535  patt.push_back({ 0, -2, 3, 0});
536  patt.push_back({ 0, -3, 1, 0});
537  patt.push_back({ 0, -3, 2, 0});
538  patt.push_back({ 0, -3, 3, 0});
539  patt.push_back({ 0, -4, 2, 0});
540  patt.push_back({ 0, -4, 3, 0});
541  patt.push_back({ 0, 0, 0, 1});
542  patt.push_back({ 0, 0, 1, 1});
543  patt.push_back({ 0, -1, 0, 1});
544  patt.push_back({ 0, -1, 1, 1});
545  patt.push_back({ 0, -1, 2, 1});
546  patt.push_back({ 0, -2, 2, 1});
547  patt.push_back({ 0, -2, 3, 1});
548  patt.push_back({ 0, -3, 2, 1});
549  patt.push_back({ 0, -3, 3, 1});
550  patt.push_back({ 0, 0, 0, -1});
551  patt.push_back({ 0, 0, -1, -1});
552  patt.push_back({ 0, -1, 0, -1});
553  patt.push_back({ 0, -1, -1, -1});
554  patt.push_back({ 0, -2, 0, -1});
555  patt.push_back({ 0, -2, 1, -1});
556  patt.push_back({ 0, -3, 1, -1});
557  patt.push_back({ 0, -3, 2, -1});
558  patt.push_back({ -1, -1, 0, 0});
559  patt.push_back({ -1, -1, 1, 0});
560  patt.push_back({ -1, -2, 0, 0});
561  patt.push_back({ -1, -2, 1, 0});
562  patt.push_back({ -1, -3, 1, 0});
563  patt.push_back({ -1, -3, 2, 0});
564  patt.push_back({ -1, -3, 3, 0});
565  patt.push_back({ -1, -4, 2, 0});
566  patt.push_back({ -1, -4, 3, 0});
567  patt.push_back({ 1, 0, 1, 0});
568  patt.push_back({ 1, 0, 0, 0});
569  patt.push_back({ 1, 0, -1, 0});
570  patt.push_back({ 1, -1, 0, 0});
571  patt.push_back({ 1, -1, 1, 0});
572  patt.push_back({ 1, -2, 2, 0});
573  patt.push_back({ 1, -2, 3, 0});
574  patt.push_back({ 1, -3, 2, 0});
575  patt.push_back({ 1, -3, 3, 0});
576  patt.push_back({ -1, -1, 0, 1});
577  patt.push_back({ -1, -1, 1, 1});
578  patt.push_back({ -1, -2, 0, 1});
579  patt.push_back({ -1, -2, 1, 1});
580  patt.push_back({ -1, -2, 2, 1});
581  patt.push_back({ -1, -3, 1, 1});
582  patt.push_back({ -1, -3, 2, 1});
583  patt.push_back({ -1, -3, 3, 1});
584  patt.push_back({ -1, -4, 2, 1});
585  patt.push_back({ -1, -4, 3, 1});
586  patt.push_back({ 1, 0, -1, -1});
587  patt.push_back({ 1, 0, 0, -1});
588  patt.push_back({ 1, -1, -1, -1});
589  patt.push_back({ 1, -1, 0, -1});
590  patt.push_back({ 1, -1, 1, -1});
591  patt.push_back({ 1, -2, 1, -1});
592  patt.push_back({ 1, -2, 2, -1});
593  patt.push_back({ 1, -3, 1, -1});
594  patt.push_back({ 1, -3, 2, -1});
595  patt.push_back({ -1, -1, 1, 2});
596  patt.push_back({ -1, -1, 2, 2});
597  patt.push_back({ -1, -2, 1, 2});
598  patt.push_back({ -1, -2, 2, 2});
599  patt.push_back({ -1, -2, 3, 2});
600  patt.push_back({ -1, -3, 2, 2});
601  patt.push_back({ -1, -3, 3, 2});
602  patt.push_back({ -1, -3, 4, 2});
603  patt.push_back({ 1, 0, -1, -2});
604  patt.push_back({ 1, 0, 0, -2});
605  patt.push_back({ 1, -1, 1, -2});
606  patt.push_back({ 1, -1, 0, -2});
607  patt.push_back({ 1, -1, -1, -2});
608  patt.push_back({ 1, -2, 0, -2});
609  patt.push_back({ 1, -2, 1, -2});
610  patt.push_back({ -2, -2, 0, 1});
611  patt.push_back({ -2, -2, 1, 1});
612  patt.push_back({ -2, -3, 1, 1});
613  patt.push_back({ -2, -3, 2, 1});
614  patt.push_back({ -2, -4, 2, 1});
615  patt.push_back({ -2, -4, 3, 1});
616  patt.push_back({ -2, -5, 3, 1});
617  patt.push_back({ 2, 1, 0, -1});
618  patt.push_back({ 2, 0, 1, -1});
619  patt.push_back({ 2, 0, 0, -1});
620  patt.push_back({ 2, 0, -1, -1});
621  patt.push_back({ 2, -1, 1, -1});
622  patt.push_back({ 2, -1, 0, -1});
623  patt.push_back({ 2, -2, 2, -1});
624  patt.push_back({ 2, -2, 1, -1});
625  patt.push_back({ -2, -2, 1, 2});
626  patt.push_back({ -2, -2, 2, 2});
627  patt.push_back({ -2, -3, 1, 2});
628  patt.push_back({ -2, -3, 2, 2});
629  patt.push_back({ -2, -3, 3, 2});
630  patt.push_back({ -2, -4, 2, 2});
631  patt.push_back({ -2, -4, 3, 2});
632  patt.push_back({ -2, -4, 4, 2});
633  patt.push_back({ 2, 1, 0, -2});
634  patt.push_back({ 2, 1, -1, -2});
635  patt.push_back({ 2, 0, 1, -2});
636  patt.push_back({ 2, 0, 0, -2});
637  patt.push_back({ 2, 0, -1, -2});
638  patt.push_back({ 2, 0, -2, -2});
639  patt.push_back({ 2, -1, 2, -2});
640  patt.push_back({ 2, -1, 1, -2});
641  patt.push_back({ 2, -1, 0, -2});
642  patt.push_back({ 2, -1, -1, -2});
643  patt.push_back({ 2, -2, 0, -2});
644  patt.push_back({ 2, -2, 1, -2});
645  patt.push_back({ -2, -2, 1, 3});
646  patt.push_back({ -2, -2, 2, 3});
647  patt.push_back({ -2, -3, 2, 3});
648  patt.push_back({ -2, -3, 3, 3});
649  patt.push_back({ -2, -3, 4, 3});
650  patt.push_back({ -2, -4, 3, 3});
651  patt.push_back({ -2, -4, 4, 3});
652  patt.push_back({ 2, 1, -1, -3});
653  patt.push_back({ 2, 0, -1, -3});
654  patt.push_back({ 2, 0, -2, -3});
655  patt.push_back({ 2, -1, 0, -3});
656  patt.push_back({ 2, -2, 0, -3});
657  patt.push_back({ 2, -2, 1, -3});
658  patt.push_back({ -2, -2, 2, 4});
659  patt.push_back({ -2, -3, 3, 4});
660  patt.push_back({ -2, -3, 4, 4});
661  patt.push_back({ -2, -4, 4, 4});
662  patt.push_back({ 2, -1, 0, 4});
663  patt.push_back({ 2, -1, -1, 4});
664  patt.push_back({ 2, -2, 0, 4});
665 
666 }
667 
668 void TRGGRLMatchModule::make_veto_map(StoreArray<CDCTriggerTrack> track2Dlist, std::vector<bool>& map_veto)
669 {
670  for (int i = 0; i < track2Dlist.getEntries(); i++) {
671  int _w = (int)(2271.7 * track2Dlist[i]->getOmega()) ; // omega from -33 to 33
672  if (_w >= 33) { _w = 33;}
673  else if (_w <= -33) { _w = -33;}
674  int _phi = (int)((track2Dlist[i]->getPhi0() + 2 * M_PI) / (M_PI / 32.0)); // phi_i digitized to 0 ~ 63
675 
676  int charge = 0;
677  if (_w > 0) {charge = 1;}
678  else if (_w < 0) {charge = -1;}
679  else {charge = 0;}
680 
681  _w = abs(_w);
682 
683  int L = _phi, R = _phi;
684  if (_w >= 0 && _w <= 8) { L = _phi; }
685  else if (_w >= 9 && _w <= 15) {
686  if (charge < 0) { L = _phi + 1; }
687  else { L = _phi; }
688  } else if (_w >= 16 && _w <= 24) {
689  if (charge < 0) { L = _phi + 2; }
690  else { L = _phi; }
691  } else if (_w >= 25 && _w <= 27) {
692  if (charge < 0) { L = _phi + 3; }
693  else { L = _phi; }
694  } else if (_w >= 28 && _w <= 30) {
695  if (charge < 0) { L = _phi + 3; }
696  else { L = _phi + 1; }
697  } else if (_w >= 31 && _w <= 32) {
698  if (charge < 0) { L = _phi + 4; }
699  else { L = _phi + 1; }
700  } else {
701  if (charge < 0) { L = _phi + 5; }
702  else { L = _phi + 1; }
703  }
704 
705  if (_w >= 0 && _w <= 8) { R = _phi; }
706  else if (_w >= 9 && _w <= 15) {
707  if (charge < 0) { R = _phi; }
708  else { R = _phi - 1; }
709  } else if (_w >= 16 && _w <= 24) {
710  if (charge < 0) { R = _phi; }
711  else { R = _phi - 2; }
712  } else if (_w >= 25 && _w <= 27) {
713  if (charge < 0) { R = _phi; }
714  else { R = _phi - 3; }
715  } else if (_w >= 28 && _w <= 30) {
716  if (charge < 0) { R = _phi + 1; }
717  else { R = _phi - 3; }
718  } else if (_w >= 21 && _w <= 32) {
719  if (charge < 0) { R = _phi + 1; }
720  else { R = _phi - 4; }
721  } else {
722  if (charge < 0) { R = _phi + 1; }
723  else { R = _phi - 5; }
724  }
725 
726  // L should be > R
727  for (int j = R - 1; j < L + 2; j++) {
728  map_veto[N64(j)] = true;
729  }
730  }
731 
732 }
733 
735  std::vector<bool> phimap_i,
736  std::vector< std::vector<int> >& pattern_base0, std::vector< std::vector<int> >& pattern_base2,
738  StoreObjPtr<TRGGRLInfo> trgInfo)
739 {
740  std::vector<bool> SL0(64, 0);
741  std::vector<bool> SL1(64, 0);
742  std::vector<bool> SL2(64, 0);
743  std::vector<bool> SL3(64, 0);
744  std::vector<bool> SL4(64, 0);
745  std::vector<bool> ST0(64, 0);
746  std::vector<bool> ST0_36b(36, 0);
747  std::vector<bool> ST2(64, 0);
748  std::vector<int> patt_ID(64, -1);
749 
750 // std::vector<bool> st_ec1(64, 0);
751 // std::vector<bool> st_ec1_36b(36, 0);
752 // std::vector<bool> st_ec2(64, 0);
753 // std::vector<bool> st_ec2_36b(36, 0);
754 
755 //-- collecting TSF info in SL0~4
756  for (int i = 0; i < tslist.getEntries(); i++) {
757  int id = tslist[i]->getSegmentID();
758  int sl = 0;
759  if (id >= 0 * 32 && id < 5 * 32) {sl = 0; id -= 0;}
760  else if (id >= 5 * 32 && id < 10 * 32) {sl = 1; id -= 5 * 32;}
761  else if (id >= 10 * 32 && id < 16 * 32) {sl = 2; id -= 10 * 32;}
762  else if (id >= 16 * 32 && id < 23 * 32) {sl = 3; id -= 16 * 32;}
763  else if (id >= 23 * 32 && id < 31 * 32) {sl = 4; id -= 23 * 32;}
764  else continue;
765 
766  if (sl == 0) {
767  int X = (int)(id / 5), Y = id % 5;
768  if (Y == 0 || Y == 1) { SL0[2 * X] = true; }
769  else if (Y == 3 || Y == 4) { SL0[2 * X + 1] = true; }
770  else { SL0[2 * X] = true; SL0[2 * X + 1] = true; }
771  } else if (sl == 1) {
772  int X = (int)(id / 5), Y = id % 5;
773  if (Y == 0 || Y == 1) { SL1[2 * X] = true; }
774  else if (Y == 3 || Y == 4) { SL1[2 * X + 1] = true; }
775  else { SL1[2 * X] = true; SL1[2 * X + 1] = true; }
776  } else if (sl == 2) {
777  int X = (int)(id / 3);
778  SL2[X] = true;
779  } else if (sl == 3) {
780  int X = (int)(id / 7), Y = id % 7;
781  if (Y == 0 || Y == 1 || Y == 2) { SL3[2 * X] = true; }
782  else if (Y == 4 || Y == 5 || Y == 6) { SL3[2 * X + 1] = true; }
783  else { SL3[2 * X] = true; SL3[2 * X + 1] = true; }
784  } else if (sl == 4) {
785  int X = (int)(id / 4);
786  SL4[X] = true;
787  }
788 
789  }
790 
791 //-- making veto
792  for (int i = 0; i < 64; i++) {
793  if (map_veto[i]) {SL0[i] = false; SL1[i] = false; SL2[i] = false;}
794  }
795  /*
796  for (int i = 0; i < 64; i++) { std::cout<<map_veto[63-i]; if((64-i)%10==1) std::cout<<" ";}
797  std::cout<<std::endl;
798  for (int i = 0; i < 64; i++) { std::cout<<SL4[63-i]; if((64-i)%10==1) std::cout<<" ";}
799  std::cout<<std::endl;
800  for (int i = 0; i < 64; i++) { std::cout<<SL3[63-i]; if((64-i)%10==1) std::cout<<" ";}
801  std::cout<<std::endl;
802  for (int i = 0; i < 64; i++) { std::cout<<SL2[63-i]; if((64-i)%10==1) std::cout<<" ";}
803  std::cout<<std::endl;
804  for (int i = 0; i < 64; i++) { std::cout<<SL1[63-i]; if((64-i)%10==1) std::cout<<" ";}
805  std::cout<<std::endl;
806  for (int i = 0; i < 64; i++) { std::cout<<SL0[63-i]; if((64-i)%10==1) std::cout<<" ";}
807  std::cout<<std::endl;
808  */
809 //-- doing short tracking
810 
811  std::vector< std::vector<int> > stlist_buf(0);
812 
813  // -- ST finding with SL2
814  for (int i = 0; i < 64; i++) {
815 
816  int ID0 = 0;
817  int ID1 = 0;
818  int ID2 = 0;
819  int ID3 = 0;
820  int ID4 = 0;
821  stlist_buf.push_back({0, 0, 0, 0, 0, 0});
822 
823  if (!SL2[i]) continue;
824  bool SL2_already_found = false;
825 
826  for (int p = 0; p < 137; p++) {
827 
828  // following patterns will not be used.
829  if (p == 4) continue;
830  if (p == 5) continue;
831  if (p == 17) continue;
832  if (p == 26) continue;
833  if (p == 38) continue;
834  if (p == 41) continue;
835  if (p == 42) continue;
836  if (p == 47) continue;
837  if (p == 50) continue;
838  if (p == 60) continue;
839  if (p == 63) continue;
840  if (p == 64) continue;
841  if (p == 74) continue;
842  if (p == 93) continue;
843  if (p == 94) continue;
844  if (p == 95) continue;
845  if (p == 96) continue;
846  if (p == 104) continue;
847  if (p == 113) continue;
848  if (p == 114) continue;
849  if (p == 115) continue;
850  if (p == 123) continue;
851  if (p == 134) continue;
852  if (p == 135) continue;
853  if (p == 136) continue;
854 
855  int x0 = pattern_base2[p][0];
856  int x1 = pattern_base2[p][1];
857  int x3 = pattern_base2[p][2];
858  int x4 = pattern_base2[p][3];
859 
860 
861  if (SL2[i] && SL0[N64(i + x0)] && SL1[N64(i + x1)] && SL3[N64(i + x3)] && SL4[N64(i + x4)] && !SL2_already_found) {
862  ST2[i] = true;
863  ID0 = N64(i + x0);
864  ID1 = N64(i + x1);
865  ID2 = i;
866  ID3 = N64(i + x3);
867  ID4 = N64(i + x4);
868  SL2_already_found = true; // if it has been found in previous pattern, no need to do it again.
869  }
870 
871  // if a pattern is found, no need to look for other pattern
872  if (SL2_already_found) break;
873 
874  }
875 
876  if (SL2_already_found) {
877  stlist_buf[i][0] = 1;
878  stlist_buf[i][1] = ID0;
879  stlist_buf[i][2] = ID1;
880  stlist_buf[i][3] = ID2;
881  stlist_buf[i][4] = ID3;
882  stlist_buf[i][5] = ID4;
883  }
884  }
886 //-- ST finding with SL0
887  for (int i = 0; i < 64; i++) {
888 
889  if (!SL0[i]) continue;
890  bool SL0_already_found = false;
891 
892  for (int p = 0; p < 137; p++) {
893 
894  // following patterns will not be used.
895  if (p == 4) continue;
896  if (p == 5) continue;
897  if (p == 17) continue;
898  if (p == 26) continue;
899  if (p == 38) continue;
900  if (p == 41) continue;
901  if (p == 42) continue;
902  if (p == 47) continue;
903  if (p == 50) continue;
904  if (p == 60) continue;
905  if (p == 63) continue;
906  if (p == 64) continue;
907  if (p == 74) continue;
908  if (p == 93) continue;
909  if (p == 94) continue;
910  if (p == 95) continue;
911  if (p == 96) continue;
912  if (p == 104) continue;
913  if (p == 113) continue;
914  if (p == 114) continue;
915  if (p == 115) continue;
916  if (p == 123) continue;
917  if (p == 134) continue;
918  if (p == 135) continue;
919  if (p == 136) continue;
920 
921  int y1 = pattern_base0[p][0];
922  int y2 = pattern_base0[p][1];
923  int y3 = pattern_base0[p][2];
924  int y4 = pattern_base0[p][3];
925 
926  if (SL0[i] && SL1[N64(i + y1)] && SL2[N64(i + y2)] && SL3[N64(i + y3)] && SL4[N64(i + y4)] && !SL0_already_found) {
927  ST0[i] = true;
928  if (patt_ID[i] < 0) { patt_ID[i] = p; }
929  SL0_already_found = true; // if it has been found in previous pattern, no need to do it again.
930  }
931 
932  // if a pattern is found, no need to look for other pattern
933  if (SL0_already_found) break;
934 
935  }
936 
937  }
939 //-- extrapolation
940  /*
941  for (int i = 0; i < 64; i++) {
942  if (patt_ID[i] == -1) continue;
943 
944  int ec = 0, l = 0, r = 0;
945  extrapolation(patt_ID[i], l, r, ec);
946  if (ec == 1) {
947  for (int e = l; e <= r; e++) { st_ec1[N64(i + e)] = true; }
948  }
949  if (ec == 2) {
950  for (int e = l; e <= r; e++) { st_ec2[N64(i + e)] = true; }
951  }
952 
953  }
954  */
955 //-- 64b into 36b
956  for (int i = 0; i < 4; i++) {
957  ST0_36b[0 + 9 * i] = ST0[0 + 16 * i] or ST0[1 + 16 * i];
958  ST0_36b[1 + 9 * i] = ST0[1 + 16 * i] or ST0[2 + 16 * i] or ST0[3 + 16 * i];
959  ST0_36b[2 + 9 * i] = ST0[3 + 16 * i] or ST0[4 + 16 * i] or ST0[5 + 16 * i];
960  ST0_36b[3 + 9 * i] = ST0[5 + 16 * i] or ST0[6 + 16 * i] or ST0[7 + 16 * i];
961  ST0_36b[4 + 9 * i] = ST0[7 + 16 * i] or ST0[8 + 16 * i];
962  ST0_36b[5 + 9 * i] = ST0[8 + 16 * i] or ST0[9 + 16 * i] or ST0[10 + 16 * i];
963  ST0_36b[6 + 9 * i] = ST0[10 + 16 * i] or ST0[11 + 16 * i] or ST0[12 + 16 * i];
964  ST0_36b[7 + 9 * i] = ST0[12 + 16 * i] or ST0[13 + 16 * i] or ST0[14 + 16 * i];
965  ST0_36b[8 + 9 * i] = ST0[14 + 16 * i] or ST0[15 + 16 * i];
966  /*
967  st_ec1_36b[0 + 9 * i] = st_ec1[0 + 16 * i] or st_ec1[1 + 16 * i];
968  st_ec1_36b[1 + 9 * i] = st_ec1[1 + 16 * i] or st_ec1[2 + 16 * i] or st_ec1[3 + 16 * i];
969  st_ec1_36b[2 + 9 * i] = st_ec1[3 + 16 * i] or st_ec1[4 + 16 * i] or st_ec1[5 + 16 * i];
970  st_ec1_36b[3 + 9 * i] = st_ec1[5 + 16 * i] or st_ec1[6 + 16 * i] or st_ec1[7 + 16 * i];
971  st_ec1_36b[4 + 9 * i] = st_ec1[7 + 16 * i] or st_ec1[8 + 16 * i];
972  st_ec1_36b[5 + 9 * i] = st_ec1[8 + 16 * i] or st_ec1[9 + 16 * i] or st_ec1[10 + 16 * i];
973  st_ec1_36b[6 + 9 * i] = st_ec1[10 + 16 * i] or st_ec1[11 + 16 * i] or st_ec1[12 + 16 * i];
974  st_ec1_36b[7 + 9 * i] = st_ec1[12 + 16 * i] or st_ec1[13 + 16 * i] or st_ec1[14 + 16 * i];
975  st_ec1_36b[8 + 9 * i] = st_ec1[14 + 16 * i] or st_ec1[15 + 16 * i];
976 
977  st_ec2_36b[0 + 9 * i] = st_ec2[0 + 16 * i] or st_ec2[1 + 16 * i];
978  st_ec2_36b[1 + 9 * i] = st_ec2[1 + 16 * i] or st_ec2[2 + 16 * i] or st_ec2[3 + 16 * i];
979  st_ec2_36b[2 + 9 * i] = st_ec2[3 + 16 * i] or st_ec2[4 + 16 * i] or st_ec2[5 + 16 * i];
980  st_ec2_36b[3 + 9 * i] = st_ec2[5 + 16 * i] or st_ec2[6 + 16 * i] or st_ec2[7 + 16 * i];
981  st_ec2_36b[4 + 9 * i] = st_ec2[7 + 16 * i] or st_ec2[8 + 16 * i];
982  st_ec2_36b[5 + 9 * i] = st_ec2[8 + 16 * i] or st_ec2[9 + 16 * i] or st_ec2[10 + 16 * i];
983  st_ec2_36b[6 + 9 * i] = st_ec2[10 + 16 * i] or st_ec2[11 + 16 * i] or st_ec2[12 + 16 * i];
984  st_ec2_36b[7 + 9 * i] = st_ec2[12 + 16 * i] or st_ec2[13 + 16 * i] or st_ec2[14 + 16 * i];
985  st_ec2_36b[8 + 9 * i] = st_ec2[14 + 16 * i] or st_ec2[15 + 16 * i];
986  */
987  }
988 
989 
990 //-- Summary info
991 
992  int N_ST = 0;
993  bool s2s3 = false;
994  bool s2s5 = false;
995  bool s2so = false;
996  bool s2f3 = false;
997  bool s2f5 = false;
998  bool s2fo = false;
999 
1000 //-- short track counting on ST2
1001  for (int i = 0; i < 64; i++) {
1002  if (ST2[i]) {
1003  N_ST++;
1004  ST2[i] = false;
1005  int L = i - 1, R = i + 1;
1006  while (ST2[N64(L)]) {
1007  ST2[N64(L)] = false;
1008  L--;
1009  }
1010  while (ST2[N64(R)]) {
1011  ST2[N64(R)] = false;
1012  R++;
1013  }
1014 
1015  //-- Fill the store array
1016  L++; R--;
1017  int index = N64((L + R) / 2); // fill the middle one when multiple ST is found continuously in the map
1018  TRGGRLShortTrack* st = grlst.appendNew();
1019  st->set_TS_ID(0, stlist_buf[index][1]);
1020  st->set_TS_ID(1, stlist_buf[index][2]);
1021  st->set_TS_ID(2, stlist_buf[index][3]);
1022  st->set_TS_ID(3, stlist_buf[index][4]);
1023  st->set_TS_ID(4, stlist_buf[index][5]);
1024  }
1025  }
1026 
1027 //-- b2b info with ST0 and phi_i map
1028  for (int i = 0; i < 36; i++) {
1029  s2s3 = (ST0_36b[i] and (ST0_36b[N36(i + 18)] or ST0_36b[N36(i + 17)] or ST0_36b[N36(i + 19)])) or s2s3;
1030  s2s5 = (ST0_36b[i] and (ST0_36b[N36(i + 18)] or ST0_36b[N36(i + 17)] or ST0_36b[N36(i + 19)]
1031  or ST0_36b[N36(i + 16)] or ST0_36b[N36(i + 20)])) or s2s5;
1032  s2so = (ST0_36b[i] and (ST0_36b[N36(i + 18)] or ST0_36b[N36(i + 17)] or ST0_36b[N36(i + 19)]
1033  or ST0_36b[N36(i + 16)] or ST0_36b[N36(i + 20)]
1034  or ST0_36b[N36(i + 15)] or ST0_36b[N36(i + 21)]
1035  or ST0_36b[N36(i + 14)] or ST0_36b[N36(i + 22)]
1036  or ST0_36b[N36(i + 13)] or ST0_36b[N36(i + 23)]
1037  or ST0_36b[N36(i + 12)] or ST0_36b[N36(i + 24)]
1038  or ST0_36b[N36(i + 11)] or ST0_36b[N36(i + 25)]
1039  or ST0_36b[N36(i + 10)] or ST0_36b[N36(i + 26)]
1040  or ST0_36b[N36(i + 9)] or ST0_36b[N36(i + 27)])) or s2so ;
1041 
1042  s2f3 = (phimap_i[i] and (ST0_36b[N36(i + 18)] or ST0_36b[N36(i + 17)] or ST0_36b[N36(i + 19)])) or s2f3;
1043  s2f5 = (phimap_i[i] and (ST0_36b[N36(i + 18)] or ST0_36b[N36(i + 17)] or ST0_36b[N36(i + 19)]
1044  or ST0_36b[N36(i + 16)] or ST0_36b[N36(i + 20)])) or s2f5;
1045  s2fo = (phimap_i[i] and (ST0_36b[N36(i + 18)] or ST0_36b[N36(i + 17)] or ST0_36b[N36(i + 19)]
1046  or ST0_36b[N36(i + 16)] or ST0_36b[N36(i + 20)]
1047  or ST0_36b[N36(i + 15)] or ST0_36b[N36(i + 21)]
1048  or ST0_36b[N36(i + 14)] or ST0_36b[N36(i + 22)]
1049  or ST0_36b[N36(i + 13)] or ST0_36b[N36(i + 23)]
1050  or ST0_36b[N36(i + 12)] or ST0_36b[N36(i + 24)]
1051  or ST0_36b[N36(i + 11)] or ST0_36b[N36(i + 25)]
1052  or ST0_36b[N36(i + 10)] or ST0_36b[N36(i + 26)]
1053  or ST0_36b[N36(i + 9)] or ST0_36b[N36(i + 27)])) or s2fo ;
1054  }
1055 
1056 //-- set results
1057  trgInfo->setNshorttrk(N_ST);
1058  trgInfo->sets2s3(s2s3);
1059  trgInfo->sets2s5(s2s5);
1060  trgInfo->sets2so(s2so);
1061  trgInfo->sets2f3(s2f3);
1062  trgInfo->sets2f5(s2f5);
1063  trgInfo->sets2fo(s2fo);
1064  trgInfo->setbwdsb(0);
1065  trgInfo->setbwdnb(0);
1066  trgInfo->setfwdsb(0);
1067  trgInfo->setfwdnb(0);
1068  trgInfo->setbrlfb(0);
1069  trgInfo->setbrlnb(0);
1070 
1071 }
1072 
1073 
1074 void
1075 TRGGRLMatchModule::extrapolation(int pattern, int& l, int& r, int& ec)
1076 {
1077  if (pattern == 6) {ec = 1; l = 0; r = 1;}
1078  if (pattern == 7) {ec = 1; l = 0; r = 1;}
1079  if (pattern == 8) {ec = 1; l = 0; r = 1;}
1080  if (pattern == 9) {ec = 1; l = 0; r = 2;}
1081  if (pattern == 10) {ec = 1; l = 0; r = 2;}
1082  if (pattern == 11) {ec = 1; l = 0; r = 1;}
1083  if (pattern == 12) {ec = 1; l = 0; r = 2;}
1084  if (pattern == 18) {ec = 1; l = 0; r = 2;}
1085  if (pattern == 19) {ec = 1; l = 0; r = 2;}
1086  if (pattern == 20) {ec = 1; l = 0; r = 4;}
1087  if (pattern == 21) {ec = 1; l = 0; r = 4;}
1088  if (pattern == 28) {ec = 1; l = -4; r = 0;}
1089  if (pattern == 29) {ec = 1; l = -3; r = 0;}
1090  if (pattern == 34) {ec = 1; l = 0; r = 1;}
1091  if (pattern == 35) {ec = 1; l = 0; r = 3;}
1092  if (pattern == 36) {ec = 1; l = 1; r = 3;}
1093  if (pattern == 37) {ec = 1; l = 0; r = 3;}
1094  if (pattern == 44) {ec = 1; l = -4; r = 0;}
1095  if (pattern == 45) {ec = 1; l = -2; r = 0;}
1096  if (pattern == 46) {ec = 1; l = -3; r = 0;}
1097  if (pattern == 54) {ec = 1; l = 1; r = 7;}
1098  if (pattern == 55) {ec = 1; l = 1; r = 6;}
1099  if (pattern == 56) {ec = 1; l = 1; r = 5;}
1100  if (pattern == 57) {ec = 1; l = 1; r = 5;}
1101  if (pattern == 64) {ec = 1; l = -6; r = -1;}
1102  if (pattern == 73) {ec = 1; l = 3; r = 13;}
1103  if (pattern == 81) {ec = 1; l = -10; r = -3;}
1104  if (pattern == 86) {ec = 1; l = 3; r = 12;}
1105  if (pattern == 87) {ec = 1; l = 3; r = 6;}
1106  if (pattern == 100) {ec = 1; l = 7; r = 20;}
1107  if (pattern == 101) {ec = 1; l = 5; r = 20;}
1108  if (pattern == 102) {ec = 1; l = 5; r = 20;}
1109  if (pattern == 103) {ec = 1; l = 4; r = 14;}
1110  if (pattern == 111) {ec = 1; l = -12; r = -5;}
1111  if (pattern == 112) {ec = 1; l = -18; r = -5;}
1112  if (pattern == 116) {ec = 1; l = -11; r = -6;}
1113  if (pattern == 120) {ec = 1; l = 7; r = 21;}
1114  if (pattern == 121) {ec = 1; l = 7; r = 14;}
1115  if (pattern == 122) {ec = 1; l = 7; r = 21;}
1116  if (pattern == 127) {ec = 1; l = -21; r = -8;}
1117  if (pattern == 128) {ec = 1; l = -15; r = -7;}
1118  if (pattern == 129) {ec = 1; l = -12; r = -7;}
1119  if (pattern == 132) {ec = 1; l = 10; r = 18;}
1120  if (pattern == 133) {ec = 1; l = 8; r = 18;}
1121 
1122  if (pattern == 0) {ec = 2; l = -3; r = 1;}
1123  if (pattern == 1) {ec = 2; l = -3; r = 1;}
1124  if (pattern == 3) {ec = 2; l = -3; r = 0;}
1125  if (pattern == 13) {ec = 2; l = 0; r = 3;}
1126  if (pattern == 14) {ec = 2; l = 0; r = 4;}
1127  if (pattern == 15) {ec = 2; l = 0; r = 5;}
1128  if (pattern == 22) {ec = 2; l = -4; r = -1;}
1129  if (pattern == 23) {ec = 2; l = -5; r = -1;}
1130  if (pattern == 24) {ec = 2; l = -3; r = 0;}
1131  if (pattern == 25) {ec = 2; l = -4; r = 0;}
1132  if (pattern == 30) {ec = 2; l = 1; r = 5;}
1133  if (pattern == 39) {ec = 2; l = -2; r = 0;}
1134  if (pattern == 40) {ec = 2; l = -2; r = 0;}
1135  if (pattern == 48) {ec = 2; l = 2; r = 6;}
1136  if (pattern == 49) {ec = 2; l = 3; r = 8;}
1137  if (pattern == 58) {ec = 2; l = -9; r = -3;}
1138  if (pattern == 59) {ec = 2; l = -9; r = -3;}
1139  if (pattern == 67) {ec = 2; l = 5; r = 11;}
1140  if (pattern == 75) {ec = 2; l = -13; r = -6;}
1141  if (pattern == 82) {ec = 2; l = 5; r = 9;}
1142  if (pattern == 83) {ec = 2; l = 5; r = 9;}
1143  if (pattern == 89) {ec = 2; l = -10; r = -4;}
1144  if (pattern == 92) {ec = 2; l = -10; r = -4;}
1145  if (pattern == 97) {ec = 2; l = 7; r = 19;}
1146  if (pattern == 105) {ec = 2; l = -16; r = -10;}
1147  if (pattern == 106) {ec = 2; l = -17; r = -7;}
1148  if (pattern == 109) {ec = 2; l = -17; r = -6;}
1149  if (pattern == 111) {ec = 2; l = -16; r = -7;}
1150  if (pattern == 117) {ec = 2; l = 9; r = 19;}
1151  if (pattern == 118) {ec = 2; l = 9; r = 19;}
1152  if (pattern == 124) {ec = 2; l = -17; r = -8;}
1153  if (pattern == 125) {ec = 2; l = -17; r = -8;}
1154  if (pattern == 126) {ec = 2; l = -17; r = -8;}
1155 
1156 }
1157 
1158 
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::TRGGRLMatchModule::m_phimatch_tracklist
std::string m_phimatch_tracklist
the matched 2d track list by phi matching
Definition: TRGGRLMatchModule.h:152
Belle2::TRGGRLMatchModule::fill_pattern_base2
void fill_pattern_base2(std::vector< std::vector< int > > &patt)
Fill the patterns in short tracking logic.
Definition: TRGGRLMatchModule.cc:526
Belle2::TRGGRLMatchModule::photon_cluster
bool photon_cluster(TRGECLCluster *cluster, std::vector< bool > track_phimap, double e_threshold)
determine photon from isolated cluster
Definition: TRGGRLMatchModule.cc:476
Belle2::TRGGRLPHOTON::set_e
void set_e(double e)
set energy
Definition: TRGGRLPHOTON.h:44
Belle2::KLMTriggerTrack
Store KLM TRG track information as a ROOT object.
Definition: KLMTriggerTrack.h:30
Belle2::TRGGRLMatchModule::m_clusterlist
std::string m_clusterlist
the ecl cluster list
Definition: TRGGRLMatchModule.h:142
Belle2::StoreArray::registerRelationTo
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:150
Belle2::TRGGRLMatchModule::m_klmmatch_tracklist
std::string m_klmmatch_tracklist
the matched 2d track list by KLM matching
Definition: TRGGRLMatchModule.h:156
Belle2::TRGGRLMatchModule::m_dphi_d_threshold
int m_dphi_d_threshold
max value of dphi_d to be identified as match, 1 digit = 10 degrees
Definition: TRGGRLMatchModule.h:128
Belle2::TRGECLCluster
Example Detector.
Definition: TRGECLCluster.h:25
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::TRGGRLPHOTON
a class for neutral ECL cluster in TRGGRL
Definition: TRGGRLPHOTON.h:26
Belle2::KLMTriggerTrack::getSector
int getSector() const
Get sector number.
Definition: KLMTriggerTrack.h:85
Belle2::TRGGRLMATCHKLM
a class for CDC2D-KLM Matching in TRGGRL
Definition: TRGGRLMATCHKLM.h:26
Belle2::TRGGRLMatchModule::calculationphiangle
void calculationphiangle(CDCTriggerTrack *track, TRGECLCluster *cluster, int &dphi_d, std::vector< bool > &track_phimap, std::vector< bool > &track_phimap_i)
calculate dphi_d between track and cluster
Definition: TRGGRLMatchModule.cc:368
Belle2::TRGGRLMatchModule::m_TRGGRLInfo
StoreObjPtr< TRGGRLInfo > m_TRGGRLInfo
output for TRGGRLInfo
Definition: TRGGRLMatchModule.h:110
Belle2::TRGGRLMatchModule
Match between CDC trigger track and ECL trigger cluster.
Definition: TRGGRLMatchModule.h:48
Belle2::TRGGRLMatchModule::N64
int N64(int x)
Force an int to be witnin 0 to 63.
Definition: TRGGRLMatchModule.cc:512
Belle2::CDCTriggerTrack
Track created by the CDC trigger.
Definition: CDCTriggerTrack.h:15
Belle2::TRGGRLMatchModule::N36
int N36(int x)
Force an int to be witnin 0 to 35.
Definition: TRGGRLMatchModule.cc:519
Belle2::TRGECLCluster::getPositionX
double getPositionX() const
The method to get hit average time Get Energy weighted position X.
Definition: TRGECLCluster.h:127
Belle2::TRGGRLMatchModule::m_dz_threshold
double m_dz_threshold
max value of dz to be identified as match
Definition: TRGGRLMatchModule.h:126
Belle2::RelationsInterface::addRelationTo
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
Definition: RelationsObject.h:144
Belle2::TRGGRLMatchModule::short_tracking
void short_tracking(StoreArray< CDCTriggerSegmentHit > tslist, std::vector< bool > map_veto, std::vector< bool > phimap_i, std::vector< std::vector< int > > &pattern_base0, std::vector< std::vector< int > > &pattern_base2, StoreArray< TRGGRLShortTrack > grlst, StoreObjPtr< TRGGRLInfo > trgInfo)
Short tracking logic.
Definition: TRGGRLMatchModule.cc:734
Belle2::TRGGRLMatchModule::m_grlphotonlist
std::string m_grlphotonlist
Non-matched cluster list at GRL.
Definition: TRGGRLMatchModule.h:158
Belle2::TRGGRLMATCH::set_e
void set_e(double e)
set the cluster energy
Definition: TRGGRLMATCH.h:57
Belle2::TRGGRLMatchModule::m_hitCollectionName
std::string m_hitCollectionName
Track Segment list.
Definition: TRGGRLMatchModule.h:160
Belle2::TRGGRLMatchModule::m_dr_threshold
double m_dr_threshold
max value of dr to be identified as match
Definition: TRGGRLMatchModule.h:124
Belle2::TRGGRLMatchModule::m_dphi_klm_threshold
double m_dphi_klm_threshold
max value of dphi (CDC track to KLM sector) to be identified as match (in degrees)
Definition: TRGGRLMatchModule.h:132
Belle2::TRGGRLMatchModule::terminate
virtual void terminate() override
Termination action.
Definition: TRGGRLMatchModule.cc:329
Belle2::TRGECLCluster::getPositionZ
double getPositionZ() const
Get Energy weighted position Z.
Definition: TRGECLCluster.h:131
Belle2::TRGECLCluster::getPositionY
double getPositionY() const
Get Energy weighted position Y.
Definition: TRGECLCluster.h:129
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::TRGGRLMatchModule::m_3d_tracklist
std::string m_3d_tracklist
the 3D NN track list
Definition: TRGGRLMatchModule.h:140
Belle2::TRGGRLMatchModule::initialize
virtual void initialize() override
Initialize the parameters.
Definition: TRGGRLMatchModule.cc:83
Belle2::TRGGRLShortTrack
a class for neutral ECL cluster in TRGGRL
Definition: TRGGRLShortTrack.h:27
Belle2::TRGGRLMatchModule::m_TrgGrlInformationName
std::string m_TrgGrlInformationName
Name of the StoreArray holding projects information from grl.
Definition: TRGGRLMatchModule.h:164
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::TRGGRLMatchModule::calculationdistance
void calculationdistance(CDCTriggerTrack *track, TRGECLCluster *cluster, double *ds, int _match3D)
calculate dr and dz between track and cluster
Definition: TRGGRLMatchModule.cc:333
Belle2::TRGGRLMatchModule::extrapolation
void extrapolation(int pattern, int &l, int &r, int &ec)
Short track extrapolation (to endcap) function.
Definition: TRGGRLMatchModule.cc:1075
Belle2::TRGGRLMatchModule::m_klmtracklist
std::string m_klmtracklist
the KLM track list
Definition: TRGGRLMatchModule.h:144
Belle2::TRGGRLMatchModule::endRun
virtual void endRun() override
End-of-run action.
Definition: TRGGRLMatchModule.cc:325
Belle2::TRGGRLMatchModule::m_e_threshold
double m_e_threshold
min value of isolated cluster energy
Definition: TRGGRLMatchModule.h:130
Belle2::TRGGRLMATCH::set_dphi_d
void set_dphi_d(double dphi_d)
set the dphi_d
Definition: TRGGRLMATCH.h:54
Belle2::TRGGRLMatchModule::m_2dmatch_tracklist
std::string m_2dmatch_tracklist
the distance in phi direction between track and cluster
Definition: TRGGRLMatchModule.h:150
Belle2::TRGGRLMatchModule::m_grlstCollectionName
std::string m_grlstCollectionName
GRL short track list.
Definition: TRGGRLMatchModule.h:162
Belle2::TRGGRLMatchModule::patterns_base0
std::vector< std::vector< int > > patterns_base0
Short tracking patterns based on SL0.
Definition: TRGGRLMatchModule.h:166
Belle2::TRGGRLMatchModule::track_phimap
std::vector< bool > track_phimap
36 bits phi map of all 2D tracks
Definition: TRGGRLMatchModule.h:134
Belle2::TRGGRLMatchModule::patterns_base2
std::vector< std::vector< int > > patterns_base2
Short tracking patterns based on SL2.
Definition: TRGGRLMatchModule.h:168
Belle2::TRGGRLMatchModule::track_phimap_i
std::vector< bool > track_phimap_i
36 bits phi map of all 2D tracks
Definition: TRGGRLMatchModule.h:136
Belle2::TRGGRLMatchModule::m_2d_tracklist
std::string m_2d_tracklist
the 2D finder track list
Definition: TRGGRLMatchModule.h:138
Belle2::TRGGRLMATCH::setDeltaZ
void setDeltaZ(double deltaz)
set the Delta Z
Definition: TRGGRLMATCH.h:51
Belle2::TRGGRLMatchModule::sectormatching_klm
void sectormatching_klm(CDCTriggerTrack *track, KLMTriggerTrack *klmtrack, double &dphi)
calculate dphi between 2D track and KLM track
Definition: TRGGRLMatchModule.cc:441
Belle2::TRGGRLMatchModule::make_veto_map
void make_veto_map(StoreArray< CDCTriggerTrack > track2Dlist, std::vector< bool > &map_veto)
Make the full track phi veto map for short tracking.
Definition: TRGGRLMatchModule.cc:668
Belle2::TRGGRLMatch
A class to represent a matching candidate in TRGGRL A matching candidate consists of a TRGCDCTrack an...
Definition: TRGGRLMatch.h:27
Belle2::TRGGRLMatchModule::m_3dmatch_tracklist
std::string m_3dmatch_tracklist
the matched 3d track list
Definition: TRGGRLMatchModule.h:154
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::TRGGRLMatchModule::event
virtual void event() override
Event processor.
Definition: TRGGRLMatchModule.cc:170
Belle2::TRGGRLMATCH::setDeltaR
void setDeltaR(double deltar)
set the Delta R
Definition: TRGGRLMATCH.h:48
Belle2::TRGECLCluster::getEnergyDep
double getEnergyDep() const
The method to get deposited energy.
Definition: TRGECLCluster.h:120
Belle2::TRGGRLMatchModule::beginRun
virtual void beginRun() override
Called when entering a new run.
Definition: TRGGRLMatchModule.cc:166
Belle2::TRGGRLMATCH
a class for CDC2D-ECL Matching in TRGGRL
Definition: TRGGRLMATCH.h:26
Belle2::TRGGRLMATCHKLM::set_dphi
void set_dphi(double dphi)
set the dphi
Definition: TRGGRLMATCHKLM.h:38
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::TRGGRLMatchModule::~TRGGRLMatchModule
virtual ~TRGGRLMatchModule()
Destructor.
Definition: TRGGRLMatchModule.cc:79