FitterEM.cxx

Go to the documentation of this file.

#include <map>
#include <vector>
#include "TMatrixD.h"
#include "TF2.h"
#include "Conventions/PlaneView.h"
#include "Conventions/StripEnd.h"
#include "Validity/VldContext.h"
#include "Calibrator/Calibrator.h"
#include "MessageService/MsgService.h"
#include "FitterEM.h"
#include "KeyFunc.cc"
#include "rte.cc"
#include "CalcEM.cc"
#include "TROOT.h"
#include "TMath.h"

#define DELTA_T 0.145
#define DELTA_R 0.101
#define DELTA_PHI 2*2.*TMath::Pi()/360.
#define EPSILON 1
#define CUTOFF 1e-6 //1e-7
#define MIPCUT 1.0
#define NPAR 6 

using namespace std;

CVSID("$Id: FitterEM.cxx,v 1.2 2004/09/29 00:07:29 cbs Exp $");

FitterEM::FitterEM(){
  //this constructor is just for testing function
  //do not use
  fNStp = 1;

  fNPlanes = 0;
  fMinBoxPlane = 999;
  fMaxBoxPlane = 0;
  fNStrips = 0;
  fMinBoxStrip = 999;
  fMaxBoxStrip = 0;

  fclosestPlaneToVtx = -1;
  fclosestZedDiff = 99999;
  fclosestZedView = 0;
  fclosestUStripToVtx = -1;
  fclosestVStripToVtx = -1;
  fclosestTPosUDiff = 99999;
  fclosestTPosVDiff = 99999;

  fTotMip = 0;
  fPlaneArray = new Int_t[fNStp];
  fStripArray = new Int_t[fNStp];
  fViewArray  = new Int_t[fNStp];
  fMipArray   = new Double_t[fNStp];
  fZedArray   = new Double_t[fNStp];
  fTPosArray  = new Double_t[fNStp];
  fBoxPlaneArray = new Int_t[fNStp];
  fBoxStripArray = new Int_t[fNStp];
  fMips = NULL;

  fEnergy = 0;
  fVtx = NULL;
  fdudz = 1;
  fdvdz = 1;

  fCache = new CacherEM();
  fBinCorrel = NULL;

  fPass = true; //start out optimistic
  fFittedPars = NULL; //these always constitute the "Current" fit params
  fFittedErr = NULL;
  fChi2Min = 99999;
  fNDF = 0;
  fNiter = 0;
  fStepSize = NULL;  
  fNSteps = NULL;  

  fSignStat = 0;
  fRunStat = 0;
  fKSStat = 0;

}

//**************************************************
FitterEM::FitterEM(Int_t nstp,CandStripHandleItr itr){

  fNStp = nstp;
  fDataStpItr = itr;

  fNPlanes = 0;
  fMinBoxPlane = 999;
  fMaxBoxPlane = 0;
  fNStrips = 0;
  fMinBoxStrip = 999;
  fMaxBoxStrip = 0;

  fclosestPlaneToVtx = -1;
  fclosestZedDiff = 99999;
  fclosestZedView = 0;
  fclosestUStripToVtx = -1;
  fclosestVStripToVtx = -1;
  fclosestTPosUDiff = 99999;
  fclosestTPosVDiff = 99999;

  fTotMip = 0;
  fPlaneArray = new Int_t[fNStp];
  fStripArray = new Int_t[fNStp];
  fViewArray  = new Int_t[fNStp];
  fMipArray   = new Double_t[fNStp];
  fZedArray   = new Double_t[fNStp];
  fTPosArray  = new Double_t[fNStp];
  fBoxPlaneArray = new Int_t[fNStp];
  fBoxStripArray = new Int_t[fNStp];
  fMips = NULL;

  SetUpStripInfo();

  fEnergy = 0;
  fVtx = NULL;
  fdudz = 1;
  fdvdz = 1;

  fCache = new CacherEM();
  fBinCorrel = NULL;

  fPass = true; //start out optimistic
  fFittedPars = NULL; //these always constitute the "Current" fit params
  fFittedErr = NULL;
  fChi2Min = 99999;
  fNDF = 0;
  fNiter = 0;
  fStepSize = NULL;  
  fNSteps = NULL;  

  fSignStat = 0;
  fRunStat = 0;
  fKSStat = 0;

}


