loon: BfldMapRect2d class Reference

00055    : BfldMap(grid,variant)
00056 {
00057 
00058    this->Clear();
00059 
00060    // read in a Map
00061 
00062    // strip off "-" sign if necessary
00063    Int_t varabs = TMath::Abs(variant);
00064 
00065    MSG("Bfld",Msg::kDebug) 
00066      << "BfldMapRect2d ctor map " << varabs << "." << endl;
00067 
00068    char string[255];
00069 
00070    // format and expand file name
00071    const char* search_path = ".:$BMAPPATH";
00072    sprintf(string,"bfld_%3d.dat",varabs);
00073    char* expanded = gSystem->Which(search_path,string,kReadPermission);
00074    if (!expanded) {
00075       MSG("Bfld",Msg::kError) << "can not find input file: '" 
00076                               << fDataSource <<"'"<< endl;
00077       return;
00078    }
00079    fDataSource = expanded;
00080 
00081    delete [] expanded;  // user responsible for deleting returned string
00082 
00083    // open the file for reading
00084    ifstream from;
00085    from.open(fDataSource.c_str(),ios::in);
00086    if (!from) {
00087       MSG("Bfld",Msg::kError) << "can not open input file: '" 
00088                               << fDataSource <<"'"<< endl;
00089       return;
00090    }
00091 
00092    // read in the header line
00093    from.getline(string,sizeof(string),'\n');
00094    fHeader = string;
00095 
00096    MSG("Bfld",Msg::kInfo) 
00097      << "BfldMapRect2d read file: " << fDataSource << endl;
00098    MSG("Bfld",Msg::kInfo) 
00099      << "BfldMapRect2d: " << fHeader << endl;
00100 
00101    fIsUVZ = ( fHeader.find("UVZ") != string::npos );
00102    Double_t dirScale  = 1;
00103 
00104    // search the title string for a coil current
00105    // generally a token like 15kAturns, 40kAturns
00106    string::size_type loc_aturn = fHeader.find("Aturn");
00107    // alas sometimes the header has the value string "40k A-turns"
00108    if (loc_aturn == string::npos) loc_aturn = fHeader.find("A-turn");
00109 
00110    if (loc_aturn != string::npos) {
00111      // find last space before 'Aturn'
00112      string::size_type loc_space = fHeader.find_last_of(" ",loc_aturn); 
00113      std::string sval = fHeader.substr(loc_space,loc_aturn-loc_space);
00114      MSG("Bfld",Msg::kDebug) << " found Aturns: '" << sval << "'" << endl;
00115      string::size_type loc_k = sval.find("k");
00116      if (loc_k != string::npos) {
00117        sval.replace(loc_k,1," ");
00118        MSG("Bfld",Msg::kDebug) << " found 'k' new string '" << sval << "'" << endl;
00119      }
00120      std::istringstream inputString(sval);
00121      inputString >> fGeneratedCoilCurrent;
00122      if (loc_k != string::npos) fGeneratedCoilCurrent *= 1000;
00123 
00124      if ( fHeader.find("CurrentForward") == string::npos ) {
00125        MSG("Bfld",Msg::kInfo) 
00126          << "BfldMapRect2d: maps were created focusing mu+;"
00127          << " correct this with an extra -1" << endl;
00128        dirScale = -1;
00129      }   
00130 
00131      MSG("Bfld",Msg::kInfo) 
00132        << "BfldMapRect2d: GeneratedCoilCurrent=" 
00133        << setprecision(4) << (fGeneratedCoilCurrent/1000.)
00134        << "kA*turns based on header info."
00135        << resetiosflags(ios::floatfield)  // reset to "%g" format
00136        << endl;
00137    }
00138    else {
00139      MSG("Bfld",Msg::kError) 
00140        << "BfldMapRect2d: could not find coil current in input file header "
00141        << endl;
00142    }
00143 
00144    from >> fX0 >> fDx >> fNx >> fY0 >> fDy >> fNy;
00145    // lengths read in from map files are in cm, convert to base units
00146    fX0 *= Munits::cm;
00147    fDx *= Munits::cm;
00148    fY0 *= Munits::cm;
00149    fDy *= Munits::cm;
00150 
00151    MSG("Bfld",Msg::kInfo) 
00152      << "BfldMapRect2d:"
00153      << setprecision(6)
00154      << " x0 " << fX0 << " dx " << fDx << " nx " << fNx 
00155      << " y0 " << fY0 << " dy " << fDy << " ny " << fNy 
00156      << resetiosflags(ios::floatfield)  // reset to "%g" format
00157      << endl;
00158 
00159    // next 3 elements are irrelevant (z0,dz,nz)
00160    Float_t z0, dz;
00161    Int_t   nz;
00162    from >> z0 >> dz >> nz;
00163 
00164    // read in remainder of the line for any quadflags
00165    from.getline(string,sizeof(string),'\n');
00166    TString quads = string;
00167 
00168    // replace "x" with " "
00169    quads = quads.ReplaceAll("x"," ");
00170    // quads should be of the form ' 00 00 00 00  FNAL2'
00171    MSG("Bfld",Msg::kDebug) << "quads string '" 
00172                            << quads << "'" << endl;
00173 
00174    fFileForm = 0;
00175    if (quads.Contains("FNAL"))  fFileForm = 1;
00176    if (quads.Contains("FNAL2")) fFileForm = 2;
00177 
00178    // remove any "FNAL" flag
00179    if (quads.Length() > 12) quads = quads.Remove(12);
00180 
00181    // scan the remaining part for the quad flags
00182    if (quads.Length() > 0) {
00183       sscanf(quads,"%3d%3d%3d%3d",
00184              fQuadFlags,fQuadFlags+1,fQuadFlags+2,fQuadFlags+3);
00185    } else {
00186       MSG("Bfld",Msg::kInfo) << "no quads info - use defaults" << endl;
00187       fQuadFlags[0] = 0;
00188       fQuadFlags[1] = 2;
00189       fQuadFlags[2] = 3;
00190       fQuadFlags[3] = 1;
00191    }
00192 
00193    // if x0=0,y0=0 assume map only has one quadrant 
00194    fQuadrant = ( 0 == fX0 && 0 == fY0 );
00195 
00196    MSG("Bfld",Msg::kSynopsis) 
00197      << "BfldMapRect2d: fileform " << fFileForm
00198      << "  quadflags (" 
00199      << fQuadFlags[0] << "," 
00200      << fQuadFlags[1] << "," 
00201      << fQuadFlags[2] << "," 
00202      << fQuadFlags[3] << ")  use them " 
00203      << (fQuadrant ? "'true'" : "'false'") 
00204      << endl;
00205 
00206    ResizeVectors(fNx*fNy);
00207 
00208    // build a temporary mesh so we can can use the functions
00209    // for finding the generator from x,y etc
00210    // but this isn't the one that will go into the loan pool
00211    BfldMeshRect2d tmpMesh(this); 
00212    
00213    Int_t   ii, indx;
00214    Float_t x,y,z;
00215    Float_t bx,by,bz,bmag;
00216 
00217    const Float_t GaussToKGauss = 0.001;
00218    const Float_t TeslaToKGauss = 10.;
00219 
00220    std::vector<bool> seen_gen(fNx*fNy);
00221 
00222    // start reading the file
00223    Int_t ilines=0;
00224    //original for (ilines=0; ilines < fNx*fNy; ilines++) {
00225    for (ilines=0; true ; ilines++) {
00226 
00227       switch (fFileForm) {
00228       case 0:
00229          from >> x >> y >> z >> bx >> by >> bz >> bmag;
00230          // LLNL maps were in gauss ... need kGauss
00231          bx   = bx   * GaussToKGauss;
00232          by   = by   * GaussToKGauss;
00233          bz   = bz   * GaussToKGauss;
00234          bmag = bmag * GaussToKGauss;
00235          break;
00236       case 1:
00237          from >> ii >> x >> y >> bx >> by;
00238          // FNAL (version 1) maps were in Tesla
00239          bx   = bx   * TeslaToKGauss;
00240          by   = by   * TeslaToKGauss;
00241          bz   = 0.0;
00242          // avoid unitialized values
00243          z    = 0.0;
00244          bmag = TMath::Sqrt(bx*bx+by*by+bz*bz);
00245          break;
00246       case 2:
00247          from >> x >> y >> z >> bx >> by >> bz >> bmag;
00248          // FNAL (version 2) maps were in Tesla
00249          bx   = bx   * TeslaToKGauss;
00250          by   = by   * TeslaToKGauss;
00251          bz   = bz   * TeslaToKGauss;
00252          bmag = bmag * TeslaToKGauss;
00253          break;
00254       default:
00255          MSG("Bfld",Msg::kFatal) << " bad FileForm " << fFileForm << endl;
00256          break;
00257       }
00258 
00259       if (!from.good()) {
00260         if (ilines != fNx*fNy)  //  read one too far ...
00261           MSG("Bfld",Msg::kInfo) 
00262             << "BfldMapRect2d: early EOF on " << endl 
00263             << "  "
00264             << fDataSource 
00265             << endl
00266             << "  saw " << ilines << " entries, expected " 
00267             << fNx*fNy << "; could be due to suppressed (0,0) points."
00268             << endl;
00269         break;
00270       }
00271 
00272       // The above puts field values in kGauss
00273       // convert them to base units
00274       // Also make correction for reversed current direction when generated
00275       bx   *= Munits::kilogauss * dirScale;
00276       by   *= Munits::kilogauss * dirScale;
00277       bz   *= Munits::kilogauss * dirScale;
00278       bmag *= Munits::kilogauss;
00279       // positions read in are given in cm
00280       x    *= Munits::cm;
00281       y    *= Munits::cm;
00282       z    *= Munits::cm;
00283 
00284       // all the maps aren't in consistent order (ascending/descending x,y)
00285       // so we have to actually interpret the (x,y) values to determine 
00286       // the array index.
00287       // we'll put them in the array in a way that we can lookup 
00288       // based on fX0,fDx,fNx,fY0,fDy,fNy
00289 
00290       // when reading points, find the nearest generator
00291       // not the lower left one
00292       indx = tmpMesh.NearestGenerator(x,y);
00293 
00294       if ( seen_gen[indx] ) {
00295         MAXMSG("Bfld",Msg::kWarning,20)
00296           << "Duplicate entry for indx " << indx
00297           << " x= " << x << ", y= " << y << endl;
00298       }
00299       else {
00300         seen_gen[indx] = true;
00301       }
00302 
00303       // cross check! for bad data points
00304       Int_t jxx, jyy;
00305       tmpMesh.GeneratorToIndices(indx,jxx,jyy);
00306       Float_t xx,yy;
00307       xx = fX0 + jxx*fDx;
00308       yy = fY0 + jyy*fDy;
00309       const Float_t dxy_epsilon = 0.0005;
00310       static MsgFormat f9p6("f9.6"); 
00311       static MsgFormat f5p1("f5.1"); 
00312       if ( TMath::Abs(x-xx) > dxy_epsilon*fDx ||
00313            TMath::Abs(y-yy) > dxy_epsilon*fDy    ) {
00314         MAXMSG("Bfld",Msg::kWarning,20)
00315           << " [" << setw(5) << ilines << "] off grid? "
00316           << endl 
00317           <<   "  (x,y)=(" << f9p6(x)  << "," << f9p6(y)  << ") "
00318           << " expected (" << f9p6(xx) << "," << f9p6(yy) << ") "
00319           << endl
00320           << "   error is " << f5p1(100.*TMath::Abs(x-xx)/fDx) << "% of dx "
00321           <<           ", " << f5p1(100.*TMath::Abs(y-yy)/fDy) << "% of dy "
00322           << endl;
00323       }
00324 
00325       AddGenerator(indx,bx,by,bz);
00326 
00327 #ifdef BLAHBLAH
00328       // one *might* think of using [i] rather than [indx]
00329       if (indx != ilines &&  ilines < 2*fNx ) {
00330          MSG("Bfld",Msg::kWarning) 
00331            << " i " << i << " indx " << indx
00332            << " (" << x << "," << y << ")" << endl;
00333       }
00334 #endif
00335 
00336    }
00337 
00338    MSG("Bfld",Msg::kDebug) 
00339      << "BfldMapRect2d: read " << ilines << " entries from "
00340      << fDataSource << endl;
00341 
00342    fFilledOkay = kTRUE;
00343 }
BfldMapRect2d Class Reference

Public Member Functions

Protected Attributes

Private Member Functions

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	BfldMapRect2d ()
	BfldMapRect2d (BfldGrid::Grid_t grid, Int_t variant)
virtual	~BfldMapRect2d ()
Float_t	GetX0 () const
Float_t	GetY0 () const
Float_t	GetDx () const
Float_t	GetDy () const
Int_t	GetNx () const
Int_t	GetNy () const
Bool_t	GetIsUVZ () const
const Int_t *	GetQuadFlags () const
Protected Attributes
std::string	fHeader
Float_t	fX0
Float_t	fDx
Int_t	fNx
Float_t	fY0
Float_t	fDy
Int_t	fNy
Bool_t	fIsUVZ
Bool_t	fQuadrant
Int_t	fQuadFlags [4]
Int_t	fFileForm
Private Member Functions
void	Clear (const Option_t *="")