Bug fix (Alpha)
[u/mrichter/AliRoot.git] / FASTSIM / AliFastMuonTrackingRes.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.4  2003/08/05 16:14:20  morsch
19 Some problems with too big fluctuations corrected. (A. de Falco)
20
21 Revision 1.1  2003/01/06 10:13:09  morsch
22 First commit.
23
24 */
25
26 #include "AliFastMuonTrackingRes.h"
27 #include "AliMUONFastTracking.h"
28 #include <TRandom.h>
29 #include <TF1.h>
30
31 ClassImp(AliFastMuonTrackingRes)
32
33
34 AliFastMuonTrackingRes::AliFastMuonTrackingRes() :
35     AliFastResponse("Resolution", "Muon Tracking Resolution")
36 {
37     SetBackground();
38 }
39
40 void AliFastMuonTrackingRes::Init()
41 {
42     fFastTracking = AliMUONFastTracking::Instance();
43     fFastTracking->Init(fBackground);
44 }
45
46
47
48 void AliFastMuonTrackingRes::Evaluate(Float_t   p,  Float_t  theta , Float_t   phi,
49                                          Float_t& pS,  Float_t& thetaS, Float_t&  phiS)
50 {
51
52     Double_t meanp    = fFastTracking->MeanP  (p, theta, phi, Int_t(fCharge));
53     Double_t sigmap   = fFastTracking->SigmaP (p, theta, phi, Int_t(fCharge));
54     Double_t sigma1p  = fFastTracking->Sigma1P(p, theta, phi, Int_t(fCharge));
55     Double_t normg2   = fFastTracking->NormG2 (p, theta, phi, Int_t(fCharge));
56     Double_t meang2   = fFastTracking->MeanG2 (p, theta, phi, Int_t(fCharge));
57     Double_t sigmag2  = fFastTracking->SigmaG2(p, theta, phi, Int_t(fCharge));
58     
59     Int_t ip,itheta,iphi;
60     fFastTracking->GetIpIthetaIphi(p, theta, phi, Int_t(fCharge), ip, itheta, iphi);
61     TF1* fitp = fFastTracking->GetFitP(ip,itheta,iphi);
62     
63     Float_t curmeanp   = fitp->GetParameter(0); 
64     Float_t cursigmap  = fitp->GetParameter(1); 
65     Float_t cursigma1p = fitp->GetParameter(2); 
66     Float_t curnormg2  = fitp->GetParameter(3); 
67     Float_t curmeang2  = fitp->GetParameter(4); 
68     Float_t cursigmag2 = fitp->GetParameter(5); 
69     if (curmeanp != meanp || cursigmap != sigmap || cursigma1p != sigma1p || 
70         curnormg2 != normg2 || curmeang2 != meang2 || cursigmag2 != sigmag2){ 
71       printf ("Setting new parameters for ip=%d itheta=%d iphi=%d\n",ip,itheta,iphi); 
72       fitp->SetParameters(meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
73     }
74   
75     Double_t meantheta  = fFastTracking->MeanTheta (p, theta, phi, Int_t(fCharge));
76     Double_t sigmatheta = fFastTracking->SigmaTheta(p, theta, phi, Int_t(fCharge));
77     Double_t meanphi    = fFastTracking->MeanPhi   (p, theta, phi, Int_t(fCharge));
78     Double_t sigmaphi   = fFastTracking->SigmaPhi  (p, theta, phi, Int_t(fCharge));
79     
80     if (sigmatheta<0 || sigmaphi<0) 
81       printf ("bin %d %d %d sigmatheta = %f, sigmaphi = %f\n",
82               ip,itheta,iphi,sigmatheta,sigmaphi);
83     // Components different from ip=0 have the RMS bigger than mean
84     Float_t ptp[3]  =  { 1.219576,-0.354764,-0.690117 };
85     Float_t ptph[3] =  { 0.977522, 0.016269, 0.023158 }; 
86     Float_t pphp[3] =  { 1.303256,-0.464847,-0.869322 };
87
88     // Smeared momentum
89     pS = -1.;
90     //    Float_t dpmax = 5. + ip * 2.5; 
91     //    Float_t dpmax = 5. + ip * 2; 
92     Float_t dpmax;
93     if (sigmag2<999.) dpmax = 5. * (sigmap + sigmag2); 
94     else dpmax = 5. * sigmap;
95     Float_t dp = 100;
96     while (pS<0 || TMath::Abs(dp)>dpmax) { 
97       pS = p + fitp->GetRandom();
98       dp = pS - p; 
99     }
100     // Smeared phi
101     Float_t sigmaphiold=sigmaphi;
102     if (ip==0) sigmaphi *= pphp[0] + pphp[1] * dp + pphp[2] * dp*dp; 
103     if (sigmaphi<0.5 * sigmaphiold) sigmaphi = 0.5 * sigmaphiold;
104     if (sigmaphi>2.  * sigmaphiold) sigmaphi = 2.  * sigmaphiold;
105     phiS = phi + gRandom->Gaus(meanphi, sigmaphi);
106     Float_t dphi = phiS - phi; 
107     // Smeared theta
108     Float_t sigmathetaold=sigmatheta;
109     if (ip==0) sigmatheta *= ptp[0] + ptp[1] * dp + ptp[2] * dp*dp; 
110     if (ip==0) sigmatheta *= ptph[0] + ptph[1] * dphi + ptph[2] * dphi*dphi; 
111     if (sigmatheta<0.5 * sigmathetaold) sigmatheta = 0.5 * sigmathetaold;
112     if (sigmatheta>2.  * sigmathetaold) sigmatheta = 2.  * sigmathetaold;
113     thetaS = theta + gRandom->Gaus(meantheta,sigmatheta);
114 }
115
116
117
118
119
120