[u/mrichter/AliRoot.git] / ITS / AliITSClusterParam.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  ITS cluster error and shape parameterization                             //
//                                                                           //
//  andrea.dainese@lnl.infn.it                                               //
///////////////////////////////////////////////////////////////////////////////
//#include "TFile.h"
//#include "TTree.h"
//#include <TVectorF.h>
//#include <TLinearFitter.h>
//#include <TH1F.h>
//#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"

ClassImp(AliITSClusterParam)


AliITSClusterParam* AliITSClusterParam::fgInstance = 0;


/*
  Example usage fitting parameterization (NOT YET...):
  TFile fres("resol.root");    //tree with resolution and shape 
  TTree * treeRes =(TTree*)fres.Get("Resol");
  
  AliITSClusterParam param;
  param.SetInstance(&param);
  param.FitResol(treeRes);
  param.FitRMS(treeRes);
  TFile fparam("ITSClusterParam.root","recreate");
  param.Write("Param");
  //
  //
  TFile fparam("ITSClusterParam.root");
  AliITSClusterParam *param2  =  (AliITSClusterParam *) fparam.Get("Param"); 
  param2->SetInstance(param2);
  param2->Test(treeRes);
  

  treeRes->Draw("(Resol-AliITSClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
*/

//-------------------------------------------------------------------------
AliITSClusterParam* AliITSClusterParam::Instance()
{
  //
  // Singleton implementation
  // Returns an instance of this class, it is created if neccessary
  //
  if (fgInstance == 0){
    fgInstance = new AliITSClusterParam();
  }
  return fgInstance;
}
//-------------------------------------------------------------------------
void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /*cl*/, 
				  Float_t tgl,Float_t tgphitr,
				  Float_t &ny,Float_t &nz)
{
  //
  // Get "mean shape" (original parametrization from AliITStrackerMI)
  //
  tgl = TMath::Abs(tgl);
  tgphitr = TMath::Abs(tgphitr);

  // SPD
  if (layer==0) {
    ny = 1.+tgphitr*3.2;
    nz = 1.+tgl*0.34;
    return;
  }
  if (layer==1) {
    ny = 1.+tgphitr*3.2;
    nz = 1.+tgl*0.28;
    return;
  }
  // SSD
  if (layer==4 || layer==5) {
    ny = 2.02+tgphitr*1.95;
    nz = 2.02+tgphitr*2.35;
    return;
  }
  // SDD
  ny  = 6.6-2.7*tgphitr;
  nz  = 2.8-3.11*tgphitr+0.45*tgl;
  return;
}
//--------------------------------------------------------------------------
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 &covyz,
				   Bool_t addMisalErr)
{
  //
  // Calculate cluster position error
  //
  Int_t retval=0;
  covyz=0.;
  switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
  case 0: 
    retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
    break;
  case 1: 
    retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
    break;
  case 2: 
    retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
    break;
  default: 
    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) {
    Double_t bz = (Double_t)AliTracker::GetBz();
    // 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 &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;
}
//--------------------------------------------------------------------------
Int_t AliITSClusterParam::GetErrorParamMI(Int_t layer,const AliITSRecPoint*cl,
					  Float_t tgl,Float_t tgphitr,
					  Float_t expQ,
					  Float_t &erry,Float_t &errz)
{
  //
  // Calculate cluster position error (original parametrization from 
  // AliITStrackerMI)
  //
  Float_t nz,ny;
  GetNTeor(layer, cl,tgl,tgphitr,ny,nz);  
  erry   = TMath::Sqrt(cl->GetSigmaY2()); 
  errz   = TMath::Sqrt(cl->GetSigmaZ2()); 
  //
  // PIXELS
  if (layer<2){
    if (TMath::Abs(ny-cl->GetNy())>0.6)  {
      if (ny<cl->GetNy()){
	erry*=0.4+TMath::Abs(ny-cl->GetNy());
	errz*=0.4+TMath::Abs(ny-cl->GetNy());
      }else{
	erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
	errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
      }
    }
    if (TMath::Abs(nz-cl->GetNz())>1.)  {
      erry*=TMath::Abs(nz-cl->GetNz());
      errz*=TMath::Abs(nz-cl->GetNz());	      
    }
    erry*=0.85;
    errz*=0.85;
    erry= TMath::Min(erry,float(0.0050));
    errz= TMath::Min(errz,float(0.0300));
    return 10;
  }

  //STRIPS
  if (layer>3){ 
    //factor 1.8 appears in new simulation
    //
    Float_t scale=1.8;
    if (cl->GetNy()==100||cl->GetNz()==100){
      erry = 0.004*scale;
      errz = 0.2*scale;
      return 100;
    }
    if (cl->GetNy()+cl->GetNz()>12){
      erry = 0.06*scale;
      errz = 0.57*scale;
      return 100;
    }
    Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
    Float_t chargematch = TMath::Max(double(normq/expQ),2.);
    //
    if (cl->GetType()==1 || cl->GetType()==10 ){     							       
      if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
	errz = 0.043*scale;
	erry = 0.00094*scale;
	return 101;
      }
      if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
	errz = 0.06*scale;
	erry =0.0013*scale;
	return 102;
      }
      erry = 0.0027*scale;
      errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
      return 103;
    }
    if (cl->GetType()==2 || cl->GetType()==11 ){ 
      erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
      errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
      return 104;
    }
    
    if (cl->GetType()>100 ){     							       
      if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
	errz = 0.05*scale;
	erry = 0.00096*scale;
	return 105;
      }
      if (cl->GetNy()+cl->GetNz()-nz-ny<1){
	errz = 0.10*scale;
	erry = 0.0025*scale;
	return 106;
      }

      errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
      erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
      return 107;
    }    
    Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
    if (diff<1) diff=1;
    if (diff>4) diff=4;
        
    if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
      errz = 0.14*diff;
      erry = 0.003*diff;
      return 108;
    }  
    erry = 0.04*diff;
    errz = 0.06*diff;
    return 109;
  }
  //DRIFTS
  Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
  Float_t chargematch = normq/expQ;
  chargematch/=2.4; // F. Prino Sept. 2007: SDD charge conversion keV->ADC
  Float_t factorz=1;
  Int_t   cnz = cl->GetNz()%10;
  //charge match
  if (cl->GetType()==1){
    if (chargematch<1.25){
      erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
    }
    else{
      erry = 0.003*chargematch;
      if (cl->GetNz()==3) erry*=1.5;
    }
    if (chargematch<1.0){
      errz =  0.0011*(1.+6./cl->GetQ());
    }
    else{
      errz = 0.002*(1+2*(chargematch-1.));
    }
    if (cnz>nz+0.6) {
      erry*=(cnz-nz+0.5);
      errz*=1.4*(cnz-nz+0.5);
    }
  }
  if (cl->GetType()>1){
    if (chargematch<1){
      erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
      errz =  0.0016*(1.+6./cl->GetQ());
    }
    else{
      errz = 0.0014*(1+3*(chargematch-1.));
      erry = 0.003*(1+3*(chargematch-1.));
    } 
    if (cnz>nz+0.6) {
      erry*=(cnz-nz+0.5);
      errz*=1.4*(cnz-nz+0.5);
    }
  }

  if (TMath::Abs(cl->GetY())>2.5){
    factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
  }
  if (TMath::Abs(cl->GetY())<1){
    factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
  }
  factorz= TMath::Min(factorz,float(4.));  
  errz*=factorz;

  erry= TMath::Min(erry,float(0.05));
  errz= TMath::Min(errz,float(0.05));  
  return 200;
}
//--------------------------------------------------------------------------
Int_t AliITSClusterParam::GetErrorParamAngle(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.;
  Double_t maxSigmaSDDx=100.;
  Double_t maxSigmaSDDz=400.;
  Double_t maxSigmaSSDx=100.;
  Double_t maxSigmaSSDz=1000.;
  
  Double_t paramSPDx[3]={-6.417,0.18,11.14};
  Double_t paramSPDz[2]={118.,-0.155};
  Double_t paramSDDx[2]={30.93,0.059};
  Double_t paramSDDz[2]={33.09,0.011};
  Double_t paramSSDx[2]={18.64,-0.0046};
  Double_t paramSSDz[2]={784.4,-0.828};
  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[4],rot[7]);
  Double_t angle = pvec.Angle(normvec);
  if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
  Double_t angleDeg = angle*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;
    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;
    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;
    if(sigmax > maxSigmaSSDx) sigmax = maxSigmaSSDx;
    if(sigmaz > maxSigmaSSDz) sigmax = maxSigmaSSDz;
    
  }

  // convert from micron to cm
  erry = 1.e-4*sigmax; 
  errz = 1.e-4*sigmaz;
  
  return 1;
}
//--------------------------------------------------------------------------
void AliITSClusterParam::Print(Option_t* /*option*/) const {
  //
  // Print param Information
  //

  //
  // Error parameterization
  //
  printf("NOT YET...\n");
  return;
}
Commit	Line	Data
572f41f9	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, 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	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// ITS cluster error and shape parameterization //
	21	// //
	22	// andrea.dainese@lnl.infn.it //
	23	///////////////////////////////////////////////////////////////////////////////
	24	//#include "TFile.h"
	25	//#include "TTree.h"
	26	//#include <TVectorF.h>
	27	//#include <TLinearFitter.h>
	28	//#include <TH1F.h>
	29	//#include <TProfile2D.h>
