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