StRoot: StHbtMaker/Infrastructure/StHbtCoulomb.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StHbtCoulomb.cxx,v 1.18 2006/08/15 21:41:41 jeromel Exp $
00004  *
00005  * Author: Randy Wells, Ohio State, rcwells@mps.ohio-state.edu
00006  ***************************************************************************
00007  *
00008  * Description: part of STAR HBT Framework: StHbtMaker package
00009  *    This is a Coulomb correction class which
00010  *  1. Reads in the dat from a file
00011  *  2. Performs a linear interpolation in R and creates any array of interpolations
00012  *  3. Interpolates in eta and returns the Coulomb correction to user
00013  *
00014  ***************************************************************************
00015  *
00016  * $Log: StHbtCoulomb.cxx,v $
00017  * Revision 1.18  2006/08/15 21:41:41  jeromel
00018  * Fix rhic -> rhic.bnl.gov
00019  *
00020  * Revision 1.17  2003/09/02 17:58:32  perev
00021  * gcc 3.2 updates + WarnOff
00022  *
00023  * Revision 1.16  2003/02/04 21:10:31  magestro
00024  * Cleaned up a couple functions
00025  *
00026  * Revision 1.15  2003/01/31 19:44:00  magestro
00027  * Cleared up simple compiler warnings on i386_linux24
00028  *
00029  * Revision 1.14  2000/10/26 19:48:54  rcwells
00030  * Added functionality for Coulomb correction of <qInv> in 3D correltions
00031  *
00032  * Revision 1.13  2000/07/16 21:38:22  laue
00033  * StHbtCoulomb.cxx StHbtSectoredAnalysis.cxx : updated for standalone version
00034  * StHbtV0.cc StHbtV0.hh : some cast to prevent compiling warnings
00035  * StHbtParticle.cc StHbtParticle.hh : pointers mTrack,mV0 initialized to 0
00036  * StHbtIOBinary.cc : some printouts in #ifdef STHBTDEBUG
00037  * StHbtEvent.cc : B-Field set to 0.25Tesla, we have to think about a better
00038  *                 solution
00039  *
00040  * Revision 1.12  2000/05/31 20:12:53  rcwells
00041  * Modified StHbtCoulomb for Id and Log entries
00042  *
00043  *
00044  **************************************************************************/
00045 
00046 #include "StHbtCoulomb.h"
00047 #include "Stiostream.h"
00048 #include <stdio.h>
00049 #include <cassert>
00050 #include "PhysicalConstants.h"
00051 
00052 #ifdef __ROOT__
00053 ClassImp(StHbtCoulomb)
00054 #endif
00055 
00056 StHbtCoulomb::StHbtCoulomb() {
00057   mFile = "/afs/rhic.bnl.gov/star/hbt/coul/StHbtCorrectionFiles/correctionpp.dat";
00058   if (!mFile) {
00059     cout << " No file, dummy!" << endl;
00060     assert(0);
00061   }
00062   mRadius = -1.0;
00063   mZ1Z2 = 1.0; // Default has been set up to be same sign charges
00064   cout << "You have 1 default Coulomb correction!" << endl;
00065 }
00066 
00067 StHbtCoulomb::StHbtCoulomb(const char* readFile, const double& radius, const double& charge) {
00068   mFile = readFile;
00069   mRadius = radius;
00070   CreateLookupTable(mRadius);
00071   mZ1Z2 = charge;
00072   cout << "You have 1 Coulomb correction!" << endl;
00073 }
00074 
00075 StHbtCoulomb::~StHbtCoulomb() {
00076 
00077 }
00078 
00079 void StHbtCoulomb::SetRadius(const double& radius) {
00080   cout << " StHbtCoulomb::setRadius() " << endl;
00081   mRadius = radius;
00082   CreateLookupTable(mRadius);
00083 }
00084 
00085 double StHbtCoulomb::GetRadius() {
00086   return (mRadius);
00087 }
00088 
00089 void StHbtCoulomb::SetFile(const char* readFile) {
00090   cout << " StHbtCoulomb::SetFile() " << endl;
00091   mFile = readFile;
00092   // Create new lookup table since file has changed
00093   if (mRadius>0.0) {
00094     CreateLookupTable(mRadius);
00095   }
00096 }
00097 
00098 void StHbtCoulomb::SetChargeProduct(const double& charge) {
00099   cout << " StHbtCoulomb::SetChargeProduct() " << endl;
00100   if ( mZ1Z2!=charge ) { 
00101     mZ1Z2 = charge;
00102     if ( mZ1Z2>0 ) {
00103       mFile = "/afs/rhic.bnl.gov/star/hbt/coul/StHbtCorrectionFiles/correctionpp.dat";
00104     }
00105     else {
00106       mFile = "/afs/rhic.bnl.gov/star/hbt/coul/StHbtCorrectionFiles/correctionpm.dat";
00107     }
00108     CreateLookupTable(mRadius);
00109   }
00110 }
00111 
00112 void StHbtCoulomb::CreateLookupTable(const double& radius) {
00113   cout << " StHbtCoulomb::CreateLookupTable() " << endl;
00114   // Read radii from mFile
00115   // Create array(pair) of linear interpolation between radii
00116 
00117   if (radius<0.0) {
00118     cout << " StHbtCoulomb::CreateLookupTable -> NEGATIVE RADIUS " << endl;
00119     cout << "  call StHbtCoulomb::SetRadius(r) with positive r " << endl;
00120     cerr << " StHbtCoulomb::CreateLookupTable -> NEGATIVE RADIUS " << endl;
00121     cerr << "  call StHbtCoulomb::SetRadius(r) with positive r " << endl;
00122     assert(0);
00123   }
00124   ifstream mystream(mFile);
00125   if (!mystream) {
00126     cout << "Could not open file" << endl;
00127     assert(0);
00128   }
00129   else {
00130     cout << "Input correction file opened" << endl;
00131   }
00132 
00133   static char tempstring[2001];
00134   static float radii[2000];
00135   static int NRadii = 0;
00136   NRadii = 0;
00137   if (!mystream.getline(tempstring,2000)) {
00138     cout << "Could not read radii from file" << endl;
00139     assert(0);
00140   }
00141   for (unsigned int ii=0; ii<strlen(tempstring); ii++) {
00142     while (tempstring[ii]==' ') ii++;
00143     sscanf(&tempstring[ii++],"%f",&radii[++NRadii]);
00144     while ( tempstring[ii]!=' ' && (ii)<strlen(tempstring) )ii++;
00145   }
00146   cout << " Read " << NRadii << " radii from file" << endl;
00147 
00148   static double LowRadius = -1.0;
00149   static double HighRadius = -1.0;
00150   static int LowIndex = 0;
00151   LowRadius = -1.0;
00152   HighRadius = -1.0;
00153   LowIndex = 0;
00154   for(int iii=1; iii<=NRadii-1; iii++) { // Loop to one less than #radii
00155     if ( radius >= radii[iii] && radius <= radii[iii+1] ) {
00156       LowRadius = radii[iii];
00157       HighRadius = radii[iii+1];
00158       LowIndex = iii;
00159     }
00160   }
00161   if ( (LowRadius < 0.0) || (HighRadius < 0.0) ) {
00162     cout << "StHbtCoulomb::CreateLookupTable --> Problem interpolating radius" << endl;
00163     cout << "  Check range of radii in lookup file...." << endl;
00164     cerr << "StHbtCoulomb::CreateLookupTable --> Problem interpolating radius" << endl;
00165     cerr << "  Check range of radii in lookup file...." << endl;
00166     assert(0);
00167   }
00168 
00169   static double corr[100];           // array of corrections ... must be > NRadii
00170   mNLines = 0;
00171   static double tempEta = 0;
00172   tempEta = 0;
00173   while (mystream >> tempEta) {
00174     for (int i=1; i<=NRadii; i++) {
00175       mystream >> corr[i];
00176     }
00177     static double LowCoulomb = 0;
00178     static double HighCoulomb = 0;
00179     static double nCorr = 0;
00180     LowCoulomb = corr[LowIndex];
00181     HighCoulomb = corr[LowIndex+1];
00182     nCorr = ( (radius-LowRadius)*HighCoulomb+(HighRadius-radius)*LowCoulomb )/(HighRadius-LowRadius);
00183       mEta[mNLines] = tempEta;     // Eta
00184       mCoulomb[mNLines] = nCorr;   // Interpolated Coulomb correction for radius
00185       mNLines++;
00186   }
00187   mystream.close();
00188   cout << "Lookup Table is created with " << mNLines << " points" << endl;
00189 }
00190 
00191 double StHbtCoulomb::CoulombCorrect(const double& eta) {
00192   // Interpolates in eta
00193   if (mRadius < 0.0) {
00194     cout << "StHbtCoulomb::CoulombCorrect(eta) --> Trying to correct for negative radius!" << endl;
00195     cerr << "StHbtCoulomb::CoulombCorrect(eta) --> Trying to correct for negative radius!" << endl;
00196     assert(0);
00197   }
00198   static int middle=0;
00199   middle=int( (mNLines-1)/2 );
00200   if (eta*mEta[middle]<0.0) {
00201     cout << "StHbtCoulomb::CoulombCorrect(eta) --> eta: " << eta << " has wrong sign for data file! " << endl;
00202     cerr << "StHbtCoulomb::CoulombCorrect(eta) --> eta: " << eta << " has wrong sign for data file! " << endl;
00203     assert(0);
00204   }
00205 
00206   static double Corr = 0;
00207   Corr = -1.0;
00208   
00209   if ( (eta>mEta[0]) && (mEta[0]>0.0) ) {
00210     Corr = mCoulomb[0];
00211     return (Corr);
00212   }
00213   if ( (eta<mEta[mNLines-1]) && (mEta[mNLines-1]<0.0) ) {
00214     Corr = mCoulomb[mNLines-1];
00215     return (Corr);
00216   }
00217   // This is a binary search for the bracketing pair of data points
00218   static int high = 0;
00219   static int low = 0;
00220   static int width = 0;
00221   high = mNLines-1;
00222   low = 0;
00223   width = high-low;
00224   middle = int(width/2.0); // Was instantiated above
00225   while (middle > 0) {
00226     if (mEta[low+middle] < eta) {
00227       // eta is in the 1st half
00228       high-=middle;
00229       width = high-low;
00230       middle = int(width/2.0);
00231     }
00232     else {
00233       // eta is in the 2nd half
00234       low+=middle;
00235       width = high-low;
00236       middle = int(width/2.0);
00237     }
00238   }
00239   // Make sure we found the right one
00240   if ( (mEta[low] >= eta) && (eta >= mEta[low+1]) ) {
00241     static double LowEta = 0;
00242     static double HighEta = 0;    
00243     static double LowCoulomb = 0;
00244     static double HighCoulomb = 0;
00245     LowEta = mEta[low];
00246     HighEta = mEta[low+1];    
00247     LowCoulomb = mCoulomb[low];
00248     HighCoulomb = mCoulomb[low+1];
00249     //      cout << LowEta << " *** Eta *** " << HighEta << endl;
00250     //      cout << LowCoulomb << " *** Coulomb *** " << HighCoulomb << endl;
00251     Corr = ( (eta-LowEta)*HighCoulomb+(HighEta-eta)*LowCoulomb )/(HighEta-LowEta);
00252   }
00253   if (Corr<0.0) {
00254     cout << "StHbtCoulomb::CoulombCorrect(eta) --> No correction" << endl;
00255     cout << "  Check range of eta in file: Input eta  " << eta << endl;
00256     cerr << "StHbtCoulomb::CoulombCorrect(eta) --> No correction" << endl;
00257     cerr << "  Check range of eta in file: Input eta  " << eta << endl;
00258     assert(0);
00259   } 
00260   return (Corr);
00261 
00262 }
00263 
00264 double StHbtCoulomb::CoulombCorrect(const double& eta,
00265                                 const double& radius) {
00266   // Checks radii ... input radius and mRadius
00267   // Calls createLookupTable if neccessary
00268   // Interpolate(linear) between etas in the created lookup table
00269 
00270   if (radius < 0.0) {
00271     if (mRadius < 0.0) {
00272       // Both radii are negative
00273       cout << "StHbtCoulomb::CoulombCorrect(eta,r) --> input and member radii are negative!" << endl;
00274       cerr << "StHbtCoulomb::CoulombCorrect(eta,r) --> input and member radii are negative!" << endl;
00275       assert(0);
00276     }
00277   }
00278   else {
00279     // radius > 0.0
00280     if (radius == mRadius) {
00281       // Both radii are positive and equal
00282       //      cout << "Radii are the same!!!" << endl;
00283     }
00284     else {
00285       // Both radii are positive but not equal
00286       mRadius = radius;
00287       CreateLookupTable(mRadius);
00288     }
00289   }
00290 
00291   // Interpolate in eta
00292   return ( CoulombCorrect(eta) );
00293 }
00294 
00295 double StHbtCoulomb::CoulombCorrect(const StHbtPair* pair) {
00296   return ( CoulombCorrect( Eta(pair) ) );;
00297 }
00298 
00299 double StHbtCoulomb::CoulombCorrect(const StHbtPair* pair, const double& radius) {
00300   return ( CoulombCorrect( Eta(pair),radius ) );
00301 }
00302 
00303 double StHbtCoulomb::Eta(const StHbtPair* pair) {
00304   static double px1,py1,pz1,px2,py2,pz2;
00305   static double px1new,py1new,pz1new;
00306   static double px2new,py2new,pz2new;
00307   static double vx1cms,vy1cms,vz1cms;
00308   static double vx2cms,vy2cms,vz2cms;
00309   static double VcmsX,VcmsY,VcmsZ;
00310   static double dv = 0.0;
00311   static double e1,e2,e1new,e2new;
00312   static double psi,theta;
00313   static double beta,gamma;
00314   static double VcmsXnew;
00315 
00316   px1 = pair->track1()->FourMomentum().px();
00317   py1 = pair->track1()->FourMomentum().py();
00318   pz1 = pair->track1()->FourMomentum().pz();
00319   e1 = pair->track1()->FourMomentum().e();
00320   px2 = pair->track2()->FourMomentum().px();
00321   py2 = pair->track2()->FourMomentum().py();
00322   pz2 = pair->track2()->FourMomentum().pz();
00323   e2 = pair->track2()->FourMomentum().e();
00324   
00325   VcmsX = ( px1+px2 )/( e1+e2 );
00326   VcmsY = ( py1+py2 )/( e1+e2 );
00327   VcmsZ = ( pz1+pz2 )/( e1+e2 );
00328   // Rotate Vcms to x-direction
00329   psi = atan(VcmsY/VcmsX);
00330   VcmsXnew = VcmsX*cos(psi)+VcmsY*sin(psi);
00331   VcmsX = VcmsXnew;
00332   theta = atan(VcmsZ/VcmsX);
00333   VcmsXnew = VcmsX*cos(theta)+VcmsZ*sin(theta);
00334   VcmsX = VcmsXnew;
00335   // Gamma and Beta
00336   beta = VcmsX;
00337   gamma = 1.0/::sqrt( 1.0-beta*beta );
00338 
00339   // Rotate p1 and p2 to new frame
00340   px1new = px1*cos(psi)+py1*sin(psi);
00341   py1new = -px1*sin(psi)+py1*cos(psi);
00342   px1 = px1new;
00343   px1new = px1*cos(theta)+pz1*sin(theta);
00344   pz1new = -px1*sin(theta)+pz1*cos(theta);
00345   px1 = px1new;
00346   py1 = py1new;
00347   pz1 = pz1new;
00348 
00349   px2new = px2*cos(psi)+py2*sin(psi);
00350   py2new = -px2*sin(psi)+py2*cos(psi);
00351   px2 = px2new;
00352   px2new = px2*cos(theta)+pz2*sin(theta);
00353   pz2new = -px2*sin(theta)+pz2*cos(theta);
00354   px2 = px2new;
00355   py2 = py2new;
00356   pz2 = pz2new;
00357 
00358   // Lorentz transform the x component and energy
00359   e1new = gamma*e1 - gamma*beta*px1;
00360   px1new = -gamma*beta*e1 + gamma*px1;
00361   e2new = gamma*e2 - gamma*beta*px2;
00362   px2new = -gamma*beta*e2 + gamma*px2;
00363   px1 = px1new;
00364   px2 = px2new;
00365 
00366   // New velocities
00367   vx1cms = px1/e1new;
00368   vy1cms = py1/e1new;
00369   vz1cms = pz1/e1new;
00370   vx2cms = px2/e2new;
00371   vy2cms = py2/e2new;
00372   vz2cms = pz2/e2new;
00373 
00374   // Velocity difference in CMS frame
00375   dv = ::sqrt( (vx1cms-vx2cms)*(vx1cms-vx2cms) +
00376              (vy1cms-vy2cms)*(vy1cms-vy2cms) +
00377              (vz1cms-vz2cms)*(vz1cms-vz2cms) );
00378 
00379   return ( mZ1Z2*fine_structure_const/(dv) );
00380 }
00381 
00382 StHbt1DHisto* StHbtCoulomb::CorrectionHistogram(const double& mass1, const double& mass2, const int& nBins, 
00383                                                 const double& low, const double& high) {
00384   if ( mass1!=mass2 ) {
00385     cout << "Masses not equal ... try again.  No histogram created." << endl;
00386     assert(0);
00387   }
00388   StHbt1DHisto* correction = new StHbt1DHisto("correction","Coulomb correction",nBins,low,high);
00389   const double reducedMass = mass1*mass2/(mass1+mass2);
00390   double qInv = low;
00391   //double dQinv = (high-low)/( (double)nBins );
00392   double eta;
00393   for (int ii=1; ii<=nBins; ii++) 
00394     {
00395       qInv = correction->GetBinCenter(ii);
00396       eta = 2.0*mZ1Z2*reducedMass*fine_structure_const/( qInv );
00397       CoulombCorrect( eta );
00398       correction->Fill( qInv, CoulombCorrect(eta,mRadius) );
00399     }
00400 
00401   return (correction);
00402 }
00403 
00404 #ifdef __ROOT__
00405 StHbt1DHisto* StHbtCoulomb::CorrectionHistogram(const StHbt1DHisto* histo, const double mass) {
00406 
00407   StHbt1DHisto* correction = (StHbt1DHisto*) ((StHbt1DHisto*)histo)->Clone();
00408   correction->Reset();
00409   correction->SetDirectory(0);
00410   int    nBins = correction->GetXaxis()->GetNbins();
00411   const double reducedMass = 0.5*mass;
00412   double qInv;
00413   double eta;
00414   for (int ii=1; ii<=nBins; ii++) 
00415     {
00416       qInv = correction->GetBinCenter(ii);
00417       eta = 2.0*mZ1Z2*reducedMass*fine_structure_const/( qInv );
00418       correction->Fill( qInv, CoulombCorrect(eta,mRadius) );
00419     }
00420 
00421   return (correction);
00422 }
00423 
00424 StHbt3DHisto* StHbtCoulomb::CorrectionHistogram(const StHbt3DHisto* histo, const double mass) {
00425 
00426   StHbt3DHisto* correction = (StHbt3DHisto*) ((StHbt3DHisto*)histo)->Clone();
00427   correction->Reset();
00428   correction->SetDirectory(0);
00429   int    nBinsX = correction->GetXaxis()->GetNbins();
00430   int    nBinsY = correction->GetYaxis()->GetNbins();
00431   int    nBinsZ = correction->GetZaxis()->GetNbins();
00432   const double reducedMass = 0.5*mass;
00433   double eta;
00434   double qInv;
00435   int binNumber;
00436   for (int ii=1; ii<=nBinsX; ii++) { 
00437     for (int iii=1; iii<=nBinsY; iii++) {
00438       for (int iv=1; iv<=nBinsZ; iv++) {
00439         binNumber = histo->GetBin(ii,iii,iv);
00440         qInv = histo->GetBinContent(binNumber);
00441         eta = 2.0*mZ1Z2*reducedMass*fine_structure_const/( qInv );
00442         correction->SetBinContent(binNumber, CoulombCorrect(eta,mRadius) );
00443       }
00444     }
00445   }
00446   return (correction);
00447 }
00448 #endif
00449 
00450 double StHbtCoulomb::CoulombCorrect(const double& mass, const double& charge,
00451                                     const double& radius, const double& qInv) {
00452   mRadius = radius;
00453   mZ1Z2 = charge;
00454   const double reducedMass = 0.5*mass; // must be same mass particles
00455   double eta = 2.0*mZ1Z2*reducedMass*fine_structure_const/( qInv );
00456   return ( CoulombCorrect(eta,mRadius) );
00457 }