BDProfMon Class Reference
[BeamDataUtil]

A BDSwicDevice specific for Profile Monitor data. More...

#include <BDProfMon.h>

Inheritance diagram for BDProfMon:

BDSwicDevice RawBeamSwicData List of all members.

Public Member Functions

 BDProfMon ()
 BDProfMon (const RawBeamData &swic_data)
virtual ~BDProfMon ()
virtual void SetData (const RawBeamData &swic_data)
 Set the RawBeamData from a swic scanner device.
virtual bool SetMask (const std::vector< double > &mask)
 Mask off any bad channels.
double WirePosition (int channel) const
 Return the distance from the origin to the wire on the given profile monitor channel [1-44].
double GetStats (double &xmean, double &ymean, double &xrms, double &yrms) const
 Calculate statistical values.
double GetGaussFit (double &xmean, double &ymean, double &xsigma, double &ysigma) const
 Calculate values from a Gaussian fit.

Static Public Member Functions

static int Xindex (int channel)
 Convert a horizontal channel [1-44] to an index in to the basic SWIC channels (0-95).
static int Yindex (int channel)
 Convert a vertical channel [1-44] to an index in to the basic SWIC channels (0-95).

Private Member Functions

virtual int Index (const int channel)
 Do not use. Use Xindex() or Yindex() instead.
void SuppressDeadHot (TGraphErrors *prof) const
void UpdateCache ()
double GetStats (TGraphErrors *prof, double &mean, double &rms) const

Private Attributes

double fSpacing

Detailed Description

A BDSwicDevice specific for Profile Monitor data.

This provides access to the underlying SWIC data assuming it came from a profile monitor.

SWIC scanner have 96 channels. For profile monitors these are partitioned into a horizontal and vertical section, each with 1/2 the channels. Of these 48 only 44 channels are used.

Author:
(last to touch it)
Author
bv
Version:
Revision
1.4
Date:
Date
2005/12/14 21:03:20
Contact: bv@bnl.gov

Created on: Fri Apr 15 13:09:38 2005

Id
BDProfMon.h,v 1.4 2005/12/14 21:03:20 bv Exp

Definition at line 40 of file BDProfMon.h.

Constructor & Destructor Documentation

BDProfMon::BDProfMon (  ) 

Definition at line 17 of file BDProfMon.cxx.

00018     : BDSwicDevice(), fSpacing(0)
00019 {
00020 }

BDProfMon::BDProfMon ( const RawBeamData &  swic_data  ) 

Definition at line 21 of file BDProfMon.cxx.

References UpdateCache().

00022     : BDSwicDevice(swic_data), fSpacing(0)
00023 {
00024     this->UpdateCache();
00025 }

BDProfMon::~BDProfMon (  )  [virtual]

Definition at line 26 of file BDProfMon.cxx.

00027 {
00028 }

Member Function Documentation

double BDProfMon::GetGaussFit ( double &  xmean,
double &  ymean,
double &  xsigma,
double &  ysigma 
) const

Calculate values from a Gaussian fit.

Statistics are the mean and sigma from a fit to a single Gaussian

Returns:
a double giving total fit area in both X and Y views.
Units are in Munits.

Definition at line 152 of file BDProfMon.cxx.

References fSpacing, BDSwicDevice::GetNoise(), GetStats(), BDSwicDevice::GetVoltage(), Msg::kVerbose, Munits::millivolt, Munits::mm, MSG, SuppressDeadHot(), Xindex(), and Yindex().

Referenced by BDTarget::ProfileProjection().

