made const array arguments of GetPredictedChi2 methods
[u/mrichter/AliRoot.git] / STRUCT / EvalAbso.C
CommitLineData
fd2bbc37 1void EvalAbso()
2{
3// Material numbers
4 enum {kC=6, kAl=9, kFe=10, kCu=11, kW=12, kPb=13,
5 kNiCuW=21, kVacuum=16, kAir=15, kConcrete=17,
6 kPolyCH2=18, kSteel=10, kInsulation=14, kPolyCc=20};
7
8 AliABSO *pABSO = (AliABSO*) gAlice->GetModule("ABSO");
9
10 Int_t nl[2];
11 nl[0] = pABSO->NumberOfLayers(0);
12 nl[1] = pABSO->NumberOfLayers(1);
13
14 Int_t mat[2][15], mat[2][15];
15 Float_t zmin[2][15], zmin[2][15], zmax[2][15], zmax[2][15];
16 Int_t i, j;
17
18 for (j=0; j< 2; j++) {
19 for (i=0; i< nl[j]; i++) {
20 mat[j][i] = pABSO->MaterialOfLayer(j,i)-1599;
21 zmax[j][i] = pABSO->ZPositionOfLayer(j,i);
22 if (i+1 < nl[j]) zmin[j][i+1] = pABSO->ZPositionOfLayer(j,i);
23 }
24 zmin[j][0]=0.;
25 }
26 Float_t l = zmax[0][nl[0]-1];
27
28
29//
30
31//
32// 1. Limiting angular resolution in the Branson formalism
33//
34 Float_t f0,f1,f2;
35 Float_t a, z, dens, radl, absl;
36 char name[21];
37
38 for (j=0; j< 2; j++) {
92664bc8 39 printf("\n A Z ZPos DZ dens radL absL");
fd2bbc37 40 printf("\n==================================================================================");
41 f0=f1=f2=0;
42 for (i=0; i< nl[j]; i++) {
43 pABSO->AliGetMaterial(mat[j][i], name, a, z, dens, radl, absl);
44 Float_t dz = zmax[j][i]-zmin[j][i];
45 f0 += dz/radl;
46 f1 += dz/(2.*radl)*(2.*zmin[j][i]+dz);
47 f2 += dz/(3.*radl)*(3.*zmin[j][i]*zmin[j][i]+3.*dz*zmin[j][i]+dz*dz);
92664bc8 48 printf("\n %3d %14s %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f %12.3f",
49 i+1, name, a, z, (zmax[j][i]+zmin[j][i])/2., dz, dens, radl, absl);
fd2bbc37 50 }
51//
52 Float_t deltaThetaB = 0.0136 * TMath::Sqrt(f0-f1*f1/f2);
53 Float_t lBranson = f2/f1;
54 printf("\n==================================================================================");
55 printf("\n Branson plane: %12.3f (cm)", lBranson);
56 printf("\n Limiting resolution (Branson) : %12.3f (mrad)", deltaThetaB*1000.);
57 }
58//
59// 2. Limiting Resolution as a function of thickness of extra Fe
60//
61 Float_t zmin0[15], zmax0[15];
62 Float_t x[100], y[100];
63
64 for (i=0; i< nl[0]; i++) {
65 zmin0[i] = zmin[0][i];
66 zmax0[i] = zmax[0][i];
67 }
68
69 Float_t ds = (zmax[0][3]-zmin[0][2])/100.;
70 Int_t j;
71
72 for (j=0; j< 100; j++) {
73 f0=f1=f2=0;
74 Float_t zFe = zmin[0][2]+Float_t(j)*ds;
75 zmax0[2]=zFe;
76 zmin0[3]=zFe;
77 for (i=0; i< nl[0]; i++) {
78 Float_t a, z, dens, radl, absl;
79 pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
80 Float_t dz = zmax0[i]-zmin0[i];
81 f0 += dz/radl;
82 f1 += dz/(2.*radl)*(2.*zmin0[i]+dz);
83 f2 += dz/(3.*radl)*(3.*zmin0[i]*zmin0[i]+3.*dz*zmin0[i]+dz*dz);
84 }
85 Float_t deltaThetaB = 0.0136 * TMath::Sqrt(f0-f1*f1/f2);
86 x[j] = zmax[0][3]-zFe;
87 y[j] = deltaThetaB*1000.;
88 }
89 TGraph *sGraph = new TGraph(100, x, y);
90 TCanvas *c = new TCanvas("c"," ",400,10,600,700);
91 c->Divide(2,2);
92 c->cd(1);
93 sGraph->Draw("AC");
94 c->Update();
95//
96// 3. Limiting Resolution as a function of density of Carbon
97//
98 Float_t dRho = (2.5-1.5)/100.;
99
100 for (j=0; j< 100; j++) {
101 f0=f1=f2=0;
102 Float_t rho = 1.5 + Float_t(j)*dRho;
103 for (i=0; i< nl[0]; i++) {
104 Float_t a, z, dens, radl, absl;
105 pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
106 if (i==1) radl*=1.75/rho;
107 Float_t dz = zmax[0][i]-zmin[0][i];
108 f0 += dz/radl;
109 f1 += dz/(2.*radl)*(2.*zmin[0][i]+dz);
110 f2 += dz/(3.*radl)*(3.*zmin[0][i]*zmin[0][i]+3.*dz*zmin[0][i]+dz*dz);
111 }
112 deltaThetaB = 0.0136 * TMath::Sqrt(f0-f1*f1/f2);
113 x[j] = rho;
114 y[j] = deltaThetaB*1000.;
115 }
116 TGraph *dGraph = new TGraph(100, x, y);
117 c->cd(2);
118 dGraph->Draw("AC");
119 c->Update();
120
121//
122// 4. Energy Loss
123//
124//
125 TCanvas *c1 = new TCanvas("c1"," ",400,10,600,700);
126 c1->Divide(2,2);
127
128 const Float_t kAvo=0.60221367;
129
130 char* tmp;
131 tmp = new char[5];
132 strncpy(tmp, "LOSS", 4);
133 tmp[4]='\0';
134
135 Int_t ixst;
136 Float_t ekin[100], dedx[100], de[100], pcut[100], deRad[100];
137 Float_t emin = 2;
138 Float_t emax = 200;
139 Float_t eps = (emax-emin)/100.;
140// 3 < theta < 9
141 for (j=0; j< 100; j++) {
142 ekin[j] = emin + Float_t(j)*eps;
143 de[j]=0;
144 deRad[j]=0;
145 }
146// all losses
147 for (i=0; i< nl[0]; i++) {
148 ((TGeant3*) gMC)->Gftmat(pABSO->GetMatId(mat[0][i]), 5, tmp, 100, ekin, dedx, pcut, ixst);
149 for (j=0; j< 100; j++) de[j]+=dedx[j]*(zmax[0][i]-zmin[0][i])/1000.;
150 }
151// radative losses
152 for (i=0; i< nl[0]; i++) {
153 Float_t a, z, dens, radl, absl;
154 pABSO->AliGetMaterial(mat[0][i], name, a, z, dens, radl, absl);
155 for (Int_t j=0; j<100; j++) {
156 Float_t nor= kAvo*dens/a*(zmax[0][i]-zmin[0][i]);
157 deRad[j] += (((TGeant3*) gMC)->Gbrelm(z,ekin[j],1.e10))*nor;
158 deRad[j] += (((TGeant3*) gMC)->Gprelm(z,ekin[j],1.e10))*nor;
159 }
160 }
161
162 TH2F *hr = new TH2F("hr1","Several graphs in the same pad",2,0,15,2,0,5);
163 hr->SetXTitle("X title");
164 hr->SetYTitle("Y title");
165 c1->cd(1);
166 hr->Draw();
167
168 TGraph *eGraph1 = new TGraph(100, ekin, de);
169 TGraph *eGraphR1 = new TGraph(100, ekin, deRad);
170
171 eGraph1 ->Draw("LP");
172 eGraphR1->Draw("LP");
173
174 hr = new TH2F("hr2","Several graphs in the same pad",2,15,200,2,0,5);
175 c1->cd(2);
176 hr->Draw();
177 eGraph1 ->Draw("LP");
178 eGraphR1->Draw("LP");
179
180
181// 2 < theta < 3
182 for (j=0; j< 100; j++) {
183 ekin[j] = emin + Float_t(j)*eps;
184 de[j]=0;
185 deRad[j]=0;
186 }
187// all losses
188 for (i=0; i< nl[1]; i++) {
189 ((TGeant3*) gMC)->Gftmat(pABSO->GetMatId(mat[1][i]), 5, tmp, 100, ekin, dedx, pcut, ixst);
190 for (j=0; j< 100; j++) de[j]+=dedx[j]*(zmax[1][i]-zmin[1][i])/1000.;
191 }
192// radative losses
193 for (i=0; i< nl[1]; i++) {
194 Float_t a, z, dens, radl, absl;
195 pABSO->AliGetMaterial(mat[1][i], name, a, z, dens, radl, absl);
196 for (Int_t j=0; j<100; j++) {
197 Float_t nor= kAvo*dens/a*(zmax[1][i]-zmin[1][i]);
198 deRad[j] += (((TGeant3*) gMC)->Gbrelm(z,ekin[j],1.e10))*nor;
199 deRad[j] += (((TGeant3*) gMC)->Gprelm(z,ekin[j],1.e10))*nor;
200 }
201 }
202
203 TH2F *hr2 = new TH2F("hr3","Several graphs in the same pad",2,0,15,2,0,5);
204 hr2->SetXTitle("X title");
205 hr2->SetYTitle("Y title");
206 c1->cd(3);
207 hr2->Draw();
208
209 TGraph *eGraph2 = new TGraph(100, ekin, de);
210 TGraph *eGraphR2 = new TGraph(100, ekin, deRad);
211
212 eGraph2 ->Draw("LP");
213 eGraphR2->Draw("LP");
214
215 hr2 = new TH2F("hr4","Several graphs in the same pad",2,15,200,2,0,5);
216 c1->cd(4);
217 hr2->Draw();
218 eGraph2 ->Draw("LP");
219 eGraphR2->Draw("LP");
220
221
222 c->Update();
223
224}
225
226
227
228