bug fixed for alignment, removed alignment database access from AliPMDUtility class
[u/mrichter/AliRoot.git] / 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)*(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]);
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      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)));
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
140 void 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
150 void AliLhcProcessIBS::Record()
151 {
152   // Record monitor quantities
153     fTauxArray[fAccelerator->Nt()] = fTaux/3600.;
154     fTaueArray[fAccelerator->Nt()] = fTaue/3600.;
155 }
156
157
158 void 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
200 AliLhcProcessIBS& AliLhcProcessIBS::operator=(const  AliLhcProcessIBS & /*rhs*/)
201 {
202 // Assignment operator
203     return *this;
204 }
205
206 Double_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
212   const Double_t kbc = 0.5772;
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);
218   return (1.0-3.0*x2)*p*q*(2.0*TMath::Log(0.5*cc*(p+q))-kbc);
219 }
220