drich-dev/scripts/RequiredResolution.py

0001 from math import *
0002 import matplotlib.pyplot as plt
0003 import numpy as np
0004 import getopt
0005 import sys
0006 from matplotlib.ticker import *
0007
0008 narg = len(sys.argv)
0009 if narg == 1:
0010   print("usage :",sys.argv[0],"1 1");
0011   print("First number 0 --> pi/K 1--> e/pi");
0012   print("Second number 0 --> gas 1--> aerogel");
0013   exit();
0014 agrv1 = sys.argv[1]
0015 agrv2 = sys.argv[2]
0016
0017 file = "";
0018 mass1 = 0;
0019 mass2 = 0;
0020 refIndex =0;
0021 maxmom =0;
0022
0023 fig, ax = plt.subplots();
0024 ml = AutoMinorLocator(2);
0025 #plt.axes().xaxis.set_minor_locator(ml);
0026 ax.xaxis.set_minor_locator(ml);
0027
0028 if agrv1 == '0':
0029   mass1 = 0.139;
0030   mass2 =  0.493;
0031   file = "piK.pdf";
0032   plt.title('pi/K separation');
0033   if agrv2 == '0':
0034     refIndex = 1.0008;
0035     maxmom = 60;
0036     plt.ylim([0.1,20]);
0037   elif agrv2 == '1':
0038     refIndex = 1.02;
0039     maxmom = 20;
0040     plt.ylim([0.1,30]);
0041
0042 elif agrv1 == '1':
0043   mass1 = 0.000511;
0044   mass2 =  0.139;
0045   file = "epi.pdf";
0046   plt.title("e-Pi Separation");
0047   if agrv2 == '0':
0048     refIndex = 1.0008;
0049     maxmom = 20;
0050     plt.ylim([0.1,100]);
0051   elif agrv2 == '1':
0052     refIndex = 1.02;
0053     maxmom =10;
0054     plt.ylim([0.1,100]);
0055
0056 p=[0.0];
0057 y=[0.0];
0058 z=[0.0];
0059 m=[0.0];
0060 #refIndex = 1.0008;
0061 #mass1=[0.139,0.000511];
0062 #mass2=[0.493,0.139];
0063
0064 print(mass1,mass2,refIndex,maxmom);
0065
0066 def computebeta(zz,mm):
0067   return zz/sqrt(zz*zz+mm*mm);
0068
0069 def computetheta(beta):
0070   csth = (1/(refIndex*beta));
0071   if(csth>1) : return 0;
0072   else : return acos(csth);
0073
0074 for x in range(0,maxmom,1):
0075   p.append(x);
0076   y.append(1000*computetheta(computebeta((x+0.1),mass1)));
0077   z.append(1000*computetheta(computebeta((x+0.1),mass2)));
0078
0079 Y = np.array(y);
0080 Z = np.array(z);
0081 m_a =np.subtract(Y,Z);
0082 m_0 = np.divide(m_a,3);
0083 m_1 = np.divide(m_a,3.25);
0084 m_2 = np.divide(m_a,3.50);
0085 m_3 = np.divide(m_a,3.75);
0086 m_4 = np.divide(m_a,4);
0087 m = list(m_0);
0088 #print(p);
0089 #print(y);
0090 #plt.plot(p, y, color='green', linestyle='solid', linewidth = 1.5,
0091 #         marker='o', markerfacecolor='blue', markersize=0.5);
0092 plt.plot(p, m, color='magenta', linestyle='solid', linewidth = 1.5,
0093          marker='o', markerfacecolor='blue', markersize=1.0, label="3");
0094 del m;
0095 m = list(m_1);
0096 plt.plot(p, m, color='purple', linestyle='solid', linewidth = 1.5,
0097          marker='o', markerfacecolor='blue', markersize=1.0,label="3.25");
0098 del m;
0099 m = list(m_3);
0100 plt.plot(p, m, color='blue', linestyle='solid', linewidth = 1.5,
0101          marker='o', markerfacecolor='blue', markersize=1.0,label="3.5");
0102 del m;
0103 m = list(m_3);
0104 plt.plot(p, m, color='cyan', linestyle='solid', linewidth = 1.5,
0105          marker='o', markerfacecolor='blue', markersize=1.0,label="3.75");
0106 del m;
0107 m = list(m_4);
0108 plt.plot(p, m, color='green', linestyle='solid', linewidth = 1.5,
0109          marker='o', markerfacecolor='blue', markersize=1.0,label="4");
0110
0111
0112 plt.yscale('log');
0113 plt.xlabel('momentum (GeV/c)');
0114 plt.ylabel('required resolution (mrad)');
0115 plt.legend();
0116 plt.grid(True);
0117 plt.grid(which='minor', alpha=0.3)
0118 #plt.show();
0119 #fig, ax = plt.subplots();
0120 #ml = AutoMinorLocator(2);
0121 #plt.axes().xaxis.set_minor_locator(ml);
0122 #ax.xaxis.set_minor_locator(ml);
0123 plt.savefig(file);
0124