//**************************************************
FitterEM::~FitterEM(){
  delete [] fVtx;
  delete [] fFittedPars;
  delete [] fFittedErr;
  delete [] fPlaneArray;
  delete [] fStripArray;
  delete [] fViewArray;
  delete [] fMipArray;
  delete [] fZedArray;
  delete [] fTPosArray;
  delete [] fBoxPlaneArray;
  delete [] fBoxStripArray;
  delete [] fMips;
  delete [] fStepSize;
  delete [] fNSteps;
  delete fCache;
  if(fBinCorrel) delete fBinCorrel;
  TF2 *prof = NULL; //this should have been created in the course of fitting
  if(gROOT->FindObject("FITSHOWEREM_2DPROFILE")) { 
    prof = (TF2*) gROOT->FindObject("FITSHOWEREM_2DPROFILE");
    delete prof;   //delete it!
  }
  MSG("FitShowerEM",Msg::kDebug) << "Deleting FitterEM" << endl;
}

//**************************************************
void FitterEM::SetUpStripInfo(){

  Int_t icnt = 0;
  Calibrator& cal=Calibrator::Instance();
  
  while (CandStripHandle *stp = fDataStpItr()) {
    
    if(icnt==0) {
      const VldContext *vc = stp->GetVldContext();
      cal.Reset(*vc);
    }

    int plane = stp->GetPlane();
    int strip = stp->GetStrip();
    fPlaneArray[icnt] = plane;
    fStripArray[icnt] = strip;
    fViewArray[icnt]  = stp->GetPlaneView(); //U == 2, V == 3;
    fMipArray[icnt]   = cal.GetMIP(stp->GetCharge(CalDigitType::kSigCorr,
                                                  StripEnd::kWhole));
    fTotMip += fMipArray[icnt];
    fZedArray[icnt]   = stp->GetZPos();
    fTPosArray[icnt]  = stp->GetTPos();
    fBoxPlaneArray[icnt] = 0;
    fBoxStripArray[icnt] = 0;

    MSG("FitShowerEM", Msg::kVerbose) << "pl="<<fPlaneArray[icnt]
                                      <<" st="<<fStripArray[icnt]
                                      <<" pv="<<fViewArray[icnt]
                                      <<" mips="<<fMipArray[icnt]
                                      <<endl;
    
    icnt++;
  }

}

//**************************************************
void FitterEM::SetInputParams(Double_t energy,Double_t* vtx,
                              Double_t dudz,Double_t dvdz) {
  fEnergy = energy;
  fVtx = new Double_t[3]; 
  fVtx[0] = vtx[0]; 
  fVtx[1] = vtx[1]; 
  fVtx[2] = vtx[2];
  fdudz = dudz; 
  fdvdz = dvdz;
  SetUpFitBox();
}

//**************************************************
void FitterEM::QuickInput(Double_t energy,Double_t vtxu,Double_t vtxv,
                          Double_t vtxz,Double_t dudz,Double_t dvdz,
                          Int_t view) {
  
  fFittedPars = new Double_t[NPAR];
  fFittedPars[0] = energy;
  fFittedPars[1] = vtxu/0.041; //vertex
  fFittedPars[2] = vtxv/0.041;
  fFittedPars[3] = vtxz/0.060;
  fFittedPars[4] = TMath::ATan(dudz*6.0/4.1); //angles in rads 
  fFittedPars[5] = TMath::ATan(dvdz*6.0/4.1); //in pln-stp coords
  
  fclosestZedView = view;

  MSG("FitShowerEM", Msg::kDebug) 
    << "Fit Parameters to start with:\n" 
    << "Energy = "<<fFittedPars[0]<<"\n"
    << "UVtx   = "<<fFittedPars[1]<<"\n"
    << "VVtx   = "<<fFittedPars[2]<<"\n"
    << "ZVtx   = "<<fFittedPars[3]<<"\n"
    << "UAng   = "<<fFittedPars[4]<<"\n"
    << "VAng   = "<<fFittedPars[5]<<endl;

}

