TPC/SCDCorrections/Langevin.cxx

0001 #include "Langevin.h"
0002 #include <string>
0003
0004 #include "TH3F.h"
0005 #include "TFile.h"
0006 #include "TMath.h"
0007
0008 //=====================
0009 Langevin::Langevin() {
0010   fDebug = 0;
0011   fEr = NULL;
0012   fEp = NULL;
0013   fEz = NULL;
0014   fDeltaR = NULL;
0015   fRDeltaPHI = NULL;
0016   fInnerRadius=1.;
0017   fOutterRadius=2.;
0018   fHalfLength=1.;
0019   fNRadialSteps=1;
0020   fNAzimuthalSteps=1;
0021   fNLongitudinalSteps=1;
0022   fFileNameRoot="rho";
0023 }
0024 //=====================
0025 Langevin::~Langevin() {
0026 }
0027 //=====================
0028 void Langevin::Make() {
0029   InitMaps();
0030
0031   //------------
0032   //STEER
0033   if(fEr==NULL) {
0034     printf("Input fields not found. A B O R T I N G\n");
0035     return;
0036   }
0037   if(fDebug>0) printf("Langevin solutions are being computed... \n");
0038
0039   const float prec=1e-8;
0040   const float kE0 = 400; //V/cm
0041   //const float kB0 = 5; //kGauss ALICE / STAR
0042   const float kB0 = 15; //kGauss PHENIX
0043   //Langevin gas parameters
0044   // P9 gas
0045   //const float kT1 = 1.34; //NIM Phys Res A235,296 (1985)
0046   //const float kT2 = 1.11; //NIM Phys Res A235,296 (1985)
0047   // P10 gas
0048   //const float kT1 = 1.36; //measured by STAR
0049   //const float kT2 = 1.11; //measured by STAR
0050   // 85.7%Ne / 9.5%CO2 / 4.5%N2
0051   const float kT1 = 1.0; // measured by ALICE
0052   const float kT2 = 1.0; // measured by ALICE
0053   const float kV0 = 4.0; // [cm/us] @ E=400V/cm (need fine tunning)
0054   float wt = -10*kB0*kV0/kE0; //[kGauss] * [cm/us] / [V/cm] //0.32 != -10*5*4/400
0055   float c0 = 1.0/(1.0+kT2*kT2*wt*wt);
0056   float c1 = kT1*wt/(1.0+kT1*kT1*wt*wt);
0057   float c2 = kT2*kT2*wt*wt*c0;
0058   if(fDebug>1) printf("wt = %f \n",wt);
0059   if(fDebug>1) printf("c0 = %f \n",c0);
0060   if(fDebug>1) printf("c1 = %f \n",c1);
0061   if(fDebug>1) printf("c2 = %f \n",c2);
0062   float dz = fEz->GetZaxis()->GetBinWidth(1);
0063   for(int r=0; r!=fNRadialSteps; ++r)
0064     for(int p=0; p!=fNAzimuthalSteps; ++p)
0065       for(int z=0; z!=fNLongitudinalSteps; ++z) {
0066     float intDR = 0, intRDPHI = 0;
0067     // integrating over drift ditance
0068     int minZ = 0;
0069     int maxZ = z+1;
0070     if( fEz->GetZaxis()->GetBinCenter(z+1)>0 ) {
0071       minZ=z;
0072       maxZ=fNLongitudinalSteps;
0073     }
0074     for(int zprime=minZ; zprime!=maxZ; ++zprime) {
0075       float dR, RdPHI;
0076       float dEz = fEz->GetBinContent( r+1, p+1, z+1 ) + kE0;
0077       dR = c0*fEr->GetBinContent( r+1, p+1, zprime+1 )/dEz;
0078       dR += c1*fEp->GetBinContent( r+1, p+1, zprime+1 )/dEz;
0079       RdPHI = -c1*fEr->GetBinContent( r+1, p+1, zprime+1 )/dEz;
0080       RdPHI += c0*fEp->GetBinContent( r+1, p+1, zprime+1 )/dEz;
0081       if(!TMath::AreEqualAbs( dR, 0, prec )) intDR += dR*dz;
0082       if(!TMath::AreEqualAbs( RdPHI, 0, prec )) intRDPHI += RdPHI*dz;
0083     }
0084     fDeltaR->SetBinContent( r+1, p+1, z+1, intDR );
0085     fRDeltaPHI->SetBinContent( r+1, p+1, z+1, intRDPHI );
0086     if(fDebug>2) printf("@{Ir,Ip,Iz}={%d,%d,%d}, deltaR=%f\n",r,p,z,intDR);
0087     if(fDebug>2) printf("@{Ir,Ip,Iz}={%d,%d,%d}, RdeltaPHI=%f\n",r,p,z,intRDPHI);
0088       }
0089   if(fDebug>0) printf("[DONE]\n");
0090   //------------
0091
0092   SaveMaps();
0093 }
0094 //=====================
0095 void Langevin::ReadFile() {
0096   const char *inputfile = Form("%s_1.root",fFileNameRoot.data());
0097   if(fDebug>0) printf("Langevin is reading the input file %s... ",inputfile);
0098   TFile *ifile = new TFile(inputfile);
0099   fEr = (TH3F*) ifile->Get("Er");
0100   fEp = (TH3F*) ifile->Get("Ep");
0101   fEz = (TH3F*) ifile->Get("Ez");
0102   if(fDebug>0) printf("[DONE]\n");
0103 }
0104 //=====================
0105 void Langevin::InitMaps() {
0106   ReadFile();
0107   if(fDebug>0) printf("Langevin is initializing distortion maps... ");
0108   fDeltaR = new TH3F("mapDeltaR","#delta_{R} [cm];Radial [cm];Azimuthal [rad];Longitudinal [cm]",
0109              fNRadialSteps,fInnerRadius,fOutterRadius,
0110              fNAzimuthalSteps,0,TMath::TwoPi(),
0111              fNLongitudinalSteps,-fHalfLength,+fHalfLength);
0112   fRDeltaPHI = new TH3F("mapRDeltaPHI","R#delta_{#varphi} [cm];Radial [cm];Azimuthal [rad];Longitudinal [cm]",
0113              fNRadialSteps,fInnerRadius,fOutterRadius,
0114              fNAzimuthalSteps,0,TMath::TwoPi(),
0115              fNLongitudinalSteps,-fHalfLength,+fHalfLength);
0116   if(fDebug>0) printf("[DONE]\n");
0117 }
0118 //=====================
0119 void Langevin::SaveMaps() {
0120   const char *outputfile= Form("%s_2.root",fFileNameRoot.data());
0121   if(fDebug>0) printf("Langevin saving distortion maps... ");
0122   TFile *ofile = new TFile(outputfile,"RECREATE");
0123   ofile->WriteObject(fDeltaR,"mapDeltaR");
0124   ofile->WriteObject(fRDeltaPHI,"mapRDeltaPHI");
0125   ofile->Close();
0126   if(fDebug>0) printf("[DONE]\n");
0127 }