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