CaloMiscalibTools/test/make_ecal_miscaplibration.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.02
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 miscalib(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             return random.gauss(1,MINRES)
0030 #            return 1.0
0031         else:
0032             return MINRES
0033     if endcap != 1:
0034         # from AN
0035         #0.000 0.261         6.19 0.12
0036         #0.783 0.957        10.7 0.27
0037         #1.305 1.479        15.0 0.42
0038         # get eta from etaindex
0039         etacur=(etaindex-0.5)*0.0174
0040         eta_ev.append(300./5*lumi)
0041         eta_ev.append(270./5*lumi)
0042         eta_ev.append(170./5*lumi)
0043         res.append(math.sqrt(6.19**2/eta_ev[0]+0.12**2))
0044         res.append(math.sqrt(10.7**2/eta_ev[1]+0.27**2))
0045         res.append(math.sqrt(15**2/eta_ev[2]+0.42**2))
0046         # extrapolation
0047
0048         for ieta in range(0,NPOINT):
0049             if eta[ieta]> etacur:
0050                 break
0051         indexmin=ieta-1
0052         indexmax=ieta
0053         # in case there are no more points left
0054         if(indexmin<0):
0055             indexmin=indexmin+1
0056             indexmax=indexmax+1
0057         # now real extrapolation
0058         real_res=res[indexmin]+(etacur-eta[indexmin])*(res[indexmax]-res[indexmin])/(eta[indexmax]-eta[indexmin])
0059         real_res=real_res/100
0060         if real_res>MINRES:
0061             real_res=MINRES
0062         if randval==1:
0063             return random.gauss(1,real_res)
0064         else:
0065             return real_res
0066     if endcap:
0067         if ENDCAPSINGLE==0:
0068             # first time called, we have to derive it from last barrel
0069             ENDCAPSINGLE=miscalib(lumi,0,1,85,1,0)
0070         if randval==1:
0071             return random.gauss(1,ENDCAPSINGLE)
0072         else:
0073             return ENDCAPSINGLE
0074
0075
0076
0077 #main
0078 if len(sys.argv)==1:
0079     print('Usage: '+sys.argv[0]+' <barrel|endcap> <lumi> <filename> [MINRES=0.02] [SEED=random]')
0080     print('       put lumi=0 for precalibration values (random at MINRES)')
0081     sys.exit(1)
0082
0083 if sys.argv[1]=='barrel':
0084     endcap=0
0085 elif sys.argv[1]=='endcap':
0086     endcap=1
0087 else:
0088     print('please specify one of <barrel|endcap>')
0089     sys.exit(1)
0090
0091 if endcap==1:
0092     # read non existing cells file
0093     ne_z=[]
0094     ne_mod=[]
0095     ne_xtal=[]
0096     necells=0
0097
0098 #    necell=open('non_existing_cell_endcap','r')
0099
0100 #    for line in necell:
0101 #        necells=necells+1
0102 #        curr_list=line.split()
0103 #        ne_z.append(string.atoi(curr_list[0]))
0104 #        ne_mod.append(string.atoi(curr_list[1]))
0105 #        ne_xtal.append(string.atoi(curr_list[2]))
0106 #    necell.close()
0107 #    print 'Read ',necells,' non-existing cells for endcap'
0108
0109
0110 lumi=string.atof(sys.argv[2])
0111
0112 if lumi<0:
0113     print('lumi = '+str(lumi)+' not valid')
0114     sys.exit(1)
0115
0116 fileout=sys.argv[3]
0117
0118 if len(sys.argv)>=5:
0119     MINRES=string.atof(sys.argv[4])
0120
0121 print('Using minimal resolution: '+str(MINRES))
0122
0123 if len(sys.argv)>=6:
0124     SEED=string.atoi(sys.argv[5])
0125     print('Using fixed seed for random generation: '+str(SEED))
0126     random.seed(SEED)
0127
0128 # now open file
0129 xmlfile=open(fileout,'w')
0130 # write header
0131 xmlfile.write('<?xml version="1.0" ?>\n')
0132 xmlfile.write('\n')
0133 xmlfile.write('<CalibrationConstants>\n')
0134 if endcap==0:
0135     xmlfile.write('  <EcalBarrel>\n')
0136 else:
0137     xmlfile.write('  <EcalEndcap>\n')
0138 # define ranges
0139 mineta=1
0140 minphi=1
0141 if endcap==0:
0142     maxeta=85
0143     maxphi=360
0144 else:
0145     maxeta=100
0146     maxphi=100
0147
0148 count=0
0149 for zindex in (-1,1):
0150     for etaindex in range(mineta,maxeta+1):
0151         for phiindex in range(minphi,maxphi+1):
0152             miscal_fac=miscalib(lumi,endcap,zindex,etaindex,phiindex,1)
0153             # create line:
0154             if endcap==0:
0155            if zindex==-1:
0156                   line='        <Cell eta_index="'+str(-etaindex)+'" phi_index="'+str(phiindex)+'" scale_factor="'+str(miscal_fac)+'"/>\n'
0157                   xmlfile.write(line)
0158                   count=count+1
0159            else:
0160                   line='        <Cell eta_index="'+str(+etaindex)+'" phi_index="'+str(phiindex)+'" scale_factor="'+str(miscal_fac)+'"/>\n'
0161                   xmlfile.write(line)
0162                   count=count+1
0163             else:
0164                 goodxtal=1
0165                 # check if it exists
0166 #                for bad in range(0,necells):
0167 #                    if ne_xtal[bad]==phiindex:
0168 #                        if ne_mod[bad]==etaindex:
0169 #                            if ne_z[bad]==zindex:
0170 #                                goodxtal=0
0171 #                                break
0172 #                if goodxtal==1:
0173 #          if zindex==-1:
0174
0175 #                line='        <Cell module_index="'+str(etaindex)+'" crystal_index="'+str(phiindex)+'"  scale_factor="'+str(miscal_fac)+'"/>\n'
0176                 line='        <Cell x_index="'+str(etaindex)+'" y_index="'+str(phiindex)+'" z_index="'+str(zindex)+'" scale_factor="'+str(miscal_fac)+'"/>\n'
0177
0178                 xmlfile.write(line)
0179                 count=count+1
0180                 #                   else :
0181                 #                      line='        <Cell module_index="'+str(etaindex)+'" crystal_index="'+str(phiindex)+'"  scale_factor="'+str(miscal_fac)+'"/>\n'
0182                 #             xmlfile.write(line)
0183                 #                      count=count+1
0184
0185 # write footer
0186 if endcap==0:
0187     xmlfile.write('  </EcalBarrel>\n')
0188 else:
0189     xmlfile.write('  </EcalEndcap>\n')
0190 xmlfile.write('</CalibrationConstants>\n')
0191 xmlfile.close()
0192 print('File '+fileout+' written with '+str(count)+' lines')
0193 #print miscalib(5,0,1,85,1,0)
0194 #print miscalib(5,1,-1,10,1,0)
0195 #print miscalib(5,1,-1,10,1,1)
0196 #print miscalib(5,1,-1,10,1,1)