BfldDebugData.cc

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <vector>
#include <unistd.h>
#include <string.h>
#include <stdlib.h>

#include "TROOT.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TVector3.h"

#include "MessageService/MsgService.h"

#include "BField/BfldValidate.h"
#include "BField/BfldLoanPool.h"
#include "BField/BFLHandler.h"

CVSID("$Id: BfldDebugData.cc,v 1.4 2007/11/11 05:20:49 rhatcher Exp $");

// --- Constants ---

const int kVorCoarseness = 0; //Voronoi... coarseness....
const bool kVorIsUseEverywhere = true;

const float kVorGlitchCutoff = 5.0;

const Int_t kRectMap = 142;
//The file used in BfldValidate is 142 (more precisely, it is
//the file /home/minos/software/labyrinth/ephemera/bfield/bfld_142.dat

//From BfldValidate.cxx
const float kSpacing = 5 * Munits::cm;
const float kHalfLength = 375 * Munits::cm;
const int kDivs = 2* (int) (kHalfLength/kSpacing) + 1;

const BfldInterpMethod::InterpMethod_t kVorInterpMethod =
BfldInterpMethod::kPlanar;

const char kTrailFileName[30] = "../NuMuEventHits.dat";
const char kGeneratorsFileName[40] = "../NearDetBField.default";

const char kGridOutFileName[20] = "bfldgriddata.root";
const char kTrailOutFileName[20] = "bfldtraildata.root";

// --- Custom data types ---

class GenTester {
public:
  GenTester(const char *genFile);
  ~GenTester(void);
  TVector3 GetBFieldAtNearestGenerator(TVector3 &v);

private:
  float SumSquares(float x, float y) {
    return x * x + y * y;
  }

  static const int kDefaultNumGens = 3072;

  vector<TVector3> *fGens;
  vector<TVector3> *fBfld;
};


// --- Global variables ---

TROOT root("rootie","tootie");
GenTester *gGenTester = 0;

// --- function prototypes ---

void MakeGridData(const char *outFileName,bool isGen);
void MakeTrailData(const char *outFileName,const char *inFileName,
                   bool isGen);
void StoreDataPoint(TVector3 &r,BField &rectField,BField &vorField,
                    TNtuple *nt,bool isGen);

// --- main and other functions ---

int main(int argc,char **argv) {
  bool isGrid = false,
    isTrail = false,
    isGen = false,
    isDebug = false,
    isUseDefaultTrailFile = true;

  char inFileName[255];

  char opt = getopt(argc,argv,"rtghdf:");
  if(opt == EOF) { //defaults
    isGrid = true;
    isTrail = false;
    isGen = false;
    isDebug = false;
  }
  else {
    while(opt != EOF) {
      switch(opt) {
      case EOF:
        break;
      case '?':
        cout << "Invalid option. Use BfldDebugData -h for help." << endl;
        exit(1);
        break;
      case 'h':
        cout << "BfldDebugData [-rtghfd] -- generate some BField values."
             << endl;
        cout << "By default, generates Bfield values on a rectangular grid."
             << endl;
        cout << "-r: Take values on a rectangular grid. Output to file" << endl
             << "    bfldgriddata.root" << endl;
        cout << "-t: Take values on a muon trail. Can be used with" << endl;
        cout << "    -r or alone. Output to file bfldtraildata.root" << endl;
        cout << "-g: Generate nearest generator data" << endl;
        cout << "    Can be used with any of the other options." << endl;
        cout << "-f filename: With -t, use filename as the test point" << endl
             << "    file instead of the default." << endl;
        cout << "-d: Print out lots of debugging information." << endl;
        cout << "-h: show this help message." << endl;
        exit(0);
        break;
      case 'r':
        isGrid = true;
        break;
      case 't':
        isTrail = true;
        break;
      case 'g':
        isGen = true;
        break;
      case 'd':
        isDebug = true;
        break;
      case 'f':
        if(optarg) {
          strcpy(inFileName,optarg);
          isUseDefaultTrailFile = false;
        }
        break;
      default:
        cout << "Error processing arguments. Exiting." << endl;
        exit(1);
      } //switch(opt)
      opt = getopt(argc,argv,"rtghdf:");
    } //while(opt != EOF)
  } //if(opt == EOF) {...} else

  //We don't want all the information Msg normally prints out.
  MsgService::Instance()->GetStream("Bfldtest")->SetFormat(Msg::kDebug,0);
  MsgService::Instance()->GetStream("Bfldtest")->SetFormat(Msg::kWarning,0);

  if(isDebug) {
    MsgService::Instance()->GetStream("Bfldtest")->SetLogLevel(Msg::kDebug);
    MSG("Bfldtest",Msg::kDebug) << "Debugging info enabled." << endl;
  }

  if(isGen) {
    MSG("Bfldtest",Msg::kInfo)
      << "This run will generate nearest-generator data." << endl;
    gGenTester = new GenTester(kGeneratorsFileName);
  }
  else
    MSG("Bfldtest",Msg::kInfo)
      << "This run will NOT generate nearest-generator data." << endl;

  if(isGrid) {
    MSG("Bfldtest",Msg::kInfo) << "Using rectangular grid." << endl;      
    MakeGridData(kGridOutFileName,isGen);
  }
   
  if(isTrail) {
    MSG("Bfldtest",Msg::kInfo) << "Using realistic muon trail." << endl;
    if(isUseDefaultTrailFile)
      strcpy(inFileName,kTrailFileName);

    MakeTrailData(kTrailOutFileName,inFileName,isGen);
  }
  
  MSG("Bfldtest",Msg::kInfo) << "Done." << endl;

  return 0;
}

