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