File indexing completed on 2024-11-25 02:29:07
0001
0002 import string,time,thread,random,math,sys
0003
0004
0005 MINRES=0.04
0006 SEED=-1
0007 ENDCAPSINGLE=0
0008 NPOINT=3
0009 eta=[]
0010 eta_ev=[]
0011 res=[]
0012
0013 eta.append((0.+0.261)/2)
0014 eta.append((0.957+0.783)/2)
0015 eta.append((1.479+1.305)/2)
0016
0017
0018
0019
0020 def miscalibecal(lumi,endcap,z,etaindex,phiindex,randval):
0021 global MINRES,ENDCAPSINGLE,NPOINT
0022 if endcap:
0023 modindex=etaindex
0024 crysindex=phiindex
0025
0026 if lumi==0:
0027 if randval==1:
0028
0029 if endcap != 1:
0030 gauss_smear=0.014+(etaindex-1)*(0.022-0.014)/(85.-1.)
0031 else:
0032 gauss_smear=MINRES
0033 return random.gauss(1,gauss_smear)
0034
0035 else:
0036 return MINRES
0037 if endcap != 1:
0038
0039
0040
0041
0042
0043 etacur=(etaindex-0.5)*0.0174
0044 eta_ev.append(300./5*lumi)
0045 eta_ev.append(270./5*lumi)
0046 eta_ev.append(170./5*lumi)
0047 res.append(math.sqrt(6.19**2/eta_ev[0]+0.12**2))
0048 res.append(math.sqrt(10.7**2/eta_ev[1]+0.27**2))
0049 res.append(math.sqrt(15**2/eta_ev[2]+0.42**2))
0050
0051
0052 for ieta in range(0,NPOINT):
0053 if eta[ieta]> etacur:
0054 break
0055 indexmin=ieta-1
0056 indexmax=ieta
0057
0058 if(indexmin<0):
0059 indexmin=indexmin+1
0060 indexmax=indexmax+1
0061
0062 real_res=res[indexmin]+(etacur-eta[indexmin])*(res[indexmax]-res[indexmin])/(eta[indexmax]-eta[indexmin])
0063 real_res=real_res/100
0064 if real_res>MINRES:
0065 real_res=MINRES
0066 if randval==1:
0067 return random.gauss(1,real_res)
0068 else:
0069 return real_res
0070 if endcap:
0071 if ENDCAPSINGLE==0:
0072
0073 ENDCAPSINGLE=miscalibecal(lumi,0,1,85,1,0)
0074 if randval==1:
0075 return random.gauss(1,ENDCAPSINGLE)
0076 else:
0077 return ENDCAPSINGLE
0078
0079
0080 def miscalibhcal(lumi,det,etaindex,phiindex,depth,randval):
0081 global MINRES,NPOINT
0082
0083 if lumi==0:
0084 if randval==1:
0085 return random.gauss(1.1,MINRES)
0086
0087 else:
0088 return MINRES
0089 elif lumi==1:
0090 if abs(etaindex)<=11:
0091 if randval==1:
0092 return random.gauss(1,0.04)
0093 else:
0094 return 0.04
0095 elif abs(etaindex)<=16:
0096 if randval==1:
0097 return random.gauss(1,0.03)
0098 else:
0099 return 0.03
0100 else:
0101 if randval==1:
0102 return random.gauss(1,0.02)
0103 else:
0104 return 0.02
0105 elif lumi==10:
0106 if abs(etaindex)<=5:
0107 if randval==1:
0108 return random.gauss(1,0.03)
0109 else:
0110 return 0.03
0111 if abs(etaindex)<=11:
0112 if randval==1:
0113 return random.gauss(1,0.025)
0114 else:
0115 return 0.025
0116 elif abs(etaindex)<=16:
0117 if randval==1:
0118 return random.gauss(1,0.014)
0119 else:
0120 return 0.014
0121 elif abs(etaindex)<=20:
0122 if randval==1:
0123 return random.gauss(1,0.01)
0124 else:
0125 return 0.01
0126 else:
0127 if randval==1:
0128 return random.gauss(1,0.02)
0129 else:
0130 return 0.02
0131 else:
0132 print('lumi must be 0, 1 or 10')
0133
0134
0135
0136 if len(sys.argv)==1:
0137 print('Usage: '+sys.argv[0]+' <ecalbarrel|ecalendcap|HCAL> <lumi> <filename> <validcell for hcal> [MINRES=0.02] [SEED=random]')
0138 print(' put lumi=0 for precalibration values (CSA07)')
0139 sys.exit(1)
0140
0141 ecalbarrel=0
0142 ecalendcap=0
0143 hcal=0
0144
0145 if sys.argv[1]=='ecalbarrel':
0146
0147 ecalbarrel=1
0148 elif sys.argv[1]=='ecalendcap':
0149
0150 ecalendcap=1
0151 elif sys.argv[1]=='HCAL':
0152 hcal=1
0153 else:
0154 print('please specify one of <barrel|endcap|HCAL>')
0155 sys.exit(1)
0156
0157 lumi=string.atof(sys.argv[2])
0158
0159 if lumi<0:
0160 print('lumi = '+str(lumi)+' not valid')
0161 sys.exit(1)
0162
0163 fileout=sys.argv[3]
0164
0165 if len(sys.argv)>=5:
0166 validcell=sys.argv[4]
0167
0168
0169
0170
0171
0172
0173
0174 if hcal==1:
0175
0176
0177 hcalcell=open(validcell,'r')
0178 hcal_det=[]
0179 hcal_eta=[]
0180 hcal_phi=[]
0181 hcal_depth=[]
0182 hcalcells=0
0183
0184 for line in hcalcell:
0185 curr_list=line.split()
0186 if curr_list[0]=="#":
0187
0188 continue
0189 hcal_eta.append(string.atoi(curr_list[0]))
0190 hcal_phi.append(string.atoi(curr_list[1]))
0191 hcal_depth.append(string.atoi(curr_list[2]))
0192 hcal_det.append(curr_list[3])
0193 hcalcells=hcalcells+1
0194
0195 hcalcell.close()
0196 print('Read ',hcalcells,' valid cells for hcal')
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212 if len(sys.argv)>=6:
0213 MINRES=string.atof(sys.argv[5])
0214
0215 if (hcal==0) or (lumi==0):
0216 print('Using minimal resolution: '+str(MINRES))
0217
0218 if len(sys.argv)>=7:
0219 SEED=string.atoi(sys.argv[6])
0220 print('Using fixed seed for random generation: '+str(SEED))
0221 random.seed(SEED)
0222
0223
0224 txtfile=open(fileout+'.txt','w')
0225 xmlfile=open(fileout+'.xml','w')
0226 xmlfileinv=open('inv_'+fileout+'.xml','w')
0227
0228 xmlfile.write('<?xml version="1.0" ?>\n')
0229 xmlfile.write('\n')
0230 xmlfile.write('<CalibrationConstants>\n')
0231
0232 xmlfileinv.write('<?xml version="1.0" ?>\n')
0233 xmlfileinv.write('\n')
0234 xmlfileinv.write('<CalibrationConstants>\n')
0235
0236 if ecalbarrel==1:
0237 xmlfile.write(' <EcalBarrel>\n')
0238 xmlfileinv.write(' <EcalBarrel>\n')
0239 elif ecalendcap==1:
0240 xmlfile.write(' <EcalEndcap>\n')
0241 xmlfileinv.write(' <EcalEndcap>\n')
0242 elif hcal==1:
0243 xmlfile.write(' <Hcal>\n')
0244 xmlfileinv.write(' <Hcal>\n')
0245
0246
0247 count=0
0248 if ecalbarrel==1 or ecalendcap==1:
0249
0250
0251 mineta=1
0252 minphi=1
0253 if ecalbarrel==1:
0254 maxeta=85
0255 maxphi=360
0256 else:
0257 maxeta=100
0258 maxphi=100
0259
0260 for zindex in (-1,1):
0261 for etaindex in range(mineta,maxeta+1):
0262 for phiindex in range(minphi,maxphi+1):
0263 miscal_fac=miscalibecal(lumi,ecalendcap,zindex,etaindex,phiindex,1)
0264
0265 if ecalbarrel==1:
0266 if zindex==-1:
0267 line=' <Cell eta_index="'+str(-etaindex)+'" phi_index="'+str(phiindex)+'" scale_factor="'+str(miscal_fac)+'"/>\n'
0268 lineinv=' <Cell eta_index="'+str(-etaindex)+'" phi_index="'+str(phiindex)+'" scale_factor="'+str(1./miscal_fac)+'"/>\n'
0269 txtline=str(etaindex)+' '+str(phiindex)+' '+str(miscal_fac)+'\n'
0270 xmlfile.write(line)
0271 xmlfileinv.write(lineinv)
0272 txtfile.write(txtline)
0273 count=count+1
0274 else:
0275 line=' <Cell eta_index="'+str(+etaindex)+'" phi_index="'+str(phiindex)+'" scale_factor="'+str(miscal_fac)+'"/>\n'
0276 lineinv=' <Cell eta_index="'+str(+etaindex)+'" phi_index="'+str(phiindex)+'" scale_factor="'+str(1./miscal_fac)+'"/>\n'
0277 txtline=str(-etaindex)+' '+str(phiindex)+' '+str(miscal_fac)+'\n'
0278 xmlfile.write(line)
0279 xmlfileinv.write(lineinv)
0280 txtfile.write(txtline)
0281 count=count+1
0282 else:
0283 goodxtal=1
0284 line=' <Cell x_index="'+str(etaindex)+'" y_index="'+str(phiindex)+'" z_index="'+str(zindex)+'" scale_factor="'+str(miscal_fac)+'"/>\n'
0285 lineinv=' <Cell x_index="'+str(etaindex)+'" y_index="'+str(phiindex)+'" z_index="'+str(zindex)+'" scale_factor="'+str(1./miscal_fac)+'"/>\n'
0286 txtline=str(etaindex)+' '+str(phiindex)+' '+str(zindex)+' '+str(miscal_fac)+'\n'
0287 xmlfile.write(line)
0288 xmlfileinv.write(lineinv)
0289 txtfile.write(txtline)
0290 count=count+1
0291
0292
0293 if hcal==1:
0294
0295
0296
0297
0298
0299
0300
0301
0302 txtline='# eta phi dep det value \n'
0303 txtfile.write(txtline)
0304 for cell in range(0,hcalcells):
0305 etaindex=hcal_eta[cell]
0306 phiindex=hcal_phi[cell]
0307 depthindex=hcal_depth[cell]
0308 detindex=-1
0309 if(hcal_det[cell]=="HB"):
0310 detindex=1;
0311 if(hcal_det[cell]=="HE"):
0312 detindex=2;
0313 if(hcal_det[cell]=="HO"):
0314 detindex=3;
0315 if(hcal_det[cell]=="HF"):
0316 detindex=4;
0317
0318 if(detindex>0):
0319
0320
0321
0322
0323 miscal_fac=miscalibhcal(lumi,detindex,etaindex,phiindex,depthindex,1)
0324 line=' <Cell det_index="'+str(detindex)+'" eta_index="'+str(etaindex)+'" phi_index="'+str(phiindex)+'" depth_index="'+str(depthindex)+'" scale_factor="'+str(miscal_fac)+'"/>\n'
0325 lineinv=' <Cell det_index="'+str(detindex)+'" eta_index="'+str(etaindex)+'" phi_index="'+str(phiindex)+'" depth_index="'+str(depthindex)+'" scale_factor="'+str(1./miscal_fac)+'"/>\n'
0326 xmlfile.write(line)
0327 xmlfileinv.write(lineinv)
0328 txtline=" "+str(etaindex)+" "+str(phiindex)+" "+str(depthindex)+" "+hcal_det[cell]+" "+str(miscal_fac)+"\n"
0329 txtfile.write(txtline)
0330 count=count+1
0331 else:
0332 txtline=" "+str(etaindex)+" "+str(phiindex)+" "+str(depthindex)+" "+hcal_det[cell]+" 1.0\n"
0333 txtfile.write(txtline)
0334
0335
0336
0337
0338 if ecalbarrel==1:
0339 xmlfile.write(' </EcalBarrel>\n')
0340 xmlfileinv.write(' </EcalBarrel>\n')
0341 elif ecalendcap==1:
0342 xmlfile.write(' </EcalEndcap>\n')
0343 xmlfileinv.write(' </EcalEndcap>\n')
0344 elif hcal==1:
0345 xmlfile.write(' </Hcal>\n')
0346 xmlfileinv.write(' </Hcal>\n')
0347
0348 xmlfile.write('</CalibrationConstants>\n')
0349 xmlfileinv.write('</CalibrationConstants>\n')
0350 xmlfile.close()
0351 xmlfileinv.close()
0352 txtfile.close()
0353 print('File '+fileout+'.xml'+' written with '+str(count)+' lines')
0354 print('File inv_'+fileout+'.xml'+' written with '+str(count)+' lines')
0355 print('File '+fileout+'.txt'+' written with '+str(count)+' lines')
0356
0357
0358
0359