TestAcc.cc

Go to the documentation of this file.

#include <iostream> //System includes
#include <fstream>
#include <vector>
#include <string.h>
#include <stdio.h>

#include "TVector3.h" //ROOT includes
#include "TNtuple.h"
#include "TFile.h"

#include "MessageService/MsgService.h"
#include "BField/BfldValidate.h" //BField headers
#include "BField/BfldLoanPool.h"
#include "BField/BField.h"

#include "TestConstants.h" //My own headers

CVSID("$Id: TestAcc.cc,v 1.2 2001/08/15 18:57:27 agoldst Exp $");

const int kDefaultNumAccTestPts = 3012;

const char kRectDataName[] = "rectData";
const char kVorDataName[] = "vorData";

void BenchmarkAccuracy(const char *outFileName,const char *rectMapFileName,
                       const char *vorMapFileName,const char *vorMeshFileName,
                       const char *vorBFldFileName,
                       const char *rectAccTestPtsFileName,
                       const char *vorAccTestPtsFileName);
bool GetAccTestPts(vector<TVector3> &pts,vector<TVector3> &bfld,
                   const char *fileName);
void AccTest(const vector<TVector3> &pts,const vector<TVector3> &trueB,
             BField *field,TNtuple *outData,int algorithm);

int GetNumPointsRect(const char *fileName);
int GetNumPointsVor(const char *fileName);
BField *NewBFieldRect(const char *mapFileName);
BField *NewBFieldVor(const char *vorMapFileName,const char *meshFileName,
                     const char *vorBFieldFileName);

//define WRITE_OUT_BIG_ERR_PTS

#ifdef WRITE_OUT_BIG_ERR_PTS
const double kSmallestBigErr = 0.10; //10% is big!

ofstream *gBigErrPtsFile;
#endif

void BenchmarkAccuracy(const char *outFileName,const char *rectMapFileName,
                       const char *vorMapFileName,const char *vorMeshFileName,
                       const char *vorBFldFileName,
                       const char *rectAccTestPtsFileName,
                       const char *vorAccTestPtsFileName) {

  MSG("Bfldtest",Msg::kDebug) << "Running accuracy test..." << endl;

  TFile outFile(outFileName,"RECREATE"); //Root datafile to store results.

  vector<TVector3> testPts(kDefaultNumAccTestPts);
  vector<TVector3> testField(kDefaultNumAccTestPts);
  
  BField *theRectField, *theVorField;
  TNtuple *theData; //Ntuple to store data as we go

#ifdef WRITE_OUT_BIG_ERR_PTS
  gBigErrPtsFile = new ofstream("bigErrPts.bfb",ios::out);
#endif

  MSG("Bfldtest",Msg::kDebug) << "Building rect2d BField using file "
                              << rectMapFileName << endl;

  theRectField = NewBFieldRect(rectMapFileName);
  
  if(!GetAccTestPts(testPts,testField,rectAccTestPtsFileName)) {
    MSG("Bfldtest",Msg::kWarning)
      << "Problem reading in rect2d test point file. Bailing." << endl;
    return;
  }

  theData = new TNtuple(kRectDataName,kRectDataName,
                        "x:y:z:btx:bty:btz:brcx:brcy:brcz:brbx:brby:brbz:brpx:brpy:brpz:brvx:brvy:brvz");

  //Translation:
  //First letter is 'b' for magnetic field.
  //Second letter is algorithm, 'r' for rect and 'v' for vor
  //Third letter is interpolation method's first letter exc.
  //PlanarVec, which is 'v'.
  //Fourth letter is component of the vector.

  AccTest(testPts,testField,theRectField,theData,kAlgRect);

  MSG("Bfldtest",Msg::kDebug) << "Saving Ntuple " << kRectDataName << endl;

  //Because TFile::Write will write whatever ntuples are in memory to disk,
  //we have to remove the rect2d data ntuple from memory before we build
  //and write the Voronoi data ntuple.

  outFile.Write();
  delete theData;
  delete theRectField;

  MSG("Bfldtest",Msg::kDebug)
    << "Building Voronoi BField using the following: " << endl;
  MSG("Bfldtest",Msg::kDebug) << "Voronoi diagram: " << vorMapFileName << endl;
  MSG("Bfldtest",Msg::kDebug) << "Voronoi mesh: " << vorMeshFileName << endl;
  MSG("Bfldtest",Msg::kDebug) << "Generator BField data: "
                              << vorBFldFileName << endl;
  
  theVorField = NewBFieldVor(vorMapFileName,vorMeshFileName,
                             vorBFldFileName);
  
  if(!GetAccTestPts(testPts,testField,vorAccTestPtsFileName)) {
    MSG("Bfldtest",Msg::kWarning)
      << "Problem reading in Voronoi test point file. Bailing." << endl;
    return;
  }

  theData = new TNtuple(kVorDataName,kVorDataName,
                        "x:y:z:btx:bty:btz:bvcx:bvcy:bvcz:bvpx:bvpy:bvpz");

  AccTest(testPts,testField,theVorField,theData,kAlgVor);
  
  MSG("Bfldtest",Msg::kDebug) << "Saving Ntuple " << kVorDataName << endl;

  outFile.Write();
  delete theData;
  delete theVorField;

  MSG("Bfldtest",Msg::kInfo) << "Done with accuracy test." << endl;
  MSG("Bfldtest",Msg::kInfo) << "Output written to " << outFileName << endl;
}

