ParameterizeBField_codegen.C

Go to the documentation of this file.

//
//
// This macro is used to generate ASCII files containing a parameterized 
// magnetic field map in a format suitable for use with the JANA
// JCalibrationFile class. It takes as input the root files 
// BfieldParameters_Bz.root and BfieldParameters_Bx.root which themselves
// are produced by the ParameterizeBField.C macro using bfield.root
// file as input. (The bfield.root file is generated with the bfield2root
// utility from a full map already in the calibDB).
//

//------------------------------
// ParameterizeBField_codegen
//------------------------------
void ParameterizeBField_codegen(const char *based_on="_1500_poisson_20090814_01")
{
   gROOT->Reset();

   // Create master params file
   stringstream fname;
   fname<<"solenoid"<<based_on<<"_params";
   ofstream fout_master(fname.str().c_str());

   // Write header info to file
   time_t now = time(NULL);
   fout_master<<"#"<<endl;
   fout_master<<"# Solenoidal field parameterization"<<endl;
   fout_master<<"#"<<endl;
   fout_master<<"# This file auto-generated by ParameterizeBField_codegen.C macro."<<endl;
   fout_master<<"# "<<ctime(&now);
   fout_master<<"# Produced by "<<getenv("USER")<<" on "<<getenv("HOST")<<endl;
   fout_master<<"#"<<endl;
   fout_master<<"#      input based on: "<<based_on<<endl;
   fout_master<<"#"<<endl;    
   fout_master<<"# This is the master parameters list. The parameterization of the"<<endl;      
   fout_master<<"# field is broken up into several sections in z. Each section has"<<endl;
   fout_master<<"# its own set of parameters kept in its own namepath."<<endl;
   fout_master<<"#"<<endl;
   fout_master<<"# Each line below represents a section. Sections are either representing"<<endl;
   fout_master<<"# Bx or Bz. Given in the table are the limits over with the section"<<endl;
   fout_master<<"# covers."<<endl;
   fout_master<<"#"<<endl;
   fout_master<<"#"<<endl;
   fout_master<<"#% Bi  sec  zmin  zmax  zmid  znorm  rmin  rmax  rmid  rnorm  namepath"<<endl;

   // Write indivdual params files, adding entries to master file
   WriteFile("BfieldParameters_Bz.root", "Bz", fout_master, fname);
   WriteFile("BfieldParameters_Bx.root", "Bx", fout_master, fname);
   
   // Close master params file
   fout_master.close();
}

