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