mag. field and geometry initialised from GRP
[u/mrichter/AliRoot.git] / LHC / multi.C
CommitLineData
db650755 1Double_t fs1(Double_t *x, Double_t* dum);
2Double_t fs2(Double_t *x, Double_t* dum);
3Double_t fs3(Double_t *x, Double_t* dum);
4
5static Int_t n; // array size
6static Double_t* as1; // dsigma/db
7static Double_t* as2; // <N>(b)
8static Double_t* ab; // b
9// | PbPb | geom. | hard
10// | 0 - 5 fm | 0 - 10 % | 0 - 42.8 %
11// | 5 - 8.6 fm | 10 - 30 % | 42.8 - 83.5 %
12// | 8.6 - 11.2 fm | 30 - 50 % | 83.5 - 96.8 %
13// | 11.2 - 13.2 fm | 50 - 70 % | 96.8 - 99.6 %
14// | 13.2 - 15.0 fm | 70 - 90 % | 99.6 - 99.9 %
15// | 15.0 - oo fm | 90 -100 % | 99.9 - 100.0 %
16
17void multi()
18{
19 const Float_t drift = 90;
20
21// open file and get graph objects
65fa119c 22 TFile* file = new TFile("DsigmaDb.root", "read");
db650755 23 TGraph* gs1 = gAlice = (TGraph*)(file->Get("Graph;1"));
24 TGraph* gs2 = gAlice = (TGraph*)(file->Get("Graph;2"));
25// initialize arrays
26 n = gs1->GetN();
27 ab = gs1->GetX();
28 as1 = gs1->GetY();
29 as2 = gs2->GetY();
30//
31// function
32//
33// dsigma_mb/db
34 TF1* tf1 = new TF1("tf1", fs1, 0., ab[n-1], 0);
35// dsigma_hard/db
36 TF1* tf2 = new TF1("tf2", fs3, 0., ab[n-1], 0);
37
38 Float_t db = 20./100.;
39 for (Int_t i = 0; i< 100; i++)
40 {
41 Float_t bb = Float_t(i)*db;
42 Float_t s = tf1->Integral(0,bb);
43 Float_t sh = tf2->Integral(0,bb);
44 printf("\n %d %f %f %f %f" , i, bb, s, s/7.769*100., sh/4.327*100.);
45 }
46
47
48// histogram for random distribution
49 TH1F* hs1 = new TH1F("hs1", "dSigma/db", n, 0,ab[n-1]);
50 hs1->Eval(tf1);
51// hs1->Draw();
52 TH1F* hmult = new TH1F("hmult", "Multiplicity", 100, 0, 10000);
53 TH1F* hmulte = new TH1F("hmulte", "Extra Multiplicity", 100, 0, 10000);
54//
55 Float_t t1, t2;
56 Double_t b, mult, dm, multe;
57 Int_t i;
58 Double_t *par;
59 Double_t scale = 8000./as2[0];
60 Double_t dT = 125./2.;
61
62 for (Int_t i=0; i < 1000000; i++)
63 {
64 // impact parameters
65 b = hs1->GetRandom();
66 // multiplicity
67 mult = scale*(Float_t) fs2(&b,par);
68 mult += gRandom->Gaus(0, TMath::Sqrt(mult));
69
70 hmult->Fill(Float_t(mult));
71 t1 = 0;
72 t2 = 0;
73 multe = 0;
74 Double_t arg;
75/*
76 while (1) {
77 arg=0.;
78 while(arg==0.) arg = gRandom->Rndm();
79 t1 -= dT*TMath::Log(arg); // (musec)
80 if (t1 > drift) break;
81 b = hs1->GetRandom();
82 mult = scale*(Float_t) fs2(&b,par);
83 mult += gRandom->Gaus(0, TMath::Sqrt(mult));
84// multe+=mult*(drift-t1)/drift;
85 multe+=mult;
86 }
87
88 while (1) {
89 arg=0.;
90 while(arg==0.) arg = gRandom->Rndm();
91 t2 -= dT*TMath::Log(arg); // (musec)
92 if (t2 > drift) break;
93 b = hs1->GetRandom();
94 mult = scale*(Float_t) fs2(&b,par);
95 mult += gRandom->Gaus(0, TMath::Sqrt(mult));
96// multe+=mult*(drift-t2)/drift;
97 multe+=mult;
98 }
99*/
100 hmulte->Fill(Float_t(multe));
101 }
102 Double_t scale = 1./hmult->Integral();
103 hmult->Scale(scale/100.);
104 printf("\n scale %f \n ", scale);
105
106
107 TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
108// c1->Divide(1,2);
109// c1->cd(1);
110 hmult->Draw();
111// c1->cd(2);
112// hmulte->Draw();
113
114
115 TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
116 c2->Divide(1,2);
117 c2->cd(1);
118 tf1->Draw();
119 tf1->GetHistogram()->SetXTitle("b (fm)");
120 tf1->GetHistogram()->SetYTitle("d#sigma /db (barn/fm)");
121 c2->cd(2);
122 tf2->Draw();
123 tf2->GetHistogram()->SetXTitle("b (fm)");
124 tf2->GetHistogram()->SetYTitle("d#sigma /db (a.u.)");
125 c2->Update();
126 out = new TFile("mult.root", "recreate");
127 hmult->Write();
128 out->Close();
129
130}
131
132
133Double_t fs1(Double_t* x, Double_t* dum)
134{
135 Double_t xx = x[0];
136 Int_t i;
137 Double_t y;
138
139 for (i=0; i<n; i++) {
140
141 if (xx <ab[i]) {
142 y = as1[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as1[i]-as1[i-1]);
143 break;
144 }
145 }
146 return y;
147}
148
149
150Double_t fs2(Double_t *x, Double_t* dum)
151{
152 Double_t xx = x[0];
153 Int_t i;
154 for (i=0; i<n; i++) {
155 if (xx <ab[i]) {
156 Double_t y = as2[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as2[i]-as2[i-1]);
157 break;
158 }
159 }
160 return y;
161}
162
163Double_t fs3(Double_t *x, Double_t* dum)
164{
165 Double_t xx = x[0];
166 Int_t i;
167 Double_t y1, y2;
168
169 for (i=0; i<n; i++) {
170 if (xx <ab[i]) {
171 y1 = as1[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as1[i]-as1[i-1]);
172 break;
173 }
174 }
175
176 for (i=0; i<n; i++) {
177 if (xx <ab[i]) {
178 Double_t y2 = as2[i-1]+(xx-ab[i-1])/(ab[i]-ab[i-1])*(as2[i]-as2[i-1]);
179 break;
180 }
181 }
182
183
184 return y1*y2;
185}
186
187