00153 {
00154     TGraphErrors xProf,yProf;
00155   
00156     xProf.Clear("");
00157     yProf.Clear("");
00158     int nxpoints = 0, nypoints = 0;
00159   
00160     for (int ch = 1; ch <= 44; ++ch) {
00161     
00162         double X = this->WirePosition(ch);
00163         double Y = X;
00164     
00165         int x_index = this->Xindex(ch);
00166         int y_index = this->Yindex(ch);
00167     
00168         double Qx = this->GetVoltage(x_index);
00169         double Qy = this->GetVoltage(y_index);
00170 
00171     
00172         Qx = Qx > 0 ? Qx : 0;
00173         Qy = Qy > 0 ? Qy : 0;
00174     
00175     
00176         double xNoise = this->GetNoise(x_index);
00177         double yNoise = this->GetNoise(y_index);
00178 
00179         xNoise = xNoise > 0 ? xNoise : 1.0*Munits::millivolt;
00180         yNoise = yNoise > 0 ? yNoise : 1.0*Munits::millivolt;
00181 
00182         // The reason I use mm as the units to fit is because MINUIT doesnt like small numbers
00183         // when fitting - M.B.
00184 
00185         // Reminder to myself:
00186         // The silly profile is not actually a Gaussian, by setting the x axis error
00187         // to zero, the fit seems to do a better job of fitting the core
00188         // - M.B.
00189 
00190 
00191         xProf.SetPoint(nxpoints,X/Munits::mm,Qx/Munits::millivolt);
00192         //      xProf.SetPointError(nxpoints,fSpacing/(2.0*Munits::mm),xNoise/Munits::millivolt);
00193         xProf.SetPointError(nxpoints,0.0,xNoise/Munits::millivolt);
00194         nxpoints++;
00195 
00196         yProf.SetPoint(nypoints,Y/Munits::mm,Qy/Munits::millivolt);
00197         //      yProf.SetPointError(nypoints,fSpacing/(2.0*Munits::mm),yNoise/Munits::millivolt);
00198         yProf.SetPointError(nypoints,0.0,yNoise/Munits::millivolt);
00199         nypoints++;
00200 
00201     }
00202 
00203     // Remove dead and hot channels from the profiles being fitted.
00204     SuppressDeadHot(&xProf);
00205     SuppressDeadHot(&yProf);
00206 
00207     // TF1 *gausF = new TF1("gausF","[0]*TMath::Gaus(x,[1],[2],1)+[3]+[4]*x");
00208     // For some reason this is faster than if I use a TMath::Gaus:
00209     //  TF1 *gausF = new TF1("gausF","[0]*exp((-1*(x-[1])*(x-[1]))/(2*[2]*[2]))/(sqrt(2*3.142)*[2])+[3]+[4]*x");
00210     TF1 gausF("gausF","[0]*exp((-1*(x-[1])*(x-[1]))/(2*[2]*[2]))/(sqrt(2*3.142)*[2])");
00211 
00212     gausF.SetParName(0,"Area");
00213     gausF.SetParName(1,"Mean");
00214     gausF.SetParName(2,"Sigma");
00215     //  gausF.SetParName(3,"norm");
00216     //  gausF.SetParName(3,"slope");
00217   
00218     gausF.SetParameter(0,2000);
00219     gausF.SetParameter(1,0.0);
00220     gausF.SetParameter(2,0.5);
00221     //  gausF.SetParameter(3,0.5);
00222     //  gausF.SetParameter(4,0.0);
00223     gausF.SetParLimits(0,0.00,100000000);
00224     gausF.SetParLimits(2,0.0001,1000);
00225 
00226     // xProf.Draw("AP");
00227     // gSystem->ProcessEvents();
00228 
00229     double xarea,yarea;
00230     double xStatMean,yStatMean,xrms,yrms;
00231   
00232     xarea  = -999;
00233     xmean  = -999;
00234     xsigma = -999;
00235     xarea = -999;
00236     ymean  = -999;
00237     ysigma = -999;
00238 
00239     // get the statistical mean and rms of the profile distribution
00240     //
00241     xarea = GetStats(&xProf,xStatMean,xrms);
00242     xarea*=Munits::millivolt;
00243   
00244     MSG("BD",Msg::kVerbose) << "X profile data" << endl;
00245     for (int i = 0; i<xProf.GetN(); i++)
00246         MSG("BD",Msg::kVerbose) << "x[" << i << "] = " <<xProf.GetX()[i] 
00247                                 << " y[" << i << "] = " <<xProf.GetY()[i] 
00248                                 << " ex[" << i << "] = " <<xProf.GetErrorX(i)
00249                                 << " ey[" << i << "] = " <<xProf.GetErrorY(i)
00250                                 << endl;
00251       
00252 
00253     //  xProf.Print("");
00254     MSG("BD",Msg::kVerbose) << endl;
00255     MSG("BD",Msg::kVerbose) << "X charge (V) : " << xarea << endl;
00256     MSG("BD",Msg::kVerbose) << "X stat mean  : " << xStatMean << endl;
00257     MSG("BD",Msg::kVerbose) << "X rms  : " << xrms << endl;
00258 
00259     
00260     // Check we have enough points and enough charge for a reasonable fit 
00261   
00262     if (xProf.GetN() > 5 &&
00263         xarea > 100*Munits::millivolt &&
00264         fabs(xStatMean) < 20*fSpacing/Munits::mm &&
00265         fabs(xrms) < 11*fSpacing/Munits::mm) {
00266         gausF.SetParameter(0,xarea/Munits::millivolt);
00267         gausF.SetParameter(1,xStatMean);
00268         gausF.SetParameter(2,xrms*0.6);
00269         gausF.SetParameter(3,0.5);
00270         gausF.SetParameter(4,0.0);
00271         if (xProf.Fit("gausF","q","",-22*fSpacing/Munits::mm,22*fSpacing/Munits::mm) == 0) {
00272             gMinuit->Command("set str 2");
00273             // Check the gaussian values are reasonable then 
00274             // run a more accurate minimization with MINOS and
00275             // strategy 2      
00276             bool goodfit = true;
00277 
00278             // check if fit results are physical
00279             if (gausF.GetParameter(0) < 10)
00280                 goodfit = false;
00281             if (fabs(gausF.GetParameter(1)) > 22*fSpacing/Munits::mm)
00282                 goodfit = false;
00283             if (gausF.GetParameter(2) > 22*fSpacing/Munits::mm ||
00284                 gausF.GetParameter(2) < 0.0)
00285                 goodfit = false;
00286       
00287             if (goodfit) {
00288                 xarea = gausF.GetParameter(0)/fSpacing*Munits::mm*Munits::millivolt;
00289                 xmean = gausF.GetParameter(1)*Munits::mm;
00290                 xsigma = gausF.GetParameter(2)*Munits::mm;
00291                 if (xProf.Fit("gausF","q E M","",-22*fSpacing/Munits::mm,22*fSpacing/Munits::mm) == 0){
00292                     gMinuit->Command("set str 1");
00293                     // remember to correct the fit result for area
00294                     // by dividing by the "bin size" and correcting the units
00295                     xarea = gausF.GetParameter(0)/fSpacing*Munits::mm*Munits::millivolt;
00296                     xmean = gausF.GetParameter(1)*Munits::mm;
00297                     xsigma = gausF.GetParameter(2)*Munits::mm;
00298                 }
00299             }
00300             // gausF.Draw("same");
00301             // gSystem->ProcessEvents();
00302         }
00303     }
00304   
00305     // yProf.Draw("AP");
00306     // gSystem->ProcessEvents();
00307 
00308 
00309     yarea = GetStats(&yProf,yStatMean,yrms);
00310     yarea*=Munits::millivolt;
00311 
00312   
00313     MSG("BD",Msg::kVerbose) << "Y profile data" << endl;
00314     for (int i = 0; i<yProf.GetN(); i++)
00315         MSG("BD",Msg::kVerbose) << "x[" << i << "] = " <<yProf.GetX()[i] 
00316                                 << " y[" << i << "] = " <<yProf.GetY()[i] 
00317                                 << " ex[" << i << "] = " <<yProf.GetErrorX(i)
00318                                 << " ey[" << i << "] = " <<yProf.GetErrorY(i)
00319                                 << endl;
00320       
00321     //  yProf.Print("");
00322     MSG("BD",Msg::kVerbose) << endl;     
00323     MSG("BD",Msg::kVerbose) << "Y charge (V) : " << yarea << endl;
00324     MSG("BD",Msg::kVerbose) << "Y stat mean  : " << yStatMean << endl;
00325     MSG("BD",Msg::kVerbose) << "Y rms  : " << yrms << endl;
00326 
00327   
00328     if (yProf.GetN() > 5 &&
00329         yarea > 100*Munits::millivolt &&
00330         fabs(yStatMean) < 20*fSpacing/Munits::mm &&
00331         fabs(yrms) < 11*fSpacing/Munits::mm) {
00332     
00333         gausF.SetParameter(0,yarea/Munits::millivolt);
00334         gausF.SetParameter(1,yStatMean);
00335         gausF.SetParameter(2,yrms*0.6);
00336         gausF.SetParameter(3,0.5);
00337         gausF.SetParameter(4,0.0);
00338    
00339         if (yProf.Fit("gausF","q","",-22*fSpacing/Munits::mm,22*fSpacing/Munits::mm) == 0) {
00340             gMinuit->Command("set str 2");
00341             // Check the gaussian values are reasonable then 
00342             // run a more accurate minimization with MINOS and
00343             // strategy 2
00344             bool goodfit=true;
00345 
00346             // check if fit results are physical
00347             if (gausF.GetParameter(0) < 10)
00348                 goodfit = false;
00349             if (fabs(gausF.GetParameter(1)) > 22*fSpacing/Munits::mm)
00350                 goodfit = false;
00351             if (gausF.GetParameter(2) > 22*fSpacing/Munits::mm ||
00352                 gausF.GetParameter(2) < 0.0)
00353                 goodfit = false;
00354             if (goodfit) {
00355                 yarea = gausF.GetParameter(0)/fSpacing*Munits::mm*Munits::millivolt;
00356                 ymean = gausF.GetParameter(1)*Munits::mm;
00357                 ysigma = gausF.GetParameter(2)*Munits::mm;
00358                 // run a more accurate minimization with MINOS      
00359                 if (yProf.Fit("gausF","q E M","",-22*fSpacing/Munits::mm,22*fSpacing/Munits::mm) == 0){
00360                     gMinuit->Command("set str 1");
00361                     yarea = gausF.GetParameter(0)/fSpacing*Munits::mm*Munits::millivolt;
00362                     ymean = gausF.GetParameter(1)*Munits::mm;
00363                     ysigma = gausF.GetParameter(2)*Munits::mm;
00364                 }
00365             }
00366             // gausF.Draw("same");
00367             // gSystem->ProcessEvents();
00368         }
00369     }
00370   
00371     //  string name = this->GetData().GetName();
00372 
00373     // get the statistics from all channels
00374     // and total charge on the profile monitors
00375     // to use as as a check
00376     double totarea=0;
00377 
00378     totarea=this->GetStats(xStatMean,yStatMean,xrms,yrms);
00379 
00380     MSG("BD",Msg::kVerbose) << "Total area "
00381                             << totarea 
00382                             << " Prof Means = "
00383                             << xStatMean << ","
00384                             << yStatMean
00385                             << " RMS = "
00386                             << xrms << ","
00387                             << yrms << endl;
00388 
00389   
00390     MSG("BD",Msg::kVerbose) << "Fitted area "
00391                             << xarea+yarea 
00392                             << " Gauss Means = "
00393                             << xmean << ","
00394                             << ymean
00395                             << " Gaus sigmas = "
00396                             << xsigma << ","
00397                             << ysigma << endl;
00398     
00399     return xarea+yarea;
00400 }

