TPC/SCDCorrections/FieldMapsLaplace.cxx

0001 #include "FieldMapsLaplace.h"
0002 #include "LaplaceSolution.h"
0003 #include <string>
0004 #include <iostream>
0005 #include "TH3F.h"
0006 #include "TAxis.h"
0007 #include "TFile.h"
0008 #include "TMath.h"
0009
0010
0011 //=====================
0012 void FieldMapsLaplace::ComputeE() {
0013   const float prec=1e-8;
0014
0015   //------------
0016   //STEER
0017   if(fDebug>0) printf("FieldMaps is being computed from Laplace solutions... \n");
0018   LaplaceSolution *laplace_green_solution;
0019   if(0) {
0020     laplace_green_solution = new LaplaceSolution(fLSNameRoot);
0021   } else {
0022     laplace_green_solution = new LaplaceSolution(fInnerRadius/100/*[m]*/,fOutterRadius/100/*[m]*/,fHalfLength/100/*[m]*/);
0023   }
0024   int brprimeINI=0;
0025   int brprimeEND=fNRadialSteps;
0026   if(fRadialBin>-0.5) {
0027     brprimeINI = fRadialBin;
0028     brprimeEND = fRadialBin+1;
0029   }
0030   // constructing volume element
0031   float dr = (fOutterRadius-fInnerRadius)/fNRadialSteps;
0032   float dphi = TMath::TwoPi()/fNAzimuthalSteps;
0033   float dz = 2.0*fHalfLength/fNLongitudinalSteps;
0034   float dV = dr*dphi*dz;
0035
0036   float e0 = 8.854187817e+1; //[fC]/[V.cm]
0037   // looping over volume
0038   for(int br=0; br!=fNRadialSteps; ++br) {
0039     float r = fEr->GetXaxis()->GetBinCenter( br+1 ); //[cm]
0040     if(fDebug>0) printf("FieldMaps: %.2f -> 1\n",float(br)/fNRadialSteps);
0041     for(int bp=0; bp!=fNAzimuthalSteps; ++bp) {
0042       float phi = fEr->GetYaxis()->GetBinCenter( bp+1 ); //[rad]
0043       for(int bz=0; bz!=fNLongitudinalSteps; ++bz) {
0044     float z = fEr->GetZaxis()->GetBinCenter( bz+1 ); //[cm]
0045     //std::cout << " z :" << z << std::endl;
0046     if(fDebug>1) printf("FieldMaps: Integral( dvprime Q(rprime) G_E(rprime,[%f,%f,%f]) )",r,phi,z);
0047     float intEr = 0, intEp = 0, intEz = 0;
0048     if(fMirrorZ) if(z>0) continue;
0049     //std::cout << " z :" << z << std::endl;
0050     // looping over volume prime
0051     for(int bzprime=0; bzprime!=fNLongitudinalSteps; ++bzprime) {
0052       float zprime = fEr->GetZaxis()->GetBinCenter( bzprime+1 ); //[cm]
0053       //std::cout << " z*z' :" << z*zprime << std::endl;
0054       if(z*zprime < 0) continue; // only integrates same side along Z
0055       //std::cout << " z :"<< z << " z' :" << zprime << " z*z' :" << z*zprime << std::endl;
0056       for(int brprime=brprimeINI; brprime!=brprimeEND; ++brprime) {
0057         float rprime = fEr->GetXaxis()->GetBinCenter( brprime+1 ); //[cm]
0058         for(int bphiprime=0; bphiprime!=fNAzimuthalSteps; ++bphiprime) {
0059           float phiprime = fEr->GetYaxis()->GetBinCenter( bphiprime+1 );
0060           //float phiprime = fEr->GetYaxis()->GetBinCenter(1);
0061           float charge = ReadCharge(rprime,phiprime,zprime,dr,dphi,dz); // fC
0062           if(TMath::AreEqualAbs( charge, 0, prec )) continue;
0063           float GR = laplace_green_solution->Er(r/100,phi,TMath::Abs(z)/100, rprime/100, phiprime, TMath::Abs(zprime)/100); // GRad
0064           float GP = 0;//laplace_green_solution->Ephi(r/100,phi,TMath::Abs(z)/100, rprime/100, phi, TMath::Abs(zprime)/100); // GPhi
0065           float GZ = laplace_green_solution->Ez(r/100,phi,TMath::Abs(z)/100, rprime/100, phiprime, TMath::Abs(zprime)/100); // GPhi
0066           if(!TMath::AreEqualAbs( GR, 0, prec )) intEr += GR*charge; // [fC/m^2]
0067           //std::cout << " z :"<< z << " z' :" << zprime << " GR :" << GR << " q :" << charge <<" intEr :" << intEr<< std::endl;
0068           if(!TMath::AreEqualAbs( GP, 0, prec )) intEp += GP*charge; // [fC/m^2]
0069           if(!TMath::AreEqualAbs( GZ, 0, prec )) intEz += GZ*charge; // [fC/m^2]
0070         }//loop over phi'
0071       }// loop over r'
0072     }//loop over z'
0073     if(fDebug>1) printf(" => Er=%f ( similar for Ep=%f, Ez=%f )\n",intEr,intEp,intEz);
0074     float toVperCM = 0.155*1e-4/e0;
0075     fEr->Fill(r,phi,z,intEr*toVperCM); //[V/cm]
0076     fEp->Fill(r,phi,z,intEp*toVperCM); //[V/cm]
0077     fEz->Fill(r,phi,z,intEz*toVperCM); //[V/cm]
0078       } // loop over z
0079     } // loop over phi
0080   } //loop over r
0081   if(fDebug>0) printf("[DONE]\n");
0082   //------------
0083 }