Another histos for lumi
[u/mrichter/AliRoot.git] / ITS / AliITSClusterParam.cxx
index d96a77d..1babafb 100644 (file)
 //#include <TProfile2D.h>
 #include <TVector3.h>
 #include "TMath.h"
+#include "AliTracker.h"
 #include "AliGeomManager.h"
 #include "AliITSRecPoint.h"
 #include "AliITSClusterParam.h"
 #include "AliITSReconstructor.h"
 #include "AliExternalTrackParam.h"
+#include "AliCheb3DCalc.h"
 
 ClassImp(AliITSClusterParam)
 
@@ -112,36 +114,64 @@ void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /*cl*/,
 Int_t AliITSClusterParam::GetError(Int_t layer,
                                   const AliITSRecPoint *cl,
                                   Float_t tgl,Float_t tgphitr,Float_t expQ,
-                                  Float_t &erry,Float_t &errz)
+                                  Float_t &erry,Float_t &errz,Float_t &covyz,
+                                  Bool_t addMisalErr)
 {
   //
   // Calculate cluster position error
   //
+  static Double_t bz = (Double_t)AliTracker::GetBz();
+  Int_t retval=0;
+  covyz=0.;
   switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
   case 0: 
-    return GetErrorOrigRecPoint(cl,erry,errz);
+    retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
     break;
   case 1: 
-    return GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
+    retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
     break;
   case 2: 
-    return GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
+    retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
+    break;
+  case 3: 
+    retval = GetErrorParamAngleOld(layer,cl,tgl,tgphitr,erry,errz);
     break;
   default: 
-    return GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
+    retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
     break;
   }
 
+  // for SSD use the error provided by the cluster finder 
+  // if single-sided clusters are enabled
+  if(layer>=4 && AliITSReconstructor::GetRecoParam()->GetUseBadChannelsInClusterFinderSSD()) { 
+    //printf("error 1 erry errz covyz %10.7f %10.7f %15.13f\n",erry,errz,covyz);
+    retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
+    //printf("type %d erry errz covyz %10.7f %10.7f %15.13f\n",cl->GetType(),erry,errz,covyz);
+  }
+  
+  if(addMisalErr) {
+    // add error due to misalignment (to be improved)
+    Float_t errmisalY2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer,bz)
+      *AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer,bz);
+    Float_t errmisalZ2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer,bz)
+      *AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer,bz);
+    erry = TMath::Sqrt(erry*erry+errmisalY2);
+    errz = TMath::Sqrt(errz*errz+errmisalZ2);
+  }
+
+  return retval;
+
 }
 //--------------------------------------------------------------------------
 Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
-                                  Float_t &erry,Float_t &errz)
+                                  Float_t &erry,Float_t &errz,Float_t &covyz)
 {
   //
   // Calculate cluster position error (just take error from AliITSRecPoint)
   //
   erry   = TMath::Sqrt(cl->GetSigmaY2()); 
   errz   = TMath::Sqrt(cl->GetSigmaZ2()); 
+  covyz  = cl->GetSigmaYZ();
 
   return 1;
 }
