Assymmetry due TDI taken into account.
[u/mrichter/AliRoot.git] / LHC / AliLhcBeam.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 /*
17 $Log$
18 */
19
20 #include "AliLhcBeam.h"
21 #include "AliLHC.h"
22 #include <TMath.h>
23 #include <TCanvas.h>
24 #include <TGraph.h>
25 #include <TMultiGraph.h>
26
27 ClassImp(AliLhcBeam)
28
29 AliLhcBeam::AliLhcBeam(AliLHC* lhc)
30 {
31 // Constructor
32   fAccelerator = lhc;
33 }
34
35 AliLhcBeam::AliLhcBeam(const AliLhcBeam& beam)
36 {
37 // copy constructor
38 }
39
40 AliLhcBeam::~AliLhcBeam()
41 {
42 // Destructor
43
44 }
45
46 void AliLhcBeam::Init()
47 {
48   // Initialization
49   printf("\n Initializing Beam");
50   printf("\n ^^^^^^^^^^^^^^^^^");
51   // Scale energy with regidity 
52   fEnergy     *= fZ/fA;
53   fGamma       = fEnergy/0.938272;
54   fEmittance   = fNEmittance/fGamma;
55   fEmittance0  = fEmittance;
56   fEmittanceL *= fZ;
57   fEmittanceL0 = fEmittanceL; 
58   fN0=fN;
59
60   printf("\n Beam Energy                :%10.3e GeV", fEnergy);
61   printf("\n Beam Normalized Emittance  :%10.3e cm ", fNEmittance);
62   printf("\n Beam Particles per Bunch   :%10.3e    ", fN);
63 }
64
65 void AliLhcBeam::RemoveParticles(Float_t loss)
66 {
67   fN-=loss;
68 }
69
70 void AliLhcBeam::IncreaseEmittance(Float_t de, Float_t del)
71 {
72   fEmittance    *= (1.+de);
73   fEmittanceL   *= (1.+del);
74   fEnergySpread *= (1.+del);
75 }
76
77 AliLhcBeam& AliLhcBeam::operator=(const  AliLhcBeam & rhs)
78 {
79 // Assignment operator
80     return *this;
81 }
82 void AliLhcBeam::SetMonitor(Int_t n)
83 {
84   fNmax = n;
85   if (fEmittanceArray)  delete fEmittanceArray;
86   if (fEmittanceLArray) delete fEmittanceLArray;
87
88
89   fEmittanceArray  = new Float_t[n];
90   fEmittanceLArray = new Float_t[n];
91 }
92
93 void AliLhcBeam::Record()
94 {
95     fEmittanceArray [fAccelerator->Nt()] = fEmittance/fEmittance0;
96     fEmittanceLArray[fAccelerator->Nt()] = fEmittanceL/fEmittanceL0;
97 }
98
99
100 void AliLhcBeam::DrawPlots()
101 {
102   // Draw monitor plots
103   Float_t* t =  fAccelerator->TimeA();
104   
105   TH1 *e1 = new TH1F("e1","Hor. Emittance",fNmax,0,t[fNmax]);
106   e1->SetMinimum(1);
107   e1->SetMaximum(fEmittanceArray[fNmax]*1.1);
108   e1->SetStats(0);
109   e1->GetXaxis()->SetTitle("t (h)");
110   e1->GetYaxis()->SetTitle("rel. Emittance (t)");
111
112   TH1 *e2 = new TH1F("e2","Long. Emittance",fNmax,0,t[fNmax]);
113   e2->SetMinimum(1);
114   e2->SetMaximum(fEmittanceLArray[fNmax]*1.1);
115   e2->SetStats(0);
116   e2->GetXaxis()->SetTitle("t (h)");
117   e2->GetYaxis()->SetTitle("rel. Emittance (t)");
118
119
120   TGraph* grE  = new TGraph(fNmax, t, fEmittanceArray);
121   grE->SetHistogram(e1);
122   TGraph* grEl = new TGraph(fNmax, t, fEmittanceLArray);
123   grEl->SetHistogram(e2);
124   grEl->SetLineStyle(2);
125  
126   TMultiGraph* mg = new TMultiGraph();
127   mg->Add(grE);
128   mg->Add(grEl);
129
130   TCanvas *c2 = new TCanvas("c2","Emittance Increase", 200, 10, 700, 500);
131   c2->SetGrid();
132   mg->Draw("AC");  
133   mg->GetXaxis()->SetTitle("t (h)");
134   mg->GetYaxis()->SetTitle("rel. Emittance(t)");
135   mg->Draw("AC");
136 }
137
138
139
140