PWGCF/FEMTOSCOPY/AliFemto/AliFemtoCoulomb.cxx

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