1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
31 #include <TMultiGraph.h>
33 #include "AliLhcProcessIBS.h"
35 #include "AliLhcIRegion.h"
36 #include "AliLhcBeam.h"
38 ClassImp(AliLhcProcessIBS)
40 Double_t func(Double_t *x, Double_t *par);
42 AliLhcProcessIBS::AliLhcProcessIBS(AliLHC* lhc, const char* name, const char* title)
43 :AliLhcProcess(lhc,name,title),
52 for (Int_t i = 0; i < 2; i++) {
59 AliLhcProcessIBS::AliLhcProcessIBS(const AliLhcProcessIBS& ibs):
69 for (Int_t i = 0; i < 2; i++) {
77 AliLhcProcessIBS::~AliLhcProcessIBS()
83 void AliLhcProcessIBS::Init()
86 const Float_t r0=1.535e-16;
88 printf("\n Initializing Process %s", GetName());
89 printf("\n ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^");
91 fIRegions = fAccelerator->IRegions();
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();
100 void AliLhcProcessIBS::Evolve(Float_t dt)
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++) {
107 Float_t sige = fBeam[i]->EnergySpread();
108 Float_t avbeta = fAccelerator->AverageBeta();
109 Float_t avd = fAccelerator->AverageDisp();
111 Float_t gamma = fBeam[i]->Gamma();
112 Float_t epsx = fBeam[i]->Emittance()*gamma;
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);
119 Float_t ssigy = ssigx;
121 Float_t asd = fBeam[i]->N()*fR[i]*fR[i]*fE[i]/
122 (16.*TMath::Pi()*gamma*epsx*epsy*epse);
125 Float_t d = sige*avd/sigx;
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]);
135 TF1 *fct = new TF1("func",func, 0., 1., 3);
137 Float_t f = (Float_t) 8.0*TMath::Pi()*
138 fct->Integral(0., 1., par, 1.e-5);
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);
150 void AliLhcProcessIBS::SetMonitor(Int_t n)
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];
160 void AliLhcProcessIBS::Record()
162 // Record monitor quantities
163 fTauxArray[fAccelerator->Nt()] = fTaux/3600.;
164 fTaueArray[fAccelerator->Nt()] = fTaue/3600.;
168 void AliLhcProcessIBS::DrawPlots()
170 // Draw monitor plots
171 Float_t* t = fAccelerator->TimeA();
173 TH1 *t1 = new TH1F("t1","Hor. IBS growth time",fNmax,0,t[fNmax]);
175 t1->SetMaximum(fTauxArray[fNmax]*1.1);
177 t1->GetXaxis()->SetTitle("t (h)");
178 t1->GetYaxis()->SetTitle("tau_x (t)");
180 TH1 *t2 = new TH1F("t2","Long. IBS growth time",fNmax,0,t[fNmax]);
182 t2->SetMaximum(fTaueArray[fNmax]*1.1);
184 t2->GetXaxis()->SetTitle("t (h)");
185 t2->GetYaxis()->SetTitle("tau_l (t)");
187 TGraph* grTaux = new TGraph(fNmax, fAccelerator->TimeA(), fTauxArray);
188 grTaux->SetHistogram(t1);
190 TGraph* grTaue = new TGraph(fNmax, fAccelerator->TimeA(), fTaueArray);
191 grTaue->SetHistogram(t2);
192 grTaue->SetLineStyle(2);
194 TMultiGraph* mg = new TMultiGraph();
198 TCanvas *c3 = new TCanvas("c3","IBS", 200, 10, 700, 500);
201 mg->GetXaxis()->SetTitle("t (h)");
202 mg->GetYaxis()->SetTitle("IBS Growth Time (h)");
210 AliLhcProcessIBS& AliLhcProcessIBS::operator=(const AliLhcProcessIBS & /*rhs*/)
212 // Assignment operator
216 Double_t func(Double_t *x, Double_t *par)
220 Double_t cc = par[2];
222 const Double_t kbc = 0.5772;
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);