1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.12 2006/08/28 10:31:17 morsch
19 Coding rule violations corrected (A. de Falco)
21 Revision 1.10 2006/01/27 09:51:37 morsch
22 Some small corrections to avoid infinite loops at high momenta.
25 Revision 1.9 2004/02/03 16:45:20 morsch
26 Unique name for functions (TF1).
28 Revision 1.8 2004/01/27 18:02:38 hristov
29 Removing some warning (Sun)
31 Revision 1.7 2003/11/13 14:21:57 morsch
32 Coding Rule violation corrections.
34 Revision 1.6 2003/08/12 15:16:25 morsch
35 Saver initialisation of fFitp array. (Lenaic COUEDEL)
37 Revision 1.5 2003/08/05 16:14:20 morsch
38 Some problems with too big fluctuations corrected. (A. de Falco)
40 Revision 1.2 2003/01/08 10:29:33 morsch
41 Path to data file changed.
43 Revision 1.1 2003/01/06 10:13:33 morsch
48 //-------------------------------------------------------------------------
49 // Class AliMUONFastTracking
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 //-------------------------------------------------------------------------
59 #include "AliMUONFastTracking.h"
60 #include "AliMUONFastTrackingEntry.h"
69 #include <Riostream.h>
71 ClassImp(AliMUONFastTracking)
74 AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
76 static 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);
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));
87 Double_t value = fasymm + par[3] * fgauss;
88 return TMath::Abs(value);
91 AliMUONFastTracking::AliMUONFastTracking(const AliMUONFastTracking & ft):
96 fDeltaP((fPmax-fPmin)/fNbinp),
100 fDeltaTheta((fThetamax-fThetamin)/fNbintheta),
104 fDeltaPhi((fPhimax-fPhimin)/fNbinphi),
115 AliMUONFastTracking* AliMUONFastTracking::Instance()
117 // Set random number generator
118 if (fgMUONFastTracking) {
119 return fgMUONFastTracking;
121 fgMUONFastTracking = new AliMUONFastTracking();
122 return fgMUONFastTracking;
126 AliMUONFastTracking::AliMUONFastTracking():
130 fDeltaP((fPmax-fPmin)/fNbinp),
134 fDeltaTheta((fThetamax-fThetamin)/fNbintheta),
138 fDeltaPhi((fPhimax-fPhimin)/fNbinphi),
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;
156 void AliMUONFastTracking::Init(Float_t bkg)
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;
173 if (fClusterFinder==kOld) sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root");
174 else sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT-AZ.root");
176 TFile *file = new TFile(filename);
183 void AliMUONFastTracking::ReadLUT(TFile* file)
186 // read the lookup tables from file
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];
193 printf ("Reading parameters from LUT file %s...\n",file->GetName());
195 const Float_t kBkg[4] = {0, 0.5, 1, 2};
196 for (Int_t ibkg=0; ibkg<4; ibkg++) {
197 sprintf (tag,"BKG%g",kBkg[ibkg]);
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);
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");
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);
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));
245 for (Int_t i=0; i<kSplitP; i++) {
246 for (Int_t j=0; j<kSplitTheta; j++) {
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));
256 TGraph *graph = new TGraph(3);
257 TF1 *f = new TF1("f","[0]+[1]*x");
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++){
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());
266 fEntry[ip][itheta][iphi][0]->SetSigmaG2(f->Eval(0));
272 printf ("parameters read. \n");
275 void AliMUONFastTracking::GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
276 Int_t &nbintheta, Float_t &thetamin,
278 Int_t &nbinphi, Float_t &phimin, Float_t &phimax) const
281 // gets the binning for the discrete parametrizations in the lookup table
286 nbintheta = fNbintheta;
287 thetamin = fThetamin;
288 thetamax = fThetamax;
295 void AliMUONFastTracking::GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi,
296 Int_t charge, Int_t &ip, Int_t &itheta,
300 // gets the id of the cells in the LUT for a given (p,theta,phi, charge)
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);
309 if (ip>= fNbinp) ip = fNbinp-1;
310 if (itheta< 0) itheta = 0;
311 if (itheta>= fNbintheta) itheta = fNbintheta-1;
313 if (iphi< 0) iphi = 0;
314 if (iphi>= fNbinphi) iphi = fNbinphi-1;
317 void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta,
318 Int_t &nSplitP, Int_t &nSplitTheta) const
321 // the first cell is splitted in more bins for theta and momentum
322 // parameterizations. Get the number of divisions for the splitted bins
324 if (ip==0) nSplitP = 5;
326 if (itheta==0) nSplitTheta = 3;
327 else nSplitTheta = 1;
330 Float_t AliMUONFastTracking::Efficiency(Float_t p, Float_t theta,
331 Float_t phi, Int_t charge){
333 // gets the tracking efficiency
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);
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);
344 Float_t eff = fCurrentEntry[ip][itheta][iphi]->GetEff(ibinp,ibintheta);
348 Float_t AliMUONFastTracking::Acceptance(Float_t p, Float_t theta,
349 Float_t phi, Int_t charge){
351 // gets the geometrical acceptance
353 if (theta<fThetamin || theta>fThetamax) return 0;
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
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);
365 Float_t acc = fCurrentEntry[ip][itheta][iphi]->GetAcc(ibinp,ibintheta);
369 Float_t AliMUONFastTracking::MeanP(Float_t p, Float_t theta,
370 Float_t phi, Int_t charge) const
373 // gets the mean value of the prec-pgen distribution
375 Int_t ip=0, itheta=0, iphi=0;
376 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
377 return fCurrentEntry[ip][itheta][iphi]->GetMeanp();
380 Float_t AliMUONFastTracking::SigmaP(Float_t p, Float_t theta,
381 Float_t phi, Int_t charge) const
384 // gets the width of the prec-pgen distribution
386 Int_t ip=0, itheta=0, iphi=0;
388 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
389 // central value and corrections with spline
390 Float_t sigmap = fCurrentEntry[ip][itheta][iphi]->GetSigmap();
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;
397 if (p>fPmax-fDeltaP/2.) {
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();
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;
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;
428 Float_t AliMUONFastTracking::Sigma1P(Float_t p, Float_t theta,
429 Float_t phi, Int_t charge) const
432 // gets the width correction of the prec-pgen distribution (see FitP)
434 Int_t ip=0, itheta=0, iphi=0;
435 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
437 // linear extrapolation of sigmap for p out of range
438 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigma1p();
439 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigma1p();
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) );
445 else return fCurrentEntry[ip][itheta][iphi]->GetSigma1p();
448 Float_t AliMUONFastTracking::NormG2(Float_t p, Float_t theta,
449 Float_t phi, Int_t charge) const
452 // gets the relative normalization of the background
453 // (gaussian) component in the prec-pgen distribution
455 Int_t ip=0, itheta=0, iphi=0;
456 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
458 // linear extrapolation of sigmap for p out of range
459 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetNormG2();
460 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetNormG2();
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) );
466 else return fCurrentEntry[ip][itheta][iphi]->GetNormG2();
469 Float_t AliMUONFastTracking::MeanG2(Float_t p, Float_t theta,
470 Float_t phi, Int_t charge) const
473 // gets the mean value of the background
474 // (gaussian) component in the prec-pgen distribution
476 Int_t ip=0, itheta=0, iphi=0;
477 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
479 // linear extrapolation of sigmap for p out of range
480 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetMeanG2();
481 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetMeanG2();
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) );
487 else return fCurrentEntry[ip][itheta][iphi]->GetMeanG2();
490 Float_t AliMUONFastTracking::SigmaG2(Float_t p, Float_t theta,
491 Float_t phi, Int_t charge) const
494 // gets the width of the background
495 // (gaussian) component in the prec-pgen distribution
497 Int_t ip=0, itheta=0, iphi=0;
498 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
500 // linear extrapolation of sigmap for p out of range
501 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigmaG2();
502 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigmaG2();
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) );
508 else return fCurrentEntry[ip][itheta][iphi]->GetSigmaG2();
512 Float_t AliMUONFastTracking::MeanTheta(Float_t p, Float_t theta,
513 Float_t phi, Int_t charge) const
516 // gets the mean value of the thetarec-thetagen distribution
518 Int_t ip=0, itheta=0, iphi=0;
519 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
520 return fCurrentEntry[ip][itheta][iphi]->GetMeantheta();
523 Float_t AliMUONFastTracking::SigmaTheta(Float_t p, Float_t theta,
524 Float_t phi, Int_t charge) const
527 // gets the width of the thetarec-thetagen distribution
529 Int_t ip=0, itheta=0, iphi=0;
531 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
532 // central value and corrections with spline
533 Float_t sigmatheta = fCurrentEntry[ip][itheta][iphi]->GetSigmatheta();
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
541 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigmatheta();
542 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigmatheta();
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;
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;
562 Float_t AliMUONFastTracking::MeanPhi(Float_t p, Float_t theta,
563 Float_t phi, Int_t charge) const
566 // gets the mean value of the phirec-phigen distribution
568 Int_t ip=0, itheta=0, iphi=0;
569 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
570 return fCurrentEntry[ip][itheta][iphi]->GetMeanphi();
573 Float_t AliMUONFastTracking::SigmaPhi(Float_t p, Float_t theta,
574 Float_t phi, Int_t charge){
576 // gets the width of the phirec-phigen distribution
578 Int_t ip=0, itheta=0, iphi=0;
580 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
581 // central value and corrections with spline
582 Float_t sigmaphi = fCurrentEntry[ip][itheta][iphi]->GetSigmaphi();
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.) {
589 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->GetSigmaphi();
590 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->GetSigmaphi();
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;
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;
610 void AliMUONFastTracking::SetSpline(){
612 // sets the spline functions for a smooth behaviour of the parameters
613 // when going from one cell to another
615 printf ("Setting spline functions...");
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];
623 // let's calculate the x axis for p, theta, phi
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);
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.;
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];
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++) {
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();
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]);
665 for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
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
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();
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]);
694 for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
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
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();
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]);
721 printf ("...done\n");
724 void AliMUONFastTracking::SetBackground(Float_t bkg){
726 // linear interpolation of the parameters in the LUT between 2 values where
727 // the background has been actually calculated
729 if (bkg>2) printf ("WARNING: unsafe extrapolation!\n");
732 Float_t bkgLevel[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
734 for (ibkg=0; ibkg<4; ibkg++) if ( bkg < bkgLevel[ibkg]) break;
735 if (ibkg == 4) ibkg--;
736 if (ibkg == 0) ibkg++;
738 Float_t x0 = bkgLevel[ibkg-1];
739 Float_t x1 = bkgLevel[ibkg];
740 Float_t x = (bkg - x0) / (x1 - x0);
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++){
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);
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);
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);
785 for (Int_t i=0; i<kSplitP; i++) {
786 for (Int_t j=0; j<kSplitTheta; j++) {
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);
799 TF1* AliMUONFastTracking::GetFitP(Int_t ip,Int_t itheta,Int_t iphi) {
800 // gets the correct prec-pgen distribution for a given LUT cell
801 if (!fFitp[ip][itheta][iphi]) {
803 sprintf(name, "fit_%d_%d_%d", ip, itheta, iphi);
804 fFitp[ip][itheta][iphi] = new TF1(name ,FitP,-20.,20.,6);
805 fFitp[ip][itheta][iphi]->SetNpx(500);
806 fFitp[ip][itheta][iphi]->SetParameters(0.,0.,0.,0.,0.,0.);
808 return fFitp[ip][itheta][iphi];
811 AliMUONFastTracking& AliMUONFastTracking::operator=(const AliMUONFastTracking& rhs)
813 // Assignment operator
818 void AliMUONFastTracking::Copy(TObject&) const
823 Fatal("Copy","Not implemented!\n");