//**************************************************
Double_t *FitterEM::GetFittedPars(){
  //convert back from internal coords to CandFitShowerEM coords
  fFittedPars[0] = fFittedPars[0]; //energy is ok
  fFittedPars[1] = fFittedPars[1]*0.041+fVtx[0]-fclosestTPosUDiff;//u vertex
  fFittedPars[2] = fFittedPars[2]*0.041+fVtx[1]-fclosestTPosVDiff;//v vertex
  fFittedPars[3] = fFittedPars[3]*0.060+fVtx[2]-fclosestZedDiff;//z vertex
  fFittedPars[4] = TMath::Tan(fFittedPars[4])*4.1/6.0; //dudz
  fFittedPars[5] = TMath::Tan(fFittedPars[5])*4.1/6.0;//dvdz
  return fFittedPars; 
}

//**************************************************
void FitterEM::SetStepSize(Double_t *stepSize) {
  
  fStepSize = new Double_t[NPAR];
  for(int i=0;i<NPAR;i++) fStepSize[i] = stepSize[i];
  
  fStepSize[0] *= fEnergy;           //assume a %    --> GeV
  fStepSize[1] /= 0.041;             //assume in m   --> units of strip
  fStepSize[2] /= 0.041;             //assume in m   --> units of strip
  fStepSize[3] /= 0.060;             //assume in m   --> units of plane
  
  //assume in deg in metre-metre --> rads in plane-strip
  fStepSize[4]=TMath::ATan(TMath::Tan(fStepSize[4]*TMath::Pi()/180.)*6.0/4.1);
  fStepSize[5]=TMath::ATan(TMath::Tan(fStepSize[5]*TMath::Pi()/180.)*6.0/4.1);
  
  MSG("FitShowerEM",Msg::kDebug) << "Step Sizes: " << fStepSize[0] 
                                 << " " << fStepSize[1] << " " << fStepSize[2]
                                 << " " << fStepSize[3] << " " << fStepSize[4]
                                 << " " << fStepSize[5] << endl;
  
}

//**************************************************
void FitterEM::SetNSteps(Int_t *NSteps) {
  
  fNSteps = new Int_t[NPAR];
  for(int i=0;i<NPAR;i++) fNSteps[i] = NSteps[i];

  MSG("FitShowerEM",Msg::kDebug) << "Number of Steps: " << fNSteps[0] 
                                 << " " << fNSteps[1] << " " << fNSteps[2]
                                 << " " << fNSteps[3] << " " << fNSteps[4]
                                 << " " << fNSteps[5] << endl;
}

