7 from ROOT
import Belle2
8 from modularAnalysis
import *
9 from array
import array
11 from heapq
import nlargest
12 from effCalculation
import EffCalculation
22 sys.exit(
"ztsim02.py> # of arg is strange. Exit.")
25 def Count_mcpart(parts):
30 pdglist = [11, 13, 2212, 321, 211]
31 pdg = math.fabs(part.getPDG())
32 charge = part.getCharge()
33 theta = part.getMomentum().Theta() * Fac
34 if theta > dec_cdc[0]
and theta < dec_cdc[1]:
38 if math.fabs(pdg) == 22:
39 if theta > dec_ecl[0]
and theta < dec_ecl[1]:
44 def Count_part_inDetector(parts):
48 charge = part.getCharge()
51 if math.fabs(charge) > 0:
53 elif math.fabs(pdg) == 22:
58 def Theta_MCTruth(part):
59 theta = part.getMomentum().Theta() * Fac
63 def NeutralCluster(clusters):
65 for cluster
in clusters:
66 if cluster.getRelatedFrom(
'TRGNNTracks'):
69 neutral_cluster.append(cluster)
70 return neutral_cluster
73 def Vec_Cluster(cluster, CMS=False):
74 x = cluster.getPositionX()
75 y = cluster.getPositionY()
76 z = cluster.getPositionZ()
77 e = cluster.getEnergyDep()
78 vec = ROOT.TVector3(x, y, z)
81 v_mom = ROOT.TLorentzVector(vec_unit, e)
83 new_vec = trans.rotateLabToCms() * v_mom
86 new_theta = new_vec.Theta() * Fac
87 new_phi = new_vec.Phi() * Fac
91 NVC = [new_e, new_theta, new_phi]
95 def MatchTrk_Info(matchtrk, clusters):
96 RV = [0., 0., 0., 0., 0., 0.]
97 RV[0] = matchtrk.getPt()
98 RV[1] = matchtrk.getPz()
99 RV[2] = (PI - matchtrk.getTheta()) * Fac
100 RV[3] = matchtrk.getPhi() * Fac
101 clusterId = matchtrk.getECLClusterId()
102 for cluster
in clusters:
103 clusterId2 = cluster.getClusterId()
104 if clusterId == clusterId2:
105 vec_ = Vec_Cluster(cluster,
False)
106 vec_cms = Vec_Cluster(cluster,
False)
112 def Max_Cluster(neuclus, CMS, n):
114 for neuclu
in neuclus:
115 NV = Vec_Cluster(neuclu, CMS)
117 if E
and len(E) >= (n + 1):
118 seqen = nlargest(n + 1, E, key=
lambda item: item[0])
121 return [-1, -999, -999]
124 def Etot_Cluster(eclcluster):
126 for clu
in eclcluster:
127 e = clu.getEnergyDep()
132 def Max_DeltPhi_trk(trk_2d_list):
134 for i, trk
in enumerate(trk_2d_list):
135 phi1 = (trk.getPhi0()) * Fac
136 for j, trk2
in enumerate(trk_2d_list):
140 phi2 = (trk2.getPhi0()) * Fac
141 Delt_phi = math.fabs(phi1 - phi2)
142 Delt.append(Delt_phi)
143 return max(Delt
or [-100.0])
146 def BhabhaVeto1(matchtrks):
148 if len(matchtrks) <= 1:
150 for i, matchtrk1
in enumerate(matchtrks):
151 if i == len(matchtrks) - 1:
153 cdctrk1 = matchtrk1.getRelatedTo(
'TRGNNTracks')
154 cluster1 = matchtrk1.getRelatedTo(
'TRGECLClusters')
155 if cdctrk1
and cluster1:
156 tanLam1 = cdctrk1.getTanLambda()
157 theta1 = math.acos(tanLam1 / math.sqrt(1. + tanLam1 * tanLam1)) * Fac
158 phi1 = (cdctrk1.getPhi0()) * Fac
161 e1 = cluster1.getEnergyDep()
162 for j, matchtrk2
in enumerate(matchtrks):
166 cdctrk2 = matchtrk2.getRelatedTo(
'TRGNNTracks')
167 cluster2 = matchtrk2.getRelatedTo(
'TRGECLClusters')
168 if cdctrk2
and cluster2:
169 tanLam2 = cdctrk2.getTanLambda()
170 theta2 = math.acos(tanLam2 / math.sqrt(1. + tanLam2 * tanLam2)) * Fac
171 phi2 = (cdctrk2.getPhi0()) * Fac
174 e2 = cluster2.getEnergyDep()
175 Delt_theta = theta1 + theta2 - 180
176 Delt_phi = math.fabs(phi1 - phi2) - 180
177 DeltArray = [Delt_theta, Delt_phi, e1, e2, e1 + e2]
178 Delt.append(DeltArray)
182 return [-999., -999., -999., -999., -999]
185 def SBhabhaVeto(matchtrks):
187 if len(matchtrks) <= 1:
189 for i, matchtrk1
in enumerate(matchtrks):
190 cdctrk1 = matchtrk1.getRelatedTo(
'TRGNNTracks')
191 cluster1 = matchtrk1.getRelatedTo(
'TRGECLClusters')
192 if cdctrk1
and cluster1:
193 e1 = cluster1.getEnergyDep()
201 def eclBhabhaVeto(eclclusters):
203 for i, cluster1
in enumerate(eclclusters):
204 if i == len(eclclusters) - 1:
206 e1 = cluster1.getEnergyDep()
207 x1 = cluster1.getPositionX()
208 y1 = cluster1.getPositionY()
209 z1 = cluster1.getPositionZ()
210 vec1 = ROOT.TVector3(x1, y1, z1)
211 theta1 = vec1.Theta()
215 for j, cluster2
in enumerate(eclclusters):
216 e2 = cluster2.getEnergyDep()
219 x2 = cluster2.getPositionX()
220 y2 = cluster2.getPositionY()
221 z2 = cluster2.getPositionZ()
222 vec2 = ROOT.TVector3(x2, y2, z2)
223 theta2 = vec2.Theta()
227 delt_theta = math.fabs((theta1 + theta2) * Fac - 180.0)
228 delt_phi = math.fabs(phi1 - phi2) * Fac - 180.0
235 bhacl = [delt_theta, delt_phi, efr, ebr, efr + ebr]
236 neclBhabha.append(bhacl)
240 def Max_DeltPhi_cluster(eclclusters):
242 for i, cluster1
in enumerate(eclclusters):
243 e1 = cluster1.getEnergyDep()
246 x1 = cluster1.getPositionX()
247 y1 = cluster1.getPositionY()
248 z1 = cluster1.getPositionZ()
249 vec1 = ROOT.TVector3(x1, y1, z1)
253 for j, cluster2
in enumerate(eclclusters):
254 e2 = cluster2.getEnergyDep()
259 x2 = cluster2.getPositionX()
260 y2 = cluster2.getPositionY()
261 z2 = cluster2.getPositionZ()
262 vec2 = ROOT.TVector3(x2, y2, z2)
266 delt_phi = math.fabs(phi1 - phi2) * 180.0 / 3.1415926
267 eclphi_col.append(delt_phi)
268 return max(eclphi_col
or [-1.0])
271 def Time_Window(clusters, eventtime):
275 event_time = ev.m_eventtiming
276 event_tot = ev.m_etot
277 energy.append([event_time, event_tot])
279 tmp = max(energy, key=
lambda item: item[1])
282 event_time = -999999.
283 for cluster
in clusters:
284 ctime_ori = cluster.getTimeAve()
285 ctime = ctime_ori - event_time
286 if math.fabs(ctime) < 100:
287 new_clusters.append(cluster)
291 def Cluster_Threshold(clusters, threshold, CMS):
293 for cluster
in clusters:
296 newv = Vec_Cluster(cluster, CMS)
299 eng = cluster.getEnergyDep()
301 new_clusters.append(cluster)
305 def Back_to_Back(clusters1, clusters2):
307 for cluster1
in clusters1:
308 cid1 = cluster1.getClusterId()
309 vec1 = Vec_Cluster(cluster1,
False)
312 for cluster2
in clusters2:
313 cid2 = cluster2.getClusterId()
316 vec2 = Vec_Cluster(cluster2,
False)
319 delttheta = math.fabs(theta1 + theta2 - 180)
320 deltphi = math.fabs(math.fabs(phi1 - phi2) - 180)
321 if delttheta < 100
and deltphi < 100:
326 def Time_Cluster(clusters, eventtime):
331 event_time = ev.m_eventtiming
332 event_tot = ev.m_etot
333 energy.append([event_time, event_tot])
335 tmp = max(energy, key=
lambda item: item[1])
337 for cluster
in clusters:
338 ctime_ori = cluster.getTimeAve()
339 ctime = ctime_ori - event_time
340 time_list.append(ctime)
344 def PrintBranchDef():
346 print(
'ntrk_2dfinder: # 2d finder track')
347 print(
'ntrk_2dfitter: # 2d fitter track')
348 print(
'ntrk_3dfitter: # 3d fitter track')
349 print(
'ntrk_NN : # Neuro network track')
350 print(
'ntrk_2Dmatch : # 2d finder track w/ associated ecl cluster')
351 print(
'ntrk_3Dmatch : # 3d finder track w/ associated ecl cluster')
352 print(
'ntrk_klm : # KLM track')
353 print(
'nhit_klm : # KLM hit',
'\n')
355 print(
'ncluster: # ecl cluster')
356 print(
'ncluster_1000b: # ecl cluster with threshold >1.0GeV, exclude TC ID 1,2, 17')
357 print(
'ncluster_2000e: # ecl cluster with threshold >2.0GeV in TC ID 1, 17')
358 print(
'max_cluster[3]: [energy, theta, phi] of the largest energetic ecl cluster')
359 print(
'smax_cluster[3]: [energy, theta, phi] of the secondary energetic ecl cluster')
360 print(
'ncluster_neutral: # ecl cluster w/o associated cdc track')
361 print(
'max_cluster_neutral[3]: [energy, theta, phi] of the largest energetic ecl neutral cluster')
362 print(
'smax_cluster_neutral[3]: [energy, theta, phi] of the secondary energetic ecl neutral cluster')
363 print(
'time_cluster: the cluster timing obtained with cluster.timing-event.timing',
'\n')
365 print(
'nbbc: # back to back cluster pairs')
366 print(
'nbbtc: # back to back track and cluster pairs')
367 print(
'bhabha: bhabha veto logic, 1: bhabha, 0: non bhabha')
368 print(
'sbhabha: bhabha with single track veto logic, 1: bhabha, 0: non bhabha')
369 print(
'eclbhabha: eclbhabha veto logic, 1: bhabha, 0: non bhabha')
370 print(
'bhabha_var[5]: variables used in bhabha logic, [Delt_theta, Delt_phi, E1, E2, E1+E2]')
371 print(
' For two tracks: Delt_theta: theta1+theta2-180, Delt_phi: |phi1-phi2|-180')
372 print(
' E1, E2 are the ecl clusters energy accociated with the two tracks')
373 print(
'eclbhabha_var[5]: variables used in eclbhabha logic, [Delt_theta, Delt_phi, E1, E2, E1+E2]')
374 print(
' For two eclclusers: Delt_theta: theta1+theta2-180, Delt_phi: |phi1-phi2|-180')
375 print(
' E1, E2 are the ecl clusters energy')
380 """This module is to calculate some variables which are useful for the trigger development"""
382 file = ROOT.TFile(argvs[2],
'recreate')
384 tgrl = ROOT.TTree(
'tgrl',
'tree with GRL_Logic')
386 ntrk_2dfinder_t = array(
'i', [-1])
388 ntrk_2dfitter_t = array(
'i', [-1])
390 ntrk_3dfitter_t = array(
'i', [-1])
392 ntrk_NN_t = array(
'i', [-1])
394 ntrk_2Dmatch_t = array(
'i', [-1])
396 ntrk_3Dmatch_t = array(
'i', [-1])
398 max_deltphi_2dfinder_t = array(
'f', [0.0])
400 cpair_t = array(
'i', 8 * [0])
404 etot_t = array(
'f', [0.0])
406 ncluster_1000b_t = array(
'i', [-1])
408 ncluster_2000e_t = array(
'i', [-1])
410 ncluster_t = array(
'i', [-1])
412 ncluster_neutral_t = array(
'i', [-1])
414 max_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
416 max_cms_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
418 max_cluster_t = array(
'f', ncomp_clu * [0.0])
420 max_cms_cluster_t = array(
'f', ncomp_clu * [0.0])
422 smax_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
424 smax_cms_cluster_neutral_t = array(
'f', ncomp_clu * [0.0])
426 smax_cluster_t = array(
'f', ncomp_clu * [0.0])
428 smax_cms_cluster_t = array(
'f', ncomp_clu * [0.0])
430 max_deltphi_cluster_t = array(
'f', [0.0])
432 time_cluster_t = array(
'f', 100 * [-99999.])
434 ntrk_klm_t = array(
'i', [-1])
436 nhit_klm_t = array(
'i', [-1])
438 npair_tc_t = array(
'i', [-1])
440 npair_cc_t = array(
'i', [-1])
444 bhabha_t = array(
'i', [0])
446 sbhabha_t = array(
'i', [0])
448 eclbhabha_t = array(
'i', [0])
450 bhabha_var_t = array(
'f', nbha_var * [0.0])
452 eclbhabha_var_t = array(
'f', nbha_var * [0.0])
455 tgrl.Branch(
'ntrk_2dfinder', ntrk_2dfinder_t,
'ntrk_2dfinder/i')
456 tgrl.Branch(
'ntrk_2dfitter', ntrk_2dfitter_t,
'ntrk_2dfitter/i')
457 tgrl.Branch(
'ntrk_3dfitter', ntrk_3dfitter_t,
'ntrk_3dfitter/i')
458 tgrl.Branch(
'ntrk_NN', ntrk_NN_t,
'ntrk_NN/i')
459 tgrl.Branch(
'ntrk_2Dmatch', ntrk_2Dmatch_t,
'ntrk_2Dmatch/i')
460 tgrl.Branch(
'ntrk_3Dmatch', ntrk_3Dmatch_t,
'ntrk_3Dmatch/i')
461 tgrl.Branch(
'npair_tc', npair_tc_t,
'npair_tc/i')
462 tgrl.Branch(
'npair_cc', npair_cc_t,
'npair_cc/i')
463 tgrl.Branch(
'max_deltphi_2dfinder', max_deltphi_2dfinder_t,
'max_deltphi_2dfinder/f')
466 tgrl.Branch(
'etot', etot_t,
'etot/f')
467 tgrl.Branch(
'ncluster', ncluster_t,
'ncluster/i')
468 tgrl.Branch(
'ncluster_1000b', ncluster_1000b_t,
'ncluster_1000b/i')
469 tgrl.Branch(
'ncluster_2000e', ncluster_2000e_t,
'ncluster_2000e/i')
470 tgrl.Branch(
'ncluster_neutral', ncluster_neutral_t,
'ncluster_neutral/i')
471 tgrl.Branch(
'max_cluster_neutral', max_cluster_neutral_t,
'max_cluster_neutral[3]/f')
473 tgrl.Branch(
'max_cluster', max_cluster_t,
'max_cluster[3]/f')
475 tgrl.Branch(
'smax_cluster_neutral', smax_cluster_neutral_t,
'smax_cluster_neutral[3]/f')
477 tgrl.Branch(
'smax_cluster', smax_cluster_t,
'smax_cluster[3]/f')
479 tgrl.Branch(
'max_deltphi_cluster', max_deltphi_cluster_t,
'max_deltphi_cluster/f')
480 tgrl.Branch(
'time_cluster', time_cluster_t,
'time_cluster[100]/f')
483 tgrl.Branch(
'ntrk_klm', ntrk_klm_t,
'ntrk_klm/i')
484 tgrl.Branch(
'nhit_klm', nhit_klm_t,
'nhit_klm/i')
487 tgrl.Branch(
'bhabha', bhabha_t,
'bhabha/i')
488 tgrl.Branch(
'sbhabha', sbhabha_t,
'sbhabha/i')
489 tgrl.Branch(
'eclbhabha', eclbhabha_t,
'eclbhabha/i')
490 tgrl.Branch(
'bhabha_var', bhabha_var_t,
'bhabha_var[5]/f')
491 tgrl.Branch(
'eclbhabha_var', eclbhabha_var_t,
'eclbhabha_var[5]/f')
493 tgrl.Branch(
'cpair', cpair_t,
'cpair[8]/i')
497 MCPart = Belle2.PyStoreArray('MCParticles')
499 self.eplus_t[0]=Theta_MCTruth(MCPart[2])
500 self.eminus_t[0]=Theta_MCTruth(MCPart[3])
501 if self.eplus_t[0] > self.eminus_t[0]:
502 self.efrd_t[0]=self.eminus_t[0]
503 self.ebkd_t[0]=self.eplus_t[0]
505 self.efrd_t[0]=self.eplus_t[0]
506 self.ebkd_t[0]=self.eminus_t[0]
507 count_part = Count_part_inDetector(MCPart)
508 self.n_par_t[0]=count_part[0]
509 self.n_par_t[1]=count_part[1]
510 count_mcpart = Count_mcpart(MCPart)
511 self.n_mcpar_t[0]=count_mcpart[0]
512 self.n_mcpar_t[1]=count_mcpart[1]
527 clusters = Time_Window(clusters_original, eventtime)
530 time_clu_list = Time_Cluster(clusters_original, eventtime)
531 for i
in range(len(time_clu_list)):
535 clusters_100 = Cluster_Threshold(clusters, 0.1,
False)
536 clusters_100_cms = Cluster_Threshold(clusters, 0.1,
False)
537 clusters_300_cms = Cluster_Threshold(clusters, 0.3,
False)
538 clusters_400_cms = Cluster_Threshold(clusters, 0.4,
False)
539 clusters_500_cms = Cluster_Threshold(clusters, 0.5,
False)
540 clusters_700_cms = Cluster_Threshold(clusters, 0.7,
False)
541 clusters_1000_cms = Cluster_Threshold(clusters, 1.0,
False)
542 clusters_2000_cms = Cluster_Threshold(clusters, 2.0,
False)
543 clusters_2500_cms = Cluster_Threshold(clusters, 2.5,
False)
544 self.
cpair_t[0] = Back_to_Back(clusters, clusters)
545 self.
cpair_t[1] = Back_to_Back(clusters_300_cms, clusters)
546 self.
cpair_t[2] = Back_to_Back(clusters_400_cms, clusters)
547 self.
cpair_t[3] = Back_to_Back(clusters_500_cms, clusters)
548 self.
cpair_t[4] = Back_to_Back(clusters_700_cms, clusters)
549 self.
cpair_t[5] = Back_to_Back(clusters_1000_cms, clusters)
550 self.
cpair_t[6] = Back_to_Back(clusters_2000_cms, clusters)
551 self.
cpair_t[7] = Back_to_Back(clusters_2500_cms, clusters)
562 self.
npair_tc_t[0] = trginfo.getNbbTrkCluster()
568 neutral_clusters = NeutralCluster(clusters)
570 max_cluster_neu = Max_Cluster(neutral_clusters,
False, 0)
571 max_cluster = Max_Cluster(clusters,
False, 0)
574 smax_cluster_neu = Max_Cluster(neutral_clusters,
False, 1)
575 smax_cluster = Max_Cluster(clusters,
False, 1)
588 self.
etot_t[0] = Etot_Cluster(clusters)
591 bhabhaveto_1 = BhabhaVeto1(matchlist)
593 for bha
in bhabhaveto_1:
594 if math.fabs(bha[0]) < 50
and math.fabs(bha[0]) > 10
and math.fabs(bha[1]) < 20:
595 if bha[2] > 2.0
and bha[3] > 2.0
and bha[4] > 6.0:
596 if bha[2] > 3.0
or bha[3] > 3.0:
597 if len(trk_2d_finder) == 2:
601 if len(bhabhaveto_1) >= 1:
606 eclbhabhaveto = eclBhabhaVeto(clusters)
608 for eclbha
in eclbhabhaveto:
609 if math.fabs(eclbha[0]) < 50
and math.fabs(eclbha[1]) < 50:
610 if eclbha[2] > 2.0
and eclbha[3] > 2.0
and eclbha[4] > 6.0:
611 if eclbha[2] > 3.0
or eclbha[3] > 3.0:
617 if len(trk_2d_finder) == 1:
618 if eclbha_logic == 1:
619 ecol = SBhabhaVeto(matchlist)
626 if len(eclbhabhaveto) >= 1:
632 """Write and close the file"""
639 if __name__ ==
"__main__":
641 inputMdst(
'default', argvs[1], main)