#include <BDProfMon.h>
Inheritance diagram for BDProfMon:

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 |
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.
Created on: Fri Apr 15 13:09:38 2005
Definition at line 40 of file BDProfMon.h.
| BDProfMon::BDProfMon | ( | ) |
| 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] |
| 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
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.
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] |
| 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).
Definition at line 64 of file BDProfMon.cxx.
Referenced by GetGaussFit(), and GetStats().
| int BDProfMon::Yindex | ( | int | channel | ) | [static] |
Convert a vertical channel [1-44] to an index in to the basic SWIC channels (0-95).
Definition at line 69 of file BDProfMon.cxx.
Referenced by GetGaussFit(), and GetStats().
double BDProfMon::fSpacing [private] |
Definition at line 112 of file BDProfMon.h.
Referenced by GetGaussFit(), UpdateCache(), and WirePosition().
1.4.7