double BDProfMon::GetStats ( TGraphErrors *  prof,
double &  mean,
double &  rms 
) const [private]

Definition at line 128 of file BDProfMon.cxx.

00128                                                                               {
00129 
00130     double qx=0,qx2=0,qtotx=0;
00131 
00132     for (int i = 0; i < prof->GetN(); i++) {
00133         double X = prof->GetX()[i];
00134         double Qx = prof->GetY()[i];
00135         qtotx += Qx;
00136         qx += Qx*X;
00137         qx2 += Qx*X*X;
00138     }
00139     if (qtotx == 0.0) {
00140         mean = 0;
00141         rms = 0;
00142     }
00143     else {
00144         mean = qx/qtotx;
00145         double tmp = qx2/qtotx-mean*mean;
00146         if (tmp>0.0) rms = sqrt(tmp);
00147         else rms=0;
00148     }
00149     return qtotx;       
00150 }

double BDProfMon::GetStats ( double &  xmean,
double &  ymean,
double &  xrms,
double &  yrms 
) const

Calculate statistical values.

Statistics are mean and RMS using a charge weighted distribution.

Returns:
a double giving total intensity in both X and Y views.
Units are in Munits.

Definition at line 78 of file BDProfMon.cxx.

References BDSwicDevice::GetVoltage(), Xindex(), and Yindex().

