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