TPC/SCDCorrections/FieldMaps.cxx

0001 #include "FieldMaps.h"
0002 #include <string>
0003
0004 #include "TH3F.h"
0005 #include "TFile.h"
0006 #include "TMath.h"
0007
0008 //=====================
0009 FieldMaps::FieldMaps() {
0010   fDebug = 0;
0011   fRadialBin = -1;
0012   fEr = NULL;
0013   fEp = NULL;
0014   fEz = NULL;
0015   fRho = NULL;
0016   fInnerRadius=1.;
0017   fOutterRadius=2.;
0018   fHalfLength=1.;
0019   fNRadialSteps=1;
0020   fNAzimuthalSteps=1;
0021   fNLongitudinalSteps=1;
0022   fFileNameRoot="rho";
0023   fLSNameRoot="none";
0024   fMirrorZ=false;
0025 }
0026 //=====================
0027 FieldMaps::~FieldMaps() {
0028 }
0029 //=====================
0030 void FieldMaps::Make(int n) {
0031   fRadialBin = n;
0032   InitMaps();
0033   ComputeE();
0034   SaveMaps();
0035 }
0036 //=====================
0037 void FieldMaps::InitMaps() {
0038   if(fDebug>0) printf("FieldMaps is initializing field maps... ");
0039   fEr = new TH3F("Er","Er [V/cm];Radial [cm];Azimuthal [rad];Longitudinal [cm]",
0040          fNRadialSteps,fInnerRadius,fOutterRadius,
0041          fNAzimuthalSteps,0,TMath::TwoPi(),
0042          fNLongitudinalSteps,-fHalfLength,+fHalfLength);
0043   fEp = new TH3F("Ep","Ephi [V/cm];Radial [cm];Azimuthal [rad];Longitudinal [cm]",
0044          fNRadialSteps,fInnerRadius,fOutterRadius,
0045          fNAzimuthalSteps,0,TMath::TwoPi(),
0046          fNLongitudinalSteps,-fHalfLength,+fHalfLength);
0047   fEz = new TH3F("Ez","Ez [V/cm];Radial [cm];Azimuthal [rad];Longitudinal [cm]",
0048          fNRadialSteps,fInnerRadius,fOutterRadius,
0049          fNAzimuthalSteps,0,TMath::TwoPi(),
0050          fNLongitudinalSteps,-fHalfLength,+fHalfLength);
0051   if(fDebug>0) printf("[DONE]\n");
0052   ReadFile();
0053 }
0054 //=====================
0055 void FieldMaps::ReadFile() {
0056   const char *inputfile= Form("%s_0.root",fFileNameRoot.data());
0057   if(fDebug>0) printf("FieldMaps is reading from %s... ",inputfile);
0058   TFile *ifile = new TFile(inputfile);
0059   if(ifile->IsZombie()) ifile=NULL;
0060   if(!ifile) exit(1);
0061   fRho = (TH3F*) ifile->Get("rho");
0062   if(fDebug>0) printf("[DONE]\n");
0063   // PLUG Gr Gz HERE
0064 }
0065 //=====================
0066 void FieldMaps::SaveMaps() {
0067   TString filename;
0068   if(fRadialBin<0) filename = Form("%s_1.root",fFileNameRoot.data());
0069   else filename = Form("%s_1_%d.root",fFileNameRoot.data(),fRadialBin);
0070   if(fDebug>1) printf("FieldMaps saving field maps into %s... ",filename.Data());
0071   TFile *ofile = new TFile(filename.Data(),"RECREATE");
0072   ofile->WriteObject(fEr,"Er");
0073   ofile->WriteObject(fEp,"Ep");
0074   ofile->WriteObject(fEz,"Ez");
0075   ofile->Close();
0076   if(fDebug>1) printf("[DONE]\n");
0077 }
0078 //=====================
0079 float FieldMaps::ReadCharge(float rprime, float phiprime, float zprime, float dr, float dphi,float dz) {
0080   int brmin = fRho->GetXaxis()->FindBin(rprime - dr/2 + 1e-8);
0081   int brmax = fRho->GetXaxis()->FindBin(rprime + dr/2 - 1e-8);
0082   int bpmin = fRho->GetYaxis()->FindBin(phiprime - dphi/2 + 1e-8);
0083   int bpmax = fRho->GetYaxis()->FindBin(phiprime + dphi/2 - 1e-8);
0084   int bzmin = fRho->GetZaxis()->FindBin(zprime - dz/2 + 1e-8);
0085   int bzmax = fRho->GetZaxis()->FindBin(zprime + dz/2 - 1e-8);
0086   if(brmin<1||bpmin<1||bzmin<1) return 0;
0087   if(brmax<brmin) brmax = brmin +1;
0088   if(bpmax<bpmin) bpmax = bpmin +1;
0089   if(bzmax<bzmin) bzmax = bzmin +1;
0090
0091   float ddr = fRho->GetXaxis()->GetBinWidth(1);
0092   float ddp = fRho->GetYaxis()->GetBinWidth(1);
0093   float ddz = fRho->GetZaxis()->GetBinWidth(1);
0094   float totalQ = 0;
0095   for(int br=brmin; br!=brmax+1; ++br)
0096     for(int bp=bpmin; bp!=bpmax+1; ++bp)
0097       for(int bz=bzmin; bz!=bzmax+1; ++bz) {
0098     float rho = fRho->GetBinContent(br,bp,bz);
0099     float r = fRho->GetXaxis()->GetBinCenter( br );
0100     float dv = r*ddp*ddr*ddz;
0101     totalQ += rho*dv;
0102       }
0103   //return rprime*fRho->Integral( brmin,brmax , bpmin,bpmax , bzmin,bzmax );
0104   return totalQ;
0105 }