e50912db	30	#include <TVector3.h>
572f41f9	31	#include "TMath.h"
d9ead1a0	32	#include "AliTracker.h"
e50912db	33	#include "AliGeomManager.h"
572f41f9	34	#include "AliITSRecPoint.h"
572f41f9	35	#include "AliITSClusterParam.h"
e50912db	36	#include "AliITSReconstructor.h"
e50912db	37	#include "AliExternalTrackParam.h"
572f41f9	38
	39	ClassImp(AliITSClusterParam)
	40
	41
	42	AliITSClusterParam* AliITSClusterParam::fgInstance = 0;
	43
	44
	45	/*
	46	Example usage fitting parameterization (NOT YET...):
	47	TFile fres("resol.root"); //tree with resolution and shape
	48	TTree * treeRes =(TTree*)fres.Get("Resol");
	49
	50	AliITSClusterParam param;
	51	param.SetInstance(&param);
	52	param.FitResol(treeRes);
	53	param.FitRMS(treeRes);
	54	TFile fparam("ITSClusterParam.root","recreate");
	55	param.Write("Param");
	56	//
	57	//
	58	TFile fparam("ITSClusterParam.root");
	59	AliITSClusterParam param2 = (AliITSClusterParam ) fparam.Get("Param");
	60	param2->SetInstance(param2);
	61	param2->Test(treeRes);
	62
	63
	64	treeRes->Draw("(Resol-AliITSClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
	65	*/
	66
	67	//-------------------------------------------------------------------------
	68	AliITSClusterParam* AliITSClusterParam::Instance()
	69	{
	70	//
	71	// Singleton implementation
	72	// Returns an instance of this class, it is created if neccessary
	73	//
	74	if (fgInstance == 0){
	75	fgInstance = new AliITSClusterParam();
	76	}
	77	return fgInstance;
	78	}
	79	//-------------------------------------------------------------------------
	80	void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /cl/,
