ITS/AliITSClusterParam.cxx

   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>
  30 #include <TVector3.h>
  31 #include "TMath.h"
  32 #include "AliTracker.h"
  33 #include "AliGeomManager.h"
  34 #include "AliITSRecPoint.h"
  35 #include "AliITSClusterParam.h"
  36 #include "AliITSReconstructor.h"
  37 #include "AliExternalTrackParam.h"
  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*/,
  81                                   Float_t tgl,Float_t tgphitr,
  82                                   Float_t &ny,Float_t &nz)
  83 {
  84   //
  85   // Get "mean shape" (original parametrization from AliITStrackerMI)
  86   //
  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;
  94     return;
  95   }
  96   if (layer==1) {
  97     ny = 1.+tgphitr*3.2;
  98     nz = 1.+tgl*0.28;
  99     return;
 100   }
 101   // SSD
 102   if (layer==4 || layer==5) {
 103     ny = 2.02+tgphitr*1.95;
 104     nz = 2.02+tgphitr*2.35;
 105     return;
 106   }
 107   // SDD
 108   ny  = 6.6-2.7*tgphitr;
 109   nz  = 2.8-3.11*tgphitr+0.45*tgl;
 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,
 116                                    Float_t &erry,Float_t &errz,Float_t &covyz,
 117                                    Bool_t addMisalErr)
 118 {
 119   //
 120   // Calculate cluster position error
 121   //
 122   Int_t retval=0;
 123   covyz=0.;
 124   switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
 125   case 0:
 126     retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
 127     break;
 128   case 1:
 129     retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
 130     break;
 131   case 2:
 132     retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
 133     break;
 134   default:
 135     retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
 136     break;
 137   }
 138
 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
 147   if(addMisalErr) {
 148     Double_t bz = (Double_t)AliTracker::GetBz();
 149     // add error due to misalignment (to be improved)
 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);
 154     erry = TMath::Sqrt(erry*erry+errmisalY2);
 155     errz = TMath::Sqrt(errz*errz+errmisalZ2);
 156   }
 157
 158   return retval;
 159
 160 }
 161 //--------------------------------------------------------------------------
 162 Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
 163                                    Float_t &erry,Float_t &errz,Float_t &covyz)
 164 {
 165   //
 166   // Calculate cluster position error (just take error from AliITSRecPoint)
 167   //
 168   erry   = TMath::Sqrt(cl->GetSigmaY2());
 169   errz   = TMath::Sqrt(cl->GetSigmaZ2());
 170   covyz  = cl->GetSigmaYZ();
 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)
 183   //
 184   Float_t nz,ny;
 185   GetNTeor(layer, cl,tgl,tgphitr,ny,nz);
 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.5*TMath::Abs(ny-cl->GetNy());
 197         errz*=0.7+0.5*TMath::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     }
 226     Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
 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
 280   Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
 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.5*TMath::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 //--------------------------------------------------------------------------
 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]);
 373   TVector3 normvec(rot[1],rot[4],rot[7]);
 374   Double_t angle = pvec.Angle(normvec);
 375   if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
 376   Double_t angleDeg = angle*180./TMath::Pi();
 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 //--------------------------------------------------------------------------
 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