ATLAS Offline Software
TRTCalib_cfilter.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 # Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
4 
5 import os,sys, math
6 
7 def tshift_poly(dt,p0,p1,p2,p3):
8  p0n = p0 + p1*dt + p2*dt*dt + p3*dt*dt*dt
9  p1n = p1 + 2*p2*dt + 3*p3*dt*dt
10  p2n = p2 + 3*p3*dt
11  p3n = p3
12  return p0n,p1n,p2n,p3n
13 
14 def p3root(a, b, c, d): # this is to replace the numpy root finder, should numpy not be available
15 
16  r = 18.0
17  pi = 3.14159265359
18  if a*a < 1.0e-48 :
19  return r
20  p = (3*a*c-b*b)/(3*a*a)
21  q = (2*b*b*b-9*a*b*c+27*a*a*d)/(27*a*a*a)
22 
23  offset = -b/(3*a)
24  discriminant = 4*p*p*p+27*q*q
25  if discriminant > 0:
26  v = math.sqrt(q*q/4 + p*p*p/27 )
27  if (-q/2 + v) >0 and (-q/2-v) > 0 :
28  r = pow(-q/2 + v,0.3333333)+pow(-q/2-v,0.3333333)+offset
29  if discriminant < 0 and p < 0:
30  cphi = (3.0/2.0)*(q/p)*math.sqrt(-3/p)
31  phi = math.acos(cphi)
32  z1 = 2*math.sqrt(-p/3)*math.cos(phi/3) + offset
33  z2 = 2*math.sqrt(-p/3)*math.cos((phi+2*pi)/3) + offset
34  z3 = 2*math.sqrt(-p/3)*math.cos((phi+4*pi)/3) + offset
35  if abs(z1-18)<abs(z2-18) and abs(z1-18)<abs(z3-18):
36  r = z1
37  elif abs(z2-18)<abs(z3-18):
38  r = z2
39  else:
40  r = z3
41 
42  return r
43 
44 keeprt=False
45 keept0=False
46 shiftrt=False
47 shiftt0=False
48 do_global=False
49 
50 if sys.argv[1]=="t0shift":
51  do_global=True
52  keeprt=True
53  keept0=True
54  shiftt0=True
55  Dt0=float(sys.argv[2])
56  olddbs=open(sys.argv[3]).readlines()
57  newdbs=open(sys.argv[3]).readlines()
58 else:
59  olddbs=open(sys.argv[1]).readlines()
60  newdbs=open(sys.argv[2]).readlines()
61 
62  if len(sys.argv)>4:
63  olddicts=open(sys.argv[3]).readlines()
64  keeprt=sys.argv[4].find("keeprt")>=0
65  keept0=sys.argv[4].find("keept0")>=0
66  shiftrt=sys.argv[4].find("shiftrt")>=0
67  shiftt0=sys.argv[4].find("shiftt0")>=0
68  if keeprt:
69  print ("cfilter called with 4 arguments. Option keeprt")
70  if keept0:
71  print ("cfilter called with 4 arguments. Option keept0")
72  if shiftrt:
73  print ("cfilter called with 4 arguments. Option shiftrt")
74  if shiftt0:
75  print ("cfilter called with 4 arguments. Option shiftt0")
76  else:
77  print ("cfilter WARNING: No oldt0s given. The two first arguments must be calibout.txt calibout.txt.")
78  keeprt=sys.argv[3].find("keeprt")>=0
79  keept0=sys.argv[3].find("keept0")>=0
80  shiftrt=sys.argv[3].find("shiftrt")>=0
81  shiftt0=sys.argv[3].find("shiftt0")>=0
82 
83 
84 
85 oldrts=[]
86 oldt0s=[]
87 olderrors=[]
88 for olddb in olddbs:
89  if not olddb.find('#')==0:
90  if len(olddb.split(":")[-1].split())>10:
91  olderrors.append(olddb.strip())
92  elif len(olddb.split(":")[-1].split())>4:
93  oldrts.append(olddb.strip())
94  else:
95  oldt0s.append(olddb.strip())
96 
97 newrts =[]
98 newt0s =[]
99 newerrors =[]
100 for newdb in newdbs:
101  if not newdb.find('#')==0:
102  if len(newdb.split(":")[-1].split())>10:
103  newerrors.append(newdb.strip())
104  elif len(newdb.split(":")[-1].split())>4:
105  newrts.append(newdb.strip())
106  else:
107  newt0s.append(newdb.strip())
108 
109 rts=newrts
110 t0s=newt0s
111 if keeprt: rts=oldrts
112 if keept0: t0s=oldt0s
113 
114 #for rt in rts:
115 # print (rt)
116 
117 #rtfix=[0,0]
118 #rtfix=[1,17.7]
119 rtfix=[1,18]
120 rtconsts={}
121 rtshifts={}
122 nprint = 0
123 for rt in rts:
124  #print (rt.strip())
125  rtkeytokens=rt.strip().split(':')[0].strip().split()
126  rtdbkey="%s %s %s %s %s"%(rtkeytokens[0],rtkeytokens[1],rtkeytokens[2],rtkeytokens[3],rtkeytokens[4])
127  rtdictkey="_%s_%s_%s_%s_%s"%(rtkeytokens[0],rtkeytokens[1],rtkeytokens[2],rtkeytokens[3],rtkeytokens[4])
128 
129  p0=float(rt.strip().split(':')[-1].split()[1].strip())
130  p1=float(rt.strip().split(':')[-1].split()[2].strip())
131  p2=float(rt.strip().split(':')[-1].split()[3].strip())
132  p3=float(rt.strip().split(':')[-1].split()[4].strip())
133 
134  # create polynomial shifted to get the t where r(t)=1
135  # find root (replacing the numpy root finder)
136  poly_orig = [p3, p2, p1, p0 - rtfix[0]]
137  r_orig = p3root(p3, p2, p1, p0 - rtfix[0])
138  shiftval = ( r_orig - rtfix[1] ) # calculate first shift along t so that r(18)=1
139 
140  p0n, p1n, p2n, p3n = tshift_poly(shiftval, p0, p1, p2, p3)
141  poly_new = [p3n, p2n, p1n, p0n]
142  r_new = p3root(p3n, p2n, p1n, p0n)
143 
144  shiftvaln = 0.0 #poly_new.r[2] #calculate second shift along t so that r(0)=0
145  p0nn, p1nn, p2nn, p3nn = tshift_poly(shiftvaln, p0n, p1n, p2n, p3n)
146  poly_newnew = [p3nn, p2nn, p1nn, p0nn]
147  # end replace numpy
148 
149  rtconsts[rtdbkey]=[]
150  if shiftrt:
151  if nprint<10 :
152  print ("%16s ... original polynomial: r=%.1f mm @ %f ns, shifting %f ns in t, new polynomial: r(%.1f ns) = %f mm, dt0 = %f" % (
153  rtdbkey,
154  rtfix[0],
155  r_orig,
156  shiftval,
157  rtfix[1],
158  p3nn*rtfix[1]*rtfix[1]*rtfix[1] + p2nn*rtfix[1]*rtfix[1] + p1nn*rtfix[1] + p0nn,
159  shiftvaln
160  ))
161  print ("original polynomial: ",p0,p1,p2,p3)
162  print ("new polynomial: ",p0nn,p1nn,p2nn,p3nn)
163  nprint = nprint +1
164  #rtconsts[rtdbkey].append(0.0)
165  rtconsts[rtdbkey].append(p0nn)
166  rtconsts[rtdbkey].append(p1nn)
167  rtconsts[rtdbkey].append(p2nn)
168  rtconsts[rtdbkey].append(p3nn)
169  rtshifts[rtkeytokens[0]]=shiftvaln #only compensate t0s for the second shift
170  else:
171  rtconsts[rtdbkey].append(p0)
172  rtconsts[rtdbkey].append(p1)
173  rtconsts[rtdbkey].append(p2)
174  rtconsts[rtdbkey].append(p3)
175  rtshifts[rtkeytokens[0]]=p0
176 
177 t0consts={}
178 for t0 in t0s:
179  t0keytokens=t0.strip().split(':')[0].strip().split()
180  t0dbkey="%s %s %s %s %s"%(t0keytokens[0],t0keytokens[1],t0keytokens[2],t0keytokens[3],t0keytokens[4])
181  t0consts[t0dbkey]=[]
182  t0consts[t0dbkey].append(float(t0.strip().split(':')[-1].split()[0].strip()))
183  t0consts[t0dbkey].append(float(t0.strip().split(':')[-1].split()[1].strip()))
184 
185 
186 #shift rts
187 dbrtouts=[]
188 for part in rtconsts:
189  dbrtouts.append("%s : 0 %e %e %e %e"%(part,rtconsts[part][0],rtconsts[part][1],rtconsts[part][2],rtconsts[part][3]))
190 
191 #shift t0s
192 dbt0outs=[]
193 for straw in t0consts:
194  if not do_global:
195  if shiftt0 and rtshifts.has_key(straw.split()[0]):
196  dbt0outs.append("%s : %f %f"%(straw,t0consts[straw][0]+rtshifts[straw.split()[0]],t0consts[straw][1]))
197  #dbt0outs.append("%s : %f %f SHIFTED %f"%(straw,t0consts[straw][0]+rtshifts[straw.split()[0]],t0consts[straw][1], rtshifts[straw.split()[0]]))
198  else:
199  dbt0outs.append("%s : %f %f"%(straw,t0consts[straw][0],t0consts[straw][1]))
200  else:
201  dbt0outs.append("%s : %f %f"%(straw,t0consts[straw][0]+Dt0,t0consts[straw][1]))
202 
203 
204 dbrtouts.sort()
205 dbt0outs.sort()
206 
207 dbfile=open("dbconst.txt","w")
208 # In case there are not errors:
209 if len (newerrors)==0:
210  dbfile.write("# Fileformat=1\n")
211  dbfile.write("# RtRelation\n")
212  for dbrtout in dbrtouts:
213  dbfile.write(dbrtout + '\n')
214  dbfile.write("# StrawT0\n")
215  for dbt0out in dbt0outs:
216  dbfile.write(dbt0out + '\n')
217  dbfile.write("#GLOBALOFFSET 0.0000\n")
218 elif len (newerrors)>1:
219  dbfile.write("# Fileformat=2\n")
220  dbfile.write("# RtRelation\n")
221  for dbrtout in dbrtouts:
222  dbfile.write(dbrtout + '\n')
223  dbfile.write("# errors\n")
224  for errorout in newerrors:
225  dbfile.write(errorout + '\n')
226  dbfile.write("# StrawT0\n")
227  for dbt0out in dbt0outs:
228  dbfile.write(dbt0out + '\n')
229  dbfile.write("#GLOBALOFFSET 0.0000\n")
230 
231 
232 if not do_global:
233 #if do_global:
234  t0dicts={}
235  for olddict in olddicts:
236  t0dicts[olddict.split()[0]]=float(olddict.split()[1])
237 
238  dictouts=[]
239  for straw in t0dicts:
240  keytokens=straw.split('_')
241  if keytokens[1]=="all": detector='-3'
242  else: detector=keytokens[1].strip()
243  if shiftt0 and rtshifts.has_key(detector):
244  dictouts.append('%s %f'%(straw,t0dicts[straw]+rtshifts[detector]))
245  #dictouts.append('%s %f SHIFTED %f'%(straw,t0dicts[straw]+rtshifts[detector],rtshifts[detector]))
246  else:
247  dictouts.append('%s %f'%(straw,t0dicts[straw]))
248 
249  dictouts.sort()
250 
251  dictfile=open("dictconst.txt","w")
252  for dictout in dictouts:
253  dictfile.write(dictout + '\n')
254 
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
dumpHVPathFromNtuple.append
bool append
Definition: dumpHVPathFromNtuple.py:91
TRTCalib_cfilter.p3root
def p3root(a, b, c, d)
Definition: TRTCalib_cfilter.py:14
TRTCalib_cfilter.tshift_poly
def tshift_poly(dt, p0, p1, p2, p3)
Definition: TRTCalib_cfilter.py:7
Trk::open
@ open
Definition: BinningType.h:40
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Trk::split
@ split
Definition: LayerMaterialProperties.h:38
python.LArMinBiasAlgConfig.float
float
Definition: LArMinBiasAlgConfig.py:65
Generated on Fri Apr 25 2025 21:22:35 for ATLAS Offline Software by   doxygen 1.8.18