2 10mu extrusions corrected.
[u/mrichter/AliRoot.git] / STRUCT / EvalAbso.C
1 void 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++) {
39         printf("\n                        A            Z            dens        radL        absL");
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);       
48             printf("\n %3d %14s %12.3f %12.3f %12.3f %12.3f %12.3f",
49                    i+1, name, a, z, dens, radl, absl);
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