// requires ROOT: http://root.cern.ch/ // // usage: // within root: .x make_pcd_plots+ ("ltt" , "transfer tube length (cm)" , 1 , "technote_parameters.txt" , 2.5 , 3.5 , 260.0 , 30.0) // batch mode: root -b -l -q 'make_pcd_plots.macro+ ("ltt","length of transfer tube (cm)",1,"technote_parameters.txt",2.5,3.5,260.0,30.0)' // // !!!!!IMPORTANT!!!!!: the "+" is very imporant: it tells ROOT to compile the code first. // this saves orders of magnitude in time for the calculations // // name: Polarization Dynamics in a Two Chambered Cell as a Function of Pumping Chamber Diameter Plotter // by: Jaideep Singh (singhj AT jlab DOT org) // date: Janurary 17, 2007 // // description: plots various parameters as a function of pumping chamber diameter // plot parameters are passed through as arguments // requires a parameters file to specify experimental conditions // PG-[num] refers to the "Polarization Gradients..." technote equation number [num] // // references: "Polarization Gradients in A Two Chambered Cell" // http://www.jlab.org/~singhj/docs/polgrad137.ps.gz // "Some Considerations Regarding Pumping Chamber Diameter" // http://www.jlab.org/~singhj/docs/pc_size_note_v1.ps.gz // // description type default // arguments: (1) variable to plot string "ltt" // (2) y axis title string "transfer tube length (cm)" // (3) number of pc operating integer 1 // temperatures to plot // (4) parameter list filename string "technote_parameters.txt" // (5) min pc diameter (inches) double 2.5 // (6) max pc diameter (inches) double 3.5 // (7) mid pc operating temp (C) double 260.0 // (7) delta temp (C) double 30.0 // // if the number of temperatures to plot is 1, then only the middle temperature will be plotted // if the number of temperatures to plot is not 1, then the middle temperature and (middle +/- delta) temps will be plotted // for example: num = 3, mid = 260 , delta = 30 => T_pc = 230,260,290 will be plotted // num = 1, mid = 260 , delta = 30 => T_pc = 260 will be plotted // // format of parameter file list: parameter typical value and units // alkali polarization 0.75 - // maximum X factor when S/V = 1 cm^-1 0.4 - // pumping chamber center to target chamber center 14.2875 cm // inner transfer tube area 0.5 cm^2 // inner target chamber area 2.0 cm^2 // inner target chamber length 40.0 cm // target chamber temperature 40 Celsius // operating target chamber density 11 amg // K mole fraction in alloy 0.900 - // Rb mole fraction in alloy 0.045 - // Na mole fraction in alloy 0.0 - // beam current 15 microamps // beam energy 6000 MeV // room temperature 23 Celsius // wall relaxation time constant 75 hrs // N2 to He density ratio 0.01 - // // EXAMPLE: // ---------------------------- // technote_parameters.txt // ---------------------------- // 0.75 // 0.4 // 14.2875 // 0.5 // 2.0 // 40.0 // 40.0 // 11.0 // 0.90 // 0.045 // 0.0 // 15.0 // 6000.0 // 23.0 // 75.0 // 0.01 // ---------------------------- // // command: root -b -l -q 'make_pcd_plots.macro+ ("ltt","length of transfer tube (cm)",1,"technote_parameters.txt",2.5,3.5,260.0,30.0)' // // output: // root [0] // Processing make_pcd_plots.macro+... // Info in : creating shared library (stuff)./make_pcd_plots_macro.so // START // Reading Data from: technote_parameters.txt // Read 16 parameters and closed file: technote_parameters.txt // PA = 0.75 // X_max = 0.4 // center_2_center (cm) = 14.2875 // att (cm^2) = 0.5 // atc (cm^2) = 2 // ltc (cm) = 40 // ttc (K) = 313.15 // ntc (amg) = 11 // mf_k = 0.9 // mf_rb = 0.045 // mf_na = 0 // I (microamps) = 15 // ebeam (MeV) = 6000 // troom (K) = 296.15 // gwall (1/hrs) = 0.0133333 // n2_to_he = 0.01 // Found ymin = 9.43462 and ymax = 10.7046 // done drawing plot frame! // did only 260 C // Info in : eps file ltt.eps has been created // Real time 0:00:00, CP time 0.160 // END // outputs: an eps with the name "[par].eps" where [par] is the name of the parameter plotted vs pumping chamber diameter // dashed blue line = middle temp - deltaT // solid black line = middle temperature // dotted red lin = middle temp + deltaT // // head files needed to compile macro #include "TComplex.h" #include "TMath.h" #include "TString.h" #include "TRegexp.h" #include "TCanvas.h" #include "TPad.h" #include "TF1.h" #include "TGaxis.h" #include "TH2F.h" #include "TLatex.h" #include "TLine.h" #include "TStyle.h" #include "Riostream.h" #include "TStopwatch.h" // calculates mass of alkali metal in micrograms Double_t falkmass(TString element , Double_t nd , Double_t vol) { Double_t mass; if ((element != "K") && (element != "Na")) { // "Rb" or anything else mass = 85.4678; } else if (element = "K") { // "K" only mass = 39.0983; } else { // "Na" only mass = 22.989770; } Double_t NA = 6.0221415E+23; return 1000000.0*(nd*(1.0E+14)*vol/NA)*mass; } // calculates pure alkali vapor density from 1995 CRC vapor pressure constants // return denisty in 10^14 parts per cm^3 Double_t falknd(TString element , Double_t temp) { Double_t A = 21.45; Double_t B = 9302.0; if ((element != "K") && (element != "Na")) { // "Rb" or anything else if (temp < 312.46) { A = 22.71; B = 9705.0; } } else if (element = "K") { // "K" only A = 21.66; B = 10253.0; if (temp < 336.53) { A = 22.95; B = 10698.0; } } else { // "Na" only A = 22.36; B = 12381.0; if (temp < 370.87) { A = 23.73; B = 12901.0; } } Double_t vp = TMath::Exp(A-B/temp); Double_t R = 8.314472; Double_t NA = 6.0221415E+23; Double_t nd = vp/R/temp; nd = nd*NA; nd = nd/100.0/100.0/100.0; nd = nd/1.0E+14; return nd; } // approximate mean number of nuclear spins flips due to atomic ions // approximation valid for helium densities from 7 to 14 amagats and N2 to He density ratios of 1/3% to 3% //PG-74 Double_t fna_fit(Double_t h , Double_t p) { Double_t p2 = TMath::Power(p,2); Double_t p3 = TMath::Power(p,3); Double_t h2 = TMath::Power(h,2); Double_t h3 = TMath::Power(h,3); Double_t c0 = ( +1.02829 ) + ( +0.22063 )*p + ( -0.054940 )*p2 + ( +0.004915914 )*p3; Double_t c1 = ( +0.19515 ) + ( -0.90280 )*p + ( +0.163844 )*p2 + ( -0.012998767 )*p3; Double_t c2 = ( -0.98146 ) + ( +0.76918 )*p + ( -0.118209 )*p2 + ( +0.008008600 )*p3; Double_t c3 = ( +0.35643 ) + ( -0.19663 )*p + ( +0.025395 )*p2 + ( -0.001295989 )*p3; return c0 + c1*h + c2*h2 + c3*h3; } //PG-75 Double_t fh(Double_t he) { return he/10.0; } Double_t fp(Double_t he , Double_t n2) { return 100*n2/he; } // full calculation of mean number of nuclear spins flips due to atomic ions // from Bonin et al PRA 38, pp4481-7 //PG-73 Double_t fna_full(Double_t r , Double_t Omega) { Double_t Omega2 = TMath::Power(Omega,2.0); Double_t Omega4 = TMath::Power(Omega,4.0); Double_t Q = (4.0+9.0*Omega2)/108.0; Double_t R = Omega*TMath::Sqrt(8.0-13.0*Omega2+16.0*Omega4)/12.0/TMath::Sqrt(3.0); Double_t S = TMath::Power(Q+R,1.0/3.0); Double_t T = TMath::Power(TMath::Abs(Q-R),1.0/3.0); if (Q-R < 0.0) { T = -T; } Double_t gamma1_real = S+T-2.0/3.0; Double_t gamma1_imag = 0.0; Double_t gamma2_real = (-1.0/2.0)*(S+T)-2.0/3.0; Double_t gamma2_imag = (TMath::Sqrt(3.0)/2.0)*(S-T); TComplex gamma1(gamma1_real,gamma1_imag); TComplex gamma2(gamma2_real,gamma2_imag); TComplex gamma2_cc = TComplex::Conjugate(gamma2); Double_t gamma2_2 = TMath::Power(TComplex::Abs(gamma2),2.0); Double_t denom = gamma2_2 + 2.0*gamma1_real*(1.0+gamma1_real); Double_t a1 = (-gamma2_2 + Omega2/2.0)/denom; TComplex denom_c(denom,0.0); TComplex Omega2_c(Omega2,0.0); TComplex two(2.0,0.0); TComplex a2 = (Omega2_c-two*gamma1*gamma2_cc)*(gamma1-gamma2_cc)/(gamma2_cc-gamma2)/denom_c; Double_t na1 = 1.0 + a1/(1-r*gamma1_real); TComplex r_c(r,0.0); TComplex na2_c = a2/(TComplex::One()-r_c*gamma2); Double_t na2_real = na2_c.Re(); return na1 + na2_real; } //PG-73 Double_t fomega(Double_t he) { Double_t kex = 15.0E9; Double_t Aah = 8.66E9; Double_t tex = 1.0/(kex*he); return 2.0*TMath::Pi()*Aah*tex; } //PG-73 Double_t fr(Double_t he , Double_t n2) { Double_t kn2 = 27.0E9; Double_t km = 0.060E9; Double_t kex = 15.0E9; // Double_t Aah = 8.66E9; Double_t ta = 1.0/(kn2*n2+km*TMath::Power(he,2.0)); Double_t tex = 1.0/(kex*he); return ta/tex; } Double_t fna(Double_t he , Double_t n2 , Double_t choice) { Double_t r = fr(he,n2); Double_t Omega = fomega(he); if (choice >= 0.0) { return fna_full(r,Omega); } else { return fna_fit(r,Omega); } } // He-He dipolar "self" relaxation // functional fit to Newbury et al, PRA 48, pp4411-20 //PG-35,36 Double_t fgdip(Double_t he , Double_t temp) { Double_t c0 = +1.2319E+0; Double_t c1 = +2.8591E-1; Double_t c2 = -2.1793E-1; Double_t c3 = -1.4426E-2; Double_t c4 = +5.3315E-1; Double_t c5 = +1.2376E+3; Double_t T0 = +296.15; Double_t t = temp/T0; return he/(744.0*(c0*TMath::Power(t,c1) + c2 + c3*t + c4/(1+c5*t))); } //PG-94 Double_t fups(Double_t tpc , Double_t ttc) { Double_t t = tpc/ttc; Double_t top = t - 1.0; Double_t bot = TMath::Power(t,0.3) - 1.0; Double_t right = TMath::Power(ttc/273.15,0.7); return 0.3*top/bot*right; } // diffusion rate into/out of target chamber in 1/hrs //PG-93 Double_t fdtc(Double_t att , Double_t ltt , Double_t vtc , Double_t ntc , Double_t t , Double_t ttc) { Double_t tpc = t*ttc; Double_t ups = fups(tpc,ttc); return 6488.21*att*ups/vtc/ltt/ntc; } // diffusion rate into/out of target chamber in 1/hrs //PG-89 Double_t fdpc(Double_t att , Double_t ltt , Double_t vtc , Double_t ntc , Double_t t , Double_t ttc, Double_t v) { Double_t dtc = fdtc(att , ltt , vtc , ntc , t , ttc); return (t/v)*dtc; } // all inputs in 10^14 cm^-3 // alkali spin exchange rate in 1/hrs //PG TABLE 1 Double_t fgse(Double_t rb , Double_t k , Double_t na) { return rb/40.9 + k/45.5 + na/45.5; } // pumping chamber relaxation rate in 1/hrs Double_t fgpc(Double_t gse , Double_t X , Double_t gwall , Double_t gdip) { return gse*X + gwall + gdip; } // target chamber relaxation rate in 1/hrs Double_t fgtc(Double_t gse , Double_t X , Double_t gwall , Double_t gdip , Double_t gbeam) { return gse*X + gwall + gdip + gbeam; } // transfer tube spin exchange rate in 1/hrs //PG-120 Double_t fgsett(Double_t tpc , Double_t rb_pc , Double_t k_pc , Double_t na_pc , Double_t ttc , Double_t rb_tc , Double_t k_tc , Double_t na_tc) { Double_t lt = TMath::Log(tpc/ttc); Double_t gse_rb = ( tpc*fgse(rb_pc , 0.0 , 0.0 ) - ttc*fgse(rb_tc , 0.0 , 0.0 ) )/9302.0/lt; Double_t gse_k = ( tpc*fgse(0.0 , k_pc , 0.0 ) - ttc*fgse(0.0 , k_tc , 0.0 ) )/10253.0/lt; Double_t gse_na = ( tpc*fgse(0.0 , 0.0 , na_pc) - ttc*fgse(0.0 , 0.0 , na_tc) )/12381.0/lt; return gse_rb + gse_k + gse_na; } // transfer tube relaxation rate in 1/hrs //PG-121 Double_t fgtt(Double_t gsett , Double_t X , Double_t gwall , Double_t gdip_pc , Double_t gdip_tc) { return gsett*(1.0+X) + gwall + TMath::Sqrt(gdip_pc*gdip_tc); } // slow polarization rate in 1/hrs, full calc //PG-12 Double_t fgslow(Double_t dtc , Double_t dpc , Double_t gse , Double_t gpc , Double_t gtc) { Double_t u = (gse+gpc-gtc)/(dpc+dtc); Double_t ftcmfpc = (dpc-dtc)/(dpc+dtc); Double_t gs = TMath::Sqrt(1.0+(2.0*ftcmfpc+u)*u); gs = 1.0 - gs; gs = gs*(dtc+dpc); gs = gs + gse + gpc + gtc; gs = gs/2.0; return gs; } // slow polarization rate in 1/hrs, 0th order approx //PG-21 with delta gamma = 0 Double_t fgslow0(Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc) { return fpc*(gse+gpc)+ftc*gtc; } // slow polarization rate in 1/hrs, 1st order approx //PG-21 Double_t fgslow1(Double_t dtc , Double_t dpc , Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc) { Double_t u = (gse+gpc-gtc)/(dpc+dtc); return fpc*(gse+gpc) + ftc*gtc - fpc*ftc*(dpc+dtc)*u*u; } // fast polarization rate in 1/hrs, full calc //PG-13 Double_t fgfast(Double_t dtc , Double_t dpc , Double_t gse , Double_t gpc , Double_t gtc) { Double_t u = (gse+gpc-gtc)/(dpc+dtc); Double_t ftcmfpc = (dpc-dtc)/(dpc+dtc); Double_t gf = TMath::Sqrt(1.0+(2.0*ftcmfpc+u)*u); gf = 1.0 + gf; gf = gf*(dtc+dpc); gf = gf + gse + gpc + gtc; gf = gf/2.0; return gf; } // fast polarization rate in 1/hrs, 0th order approx //PG-22 with delta gamma = 0 Double_t fgfast0(Double_t dtc , Double_t dpc , Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc) { return dtc + dpc + gse + gpc + gtc - fpc*(gse+gpc) - ftc*gtc; } // fast polarization rate in 1/hrs, 1st order approx //PG-22 Double_t fgfast1(Double_t dtc , Double_t dpc , Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc) { Double_t u = (gse+gpc-gtc)/(dpc+dtc); return dtc + dpc + gse + gpc + gtc - fpc*(gse+gpc) - ftc*gtc + fpc*ftc*(dpc+dtc)*u*u; } // equilibrium polarization in pumping chamber //PG-17 Double_t fppceq2(Double_t PA , Double_t dtc , Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc) { Double_t top = PA*gse*fpc; Double_t bot = gse*fpc+gpc*fpc+ftc*gtc/(1.0+gtc/dtc); return top/bot; } // equilibrium polarization in target chamber //PG-18 Double_t fptceq2(Double_t PA , Double_t dtc , Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc) { Double_t ppceq2 = fppceq2(PA , dtc , gse , gpc , gtc , fpc , ftc); Double_t bot = (1.0+gtc/dtc); return ppceq2/bot; } // polarization in pumping chamber at time t (in hours) , single value //PG-14 Double_t fppc2(Double_t t, Double_t PA , Double_t dtc , Double_t dpc , Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc, Double_t ppc0 , Double_t ptc0) { Double_t ppceq2 = fppceq2(PA , dtc , gse , gpc , gtc , fpc , ftc); Double_t ptceq2 = fptceq2(PA , dtc , gse , gpc , gtc , fpc , ftc); Double_t gslow = fgslow(dtc , dpc , gse , gpc , gtc); Double_t gfast = fgfast(dtc , dpc , gse , gpc , gtc); Double_t a = -(gse+gpc+dpc); Double_t b = dpc; Double_t c = dtc; Double_t d = -(gtc+dtc); Double_t B = gse*PA; Double_t cpc = gslow*(ppceq2-ppc0) - b*ptc0 - a*ppc0 - B; cpc = cpc/(gfast-gslow); Double_t ctc = gslow*(ptceq2-ptc0) - d*ptc0 - c*ppc0; ctc = ctc/(gfast-gslow); return ppceq2 + (ppc0-ppceq2-cpc)*TMath::Exp(-gslow*t) + cpc*TMath::Exp(-gfast*t); } // polarization in pumping chamber at time t (in hours) , functional form //PG-14 Double_t t_fppc2(Double_t *x , Double_t *par) { Double_t t = x[0]; Double_t PA = par[0]; Double_t dtc = par[1]; Double_t dpc = par[2]; Double_t gse = par[3]; Double_t gpc = par[4]; Double_t gtc = par[5]; Double_t fpc = par[6]; Double_t ftc = par[7]; Double_t ppc0 = par[8]; Double_t ptc0 = par[9]; return fppc2(t,PA,dtc,dpc,gse,gpc,gtc,fpc,ftc,ppc0,ptc0); } // polarization in target chamber at time t (in hours) , single value //PG-15 Double_t fptc2(Double_t t, Double_t PA , Double_t dtc , Double_t dpc , Double_t gse , Double_t gpc , Double_t gtc , Double_t fpc , Double_t ftc, Double_t ppc0 , Double_t ptc0) { Double_t ppceq2 = fppceq2(PA , dtc , gse , gpc , gtc , fpc , ftc); Double_t ptceq2 = fptceq2(PA , dtc , gse , gpc , gtc , fpc , ftc); Double_t gslow = fgslow(dtc , dpc , gse , gpc , gtc); Double_t gfast = fgfast(dtc , dpc , gse , gpc , gtc); Double_t a = -(gse+gpc+dpc); Double_t b = dpc; Double_t c = dtc; Double_t d = -(gtc+dtc); Double_t B = gse*PA; Double_t cpc = gslow*(ppceq2-ppc0) - b*ptc0 - a*ppc0 - B; cpc = cpc/(gfast-gslow); Double_t ctc = gslow*(ptceq2-ptc0) - d*ptc0 - c*ppc0; ctc = ctc/(gfast-gslow); return ptceq2 + (ptc0-ptceq2-ctc)*TMath::Exp(-gslow*t) + ctc*TMath::Exp(-gfast*t); } // polarization in target chamber at time t (in hours) , functional form //PG-15 Double_t t_fptc2(Double_t *x , Double_t *par) { Double_t t = x[0]; Double_t PA = par[0]; Double_t dtc = par[1]; Double_t dpc = par[2]; Double_t gse = par[3]; Double_t gpc = par[4]; Double_t gtc = par[5]; Double_t fpc = par[6]; Double_t ftc = par[7]; Double_t ppc0 = par[8]; Double_t ptc0 = par[9]; return fptc2(t,PA,dtc,dpc,gse,gpc,gtc,fpc,ftc,ppc0,ptc0); } // energy loss due to ionizing colliions //PG-44,45,46,47,48,49, TABLE 3 Double_t fdEdx(Double_t ebeam = 2000.0) { Double_t me = 0.510998918; Double_t g = ebeam/me; Double_t g2 = TMath::Power(g,2); Double_t fg = (1.0+2.0/g-2.0/g2)*TMath::Log(2) - 1.0/8.0*TMath::Power(1.0-1.0/g,2) - 1.0/g2; Double_t Z = 2.0; Double_t ibb = 41.8/1.0E6; Double_t cdelta = -11.1393; Double_t x0 = 2.2017; Double_t x1 = 3.6122; Double_t a = 0.13443; Double_t m = 5.8347; Double_t beta = TMath::Sqrt(1.0-1.0/g2); Double_t beta2 = TMath::Power(beta,2); Double_t x = TMath::Log10(beta*g); Double_t delta = 0.0; if (x < x0) { delta = 0.0; } else if (x < x1) { delta = 4.6052*x+cdelta+a*TMath::Power(x1-x,m); } else { delta = 4.6052*x+cdelta; } Double_t dedx = 2.0*TMath::Log(g-1); dedx += TMath::Log(g+1.0); dedx += 2.0*TMath::Log(me/ibb); dedx -= fg; dedx -= delta; dedx = dedx*Z/beta2; Double_t c0 = 510000.0*(1.0E-24)/2.0; return dedx*c0; } // maximum fraction of ions that are molecular ions //PG-71 Double_t fnmax(Double_t he , Double_t n2) { Double_t he2 = TMath::Power(he,2); Double_t gmnh = 29.0E+6; // Double_t mi2ctn2 = (30.0E+9)*n2; Double_t mi3ctn2 = (9.8E+9)*n2*he; Double_t ai2ctn2 = (27.0E+9)*n2; Double_t miform = (60.0E+6)*he2; return gmnh/(mi3ctn2+mi3ctn2)/(1.0+ai2ctn2/miform); } // beam relaxation rate in 1/hours //PG-42 Double_t fgbeam(Double_t na , Double_t nmax , Double_t I , Double_t atc , Double_t ebeam = 2000.0) { Double_t dedx = fdEdx(ebeam); Double_t e = 1.60217653E-19; Double_t Ei = 43.2; Double_t c0 = 3600.0; c0 = c0/(1.0E+6); return (dedx/e/Ei)*(I/atc)*(na+nmax/2.0)*c0; } // pumping to target chamber temperature ratio //PG-86 Double_t ft(Double_t tpc , Double_t ttc) { return tpc/ttc; } // pumping to target chamber volume ratio //PG-86 Double_t fv(Double_t vpc , Double_t vtc) { return vpc/vtc; } // pumping chamber operating density in amg //PG-86 Double_t fnpc(Double_t nfill , Double_t t , Double_t v) { return nfill*(1.0+v)/(t+v); } // target chamber operating density in amg //PG-86 Double_t fntc(Double_t nfill , Double_t t , Double_t v) { return t*nfill*(1.0+v)/(t+v); } // fraction of nuclei in target chamber //PG-88 Double_t fftc(Double_t t , Double_t v) { return t/(t+v); } // fraction of nuclei in pumping chamber //PG-87 Double_t ffpc(Double_t t , Double_t v) { return v/(t+v); } // approx fraction of nuclei in transfer tube //PG-100 Double_t fftt(Double_t t , Double_t fpc , Double_t ftc , Double_t vtt , Double_t vtc) { return vtt/vtc/(1.0+fpc/ftc)/(t-1.0)*TMath::Log(t); } // volume of pumping chamber assuming spherical shape Double_t fvpc(Double_t rpc) { return 4.0*TMath::Pi()/3.0*TMath::Power(rpc,3); } //volume of transfer tube cylinder Double_t fvtt(Double_t att , Double_t ltt) { return att*ltt; } //volume of target chamber cylinder Double_t fvtc(Double_t atc , Double_t ltc) { return atc*ltc; } // pressure in atm from a given density and temp of an ideal gas Double_t fpress(Double_t den , Double_t temp) { return den*temp/273.15; } // relative polarization gradient in percent, 0th order cal //PG-129 ignoring all terms beyond unity Double_t fpg_0(Double_t gtc , Double_t dtc) { return 100.0*gtc/dtc; } // relative polarization gradient in percent, 1st order cal //PG-129 Double_t fpgnb_1(Double_t gtc0 , Double_t dtc , Double_t gsett , Double_t gpc , Double_t ftt , Double_t ftc) { Double_t scale = (gsett/gtc0+TMath::Sqrt(gpc/gtc0))*ftt/ftc/2.0-gtc0/dtc; return 100.0*gtc0/dtc*(1.0+scale); } // relative polarization gradient due to beam only in percent, 1st order cal //PG-136 Double_t fpgb_1(Double_t gbeam , Double_t dtc , Double_t gtc0) { Double_t scale = (2.0*gtc0+gbeam)/dtc; return 100.0*gbeam/dtc*(1.0-scale); } // relative polarization gradient in percent, full calc for a two chambered cell //PG-125 Double_t fpg2_full(Double_t gtc , Double_t dtc) { return 100.0*(1.0-1.0/(1.0+gtc/dtc)); } // relative polarization gradient in percent, full calc for a three chambered cell //PG-124 Double_t fpg3_full(Double_t gtc , Double_t dtc , Double_t ftt , Double_t ftc , Double_t gtt) { Double_t u = gtc/dtc; Double_t ratio = 1.0/(1.0+u+(1+u/2.0)*ftt/ftc*gtt/2.0/dtc); return 100.0*(1.0 - ratio); } //Double_t expfit(Double_t *x , Double_t *par) //{ // Double_t t = x[0]; // Double_t A = par[0]; // Double_t tau = par[1]; // return A*(1.0-TMath::Exp(-t/tau)); //} //Double_t rnum(Double_t num , Int_t digits) //{ // Double_t scale = TMath::Power(10.0,digits); // Int_t temp = TMath::FloorNint(num*scale); // return 1.0*temp/scale; //} // calculates all parameters and returns the selected parameter // essentially the "heart" of this program Double_t mp_yaxis(Double_t *x, Double_t *par) { // inputs Double_t in_pcd = x[0]; Double_t rpc = in_pcd/2.0*2.54; Double_t tpc = par[0]+273.15; Double_t choice = par[1]; Double_t PA = par[2]; Double_t X = 3.0*par[3]/rpc; Double_t center_2_center = par[4]; Double_t att = par[5]; Double_t atc = par[6]; Double_t ltc = par[7]; Double_t ttc = par[8] + 273.15; Double_t ntc = par[9]; Double_t mf_k = par[10]; Double_t mf_rb = par[11]; Double_t mf_na = par[12]; Double_t I = par[13]; Double_t ebeam_MeV = par[14]; Double_t troom = par[15] + 273.15; Double_t gwall = 1.0/par[16]; Double_t n2_to_he = par[17]; // calculated outputs // to choose a particular output, specify the variable name as the first argument // for example, for the transfer tube length, the variable is ltt => so the argument would be "ltt" Double_t dedx = fdEdx(ebeam_MeV); Double_t rtc = TMath::Sqrt(atc/TMath::Pi()); Double_t ltt = center_2_center-rtc-rpc; Double_t vpc = fvpc(rpc); Double_t vtt = fvtt(att,ltt); Double_t vtc = fvtc(atc,ltc); Double_t t = ft(tpc,ttc); Double_t v = fv(vpc,vtc); Double_t he = ntc/t*(t+v)/(1.0+v); Double_t n2 = n2_to_he*he; Double_t fpc = ffpc(t,v); Double_t ftc = fftc(t,v); Double_t ftt = fftt(t,fpc,ftc,vtt,vtc); Double_t npc = fnpc(he,t,v); Double_t p_room = fpress(he,troom); Double_t p_op = fpress(npc,tpc); Double_t ups = fups(tpc,ttc); Double_t rb_pc = mf_rb*falknd("Rb",tpc); Double_t rb_tc = mf_rb*falknd("Rb",ttc); Double_t rb_pcm = falkmass("Rb",rb_pc,vpc); Double_t rb_tcm = falkmass("Rb",rb_tc,vtc); Double_t k_pc = mf_k*falknd("K",tpc); Double_t k_tc = mf_k*falknd("K",ttc); Double_t k_pcm = falkmass("K",k_pc,vpc); Double_t k_tcm = falkmass("K",k_tc,vtc); Double_t na_pc = mf_na*falknd("Na",tpc); Double_t na_tc = mf_na*falknd("Na",ttc); Double_t na_pcm = falkmass("Na",na_pc,vpc); Double_t na_tcm = falkmass("Na",na_tc,vtc); Double_t gsepc = fgse(rb_pc,k_pc,na_pc); Double_t gsetc = fgse(rb_tc,k_tc,na_tc); Double_t gsett = fgsett(tpc,rb_pc,k_pc,na_pc,ttc,rb_tc,k_tc,na_tc); Double_t gdip_pc = fgdip(npc,tpc); Double_t gdip_tc = fgdip(ntc,ttc); Double_t tau = 1.0/(gwall + fgdip(he,troom)); Double_t h = fh(ntc); Double_t p = fp(ntc,n2*ntc/he); Double_t na_fit = fna_fit(h,p); Double_t Omega = fomega(ntc); Double_t r = fr(ntc,n2*ntc/he); Double_t na_full = fna_full(r,Omega); Double_t nmax = fnmax(ntc,n2*ntc/he); Double_t gbeam = fgbeam(na_full,nmax,I,atc,ebeam_MeV); Double_t dtc = fdtc(att,ltt,vtc,ntc,t,ttc); Double_t dpc = fdpc(att,ltt,vtc,ntc,t,ttc,v); Double_t gpc = fgpc(gsepc,X,gwall,gdip_pc); Double_t gtc0 = fgtc(gsetc,X,gwall,gdip_tc,0.0); Double_t gtc = fgtc(gsetc,X,gwall,gdip_tc,gbeam); Double_t gtt = fgtt(gsett,X,gwall,gdip_pc,gdip_tc); Double_t gfast0 = fgfast0(dtc,dpc,gsepc,gpc,gtc,fpc,ftc); Double_t gfast1 = fgfast1(dtc,dpc,gsepc,gpc,gtc,fpc,ftc); Double_t gfast = fgfast(dtc,dpc,gsepc,gpc,gtc); Double_t gslow0 = fgslow0(gsepc,gpc,gtc,fpc,ftc); Double_t gslow1 = fgslow1(dtc,dpc,gsepc,gpc,gtc,fpc,ftc); Double_t gslow = fgslow(dtc,dpc,gsepc,gpc,gtc); Double_t p_pc_eq2 = fppceq2(PA,dtc,gsepc,gpc,gtc0,fpc,ftc); Double_t p_tc_eq2 = fptceq2(PA,dtc,gsepc,gpc,gtc0,fpc,ftc); Double_t pgnb_0 = fpg_0(gtc0,dtc); Double_t pgb_0 = fpg_0(gbeam,dtc); Double_t pgtot_0 = pgnb_0+pgb_0; Double_t pgnb_1 = fpgnb_1(gtc0,dtc,gsett,gpc,ftt,ftc); Double_t pgb_1 = fpgb_1(gbeam,dtc,gtc0); Double_t pgtot_1 = pgnb_1+pgb_1; Double_t pg2_full = fpg2_full(gtc,dtc); Double_t pg3_full = fpg3_full(gtc,dtc,ftt,ftc,gtt); Double_t fpc_ftc = fpc/ftc; Double_t p_pc_eq2b = fppceq2(PA,dtc,gsepc,gpc,gtc,fpc,ftc); Double_t p_tc_eq2b = fptceq2(PA,dtc,gsepc,gpc,gtc,fpc,ftc); Double_t pc_b2n = p_pc_eq2b/p_pc_eq2; Double_t tc_b2n = p_tc_eq2b/p_tc_eq2; Double_t gse_eff = fpc*gsepc; Double_t g_eff = fpc*gpc + ftc*gtc; Double_t output = 0; // insert "code2.txt" here if (choice == 0) { output = in_pcd; } else if (choice == 1) { output = rpc; } else if (choice == 2) { output = tpc; } else if (choice == 3) { output = choice; } else if (choice == 4) { output = PA; } else if (choice == 5) { output = X; } else if (choice == 6) { output = center_2_center; } else if (choice == 7) { output = att; } else if (choice == 8) { output = atc; } else if (choice == 9) { output = ltc; } else if (choice == 10) { output = ttc; } else if (choice == 11) { output = ntc; } else if (choice == 12) { output = mf_k; } else if (choice == 13) { output = mf_rb; } else if (choice == 14) { output = mf_na; } else if (choice == 15) { output = I; } else if (choice == 16) { output = ebeam_MeV; } else if (choice == 17) { output = troom; } else if (choice == 18) { output = gwall; } else if (choice == 19) { output = n2_to_he; } else if (choice == 20) { output = dedx; } else if (choice == 21) { output = rtc; } else if (choice == 22) { output = ltt; } else if (choice == 23) { output = vpc; } else if (choice == 24) { output = vtt; } else if (choice == 25) { output = vtc; } else if (choice == 26) { output = t; } else if (choice == 27) { output = v; } else if (choice == 28) { output = he; } else if (choice == 29) { output = n2; } else if (choice == 30) { output = fpc; } else if (choice == 31) { output = ftc; } else if (choice == 32) { output = ftt; } else if (choice == 33) { output = npc; } else if (choice == 34) { output = p_room; } else if (choice == 35) { output = p_op; } else if (choice == 36) { output = ups; } else if (choice == 37) { output = rb_pc; } else if (choice == 38) { output = rb_tc; } else if (choice == 39) { output = rb_pcm; } else if (choice == 40) { output = rb_tcm; } else if (choice == 41) { output = k_pc; } else if (choice == 42) { output = k_tc; } else if (choice == 43) { output = k_pcm; } else if (choice == 44) { output = k_tcm; } else if (choice == 45) { output = na_pc; } else if (choice == 46) { output = na_tc; } else if (choice == 47) { output = na_pcm; } else if (choice == 48) { output = na_tcm; } else if (choice == 49) { output = gsepc; } else if (choice == 50) { output = gsetc; } else if (choice == 51) { output = gsett; } else if (choice == 52) { output = gdip_pc; } else if (choice == 53) { output = gdip_tc; } else if (choice == 54) { output = tau; } else if (choice == 55) { output = h; } else if (choice == 56) { output = p; } else if (choice == 57) { output = na_fit; } else if (choice == 58) { output = Omega; } else if (choice == 59) { output = r; } else if (choice == 60) { output = na_full; } else if (choice == 61) { output = nmax; } else if (choice == 62) { output = gbeam; } else if (choice == 63) { output = dtc; } else if (choice == 64) { output = dpc; } else if (choice == 65) { output = gpc; } else if (choice == 66) { output = gtc0; } else if (choice == 67) { output = gtc; } else if (choice == 68) { output = gtt; } else if (choice == 69) { output = gfast0; } else if (choice == 70) { output = gfast1; } else if (choice == 71) { output = gfast; } else if (choice == 72) { output = gslow0; } else if (choice == 73) { output = gslow1; } else if (choice == 74) { output = gslow; } else if (choice == 75) { output = p_pc_eq2; } else if (choice == 76) { output = p_tc_eq2; } else if (choice == 77) { output = pgnb_0; } else if (choice == 78) { output = pgb_0; } else if (choice == 79) { output = pgtot_0; } else if (choice == 80) { output = pgnb_1; } else if (choice == 81) { output = pgb_1; } else if (choice == 82) { output = pgtot_1; } else if (choice == 83) { output = pg2_full; } else if (choice == 84) { output = pg3_full; } else if (choice == 85) { output = fpc_ftc; } else if (choice == 86) { output = p_pc_eq2b; } else if (choice == 87) { output = p_tc_eq2b; } else if (choice == 88) { output = pc_b2n; } else if (choice == 89) { output = tc_b2n; } else if (choice == 90) { output = gse_eff; } else if (choice == 91) { output = g_eff; } else if (choice == 92) { output = output; } // insert "code2.txt" here return output; } void make_pcd_plots(TString in_choice_str = "ltt", TString in_yaxis = "transfer tube length (cm)" , Int_t num_temps = 1, TString filename = "technote_parameters.txt", Double_t in_xmin = 2.5, Double_t in_xmax = 3.5 , Double_t in_tpc_mid = 260.0 , Double_t in_deltaT = 30.0) { gROOT->Reset(); gROOT->DeleteAll(); cout << "START" << endl; TStopwatch timer; timer.Start(); // reading in parameters ifstream datafile; datafile.open(filename.Data()); cout << "Reading Data from: " << filename << endl; Int_t count = 2; const Int_t num_par = 18; Double_t in_par[num_par]; Double_t in_temp; while ( !(datafile.eof()) && (count <= num_par-1) ) { if (datafile.good()) { datafile >> in_temp; if (!(datafile.eof())) { in_par[count] = in_temp; count++; } } else { cout << "Trouble reading file" << endl << "Ending file: " << filename << endl; exit(0); } } datafile.close(); cout << "Read " << count-2 << " parameters and closed file: " << filename << endl; if (count-2 != num_par-2) { cout << "Number of parameters read (" << count-2 << ") does not equal required number of parameters (" << num_par-2 << ")" << endl; cout << "Ending file: " << filename << endl; exit(0); } // outputting values cout << "PA = " << in_par[2] << endl; // unitless (0-1) cout << "X_max = " << in_par[3] << endl; // unitless (0-1) cout << "center_2_center (cm) = " << in_par[4] << endl; // cm (>0) cout << "att (cm^2) = " << in_par[5] << endl; // cm^2 (>0) cout << "atc (cm^2) = " << in_par[6] << endl; // cm^2 (>0) cout << "ltc (cm) = " << in_par[7] << endl; // cm (>0) cout << "ttc (K) = " << in_par[8] + 273.15 << endl; // C cout << "ntc (amg) = " << in_par[9] << endl; // amg (>0) cout << "mf_k = " << in_par[10] << endl; // unitless (0-1) cout << "mf_rb = " << in_par[11] << endl; // unitless (0-1) cout << "mf_na = " << in_par[12] << endl; // unitless (0-1) cout << "I (microamps) = " << in_par[13] << endl; // microamps (>0) cout << "ebeam (MeV) = " << in_par[14] << endl; // MeV (>0) cout << "troom (K) = " << in_par[15] + 273.15 << endl; // C (>0) cout << "gwall (1/hrs) = " << 1.0/in_par[16] << endl; // hours (>0) cout << "n2_to_he = " << in_par[17] << endl; // unitless (>0) // common to all plots gROOT->SetStyle("Plain"); gROOT->ForceStyle(); TStyle *basestyle = gROOT->GetStyle("Plain"); basestyle->SetOptTitle(0); basestyle->SetOptStat(0); basestyle->SetOptTitle(0); basestyle->SetOptStat(0); basestyle->SetLabelSize(0.0,"X"); basestyle->SetLabelSize(0.0,"Y"); basestyle->SetTickLength(0.0,"X"); basestyle->SetTickLength(0.0,"Y"); basestyle->SetLabelOffset(0.01,"X"); basestyle->SetLabelOffset(0.01,"Y"); // defining plot sizes Double_t a = 1.0; Double_t h0 = 800; Double_t w0 = a*h0; Double_t pad1l = h0/5; Double_t pad1b = h0/8; Double_t pad1t = h0/20; Double_t pad1r = h0/20; Double_t pad1w = pad1l+w0+pad1r; Double_t pad1h = pad1t+h0+pad1b; pad1r = pad1r/pad1w; pad1l = pad1l/pad1w; pad1t = pad1t/pad1h; pad1b = pad1b/pad1h; Int_t canh = TMath::Nint(pad1h)+1; Int_t canw = TMath::Nint(pad1w)+1; //creating canvas and main pad TCanvas *canvas = new TCanvas("canvas","Canvas Title",canw,canh); canvas->cd(); basestyle->SetPadLeftMargin(pad1l); basestyle->SetPadRightMargin(pad1r); basestyle->SetPadTopMargin(pad1t); basestyle->SetPadBottomMargin(pad1b); TPad *pad1 = new TPad("pad1","AllP",0.000,0.0,pad1w/canw,1.0); pad1->SetLineWidth(2); pad1->Draw(); pad1->cd(); // setting up function Int_t in_choice = 0; // insert "code1.txt" here if (in_choice_str == "in_pcd") { in_choice = 0; } else if (in_choice_str == "rpc") { in_choice = 1; } else if (in_choice_str == "tpc") { in_choice = 2; } else if (in_choice_str == "choice") { in_choice = 3; } else if (in_choice_str == "PA") { in_choice = 4; } else if (in_choice_str == "X") { in_choice = 5; } else if (in_choice_str == "center_2_center") { in_choice = 6; } else if (in_choice_str == "att") { in_choice = 7; } else if (in_choice_str == "atc") { in_choice = 8; } else if (in_choice_str == "ltc") { in_choice = 9; } else if (in_choice_str == "ttc") { in_choice = 10; } else if (in_choice_str == "ntc") { in_choice = 11; } else if (in_choice_str == "mf_k") { in_choice = 12; } else if (in_choice_str == "mf_rb") { in_choice = 13; } else if (in_choice_str == "mf_na") { in_choice = 14; } else if (in_choice_str == "I") { in_choice = 15; } else if (in_choice_str == "ebeam_MeV") { in_choice = 16; } else if (in_choice_str == "troom") { in_choice = 17; } else if (in_choice_str == "gwall") { in_choice = 18; } else if (in_choice_str == "n2_to_he") { in_choice = 19; } else if (in_choice_str == "dedx") { in_choice = 20; } else if (in_choice_str == "rtc") { in_choice = 21; } else if (in_choice_str == "ltt") { in_choice = 22; } else if (in_choice_str == "vpc") { in_choice = 23; } else if (in_choice_str == "vtt") { in_choice = 24; } else if (in_choice_str == "vtc") { in_choice = 25; } else if (in_choice_str == "t") { in_choice = 26; } else if (in_choice_str == "v") { in_choice = 27; } else if (in_choice_str == "he") { in_choice = 28; } else if (in_choice_str == "n2") { in_choice = 29; } else if (in_choice_str == "fpc") { in_choice = 30; } else if (in_choice_str == "ftc") { in_choice = 31; } else if (in_choice_str == "ftt") { in_choice = 32; } else if (in_choice_str == "npc") { in_choice = 33; } else if (in_choice_str == "p_room") { in_choice = 34; } else if (in_choice_str == "p_op") { in_choice = 35; } else if (in_choice_str == "ups") { in_choice = 36; } else if (in_choice_str == "rb_pc") { in_choice = 37; } else if (in_choice_str == "rb_tc") { in_choice = 38; } else if (in_choice_str == "rb_pcm") { in_choice = 39; } else if (in_choice_str == "rb_tcm") { in_choice = 40; } else if (in_choice_str == "k_pc") { in_choice = 41; } else if (in_choice_str == "k_tc") { in_choice = 42; } else if (in_choice_str == "k_pcm") { in_choice = 43; } else if (in_choice_str == "k_tcm") { in_choice = 44; } else if (in_choice_str == "na_pc") { in_choice = 45; } else if (in_choice_str == "na_tc") { in_choice = 46; } else if (in_choice_str == "na_pcm") { in_choice = 47; } else if (in_choice_str == "na_tcm") { in_choice = 48; } else if (in_choice_str == "gsepc") { in_choice = 49; } else if (in_choice_str == "gsetc") { in_choice = 50; } else if (in_choice_str == "gsett") { in_choice = 51; } else if (in_choice_str == "gdip_pc") { in_choice = 52; } else if (in_choice_str == "gdip_tc") { in_choice = 53; } else if (in_choice_str == "tau") { in_choice = 54; } else if (in_choice_str == "h") { in_choice = 55; } else if (in_choice_str == "p") { in_choice = 56; } else if (in_choice_str == "na_fit") { in_choice = 57; } else if (in_choice_str == "Omega") { in_choice = 58; } else if (in_choice_str == "r") { in_choice = 59; } else if (in_choice_str == "na_full") { in_choice = 60; } else if (in_choice_str == "nmax") { in_choice = 61; } else if (in_choice_str == "gbeam") { in_choice = 62; } else if (in_choice_str == "dtc") { in_choice = 63; } else if (in_choice_str == "dpc") { in_choice = 64; } else if (in_choice_str == "gpc") { in_choice = 65; } else if (in_choice_str == "gtc0") { in_choice = 66; } else if (in_choice_str == "gtc") { in_choice = 67; } else if (in_choice_str == "gtt") { in_choice = 68; } else if (in_choice_str == "gfast0") { in_choice = 69; } else if (in_choice_str == "gfast1") { in_choice = 70; } else if (in_choice_str == "gfast") { in_choice = 71; } else if (in_choice_str == "gslow0") { in_choice = 72; } else if (in_choice_str == "gslow1") { in_choice = 73; } else if (in_choice_str == "gslow") { in_choice = 74; } else if (in_choice_str == "p_pc_eq2") { in_choice = 75; } else if (in_choice_str == "p_tc_eq2") { in_choice = 76; } else if (in_choice_str == "pgnb_0") { in_choice = 77; } else if (in_choice_str == "pgb_0") { in_choice = 78; } else if (in_choice_str == "pgtot_0") { in_choice = 79; } else if (in_choice_str == "pgnb_1") { in_choice = 80; } else if (in_choice_str == "pgb_1") { in_choice = 81; } else if (in_choice_str == "pgtot_1") { in_choice = 82; } else if (in_choice_str == "pg2_full") { in_choice = 83; } else if (in_choice_str == "pg3_full") { in_choice = 84; } else if (in_choice_str == "fpc_ftc") { in_choice = 85; } else if (in_choice_str == "p_pc_eq2b") { in_choice = 86; } else if (in_choice_str == "p_tc_eq2b") { in_choice = 87; } else if (in_choice_str == "pc_b2n") { in_choice = 88; } else if (in_choice_str == "tc_b2n") { in_choice = 89; } else if (in_choice_str == "gse_eff") { in_choice = 90; } else if (in_choice_str == "g_eff") { in_choice = 91; } else if (in_choice_str == "output") { in_choice = 92; } // end "code1.txt" here Double_t xmin = in_xmin; Double_t xmax = in_xmax; in_par[0] = in_tpc_mid; in_par[1] = in_choice; TF1 *f0 = new TF1("f0",mp_yaxis,xmin,xmax,num_par); f0->SetParameters(in_par); Double_t ymin = 0.0; Double_t ymax = 0.0; Double_t y_xmin,y_xmax,ynmin,ynmax; // finding min and max y values Double_t tstart = in_tpc_mid-in_deltaT; if (num_temps == 1) { tstart = in_tpc_mid; } for (Int_t count = 0; count < num_temps; count++) { f0->SetParameter(0,tstart+count*in_deltaT); y_xmin = f0->Eval(xmin); y_xmax = f0->Eval(xmax); ynmin = TMath::Min(y_xmin,y_xmax); ynmax = TMath::Max(y_xmin,y_xmax); if (count == 0) { ymin = ynmin; ymax = ynmax; } else { if (ynmin < ymin) { ymin = ynmin; } if (ynmax > ymax) { ymax = ynmax; } } } cout << "Found ymin = " << ymin << " and ymax = " << ymax << endl; //creating main pad graph window and axes Double_t dx = (xmax-xmin)/10.0; Double_t dy = (ymax-ymin)/10.0; ymin -= dy; ymax += dy; xmin -= dx; xmax += dx; TH2F *frame = new TH2F("frame","",10,xmin,xmax,10,ymin,ymax); frame->SetLineWidth(2); frame->SetAxisRange(xmin,xmax,"X"); frame->SetAxisRange(ymin,ymax,"Y"); frame->GetXaxis()->SetTicks("-"); frame->GetYaxis()->SetTicks("+"); frame->Draw(); //the ususal axes TGaxis *x_axis = new TGaxis(xmin,ymin,xmax,ymin,xmin,xmax); TGaxis *y_axis = new TGaxis(xmin,ymin,xmin,ymax,ymin,ymax); x_axis->SetOption("-"); y_axis->SetOption("+"); x_axis->SetLabelOffset(-0.075); y_axis->SetLabelOffset(-0.035); x_axis->SetLabelSize(0.0275); y_axis->SetLabelSize(0.0275); x_axis->SetLineWidth(2); y_axis->SetLineWidth(2); x_axis->Draw(); y_axis->Draw(); //the other and unlabeled axes TGaxis *x2_axis = new TGaxis(xmin,ymax,xmax,ymax,xmin,xmax); x2_axis->SetOption("+U"); x2_axis->SetLabelOffset(-0.075); x2_axis->SetLabelSize(0.025); x2_axis->SetLineWidth(2); x2_axis->Draw(); TGaxis *y2_axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax); y2_axis->SetOption("-U"); y2_axis->SetLabelOffset(-0.035); y2_axis->SetLabelSize(0.025); y2_axis->SetLineWidth(2); y2_axis->Draw(); // drawing axis titles TLatex x_title; x_title.SetTextAlign(22); x_title.SetTextAngle(0); x_title.SetTextSize(0.025); x_title.DrawLatex((xmax+xmin)/2.0,ymin-(1.125)*dy,"pumping chamber diameter (in)"); TLatex y_title; y_title.SetTextAlign(22); y_title.SetTextAngle(90); y_title.SetTextSize(0.040); y_title.DrawLatex(xmin-(1.8)*dx,(ymax+ymin)/2.0,in_yaxis); // defining and plotting function cout << "done drawing plot frame!" << endl; f0->SetLineColor(1); f0->SetLineWidth(2); f0->SetParameter(0,in_tpc_mid); f0->SetLineStyle(1); f0->DrawCopy("LSAME"); if (num_temps > 1) { f0->SetParameter(0,in_tpc_mid-in_deltaT); f0->SetLineColor(4); f0->SetLineStyle(2); f0->DrawCopy("LSAME"); f0->SetParameter(0,in_tpc_mid+in_deltaT); f0->SetLineStyle(3); f0->SetLineColor(2); f0->DrawCopy("LSAME"); cout << "did all temps" << endl; } else { cout << "did only 260 C" << endl; } TString outeps = in_choice_str; outeps += ".eps"; canvas->Print(outeps.Data()); timer.Stop(); timer.Print(); cout << "END" << endl; }