PHGenFitPkg/GenFitExp/Field.cc

0001 /*
0002    Authors: Haiwang Yu
0003 */
0004 #include "Field.h"
0005
0006 #include <phfield/PHField.h>
0007
0008 #include <TVector3.h>                   // for TVector3
0009
0010 #include <CLHEP/Units/SystemOfUnits.h>
0011
0012 #include <cassert>
0013 #include <limits>
0014
0015 namespace genfit
0016 {
0017 Field::Field(const PHField* field)
0018   : field_(field)
0019 {
0020   assert(field_);
0021 }
0022
0023 //
0024 //bool Field::initialize(std::string inname)
0025 //{
0026 //  field_map_r_ = new TH2D("field_map_r_", "B_{r} [kGauss]; z [cm]; r [cm]", 401, -401, 401, 151, -1, 301);
0027 //  field_map_r_->SetDirectory(0);
0028 //  field_map_z_ = new TH2D("field_map_z_", "B_{z} [kGauss]; z [cm]; r [cm]", 401, -401, 401, 151, -1, 301);
0029 //  field_map_z_->SetDirectory(0);
0030 //  TFile* fin = TFile::Open(inname.data(), "READ");
0031 //  if (!fin)
0032 //  {
0033 //    LogERROR("Input TFile is invalid!");
0034 //    return false;
0035 //  }
0036 //  TNtuple* T = (TNtuple*) fin->Get("fieldmap");
0037 //  if (!T)
0038 //  {
0039 //    LogERROR("Input filed map NTuple not found!");
0040 //    ;
0041 //    return false;
0042 //  }
0043 //
0044 //  Float_t r;
0045 //  Float_t z;
0046 //  Float_t br;
0047 //  Float_t bz;
0048 //  T->SetBranchAddress("r", &r);
0049 //  T->SetBranchAddress("z", &z);
0050 //  T->SetBranchAddress("br", &br);
0051 //  T->SetBranchAddress("bz", &bz);
0052 //
0053 //  for (long ientry = 0; ientry < T->GetEntries(); ientry++)
0054 //  {
0055 //    //if (T->GetEntry(ientry) < 0) continue;
0056 //    T->GetEntry(ientry);
0057 //    //std::cout<<r <<" ,"<<z <<" ,"<<br<<" ,"<<bz<<"\n";
0058 //    field_map_r_->SetBinContent(
0059 //        field_map_r_->GetXaxis()->FindBin(z),
0060 //        field_map_r_->GetYaxis()->FindBin(r),
0061 //        br * 10);
0062 //    field_map_z_->SetBinContent(
0063 //        field_map_z_->GetXaxis()->FindBin(z),
0064 //        field_map_z_->GetYaxis()->FindBin(r),
0065 //        bz * 10);
0066 //  }
0067 //  fin->Close();
0068 //
0069 //  //field_map_z_->Print("all");
0070 //
0071 //  return true;
0072 //}
0073
0074 //void Field::plot(std::string option){
0075 //
0076 ////    TH2D *hbr = new TH2D("hbr","|B_{r}| [kGauss]; z [cm]; r [cm]",401, -401, 401, 151, -1, 301);
0077 ////    TH2D *hbz = new TH2D("hbz","|B_{z}| [kGauss]; z [cm]; r [cm]",401, -401, 401, 151, -1, 301);
0078 //  TH2D *hbr = new TH2D("hbr","|B_{r}| [kGauss]; z [cm]; r [cm]",300, -350, 350, 200, 0, 290);
0079 //  TH2D *hbz = new TH2D("hbz","|B_{z}| [kGauss]; z [cm]; r [cm]",300, -350, 350, 200, 0, 290);
0080 //  for (double r = 0; r < 300; r+=1)
0081 //      for (double z = 400; z > -400; z-=1) {
0082 //          double bx;
0083 //          double by;
0084 //          double bz;
0085 //          get(r, 0, z, bx, by, bz);
0086 //
0087 ////            double br = sqrt(bx * bx + by * by);
0088 //
0089 //          //std::cout <<"DEBUG: "<<__LINE__<<":"<< r << " ," << z << " ," << br << " ," << bz << "\n";
0090 //
0091 //          hbr->SetBinContent(
0092 //                  hbr->GetXaxis()->FindBin(z),
0093 //                  hbr->GetYaxis()->FindBin(r),
0094 //                  bz*10);
0095 //          hbz->SetBinContent(
0096 //                  hbz->GetXaxis()->FindBin(z),
0097 //                  hbz->GetYaxis()->FindBin(r),
0098 //                  bz*10);
0099 //      }
0100 //
0101 //  TCanvas *c0 = new TCanvas("c0","c0");
0102 //  c0->Divide(1,2);
0103 //  c0->cd(1);
0104 //  hbr->SetStats(0);
0105 //  hbr->Draw("colz");
0106 //  c0->cd(2);
0107 //  hbz->SetStats(0);
0108 //  hbz->Draw("colz");
0109 //  c0->Update();
0110 //
0111 //  c0->SaveAs("Field_plot.root");
0112 //  c0->SaveAs("Field_plot.pdf");
0113 //
0114 //  delete c0;
0115 //  delete hbr;
0116 //  delete hbz;
0117 //
0118 //  return;
0119 //}
0120
0121 TVector3 Field::get(const TVector3& v) const
0122 {
0123   double x = v.x();
0124   double y = v.y();
0125   double z = v.z();
0126   double Bx;
0127   double By;
0128   double Bz;
0129   get(x, y, z, Bx, By, Bz);
0130   return TVector3(Bx, By, Bz);
0131 }
0132
0133 void Field::get(const double& x, const double& y, const double& z, double& Bx, double& By, double& Bz) const
0134 {
0135   assert(field_);
0136
0137   const double Point[] = {x * CLHEP::cm, y * CLHEP::cm, z * CLHEP::cm, 0};
0138   double Bfield[] = {std::numeric_limits<double>::signaling_NaN(),
0139                      std::numeric_limits<double>::signaling_NaN(),
0140                      std::numeric_limits<double>::signaling_NaN(),
0141                      std::numeric_limits<double>::signaling_NaN(),
0142                      std::numeric_limits<double>::signaling_NaN(),
0143                      std::numeric_limits<double>::signaling_NaN()};
0144
0145   field_->GetFieldValue(Point, Bfield);
0146
0147   Bx = Bfield[0] / CLHEP::kilogauss;
0148   By = Bfield[1] / CLHEP::kilogauss;
0149   Bz = Bfield[2] / CLHEP::kilogauss;
0150
0151   //  double r = sqrt(x * x + y * y);
0152   //
0153   //  if (fabs(r) > 300 || fabs(z) > 400)
0154   //  {
0155   //    Bx = 0;
0156   //    By = 0;
0157   //    Bz = 0;
0158   //    return;
0159   //  }
0160   //
0161   //  int bin_z = field_map_r_->GetXaxis()->FindBin(z);
0162   //  int bin_r = field_map_r_->GetYaxis()->FindBin(r);
0163   //
0164   //  Bz = field_map_z_->GetBinContent(bin_z, bin_r);
0165   //  double Br = field_map_r_->GetBinContent(bin_z, bin_r);
0166   //
0167   //  Bx = x / r * Br;
0168   //  By = y / r * Br;
0169   //
0170   //  if (r == 0)
0171   //  {
0172   //    Bx = 0;
0173   //    By = 0;
0174   //  }
0175
0176   //std::cout<<"DEBUG: "<<__LINE__<<": "<<z<<","<<r<<","<<bin_z<<","<<bin_r<<": "<<Br<<","<<Bz<<"\n";
0177 }
0178
0179 } /* End of namespace genfit */