CaloMiscalibTools/test/make_miscalibration.py

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