15 from ROOT
import Belle2
16 import modularAnalysis
as ma
17 from array
import array
20 from heapq
import nlargest
21 from effCalculation
import EffCalculation
30 sys.exit(
"ztsim02.py> # of arg is strange. Exit.")
33 def Count_mcpart(parts):
38 pdglist = [11, 13, 2212, 321, 211]
39 pdg = math.fabs(part.getPDG())
40 charge = part.getCharge()
41 theta = part.getMomentum().Theta() * Fac
42 if theta > dec_cdc[0]
and theta < dec_cdc[1]:
46 if math.fabs(pdg) == 22:
47 if theta > dec_ecl[0]
and theta < dec_ecl[1]:
52 def Count_part_inDetector(parts):
56 charge = part.getCharge()
59 if math.fabs(charge) > 0:
61 elif math.fabs(pdg) == 22:
66 def Theta_MCTruth(part):
67 theta = part.getMomentum().Theta() * Fac
71 def NeutralCluster(clusters):
73 for cluster
in clusters:
74 if cluster.getRelatedFrom(
'TRGNNTracks'):
77 neutral_cluster.append(cluster)
78 return neutral_cluster
81 def Vec_Cluster(cluster, CMS=False):
82 x = cluster.getPositionX()
83 y = cluster.getPositionY()
84 z = cluster.getPositionZ()
85 e = cluster.getEnergyDep()
86 vec = ROOT.TVector3(x, y, z)
89 v_mom = ROOT.TLorentzVector(vec_unit, e)
91 new_vec = trans.rotateLabToCms() * v_mom
94 new_theta = new_vec.Theta() * Fac
95 new_phi = new_vec.Phi() * Fac
99 NVC = [new_e, new_theta, new_phi]
103 def MatchTrk_Info(matchtrk, clusters):
104 RV = [0., 0., 0., 0., 0., 0.]
105 RV[0] = matchtrk.getPt()
106 RV[1] = matchtrk.getPz()
107 RV[2] = (PI - matchtrk.getTheta()) * Fac
108 RV[3] = matchtrk.getPhi() * Fac
109 clusterId = matchtrk.getECLClusterId()
110 for cluster
in clusters:
111 clusterId2 = cluster.getClusterId()
112 if clusterId == clusterId2:
113 vec_ = Vec_Cluster(cluster,
False)
114 vec_cms = Vec_Cluster(cluster,
False)
120 def Max_Cluster(neuclus, CMS, n):
122 for neuclu
in neuclus:
123 NV = Vec_Cluster(neuclu, CMS)
125 if E
and len(E) >= (n + 1):
126 seqen = nlargest(n + 1, E, key=
lambda item: item[0])
129 return [-1, -999, -999]
132 def Etot_Cluster(eclcluster):
134 for clu
in eclcluster:
135 e = clu.getEnergyDep()
140 def Max_DeltPhi_trk(trk_2d_list):
142 for i, trk
in enumerate(trk_2d_list):
143 phi1 = (trk.getPhi0()) * Fac
144 for j, trk2
in enumerate(trk_2d_list):
148 phi2 = (trk2.getPhi0()) * Fac
149 Delt_phi = math.fabs(phi1 - phi2)
150 Delt.append(Delt_phi)
151 return max(Delt
or [-100.0])
154 def BhabhaVeto1(matchtrks):
156 if len(matchtrks) <= 1:
158 for i, matchtrk1
in enumerate(matchtrks):
159 if i == len(matchtrks) - 1:
161 cdctrk1 = matchtrk1.getRelatedTo(
'TRGNNTracks')
162 cluster1 = matchtrk1.getRelatedTo(
'TRGECLClusters')
163 if cdctrk1
and cluster1:
164 tanLam1 = cdctrk1.getTanLambda()
165 theta1 = math.acos(tanLam1 / math.sqrt(1. + tanLam1 * tanLam1)) * Fac
166 phi1 = (cdctrk1.getPhi0()) * Fac
169 e1 = cluster1.getEnergyDep()
170 for j, matchtrk2
in enumerate(matchtrks):
174 cdctrk2 = matchtrk2.getRelatedTo(
'TRGNNTracks')
175 cluster2 = matchtrk2.getRelatedTo(
'TRGECLClusters')
176 if cdctrk2
and cluster2:
177 tanLam2 = cdctrk2.getTanLambda()
178 theta2 = math.acos(tanLam2 / math.sqrt(1. + tanLam2 * tanLam2)) * Fac
179 phi2 = (cdctrk2.getPhi0()) * Fac
182 e2 = cluster2.getEnergyDep()
183 Delt_theta = theta1 + theta2 - 180
184 Delt_phi = math.fabs(phi1 - phi2) - 180
185 DeltArray = [Delt_theta, Delt_phi, e1, e2, e1 + e2]
186 Delt.append(DeltArray)
190 return [-999., -999., -999., -999., -999]
193 def SBhabhaVeto(matchtrks):
195 if len(matchtrks) <= 1:
197 for i, matchtrk1
in enumerate(matchtrks):
198 cdctrk1 = matchtrk1.getRelatedTo(
'TRGNNTracks')
199 cluster1 = matchtrk1.getRelatedTo(
'TRGECLClusters')
200 if cdctrk1
and cluster1:
201 e1 = cluster1.getEnergyDep()
209 def eclBhabhaVeto(eclclusters):
211 for i, cluster1
in enumerate(eclclusters):
212 if i == len(eclclusters) - 1:
214 e1 = cluster1.getEnergyDep()
215 x1 = cluster1.getPositionX()
216 y1 = cluster1.getPositionY()
217 z1 = cluster1.getPositionZ()
218 vec1 = ROOT.TVector3(x1, y1, z1)
219 theta1 = vec1.Theta()
223 for j, cluster2
in enumerate(eclclusters):
224 e2 = cluster2.getEnergyDep()
227 x2 = cluster2.getPositionX()
228 y2 = cluster2.getPositionY()
229 z2 = cluster2.getPositionZ()
230 vec2 = ROOT.TVector3(x2, y2, z2)
231 theta2 = vec2.Theta()
235 delt_theta = math.fabs((theta1 + theta2) * Fac - 180.0)
236 delt_phi = math.fabs(phi1 - phi2) * Fac - 180.0
243 bhacl = [delt_theta, delt_phi, efr, ebr, efr + ebr]
244 neclBhabha.append(bhacl)
248 def Max_DeltPhi_cluster(eclclusters):
250 for i, cluster1
in enumerate(eclclusters):
251 e1 = cluster1.getEnergyDep()
254 x1 = cluster1.getPositionX()
255 y1 = cluster1.getPositionY()
256 z1 = cluster1.getPositionZ()
257 vec1 = ROOT.TVector3(x1, y1, z1)
261 for j, cluster2
in enumerate(eclclusters):
262 e2 = cluster2.getEnergyDep()
267 x2 = cluster2.getPositionX()
268 y2 = cluster2.getPositionY()
269 z2 = cluster2.getPositionZ()
270 vec2 = ROOT.TVector3(x2, y2, z2)
274 delt_phi = math.fabs(phi1 - phi2) * 180.0 / 3.1415926
275 eclphi_col.append(delt_phi)
276 return max(eclphi_col
or [-1.0])
279 def Time_Window(clusters, eventtime):
283 event_time = ev.m_eventtiming
284 event_tot = ev.m_etot
285 energy.append([event_time, event_tot])
287 tmp = max(energy, key=
lambda item: item[1])
290 event_time = -999999.
291 for cluster
in clusters:
292 ctime_ori = cluster.getTimeAve()
293 ctime = ctime_ori - event_time
294 if math.fabs(ctime) < 100:
295 new_clusters.append(cluster)
299 def Cluster_Threshold(clusters, threshold, CMS):
301 for cluster
in clusters:
304 newv = Vec_Cluster(cluster, CMS)
307 eng = cluster.getEnergyDep()
309 new_clusters.append(cluster)
313 def Back_to_Back(clusters1, clusters2):
315 for cluster1
in clusters1:
316 cid1 = cluster1.getClusterId()
317 vec1 = Vec_Cluster(cluster1,
False)
320 for cluster2
in clusters2:
321 cid2 = cluster2.getClusterId()
324 vec2 = Vec_Cluster(cluster2,
False)
327 delttheta = math.fabs(theta1 + theta2 - 180)
328 deltphi = math.fabs(math.fabs(phi1 - phi2) - 180)
329 if delttheta < 100
and deltphi < 100:
334 def Time_Cluster(clusters, eventtime):
339 event_time = ev.m_eventtiming
340 event_tot = ev.m_etot
341 energy.append([event_time, event_tot])
343 tmp = max(energy, key=
lambda item: item[1])
345 for cluster
in clusters:
346 ctime_ori = cluster.getTimeAve()
347 ctime = ctime_ori - event_time
348 time_list.append(ctime)
352 def PrintBranchDef():
354 print(
'ntrk_2dfinder: # 2d finder track')
355 print(
'ntrk_2dfitter: # 2d fitter track')
356 print(
'ntrk_3dfitter: # 3d fitter track')
357 print(
'ntrk_NN : # Neuro network track')
358 print(
'ntrk_2Dmatch : # 2d finder track w/ associated ecl cluster')
359 print(
'ntrk_3Dmatch : # 3d finder track w/ associated ecl cluster')
360 print(
'ntrk_klm : # KLM track')
361 print(
'nhit_klm : # KLM hit',
'\n')
363 print(
'ncluster: # ecl cluster')
364 print(
'ncluster_1000b: # ecl cluster with threshold >1.0GeV, exclude TC ID 1,2, 17')
365 print(
'ncluster_2000e: # ecl cluster with threshold >2.0GeV in TC ID 1, 17')
366 print(
'max_cluster[3]: [energy, theta, phi] of the largest energetic ecl cluster')
367 print(
'smax_cluster[3]: [energy, theta, phi] of the secondary energetic ecl cluster')
368 print(
'ncluster_neutral: # ecl cluster w/o associated cdc track')
369 print(
'max_cluster_neutral[3]: [energy, theta, phi] of the largest energetic ecl neutral cluster')
370 print(
'smax_cluster_neutral[3]: [energy, theta, phi] of the secondary energetic ecl neutral cluster')
371 print(
'time_cluster: the cluster timing obtained with cluster.timing-event.timing',
'\n')
373 print(
'nbbc: # back to back cluster pairs')
374 print(
'nbbtc: # back to back track and cluster pairs')
375 print(
'bhabha: bhabha veto logic, 1: bhabha, 0: non bhabha')
376 print(
'sbhabha: bhabha with single track veto logic, 1: bhabha, 0: non bhabha')
377 print(
'eclbhabha: eclbhabha veto logic, 1: bhabha, 0: non bhabha')
378 print(
'bhabha_var[5]: variables used in bhabha logic, [Delt_theta, Delt_phi, E1, E2, E1+E2]')
379 print(
' For two tracks: Delt_theta: theta1+theta2-180, Delt_phi: |phi1-phi2|-180')
380 print(
' E1, E2 are the ecl clusters energy accociated with the two tracks')
381 print(
'eclbhabha_var[5]: variables used in eclbhabha logic, [Delt_theta, Delt_phi, E1, E2, E1+E2]')
382 print(
' For two eclclusers: Delt_theta: theta1+theta2-180, Delt_phi: |phi1-phi2|-180')
383 print(
' E1, E2 are the ecl clusters energy')
388 """This module is to calculate some variables which are useful for the trigger development"""
390 file = ROOT.TFile(argvs[2],
'recreate')
392 tgrl = ROOT.TTree(
'tgrl',
'tree with GRL_Logic')
394 ntrk_2dfinder_t = array(
'i', [-1])
396 ntrk_2dfitter_t = array(
'i', [-1])
398 ntrk_3dfitter_t = array(
'i', [-1])
400 ntrk_NN_t = array(
'i', [-1])
402 ntrk_2Dmatch_t = array(
'i', [-1])
404 ntrk_3Dmatch_t = array(
'i', [-1])
406 max_deltphi_2dfinder_t = array(
'f', [0.0])
408 cpair_t = array(
'i', 8 * [0])
412 etot_t = array(
'f', [0.0])
414 ncluster_1000b_t = array(
'i', [-1])
416 ncluster_2000e_t = array(
'i', [-1])
418 ncluster_t = array(
'i', [-1])
420 ncluster_neutral_t = array(
'i', [-1])
422 max_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
424 max_cms_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
426 max_cluster_t = array(
'f', ncomp_clu * [0.0])
428 max_cms_cluster_t = array(
'f', ncomp_clu * [0.0])
430 smax_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
432 smax_cms_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
434 smax_cluster_t = array(
'f', ncomp_clu * [0.0])
436 smax_cms_cluster_t = array(
'f', ncomp_clu * [0.0])
438 max_deltphi_cluster_t = array(
'f', [0.0])
440 time_cluster_t = array(
'f', 100 * [-99999.])
442 ntrk_klm_t = array(
'i', [-1])
444 nhit_klm_t = array(
'i', [-1])
446 npair_tc_t = array(
'i', [-1])
448 npair_cc_t = array(
'i', [-1])
452 bhabha_t = array(
'i', [0])
454 sbhabha_t = array(
'i', [0])
456 eclbhabha_t = array(
'i', [0])
458 bhabha_var_t = array(
'f', nbha_var * [0.0])
460 eclbhabha_var_t = array(
'f', nbha_var * [0.0])
463 tgrl.Branch(
'ntrk_2dfinder', ntrk_2dfinder_t,
'ntrk_2dfinder/i')
464 tgrl.Branch(
'ntrk_2dfitter', ntrk_2dfitter_t,
'ntrk_2dfitter/i')
465 tgrl.Branch(
'ntrk_3dfitter', ntrk_3dfitter_t,
'ntrk_3dfitter/i')
466 tgrl.Branch(
'ntrk_NN', ntrk_NN_t,
'ntrk_NN/i')
467 tgrl.Branch(
'ntrk_2Dmatch', ntrk_2Dmatch_t,
'ntrk_2Dmatch/i')
468 tgrl.Branch(
'ntrk_3Dmatch', ntrk_3Dmatch_t,
'ntrk_3Dmatch/i')
469 tgrl.Branch(
'npair_tc', npair_tc_t,
'npair_tc/i')
470 tgrl.Branch(
'npair_cc', npair_cc_t,
'npair_cc/i')
471 tgrl.Branch(
'max_deltphi_2dfinder', max_deltphi_2dfinder_t,
'max_deltphi_2dfinder/f')
474 tgrl.Branch(
'etot', etot_t,
'etot/f')
475 tgrl.Branch(
'ncluster', ncluster_t,
'ncluster/i')
476 tgrl.Branch(
'ncluster_1000b', ncluster_1000b_t,
'ncluster_1000b/i')
477 tgrl.Branch(
'ncluster_2000e', ncluster_2000e_t,
'ncluster_2000e/i')
478 tgrl.Branch(
'ncluster_neutral', ncluster_neutral_t,
'ncluster_neutral/i')
479 tgrl.Branch(
'max_cluster_neutral', max_cluster_neutral_t,
'max_cluster_neutral[3]/f')
481 tgrl.Branch(
'max_cluster', max_cluster_t,
'max_cluster[3]/f')
483 tgrl.Branch(
'smax_cluster_neutral', smax_cluster_neutral_t,
'smax_cluster_neutral[3]/f')
485 tgrl.Branch(
'smax_cluster', smax_cluster_t,
'smax_cluster[3]/f')
487 tgrl.Branch(
'max_deltphi_cluster', max_deltphi_cluster_t,
'max_deltphi_cluster/f')
488 tgrl.Branch(
'time_cluster', time_cluster_t,
'time_cluster[100]/f')
491 tgrl.Branch(
'ntrk_klm', ntrk_klm_t,
'ntrk_klm/i')
492 tgrl.Branch(
'nhit_klm', nhit_klm_t,
'nhit_klm/i')
495 tgrl.Branch(
'bhabha', bhabha_t,
'bhabha/i')
496 tgrl.Branch(
'sbhabha', sbhabha_t,
'sbhabha/i')
497 tgrl.Branch(
'eclbhabha', eclbhabha_t,
'eclbhabha/i')
498 tgrl.Branch(
'bhabha_var', bhabha_var_t,
'bhabha_var[5]/f')
499 tgrl.Branch(
'eclbhabha_var', eclbhabha_var_t,
'eclbhabha_var[5]/f')
501 tgrl.Branch(
'cpair', cpair_t,
'cpair[8]/i')
505 MCPart = Belle2.PyStoreArray('MCParticles')
507 self.eplus_t[0]=Theta_MCTruth(MCPart[2])
508 self.eminus_t[0]=Theta_MCTruth(MCPart[3])
509 if self.eplus_t[0] > self.eminus_t[0]:
510 self.efrd_t[0]=self.eminus_t[0]
511 self.ebkd_t[0]=self.eplus_t[0]
513 self.efrd_t[0]=self.eplus_t[0]
514 self.ebkd_t[0]=self.eminus_t[0]
515 count_part = Count_part_inDetector(MCPart)
516 self.n_par_t[0]=count_part[0]
517 self.n_par_t[1]=count_part[1]
518 count_mcpart = Count_mcpart(MCPart)
519 self.n_mcpar_t[0]=count_mcpart[0]
520 self.n_mcpar_t[1]=count_mcpart[1]
535 clusters = Time_Window(clusters_original, eventtime)
538 time_clu_list = Time_Cluster(clusters_original, eventtime)
539 for i
in range(len(time_clu_list)):
543 clusters_100 = Cluster_Threshold(clusters, 0.1,
False)
544 clusters_100_cms = Cluster_Threshold(clusters, 0.1,
False)
545 clusters_300_cms = Cluster_Threshold(clusters, 0.3,
False)
546 clusters_400_cms = Cluster_Threshold(clusters, 0.4,
False)
547 clusters_500_cms = Cluster_Threshold(clusters, 0.5,
False)
548 clusters_700_cms = Cluster_Threshold(clusters, 0.7,
False)
549 clusters_1000_cms = Cluster_Threshold(clusters, 1.0,
False)
550 clusters_2000_cms = Cluster_Threshold(clusters, 2.0,
False)
551 clusters_2500_cms = Cluster_Threshold(clusters, 2.5,
False)
552 self.
cpair_tcpair_t[0] = Back_to_Back(clusters, clusters)
553 self.
cpair_tcpair_t[1] = Back_to_Back(clusters_300_cms, clusters)
554 self.
cpair_tcpair_t[2] = Back_to_Back(clusters_400_cms, clusters)
555 self.
cpair_tcpair_t[3] = Back_to_Back(clusters_500_cms, clusters)
556 self.
cpair_tcpair_t[4] = Back_to_Back(clusters_700_cms, clusters)
557 self.
cpair_tcpair_t[5] = Back_to_Back(clusters_1000_cms, clusters)
558 self.
cpair_tcpair_t[6] = Back_to_Back(clusters_2000_cms, clusters)
559 self.
cpair_tcpair_t[7] = Back_to_Back(clusters_2500_cms, clusters)
570 self.
npair_tc_tnpair_tc_t[0] = trginfo.getNbbTrkCluster()
571 self.
npair_cc_tnpair_cc_t[0] = trginfo.getNbbCluster()
576 neutral_clusters = NeutralCluster(clusters)
578 max_cluster_neu = Max_Cluster(neutral_clusters,
False, 0)
579 max_cluster = Max_Cluster(clusters,
False, 0)
582 smax_cluster_neu = Max_Cluster(neutral_clusters,
False, 1)
583 smax_cluster = Max_Cluster(clusters,
False, 1)
596 self.
etot_tetot_t[0] = Etot_Cluster(clusters)
599 bhabhaveto_1 = BhabhaVeto1(matchlist)
601 for bha
in bhabhaveto_1:
602 if math.fabs(bha[0]) < 50
and math.fabs(bha[0]) > 10
and math.fabs(bha[1]) < 20:
603 if bha[2] > 2.0
and bha[3] > 2.0
and bha[4] > 6.0:
604 if bha[2] > 3.0
or bha[3] > 3.0:
605 if len(trk_2d_finder) == 2:
607 self.
bhabha_tbhabha_t[0] = bha_logic
609 if len(bhabhaveto_1) >= 1:
610 for i
in range(self.
nbha_varnbha_var):
614 eclbhabhaveto = eclBhabhaVeto(clusters)
616 for eclbha
in eclbhabhaveto:
617 if math.fabs(eclbha[0]) < 50
and math.fabs(eclbha[1]) < 50:
618 if eclbha[2] > 2.0
and eclbha[3] > 2.0
and eclbha[4] > 6.0:
619 if eclbha[2] > 3.0
or eclbha[3] > 3.0:
625 if len(trk_2d_finder) == 1:
626 if eclbha_logic == 1:
627 ecol = SBhabhaVeto(matchlist)
634 if len(eclbhabhaveto) >= 1:
635 for i
in range(self.
nbha_varnbha_var):
640 """Write and close the file"""
642 self.
filefile.Write()
643 self.
filefile.Close()
647 if __name__ ==
"__main__":
648 main = b2.create_path()
649 ma.inputMdst(
'default', argvs[1], main)
The DetectorSet class for sets of detector IDs in the form of EDetector values.
a (simplified) python wrapper for StoreArray.
a (simplified) python wrapper for StoreObjPtr.
cpair_t
#cluster pairs with different energy threshold
eclbhabha_var_t
variables used in eclbhabha logic, [Delt_theta, Delt_phi, E1, E2, E1+E2]
ncluster_1000b_t
#ecl cluster with threshold >1.0GeV, exclude TC ID 1,2, 17
max_cluster_neutral_t
[energy, theta, phi] of the largest energetic ecl neutral cluster
ntrk_3dfitter_t
#3d fitter tracks
eclbhabha_t
eclbhabha veto logic, 1: bhabha, 0: non bhabha
etot_t
the total deposited cluster energy in ecl
smax_cluster_neutral_t
[energy, theta, phi] of the secondary energetic ecl neutral cluster
int ncomp_clu
#array components
smax_cluster_t
[energy, theta, phi] of the secondary energetic ecl cluster
tgrl
the tree in the output file
ntrk_2dfinder_t
#2d finder tracks
ntrk_3Dmatch_t
#3d matched tracks
time_cluster_t
the timing of clusters
bhabha_t
bhabha veto logic, 1: bhabha, 0: non bhabha
ntrk_2dfitter_t
#2d fitter tracks
max_cluster_t
[energy, theta, phi] of the largest energetic ecl cluster
int nbha_var
#array components
sbhabha_t
bhabha with single track veto logic, 1: bhabha, 0: non bhabha
max_deltphi_cluster_t
max delt phi angle between two clusters
bhabha_var_t
variables used in bhabha logic, [Delt_theta, Delt_phi, E1, E2, E1+E2]
max_deltphi_2dfinder_t
#max phi angle between two 2d finder tracks