Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-25 02:29:07

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