e50912db	81	Float_t tgl,Float_t tgphitr,
572f41f9	82	Float_t &ny,Float_t &nz)
	83	{
	84	//
e50912db	85	// Get "mean shape" (original parametrization from AliITStrackerMI)
572f41f9	86	//
e50912db	87	tgl = TMath::Abs(tgl);
	88	tgphitr = TMath::Abs(tgphitr);
	89
	90	// SPD
	91	if (layer==0) {
	92	ny = 1.+tgphitr*3.2;
	93	nz = 1.+tgl*0.34;
572f41f9	94	return;
572f41f9	95	}
e50912db	96	if (layer==1) {
	97	ny = 1.+tgphitr*3.2;
	98	nz = 1.+tgl*0.28;
572f41f9	99	return;
572f41f9	100	}
e50912db	101	// SSD
	102	if (layer==4 \|\| layer==5) {
	103	ny = 2.02+tgphitr*1.95;
	104	nz = 2.02+tgphitr*2.35;
572f41f9	105	return;
572f41f9	106	}
e50912db	107	// SDD
	108	ny = 6.6-2.7*tgphitr;
	109	nz = 2.8-3.11tgphitr+0.45tgl;
	110	return;
	111	}
	112	//--------------------------------------------------------------------------
	113	Int_t AliITSClusterParam::GetError(Int_t layer,
	114	const AliITSRecPoint *cl,
	115	Float_t tgl,Float_t tgphitr,Float_t expQ,
d9ead1a0	116	Float_t &erry,Float_t &errz,Float_t &covyz,
8c139cf3	117	Bool_t addMisalErr)
e50912db	118	{
	119	//
	120	// Calculate cluster position error
	121	//
401eff16	122	Int_t retval=0;
d9ead1a0	123	covyz=0.;
e50912db	124	switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
e50912db	125	case 0:
d9ead1a0	126	retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
e50912db	127	break;
e50912db	128	case 1:
401eff16	129	retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
e50912db	130	break;
e50912db	131	case 2:
401eff16	132	retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
e50912db	133	break;
e50912db	134	default:
401eff16	135	retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
e50912db	136	break;
	137	}
	138
