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