Line data Source code
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 : /* $Id$ */
17 :
18 : //-------------------------------------------------------------------------
19 : // Class AliMUONFastTracking
20 : //
21 : // Manager for the fast simulation of tracking in the muon spectrometer
22 : // This class reads the lookup tables containing the parameterization
23 : // of the deltap, deltatheta, deltaphi for different background levels
24 : // and provides the related smeared parameters.
25 : // Used by AliFastMuonTrackingEff, AliFastMuonTrackingAcc,
26 : // AliFastMuonTrackingRes.
27 : //-------------------------------------------------------------------------
28 :
29 : #include <stdio.h>
30 : #include <stdlib.h>
31 : #include <string.h>
32 :
33 : #include <Riostream.h>
34 : #include <TF1.h>
35 : #include <TFile.h>
36 : #include <TH3.h>
37 : #include <TMath.h>
38 : #include <TRandom.h>
39 : #include <TSpline.h>
40 :
41 : #include "AliMUONFastTracking.h"
42 : #include "AliMUONFastTrackingEntry.h"
43 :
44 : using std::cout;
45 : using std::endl;
46 12 : ClassImp(AliMUONFastTracking)
47 :
48 :
49 : AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
50 :
51 : static Double_t FitP(Double_t *x, Double_t *par){
52 : // Fit function
53 0 : Double_t dx = x[0] - par[0];
54 0 : Double_t dx2 = x[0] - par[4];
55 0 : Double_t sigma = par[1] * ( 1 + par[2] * dx);
56 0 : if (sigma == 0) {
57 :
58 0 : return 0.;
59 : }
60 0 : Double_t fasymm = TMath::Exp(-0.5 * dx * dx / (sigma * sigma));
61 0 : Double_t sigma2 = par[1] * par[5];
62 0 : Double_t fgauss = TMath::Exp(-0.5 * dx2 * dx2 / (sigma2 * sigma2));
63 0 : Double_t value = fasymm + par[3] * fgauss;
64 0 : return TMath::Abs(value);
65 0 : }
66 :
67 : AliMUONFastTracking::AliMUONFastTracking(const AliMUONFastTracking & ft):
68 0 : TObject(),
69 0 : fNbinp(10),
70 0 : fPmin(0.),
71 0 : fPmax(200.),
72 0 : fDeltaP((fPmax-fPmin)/fNbinp),
73 0 : fNbintheta(10),
74 0 : fThetamin(2.),
75 0 : fThetamax(9.),
76 0 : fDeltaTheta((fThetamax-fThetamin)/fNbintheta),
77 0 : fNbinphi(10),
78 0 : fPhimin(-180.),
79 0 : fPhimax(180.),
80 0 : fDeltaPhi((fPhimax-fPhimin)/fNbinphi),
81 0 : fPrintLevel(1),
82 0 : fBkg(0.),
83 0 : fSpline(0),
84 0 : fClusterFinder(kOld)
85 0 : {
86 : // Copy constructor
87 0 : ft.Copy(*this);
88 0 : }
89 :
90 :
91 : AliMUONFastTracking* AliMUONFastTracking::Instance()
92 : {
93 : // Set random number generator
94 0 : if (fgMUONFastTracking) {
95 0 : return fgMUONFastTracking;
96 : } else {
97 0 : fgMUONFastTracking = new AliMUONFastTracking();
98 0 : return fgMUONFastTracking;
99 : }
100 0 : }
101 :
102 0 : AliMUONFastTracking::AliMUONFastTracking():
103 0 : fNbinp(10),
104 0 : fPmin(0.),
105 0 : fPmax(200.),
106 0 : fDeltaP((fPmax-fPmin)/fNbinp),
107 0 : fNbintheta(10),
108 0 : fThetamin(2.),
109 0 : fThetamax(9.),
110 0 : fDeltaTheta((fThetamax-fThetamin)/fNbintheta),
111 0 : fNbinphi(10),
112 0 : fPhimin(-180.),
113 0 : fPhimax(180.),
114 0 : fDeltaPhi((fPhimax-fPhimin)/fNbinphi),
115 0 : fPrintLevel(1),
116 0 : fBkg(0.),
117 0 : fSpline(0),
118 0 : fClusterFinder(kOld)
119 0 : {
120 : //
121 : // constructor
122 : //
123 0 : for (Int_t i = 0; i<20;i++) {
124 0 : for (Int_t j = 0; j<20; j++) {
125 0 : for (Int_t k = 0; k<20; k++) {
126 0 : fFitp[i][j][k] = 0x0;
127 : }
128 : }
129 : }
130 0 : }
131 :
132 : void AliMUONFastTracking::Init(Float_t bkg)
133 : {
134 : //
135 : // Initialization
136 : //
137 0 : for (Int_t ip=0; ip< fNbinp; ip++){
138 0 : for (Int_t itheta=0; itheta< fNbintheta; itheta++){
139 0 : for (Int_t iphi=0; iphi< fNbinphi; iphi++){
140 0 : fCurrentEntry[ip][itheta][iphi] = new AliMUONFastTrackingEntry;
141 0 : for (Int_t ibkg=0; ibkg<4; ibkg++){
142 0 : fEntry[ip][itheta][iphi][ibkg] = new AliMUONFastTrackingEntry;
143 : }
144 : }
145 : }
146 : }
147 :
148 0 : char filename [100];
149 0 : if (fClusterFinder==kOld) snprintf (filename, 100, "$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root");
150 0 : else snprintf (filename, 100, "$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT-AZ.root");
151 :
152 0 : TFile *file = new TFile(filename);
153 0 : ReadLUT(file);
154 0 : SetBackground(bkg);
155 0 : UseSpline(0);
156 0 : }
157 :
158 :
159 : void AliMUONFastTracking::ReadLUT(TFile* file)
160 : {
161 : //
162 : // read the lookup tables from file
163 : //
164 0 : TH3F *heff[5][3], *hacc[5][3], *hmeanp, *hsigmap, *hsigma1p, *hchi2p;
165 : TH3F *hnormg2, *hmeang2, *hsigmag2, *hmeantheta, *hsigmatheta, *hchi2theta;
166 : TH3F *hmeanphi, *hsigmaphi, *hchi2phi;
167 0 : char tag[40], tag2[40];
168 :
169 0 : printf ("Reading parameters from LUT file %s...\n",file->GetName());
170 :
171 : const Float_t kBkg[4] = {0, 0.5, 1, 2};
172 0 : for (Int_t ibkg=0; ibkg<4; ibkg++) {
173 0 : snprintf (tag, 40, "BKG%g",kBkg[ibkg]);
174 0 : file->cd(tag);
175 0 : for (Int_t isplp = 0; isplp<kSplitP; isplp++) {
176 0 : for (Int_t ispltheta = 0; ispltheta<kSplitTheta; ispltheta++) {
177 0 : snprintf (tag2, 40, "heff[%d][%d]",isplp,ispltheta);
178 0 : heff[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
179 0 : snprintf (tag2, 40, "hacc[%d][%d]",isplp,ispltheta);
180 0 : hacc[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
181 : }
182 : }
183 0 : hmeanp = (TH3F*)gDirectory->Get("hmeanp");
184 0 : hsigmap = (TH3F*)gDirectory->Get("hsigmap");
185 0 : hsigma1p = (TH3F*)gDirectory->Get("hsigma1p");
186 0 : hchi2p = (TH3F*)gDirectory->Get("hchi2p");
187 0 : hnormg2 = (TH3F*)gDirectory->Get("hnormg2");
188 0 : hmeang2 = (TH3F*)gDirectory->Get("hmeang2");
189 0 : hsigmag2 = (TH3F*)gDirectory->Get("hsigmag2");
190 0 : hmeantheta = (TH3F*)gDirectory->Get("hmeantheta");
191 0 : hsigmatheta = (TH3F*)gDirectory->Get("hsigmatheta");
192 0 : hchi2theta = (TH3F*)gDirectory->Get("hchi2theta");
193 0 : hmeanphi = (TH3F*)gDirectory->Get("hmeanphi");
194 0 : hsigmaphi = (TH3F*)gDirectory->Get("hsigmaphi");
195 0 : hchi2phi = (TH3F*)gDirectory->Get("hchi2phi");
196 :
197 0 : for (Int_t ip=0; ip<fNbinp ;ip++) {
198 0 : for (Int_t itheta=0; itheta<fNbintheta ;itheta++) {
199 0 : for (Int_t iphi=0; iphi<fNbinphi ;iphi++) {
200 0 : Float_t p = fPmin + fDeltaP * (ip + 0.5);
201 0 : Float_t theta = fThetamin + fDeltaTheta * (itheta + 0.5);
202 0 : Float_t phi = fPhimin + fDeltaPhi * (iphi + 0.5);
203 :
204 0 : fEntry[ip][itheta][iphi][ibkg]->SetP(p);
205 0 : fEntry[ip][itheta][iphi][ibkg]->SetMeanp(hmeanp->GetBinContent(ip+1,itheta+1,iphi+1));
206 0 : fEntry[ip][itheta][iphi][ibkg]->SetSigmap(TMath::Abs(hsigmap->GetBinContent(ip+1,itheta+1,iphi+1)));
207 0 : fEntry[ip][itheta][iphi][ibkg]->SetSigma1p(hsigma1p->GetBinContent(ip+1,itheta+1,iphi+1));
208 0 : fEntry[ip][itheta][iphi][ibkg]->SetChi2p(hchi2p->GetBinContent(ip+1,itheta+1,iphi+1));
209 0 : fEntry[ip][itheta][iphi][ibkg]->SetNormG2(hnormg2->GetBinContent(ip+1,itheta+1,iphi+1));
210 0 : fEntry[ip][itheta][iphi][ibkg]->SetMeanG2(hmeang2->GetBinContent(ip+1,itheta+1,iphi+1));
211 0 : if (ibkg == 0) fEntry[ip][itheta][iphi][ibkg]->SetSigmaG2(9999);
212 0 : else fEntry[ip][itheta][iphi][ibkg]->SetSigmaG2(hsigmag2->GetBinContent(ip+1,itheta+1,iphi+1));
213 0 : fEntry[ip][itheta][iphi][ibkg]->SetTheta(theta);
214 0 : fEntry[ip][itheta][iphi][ibkg]->SetMeantheta(hmeantheta->GetBinContent(ip+1,itheta+1,iphi+1));
215 0 : fEntry[ip][itheta][iphi][ibkg]->SetSigmatheta(TMath::Abs(hsigmatheta->GetBinContent(ip+1,itheta+1,iphi+1)));
216 0 : fEntry[ip][itheta][iphi][ibkg]->SetChi2theta(hchi2theta->GetBinContent(ip+1,itheta+1,iphi+1));
217 0 : fEntry[ip][itheta][iphi][ibkg]->SetPhi(phi);
218 0 : fEntry[ip][itheta][iphi][ibkg]->SetMeanphi(hmeanphi->GetBinContent(ip+1,itheta+1,iphi+1));
219 0 : fEntry[ip][itheta][iphi][ibkg]->SetSigmaphi(TMath::Abs(hsigmaphi->GetBinContent(ip+1,itheta+1,iphi+1)));
220 0 : fEntry[ip][itheta][iphi][ibkg]->SetChi2phi(hchi2phi->GetBinContent(ip+1,itheta+1,iphi+1));
221 0 : for (Int_t i=0; i<kSplitP; i++) {
222 0 : for (Int_t j=0; j<kSplitTheta; j++) {
223 0 : fEntry[ip][itheta][iphi][ibkg]->SetAcc(i,j,hacc[i][j]->GetBinContent(ip+1,itheta+1,iphi+1));
224 0 : fEntry[ip][itheta][iphi][ibkg]->SetEff(i,j,heff[i][j]->GetBinContent(ip+1,itheta+1,iphi+1));
225 : }
226 : }
227 : } // iphi
228 : } // itheta
229 : } // ip
230 : } // ibkg
231 :
232 0 : TGraph *graph = new TGraph(3);
233 0 : TF1 *f = new TF1("f","[0]+[1]*x");
234 :
235 0 : for (Int_t ip=0; ip< fNbinp; ip++){
236 0 : for (Int_t itheta=0; itheta< fNbintheta; itheta++){
237 0 : for (Int_t iphi=0; iphi< fNbinphi; iphi++){
238 0 : graph->SetPoint(0,0.5,fEntry[ip][itheta][iphi][1]->GetSigmaG2());
239 0 : graph->SetPoint(1,1,fEntry[ip][itheta][iphi][2]->GetSigmaG2());
240 0 : graph->SetPoint(2,2,fEntry[ip][itheta][iphi][3]->GetSigmaG2());
241 0 : graph->Fit("f","q");
242 0 : fEntry[ip][itheta][iphi][0]->SetSigmaG2(f->Eval(0));
243 : }
244 : }
245 : }
246 0 : f->Delete();
247 0 : graph->Delete();
248 0 : printf ("parameters read. \n");
249 0 : }
250 :
251 : void AliMUONFastTracking::GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
252 : Int_t &nbintheta, Float_t &thetamin,
253 : Float_t &thetamax,
254 : Int_t &nbinphi, Float_t &phimin, Float_t &phimax) const
255 : {
256 : //
257 : // gets the binning for the discrete parametrizations in the lookup table
258 : //
259 0 : nbinp = fNbinp;
260 0 : pmin = fPmin;
261 0 : pmax = fPmax;
262 0 : nbintheta = fNbintheta;
263 0 : thetamin = fThetamin;
264 0 : thetamax = fThetamax;
265 0 : nbinphi = fNbinphi;
266 0 : phimin = fPhimin;
267 0 : phimax = fPhimax;
268 0 : }
269 :
270 :
271 : void AliMUONFastTracking::GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi,
272 : Int_t charge, Int_t &ip, Int_t &itheta,
273 : Int_t &iphi) const
274 : {
275 : //
276 : // gets the id of the cells in the LUT for a given (p,theta,phi, charge)
277 : //
278 0 : if (charge < 0) phi = -phi;
279 0 : ip = Int_t (( p - fPmin ) / fDeltaP);
280 0 : itheta = Int_t (( theta - fThetamin ) / fDeltaTheta);
281 0 : iphi = Int_t (( phi - fPhimin ) / fDeltaPhi);
282 :
283 :
284 0 : if (ip< 0) ip = 0;
285 0 : if (ip>= fNbinp) ip = fNbinp-1;
286 0 : if (itheta< 0) itheta = 0;
287 0 : if (itheta>= fNbintheta) itheta = fNbintheta-1;
288 :
289 0 : if (iphi< 0) iphi = 0;
290 0 : if (iphi>= fNbinphi) iphi = fNbinphi-1;
291 0 : }
292 :
293 : void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta,
294 : Int_t &nSplitP, Int_t &nSplitTheta) const
295 : {
296 : //
297 : // the first cell is splitted in more bins for theta and momentum
298 : // parameterizations. Get the number of divisions for the splitted bins
299 : //
300 0 : if (ip==0) nSplitP = 5;
301 0 : else nSplitP = 2;
302 0 : if (itheta==0) nSplitTheta = 3;
303 0 : else nSplitTheta = 1;
304 0 : }
305 :
306 : Float_t AliMUONFastTracking::Efficiency(Float_t p, Float_t theta,
307 : Float_t phi, Int_t charge){
308 : //
309 : // gets the tracking efficiency
310 : //
311 0 : Int_t ip=0, itheta=0, iphi=0;
312 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
313 0 : Int_t nSplitP, nSplitTheta;
314 0 : GetSplit(ip,itheta,nSplitP,nSplitTheta);
315 :
316 0 : Float_t dp = p - fPmin;
317 0 : Int_t ibinp = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP);
318 0 : Float_t dtheta = theta - fThetamin;
319 0 : Int_t ibintheta = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta);
320 0 : Float_t eff = fCurrentEntry[ip][itheta][iphi]->GetEff(ibinp,ibintheta);
321 0 : return eff;
322 0 : }
323 :
324 : Float_t AliMUONFastTracking::Acceptance(Float_t p, Float_t theta,
325 : Float_t phi, Int_t charge){
326 : //
327 : // gets the geometrical acceptance
328 : //
329 0 : if (theta<fThetamin || theta>fThetamax) return 0;
330 :
331 0 : Int_t ip=0, itheta=0, iphi=0;
332 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
333 0 : Int_t nSplitP, nSplitTheta;
334 0 : GetSplit(ip,itheta,nSplitP,nSplitTheta);
335 : // central value and corrections with spline
336 :
337 0 : Float_t dp = p - fPmin;
338 0 : Int_t ibinp = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP);
339 0 : Float_t dtheta = theta - fThetamin;
340 0 : Int_t ibintheta = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta);
341 0 : Float_t acc = fCurrentEntry[ip][itheta][iphi]->GetAcc(ibinp,ibintheta);
342 : return acc;
343 0 : }
344 :
345 : Float_t AliMUONFastTracking::MeanP(Float_t p, Float_t theta,
346 : Float_t phi, Int_t charge) const
347 : {
348 : //
349 : // gets the mean value of the prec-pgen distribution
350 : //
351 0 : Int_t ip=0, itheta=0, iphi=0;
352 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
353 0 : return fCurrentEntry[ip][itheta][iphi]->GetMeanp();
354 0 : }
355 :
356 : Float_t AliMUONFastTracking::SigmaP(Float_t p, Float_t theta,
357 : Float_t phi, Int_t charge) const
358 : {
359 : //
360 : // gets the width of the prec-pgen distribution
361 : //
362 0 : Int_t ip=0, itheta=0, iphi=0;
363 : Int_t index;
364 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
365 : // central value and corrections with spline
366 0 : Float_t sigmap = fCurrentEntry[ip][itheta][iphi]->GetSigmap();
367 0 : if (!fSpline) return sigmap;
368 : // corrections vs p, theta, phi
369 0 : index = iphi + fNbinphi * itheta;
370 0 : Double_t xmin,ymin,xmax,ymax;
371 0 : Float_t frac1 = fSplineSigmap[index][0]->Eval(p)/sigmap;
372 :
373 0 : if (p>fPmax-fDeltaP/2.) {
374 0 : Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigmap();
375 0 : Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigmap();
376 0 : Float_t s3 = fCurrentEntry[fNbinp-3][itheta][iphi]->GetSigmap();
377 0 : Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
378 0 : Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
379 0 : Float_t p3 = fDeltaP * (fNbinp - 3 + 0.5) + fPmin;
380 0 : Float_t p12 = p1 * p1, p22 = p2 * p2, p32 = p3 * p3;
381 0 : Float_t d = p12*p2 + p1*p32 + p22*p3 - p32*p2 - p3*p12 - p22*p1;
382 0 : Float_t a = (s1*p2 + p1*s3 + s2*p3 - s3*p2 - p3*s1 - s2*p1) / d;
383 0 : Float_t b = (p12*s2 + s1*p32 + p22*s3 - p32*s2 - s3*p12 - p22*s1)/d;
384 0 : Float_t c = (p12*p2*s3 + p1*p32*s2 + p22*p3*s1
385 0 : - p32*p2*s1 - p3*p12*s2 - p22*p1*s3) / d;
386 0 : Float_t sigma = a * p * p + b * p + c;
387 0 : frac1 = sigma/sigmap;
388 0 : }
389 0 : index = iphi + fNbinphi * ip;
390 0 : fSplineEff[index][1]->GetKnot(0,xmin,ymin);
391 0 : fSplineEff[index][1]->GetKnot(9,xmax,ymax);
392 0 : if (theta>xmax) theta = xmax;
393 0 : Float_t frac2 = fSplineSigmap[index][1]->Eval(theta)/sigmap;
394 0 : index = itheta + fNbintheta * ip;
395 0 : fSplineEff[index][2]->GetKnot(0,xmin,ymin);
396 0 : fSplineEff[index][2]->GetKnot(9,xmax,ymax);
397 0 : if (phi>xmax) phi = xmax;
398 0 : Float_t frac3 = fSplineSigmap[index][2]->Eval(phi)/sigmap;
399 0 : Float_t sigmatot = sigmap * frac1 * frac2 * frac3;
400 0 : if (sigmatot<0) sigmatot = sigmap;
401 : return sigmatot;
402 0 : }
403 :
404 : Float_t AliMUONFastTracking::Sigma1P(Float_t p, Float_t theta,
405 : Float_t phi, Int_t charge) const
406 : {
407 : //
408 : // gets the width correction of the prec-pgen distribution (see FitP)
409 : //
410 0 : Int_t ip=0, itheta=0, iphi=0;
411 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
412 0 : if (p>fPmax) {
413 : // linear extrapolation of sigmap for p out of range
414 0 : Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigma1p();
415 0 : Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigma1p();
416 0 : Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
417 0 : Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
418 0 : Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
419 : return sigma;
420 : }
421 0 : else return fCurrentEntry[ip][itheta][iphi]->GetSigma1p();
422 0 : }
423 :
424 : Float_t AliMUONFastTracking::NormG2(Float_t p, Float_t theta,
425 : Float_t phi, Int_t charge) const
426 : {
427 : //
428 : // gets the relative normalization of the background
429 : // (gaussian) component in the prec-pgen distribution
430 : //
431 0 : Int_t ip=0, itheta=0, iphi=0;
432 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
433 0 : if (p>fPmax) {
434 : // linear extrapolation of sigmap for p out of range
435 0 : Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetNormG2();
436 0 : Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetNormG2();
437 0 : Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
438 0 : Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
439 0 : Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
440 : return norm;
441 : }
442 0 : else return fCurrentEntry[ip][itheta][iphi]->GetNormG2();
443 0 : }
444 :
445 : Float_t AliMUONFastTracking::MeanG2(Float_t p, Float_t theta,
446 : Float_t phi, Int_t charge) const
447 : {
448 : //
449 : // gets the mean value of the background
450 : // (gaussian) component in the prec-pgen distribution
451 : //
452 0 : Int_t ip=0, itheta=0, iphi=0;
453 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
454 0 : if (p>fPmax) {
455 : // linear extrapolation of sigmap for p out of range
456 0 : Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetMeanG2();
457 0 : Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetMeanG2();
458 0 : Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
459 0 : Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
460 0 : Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
461 : return norm;
462 : }
463 0 : else return fCurrentEntry[ip][itheta][iphi]->GetMeanG2();
464 0 : }
465 :
466 : Float_t AliMUONFastTracking::SigmaG2(Float_t p, Float_t theta,
467 : Float_t phi, Int_t charge) const
468 : {
469 : //
470 : // gets the width of the background
471 : // (gaussian) component in the prec-pgen distribution
472 : //
473 0 : Int_t ip=0, itheta=0, iphi=0;
474 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
475 0 : if (p>fPmax) {
476 : // linear extrapolation of sigmap for p out of range
477 0 : Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigmaG2();
478 0 : Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigmaG2();
479 0 : Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
480 0 : Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
481 0 : Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
482 : return sigma;
483 : }
484 0 : else return fCurrentEntry[ip][itheta][iphi]->GetSigmaG2();
485 0 : }
486 :
487 :
488 : Float_t AliMUONFastTracking::MeanTheta(Float_t p, Float_t theta,
489 : Float_t phi, Int_t charge) const
490 : {
491 : //
492 : // gets the mean value of the thetarec-thetagen distribution
493 : //
494 0 : Int_t ip=0, itheta=0, iphi=0;
495 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
496 0 : return fCurrentEntry[ip][itheta][iphi]->GetMeantheta();
497 0 : }
498 :
499 : Float_t AliMUONFastTracking::SigmaTheta(Float_t p, Float_t theta,
500 : Float_t phi, Int_t charge) const
501 : {
502 : //
503 : // gets the width of the thetarec-thetagen distribution
504 : //
505 0 : Int_t ip=0, itheta=0, iphi=0;
506 : Int_t index;
507 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
508 : // central value and corrections with spline
509 0 : Float_t sigmatheta = fCurrentEntry[ip][itheta][iphi]->GetSigmatheta();
510 0 : if (!fSpline) return sigmatheta;
511 : // corrections vs p, theta, phi
512 0 : index = iphi + fNbinphi * itheta;
513 0 : Double_t xmin,ymin,xmax,ymax;
514 0 : Float_t frac1 = fSplineSigmatheta[index][0]->Eval(p)/sigmatheta;
515 0 : if (p>fPmax-fDeltaP/2.) {
516 : // linear extrapolation of sigmap for p out of range
517 0 : Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigmatheta();
518 0 : Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigmatheta();
519 0 : Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
520 0 : Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
521 0 : Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
522 0 : frac1=sigma/sigmatheta;
523 0 : }
524 0 : index = iphi + fNbinphi * ip;
525 0 : fSplineEff[index][1]->GetKnot(0,xmin,ymin);
526 0 : fSplineEff[index][1]->GetKnot(9,xmax,ymax);
527 0 : if (theta>xmax) theta = xmax;
528 0 : Float_t frac2 = fSplineSigmatheta[index][1]->Eval(theta)/sigmatheta;
529 0 : index = itheta + fNbintheta * ip;
530 0 : fSplineEff[index][2]->GetKnot(0,xmin,ymin);
531 0 : fSplineEff[index][2]->GetKnot(9,xmax,ymax);
532 0 : if (phi>xmax) phi = xmax;
533 0 : Float_t frac3 = fSplineSigmatheta[index][2]->Eval(phi)/sigmatheta;
534 0 : return sigmatheta * frac1 * frac2 * frac3;
535 0 : }
536 :
537 :
538 : Float_t AliMUONFastTracking::MeanPhi(Float_t p, Float_t theta,
539 : Float_t phi, Int_t charge) const
540 : {
541 : //
542 : // gets the mean value of the phirec-phigen distribution
543 : //
544 0 : Int_t ip=0, itheta=0, iphi=0;
545 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
546 0 : return fCurrentEntry[ip][itheta][iphi]->GetMeanphi();
547 0 : }
548 :
549 : Float_t AliMUONFastTracking::SigmaPhi(Float_t p, Float_t theta,
550 : Float_t phi, Int_t charge){
551 : //
552 : // gets the width of the phirec-phigen distribution
553 : //
554 0 : Int_t ip=0, itheta=0, iphi=0;
555 : Int_t index;
556 0 : GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
557 : // central value and corrections with spline
558 0 : Float_t sigmaphi = fCurrentEntry[ip][itheta][iphi]->GetSigmaphi();
559 0 : if (!fSpline) return sigmaphi;
560 : // corrections vs p, theta, phi
561 0 : index = iphi + fNbinphi * itheta;
562 0 : Float_t frac1 = fSplineSigmaphi[index][0]->Eval(p)/sigmaphi;
563 0 : Double_t xmin,ymin,xmax,ymax;
564 0 : if (p>fPmax-fDeltaP/2.) {
565 0 : Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigmaphi();
566 0 : Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigmaphi();
567 0 : Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
568 0 : Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
569 0 : Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
570 0 : frac1 = sigma/sigmaphi;
571 0 : }
572 :
573 0 : index = iphi + fNbinphi * ip;
574 0 : fSplineEff[index][1]->GetKnot(0,xmin,ymin);
575 0 : fSplineEff[index][1]->GetKnot(9,xmax,ymax);
576 0 : if (theta>xmax) theta = xmax;
577 0 : Float_t frac2 = fSplineSigmaphi[index][1]->Eval(theta)/sigmaphi;
578 0 : index = itheta + fNbintheta * ip;
579 0 : fSplineEff[index][2]->GetKnot(0,xmin,ymin);
580 0 : fSplineEff[index][2]->GetKnot(9,xmax,ymax);
581 0 : if (phi>xmax) phi = xmax;
582 0 : Float_t frac3 = fSplineSigmaphi[index][2]->Eval(phi)/sigmaphi;
583 0 : return sigmaphi * frac1 * frac2 * frac3;
584 0 : }
585 :
586 : void AliMUONFastTracking::SetSpline(){
587 : //
588 : // sets the spline functions for a smooth behaviour of the parameters
589 : // when going from one cell to another
590 : //
591 0 : printf ("Setting spline functions...");
592 0 : char splname[40];
593 0 : Double_t x[20][3];
594 0 : Double_t x2[50][3];
595 0 : Int_t nbins[3] = {fNbinp, fNbintheta, fNbinphi};
596 0 : Double_t xspl[20],yeff[50],ysigmap[20],ysigma1p[20];
597 0 : Double_t yacc[50], ysigmatheta[20],ysigmaphi[20];
598 0 : Double_t xsp2[50];
599 : // let's calculate the x axis for p, theta, phi
600 :
601 : Int_t i, ispline, ivar;
602 0 : for (i=0; i< fNbinp; i++) x[i][0] = fPmin + fDeltaP * (i + 0.5);
603 0 : for (i=0; i< fNbintheta; i++) x[i][1] = fThetamin + fDeltaTheta * (i + 0.5);
604 0 : for (i=0; i< fNbinphi; i++) x[i][2] = fPhimin + fDeltaPhi * (i + 0.5);
605 :
606 0 : for (i=0; i< 5 * fNbinp; i++) x2[i][0] = fPmin + fDeltaP * (i + 0.5)/5.;
607 0 : for (i=0; i< 5 * fNbintheta; i++) x2[i][1] = fThetamin + fDeltaTheta * (i + 0.5)/5.;
608 0 : for (i=0; i< 5 * fNbinphi; i++) x2[i][2] = fPhimin + fDeltaPhi * (i + 0.5)/5.;
609 :
610 : // splines in p
611 : ivar = 0;
612 0 : for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
613 0 : for (i=0; i<5 * nbins[ivar]; i++) xsp2[i] = x2[i][ivar];
614 : ispline=0;
615 0 : for (Int_t itheta=0; itheta< fNbintheta; itheta++){
616 0 : for (Int_t iphi=0; iphi< fNbinphi; iphi++){
617 0 : for (Int_t ip=0; ip<fNbinp; ip++) {
618 0 : ysigmap[ip] = fCurrentEntry[ip][itheta][iphi]->GetSigmap();
619 0 : ysigma1p[ip] = fCurrentEntry[ip][itheta][iphi]->GetSigma1p();
620 0 : ysigmatheta[ip] = fCurrentEntry[ip][itheta][iphi]->GetSigmatheta();
621 0 : ysigmaphi[ip] = fCurrentEntry[ip][itheta][iphi]->GetSigmaphi();
622 : }
623 0 : if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
624 0 : snprintf (splname, 40, "fSplineEff[%d][%d]",ispline,ivar);
625 0 : fSplineEff[ispline][ivar] = new TSpline3(splname,xsp2,yeff,5 * nbins[ivar]);
626 0 : snprintf (splname, 40, "fSplineAcc[%d][%d]",ispline,ivar);
627 0 : fSplineAcc[ispline][ivar] = new TSpline3(splname,xsp2,yacc,5 * nbins[ivar]);
628 0 : snprintf (splname, 40, "fSplineSigmap[%d][%d]",ispline,ivar);
629 0 : fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
630 0 : snprintf (splname, 40, "fSplineSigma1p[%d][%d]",ispline,ivar);
631 0 : fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
632 0 : snprintf (splname, 40, "fSplineSigmatheta[%d][%d]",ispline,ivar);
633 0 : fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
634 0 : snprintf (splname, 40, "fSplineSigmaphi[%d][%d]",ispline,ivar);
635 0 : fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
636 0 : ispline++;
637 : }
638 : }
639 :
640 : ivar = 1;
641 0 : for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
642 : ispline=0;
643 0 : for (Int_t ip=0; ip<fNbinp; ip++) {
644 0 : for (Int_t iphi=0; iphi< fNbinphi; iphi++){
645 0 : for (Int_t itheta=0; itheta< fNbintheta; itheta++){
646 : // for efficiency and acceptance let's take the central value
647 0 : ysigmap[itheta] = fCurrentEntry[ip][itheta][iphi]->GetSigmap();
648 0 : ysigma1p[itheta] = fCurrentEntry[ip][itheta][iphi]->GetSigma1p();
649 0 : ysigmatheta[itheta] = fCurrentEntry[ip][itheta][iphi]->GetSigmatheta();
650 0 : ysigmaphi[itheta] = fCurrentEntry[ip][itheta][iphi]->GetSigmaphi();
651 : }
652 0 : if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
653 0 : snprintf (splname, 40, "fSplineEff[%d][%d]",ispline,ivar);
654 0 : fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
655 0 : snprintf (splname, 40, "fSplineAcc[%d][%d]",ispline,ivar);
656 0 : fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
657 0 : snprintf (splname, 40, "fSplineSigmap[%d][%d]",ispline,ivar);
658 0 : fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
659 0 : snprintf (splname, 40, "fSplineSigma1p[%d][%d]",ispline,ivar);
660 0 : fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
661 0 : snprintf (splname, 40, "fSplineSigmatheta[%d][%d]",ispline,ivar);
662 0 : fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
663 0 : snprintf (splname, 40, "fSplineSigmaphi[%d][%d]",ispline,ivar);
664 0 : fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
665 0 : ispline++;
666 : }
667 : }
668 :
669 : ivar = 2;
670 0 : for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
671 : ispline=0;
672 0 : for (Int_t ip=0; ip<fNbinp; ip++) {
673 0 : for (Int_t itheta=0; itheta< fNbintheta; itheta++){
674 0 : for (Int_t iphi=0; iphi< fNbinphi; iphi++){
675 : // for efficiency and acceptance let's take the central value
676 0 : ysigmap[iphi] = fCurrentEntry[ip][itheta][iphi]->GetSigmap();
677 0 : ysigma1p[iphi] = fCurrentEntry[ip][itheta][iphi]->GetSigma1p();
678 0 : ysigmatheta[iphi] = fCurrentEntry[ip][itheta][iphi]->GetSigmatheta();
679 0 : ysigmaphi[iphi] = fCurrentEntry[ip][itheta][iphi]->GetSigmaphi();
680 : }
681 0 : if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
682 0 : snprintf (splname, 40, "fSplineEff[%d][%d]",ispline,ivar);
683 0 : fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
684 0 : snprintf (splname, 40, "fSplineAcc[%d][%d]",ispline,ivar);
685 0 : fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
686 0 : snprintf (splname, 40, "fSplineSigmap[%d][%d]",ispline,ivar);
687 0 : fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
688 0 : snprintf (splname, 40, "fSplineSigma1p[%d][%d]",ispline,ivar);
689 0 : fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
690 0 : snprintf (splname, 40, "fSplineSigmatheta[%d][%d]",ispline,ivar);
691 0 : fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
692 0 : snprintf (splname, 40, "fSplineSigmaphi[%d][%d]",ispline,ivar);
693 0 : fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
694 0 : ispline++;
695 : }
696 : }
697 0 : printf ("...done\n");
698 0 : }
699 :
700 : void AliMUONFastTracking::SetBackground(Float_t bkg){
701 : //
702 : // linear interpolation of the parameters in the LUT between 2 values where
703 : // the background has been actually calculated
704 : //
705 0 : if (bkg>2) printf ("WARNING: unsafe extrapolation!\n");
706 0 : fBkg = bkg;
707 :
708 0 : Float_t bkgLevel[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
709 : Int_t ibkg;
710 0 : for (ibkg=0; ibkg<4; ibkg++) if ( bkg < bkgLevel[ibkg]) break;
711 0 : if (ibkg == 4) ibkg--;
712 0 : if (ibkg == 0) ibkg++;
713 :
714 0 : Float_t x0 = bkgLevel[ibkg-1];
715 0 : Float_t x1 = bkgLevel[ibkg];
716 0 : Float_t x = (bkg - x0) / (x1 - x0);
717 :
718 : Float_t y0, y1;
719 :
720 0 : for (Int_t ip=0; ip< fNbinp; ip++){
721 0 : for (Int_t itheta=0; itheta< fNbintheta; itheta++){
722 0 : for (Int_t iphi=0; iphi< fNbinphi; iphi++){
723 0 : fCurrentEntry[ip][itheta][iphi]->SetP(fEntry[ip][itheta][iphi][ibkg]->GetP());
724 0 : fCurrentEntry[ip][itheta][iphi]->SetTheta(fEntry[ip][itheta][iphi][ibkg]->GetTheta());
725 0 : fCurrentEntry[ip][itheta][iphi]->SetPhi(fEntry[ip][itheta][iphi][ibkg]->GetPhi());
726 0 : fCurrentEntry[ip][itheta][iphi]->SetChi2p(-1);
727 0 : fCurrentEntry[ip][itheta][iphi]->SetChi2theta(-1);
728 0 : fCurrentEntry[ip][itheta][iphi]->SetChi2phi(-1);
729 :
730 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetMeanp();
731 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetMeanp();
732 0 : fCurrentEntry[ip][itheta][iphi] ->SetMeanp((y1 - y0) * x + y0);
733 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetMeantheta();
734 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetMeantheta();
735 0 : fCurrentEntry[ip][itheta][iphi] ->SetMeantheta((y1 - y0) * x +y0);
736 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetMeanphi();
737 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetMeanphi();
738 0 : fCurrentEntry[ip][itheta][iphi] ->SetMeanphi((y1 - y0) * x + y0);
739 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetSigmap();
740 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetSigmap();
741 0 : fCurrentEntry[ip][itheta][iphi] ->SetSigmap((y1 - y0) * x + y0);
742 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetSigmatheta();
743 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetSigmatheta();
744 0 : fCurrentEntry[ip][itheta][iphi] ->SetSigmatheta((y1 - y0) * x+y0);
745 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetSigmaphi();
746 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetSigmaphi();
747 0 : fCurrentEntry[ip][itheta][iphi] ->SetSigmaphi((y1 - y0) * x + y0);
748 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetSigma1p();
749 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetSigma1p();
750 0 : fCurrentEntry[ip][itheta][iphi] ->SetSigma1p((y1 - y0) * x + y0);
751 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetNormG2();
752 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetNormG2();
753 0 : fCurrentEntry[ip][itheta][iphi] ->SetNormG2((y1 - y0) * x + y0);
754 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetMeanG2();
755 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetMeanG2();
756 0 : fCurrentEntry[ip][itheta][iphi] ->SetMeanG2((y1 - y0) * x + y0);
757 :
758 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetSigmaG2();
759 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetSigmaG2();
760 0 : fCurrentEntry[ip][itheta][iphi] ->SetSigmaG2((y1 - y0) * x + y0);
761 0 : for (Int_t i=0; i<kSplitP; i++) {
762 0 : for (Int_t j=0; j<kSplitTheta; j++) {
763 0 : fCurrentEntry[ip][itheta][iphi]->SetAcc(i,j,fEntry[ip][itheta][iphi][ibkg]->GetAcc(i,j));
764 0 : y0 = fEntry[ip][itheta][iphi][ibkg-1]->GetEff(i,j);
765 0 : y1 = fEntry[ip][itheta][iphi][ibkg]->GetEff(i,j);
766 0 : fCurrentEntry[ip][itheta][iphi]->SetEff(i,j, (y1 - y0) * x + y0);
767 : }
768 : }
769 : }
770 : }
771 : }
772 0 : SetSpline();
773 0 : }
774 :
775 : TF1* AliMUONFastTracking::GetFitP(Int_t ip,Int_t itheta,Int_t iphi) {
776 : // gets the correct prec-pgen distribution for a given LUT cell
777 0 : if (!fFitp[ip][itheta][iphi]) {
778 0 : char name[256];
779 0 : snprintf(name, 256, "fit_%d_%d_%d", ip, itheta, iphi);
780 0 : fFitp[ip][itheta][iphi] = new TF1(name ,FitP,-20.,20.,6);
781 0 : fFitp[ip][itheta][iphi]->SetNpx(500);
782 0 : fFitp[ip][itheta][iphi]->SetParameters(0.,0.,0.,0.,0.,0.);
783 0 : }
784 0 : return fFitp[ip][itheta][iphi];
785 0 : }
786 :
787 : AliMUONFastTracking& AliMUONFastTracking::operator=(const AliMUONFastTracking& rhs)
788 : {
789 : // Assignment operator
790 0 : rhs.Copy(*this);
791 0 : return *this;
792 : }
793 :
794 : void AliMUONFastTracking::Copy(TObject&) const
795 : {
796 : //
797 : // Copy
798 : //
799 0 : Fatal("Copy","Not implemented!\n");
800 0 : }
801 :
802 :
803 :
804 :
805 :
806 :
807 :
808 :
|