Referenced by GetGaussFit().

00079 {
00080     double qx=0,qx2=0,qtotx=0;
00081     double qy=0,qy2=0,qtoty=0;
00082     for (int ch = 1; ch <= 44; ++ch) {
00083         double X = this->WirePosition(ch);
00084         double Y = X;
00085 
00086         int x_index = this->Xindex(ch);
00087         int y_index = this->Yindex(ch);
00088 
00089         double Qx = this->GetVoltage(x_index);
00090         double Qy = this->GetVoltage(y_index);
00091 
00092         Qx = Qx > 0 ? Qx : 0;
00093         Qy = Qy > 0 ? Qy : 0;
00094 
00095         qtotx += Qx;
00096         qtoty += Qy;
00097 
00098         qx += Qx*X;
00099         qy += Qy*Y;
00100 
00101         qx2 += Qx*X*X;
00102         qy2 += Qy*Y*Y;
00103     }
00104     if (qtotx == 0.0) {
00105         xmean = 0;
00106         xrms = 0;
00107     }
00108     else {
00109         xmean = qx/qtotx;
00110         double tmp = qx2/qtotx-xmean*xmean;
00111         if (tmp>0.0) xrms = sqrt(tmp);
00112         else xrms=0;
00113     }
00114     if (qtoty == 0.0) {
00115         ymean = 0;
00116         yrms = 0;
00117     }
00118     else {
00119         ymean = qy/qtoty;
00120         double tmp = qy2/qtoty-ymean*ymean;
00121         if (tmp > 0.0) yrms = sqrt(tmp);
00122         else yrms = 0;
00123     }
00124     return qtotx+qtoty;
00125 }

