loon: BDProfMon.cxx Source File

00001 #include "BDProfMon.h"
00002 
00003 #include <Conventions/Munits.h>
00004 
00005 #include <cmath>
00006 #include <string>
00007 #include <cassert>
00008 #include <TMinuit.h>
00009 #include <TSystem.h>
00010 #include <TF1.h>
00011 #include <MessageService/MsgService.h>
00012 CVSID("$Id: BDProfMon.cxx,v 1.10 2008/11/07 16:46:04 bv Exp $");
00013 
00014 using namespace std;
00015 
00016 
00017 BDProfMon::BDProfMon() 
00018     : BDSwicDevice(), fSpacing(0)
00019 {
00020 }
00021 BDProfMon::BDProfMon(const RawBeamData& swic_data)
00022     : BDSwicDevice(swic_data), fSpacing(0)
00023 {
00024     this->UpdateCache();
00025 }
00026 BDProfMon::~BDProfMon()
00027 {
00028 }
00029 
00030 void BDProfMon::UpdateCache()
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 }
00040 
00041 void BDProfMon::SetData(const RawBeamData& swic_data)
00042 {
00043     this->BDSwicDevice::SetData(swic_data);
00044     this->UpdateCache();
00045 }
00046 /* Channels are aranged like:
00047 
00048 entire buffer:     0   1   2 ...  45  46  47 |  48  49  50 ...  93  94  95
00049 wire data block: 104 105 106 ... 149 150 151 | 152 153 154 ... 197 198 199
00050 channel (fortran): -   -   1 ...  44   -   - |   -   -   1 ...  44   -   -
00051 channel (C++):     -   -   0 ...  43   -   - |   -   -   0 ...  43   -   -
00052 view:                         X              |              Y
00053 
00054 */
00055 
00056 bool BDProfMon::SetMask(const std::vector<double>& cmask)
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 }
00064 int BDProfMon::Xindex(int channel)
00065 {
00066     if (channel<1 || channel > 44) return -1;
00067     return channel+1;
00068 }
00069 int BDProfMon::Yindex(int channel)
00070 {
00071     if (channel<1 || channel > 44) return -1;
00072     return channel+49;
00073 }
00074 double BDProfMon::WirePosition(int channel) const
00075 {
00076     return (channel-22.5)*fSpacing;
00077 }
00078 double BDProfMon::GetStats(double& xmean, double &ymean, double &xrms, double &yrms) const
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 }
00126 
00127 
00128 double BDProfMon::GetStats(TGraphErrors *prof, double& mean, double& rms) const {
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 }
00151                            
00152 double BDProfMon::GetGaussFit(double& xmean, double &ymean, double &xsigma, double &ysigma) const
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 }
00401 
00402 void BDProfMon::SuppressDeadHot(TGraphErrors *prof) const
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 }
00424 
00425 
00426 
00427 int BDProfMon::Index(int /*channel*/) { const int do_not_call=0; assert(do_not_call); }