@@ -315,6 +345,99 @@ Int_t AliITSClusterParam::GetErrorParamAngle(Int_t layer,
   // residuals, as a function of angle between track and det module plane.
   // Origin: M.Lunardon, S.Moretto)
   //
+  const int   kNcfSPDResX = 21;
+  const float kCfSPDResX[kNcfSPDResX] = {+1.1201e+01,+2.0903e+00,-2.2909e-01,-2.6413e-01,+4.2135e-01,-3.7190e-01,
+                                        +4.2339e-01,+1.8679e-01,-5.1249e-01,+1.8421e-01,+4.8849e-02,-4.3127e-01,
+                                        -1.1148e-01,+3.1984e-03,-2.5743e-01,-6.6408e-02,+3.0756e-01,+2.6809e-01,
+                                        -5.0339e-03,-1.4964e-01,-1.1001e-01};
+  const float kSPDazMax=56.000000;
+  //
+  const int   kNcfSPDMeanX = 16;
+  const float kCfSPDMeanX[kNcfSPDMeanX] = {-1.2532e+00,-3.8185e-01,-8.9039e-01,+2.6648e+00,+7.0361e-01,+1.2298e+00,
+                                          +3.2871e-01,+7.8487e-02,-1.6792e-01,-1.3966e-01,-3.1670e-01,-2.1795e-01,
+                                          -1.9451e-01,-4.9347e-02,-1.9186e-01,-1.9195e-01};
+  //
+  const int   kNcfSPDResZ = 5;
+  const float kCfSPDResZ[kNcfSPDResZ] = {+9.2384e+01,+3.4352e-01,-2.7317e+01,-1.4642e-01,+2.0868e+00};
+  const float kSPDpolMin=34.358002, kSPDpolMax=145.000000;
+  //
+  const Double_t kMaxSigmaSDDx=100.;
+  const Double_t kMaxSigmaSDDz=400.;
+  const Double_t kMaxSigmaSSDx=100.;
+  const Double_t kMaxSigmaSSDz=1000.;
+  //  
+  const Double_t kParamSDDx[2]={30.93,0.059};
+  const Double_t kParamSDDz[2]={33.09,0.011};
+  const Double_t kParamSSDx[2]={18.64,-0.0046};
+  const Double_t kParamSSDz[2]={784.4,-0.828};
+  Double_t sigmax=1000.0,sigmaz=1000.0;
+  Double_t biasx = 0.0;
+
+  Int_t volId = (Int_t)cl->GetVolumeId();
+  Double_t rotMA[9]; AliGeomManager::GetRotation(volId,rotMA);      // misaligned rotation
+  Double_t rotOR[9]; AliGeomManager::GetOrigRotation(volId,rotOR);  // original rotation
+  // difference in phi of original and misaligned sensors
+  double cross = rotOR[1]*rotMA[4]-rotOR[4]*rotMA[1];
+  cross /= TMath::Sqrt( (1.-rotOR[7]*rotOR[7]) * (1.-rotMA[7]*rotMA[7]) );
+  Double_t angleAzi = TMath::Abs(TMath::ATan(tgphitr) - TMath::ASin(cross) );
+  Double_t anglePol = TMath::Abs(TMath::ATan(tgl));
+
+  if(angleAzi>0.5*TMath::Pi()) angleAzi = TMath::Pi()-angleAzi;
+  if(anglePol>0.5*TMath::Pi()) anglePol = TMath::Pi()-anglePol;
+  Double_t angleAziDeg = angleAzi*180./TMath::Pi();
+  Double_t anglePolDeg = anglePol*180./TMath::Pi();
+  
+  if(layer==0 || layer==1) { // SPD
+    //
+    float phiInt    = angleAziDeg/kSPDazMax; // mapped to -1:1
+    if (phiInt>1) phiInt = 1; else if (phiInt<-1) phiInt = -1;
+    float phiAbsInt = (TMath::Abs(angleAziDeg+angleAziDeg) - kSPDazMax)/kSPDazMax; // mapped to -1:1
+    if (phiAbsInt>1) phiAbsInt = 1; else if (phiAbsInt<-1) phiAbsInt = -1;
+    anglePolDeg += 90; // the parameterization was provided in polar angle (90 deg - normal to sensor)
+    float polInt   = (anglePolDeg+anglePolDeg - (kSPDpolMax+kSPDpolMin))/(kSPDpolMax-kSPDpolMin); // mapped to -1:1
+    if (polInt>1) polInt = 1; else if (polInt<-1) polInt = -1;
+    //
+    sigmax = AliCheb3DCalc::ChebEval1D(phiAbsInt, kCfSPDResX , kNcfSPDResX);
+    biasx  = AliCheb3DCalc::ChebEval1D(phiInt   , kCfSPDMeanX, kNcfSPDMeanX);
+    sigmaz = AliCheb3DCalc::ChebEval1D(polInt   , kCfSPDResZ , kNcfSPDResZ);
+    //
+    // for the moment for the SPD only, need to decide where to put it
+    biasx *= 1e-4;
+    
+  } else if(layer==2 || layer==3) { // SDD
+
+    sigmax = angleAziDeg*kParamSDDx[1]+kParamSDDx[0];
+    sigmaz = kParamSDDz[0]+kParamSDDz[1]*anglePolDeg;
+    if(sigmax > kMaxSigmaSDDx) sigmax = kMaxSigmaSDDx;
+    if(sigmaz > kMaxSigmaSDDz) sigmax = kMaxSigmaSDDz;
+    
+  } else if(layer==4 || layer==5) { // SSD
+
+    sigmax = angleAziDeg*kParamSSDx[1]+kParamSSDx[0];
+    sigmaz = kParamSSDz[0]+kParamSSDz[1]*anglePolDeg;
+    if(sigmax > kMaxSigmaSSDx) sigmax = kMaxSigmaSSDx;
+    if(sigmaz > kMaxSigmaSSDz) sigmax = kMaxSigmaSSDz;
+    
+  }
+
+  // convert from micron to cm
+  erry = 1.e-4*sigmax; 
+  errz = 1.e-4*sigmaz;
+  
+
+  return 1;
+}
+//--------------------------------------------------------------------------
+Int_t AliITSClusterParam::GetErrorParamAngleOld(Int_t layer,
+                                               const AliITSRecPoint *cl,
+                                               Float_t tgl,Float_t tgphitr,
+                                               Float_t &erry,Float_t &errz)
+{
+  //
+  // Calculate cluster position error (parametrization extracted from rp-hit
+  // residuals, as a function of angle between track and det module plane.
+  // Origin: M.Lunardon, S.Moretto)
+  //
 
   Double_t maxSigmaSPDx=100.;
   Double_t maxSigmaSPDz=400.;
@@ -332,41 +455,37 @@ Int_t AliITSClusterParam::GetErrorParamAngle(Int_t layer,
   Double_t sigmax=1000.0,sigmaz=1000.0;
   
   Int_t volId = (Int_t)cl->GetVolumeId();
-  Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
-  Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
-
-
-  Double_t phitr = TMath::ATan(tgphitr);
-  Double_t alpha = TMath::ATan2(tra[1],tra[0]);
-  Double_t phiglob = alpha+phitr;
-  Double_t p[3]; 
-  p[0] = TMath::Cos(phiglob);
-  p[1] = TMath::Sin(phiglob);
-  p[2] = tgl;
-  TVector3 pvec(p[0],p[1],p[2]);
-  TVector3 normvec(rot[1],-rot[0],0.);
-  Double_t angle = pvec.Angle(normvec);
-  if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
-  Double_t angleDeg = angle*180./TMath::Pi();
+  Double_t rotMA[9]; AliGeomManager::GetRotation(volId,rotMA);      // misaligned rotation
+  Double_t rotOR[9]; AliGeomManager::GetOrigRotation(volId,rotOR);  // original rotation
+  // difference in phi of original and misaligned sensors
+  double cross = rotOR[1]*rotMA[4]-rotOR[4]*rotMA[1];
+  cross /= TMath::Sqrt( (1.-rotOR[7]*rotOR[7]) * (1.-rotMA[7]*rotMA[7]) );
+  Double_t angleAzi = TMath::Abs(TMath::ATan(tgphitr) - TMath::ASin(cross) );
+  Double_t anglePol = TMath::Abs(TMath::ATan(tgl));
+
+  if(angleAzi>0.5*TMath::Pi()) angleAzi = TMath::Pi()-angleAzi;
+  if(anglePol>0.5*TMath::Pi()) anglePol = TMath::Pi()-anglePol;
+  Double_t angleAziDeg = angleAzi*180./TMath::Pi();
+  Double_t anglePolDeg = anglePol*180./TMath::Pi();
   
   if(layer==0 || layer==1) { // SPD
 
-    sigmax = TMath::Exp(angleDeg*paramSPDx[1]+paramSPDx[0])+paramSPDx[2];
-    sigmaz = paramSPDz[0]+paramSPDz[1]*angleDeg;
+    sigmax = TMath::Exp(angleAziDeg*paramSPDx[1]+paramSPDx[0])+paramSPDx[2];
+    sigmaz = paramSPDz[0]+paramSPDz[1]*anglePolDeg;
     if(sigmax > maxSigmaSPDx) sigmax = maxSigmaSPDx;
     if(sigmaz > maxSigmaSPDz) sigmax = maxSigmaSPDz;
 
   } else if(layer==2 || layer==3) { // SDD
 
-    sigmax = angleDeg*paramSDDx[1]+paramSDDx[0];
-    sigmaz = paramSDDz[0]+paramSDDz[1]*angleDeg;
+    sigmax = angleAziDeg*paramSDDx[1]+paramSDDx[0];
+    sigmaz = paramSDDz[0]+paramSDDz[1]*anglePolDeg;
     if(sigmax > maxSigmaSDDx) sigmax = maxSigmaSDDx;
     if(sigmaz > maxSigmaSDDz) sigmax = maxSigmaSDDz;
     
   } else if(layer==4 || layer==5) { // SSD
 
-    sigmax = angleDeg*paramSSDx[1]+paramSSDx[0];
-    sigmaz = paramSSDz[0]+paramSSDz[1]*angleDeg;
+    sigmax = angleAziDeg*paramSSDx[1]+paramSSDx[0];
+    sigmaz = paramSSDz[0]+paramSSDz[1]*anglePolDeg;
     if(sigmax > maxSigmaSSDx) sigmax = maxSigmaSSDx;
     if(sigmaz > maxSigmaSSDz) sigmax = maxSigmaSSDz;