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 Mon Aug 11 01:05:30 2014 for loon by  doxygen 1.4.7