PWG2/FEMTOSCOPY/AliFemto/Infrastructure/AliFemtoCoulomb.cxx

   1 /***************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Author: Randy Wells, Ohio State, rcwells@mps.ohio-state.edu
   6  ***************************************************************************
   7  *
   8  * Description: part of STAR HBT Framework: AliFemtoMaker package
   9  *    This is a Coulomb correction class which
  10  *  1. Reads in the dat from a file
  11  *  2. Performs a linear interpolation in R and creates any array of interpolations
  12  *  3. Interpolates in eta and returns the Coulomb correction to user
  13  *
  14  ***************************************************************************
  15  *
  16  * $Log$
  17  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  18  * First version on CVS
  19  *
  20  * Revision 1.17  2003/09/02 17:58:32  perev
  21  * gcc 3.2 updates + WarnOff
  22  *
  23  * Revision 1.16  2003/02/04 21:10:31  magestro
  24  * Cleaned up a couple functions
  25  *
  26  * Revision 1.15  2003/01/31 19:44:00  magestro
  27  * Cleared up simple compiler warnings on i386_linux24
  28  *
  29  * Revision 1.14  2000/10/26 19:48:54  rcwells
  30  * Added functionality for Coulomb correction of <qInv> in 3D correltions
  31  *
  32  * Revision 1.13  2000/07/16 21:38:22  laue
  33  * AliFemtoCoulomb.cxx AliFemtoSectoredAnalysis.cxx : updated for standalone version
  34  * AliFemtoV0.cc AliFemtoV0.h : some cast to prevent compiling warnings
  35  * AliFemtoParticle.cc AliFemtoParticle.h : pointers mTrack,mV0 initialized to 0
  36  * AliFemtoIOBinary.cc : some printouts in #ifdef STHBTDEBUG
  37  * AliFemtoEvent.cc : B-Field set to 0.25Tesla, we have to think about a better
  38  *                 solution
  39  *
  40  * Revision 1.12  2000/05/31 20:12:53  rcwells
  41  * Modified AliFemtoCoulomb for Id and Log entries
  42  *
  43  *
  44  **************************************************************************/
  45
  46 #include "AliFemtoCoulomb.h"
  47 //#include "Stiostream.h"
  48 #include <stdio.h>
  49 #include <cassert>
  50 #include "PhysicalConstants.h"
  51
  52 #ifdef __ROOT__
  53 ClassImp(AliFemtoCoulomb)
  54 #endif
  55
  56 AliFemtoCoulomb::AliFemtoCoulomb() {
  57   fFile = "/afs/rhic/star/hbt/coul/AliFemtoCorrectionFiles/correctionpp.dat";
  58   if (!fFile) {
  59     cout << " No file, dummy!" << endl;
  60     assert(0);
  61   }
  62   fRadius = -1.0;
  63   fZ1Z2 = 1.0; // Default has been set up to be same sign charges
  64   cout << "You have 1 default Coulomb correction!" << endl;
  65 }
  66
  67 AliFemtoCoulomb::AliFemtoCoulomb(const char* readFile, const double& radius, const double& charge) {
  68   fFile = readFile;
  69   fRadius = radius;
  70   CreateLookupTable(fRadius);
  71   fZ1Z2 = charge;
  72   cout << "You have 1 Coulomb correction!" << endl;
  73 }
  74
  75 AliFemtoCoulomb::~AliFemtoCoulomb() {
  76
  77 }
  78
  79 void AliFemtoCoulomb::SetRadius(const double& radius) {
  80   cout << " AliFemtoCoulomb::setRadius() " << endl;
  81   fRadius = radius;
  82   CreateLookupTable(fRadius);
  83 }
  84
  85 double AliFemtoCoulomb::GetRadius() {
  86   return (fRadius);
  87 }
  88
  89 void AliFemtoCoulomb::SetFile(const char* readFile) {
  90   cout << " AliFemtoCoulomb::SetFile() " << endl;
  91   fFile = readFile;
  92   // Create new lookup table since file has changed
  93   if (fRadius>0.0) {
  94     CreateLookupTable(fRadius);
  95   }
  96 }
  97
  98 void AliFemtoCoulomb::SetChargeProduct(const double& charge) {
  99   cout << " AliFemtoCoulomb::SetChargeProduct() " << endl;
 100   if ( fZ1Z2!=charge ) {
 101     fZ1Z2 = charge;
 102     if ( fZ1Z2>0 ) {
 103       fFile = "/afs/rhic/star/hbt/coul/AliFemtoCorrectionFiles/correctionpp.dat";
 104     }
 105     else {
 106       fFile = "/afs/rhic/star/hbt/coul/AliFemtoCorrectionFiles/correctionpm.dat";
 107     }
 108     CreateLookupTable(fRadius);
 109   }
 110 }
 111
 112 void AliFemtoCoulomb::CreateLookupTable(const double& radius) {
 113   cout << " AliFemtoCoulomb::CreateLookupTable() " << endl;
 114   // Read radii from fFile
 115   // Create array(pair) of linear interpolation between radii
 116
 117   if (radius<0.0) {
 118     cout << " AliFemtoCoulomb::CreateLookupTable -> NEGATIVE RADIUS " << endl;
 119     cout << "  call AliFemtoCoulomb::SetRadius(r) with positive r " << endl;
 120     cerr << " AliFemtoCoulomb::CreateLookupTable -> NEGATIVE RADIUS " << endl;
 121     cerr << "  call AliFemtoCoulomb::SetRadius(r) with positive r " << endl;
 122     assert(0);
 123   }
 124   ifstream mystream(fFile);
 125   if (!mystream) {
 126     cout << "Could not open file" << endl;
 127     assert(0);
 128   }
 129   else {
 130     cout << "Input correction file opened" << endl;
 131   }
 132
 133   static char tempstring[2001];
 134   static float radii[2000];
 135   static int NRadii = 0;
 136   NRadii = 0;
 137   if (!mystream.getline(tempstring,2000)) {
 138     cout << "Could not read radii from file" << endl;
 139     assert(0);
 140   }
 141   for (unsigned int ii=0; ii<strlen(tempstring); ii++) {
 142     while (tempstring[ii]==' ') ii++;
 143     sscanf(&tempstring[ii++],"%f",&radii[++NRadii]);
 144     while ( tempstring[ii]!=' ' && (ii)<strlen(tempstring) )ii++;
 145   }
 146   cout << " Read " << NRadii << " radii from file" << endl;
 147
 148   static double LowRadius = -1.0;
 149   static double HighRadius = -1.0;
 150   static int LowIndex = 0;
 151   LowRadius = -1.0;
 152   HighRadius = -1.0;
 153   LowIndex = 0;
 154   for(int iii=1; iii<=NRadii-1; iii++) { // Loop to one less than #radii
 155     if ( radius >= radii[iii] && radius <= radii[iii+1] ) {
 156       LowRadius = radii[iii];
 157       HighRadius = radii[iii+1];
 158       LowIndex = iii;
 159     }
 160   }
 161   if ( (LowRadius < 0.0) || (HighRadius < 0.0) ) {
 162     cout << "AliFemtoCoulomb::CreateLookupTable --> Problem interpolating radius" << endl;
 163     cout << "  Check range of radii in lookup file...." << endl;
 164     cerr << "AliFemtoCoulomb::CreateLookupTable --> Problem interpolating radius" << endl;
 165     cerr << "  Check range of radii in lookup file...." << endl;
 166     assert(0);
 167   }
 168
 169   static double corr[100];           // array of corrections ... must be > NRadii
 170   fNLines = 0;
 171   static double tempEta = 0;
 172   tempEta = 0;
 173   while (mystream >> tempEta) {
 174     for (int i=1; i<=NRadii; i++) {
 175       mystream >> corr[i];
 176     }
 177     static double LowCoulomb = 0;
 178     static double HighCoulomb = 0;
 179     static double nCorr = 0;
 180     LowCoulomb = corr[LowIndex];
 181     HighCoulomb = corr[LowIndex+1];
 182     nCorr = ( (radius-LowRadius)*HighCoulomb+(HighRadius-radius)*LowCoulomb )/(HighRadius-LowRadius);
 183       fEta[fNLines] = tempEta;     // Eta
 184       fCoulomb[fNLines] = nCorr;   // Interpolated Coulomb correction for radius
 185       fNLines++;
 186   }
 187   mystream.close();
 188   cout << "Lookup Table is created with " << fNLines << " points" << endl;
 189 }
 190
 191 double AliFemtoCoulomb::CoulombCorrect(const double& eta) {
 192   // Interpolates in eta
 193   if (fRadius < 0.0) {
 194     cout << "AliFemtoCoulomb::CoulombCorrect(eta) --> Trying to correct for negative radius!" << endl;
 195     cerr << "AliFemtoCoulomb::CoulombCorrect(eta) --> Trying to correct for negative radius!" << endl;
 196     assert(0);
 197   }
 198   static int middle=0;
 199   middle=int( (fNLines-1)/2 );
 200   if (eta*fEta[middle]<0.0) {
 201     cout << "AliFemtoCoulomb::CoulombCorrect(eta) --> eta: " << eta << " has wrong sign for data file! " << endl;
 202     cerr << "AliFemtoCoulomb::CoulombCorrect(eta) --> eta: " << eta << " has wrong sign for data file! " << endl;
 203     assert(0);
 204   }
 205
 206   static double Corr = 0;
 207   Corr = -1.0;
 208
 209   if ( (eta>fEta[0]) && (fEta[0]>0.0) ) {
 210     Corr = fCoulomb[0];
 211     return (Corr);
 212   }
 213   if ( (eta<fEta[fNLines-1]) && (fEta[fNLines-1]<0.0) ) {
 214     Corr = fCoulomb[fNLines-1];
 215     return (Corr);
 216   }
 217   // This is a binary search for the bracketing pair of data points
 218   static int high = 0;
 219   static int low = 0;
 220   static int width = 0;
 221   high = fNLines-1;
 222   low = 0;
 223   width = high-low;
 224   middle = int(width/2.0); // Was instantiated above
 225   while (middle > 0) {
 226     if (fEta[low+middle] < eta) {
 227       // eta is in the 1st half
 228       high-=middle;
 229       width = high-low;
 230       middle = int(width/2.0);
 231     }
 232     else {
 233       // eta is in the 2nd half
 234       low+=middle;
 235       width = high-low;
 236       middle = int(width/2.0);
 237     }
 238   }
 239   // Make sure we found the right one
 240   if ( (fEta[low] >= eta) && (eta >= fEta[low+1]) ) {
 241     static double LowEta = 0;
 242     static double HighEta = 0;
 243     static double LowCoulomb = 0;
 244     static double HighCoulomb = 0;
 245     LowEta = fEta[low];
 246     HighEta = fEta[low+1];
 247     LowCoulomb = fCoulomb[low];
 248     HighCoulomb = fCoulomb[low+1];
 249     //      cout << LowEta << " *** Eta *** " << HighEta << endl;
 250     //      cout << LowCoulomb << " *** Coulomb *** " << HighCoulomb << endl;
 251     Corr = ( (eta-LowEta)*HighCoulomb+(HighEta-eta)*LowCoulomb )/(HighEta-LowEta);
 252   }
 253   if (Corr<0.0) {
 254     cout << "AliFemtoCoulomb::CoulombCorrect(eta) --> No correction" << endl;
 255     cout << "  Check range of eta in file: Input eta  " << eta << endl;
 256     cerr << "AliFemtoCoulomb::CoulombCorrect(eta) --> No correction" << endl;
 257     cerr << "  Check range of eta in file: Input eta  " << eta << endl;
 258     assert(0);
 259   }
 260   return (Corr);
 261
 262 }
 263
 264 double AliFemtoCoulomb::CoulombCorrect(const double& eta,
 265                                 const double& radius) {
 266   // Checks radii ... input radius and fRadius
 267   // Calls createLookupTable if neccessary
 268   // Interpolate(linear) between etas in the created lookup table
 269
 270   if (radius < 0.0) {
 271     if (fRadius < 0.0) {
 272       // Both radii are negative
 273       cout << "AliFemtoCoulomb::CoulombCorrect(eta,r) --> input and member radii are negative!" << endl;
 274       cerr << "AliFemtoCoulomb::CoulombCorrect(eta,r) --> input and member radii are negative!" << endl;
 275       assert(0);
 276     }
 277   }
 278   else {
 279     // radius > 0.0
 280     if (radius == fRadius) {
 281       // Both radii are positive and equal
 282       //      cout << "Radii are the same!!!" << endl;
 283     }
 284     else {
 285       // Both radii are positive but not equal
 286       fRadius = radius;
 287       CreateLookupTable(fRadius);
 288     }
 289   }
 290
 291   // Interpolate in eta
 292   return ( CoulombCorrect(eta) );
 293 }
 294
 295 double AliFemtoCoulomb::CoulombCorrect(const AliFemtoPair* pair) {
 296   return ( CoulombCorrect( Eta(pair) ) );;
 297 }
 298
 299 double AliFemtoCoulomb::CoulombCorrect(const AliFemtoPair* pair, const double& radius) {
 300   return ( CoulombCorrect( Eta(pair),radius ) );
 301 }
 302
 303 double AliFemtoCoulomb::Eta(const AliFemtoPair* pair) {
 304   static double px1,py1,pz1,px2,py2,pz2;
 305   static double px1new,py1new,pz1new;
 306   static double px2new,py2new,pz2new;
 307   static double vx1cms,vy1cms,vz1cms;
 308   static double vx2cms,vy2cms,vz2cms;
 309   static double VcmsX,VcmsY,VcmsZ;
 310   static double dv = 0.0;
 311   static double e1,e2,e1new,e2new;
 312   static double psi,theta;
 313   static double beta,gamma;
 314   static double VcmsXnew;
 315
 316   px1 = pair->track1()->FourMomentum().px();
 317   py1 = pair->track1()->FourMomentum().py();
 318   pz1 = pair->track1()->FourMomentum().pz();
 319   e1 = pair->track1()->FourMomentum().e();
 320   px2 = pair->track2()->FourMomentum().px();
 321   py2 = pair->track2()->FourMomentum().py();
 322   pz2 = pair->track2()->FourMomentum().pz();
 323   e2 = pair->track2()->FourMomentum().e();
 324
 325   VcmsX = ( px1+px2 )/( e1+e2 );
 326   VcmsY = ( py1+py2 )/( e1+e2 );
 327   VcmsZ = ( pz1+pz2 )/( e1+e2 );
 328   // Rotate Vcms to x-direction
 329   psi = atan(VcmsY/VcmsX);
 330   VcmsXnew = VcmsX*cos(psi)+VcmsY*sin(psi);
 331   VcmsX = VcmsXnew;
 332   theta = atan(VcmsZ/VcmsX);
 333   VcmsXnew = VcmsX*cos(theta)+VcmsZ*sin(theta);
 334   VcmsX = VcmsXnew;
 335   // Gamma and Beta
 336   beta = VcmsX;
 337   gamma = 1.0/::sqrt( 1.0-beta*beta );
 338
 339   // Rotate p1 and p2 to new frame
 340   px1new = px1*cos(psi)+py1*sin(psi);
 341   py1new = -px1*sin(psi)+py1*cos(psi);
 342   px1 = px1new;
 343   px1new = px1*cos(theta)+pz1*sin(theta);
 344   pz1new = -px1*sin(theta)+pz1*cos(theta);
 345   px1 = px1new;
 346   py1 = py1new;
 347   pz1 = pz1new;
 348
 349   px2new = px2*cos(psi)+py2*sin(psi);
 350   py2new = -px2*sin(psi)+py2*cos(psi);
 351   px2 = px2new;
 352   px2new = px2*cos(theta)+pz2*sin(theta);
 353   pz2new = -px2*sin(theta)+pz2*cos(theta);
 354   px2 = px2new;
 355   py2 = py2new;
 356   pz2 = pz2new;
 357
 358   // Lorentz transform the x component and energy
 359   e1new = gamma*e1 - gamma*beta*px1;
 360   px1new = -gamma*beta*e1 + gamma*px1;
 361   e2new = gamma*e2 - gamma*beta*px2;
 362   px2new = -gamma*beta*e2 + gamma*px2;
 363   px1 = px1new;
 364   px2 = px2new;
 365
 366   // New velocities
 367   vx1cms = px1/e1new;
 368   vy1cms = py1/e1new;
 369   vz1cms = pz1/e1new;
 370   vx2cms = px2/e2new;
 371   vy2cms = py2/e2new;
 372   vz2cms = pz2/e2new;
 373
 374   // Velocity difference in CMS frame
 375   dv = ::sqrt( (vx1cms-vx2cms)*(vx1cms-vx2cms) +
 376              (vy1cms-vy2cms)*(vy1cms-vy2cms) +
 377              (vz1cms-vz2cms)*(vz1cms-vz2cms) );
 378
 379   return ( fZ1Z2*fine_structure_const/(dv) );
 380 }
 381
 382 TH1D* AliFemtoCoulomb::CorrectionHistogram(const double& mass1, const double& mass2, const int& nBins,
 383                                                 const double& low, const double& high) {
 384   if ( mass1!=mass2 ) {
 385     cout << "Masses not equal ... try again.  No histogram created." << endl;
 386     assert(0);
 387   }
 388   TH1D* correction = new TH1D("correction","Coulomb correction",nBins,low,high);
 389   const double reducedMass = mass1*mass2/(mass1+mass2);
 390   double qInv = low;
 391   //double dQinv = (high-low)/( (double)nBins );
 392   double eta;
 393   for (int ii=1; ii<=nBins; ii++)
 394     {
 395       qInv = correction->GetBinCenter(ii);
 396       eta = 2.0*fZ1Z2*reducedMass*fine_structure_const/( qInv );
 397       CoulombCorrect( eta );
 398       correction->Fill( qInv, CoulombCorrect(eta,fRadius) );
 399     }
 400
 401   return (correction);
 402 }
 403
 404 #ifdef __ROOT__
 405 TH1D* AliFemtoCoulomb::CorrectionHistogram(const TH1D* histo, const double mass) {
 406
 407   TH1D* correction = (TH1D*) ((TH1D*)histo)->Clone();
 408   correction->Reset();
 409   correction->SetDirectory(0);
 410   int    nBins = correction->GetXaxis()->GetNbins();
 411   const double reducedMass = 0.5*mass;
 412   double qInv;
 413   double eta;
 414   for (int ii=1; ii<=nBins; ii++)
 415     {
 416       qInv = correction->GetBinCenter(ii);
 417       eta = 2.0*fZ1Z2*reducedMass*fine_structure_const/( qInv );
 418       correction->Fill( qInv, CoulombCorrect(eta,fRadius) );
 419     }
 420
 421   return (correction);
 422 }
 423
 424 TH3D* AliFemtoCoulomb::CorrectionHistogram(const TH3D* histo, const double mass) {
 425
 426   TH3D* correction = (TH3D*) ((TH3D*)histo)->Clone();
 427   correction->Reset();
 428   correction->SetDirectory(0);
 429   int    nBinsX = correction->GetXaxis()->GetNbins();
 430   int    nBinsY = correction->GetYaxis()->GetNbins();
 431   int    nBinsZ = correction->GetZaxis()->GetNbins();
 432   const double reducedMass = 0.5*mass;
 433   double eta;
 434   double qInv;
 435   int binNumber;
 436   for (int ii=1; ii<=nBinsX; ii++) {
 437     for (int iii=1; iii<=nBinsY; iii++) {
 438       for (int iv=1; iv<=nBinsZ; iv++) {
 439         binNumber = histo->GetBin(ii,iii,iv);
 440         qInv = histo->GetBinContent(binNumber);
 441         eta = 2.0*fZ1Z2*reducedMass*fine_structure_const/( qInv );
 442         correction->SetBinContent(binNumber, CoulombCorrect(eta,fRadius) );
 443       }
 444     }
 445   }
 446   return (correction);
 447 }
 448 #endif
 449
 450 double AliFemtoCoulomb::CoulombCorrect(const double& mass, const double& charge,
 451                                     const double& radius, const double& qInv) {
 452   fRadius = radius;
 453   fZ1Z2 = charge;
 454   const double reducedMass = 0.5*mass; // must be same mass particles
 455   double eta = 2.0*fZ1Z2*reducedMass*fine_structure_const/( qInv );
 456   return ( CoulombCorrect(eta,fRadius) );
 457 }