/* Code for generation of initial baffle configuration This uses data from GEMC Seamus Riordan sriordan@physics.umass.edu April 12, 2011 */ //#define NO_EXT_TARG #define STRAIGHT_ENDS double beamE = 11.0;// GeV double thmin = 22.0;// deg double thmax = 35.0;// deg double targ_len = 0.4;// double targ_z = 0.1;// double xmin = 0.0; double xmax = 0.8; // x = 0.8 is about the W2 < 4 cutoff double pmin = 0.8; // x = 0.8 is about the W2 < 4 cutoff double SCALE_MIN_ANG = 1.4; // Artificial scale factor to increase double SCALE_MAX_ANG = 1.4; // Artificial scale factor to increase // efficiency // double SCALE_LASTBAF = 0.20; // scale last width to block more gammas double SCALE_LASTBAF = 0.; // scale last width to block more gammas double abspmin = 0.0; // Just ignore momenta less than this #define nbaf 6 double bafspace = 0.3;// m double minsep_scale = 1.0; // Scale minsep to // help prevent more photons double width_scale = 1.00; // Scale width to // help prevent more photons double baf_thick = 0.09; double maxr = 1.42; #define MAXSLIT 30 // Radial block granularity #define NBLOCK 20 #define MASS_P 0.938 #define PI 3.1415927 #define MAXHIT 2000 #include "TGraph.h" #include "TH1F.h" #include "TLine.h" #include "TH2F.h" #include "TF1.h" #include "TChain.h" #include "TCanvas.h" #include "TROOT.h" #include "TStyle.h" double getEf(double, double); void GetPhiWidth(TH1F *, double &, double &); void GetEugene( int i, double R, double &in, double &out ); int makebaf5(){ gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); double rmin[nbaf], rmax[nbaf]; double zbaf[nbaf]; int i, j; // Calculate minimum and maximum r for each range double lastz = bafspace*nbaf + targ_z; double firstz = bafspace + targ_z; double dlmin, dlmax; double thrmin = thmin*PI/180.0; double thrmax = thmax*PI/180.0; for( i = 0; i < nbaf; i++ ){ zbaf[i] = targ_z + bafspace*(i+1); printf("%d zbaf = %f\n", i, zbaf[i] ); /* dlmin = zbaf[i] - targ_len/2.0 - baf_thick/2.0; dlmax = zbaf[i] + targ_len/2.0 + baf_thick/2.0; */ dlmin = zbaf[i] - targ_len/2.0 - baf_thick - targ_z; dlmax = zbaf[i] + targ_len/2.0 + baf_thick - targ_z; rmin[i] = dlmin*tan(thrmin); rmax[i] = dlmax*tan(thrmax); printf("rmax (%d) = %f\n", i, rmax[i]); if( rmax[i] > maxr ){ printf("reducing %f to %f\n", rmax[i], maxr); rmax[i] = maxr; } } double absminr[6] = {3.9,15.3,26.6,37.9,49.2,60.4}; // cm // Make sure this has a minimum angle enforced for( i = 0; i < nbaf; i++ ){ if( rmin[i] < absminr[i]*0.01 ){ rmin[i] = absminr[i]*0.01; } } // Break down phi profiles TChain *T = new TChain("out"); // T->Add("testcleo/solout_CLEO_LD2_120d_nobaf_11GeV.root"); // T->Add("solout_CLEO_LD2_120d_nobaf_11GeV.root"); // T->Add("solout_BaBar_nobaf_LD2_120d_11GeV.root"); //T->Add("oldcleo/solout_CLEO_nobaf_LD2_120d_11GeV.root"); T->Add("cleo/solout_CLEO_LD2_120d_nobaf_11GeV.root"); int nhit; double x[MAXHIT],y[MAXHIT],z[MAXHIT],px[MAXHIT],py[MAXHIT],pz[MAXHIT],p0x[MAXHIT],p0y[MAXHIT],p0z[MAXHIT]; double v0x[MAXHIT],v0y[MAXHIT],v0z[MAXHIT],weight[MAXHIT], trid[MAXHIT]; double th, phi, Ef, nu, xbj, Q2, phi2, W; bool isfill; T->SetBranchAddress("nhits", &nhit); T->SetBranchAddress("x", x); T->SetBranchAddress("y", y); T->SetBranchAddress("z", z); T->SetBranchAddress("px", px); T->SetBranchAddress("py", py); T->SetBranchAddress("pz", pz); T->SetBranchAddress("p0x", p0x); T->SetBranchAddress("p0y", p0y); T->SetBranchAddress("p0z", p0z); T->SetBranchAddress("v0x", v0x); T->SetBranchAddress("v0y", v0y); T->SetBranchAddress("v0z", v0z); T->SetBranchAddress("weight", weight); T->SetBranchAddress("trid", trid); TH1F *hblock[nbaf][NBLOCK]; // Histograms will store information on phi distributions for( i = 0; i < nbaf; i++ ){ for( j = 0; j < NBLOCK; j++ ){ hblock[i][j] = new TH1F(Form("hb_%d_%d", i,j), Form("hb_%d_%d", i, j), 200, 0.0, 45.0); } } for( i = 0; i < T->GetEntries(); i++ ){ if( (i % 100000)==0 ){ printf("Event %8d\n", i);} T->GetEntry(i); for( j = 0; j < nhit; j++ ){ // Convert to GeV p0x[j] *= 1e-3; p0y[j] *= 1e-3; p0z[j] *= 1e-3; px[j] *= 1e-3; py[j] *= 1e-3; pz[j] *= 1e-3; } th = fabs(atan(sqrt(p0x[0]*p0x[0]+p0y[0]*p0y[0])/p0z[0]))*180.0/3.14159; phi = atan2(p0y[0], p0x[0])*180.0/3.14159; if( weight[0] > 2e11 ) continue; Ef = sqrt(p0x[0]*p0x[0]+p0y[0]*p0y[0]+p0z[0]*p0z[0]); nu = 11.0 - Ef; Q2 = 2.0*11.0*Ef*(1.0-cos(th*3.14159/180.0)); xbj = Q2/(2.0*nu*0.938); W = sqrt(0.938*0.938 + 2.0*0.938*nu - Q2); bool hasfront = false; // Base cuts if( W*W < 4.0 ) continue; // if( xbj < xmin || xbj > xmax || Ef < pmin ) continue; /* for( j = 0; j < nhit && !hasfront; j++ ){ if( z[j] < zfirst[magidx] && fabs(trid[j]-1)<1e-3 ){ hasfront = true; } } isfill = false; */ // if( hasfront ){ for( j = 0; j < nhit; j++ ){ double r = sqrt(x[j]*x[j]+y[j]*y[j])/1000.0; // Get baffle this is associated with int thisbaf = (int) ((z[j]/1000.0-targ_z)/bafspace); if( thisbaf >= nbaf ) continue; if( thisbaf < 0 ) continue; //printf("(%f, %f)-> %d, %f < %f < %f (%f deg)\n", r, z[j]/1000.0, thisbaf, rmin[thisbaf], r, rmax[thisbaf], th); // Get block bin if( r < rmin[thisbaf] || r > rmax[thisbaf] ) continue; int thisblock = (r-rmin[thisbaf])*NBLOCK/(rmax[thisbaf]-rmin[thisbaf]); if( fabs(trid[j]-1)<1e-3 ){ phi2 = atan2(y[j],x[j])*180/3.14159; if( xbj > 0.2 ){ isfill = true; // hphiR->Fill(r, phi2-phi, weight[j]); // printf("Filling %d block %d with %f\n", thisbaf, thisblock, (phi2-phi)*PI/180.0); hblock[thisbaf][thisblock]->Fill(phi2-phi, weight[j]); // hblock[thisbaf][thisblock]->Fill(phi2-phi); } } } } // } TCanvas *can; double phimin[nbaf][NBLOCK], phimax[nbaf][NBLOCK]; for( i = 0; i < nbaf; i++ ){ for( j = 0; j < NBLOCK; j++ ){ GetPhiWidth(hblock[i][j], phimin[i][j], phimax[i][j] ); } } // This block plots the widths /* for( i = 0; i < nbaf; i++ ){ can = new TCanvas(); can->Divide(5,4); for( j = 0; j < NBLOCK; j++ ){ can->cd(j+1); // gPad->SetLogy(); hblock[i][j]->SetTitle(Form("Baffle %d, Block %d", i,j)); hblock[i][j]->GetXaxis()->SetTitle("#phi [deg]"); hblock[i][j]->GetXaxis()->CenterTitle(); hblock[i][j]->Draw(); TLine *minl = new TLine( phimin[i][j], 0, phimin[i][j], hblock[i][j]->GetMaximum()*0.9 ); TLine *maxl = new TLine( phimax[i][j], 0, phimax[i][j], hblock[i][j]->GetMaximum()*0.9 ); minl->SetLineColor(kRed); maxl->SetLineColor(kGreen); minl->Draw(); maxl->Draw(); } } */ // Plot trajectories TGraph *ginner[NBLOCK]; TGraph *gouter[NBLOCK]; TF1 *finner[NBLOCK], *fouter[NBLOCK], *fsinner[NBLOCK], *fsouter[NBLOCK]; int ntpt[NBLOCK]; double phiinner[NBLOCK][nbaf], phiouter[NBLOCK][nbaf], ztpt[NBLOCK][nbaf]; double zpts[nbaf]; double iinter[NBLOCK], blockr[NBLOCK]; double ointer[NBLOCK]; for( i = 0; i < NBLOCK; i++ ){ ntpt[i] = 0; for( j = 0; j < nbaf; j++ ){ if( phimin[j][i] > 0 ){ phiinner[i][ntpt[i]] = phimin[j][i]; phiouter[i][ntpt[i]] = phimax[j][i]; ztpt[i][ntpt[i]] = zbaf[j]; zpts[j] = zbaf[j]; ntpt[i]++; } } ginner[i] = new TGraph(ntpt[i], ztpt[i], phiinner[i]); gouter[i] = new TGraph(ntpt[i], ztpt[i], phiouter[i]); finner[i] = new TF1(Form("finner_%d", i), "pol1", 0, 2.5); fouter[i] = new TF1(Form("fouter_%d", i), "pol1", 0, 2.5); fsinner[i] = new TF1(Form("fsinner_%d", i), "pol1", 0, 2.5); fsouter[i] = new TF1(Form("fsouter_%d", i), "pol1", 0, 2.5); fouter[i]->SetLineColor(kRed); fsouter[i]->SetLineColor(kRed); finner[i]->SetLineWidth(1); fouter[i]->SetLineWidth(1); fsinner[i]->SetLineWidth(1); fsouter[i]->SetLineWidth(1); ////////////////////////////////// ginner[i]->Fit(Form("finner_%d", i), "N"); gouter[i]->Fit(Form("fouter_%d", i), "N"); fsinner[i]->SetParameter(0, finner[i]->GetParameter(0) - finner[i]->Eval(firstz)); fsinner[i]->SetParameter(1, finner[i]->GetParameter(1)); fsouter[i]->SetParameter(0, fouter[i]->GetParameter(0) - fouter[i]->Eval(firstz)); fsouter[i]->SetParameter(1, fouter[i]->GetParameter(1)); //////////////////////////////////// gouter[i]->SetMarkerColor(kRed); gouter[i]->SetMarkerStyle(7); ginner[i]->SetMarkerStyle(7); iinter[i] = finner[i]->GetParameter(0); ointer[i] = fouter[i]->GetParameter(0); blockr[i] = i; } TGraph *giinter = new TGraph(NBLOCK, blockr, iinter); TGraph *gointer = new TGraph(NBLOCK, blockr, ointer); giinter->SetMarkerStyle(7); gointer->SetMarkerStyle(7); TF1 *fiinter = new TF1("fiinter", "pol1", 0, 18 ); TF1 *fointer = new TF1("fointer", "pol3", 0, 18 ); /* // Fit intercepts new TCanvas(); giinter->Draw("AP"); giinter->Fit("fiinter", "R"); new TCanvas(); gointer->Draw("AP"); gointer->Fit("fointer", "R"); */ for( i = 0; i < NBLOCK; i++ ){ /* finner[i]->FixParameter(0, giinter->Eval(i)); fouter[i]->FixParameter(0, gointer->Eval(i)); ginner[i]->Fit(Form("finner_%d", i), "N"); gouter[i]->Fit(Form("fouter_%d", i), "N"); */ } TH2F *frm[NBLOCK]; for( i = 0; i < NBLOCK; i++ ){ frm[i] = new TH2F(Form("frm_%d", i), Form("Block %d", i ), 100, 0, 2.5, 100, -1, 25); frm[i]->GetXaxis()->SetTitle("z [m]"); frm[i]->GetXaxis()->CenterTitle(); frm[i]->GetYaxis()->SetTitle("#phi [deg]"); frm[i]->GetYaxis()->CenterTitle(); } can = new TCanvas(); can->Divide(5,4); for( i = 0; i < NBLOCK; i++ ){ can->cd(i+1); frm[i]->Draw(); /* gouter[i]->Draw("P"); ginner[i]->Draw("P"); finner[i]->Draw("lsame"); fouter[i]->Draw("lsame"); */ } /////////////////////////////////////////////////////////////////////////// int rcnt, bcnt; double theta, bover, dradius, dtheta; double theta_min, theta_max; double phi_in, phi_out, Dphi; bool clip; double width, outerang, innerang; double r[nbaf][NBLOCK+1]; double sphi[nbaf][NBLOCK+1]; double dphi[nbaf][NBLOCK+1]; double thisz, thisthmin, thisthmax; double dphimin, dphimax; int nslit = 30; // Find maximum needed angle for all R double maxwidth = 0; double thiswidth; for( i = 0; i < NBLOCK; i++ ){ dphimin = finner[i]->Eval(lastz)-finner[i]->Eval(firstz); dphimax = fouter[i]->Eval(lastz)-fouter[i]->Eval(firstz); width = dphimin*width_scale; outerang = dphimax-width; thiswidth = outerang+width-dphimin; if( thiswidth > maxwidth ){ maxwidth = thiswidth; } } // int nslit = 360.0/maxwidth; printf("maxwidth = %f, nslit = %d\n", maxwidth, nslit); double minsep; double maxsect = 360.0/nslit; for( bcnt = 0; bcnt < nbaf; bcnt++ ){ thisz = bafspace*(bcnt+1) + targ_z; for( rcnt = 0; rcnt < NBLOCK; rcnt++ ){ dphimin = finner[rcnt]->Eval(lastz)-finner[rcnt]->Eval(firstz); dphimax = fouter[rcnt]->Eval(lastz)-fouter[rcnt]->Eval(firstz); printf("dphimin = %f, dphimax = %f\n", dphimin, dphimax); minsep = minsep_scale*dphimin/(nbaf-1); // printf("minsep = %f\n", minsep*180.0/3.14159); // Any bigger than this and photons // get through width = dphimin*width_scale; outerang = dphimax-width; innerang = dphimin; // Make sure we have enough room // up front if( (width+minsep) > maxsect ){ printf("Baffle %d block %d: reducing width... from %f to %f\n", bcnt, rcnt, width, (maxsect-minsep)); width = maxsect-minsep; outerang = dphimax-width; } // printf("width = %f\n", width); // Room we have at the back bover = maxsect - (outerang-innerang+width); // printf("maxsect %f\n", maxsect); if( bover < minsep ){ printf("pulling in outer angle from %f to %f\n", outerang, maxsect - minsep - width + innerang); outerang = maxsect - minsep - width + innerang; } // Try this to fill out space at the cost of very high x /* if( bover > minsep ){ printf("increasing inner angle from %f to %f\n", dphimin, innerang + bover); innerang += bover; width += bover; } */ // Scale along z double zscale = (thisz-firstz)/(lastz-firstz); printf("zscale = %f\n", zscale ); phi_in = innerang*zscale; phi_out = outerang*zscale + width; // printf("%f %f width %f\n", dphimin, outerang, width); double photonblock_in = 0.0; if( bcnt == nbaf-1 ){ // Kill last remaining photons photonblock_in = phi_in*SCALE_LASTBAF; } double backmiddle = 0.5*( (outerang - innerang) + width ) + innerang; r[bcnt][rcnt] = (rmax[bcnt]-rmin[bcnt])*((double) rcnt)/((double) NBLOCK) + rmin[bcnt]; // printf("%d %d rmax %f rmin %f scale %f -> %f\n", bcnt, rcnt, rmax[bcnt], rmin[bcnt], ((double) rcnt)/((double) NBLOCK ), r[bcnt][rcnt] ); #ifndef STRAIGHT_ENDS sphi[bcnt][rcnt] = phi_in + photonblock_in; dphi[bcnt][rcnt] = phi_out-phi_in-photonblock_in; #else sphi[bcnt][rcnt] = phi_in + photonblock_in - backmiddle + 16.; dphi[bcnt][rcnt] = phi_out-phi_in-photonblock_in; #endif//STRAIGHT_ENDS // printf("{%d, %f, %f},\n", rcnt, phi_in, phi_out-phi_in ); } printf("\n"); r[bcnt][NBLOCK] = rmax[bcnt]; } TGraph *gbinner[NBLOCK], *gbouter[NBLOCK], *gsinner[NBLOCK], *gsouter[NBLOCK]; double ibphi[NBLOCK][nbaf], obphi[NBLOCK][nbaf], shiftphiinner[NBLOCK][nbaf], shiftphiouter[NBLOCK][nbaf]; double eout[NBLOCK][nbaf], ein[NBLOCK][nbaf]; for( rcnt = 0; rcnt < NBLOCK; rcnt++ ){ for( bcnt = 0; bcnt < nbaf; bcnt++ ){ ibphi[rcnt][bcnt] = sphi[bcnt][rcnt]; obphi[rcnt][bcnt] = sphi[bcnt][rcnt] + dphi[bcnt][rcnt]; if( bcnt < ntpt[rcnt] ){ shiftphiinner[rcnt][bcnt] = phiinner[rcnt][bcnt] - finner[rcnt]->Eval(firstz); shiftphiouter[rcnt][bcnt] = phiouter[rcnt][bcnt] - fouter[rcnt]->Eval(firstz); } } gbinner[rcnt] = new TGraph(nbaf, zpts, ibphi[rcnt]); gbouter[rcnt] = new TGraph(nbaf, zpts, obphi[rcnt]); gsinner[rcnt] = new TGraph(ntpt[rcnt], ztpt[rcnt], shiftphiinner[rcnt]); gsouter[rcnt] = new TGraph(ntpt[rcnt], ztpt[rcnt], shiftphiouter[rcnt]); gbinner[rcnt]->SetLineStyle(2); gbouter[rcnt]->SetLineStyle(2); gbinner[rcnt]->SetLineColor(kBlack); gbouter[rcnt]->SetLineColor(kRed); gsinner[rcnt]->SetMarkerStyle(7); gsouter[rcnt]->SetMarkerStyle(7); gsinner[rcnt]->SetMarkerColor(kBlack); gsouter[rcnt]->SetMarkerColor(kRed); can->cd(rcnt+1); fsinner[rcnt]->Draw("lsame"); fsouter[rcnt]->Draw("lsame"); gsinner[rcnt]->Draw("p"); gsouter[rcnt]->Draw("p"); gbinner[rcnt]->Draw("l"); gbouter[rcnt]->Draw("l"); // Lines for baffle width for( j = 0; j < nbaf; j++ ){ TLine *lbwi = new TLine( zbaf[j]-baf_thick, ibphi[rcnt][j], zbaf[j], ibphi[rcnt][j]); TLine *lbwo = new TLine( zbaf[j]-baf_thick, obphi[rcnt][j], zbaf[j], obphi[rcnt][j]); lbwi->Draw("same"); lbwo->Draw("same"); } // Lines for spacing TLine *spacelow = new TLine( firstz, ibphi[rcnt][0], firstz, ibphi[rcnt][0] + maxsect ); TLine *spacehigh = new TLine( lastz, ibphi[rcnt][nbaf-1], lastz, ibphi[rcnt][nbaf-1] + maxsect ); spacelow->SetLineColor(kBlue-10); spacehigh->SetLineColor(kBlue-10); spacelow->Draw(); spacehigh->Draw(); TLine *seplow = new TLine( firstz, ibphi[rcnt][0] + maxsect - minsep, firstz, ibphi[rcnt][0] + maxsect); TLine *sephigh = new TLine( lastz, ibphi[rcnt][nbaf-1] + maxsect - minsep, lastz, ibphi[rcnt][nbaf-1] + maxsect); seplow->SetLineColor(kGreen+4); sephigh->SetLineColor(kGreen+4); seplow->SetLineWidth(1); sephigh->SetLineWidth(1); seplow->Draw(); sephigh->Draw(); /* Draw Eugene's paths here */ /* printf("block = %d\n", rcnt ); double eugoffset, dummy; for( bcnt = 0; bcnt < 6; bcnt++ ){ GetEugene(0, r[0][rcnt], eugoffset, dummy); GetEugene(bcnt, r[bcnt][rcnt], ein[rcnt][bcnt], eout[rcnt][bcnt] ); printf("\tr = %f z = %f (%f, %f) %f\n", r[bcnt][rcnt], z[bcnt], ein[rcnt][bcnt], eout[rcnt][bcnt], eugoffset); ein[rcnt][bcnt] -= eugoffset; eout[rcnt][bcnt] -= eugoffset; } printf("\n"); TGraph *geugin = new TGraph(6, zbaf, ein[rcnt] ); TGraph *geugout = new TGraph(6, zbaf, eout[rcnt] ); geugin->SetMarkerStyle(21); geugout->SetMarkerStyle(21); geugin->Draw("LP"); geugout->Draw("LP"); */ if( rcnt == 0 ){ new TCanvas(); frm[0]->Draw(); fsinner[rcnt]->Draw("lsame"); fsouter[rcnt]->Draw("lsame"); gsinner[rcnt]->Draw("p"); gsouter[rcnt]->Draw("p"); gbinner[rcnt]->Draw("l"); gbouter[rcnt]->Draw("l"); spacelow->Draw(); spacehigh->Draw(); seplow->Draw(); sephigh->Draw(); for( j = 0; j < nbaf; j++ ){ TLine *lbwi = new TLine( zbaf[j]-baf_thick, ibphi[rcnt][j], zbaf[j], ibphi[rcnt][j]); TLine *lbwo = new TLine( zbaf[j]-baf_thick, obphi[rcnt][j], zbaf[j], obphi[rcnt][j]); lbwi->Draw("same"); lbwo->Draw("same"); } //geugin->Draw("LP"); //geugout->Draw("LP"); } } for( bcnt = 0; bcnt < nbaf; bcnt++ ){ printf("[ "); for( rcnt = 0; rcnt < NBLOCK; rcnt++ ){ printf("%f, %f, %f, %f,\n", r[bcnt][rcnt]*100.0, r[bcnt][rcnt+1]*100.0, sphi[bcnt][rcnt], dphi[bcnt][rcnt]); } printf("],\n"); } printf("( "); for( bcnt = 0; bcnt < nbaf; bcnt++ ){ printf("%5.2f, ", r[bcnt][0]*100.0); } printf(");\n"); printf("( "); for( bcnt = 0; bcnt < nbaf; bcnt++ ){ // printf("%5.2f, ", r[bcnt][NBLOCK-1]*100.0); printf("%5.2f, ", r[bcnt][NBLOCK]*100.0); } printf(");\n"); printf("nslit = %d\n", nslit); return 0; } double getEf( double x, double th ){ return x*MASS_P*beamE/(x*MASS_P + beamE*(1.0-cos(th))); } void GetPhiWidth(TH1F *h, double &min, double &max){ int bin, i; double thresh; if( h->GetEntries() < 100 ){ min = -1; max = -1; return; } // Use integral thresh = 0.02*h->Integral(); /* thresh = 0.02*h->GetMaximum(); for( i = h->GetMaximumBin(); i >= 1; i-- ){ if( h->GetBinContent(i) < thresh ){ min = h->GetBinCenter(i); break; } } */ i = 1; while( h->Integral(1,i) < thresh ){ i++; } min = h->GetBinCenter(i); // Artificially scale this min *= SCALE_MIN_ANG; /* for( i = h->GetMaximumBin(); i < h->GetNbinsX(); i++ ){ if( h->GetBinContent(i) < thresh ){ max = h->GetBinCenter(i); break; } } */ i = h->GetNbinsX(); while( h->Integral(1,i) > h->Integral()-thresh ){ i--; } max = h->GetBinCenter(i); max *= SCALE_MAX_ANG; return; } void GetEugene( int pidx, double r, double &in, double &out ){ // Eugene's Parameterization double baf[6][80] = { { 3.90 , 5.44 , 176.26 , 183.00 , 5.44 , 6.98 , 176.44 , 182.13 , 6.98 , 8.52 , 176.61 , 181.67 , 8.52 , 10.06 , 176.75 , 181.78 , 10.06 , 11.60 , 176.89 , 181.90 , 11.60 , 13.14 , 177.03 , 182.01 , 13.14 , 14.68 , 177.18 , 182.12 , 14.68 , 16.22 , 177.32 , 182.24 , 16.22 , 17.76 , 177.46 , 182.35 , 17.76 , 19.30 , 177.60 , 182.46 , 19.30 , 20.84 , 177.74 , 182.58 , 20.84 , 22.38 , 177.88 , 182.69 , 22.38 , 23.92 , 178.02 , 182.80 , 23.92 , 25.46 , 178.16 , 182.92 , 25.46 , 27.00 , 178.30 , 183.03 , 27.00 , 28.54 , 178.44 , 183.15 , 28.54 , 30.08 , 178.58 , 183.26 , 30.08 , 31.62 , 178.72 , 183.37 , 31.62 , 33.16 , 178.86 , 183.49 , 33.16 , 34.70 , 179.00 , 183.60 }, // #2 { 15.30 , 17.25 , 175.93 , 182.78 , 17.25 , 19.20 , 176.14 , 182.56 , 19.20 , 21.15 , 176.35 , 182.36 , 21.15 , 23.10 , 176.58 , 182.19 , 23.10 , 25.05 , 176.81 , 182.05 , 25.05 , 27.00 , 177.06 , 181.93 , 27.00 , 28.95 , 177.28 , 181.98 , 28.95 , 30.90 , 177.46 , 182.12 , 30.90 , 32.85 , 177.64 , 182.26 , 32.85 , 34.80 , 177.82 , 182.41 , 34.80 , 36.75 , 178.00 , 182.55 , 36.75 , 38.70 , 178.18 , 182.70 , 38.70 , 40.65 , 178.35 , 182.84 , 40.65 , 42.60 , 178.53 , 182.98 , 42.60 , 44.55 , 178.71 , 183.13 , 44.55 , 46.50 , 178.89 , 183.27 , 46.50 , 48.45 , 179.07 , 183.41 , 48.45 , 50.40 , 179.25 , 183.56 , 50.40 , 52.35 , 179.42 , 183.70 , 52.35 , 54.30 , 179.60 , 183.84 }, // #3 plate { 26.60 , 28.96 , 175.71 , 182.49 , 28.96 , 31.32 , 175.95 , 182.40 , 31.32 , 33.68 , 176.20 , 182.32 , 33.68 , 36.04 , 176.46 , 182.27 , 36.04 , 38.40 , 176.73 , 182.23 , 38.40 , 40.76 , 177.02 , 182.22 , 40.76 , 43.12 , 177.31 , 182.22 , 43.12 , 45.48 , 177.62 , 182.23 , 45.48 , 47.84 , 177.93 , 182.28 , 47.84 , 50.20 , 178.15 , 182.46 , 50.20 , 52.56 , 178.36 , 182.63 , 52.56 , 54.92 , 178.58 , 182.81 , 54.92 , 57.28 , 178.80 , 182.98 , 57.28 , 59.64 , 179.01 , 183.16 , 59.64 , 62.00 , 179.23 , 183.33 , 62.00 , 64.36 , 179.45 , 183.51 , 64.36 , 66.72 , 179.67 , 183.68 , 66.72 , 69.08 , 179.88 , 183.85 , 69.08 , 71.44 , 180.10 , 184.03 , 71.44 , 73.80 , 180.32 , 184.21 }, //#4 plate { 37.90 , 40.67 , 175.38 , 182.05 , 40.67 , 43.44 , 175.66 , 182.04 , 43.44 , 46.21 , 175.95 , 182.04 , 46.21 , 48.98 , 176.25 , 182.05 , 48.98 , 51.75 , 176.57 , 182.08 , 51.75 , 54.52 , 176.89 , 182.13 , 54.52 , 57.29 , 177.22 , 182.19 , 57.29 , 60.06 , 177.57 , 182.27 , 60.06 , 62.83 , 177.92 , 182.36 , 62.83 , 65.60 , 178.29 , 182.47 , 65.60 , 68.37 , 178.66 , 182.66 , 68.37 , 71.14 , 178.90 , 182.84 , 71.14 , 73.91 , 179.15 , 183.02 , 73.91 , 76.68 , 179.39 , 183.22 , 76.68 , 79.45 , 179.65 , 183.42 , 79.45 , 82.22 , 179.90 , 183.62 , 82.22 , 84.99 , 180.16 , 183.84 , 84.99 , 87.76 , 180.43 , 184.06 , 87.76 , 90.53 , 180.71 , 184.28 , 90.53 , 93.30 , 180.98 , 184.51 }, //#5 plate { 49.20 , 52.38 , 175.06 , 181.61 , 52.38 , 55.56 , 175.37 , 181.64 , 55.56 , 58.74 , 175.70 , 181.70 , 58.74 , 61.92 , 176.04 , 181.76 , 61.92 , 65.10 , 176.39 , 181.84 , 65.10 , 68.28 , 176.76 , 181.94 , 68.28 , 71.46 , 177.13 , 182.05 , 71.46 , 74.64 , 177.52 , 182.17 , 74.64 , 77.82 , 177.91 , 182.31 , 77.82 , 81.00 , 178.32 , 182.46 , 81.00 , 84.18 , 178.75 , 182.63 , 84.18 , 87.36 , 179.20 , 182.85 , 87.36 , 90.54 , 179.48 , 183.06 , 90.54 , 93.72 , 179.76 , 183.28 , 93.72 , 96.90 , 180.05 , 183.51 , 96.90 ,100.08 , 180.35 , 183.75 , 100.08 ,103.26 , 180.65 , 183.99 , 103.26 ,106.44 , 180.96 , 184.24 , 106.44 ,109.62 , 181.27 , 184.49 , 109.62 ,112.80 , 181.59 , 184.75 }, //#6 plate { 60.40 , 63.98 , 174.92 , 181.35 , 63.98 , 67.56 , 175.27 , 181.43 , 67.56 , 71.14 , 175.63 , 181.52 , 71.14 , 74.72 , 176.00 , 181.63 , 74.72 , 78.30 , 176.38 , 181.75 , 78.30 , 81.88 , 176.78 , 181.88 , 81.88 , 85.46 , 177.18 , 182.02 , 85.46 , 89.04 , 177.60 , 182.18 , 89.04 , 92.62 , 178.03 , 182.35 , 92.62 , 96.20 , 178.47 , 182.53 , 96.20 , 99.78 , 178.92 , 182.73 , 99.78 ,103.36 , 179.39 , 182.94 , 103.36 , 106.94 , 179.88, 183.19, 106.94 , 110.52 , 180.20, 183.43, 110.52 , 114.10 , 180.52, 183.69, 114.10 , 117.68 , 180.84, 183.94, 117.68 , 121.26 , 181.17, 184.21, 121.26 , 124.84 , 181.51, 184.48, 124.84 , 128.42 , 181.86, 184.76, 128.42 , 132.00 , 182.21, 185.04 } }; //// Find R for this int idx = 0; if( r*100.0 < baf[pidx][0] || baf[pidx][19*4+1] <= r*100 ){ in = 0.0; out = 0.0; return; } while( !( baf[pidx][idx*4] <= r*100.0 && r*100.0 < baf[pidx][idx*4+1] ) ) { idx++; } double scale = 0.0; scale = (r*100.0 - baf[pidx][idx])/(baf[pidx][idx+1] - baf[pidx][idx]); double thisin, thisout; thisin = scale*( baf[pidx][(idx+1)*4+2]-baf[pidx][(idx)*4+2] ) + baf[pidx][(idx)*4+2]; thisout = scale*( baf[pidx][(idx+1)*4+3]-baf[pidx][(idx)*4+3] ) + baf[pidx][(idx)*4+3]; // These are the block positions. The slits are complementary to them in = thisout; out = (360.0/30.0 + thisin); return; }