//**************************************************
void FitterEM::SetUpFitBox() {

  //This is the tricky part...
  //need to move into (plane,strip) Box with vertex approx at (0,0,0)
  //adjust all strip numbers by appropriate transverse vtx
   //adjust all plane numbers by z vtx
  //adjust angles to be u-z, v-z projection angles in rads 
  //calculated in units of plane and strip
  //need to set fFittedPars to be the new angles and vtx in the Box coords
  
  for(int i=0;i<fNStp;i++){
    float ZedDiff = fZedArray[i]-fVtx[2];
    if(fabs(ZedDiff)<fabs(fclosestZedDiff)){
      fclosestZedDiff = ZedDiff;
      fclosestPlaneToVtx = fPlaneArray[i];
      fclosestZedView = fViewArray[i];
    }
  }
  
  Int_t closestOtherViewPlaneToVtxPlane = 99999;
  for(int i=0;i<fNStp;i++){
    //if this is the vertex plane, use it for finding transverse vertex
    if(fPlaneArray[i]==fclosestPlaneToVtx){
      if(fViewArray[i]==2){ //U strips
        float TPosDiff = fTPosArray[i]-fVtx[0];
        if(fabs(TPosDiff)<fabs(fclosestTPosUDiff)){
          fclosestTPosUDiff = TPosDiff;
          fclosestUStripToVtx = fStripArray[i];
        }
      }
      else if(fViewArray[i]==3){ //V strips
        float TPosDiff = fTPosArray[i]-fVtx[1];
        if(fabs(TPosDiff)<fabs(fclosestTPosVDiff)){
          fclosestTPosVDiff = TPosDiff;
          fclosestVStripToVtx = fStripArray[i];
        }
      }
    }
    //otherwise, if plane-view is not same as vertex plane-view, 
    //then use it to find vertex in other view:
    else if(fViewArray[i]!=fclosestZedView){
      //look for closest plane in other view:
      if(abs(fPlaneArray[i]-fclosestPlaneToVtx)<= 
         abs(closestOtherViewPlaneToVtxPlane-fclosestPlaneToVtx)){
        closestOtherViewPlaneToVtxPlane=fPlaneArray[i];
        if(fViewArray[i]==2){ //U strips
          float TPosDiff = fTPosArray[i]-fVtx[0];
          if(fabs(TPosDiff)<fabs(fclosestTPosUDiff)){
            fclosestTPosUDiff = TPosDiff;
            fclosestUStripToVtx = fStripArray[i];
          }
        }
        else if(fViewArray[i]==3){ //V strips
          float TPosDiff = fTPosArray[i]-fVtx[1];
          if(fabs(TPosDiff)<fabs(fclosestTPosVDiff)){
            fclosestTPosVDiff = TPosDiff;
            fclosestVStripToVtx = fStripArray[i];
          }
        }
      }
    }
  }
  

  MSG("FitShowerEM", Msg::kDebug) 
    << "Box Quantities: Closest Plane To Vtx="<<fclosestPlaneToVtx
    << " with ZedDiff="<<fclosestZedDiff
    << " and in View="<<fclosestZedView<<"\n"
    << "Closest U Strip to Vtx="<<fclosestUStripToVtx
    << " with TPosDiff="<<fclosestTPosUDiff<<"\n"
    << "Closest V Strip to Vtx="<<fclosestVStripToVtx
    << " with TPosDiff="<<fclosestTPosVDiff<<"\n"
    << "(Vertex from CandShowerEM=["<<fVtx[0]<<","<<fVtx[1]<<","<<fVtx[2]<<"])"
    << endl;
  
  for(int i=0;i<fNStp;i++){
    
    fBoxPlaneArray[i] = fPlaneArray[i]-fclosestPlaneToVtx;
    if(fViewArray[i]==2) {//U strips  
      fBoxStripArray[i] = fStripArray[i]-fclosestUStripToVtx; 
    }
    if(fViewArray[i]==3) {//V strips
      fBoxStripArray[i] = fStripArray[i]-fclosestVStripToVtx; 
    }
    
    if(fBoxPlaneArray[i]<fMinBoxPlane) fMinBoxPlane = fBoxPlaneArray[i];
    if(fBoxPlaneArray[i]>fMaxBoxPlane) fMaxBoxPlane = fBoxPlaneArray[i];
    if(fBoxStripArray[i]<fMinBoxStrip) fMinBoxStrip = fBoxStripArray[i];
    if(fBoxStripArray[i]>fMaxBoxStrip) fMaxBoxStrip = fBoxStripArray[i];
  }


  fNPlanes = fMaxBoxPlane-fMinBoxPlane+1;
  fNStrips = fMaxBoxStrip-fMinBoxStrip+1;
  
  //now make fMips[plane*strip] array to simplify Mip value access in fit
  Int_t TotStp = fNPlanes*fNStrips;
  fMips = new Double_t[TotStp];
  for(int i=0;i<TotStp;i++) fMips[i] = 0;
  for(int i=0;i<fNStp;i++){    
    int n = ((fBoxPlaneArray[i]-fMinBoxPlane)*fNStrips 
             + (fBoxStripArray[i] - fMinBoxStrip));
    fMips[n] = 100.*fMipArray[i]/fTotMip; //make everything a percentage
  }
  
  MSG("FitShowerEM", Msg::kDebug) 
    << "Box Limits: "<<"\n"
    << "NPlanes="<<fNPlanes<<" with Min="<<fMinBoxPlane
    << " and Max="<<fMaxBoxPlane<<"\n"
    << "NStrips="<<fNStrips<<" with Min="<<fMinBoxStrip
    << " and Max="<<fMaxBoxStrip<<endl;

  fFittedPars = new Double_t[NPAR];
  fFittedPars[0] = fEnergy;
  fFittedPars[1] = fclosestTPosUDiff/0.041; //vertex
  fFittedPars[2] = fclosestTPosVDiff/0.041;
  fFittedPars[3] = fclosestZedDiff/0.060;  
  fFittedPars[4] = TMath::ATan(fdudz*6.0/4.1); //angles in rads 
  fFittedPars[5] = TMath::ATan(fdvdz*6.0/4.1); //in pln-stp coords
  
  MSG("FitShowerEM", Msg::kDebug) 
    << "Fit Parameters to start with:\n" 
    << "Energy = "<<fFittedPars[0]<<"\n"
    << "UVtx   = "<<fFittedPars[1]<<"\n"
    << "VVtx   = "<<fFittedPars[2]<<"\n"
    << "ZVtx   = "<<fFittedPars[3]<<"\n"
    << "UAng   = "<<fFittedPars[4]<<"\n"
    << "VAng   = "<<fFittedPars[5]<<endl;

  fBinCorrel = new BinCorrelationEM(fFittedPars);

}
  

