CaloMiscalibTools/test/make_miscalibration.py

0001 #!/usr/bin/env python3
0002 import string,time,thread,random,math,sys
0003
0004 #global variables
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 # events/crystal/fb
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     #return a random factor according to  eta
0026     if lumi==0:
0027         if randval==1:
0028             # from CSA07 discussion
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 #            return 1.0
0035         else:
0036             return MINRES
0037     if endcap != 1:
0038         # from AN
0039         #0.000 0.261         6.19 0.12
0040         #0.783 0.957        10.7 0.27
0041         #1.305 1.479        15.0 0.42
0042         # get eta from etaindex
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         # extrapolation
0051
0052         for ieta in range(0,NPOINT):
0053             if eta[ieta]> etacur:
0054                 break
0055         indexmin=ieta-1
0056         indexmax=ieta
0057         # in case there are no more points left
0058         if(indexmin<0):
0059             indexmin=indexmin+1
0060             indexmax=indexmax+1
0061         # now real extrapolation
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             # first time called, we have to derive it from last barrel
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     #return a random factor according to  eta
0083     if lumi==0:
0084         if randval==1:
0085             return random.gauss(1.1,MINRES)
0086 #            return 1.0
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 #main
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     #endcap=0
0147     ecalbarrel=1
0148 elif sys.argv[1]=='ecalendcap':
0149     #endcap=1
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 #if ecalendcap==1:
0169 #    ne_z=[]
0170 #    ne_mod=[]
0171 #    ne_xtal=[]
0172 #    necells=0
0173
0174 if hcal==1:
0175
0176 # read valid hcal cells
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             # skip comment lines
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 # read non existing cells file, not needed anymore
0198 #    necell=open('non_existing_cell_endcap','r')
0199
0200 #    for line in necell:
0201 #        necells=necells+1
0202 #        curr_list=line.split()
0203 #        ne_z.append(string.atoi(curr_list[0]))
0204 #        ne_mod.append(string.atoi(curr_list[1]))
0205 #        ne_xtal.append(string.atoi(curr_list[2]))
0206 #    necell.close()
0207 #    print 'Read ',necells,' non-existing cells for endcap'
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 # now open file
0224 txtfile=open(fileout+'.txt','w')
0225 xmlfile=open(fileout+'.xml','w')
0226 xmlfileinv=open('inv_'+fileout+'.xml','w')
0227 # write header
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     # define ranges
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                 # create line:
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 ############################## HERE HCAL
0293 if hcal==1:
0294 #    mindet=1
0295 #    maxdet=4
0296 #    mineta=-63
0297 #    maxeta=63
0298 #    minphi=0
0299 #    maxphi=127
0300 #    mindepth=1
0301 #    maxdepth=4
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 #    for detindex in range(mindet,maxdet+1):
0320 #        for etaindex in range(mineta,maxeta+1):
0321 #            for phiindex in range(minphi,maxphi+1):
0322 #                for depthindex in range(mindepth,maxdepth+1):
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 # write footer
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 #print miscalibecal(5,0,1,85,1,0)
0357 #print miscalibecal(5,1,-1,10,1,0)
0358 #print miscalibecal(5,1,-1,10,1,1)
0359 #print miscalibecal(5,1,-1,10,1,1)