int BDProfMon::Index ( const int  channel  )  [private, virtual]

Do not use. Use Xindex() or Yindex() instead.

Definition at line 427 of file BDProfMon.cxx.

00427 { const int do_not_call=0; assert(do_not_call); }

void BDProfMon::SetData ( const RawBeamData &  swic_data  )  [virtual]

Set the RawBeamData from a swic scanner device.

Reimplemented from BDSwicDevice.

Definition at line 41 of file BDProfMon.cxx.

References BDSwicDevice::SetData(), and UpdateCache().

Referenced by BDTarget::SetSpill().

00042 {
00043     this->BDSwicDevice::SetData(swic_data);
00044     this->UpdateCache();
00045 }

bool BDProfMon::SetMask ( const std::vector< double > &  mask  )  [virtual]

Mask off any bad channels.

The mask is defined on basic SWIC channels (0-95). You do not have to worry about masking unused channels, that will be done internally.

Reimplemented from BDSwicDevice.

Definition at line 56 of file BDProfMon.cxx.

References BDSwicDevice::SetMask().

00057 {
00058     vector<double> mask = cmask;
00059 
00060     const int unused[] = {0,1, 46,47, 48,49, 95, 95, -1} ;
00061     for (int ind=0; unused[ind] >= 0; ++ind) mask[ind] = 0.0;
00062     return this->BDSwicDevice::SetMask(mask);
00063 }