//**************************************************
Double_t FitterEM::PredictEMLoss(Int_t plane,Int_t strip,Double_t &flucVal)
{

  if(!fCache->ValidPS(fFittedPars)){
    fCache->ClearPSCache();
    double sumE = 0;
    vector<rte> rtmap;
    
    if(!fCache->ValidRT(fFittedPars)){
      fCache->ClearCache();
      rtmap = RTCalc(fFittedPars[0],DELTA_T,DELTA_R,EPSILON,CUTOFF,sumE);
      fCache->StoreRT(NPAR,fFittedPars,rtmap,sumE);
      fCache->PrintRTMap();
    }
    else {
      rtmap = fCache->GetRTMap();
      sumE = fCache->GetSumE();
    }
    
    map<int,spef> psmap = PSCalc(rtmap,fclosestZedView,
                                 fFittedPars[1],fFittedPars[2],
                                 fFittedPars[3],fFittedPars[4],
                                 fFittedPars[5],DELTA_PHI,sumE);
    fCache->StorePS(NPAR,fFittedPars,psmap);
    fCache->PrintPSMap();
  }
 
  map<int,spef> spmap = fCache->GetPSMap();
  int key = MakeKey(plane,strip);
  map<int,spef>::iterator theOne = spmap.find(key);
  map<int,spef>::iterator end = spmap.end();
  if(theOne==end) {
    flucVal = 0.;
    return 0;
  }
  flucVal = theOne->second.fl;
  return theOne->second.en;
  
}

//**************************************************
Double_t FitterEM::CalculateChi2(){
  
  fBinCorrel->ReInit(fFittedPars);
  Double_t *ErrArray = fBinCorrel->GetErrArray(fMinBoxPlane,fMaxBoxPlane,
                                               fMinBoxStrip,fMaxBoxStrip);
  
  Int_t TotStp = fNPlanes*fNStrips;

  //for holding all necessary function values during calculation
  Double_t *ff = new Double_t[TotStp];  //fit func values
  Double_t *fff = new Double_t[TotStp]; //fit func fluctuations
  for(Int_t i=0;i<TotStp;i++) { 
    ff[i] = -1; 
    fff[i] = 0.; 
  }
  
  Double_t lastDiagTerm = 1;
  Double_t flucVal = 0;
  //normalise covariance matrix using function: (make ff a percentage)
  for(Int_t i=0;i<TotStp;i++){
    if(ff[i]==-1) {
      ff[i] = 100.*PredictEMLoss(fMinBoxPlane+i/fNStrips,
                                 fMinBoxStrip+i%fNStrips,flucVal);    
      fff[i] = 100.*flucVal;
    }
    for(Int_t j=0;j<TotStp;j++){
      if(ff[j]==-1) {
        ff[j] = 100.*PredictEMLoss(fMinBoxPlane+j/fNStrips,
                                   fMinBoxStrip+j%fNStrips,flucVal);
        fff[j] = 100.*flucVal;
      }
      ErrArray[i*TotStp+j]*=(fff[i]*fff[j]);
      if(i==j) {
        if(ErrArray[i*TotStp+j]==0) ErrArray[i*TotStp+j] = lastDiagTerm;
        else lastDiagTerm = ErrArray[i*TotStp+j];
        MSG("FitShowerEM", Msg::kVerbose) << "ErrArray: " 
                                          << ErrArray[i*TotStp+j]
                                          << endl;      
      }
      
    }
  }
  
  //invert covariance matrix:
  Double_t *InvArray = new Double_t[TotStp*TotStp];
  TMatrixD matrix(TotStp,TotStp);
  matrix.SetMatrixArray(ErrArray);
  matrix.Invert();
  if(matrix.IsValid()) matrix.GetMatrix2Array(InvArray);
  else return -1.0;

  //loop through again to calculate chi2
  Double_t Chi2 = 0;
  Double_t signStat = 0;
  Double_t nValidPoints = 0;
  for(Int_t i=0;i<TotStp;i++){
    MSG("FitShowerEM", Msg::kVerbose) << "MIP: " << fMips[i] 
                                      << " FF: " << ff[i] << endl;
    if(fMips[i]>ff[i]) signStat+=1;
    nValidPoints+=1;
    for(Int_t j=0;j<TotStp;j++){
      Chi2+=(fMips[i]-ff[i])*InvArray[i*TotStp+j]*(fMips[j]-ff[j]);
    }
  }
  
  fSignStat = signStat/nValidPoints;
  fNDF=TotStp;
 
  delete [] InvArray;
  delete [] ff;

  return Chi2;
  
}