d9ead1a0	139	// for SSD use the error provided by the cluster finder
	140	// if single-sided clusters are enabled
	141	if(layer>=4 && AliITSReconstructor::GetRecoParam()->GetUseBadChannelsInClusterFinderSSD()) {
	142	//printf("error 1 erry errz covyz %10.7f %10.7f %15.13f\n",erry,errz,covyz);
	143	retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
	144	//printf("type %d erry errz covyz %10.7f %10.7f %15.13f\n",cl->GetType(),erry,errz,covyz);
	145	}
	146
8c139cf3	147	if(addMisalErr) {
d9ead1a0	148	Double_t bz = (Double_t)AliTracker::GetBz();
8c139cf3	149	// add error due to misalignment (to be improved)
d9ead1a0	150	Float_t errmisalY2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer,bz)
	151	*AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer,bz);
	152	Float_t errmisalZ2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer,bz)
	153	*AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer,bz);
f9119eb9	154	erry = TMath::Sqrt(erry*erry+errmisalY2);
f9119eb9	155	errz = TMath::Sqrt(errz*errz+errmisalZ2);
8c139cf3	156	}
401eff16	157
	158	return retval;
	159
572f41f9	160	}
572f41f9	161	//--------------------------------------------------------------------------
e50912db	162	Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
d9ead1a0	163	Float_t &erry,Float_t &errz,Float_t &covyz)
572f41f9	164	{
e50912db	165	//
	166	// Calculate cluster position error (just take error from AliITSRecPoint)
	167	//
	168	erry = TMath::Sqrt(cl->GetSigmaY2());
	169	errz = TMath::Sqrt(cl->GetSigmaZ2());
d9ead1a0	170	covyz = cl->GetSigmaYZ();
e50912db	171
	172	return 1;
	173	}
	174	//--------------------------------------------------------------------------
	175	Int_t AliITSClusterParam::GetErrorParamMI(Int_t layer,const AliITSRecPoint*cl,
	176	Float_t tgl,Float_t tgphitr,
	177	Float_t expQ,
	178	Float_t &erry,Float_t &errz)
	179	{
	180	//
	181	// Calculate cluster position error (original parametrization from
	182	// AliITStrackerMI)
572f41f9	183	//
572f41f9	184	Float_t nz,ny;
e50912db	185	GetNTeor(layer, cl,tgl,tgphitr,ny,nz);
572f41f9	186	erry = TMath::Sqrt(cl->GetSigmaY2());
	187	errz = TMath::Sqrt(cl->GetSigmaZ2());
	188	//
	189	// PIXELS
	190	if (layer<2){
	191	if (TMath::Abs(ny-cl->GetNy())>0.6) {
	192	if (ny<cl->GetNy()){
	193	erry*=0.4+TMath::Abs(ny-cl->GetNy());
	194	errz*=0.4+TMath::Abs(ny-cl->GetNy());
	195	}else{
	196	erry=0.7+0.5TMath::Abs(ny-cl->GetNy());
	197	errz=0.7+0.5TMath::Abs(ny-cl->GetNy());
	198	}
	199	}
	200	if (TMath::Abs(nz-cl->GetNz())>1.) {
	201	erry*=TMath::Abs(nz-cl->GetNz());
	202	errz*=TMath::Abs(nz-cl->GetNz());
	203	}
	204	erry*=0.85;
	205	errz*=0.85;
	206	erry= TMath::Min(erry,float(0.0050));
	207	errz= TMath::Min(errz,float(0.0300));
	208	return 10;
	209	}
	210
	211	//STRIPS
	212	if (layer>3){
	213	//factor 1.8 appears in new simulation
	214	//
	215	Float_t scale=1.8;
	216	if (cl->GetNy()==100\|\|cl->GetNz()==100){
	217	erry = 0.004*scale;
	218	errz = 0.2*scale;
	219	return 100;
	220	}
	221	if (cl->GetNy()+cl->GetNz()>12){
	222	erry = 0.06*scale;
	223	errz = 0.57*scale;
	224	return 100;
	225	}
e50912db	226	Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgltgl+tgphitrtgphitr));
572f41f9	227	Float_t chargematch = TMath::Max(double(normq/expQ),2.);
	228	//
	229	if (cl->GetType()==1 \|\| cl->GetType()==10 ){
	230	if (chargematch<1.0 \|\| (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
	231	errz = 0.043*scale;
	232	erry = 0.00094*scale;
	233	return 101;
	234	}
	235	if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
	236	errz = 0.06*scale;
	237	erry =0.0013*scale;
	238	return 102;
	239	}
	240	erry = 0.0027*scale;
	241	errz = TMath::Min(0.028(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)scale;
	242	return 103;
	243	}
	244	if (cl->GetType()==2 \|\| cl->GetType()==11 ){
	245	erry = TMath::Min(0.0010(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)scale;
	246	errz = TMath::Min(0.025(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)scale;
	247	return 104;
	248	}
	249
	250	if (cl->GetType()>100 ){
	251	if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
	252	errz = 0.05*scale;
	253	erry = 0.00096*scale;
	254	return 105;
	255	}
	256	if (cl->GetNy()+cl->GetNz()-nz-ny<1){
	257	errz = 0.10*scale;
	258	erry = 0.0025*scale;
	259	return 106;
	260	}
	261
	262	errz = TMath::Min(0.05(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)scale;
	263	erry = TMath::Min(0.003(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)scale;
	264	return 107;
	265	}
	266	Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
	267	if (diff<1) diff=1;
	268	if (diff>4) diff=4;
	269
	270	if (cl->GetType()==5\|\|cl->GetType()==6\|\|cl->GetType()==7\|\|cl->GetType()==8){
	271	errz = 0.14*diff;
	272	erry = 0.003*diff;
	273	return 108;
	274	}
	275	erry = 0.04*diff;
	276	errz = 0.06*diff;
	277	return 109;
	278	}
	279	//DRIFTS
e50912db	280	Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgltgl+tgphitrtgphitr));
572f41f9	281	Float_t chargematch = normq/expQ;
	282	chargematch/=2.4; // F. Prino Sept. 2007: SDD charge conversion keV->ADC
	283	Float_t factorz=1;
	284	Int_t cnz = cl->GetNz()%10;
	285	//charge match
	286	if (cl->GetType()==1){
	287	if (chargematch<1.25){
	288	erry = 0.0028*(1.+6./cl->GetQ()); // gold clusters
	289	}
	290	else{
	291	erry = 0.003*chargematch;
	292	if (cl->GetNz()==3) erry*=1.5;
	293	}
	294	if (chargematch<1.0){
	295	errz = 0.0011*(1.+6./cl->GetQ());
	296	}
	297	else{
	298	errz = 0.002(1+2(chargematch-1.));
	299	}
	300	if (cnz>nz+0.6) {
	301	erry*=(cnz-nz+0.5);
	302	errz=1.4(cnz-nz+0.5);
	303	}
	304	}
	305	if (cl->GetType()>1){
	306	if (chargematch<1){
	307	erry = 0.00385*(1.+6./cl->GetQ()); // gold clusters
	308	errz = 0.0016*(1.+6./cl->GetQ());
	309	}
	310	else{
	311	errz = 0.0014(1+3(chargematch-1.));
	312	erry = 0.003(1+3(chargematch-1.));
	313	}
	314	if (cnz>nz+0.6) {
	315	erry*=(cnz-nz+0.5);
	316	errz=1.4(cnz-nz+0.5);
	317	}
	318	}
	319
	320	if (TMath::Abs(cl->GetY())>2.5){
	321	factorz=1+2(TMath::Abs(cl->GetY())-2.5);
	322	}
	323	if (TMath::Abs(cl->GetY())<1){
	324	factorz=1.+0.5TMath::Abs(TMath::Abs(cl->GetY())-1.);
	325	}
	326	factorz= TMath::Min(factorz,float(4.));
	327	errz*=factorz;
	328
	329	erry= TMath::Min(erry,float(0.05));
	330	errz= TMath::Min(errz,float(0.05));
	331	return 200;
	332	}
	333	//--------------------------------------------------------------------------
