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