//**************************************************
void FitterEM::DoFit(){

  MSG("FitShowerEM", Msg::kDebug) << "In DoFit()" << endl;

  fFittedErr = new Double_t[NPAR];
  
  Double_t *BestPars = new Double_t[NPAR];  
  Double_t *parlowlim = new Double_t[NPAR];
  Double_t *parupplim = new Double_t[NPAR];
  Int_t nvarypar = 0;
  
  for(Int_t i=0;i<NPAR;i++){    
    if(fStepSize[i]==0||fNSteps[i]<1) {
      parlowlim[i] = fFittedPars[i];
      parupplim[i] = fFittedPars[i];
    }
    else {
      parlowlim[i] = fFittedPars[i] - fNSteps[i]*fStepSize[i];
      parupplim[i] = fFittedPars[i] + fNSteps[i]*fStepSize[i];      
      nvarypar += 1;
    }
  }

  for(Double_t a=parlowlim[0];a<=parupplim[0];a+=fStepSize[0]){
    fFittedPars[0] = a;
    for(Double_t b=parlowlim[1];b<=parupplim[1];b+=fStepSize[1]){
      fFittedPars[1] = b;
      for(Double_t c=parlowlim[2];c<=parupplim[2];c+=fStepSize[2]){
        fFittedPars[2] = c;
        for(Double_t d=parlowlim[3];d<=parupplim[3];d+=fStepSize[3]){
          fFittedPars[3] = d;
          for(Double_t e=parlowlim[4];e<=parupplim[4];e+=fStepSize[4]){
            fFittedPars[4] = e;
            for(Double_t f=parlowlim[5];f<=parupplim[5];f+=fStepSize[5]){
              fFittedPars[5] = f;

              MSG("FitShowerEM", Msg::kDebug) 
                <<a<<" "<<b<<" "<<c<<" "<<d<<" "<<e<<" "<<f<<endl;
              
              if(fFittedPars[0]<0.05) MSG("FitShowerEM", Msg::kWarning) 
                <<"Energy < 50MeV! Not calculating Chi2" << endl;
              Double_t Chi2 = CalculateChi2();
              fNiter+=1;
              if(Chi2==-1) continue; //problem with calculation
              if(Chi2<fChi2Min) {
                fChi2Min = Chi2;
                for(Int_t i=0;i<NPAR;i++) BestPars[i] = fFittedPars[i];
                fBestSignStat = fSignStat;
                //DoTests();
              }
              MSG("FitShowerEM", Msg::kDebug) << "Chi2: "
                                              << Chi2 << endl;
            }
          }
        }
      }
    }
  }
  
  for(Int_t i=0;i<NPAR;i++) fFittedPars[i] = BestPars[i];
  fNDF -= nvarypar;

  fPass = true;
  if(fChi2Min==99999||fNDF<=0) fPass = false;

}