void MakeGridData(const char *outFileName,bool isGen) {
  TFile outFile(outFileName,"RECREATE");
  TNtuple *nt;

  if(isGen)
    nt = new TNtuple("data","data",
                     "x:y:bx1:by1:bz1:bx2:by2:bz2:bgx:bgy:bgz:suspect");
  else
    nt = new TNtuple("data","data","x:y:bx1:by1:bz1:bx2:by2:bz2:suspect");

  VldContext theContext1(Detector::kNear,SimFlag::kReroot,VldTimeStamp());
  VldContext theContext2(Detector::kNear,SimFlag::kReroot,VldTimeStamp());;

  MSG("Bfldtest",Msg::kInfo) << "Building rect2d interpolator..." << endl;
  BField rectField(theContext1,-1,kRectMap);
  rectField.SetInterpMethod(BfldInterpMethod::kDefault);

  MSG("Bfldtest",Msg::kInfo) << "Building Voronoi interpolator..." << endl;
  MSG("Bfldtest",Msg::kInfo) << "(BFL/Bfld warnings disabled)" << endl;
  MsgService::Instance()->GetStream("BFL")->SetLogLevel(Msg::kFatal);
  MsgService::Instance()->GetStream("Bfld")->SetLogLevel(Msg::kFatal);

  BField vorField(theContext2,kVorCoarseness,kVorIsUseEverywhere);
  vorField.SetInterpMethod(kVorInterpMethod);

  TVector3 v(0.0,0.0,0.0),b1,b2,bGen;
  float x,y,z;

  //Scan rectangular grid:
  float theRealHalf = kDivs * kSpacing / 2.0;

  Float_t dx = theRealHalf*2/kDivs;
  Float_t dy = theRealHalf*2/kDivs;

  MSG("Bfldtest",Msg::kInfo) << "Grid parameters:" << endl;
  MSG("Bfldtest",Msg::kInfo) << "xmin = " << -theRealHalf << "m; xmax = "
       << theRealHalf << "m; xstep = " << dx << "m" << endl;
  MSG("Bfldtest",Msg::kInfo) << "ymin = " << -theRealHalf << "m; ymax = "
       << theRealHalf << "m; ystep = " << dy << "m" << endl;

  for (Int_t ix = 0; ix < kDivs; ix++) {
    x = -theRealHalf + ((float)ix + 0.5)*dx;
    for (Int_t iy = 0; iy < kDivs; iy++) {
      y = -theRealHalf + ((float)iy + 0.5)*dy;
      
      v.SetX(x); v.SetY(y);

      StoreDataPoint(v,rectField,vorField,nt,isGen);
    }
  }
  
  MSG("Bfldtest",Msg::kInfo) << "Writing file " << outFileName << endl;
  outFile.Write();
  //Since BField doesn't inherit from TObject, only the ntuple will be
  //written to outFile by TFile::Write().
  delete nt;
}

void MakeTrailData(const char *outFileName,const char *inFileName,bool isGen)
{
  TFile outFile(outFileName,"RECREATE");

  TNtuple *nt;
  if(isGen)
    nt = new TNtuple("data","data",
                     "x:y:bx1:by1:bz1:bx2:by2:bz2:bgx:bgy:bgz:suspect");
  else
    nt = new TNtuple("data","data","x:y:bx1:by1:bz1:bx2:by2:bz2:suspect");

  VldContext theContext1(Detector::kNear,SimFlag::kReroot,VldTimeStamp());
  VldContext theContext2(Detector::kNear,SimFlag::kReroot,VldTimeStamp());;
  
  MSG("Bfldtest",Msg::kInfo) << "Building rect2d interpolator..." << endl;
  BField rectField(theContext1,-1,kRectMap);
  rectField.SetInterpMethod(BfldInterpMethod::kDefault);

  MSG("Bfldtest",Msg::kInfo) << "Building Voronoi interpolator..." << endl;
  MSG("Bfldtest",Msg::kInfo) << "(BFL/Bfld Warnings disabled)" << endl;
  MsgService::Instance()->GetStream("BFL")->SetLogLevel(Msg::kFatal);
  MsgService::Instance()->GetStream("Bfld")->SetLogLevel(Msg::kFatal);

  BField vorField(theContext2,kVorCoarseness,kVorIsUseEverywhere);
  vorField.SetInterpMethod(kVorInterpMethod);

  TVector3 v(0.0,0.0,0.0); // THIS IS a very important vector
  TVector3 b1,b2,bGen;
  float x,y,z;

  //Scan the trail for points at which to evaluate the B field.
  //Note that NuMuEventHits.dat, the default, gives locations in cm.

  MSG("Bfldtest",Msg::kInfo) << "Using trail file " << inFileName << endl;

  ifstream hitsFile(inFileName,ios::in);
  if(!hitsFile.good()) {
    MSG("Bfldtest",Msg::kWarning) << "Couldn't open trail file "
                                  << inFileName << endl;
    MSG("Bfldtest",Msg::kWarning) << "Bailing out..." << endl;

    delete nt;
    return;
  }

  int id;
  
  while(!hitsFile.eof()) {
    hitsFile >> id >> x >> y >> z;
    v.SetXYZ(x*Munits::cm,y*Munits::cm,z*Munits::cm);

    StoreDataPoint(v,rectField,vorField,nt,isGen);
  }
  
  MSG("Bfldtest",Msg::kInfo) << "Writing file " << outFileName << endl;
  outFile.Write();
  //Since BField doesn't inherit from TObject, only the ntuple will be
  //written to outFile by TFile::Write().

  delete nt;
}

