ATLAS Offline Software
Loading...
Searching...
No Matches
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
5import os,sys, math
6
7def 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
14def 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
44keeprt=False
45keept0=False
46shiftrt=False
47shiftt0=False
48do_global=False
49
50if 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()
58else:
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
85oldrts=[]
86oldt0s=[]
87olderrors=[]
88for 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
97newrts =[]
98newt0s =[]
99newerrors =[]
100for 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
109rts=newrts
110t0s=newt0s
111if keeprt: rts=oldrts
112if keept0: t0s=oldt0s
113
114#for rt in rts:
115# print (rt)
116
117#rtfix=[0,0]
118#rtfix=[1,17.7]
119rtfix=[1,18]
120rtconsts={}
121rtshifts={}
122nprint = 0
123for 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
177t0consts={}
178for 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
187dbrtouts=[]
188for 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
192dbt0outs=[]
193for 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
204dbrtouts.sort()
205dbt0outs.sort()
206
207dbfile=open("dbconst.txt","w")
208# In case there are not errors:
209if 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")
218elif 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
232if 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
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15
find
std::string find(const std::string &s)
return a remapped string
Definition hcg.cxx:138
split
std::vector< std::string > split(const std::string &s, const std::string &t=":")
Definition hcg.cxx:177
TRTCalib_cfilter.p3root
p3root(a, b, c, d)
Definition TRTCalib_cfilter.py:14
TRTCalib_cfilter.tshift_poly
tshift_poly(dt, p0, p1, p2, p3)
Definition TRTCalib_cfilter.py:7
Generated on for ATLAS Offline Software by doxygen 1.14.0