e50912db	334	Int_t AliITSClusterParam::GetErrorParamAngle(Int_t layer,
	335	const AliITSRecPoint *cl,
	336	Float_t tgl,Float_t tgphitr,
	337	Float_t &erry,Float_t &errz)
	338	{
	339	//
	340	// Calculate cluster position error (parametrization extracted from rp-hit
	341	// residuals, as a function of angle between track and det module plane.
	342	// Origin: M.Lunardon, S.Moretto)
	343	//
	344
	345	Double_t maxSigmaSPDx=100.;
	346	Double_t maxSigmaSPDz=400.;
	347	Double_t maxSigmaSDDx=100.;
	348	Double_t maxSigmaSDDz=400.;
	349	Double_t maxSigmaSSDx=100.;
	350	Double_t maxSigmaSSDz=1000.;
	351
	352	Double_t paramSPDx[3]={-6.417,0.18,11.14};
	353	Double_t paramSPDz[2]={118.,-0.155};
	354	Double_t paramSDDx[2]={30.93,0.059};
	355	Double_t paramSDDz[2]={33.09,0.011};
	356	Double_t paramSSDx[2]={18.64,-0.0046};
	357	Double_t paramSSDz[2]={784.4,-0.828};
	358	Double_t sigmax=1000.0,sigmaz=1000.0;
	359
	360	Int_t volId = (Int_t)cl->GetVolumeId();
	361	Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
	362	Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
	363
	364
	365	Double_t phitr = TMath::ATan(tgphitr);
	366	Double_t alpha = TMath::ATan2(tra[1],tra[0]);
	367	Double_t phiglob = alpha+phitr;
	368	Double_t p[3];
	369	p[0] = TMath::Cos(phiglob);
	370	p[1] = TMath::Sin(phiglob);
	371	p[2] = tgl;
	372	TVector3 pvec(p[0],p[1],p[2]);
c3a24b56	373	TVector3 normvec(rot[1],rot[4],rot[7]);
e50912db	374	Double_t angle = pvec.Angle(normvec);
e50912db	375	if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
0ffe9832	376	Double_t angleDeg = angle*180./TMath::Pi();
e50912db	377
	378	if(layer==0 \|\| layer==1) { // SPD
	379
	380	sigmax = TMath::Exp(angleDeg*paramSPDx[1]+paramSPDx[0])+paramSPDx[2];
	381	sigmaz = paramSPDz[0]+paramSPDz[1]*angleDeg;
	382	if(sigmax > maxSigmaSPDx) sigmax = maxSigmaSPDx;
	383	if(sigmaz > maxSigmaSPDz) sigmax = maxSigmaSPDz;
	384
	385	} else if(layer==2 \|\| layer==3) { // SDD
	386
	387	sigmax = angleDeg*paramSDDx[1]+paramSDDx[0];
	388	sigmaz = paramSDDz[0]+paramSDDz[1]*angleDeg;
	389	if(sigmax > maxSigmaSDDx) sigmax = maxSigmaSDDx;
	390	if(sigmaz > maxSigmaSDDz) sigmax = maxSigmaSDDz;
	391
	392	} else if(layer==4 \|\| layer==5) { // SSD
	393
	394	sigmax = angleDeg*paramSSDx[1]+paramSSDx[0];
	395	sigmaz = paramSSDz[0]+paramSSDz[1]*angleDeg;
	396	if(sigmax > maxSigmaSSDx) sigmax = maxSigmaSSDx;
	397	if(sigmaz > maxSigmaSSDz) sigmax = maxSigmaSSDz;
	398
	399	}
	400
	401	// convert from micron to cm
	402	erry = 1.e-4*sigmax;
	403	errz = 1.e-4*sigmaz;
	404
	405	return 1;
	406	}
	407	//--------------------------------------------------------------------------
572f41f9	408	void AliITSClusterParam::Print(Option_t* /option/) const {
	409	//
	410	// Print param Information
	411	//
	412
	413	//
	414	// Error parameterization
	415	//
	416	printf("NOT YET...\n");
	417	return;
	418	}
	419
	420
	421
	422
	423