void StoreDataPoint(TVector3 &r,BField &rectField,BField &vorField,
                    TNtuple *nt,bool isGen) {
  TVector3 b1,b2,bGen;

  b1 = rectField.GetBField(r);
  b2 = vorField.GetBField(r);
  
  if(isGen) {
    bGen = gGenTester->GetBFieldAtNearestGenerator(r);
  }

  MSG("Bfldtest",Msg::kDebug)
    << "Voronoi interpolated b = (" << b2.X() << "," << b2.Y() << ")"
    << endl;

  if(isGen) {
    nt->Fill(r.X(),r.Y(),b1.X(),b1.Y(),b1.Z(),b2.X(),b2.Y(),b2.Z(),
             bGen.X(),bGen.Y(),bGen.Z(),b2.Mag() > kVorGlitchCutoff);
  }
  else {
    nt->Fill(r.X(),r.Y(),b1.X(),b1.Y(),b1.Z(),b2.X(),b2.Y(),b2.Z(),
             b2.Mag() > kVorGlitchCutoff);
  }
  
  if(b2.Mag() > kVorGlitchCutoff) {
    MSG("Bfldtest",Msg::kDebug)
      << "Flagged a suspicious value from Voronoi interpolator: " << endl;
    MSG("Bfldtest",Msg::kDebug) << "B(" << r.X() << "," << r.Y() << ",0) ?= "
                               << b2.Mag() << endl;
  }
  
}

GenTester::GenTester(const char *genFileName) {
  MSG("Bfldtest",Msg::kInfo) << "Reading in generator file " << genFileName
                             << endl;

  ifstream genFile(genFileName,ios::in);
  if(!genFile.good()) {
    MSG("Bfldtest",Msg::kFatal) << "Couldn't open generator file "
                                << genFileName << endl;
    MSG("Bfldtest",Msg::kFatal) << "Bailing out entirely." << endl;
    exit(1);
  }

  float id, x, y, z, bx, by, bz;
  fGens = new vector<TVector3>(kDefaultNumGens);
  fBfld = new vector<TVector3>(kDefaultNumGens);

  int i = 0,n = fGens->size();
  do {
    genFile >> id >> x >> y >> z >> bx >> by >> bz;
    if(i >= n) {
      n *= 2;
      fGens->resize(n);
      fBfld->resize(n);
    }

    (*fGens)[i].SetXYZ(x,y,z);
    (*fBfld)[i].SetXYZ(bx,by,bz);
    ++i;
  } while(!genFile.eof() && !genFile.fail());

  //For some reason, the last point read in is garbage because
  //the last line of the file is read twice. Discard it.

  fGens->resize(i - 1);
  fBfld->resize(i - 1);

  MSG("Bfldtest",Msg::kInfo) << "Read in " << i - 1 << " generators..."
                             << endl;
}


GenTester::~GenTester(void) {
  delete fGens;
  delete fBfld;
}

TVector3 GenTester::GetBFieldAtNearestGenerator(TVector3 &v) {

  float bestDistSqr = SumSquares((*fGens)[0].X() - v.X(),
                                 (*fGens)[0].Y() - v.Y());
  int bestGen = 0;
  float distSqr;

  int n = fGens->size(), i;

  for(i = 1;i < n;i++) {
    distSqr = SumSquares(v.X() - (*fGens)[i].X(),v.Y() - (*fGens)[i].Y());
    if(distSqr < bestDistSqr) {
      bestGen = i;
      bestDistSqr = distSqr;
    }
  }

  MSG("Bfldtest",Msg::kDebug)
    << "Nearest generator to (" << v.X() << "," << v.Y() << "): ("
    << (*fGens)[bestGen].X() << "," << (*fGens)[bestGen].Y() << ")"
    << endl;

  return (*fBfld)[bestGen];
}

---