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