//------------------------------
// WriteFile
//------------------------------
void WriteFile(const char *input_fname, string element, ofstream &fout_master, stringstream &fname)
{
#define MAX_SECTIONS 100

   // Open input ROOT file
   TFile *f = new TFile(input_fname);
   
   cout<<"------------------- Ignore errors in this section --------------------"<<endl;

   // Get zmin, zmax, rmin, and rmax for each section
   double zmin[MAX_SECTIONS], zmax[MAX_SECTIONS];
   double rmin[MAX_SECTIONS], rmax[MAX_SECTIONS];
   TH1D *z_bounds_hist = (TH1D*)gROOT->FindObject("z_bounds_hist");
   UInt_t Nsec = z_bounds_hist->GetNbinsX()-1;
   for(UInt_t sec=1; sec<=Nsec; sec++){
      zmin[sec-1] = z_bounds_hist->GetBinContent(z_bounds_hist->FindBin(sec));
      zmax[sec-1] = z_bounds_hist->GetBinContent(z_bounds_hist->FindBin(sec+1));

      stringstream hname;
      hname<<"sec"<<(sec)<<"_p0";   // sec1_p0, sec2_p0, ...
      TH1D *sec_p0 = (TH1D*)gROOT->FindObject(hname.str().c_str());
      if(!sec_p0)break;
      rmin[sec-1] = sec_p0->GetXaxis()->GetXmin();
      rmax[sec-1] = sec_p0->GetXaxis()->GetXmax();
   }

   // Get order of Chebyshev polynomials used for 1st level of parameterization
   UInt_t order1;
   for(order1=0; order1<20; order1++){
      stringstream hname;
      hname<<"sec1"<<"_p"<<order1;
      TH1D *sec1_p = (TH1D*)gROOT->FindObject(hname.str().c_str());
      if(!sec1_p){order1--; break;}
   }

   // Get order of Chebyshev polynomials used for 2nd level of parameterization
   TH1D *sec1_pp0 = (TH1D*)gROOT->FindObject("sec1_pp0");
   UInt_t order2 = sec1_pp0->GetNbinsX()-1;

   cout<<"----------------------------------------------------------------------"<<endl;

   // Information messages
   time_t now = time(NULL);
   cout<<endl;
   cout<<"Writing parameters for: "<<element<<endl;
   cout<<"       sections = "<<Nsec<<endl;
   cout<<"1st level order = "<<order1<<endl;
   cout<<"2nd level order = "<<order2<<endl;
   cout<<endl;

   // Parameters
   for(UInt_t sec=1; sec<=Nsec; sec++){
      
      // Open output file
      stringstream my_fname;
      my_fname<<fname.str()<<"_"<<element<<"_sec"<<sec;
      ofstream fout(my_fname.str().c_str());
      cout<<"Writing parameters to "<<fname.str()<<endl;
      
      // Add line to master file
      fout_master<<element;
      fout_master<<"\t"<<sec;
      fout_master<<"\t"<<zmin[sec-1];
      fout_master<<"\t"<<zmax[sec-1];
      fout_master<<"\t"<<(zmax[sec-1]+zmin[sec-1])/2.0;
      fout_master<<"\t"<<(zmax[sec-1]-zmin[sec-1])/2.0;
      fout_master<<"\t"<<rmin[sec-1];
      fout_master<<"\t"<<rmax[sec-1];
      fout_master<<"\t"<<(rmax[sec-1]+rmin[sec-1])/2.0;
      fout_master<<"\t"<<(rmax[sec-1]-rmin[sec-1])/2.0;
      fout_master<<"\t"<<"Magnets/Solenoid/"<<my_fname.str();
      fout_master<<endl;
      
      // Write header info to file
      fout<<"#"<<endl;
      fout<<"# Solenoidal field parameterization"<<endl;
      fout<<"#"<<endl;
      fout<<"# This file auto-generated by ParameterizeBField_codegen.C macro."<<endl;
      fout<<"# "<<ctime(&now);
      fout<<"# Produced by "<<getenv("USER")<<" on "<<getenv("HOST")<<endl;
      fout<<"#"<<endl;
      fout<<"#      input file: "<<f->GetName()<<endl;
      fout<<"#         section: "<<sec<<"  (out of "<<Nsec<<")"<<endl;
      fout<<"# 1st level order: "<<order1<<endl;
      fout<<"# 2nd level order: "<<order2<<endl;
      fout<<"#"<<endl;     
      fout<<"# In the following table, values are coefficients of Chebyshev"<<endl;    
      fout<<"# polynomials. There are 2 levels of parameterization corresponding"<<endl;     
      fout<<"# to the 2 dimensions in which the field map is parameterized (z and r)."<<endl;      
      fout<<"#"<<endl;
      fout<<"# The first level fits the field as a function of z in lab coordinates."<<endl;
      fout<<"# The second level parametrizes those coefficients as a function of r"<<endl;
      fout<<"# also in lab coordinates."<<endl;
      fout<<"#"<<endl;
      fout<<"#"<<endl;
      fout<<"#"<<endl;

      int chars_per_val = 14;
      string header((order2+1)*chars_per_val, ' ');
      for(UInt_t j=0; j<=order2; j++){
         stringstream name;
         name<<"p"<<j;
         string &s = name.str();
         header.replace((j+1)*chars_per_val-s.length()-2, s.length(), s);
      }
      fout<<"#"<<header<<endl;

      // Loop
      for(UInt_t i=0; i<=order1; i++){
         stringstream hname;
         hname<<"sec"<<sec<<"_pp"<<i;
         TH1D *h = (TH1D*)gROOT->FindObject(hname.str().c_str());
         
         string line((order2+1)*chars_per_val, ' ');
         for(UInt_t j=0; j<=order2; j++){
            stringstream val;
            val<<h->GetBinContent(j);
            string &s = val.str();
            if(s.length()>=chars_per_val){
               cerr<<"ERROR: too many characters in val! ("<<__FILE__<<":"<<__LINE__<<")"<<endl;
               return -1;
            }
            line.replace((j+1)*chars_per_val-s.length(), s.length(), s);
         }
         fout<<line<<endl;
      }
      
      fout.close();
   }
}