//**************************************************
void FitterEM::DoTests(){

  Int_t nValidData        = 0;
  Int_t nValidFunc        = 0;
  Double_t *DataVals      = new Double_t[fNPlanes*fNStrips];
  Double_t *DataVals_Asc  = new Double_t[fNPlanes*fNStrips];
  Double_t *FuncVals      = new Double_t[fNPlanes*fNStrips];
  Double_t *FuncVals_Asc  = new Double_t[fNPlanes*fNStrips];
  Double_t lowest_DataVal = 99999;
  Double_t lowest_FuncVal = 99999;
  Double_t flucVal = 0;

  for(Int_t i=0;i<fNPlanes;i++){
    for(Int_t j=1;j<=fNStrips;j++){
      Float_t bc = fMips[i*fNStrips+j];
      Float_t ff = PredictEMLoss(fMinBoxPlane+i,fMinBoxStrip+j,flucVal);

      if(ff>0){
        FuncVals[nValidFunc] = ff;
        nValidFunc+=1;
        if(ff<lowest_FuncVal) lowest_FuncVal = ff;
      }
      
      if(bc>0){
        DataVals[nValidData] = bc;
        nValidData+=1;
        if(bc<lowest_DataVal) lowest_DataVal = bc;
      }
    }
  }
  
  //put arrays into ascending order:
  DataVals_Asc[0] = lowest_DataVal;
  for(Int_t i=1;i<nValidData;i++){
    Double_t lowest_DataVal = 999999;
    for(Int_t j=0;j<nValidData;j++){
      if(DataVals[j]<lowest_DataVal&&DataVals[j]>DataVals_Asc[i-1])
        lowest_DataVal = DataVals[j];
    }
    DataVals_Asc[i] = lowest_DataVal;
  }
  delete [] DataVals;

  FuncVals_Asc[0] = lowest_FuncVal;
  for(Int_t i=1;i<nValidFunc;i++){
    Double_t lowest_FuncVal = 999999;
    for(Int_t j=0;j<nValidFunc;j++){
      if(FuncVals[j]<lowest_FuncVal&&FuncVals[j]>FuncVals_Asc[i-1])
        lowest_FuncVal = FuncVals[j];
    }
    FuncVals_Asc[i] = lowest_FuncVal;
  }
  delete [] FuncVals;

  //calculate number of runs and mean and variance of expected number of runs:
  Double_t mean_runs = 1. + 2.*nValidFunc*nValidData/(nValidData+nValidFunc);
  Double_t var_runs = (2.*nValidFunc*nValidData
                     *(2.*nValidFunc*nValidData - nValidData -nValidFunc))
    /((nValidFunc+nValidData)*(nValidFunc+nValidData)
      *(nValidFunc+nValidData-1));

  Double_t r = 0;
  Int_t nD = 0;
  Int_t nF = 0;
  Int_t currentBun = 0; // 1 => Data was last lowest; -1 => Func was last lowest
  Bool_t change = false;

  for(Int_t i=0;i<nValidData+nValidFunc;i++){
    if(DataVals_Asc[nD]<FuncVals_Asc[nF]) {
      nD+=1;
      if(currentBun!=1) {
        change = true;
        currentBun = 1;
      }
    }
    else {
      nF+=1;
      if(currentBun!=-1) {
        change = true;
        currentBun = -1;
      } 
    }
    if(change) {
      r+=1;
      change = false;
    }
    if(nD>=nValidData) break;
    if(nF>=nValidFunc) break;
  }

  fRunStat = (r - mean_runs)/sqrt(var_runs);
  
  //Now do K-S test:
  for(Int_t i=0;i<nValidData-1;i++){
    //between DataVals_Asc[i] and DataVals_Asc[i+1]
    //cumulative distribution value is (i+1)/nValidData
    //Estimate from Func by interpolating at 
    //(DataVals_Asc[i+1]-DataVals_Asc[i])/2.
    Double_t S =  Double_t(i+1)/Double_t(nValidData);
    Double_t theVal = (DataVals_Asc[i+1]+DataVals_Asc[i])/2.;
    Int_t funcIndex = -1;
    for(Int_t j=1;j<nValidFunc;j++){
      if(FuncVals_Asc[j]>theVal) {
        funcIndex = j;
        break;
      }
    }
    if(funcIndex==-1) funcIndex = nValidFunc-1;
    Double_t estimator = funcIndex-1. + (FuncVals_Asc[funcIndex]-theVal)
      /(FuncVals_Asc[funcIndex]-FuncVals_Asc[funcIndex-1]);
    estimator/=Double_t(nValidFunc);
    
    if(fabs(S-estimator)>fKSStat) fKSStat = fabs(S-estimator);
    
  }

  delete [] DataVals_Asc;
  delete [] FuncVals_Asc;
  return;

}

---