void BDProfMon::SuppressDeadHot ( TGraphErrors *  prof  )  const [private]

Definition at line 402 of file BDProfMon.cxx.

Referenced by GetGaussFit().

00403 {
00404 
00405     // remove hot channels at the edges
00406     if (prof->GetN() > 1) {
00407         if  (prof->GetY()[0] > 10*prof->GetY()[1]) 
00408             prof->RemovePoint(0);
00409         if (prof->GetY()[prof->GetN()-1] > 10*prof->GetY()[prof->GetN()-2]) 
00410             prof->RemovePoint(prof->GetN()-1);
00411     }
00412 
00413     // suppress dead/hot channels by assuming its a smoothly varying distribution
00414   
00415     for (int i = 1; i < prof->GetN()-1; i++) {
00416         double average  = (prof->GetY()[i-1]+prof->GetY()[i+1])/2.0;
00417         if (prof->GetY()[i] < 0.1*average || 
00418             prof->GetY()[i] > 10.0*average) {
00419             prof->RemovePoint(i);
00420             i--;
00421         }
00422     }
00423 }

void BDProfMon::UpdateCache (  )  [private]

Definition at line 30 of file BDProfMon.cxx.

References fSpacing, RawBeamSwicData::GetData(), RawBeamData::GetName(), and Munits::mm.

Referenced by BDProfMon(), and SetData().

00031 {
00032     if (!this->IsValid()) return;
00033     RawBeamData& rbd = this->GetData();
00034     string name = rbd.GetName();
00035     if (name == "E:M121DS" || name == "E:MTGTDS")
00036         fSpacing = 0.5*Munits::mm;
00037     else
00038         fSpacing = 1.0*Munits::mm;
00039 }

double BDProfMon::WirePosition ( int  channel  )  const

Return the distance from the origin to the wire on the given profile monitor channel [1-44].

Origin is half way between channels 22 and 23. Units in Munits.

Definition at line 74 of file BDProfMon.cxx.

References fSpacing.

00075 {
00076     return (channel-22.5)*fSpacing;
00077 }

int BDProfMon::Xindex ( int  channel  )  [static]

Convert a horizontal channel [1-44] to an index in to the basic SWIC channels (0-95).

Returns:
an int which will be negative if an illegal channel is given, o.w. the associated swic channel is given.

Definition at line 64 of file BDProfMon.cxx.

Referenced by GetGaussFit(), and GetStats().

00065 {
00066     if (channel<1 || channel > 44) return -1;
00067     return channel+1;
00068 }

int BDProfMon::Yindex ( int  channel  )  [static]

Convert a vertical channel [1-44] to an index in to the basic SWIC channels (0-95).

Returns:
an int which will be negative if an illegal channel is given, o.w. the associated swic channel is given.

Definition at line 69 of file BDProfMon.cxx.

Referenced by GetGaussFit(), and GetStats().

00070 {
00071     if (channel<1 || channel > 44) return -1;
00072     return channel+49;
00073 }

Member Data Documentation

double BDProfMon::fSpacing [private]

Definition at line 112 of file BDProfMon.h.

Referenced by GetGaussFit(), UpdateCache(), and WirePosition().

The documentation for this class was generated from the following files:
Generated on Tue Jun 18 22:00:33 2013 for loon by  doxygen 1.4.7