void AccTest(const vector<TVector3> &pts,const vector<TVector3> &trueB,
             BField *field,TNtuple *outData,int algorithm) {
  if(algorithm == kAlgRect) {
    MSG("Bfldtest",Msg::kDebug)
      << "Running accuracy test for rect2d algorithm..." << endl;
  }
  else {
    MSG("Bfldtest",Msg::kDebug)
      << "Running accuracy test for Voronoi algorithm..." << endl;
  }

  int i,j;
  TVector3 r,bt,brc,brb,brp,brv,bvc,bvp;
  float *row;

  //***Is BField smart enough to handle a change in interpolation
  //methods in mid-run?
  //I suspect not...in principle, though, it ought to work.

  //Write out the data into the ntuple, in the ugliest possible way.

  for(i = 0;i < pts.size();i++) {
    r = pts[i];
    bt = trueB[i];

    j = 0;
    if(algorithm == kAlgRect) {
      row = new float[18];
    }
    else {
      row = new float[12];
    }

    row[j] = r.X(); ++j;
    row[j] = r.Y(); ++j;
    row[j] = r.Z(); ++j;
    
    row[j] = bt.X(); ++j;
    row[j] = bt.Y(); ++j;
    row[j] = bt.Z(); ++j;


    if(algorithm == kAlgRect) {
      field->SetInterpMethod(BfldInterpMethod::kClosest);
      brc = field->GetBField(r);
      field->SetInterpMethod(BfldInterpMethod::kBilinear);
      brb = field->GetBField(r);
      field->SetInterpMethod(BfldInterpMethod::kPlanar);
      brp = field->GetBField(r);
      field->SetInterpMethod(BfldInterpMethod::kPlanarVec);
      brv = field->GetBField(r);
   
      row[j] = brc.X(); ++j;
      row[j] = brc.Y(); ++j;
      row[j] = brc.Z(); ++j;
      
      row[j] = brb.X(); ++j;
      row[j] = brb.Y(); ++j;
      row[j] = brb.Z(); ++j;
      
      row[j] = brp.X(); ++j;
      row[j] = brp.Y(); ++j;
      row[j] = brp.Z(); ++j;
      
      row[j] = brv.X(); ++j;
      row[j] = brv.Y(); ++j;
      row[j] = brv.Z();
    

#ifdef WRITE_OUT_BIG_ERR_PTS
      double err = (brb-bt).Mag() / bt.Mag();
      if(err > kSmallestBigErr) {
        *gBigErrPtsFile << "r " << r.X() << " " << r.Y()
                        << " " << err << endl;
      }
#endif

    }
    else { //so algorithm == kAlgVor
      field->SetInterpMethod(BfldInterpMethod::kClosest);
      bvc = field->GetBField(r);
      field->SetInterpMethod(BfldInterpMethod::kPlanar);
      bvp = field->GetBField(r);
    
      row[j] = bvc.X(); ++j;
      row[j] = bvc.Y(); ++j;
      row[j] = bvc.Z(); ++j;

      row[j] = bvp.X(); ++j;
      row[j] = bvp.Y(); ++j;
      row[j] = bvp.Z();

#ifdef WRITE_OUT_BIG_ERR_PTS
      double err = (bvp-bt).Mag() / bt.Mag();
      if(err > kSmallestBigErr) {
        *gBigErrPtsFile << "v " << r.X() << " " << r.Y()
                        << " " << err << endl;
      }
#endif
    } 

    outData->Fill(row);
    delete row;
  }
}

bool GetAccTestPts(vector<TVector3> &pts,vector<TVector3> &bfld,
                   const char *fileName) {
  //Read the data from fileName into pts and bfld
  //For now, assume the data will have the form of an ANSYS output file,
  //to wit, a series of lines of the form
  //id  x  y  z bx  by  bz
  //where all variables (including id) are floats

  MSG("Bfldtest",Msg::kDebug)
    << "Reading in test points file " << fileName << endl;

  ifstream inFile(fileName,ios::in);
  if(!inFile.good()) {
    MSG("Bfldtest",Msg::kWarning)
      << "Error opening test points file " << fileName << endl;
    return false;
  }

  if(pts.size() > bfld.size()) //Make sure the arrays match in size.
    bfld.resize(pts.size());
  else
    pts.resize(bfld.size());
  int n = pts.size();

  float id_ignored,x,y,z,bx,by,bz;
  int i = 0;

  while(!inFile.eof()) {
    inFile >> id_ignored >> x >> y >> z >> bx >> by >> bz;
    if(i >= n) { //The poor man's dynamically resizing array.
      n *= 2;
      pts.resize(n);
      bfld.resize(n);
    }

    pts[i].SetXYZ(x,y,z);
    bfld[i].SetXYZ(bx,by,bz);
    i++;
  }

  pts.resize(i); //Truncate the arrays to exactly the number
  bfld.resize(i); //of points read in.

  MSG("Bfldtest",Msg::kDebug)
    << "Read in " << i << " test points and B values." << endl;
  return true;
}

---