Removing the coding violation
[u/mrichter/AliRoot.git] / LHC / multi.C
1 Double_t fs1(Double_t *x, Double_t* dum);
2 Double_t fs2(Double_t *x, Double_t* dum);
3 Double_t fs3(Double_t *x, Double_t* dum);
4
5 static Int_t n;        // array size
6 static Double_t* as1;  // dsigma/db
7 static Double_t* as2;  // <N>(b)
8 static 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  
17 void multi()
18 {
19     const Float_t drift = 90;
20     
21 //  open file and get graph objects
22     TFile*  file = new TFile("DsigmaDb.root", "read");
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
133 Double_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
150 Double_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
163 Double_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