Updates for pile-up vertex (F. Prino)
[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
803d1ab0 16/* $Id$ */
11141716 17
21aa51f2 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
0b28fd57 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
11141716 33#include "AliLhcProcessIBS.h"
34#include "AliLHC.h"
35#include "AliLhcIRegion.h"
36#include "AliLhcBeam.h"
37
11141716 38ClassImp(AliLhcProcessIBS)
39
40 Double_t func(Double_t *x, Double_t *par);
41
42AliLhcProcessIBS::AliLhcProcessIBS(AliLHC* lhc, const char* name, const char* title)
64eb707f 43 :AliLhcProcess(lhc,name,title),
44 fCrossSection(0.),
45 fIRegions(0),
46 fTaux(0.),
47 fTaue(0.),
48 fTauxArray(0),
49 fTaueArray(0)
11141716 50{
51// Constructor
52}
53
64eb707f 54AliLhcProcessIBS::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
11141716 66
67AliLhcProcessIBS::~AliLhcProcessIBS()
68{
69// Destructor
70
71}
72
73void 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
90void AliLhcProcessIBS::Evolve(Float_t dt)
91{
21aa51f2 92//
93// Evolve by one time step dt
11141716 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 //
60e55aee 117 Float_t at = sige*TMath::Sqrt((1.-d)*(1.+d))/(gamma*ssigx);
118 Float_t bt = sige*TMath::Sqrt((1.-d)*(1.+d))/(gamma*ssigy);
119 Float_t ct = sige*TMath::Sqrt((1.-d)*(1.+d))*TMath::Sqrt(4.*sigy/fR[i]);
11141716 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 //
60e55aee 130 const Double_t osq2=1./TMath::Sqrt(2.);
131 fTaux = 1./(asd*f*(d-osq2*at)*(d+osq2*at));
132 fTaue = 1./(asd*f*((1.-d)*(1.+d)));
11141716 133 // printf("\n taux, taue %f %f", taux, taue);
134 // Float_t tauy = -2./at*at/asd/f;
135 fBeam[i]->IncreaseEmittance(dt/fTaux, dt/fTaue);
136 }
137
138}
139
140void AliLhcProcessIBS::SetMonitor(Int_t n)
141{
142 // Initialize Monitor
143 if (fTauxArray) delete fTauxArray;
144 if (fTaueArray) delete fTaueArray;
145 fTauxArray = new Float_t[n];
146 fTaueArray = new Float_t[n];
147 fNmax = n;
148}
149
150void AliLhcProcessIBS::Record()
151{
152 // Record monitor quantities
153 fTauxArray[fAccelerator->Nt()] = fTaux/3600.;
154 fTaueArray[fAccelerator->Nt()] = fTaue/3600.;
155}
156
157
158void AliLhcProcessIBS::DrawPlots()
159{
160 // Draw monitor plots
161 Float_t* t = fAccelerator->TimeA();
162
163 TH1 *t1 = new TH1F("t1","Hor. IBS growth time",fNmax,0,t[fNmax]);
164 t1->SetMinimum(0);
165 t1->SetMaximum(fTauxArray[fNmax]*1.1);
166 t1->SetStats(0);
167 t1->GetXaxis()->SetTitle("t (h)");
168 t1->GetYaxis()->SetTitle("tau_x (t)");
169
170 TH1 *t2 = new TH1F("t2","Long. IBS growth time",fNmax,0,t[fNmax]);
171 t2->SetMinimum(0);
172 t2->SetMaximum(fTaueArray[fNmax]*1.1);
173 t2->SetStats(0);
174 t2->GetXaxis()->SetTitle("t (h)");
175 t2->GetYaxis()->SetTitle("tau_l (t)");
176
177 TGraph* grTaux = new TGraph(fNmax, fAccelerator->TimeA(), fTauxArray);
178 grTaux->SetHistogram(t1);
179
180 TGraph* grTaue = new TGraph(fNmax, fAccelerator->TimeA(), fTaueArray);
181 grTaue->SetHistogram(t2);
182 grTaue->SetLineStyle(2);
183
184 TMultiGraph* mg = new TMultiGraph();
185 mg->Add(grTaux);
186 mg->Add(grTaue);
187
188 TCanvas *c3 = new TCanvas("c3","IBS", 200, 10, 700, 500);
189 c3->SetGrid();
190 mg->Draw("AC");
191 mg->GetXaxis()->SetTitle("t (h)");
192 mg->GetYaxis()->SetTitle("IBS Growth Time (h)");
193 mg->Draw("AC");
194}
195
196
197
198
199
e76f229f 200AliLhcProcessIBS& AliLhcProcessIBS::operator=(const AliLhcProcessIBS & /*rhs*/)
11141716 201{
202// Assignment operator
203 return *this;
204}
205
206Double_t func(Double_t *x, Double_t *par)
207{
208 Double_t a = par[0];
209 Double_t b = par[1];
210 Double_t cc = par[2];
211
21aa51f2 212 const Double_t kbc = 0.5772;
11141716 213 Double_t xx = x[0];
214 Double_t x2=xx*xx;
215 Double_t x1=1.0-x2;
216 Double_t p=1.0/TMath::Sqrt(x2+a*a*x1);
217 Double_t q=1.0/TMath::Sqrt(x2+b*b*x1);
21aa51f2 218 return (1.0-3.0*x2)*p*q*(2.0*TMath::Log(0.5*cc*(p+q))-kbc);
11141716 219}
220