LHC/AliLhcProcessIBS.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 /* $Id$ */
  17
  18 //
  19 // Realisation of AliLhcProcess for the fast simulation of the
  20 // Intra Beam Scattering process
  21 // in transverse and longitudinal direction.
  22 // Author: Andreas Morsch
  23 // andreas.morsch@cern.ch
  24 //
  25
  26 #include <TCanvas.h>
  27 #include <TF1.h>
  28 #include <TGraph.h>
  29 #include <TH1F.h>
  30 #include <TMath.h>
  31 #include <TMultiGraph.h>
  32
  33 #include "AliLhcProcessIBS.h"
  34 #include "AliLHC.h"
  35 #include "AliLhcIRegion.h"
  36 #include "AliLhcBeam.h"
  37
  38 ClassImp(AliLhcProcessIBS)
  39
  40   Double_t func(Double_t *x, Double_t *par);
  41
  42 AliLhcProcessIBS::AliLhcProcessIBS(AliLHC* lhc, const char* name, const char* title)
  43     :AliLhcProcess(lhc,name,title),
  44      fCrossSection(0.),
  45      fIRegions(0),
  46      fTaux(0.),
  47      fTaue(0.),
  48      fTauxArray(0),
  49      fTaueArray(0)
  50 {
  51 // Constructor
  52   for (Int_t i = 0; i < 2; i++) {
  53     fBeam[i] = 0;
  54     fR[i]    = 0.;
  55     fE[i]    = 0.;
  56   }
  57 }
  58
  59 AliLhcProcessIBS::AliLhcProcessIBS(const AliLhcProcessIBS& ibs):
  60     AliLhcProcess(ibs),
  61     fCrossSection(0.),
  62     fIRegions(0),
  63     fTaux(0.),
  64     fTaue(0.),
  65     fTauxArray(0),
  66     fTaueArray(0)
  67 {
  68 // Copy Constructor
  69   for (Int_t i = 0; i < 2; i++) {
  70     fBeam[i] = 0;
  71     fR[i]    = 0.;
  72     fE[i]    = 0.;
  73   }
  74 }
  75
  76
  77 AliLhcProcessIBS::~AliLhcProcessIBS()
  78 {
  79 // Destructor
  80
  81 }
  82
  83 void AliLhcProcessIBS::Init()
  84 {
  85   // Initialization
  86    const Float_t r0=1.535e-16;
  87
  88    printf("\n Initializing Process %s", GetName());
  89    printf("\n ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^");
  90
  91    fIRegions = fAccelerator->IRegions();
  92
  93    for (Int_t i = 0; i < 2; i++) {
  94      fBeam[i] = fAccelerator->Beam(i);
  95      fR[i]    = r0*fBeam[i]->Z()*fBeam[i]->Z()/fBeam[i]->A();
  96      fE[i]    = 0.938*fBeam[i]->A();
  97    }
  98 }
  99
 100 void AliLhcProcessIBS::Evolve(Float_t dt)
 101 {
 102 //
 103 // Evolve by one time step dt
 104   printf("\n Here process %s %f:", GetName(), dt);
 105    for (Int_t i=0; i<2; i++) {
 106      // Density
 107      Float_t sige    = fBeam[i]->EnergySpread();
 108      Float_t avbeta  = fAccelerator->AverageBeta();
 109      Float_t avd     = fAccelerator->AverageDisp();
 110
 111      Float_t gamma   = fBeam[i]->Gamma();
 112      Float_t epsx    = fBeam[i]->Emittance()*gamma;
 113      Float_t epsy    = epsx;
 114      Float_t epse    = fBeam[i]->LongEmittance();
 115      Float_t sigxb   = TMath::Sqrt(epsx/gamma*avbeta);
 116      Float_t sigx    = TMath::Sqrt(sigxb*sigxb+(avd*avd*sige*sige));
 117      Float_t ssigx   = TMath::Sqrt(epsx/gamma/avbeta);
 118      Float_t sigy    = sigx;
 119      Float_t ssigy   = ssigx;
 120
 121      Float_t asd = fBeam[i]->N()*fR[i]*fR[i]*fE[i]/
 122        (16.*TMath::Pi()*gamma*epsx*epsy*epse);
 123
 124      // impact parameter
 125      Float_t d    = sige*avd/sigx;
 126      //
 127      Float_t at = sige*TMath::Sqrt((1.-d)*(1.+d))/(gamma*ssigx);
 128      Float_t bt = sige*TMath::Sqrt((1.-d)*(1.+d))/(gamma*ssigy);
 129      Float_t ct = sige*TMath::Sqrt((1.-d)*(1.+d))*TMath::Sqrt(4.*sigy/fR[i]);
 130      //
 131      Double_t par[3];
 132      par[0] = at;
 133      par[1] = bt;
 134      par[2] = ct;
 135      TF1 *fct = new TF1("func",func, 0., 1., 3);
 136
 137      Float_t f = (Float_t) 8.0*TMath::Pi()*
 138          fct->Integral(0., 1., par, 1.e-5);
 139      //
 140      const Double_t osq2=1./TMath::Sqrt(2.);
 141      fTaux =  1./(asd*f*(d-osq2*at)*(d+osq2*at));
 142      fTaue =  1./(asd*f*((1.-d)*(1.+d)));
 143      //     printf("\n taux, taue %f %f", taux, taue);
 144      //     Float_t tauy = -2./at*at/asd/f;
 145      fBeam[i]->IncreaseEmittance(dt/fTaux, dt/fTaue);
 146    }
 147
 148 }
 149
 150 void AliLhcProcessIBS::SetMonitor(Int_t n)
 151 {
 152   // Initialize Monitor
 153   if (fTauxArray) delete fTauxArray;
 154   if (fTaueArray) delete fTaueArray;
 155   fTauxArray = new Float_t[n];
 156   fTaueArray = new Float_t[n];
 157   fNmax = n;
 158 }
 159
 160 void AliLhcProcessIBS::Record()
 161 {
 162   // Record monitor quantities
 163     fTauxArray[fAccelerator->Nt()] = fTaux/3600.;
 164     fTaueArray[fAccelerator->Nt()] = fTaue/3600.;
 165 }
 166
 167
 168 void AliLhcProcessIBS::DrawPlots()
 169 {
 170   // Draw monitor plots
 171   Float_t* t =  fAccelerator->TimeA();
 172
 173   TH1 *t1 = new TH1F("t1","Hor. IBS growth time",fNmax,0,t[fNmax]);
 174   t1->SetMinimum(0);
 175   t1->SetMaximum(fTauxArray[fNmax]*1.1);
 176   t1->SetStats(0);
 177   t1->GetXaxis()->SetTitle("t (h)");
 178   t1->GetYaxis()->SetTitle("tau_x (t)");
 179
 180   TH1 *t2 = new TH1F("t2","Long. IBS growth time",fNmax,0,t[fNmax]);
 181   t2->SetMinimum(0);
 182   t2->SetMaximum(fTaueArray[fNmax]*1.1);
 183   t2->SetStats(0);
 184   t2->GetXaxis()->SetTitle("t (h)");
 185   t2->GetYaxis()->SetTitle("tau_l (t)");
 186
 187   TGraph* grTaux = new TGraph(fNmax, fAccelerator->TimeA(), fTauxArray);
 188   grTaux->SetHistogram(t1);
 189
 190   TGraph* grTaue = new TGraph(fNmax, fAccelerator->TimeA(), fTaueArray);
 191   grTaue->SetHistogram(t2);
 192   grTaue->SetLineStyle(2);
 193
 194   TMultiGraph* mg = new TMultiGraph();
 195   mg->Add(grTaux);
 196   mg->Add(grTaue);
 197
 198   TCanvas *c3 = new TCanvas("c3","IBS", 200, 10, 700, 500);
 199   c3->SetGrid();
 200   mg->Draw("AC");
 201   mg->GetXaxis()->SetTitle("t (h)");
 202   mg->GetYaxis()->SetTitle("IBS Growth Time (h)");
 203   mg->Draw("AC");
 204 }
 205
 206
 207
 208
 209
 210 AliLhcProcessIBS& AliLhcProcessIBS::operator=(const  AliLhcProcessIBS & /*rhs*/)
 211 {
 212 // Assignment operator
 213     return *this;
 214 }
 215
 216 Double_t func(Double_t *x, Double_t *par)
 217 {
 218   Double_t a  = par[0];
 219   Double_t b  = par[1];
 220   Double_t cc = par[2];
 221
 222   const Double_t kbc = 0.5772;
 223   Double_t xx = x[0];
 224   Double_t x2=xx*xx;
 225   Double_t x1=1.0-x2;
 226   Double_t p=1.0/TMath::Sqrt(x2+a*a*x1);
 227   Double_t q=1.0/TMath::Sqrt(x2+b*b*x1);
 228   return (1.0-3.0*x2)*p*q*(2.0*TMath::Log(0.5*cc*(p+q))-kbc);
 229 }
 230