// --------------------------------------------------------------------------
//
// Create the field map v18a from the delivered ASCII file.
// The format of the text file is visible from the macro body.
//
// V. Friese   22/02/2019
//
// --------------------------------------------------------------------------

void createFieldMap18a()
{

  TString inFileName  = "CBMNewServ1Slice2m9WsBoxmod2_18GRD.dat";
  TString mapName     = "field_v18a";
  TString outFileName = "../" + mapName + ".root";
  
  // Uncompress ASCII file
  TString zipFileName = inFileName + ".bz2";
  TString command = "bzip2 -dk ";
  command += zipFileName;
  std::cout << "Unpacking " << zipFileName << std::endl;
  std::cout << command << std::endl;
  gSystem->Exec(command.Data());

  // Open input file
  std::ifstream inFile(inFileName);
  if ( ! inFile.is_open() ) {
    std::cout << "Could not open input file " << inFileName << std::endl;
    return;
  }
  std::cout << "\nUsing input file " << inFileName << std::endl;

  // First line: number of bins in x, y and z and step size
  Int_t nBinsX = 0;
  Int_t nBinsY = 0;
  Int_t nBinsZ = 0;
  Double_t step;
  inFile >> nBinsZ >> nBinsX >> nBinsY >> step;
  std::cout << "Number of bins: x " << nBinsX << ", y " << nBinsY << ", z "
      << nBinsZ << ", step size " << step << " cm" << std::endl;

  // Field map limits
  Double_t xMax = (nBinsX - 1) * step;
  Double_t yMax = (nBinsY - 1) * step;
  Double_t zMax = (nBinsZ - 1) * step;

  // Create temporary arrays with field values
  TArrayF bx(nBinsX*nBinsY*nBinsZ);
  TArrayF by(nBinsX*nBinsY*nBinsZ);
  TArrayF bz(nBinsX*nBinsY*nBinsZ);

  // Seven lines of description to be skipped
  std::string line;
  for (Int_t iRead = 0; iRead < 8; iRead++) getline(inFile, line);

  // Read field values
  Double_t x = 0.;
  Double_t y = 0.;
  Double_t z = 0.;
  Double_t rx = 0.;
  Double_t ry = 0.;
  Double_t rz = 0.;
  Double_t xNom = 0.;
  Double_t yNom = 0.;
  Double_t zNom = 0.;
  Int_t index = 0;
  for (Int_t ix = 0; ix < nBinsX; ix++) {
    std::cout << "\rReading file: " << Int_t(ix*100/nBinsX) << " %" << std::flush;
    xNom = ix * step;
    for(Int_t iy = 0; iy < nBinsY; iy++) {
      yNom = iy * step;
      for (Int_t iz = 0; iz < nBinsZ; iz++) {
        index = ix*nBinsY*nBinsZ + iy*nBinsZ + iz;
        zNom = iz * step;
        inFile >> x >> y >> z >> rx >> ry >> rz;
        x /= 10.; // mm->cm
        y /= 10.; // mm->cm
        z /= 10.; // mm->cm
        if ( ! inFile.good() ) {
          std::cout << "Error reading from file at ix " << ix << ", iy " << iy
              << ", iz " << iz << std::endl;
          inFile.close();
          return;
        } //? stream state
        if ( abs((x-xNom)/xNom) > 0.001 ) {
          std::cout << "Error in x: ix = " << ix << ", x = " << x << std::endl;
          inFile.close();
          return;
        } //? x ok
        if ( abs((y-yNom)/yNom) > 0.001 ) {
          std::cout << "Error in y: iy = " << iy << ", y = " << y << std::endl;
          inFile.close();
          return;
        } //? y ok
        if ( abs((z-zNom)/zNom) > 0.001 ) {
          std::cout << "Error in z: iz = " << iz << ", z = " << z << std::endl;
          inFile.close();
          return;
        } // z ok
        bx[index] = rx * 10.;  // T->kG
        by[index] = ry * 10.;  // T->kG
        bz[index] = rz * 10.;  // T->kG
      } //# iz
    } //# iy
  } //# ix
  std::cout << "\rReading file: done" << std::endl;
  inFile.close();

  auto field = new CbmFieldMapSym3();
  field->Init(nBinsX, 0., xMax, nBinsY, 0., yMax, nBinsZ, 0., zMax,
              &bx, &by, &bz);
  field->WriteRootFile(outFileName.Data(), mapName.Data());
  std::cout << "Created ROOT file " << outFileName << std::endl;
  
  gSystem->Exec(Form("rm -f %s", inFileName.Data()));

}