PWG2/UNICOR/AliUnicorCoulomb.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2048, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 //Author: Dariusz Miskowiec
  17 //Date:   2010
  18
  19 //=============================================================================
  20 // Coulomb correlation function
  21 //=============================================================================
  22
  23 #include <TRandom2.h>
  24 #include <TComplex.h>
  25 #include <TFile.h>
  26 #include <TH2.h>
  27 #include "AliUnicorCoulomb.h"
  28
  29 ClassImp(AliUnicorCoulomb)
  30
  31 //=============================================================================
  32 AliUnicorCoulomb::AliUnicorCoulomb(int zz, double mass, double R) : TGraph(1001) {
  33
  34   // constructor
  35   // zz = +1 for pi+pi+, -1 for pi-pi-, etc.
  36   // mass is the reduced mass: m1*m2/(m1+m2)
  37
  38   SetMarkerStyle(20);
  39   int nstep = 10000;
  40
  41   double rx = R;
  42   double ry = R;
  43   double rz = R;
  44   double r0 = 0;
  45   TRandom2 ran;
  46
  47   if (zz>0) SetPoint(0,0,0);
  48   else SetPoint(0,0,1e9);
  49   for (int i=1; i<1000; i++) {
  50     double qinv = i/1000.0;
  51     double k = qinv/2.0;
  52     double sum = 0;
  53     for (int j=0; j<nstep; j++) {
  54       double x0 = ran.Gaus(0.0,rx);
  55       double y0 = ran.Gaus(0.0,ry);
  56       double z0 = ran.Gaus(0.0,rz);
  57       double t0 = ran.Gaus(0.0,r0);
  58       double x1 = ran.Gaus(0.0,rx);
  59       double y1 = ran.Gaus(0.0,ry);
  60       double z1 = ran.Gaus(0.0,rz);
  61       double t1 = ran.Gaus(0.0,r0);
  62       double dx = x1-x0;
  63       double dy = y1-y0;
  64       double dz = z1-z0;
  65       double dt = t1-t0;
  66       dz += dt*k/mass;
  67       sum += WaveFunction2(zz,mass,k,dx,dy,dz);
  68     }
  69     SetPoint(i,qinv,sum/nstep);
  70   }
  71   SetPoint(1000,1000,1);
  72 }
  73 //=============================================================================
  74 double AliUnicorCoulomb::WaveFunction2(int zz, double mass, double k, double x, double y, double z) {
  75
  76   // Square of normalized Coulomb wave function.
  77   // Hypergeometrical function diverges for numerical reasons for large Qinv, z and x.
  78   // Out of the convergence limit I will return 1.0
  79
  80   double qinv12 = 1.2 * 2.0 * k;
  81   double p0 = -9.1569/qinv12;
  82   double p2 = 2.735e-2*qinv12;
  83   double wf2 = 1.0;
  84   int converg = (z>p0+p2*(x*x+y*y));
  85   if (converg) {
  86     TComplex co = WaveFunction(zz,mass,k,x,y,z);
  87     wf2 =  co.Rho2()*Gamow(zz,mass,k);
  88   }
  89   return wf2;
  90 }
  91 //=============================================================================
  92 TComplex AliUnicorCoulomb::WaveFunction(int zz, double mass, double k, double x, double y, double z) {
  93
  94   // Coulomb wave function in parabolic coordinates
  95   // Merzbacher, Quantum Mechanics, p.248
  96   // Landau-Lifshitz, Quantum Mechanics (Non-Relativistic Theory), p.518
  97
  98   double r = sqrt(x*x+y*y+z*z);
  99   //  double ksi = r+z;
 100   double eta = r-z;
 101   TComplex jjj(0,1);
 102   TComplex a = -zz*mass/k/137.036*jjj;
 103   TComplex b(1.0,0.0);
 104   TComplex c = jjj*k*eta/0.197327;
 105   TComplex wf = TComplex::Exp(jjj*k*z/0.197327)*F1(a,b,c);
 106   //printf("%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f \n",a.Re(),a.Im(),b.Re(),b.Im(),c.Re(),c.Im(),F1(a,b,c).Rho());
 107   return wf;
 108 }
 109 //=============================================================================
 110 double AliUnicorCoulomb::Gamow(int zz, double m, double k) {
 111
 112   // Squared Coulomb wave function in infinity = Gamow
 113   // zz - product of charges
 114   // m - reduced mass
 115   // k = Qinv/2 for equal masses
 116
 117   double n = -zz*m/k/137.036;
 118   double a = 2*TMath::Pi()*n;
 119   double ga = a/(1.0-exp(-a));
 120   return ga;
 121 }
 122 //=============================================================================
 123 TComplex AliUnicorCoulomb::F1(TComplex alpha, TComplex gamma, TComplex z) {
 124
 125   // confluent hypergeometric function 1F1
 126
 127   double prec = 0.0000001;
 128   TComplex term(1.0,0.0);
 129   TComplex sum(1.0,0.0);
 130   for (int n=1; n<1000; n++) {
 131     double u = (double) n;
 132     TComplex v = (alpha+u-1.0)/(gamma+u-1.0);
 133     TComplex w = z/u;
 134     term *= v;
 135     term *= w;
 136     sum += term;
 137     //    printf("%10d %10.3f %10.3f %10.3f %10.3f\n", n, v.Rho(), w.Rho(), term.Rho(), sum.Rho());
 138     if (TComplex::Abs(term/sum) < prec) return sum;
 139   }
 140   printf("F1 Maximum number of iterations reached\n");
 141   return sum;
 142 }
 143 //=============================================================================
 144  void AliUnicorCoulomb::Makehist(int zz, double m, const char *outfil) {
 145
 146   // Make a two-dim histogram coulomb(R,Q) and store it on outfil.
 147   // Later to be used via TH2::Interpolate().
 148   // zz and m are the product of charges and the reduced mass.
 149
 150   TH2D *hi = new TH2D("co","coulomb",11,-0.5,10.5,999,0.0005,0.9995);
 151   hi->SetXTitle("R (fm)");
 152   hi->SetYTitle("Q (GeV/c)");
 153   AliUnicorCoulomb *co = 0;
 154   for (int i=1; i<=hi->GetNbinsX(); i++) {
 155     double R = hi->GetXaxis()->GetBinCenter(i);
 156     printf("R = %.1f fm\n",R);
 157     co = new AliUnicorCoulomb(zz, m, R);
 158     for (int j=1; j<=hi->GetNbinsY(); j++) {
 159       double Q = hi->GetYaxis()->GetBinCenter(j);
 160       hi->SetBinContent(i,j,co->Eval(Q));
 161     }
 162     delete co;
 163   }
 164   TFile::Open(outfil,"recreate");
 165   hi->Write();
 166   gFile->Close();
 167 }
 168 //=============================================================================