loon: ProtonDist.cxx Source File

00001 #include "ProtonDist.h"
00002 
00003 #include <HistMan/HistMan.h>
00004 
00005 #include <BeamData/ana/bmnt/AcnetDevice.h>
00006 #include <BeamData/ana/bmnt/Knot.h>
00007 #include <BeamData/ana/bmnt/ProfMon.h>
00008 
00009 #include <TH2D.h>
00010 #include <TH1D.h>
00011 
00012 #include <iostream>
00013 #include <cmath>
00014 using namespace std;
00015 
00016 /* NearDet Runs with stable beam.  pseudo medium energy, ~2.5e12 POTs,
00017    263 spills once per minute.
00018 
00019    6067_0000 1106373930 1106355930 utc
00020    6067_0001
00021    6068_0000
00022    6068_0001
00023    6068_0002
00024    6068_0003 1106438416 1106420416
00025 
00026 */
00027 
00028 
00029 
00030 void ProtonDist::Book(PlotterManager& )
00031 {
00032     HistMan hm("prot");
00033 
00034     const int npos_bins = 400;
00035     const double pos_edge = 5;
00036 
00037 
00038     hm.Book<TH1D>("v1_bpm","Proton beam vertical position at E:VP121",
00039                   npos_bins,-pos_edge,pos_edge);
00040     hm.Book<TH1D>("h1_bpm","Proton beam horizontal position at E:HP121",
00041                   npos_bins,-pos_edge,pos_edge);
00042     hm.Book<TH1D>("v2_bpm","Proton beam vertical position at E:VPTGT",
00043                   npos_bins,-pos_edge,pos_edge);
00044     hm.Book<TH1D>("h2_bpm","Proton beam horizontal position at E:HPTGT",
00045                   npos_bins,-pos_edge,pos_edge);    
00046 
00047     hm.Book<TH1D>("h2_bpm_1","Proton beam horizontal position at E:HPTGT, run 6067",
00048                   npos_bins,-pos_edge,pos_edge);    
00049     hm.Book<TH1D>("h2_bpm_2","Proton beam horizontal position at E:HPTGT, run 6068",
00050                   npos_bins,-pos_edge,pos_edge);    
00051 
00052     hm.Book<TH1D>("v1_pmm","Proton beam vertical position at E:M121DS",
00053                   npos_bins,-pos_edge,pos_edge);
00054     hm.Book<TH1D>("h1_pmm","Proton beam horizontal position at E:M121DS",
00055                   npos_bins,-pos_edge,pos_edge);
00056     hm.Book<TH1D>("v2_pmm","Proton beam vertical position at E:MTGTDS",
00057                   npos_bins,-pos_edge,pos_edge);
00058     hm.Book<TH1D>("h2_pmm","Proton beam horizontal position at E:MTGTDS",
00059                   npos_bins,-pos_edge,pos_edge);
00060 
00061     hm.Book<TH1D>("h2_pmm_1","Proton beam horizontal position at E:MTGTDS, run 6067",
00062                   npos_bins,-pos_edge,pos_edge);
00063     hm.Book<TH1D>("h2_pmm_2","Proton beam horizontal position at E:MTGTDS, run 6067",
00064                   npos_bins,-pos_edge,pos_edge);
00065     
00066     hm.Book<TH1D>("v1_bpm_pmm","Proton beam horizontal difference E:VP121 - E:M121DS",
00067                   npos_bins,-pos_edge,pos_edge);
00068     hm.Book<TH1D>("h1_bpm_pmm","Proton beam horizontal difference E:HP121 - E:M121DS",
00069                   npos_bins,-pos_edge,pos_edge);
00070     hm.Book<TH1D>("v2_bpm_pmm","Proton beam horizontal difference E:VPTGT - E:MTGTDS",
00071                   npos_bins,-pos_edge,pos_edge);
00072     hm.Book<TH1D>("h2_bpm_pmm","Proton beam horizontal difference E:HPTGT - E:MTGTDS",
00073                   npos_bins,-pos_edge,pos_edge);
00074 
00075 
00076     const int nwid_bins = 300;
00077     const double wid_max = 3;
00078 
00079     hm.Book<TH1D>("v1_pms","Proton beam vertical width at E:M121DS",
00080                   nwid_bins,0,wid_max);
00081     hm.Book<TH1D>("h1_pms","Proton beam horizontal width at E:M121DS",
00082                   nwid_bins,0,wid_max);
00083     hm.Book<TH1D>("v2_pms","Proton beam vertical width at E:MTGTDS",
00084                   nwid_bins,0,wid_max);
00085     hm.Book<TH1D>("h2_pms","Proton beam horizontal width at E:MTGTDS",
00086                   nwid_bins,0,wid_max);
00087 
00088     hm.Book<TH1D>("h2_pms_1","Proton beam horizontal width at E:MTGTDS, run 6067",
00089                   nwid_bins,0,wid_max);
00090     hm.Book<TH1D>("h2_pms_2","Proton beam horizontal width at E:MTGTDS, run 6068",
00091                   nwid_bins,0,wid_max);
00092     
00093 
00094     hm.Book<TH2D>("bpm_pos","Proton Beam Center at target, measured by BPM",
00095                   80,-3,1,40,0,2);
00096     hm.Book<TH2D>("bpm_ang","Proton Beam dy/dz vs. dx/dz at target (rad), measured by BPM",
00097                   100,-100e-6,100e-6,100,-100e-6,100e-6);
00098 
00099     hm.Book<TH2D>("pm_pos","Proton Beam Center at target, mearused by ProfMon",
00100                   80,-3,1,40,0,2);
00101     //hm.Book<TH2D>("pm_wid","Proton Beam Width at target, measured by ProfMon",
00102 
00103     hm.Book<TH2D>("pm_ang","Proton Beam dy/dz vs. dx/dz at target (rad), measured by ProfMon",
00104                   100,-100e-6,100e-6,100,-100e-6,100e-6);
00105 
00106     hm.Book<TH2D>("bpm_pm_pos_diff","Difference of Proton Beam Center at target, BPM-ProfMon",
00107                   40,-1,1,40,-1,1);
00108     hm.Book<TH2D>("bpm_pm_pos_diff_1","Difference of Proton Beam Center at target, BPM-ProfMon, Run 6067",
00109                   40,-1,1,40,-1,1);
00110     hm.Book<TH2D>("bpm_pm_pos_diff_2","Difference of Proton Beam Center at target, BPM-ProfMon, Run 6068",
00111                   40,-1,1,40,-1,1);
00112     
00113     hm.Book<TH1D>("pot","Proton intensity at E:TORTGT (x10^12)",
00114                   100,0,10);
00115     hm.Book<TH1D>("pot_1","Run 6067 Proton intensity at E:TORTGT (x10^12)",
00116                   100,0,10);
00117     hm.Book<TH1D>("pot_2","Run 6068 Proton intensity at E:TORTGT (x10^12)",
00118                   100,0,10);
00119 
00120 }
00121 static bool run_6067(double ts)
00122 {
00123     const int ti_6067 = 1106355930-1; // 2005-01-22 01:05:30 UCT
00124     const int tf_6067 = 1106362789+1; // 2005-01-22 02:59:49 UTC
00125 
00126     return ts>ti_6067 && ts<tf_6067;
00127 }
00128 
00129 static bool run_6068(double ts)
00130 {
00131     const int ti_6068 = 1106409473-1; // 2005-01-22 15:57:53 UTC
00132     const int tf_6068 = 1106420416+1; // 2005-01-22 19:00:16 UTC
00133 
00134     return ts>ti_6068 && ts<tf_6068;
00135 }
00136 
00137 // Extrapolate from point 1 to 2 to point 3.  Returns the transverse
00138 // position at point 3 given longitudinal on.
00139 static double extrapolate_position(double t1, double z1, double t2, double z2, double z3)
00140 {
00141     return t1+(t2-t1)*(z3-z1)/(z2-z1);
00142 }
00143 // Return the 2d angle (atan(dx/dz)) in rad between the two points
00144 static double extrapolate_direction(double t1, double z1, double t2, double z2)
00145 {
00146     return atan2(t2-t1,z2-z1);
00147 }
00148 
00149 bool ProtonDist::Fill(PlotterManager& pm)
00150 {
00151     Knot& k = pm.GetKnot();
00152     
00153     const AcnetDevice* vp121 = k.GetDevice("E:VP121");
00154     const AcnetDevice* hp121 = k.GetDevice("E:HP121");
00155     const AcnetDevice* vptgt = k.GetDevice("E:VPTGT");
00156     const AcnetDevice* hptgt = k.GetDevice("E:HPTGT");
00157 
00158     const AcnetDevice* pm121 = k.GetDevice("E:M121DS");
00159     const AcnetDevice* pmtgt = k.GetDevice("E:MTGTDS");
00160 
00161     const AcnetDevice* tortgt = k.GetDevice("E:TORTGT");
00162 
00163     int entry_number = k.GetEntryNumber();
00164 
00165     bool ok = vp121 && hp121 && vptgt && pm121 && pmtgt && tortgt;
00166     if (!ok) {
00167         cerr << "Failed to get all acnet devices in entry " << entry_number << endl;
00168         return false;
00169     }
00170     ok = vp121->ndata && hp121->ndata && vptgt->ndata && hptgt->ndata &&
00171         pm121->ndata && pmtgt->ndata && tortgt->ndata;
00172     if (!ok) {
00173         cerr << "Not all devices have data in entry " << entry_number << endl;
00174         return false;
00175     }
00176 
00177     HistMan hm("prot");
00178 
00179     double pot = tortgt->data[0];
00180     hm.Fill1d("pot",pot);
00181 
00182     // "as found" offsets from center. (Wes Smart). 
00183 
00184     double bpm1o[2] = { 0 };
00185     double bpm2o[2] = { 0 };
00186     double pmm1o[2] = { 0 };
00187     double pmm2o[2] = { 0 };
00188     if (do_correction) {
00189         cerr << "Doing correction\n";
00190         bpm1o[0] = +0.000667*25.4;
00191         bpm1o[1] = -0.002000*25.4;
00192 
00193         bpm2o[0] = -0.001800*25.4;
00194         bpm2o[1] = -0.005667*25.4;
00195 
00196         pmm1o[0] = -0.011500*25.4;
00197         pmm1o[1] = +0.001000*25.4;
00198 
00199         pmm2o[0] = +0.005250*25.4;
00200         pmm2o[1] = -0.001000*25.4;
00201     }
00202 
00203     double bpm1[2] = { hp121->data[0]+bpm1o[0], vp121->data[0]+bpm1o[1] };
00204     double bpm2[2] = { hptgt->data[0]+bpm2o[0], vptgt->data[0]+bpm2o[1] };
00205     hm.Fill1d("h1_bpm",bpm1[0]);
00206     hm.Fill1d("v1_bpm",bpm1[1]);
00207     hm.Fill1d("h2_bpm",bpm2[0]);
00208     hm.Fill1d("v2_bpm",bpm2[1]);
00209 
00210     double pmm1[2], pmm2[2], pms1[2], pms2[2]; // (P)rofile (M)ontor (M)ean/(S)igma
00211 
00212     ProfMon prof121;
00213     prof121.SetPeds(pm.GetPeds().GetPeds("E:M121DS"));
00214     prof121.SetMask(pm.GetPeds().GetMask("E:M121DS"));
00215     prof121.SetAcnetDevice(*pm121);
00216     prof121.GetStats(pmm1[0],pmm1[1],pms1[0],pms1[1]);
00217     for (int ind=0; ind<2; ++ind) {
00218         pmm1[ind] += pmm1o[ind];
00219         pmm2[ind] += pmm2o[ind];
00220     }
00221 
00222     ProfMon proftgt;
00223     proftgt.SetPeds(pm.GetPeds().GetPeds("E:MTGTDS"));
00224     proftgt.SetMask(pm.GetPeds().GetMask("E:MTGTDS"));
00225     proftgt.SetAcnetDevice(*pmtgt);
00226     proftgt.GetStats(pmm2[0],pmm2[1],pms2[0],pms2[1]);
00227 
00228     cerr << pmm2[0] << " " << pmm2[1] << " " << pms2[0] << " " << pms2[1] << endl;
00229 
00230     hm.Fill1d("h1_pmm",pmm1[0]);
00231     hm.Fill1d("v1_pmm",pmm1[1]);
00232     hm.Fill1d("h2_pmm",pmm2[0]);
00233     hm.Fill1d("v2_pmm",pmm2[1]);
00234 
00235     hm.Fill1d("h1_pms",pms1[0]);
00236     hm.Fill1d("v1_pms",pms1[1]);
00237     hm.Fill1d("h2_pms",pms2[0]);
00238     hm.Fill1d("v2_pms",pms2[1]);
00239 
00240     hm.Fill1d("h1_bpm_pmm",bpm1[0]-pmm1[0]);
00241     hm.Fill1d("v1_bpm_pmm",bpm1[1]-pmm1[1]);
00242     hm.Fill1d("h2_bpm_pmm",bpm2[0]-pmm2[0]);
00243     hm.Fill1d("v2_bpm_pmm",bpm2[1]-pmm2[1]);
00244 
00245 
00246     HistMan ps_hm("spill");     // Per Spill HistMan
00247     TH1D* m121 = ps_hm.Book<TH1D>(Form("em121ds_%04d",entry_number),
00248                                   Form("E:M121DS spill %d (t=%10d)",
00249                                        entry_number,prof121.GetVmeSeconds()),
00250                                   96,0,96);
00251     TH1D* mtgt = ps_hm.Book<TH1D>(Form("emtgtds_%04d",entry_number),
00252                                   Form("E:MTGTDS spill %d (t=%10d)",
00253                                        entry_number,proftgt.GetVmeSeconds()),
00254                                   96,0,96);
00255     for (int ind=0; ind<96; ++ind) {
00256         m121->Fill(ind,prof121.GetVoltage(ind));
00257         mtgt->Fill(ind,proftgt.GetVoltage(ind));
00258     }
00259 
00260 
00261     // Values taken from Wes Smart's "as found" survey.  In feet.
00262     const double
00263         zh1_bpm =    -70.4923,  // HP121
00264         zv1_bpm =    -69.5756,  // VP121
00265         z1_pm   =    -68.7882,  // PM121
00266         zh2_bpm =    -29.6521,  // HPTGT
00267         zv2_bpm =    -28.7354,  // VPTGT
00268         z2_pm   =    -27.9479,  // PMTGT
00269         z_target_pme = 0.0; // 
00270 
00271 
00272 
00273     //double extrapolate_position(double t1, double z1, double t2, double z2, double z3)
00274 
00275     double th_bpm = extrapolate_position(bpm1[0],zh1_bpm,bpm2[0],zh2_bpm,z_target_pme);
00276     double tv_bpm = extrapolate_position(bpm1[1],zv1_bpm,bpm2[1],zv2_bpm,z_target_pme);
00277     double ah_bpm = extrapolate_direction(bpm1[0],zh1_bpm,bpm2[0],zh2_bpm);
00278     double av_bpm = extrapolate_direction(bpm1[1],zv1_bpm,bpm2[1],zv2_bpm);
00279 
00280     double th_pm = extrapolate_position(pmm1[0],z1_pm,pmm2[0],z2_pm,z_target_pme);
00281     double tv_pm = extrapolate_position(pmm1[1],z1_pm,pmm2[1],z2_pm,z_target_pme);
00282     double ah_pm = extrapolate_direction(pmm1[0],z1_pm,pmm2[0],z2_pm);
00283     double av_pm = extrapolate_direction(pmm1[1],z1_pm,pmm2[1],z2_pm);
00284 
00285 
00286     hm.Fill2d("bpm_pos",th_bpm,tv_bpm,pot);
00287     hm.Fill2d("bpm_ang",ah_bpm,av_bpm,pot);
00288 
00289     hm.Fill2d("pm_pos",th_pm,tv_pm,pot);
00290     hm.Fill2d("pm_ang",ah_pm,av_pm,pot);
00291     //hm.Fill2d("pm_wid");
00292 
00293     hm.Fill2d("bpm_pm_pos_diff",th_bpm-th_pm,tv_bpm-tv_pm,pot);
00294     
00295     if (run_6067(tortgt->timestamp)) {
00296         hm.Fill1d("h2_bpm_1",bpm2[0]);
00297         hm.Fill1d("h2_pmm_1",pmm2[0]);
00298         hm.Fill1d("h2_pms_1",pms2[0]);
00299         hm.Fill2d("bpm_pm_pos_diff_1",th_bpm-th_pm,tv_bpm-tv_pm,pot);
00300         hm.Fill1d("pot_1",pot);
00301     }
00302     if (run_6068(tortgt->timestamp)) {
00303         hm.Fill1d("h2_bpm_2",bpm2[0]);
00304         hm.Fill1d("h2_pmm_2",pmm2[0]);
00305         hm.Fill1d("h2_pms_2",pms2[0]);
00306         hm.Fill2d("bpm_pm_pos_diff_2",th_bpm-th_pm,tv_bpm-tv_pm,pot);
00307         hm.Fill1d("pot_2",pot);
00308     }
00309 
00310     return true;
00311 }