renaming function to avoid library conflict (discovered in test/generators/TUHKMgen)
[u/mrichter/AliRoot.git] / FASTSIM / AliFastMuonTrackingRes.cxx
index e192c4e..c3bd513 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-*/
+/* $Id: */
+
+// Implementation of AliFastResponse for the Muon Spectrometer resolution.
+// The response depends on the charge of the muon and
+// the background level.
+// The class uses the instance of an object of type AliMUONFastTracking to 
+// obtain the smearing parameters.
+// Author: andreas.morsch@cern.ch
 
 #include "AliFastMuonTrackingRes.h"
 #include "AliMUONFastTracking.h"
@@ -26,13 +31,27 @@ ClassImp(AliFastMuonTrackingRes)
 
 
 AliFastMuonTrackingRes::AliFastMuonTrackingRes() :
-    AliFastResponse("Resolution", "Muon Tracking Resolution")
+    AliFastResponse("Resolution", "Muon Tracking Resolution"),
+    fBackground(0.),
+    fCharge(1.),
+    fFastTracking(0)
+{
+// Default constructor
+}
+
+AliFastMuonTrackingRes::AliFastMuonTrackingRes(const AliFastMuonTrackingRes & res)
+    :AliFastResponse(res),
+    fBackground(0.),
+    fCharge(1.),
+    fFastTracking(0)
 {
-    SetBackground();
+// Copy constructor
+    res.Copy(*this);
 }
 
 void AliFastMuonTrackingRes::Init()
 {
+// Initialisation
     fFastTracking = AliMUONFastTracking::Instance();
     fFastTracking->Init(fBackground);
 }
@@ -42,6 +61,9 @@ void AliFastMuonTrackingRes::Init()
 void AliFastMuonTrackingRes::Evaluate(Float_t   p,  Float_t  theta , Float_t   phi,
                                         Float_t& pS,  Float_t& thetaS, Float_t&  phiS)
 {
+//
+// Evaluate Gaussian smearing from given kinematics 
+//
 
     Double_t meanp    = fFastTracking->MeanP  (p, theta, phi, Int_t(fCharge));
     Double_t sigmap   = fFastTracking->SigmaP (p, theta, phi, Int_t(fCharge));
@@ -51,42 +73,68 @@ void AliFastMuonTrackingRes::Evaluate(Float_t   p,  Float_t  theta , Float_t   p
     Double_t sigmag2  = fFastTracking->SigmaG2(p, theta, phi, Int_t(fCharge));
     
     Int_t ip,itheta,iphi;
-    TF1* fitp = fFastTracking->GetFitP();
-    
     fFastTracking->GetIpIthetaIphi(p, theta, phi, Int_t(fCharge), ip, itheta, iphi);
-    if (sigmap == 0) {
-       printf ("WARNING!!! sigmap=0:    ");
-       printf ("ip= %d itheta = %d iphi = %d   ", ip, itheta, iphi);  
-       printf ("p= %f theta = %f phi = %f\n", p, theta, phi);  
+    TF1* fitp = fFastTracking->GetFitP(ip,itheta,iphi);
+    
+    Float_t curmeanp   = fitp->GetParameter(0); 
+    Float_t cursigmap  = fitp->GetParameter(1); 
+    Float_t cursigma1p = fitp->GetParameter(2); 
+    Float_t curnormg2  = fitp->GetParameter(3); 
+    Float_t curmeang2  = fitp->GetParameter(4); 
+    Float_t cursigmag2 = fitp->GetParameter(5); 
+    if (curmeanp != meanp || cursigmap != sigmap || cursigma1p != sigma1p || 
+       curnormg2 != normg2 || curmeang2 != meang2 || cursigmag2 != sigmag2){ 
+      printf ("Setting new parameters for ip=%d itheta=%d iphi=%d\n",ip,itheta,iphi); 
+      fitp->SetParameters(meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
     }
-
-    fitp->SetParameters(meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
   
     Double_t meantheta  = fFastTracking->MeanTheta (p, theta, phi, Int_t(fCharge));
     Double_t sigmatheta = fFastTracking->SigmaTheta(p, theta, phi, Int_t(fCharge));
     Double_t meanphi    = fFastTracking->MeanPhi   (p, theta, phi, Int_t(fCharge));
     Double_t sigmaphi   = fFastTracking->SigmaPhi  (p, theta, phi, Int_t(fCharge));
-  
+    
+    if (sigmatheta<0 || sigmaphi<0) 
+      printf ("bin %d %d %d sigmatheta = %f, sigmaphi = %f\n",
+             ip,itheta,iphi,sigmatheta,sigmaphi);
     // Components different from ip=0 have the RMS bigger than mean
     Float_t ptp[3]  =  { 1.219576,-0.354764,-0.690117 };
     Float_t ptph[3] =  { 0.977522, 0.016269, 0.023158 }; 
     Float_t pphp[3] =  { 1.303256,-0.464847,-0.869322 };
 
     // Smeared momentum
-    pS   = p + fitp->GetRandom();
-    Float_t dp = pS - p; 
-
+    pS = -1.;
+    //    Float_t dpmax = 5. + ip * 2.5; 
+    //    Float_t dpmax = 5. + ip * 2; 
+    Float_t dpmax;
+    if (sigmag2<999.) dpmax = 5. * TMath::Abs(sigmap + sigmag2); 
+    else dpmax = 5. * TMath::Abs(sigmap);
+    Float_t dp = 100;
+    while (pS<0 || TMath::Abs(dp)>dpmax) { 
+      pS = p + fitp->GetRandom();
+      dp = pS - p; 
+    }
     // Smeared phi
+    Float_t sigmaphiold=sigmaphi;
     if (ip==0) sigmaphi *= pphp[0] + pphp[1] * dp + pphp[2] * dp*dp; 
+    if (sigmaphi<0.5 * sigmaphiold) sigmaphi = 0.5 * sigmaphiold;
+    if (sigmaphi>2.  * sigmaphiold) sigmaphi = 2.  * sigmaphiold;
     phiS = phi + gRandom->Gaus(meanphi, sigmaphi);
     Float_t dphi = phiS - phi; 
     // Smeared theta
+    Float_t sigmathetaold=sigmatheta;
     if (ip==0) sigmatheta *= ptp[0] + ptp[1] * dp + ptp[2] * dp*dp; 
     if (ip==0) sigmatheta *= ptph[0] + ptph[1] * dphi + ptph[2] * dphi*dphi; 
+    if (sigmatheta<0.5 * sigmathetaold) sigmatheta = 0.5 * sigmathetaold;
+    if (sigmatheta>2.  * sigmathetaold) sigmatheta = 2.  * sigmathetaold;
     thetaS = theta + gRandom->Gaus(meantheta,sigmatheta);
 }
 
-
+AliFastMuonTrackingRes& AliFastMuonTrackingRes::operator=(const  AliFastMuonTrackingRes& rhs)
+{
+// Assignment operator
+    rhs.Copy(*this);
+    return *this;
+}