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                                    Bool_t addMisalErr)
 117 {
 118   //
 119   // Calculate cluster position error
 120   //
 121   Int_t retval=0;
 122   switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
 123   case 0:
 124     retval = GetErrorOrigRecPoint(cl,erry,errz);
 125     break;
 126   case 1:
 127     retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
 128     break;
 129   case 2:
 130     retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
 131     break;
 132   default:
 133     retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
 134     break;
 135   }
 136
 137   if(addMisalErr) {
 138     // add error due to misalignment (to be improved)
 139     Float_t errmisalY2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer)
 140       *AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer);
 141     Float_t errmisalZ2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer)
 142       *AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer);
 143     erry = TMath::Sqrt(erry*erry+errmisalY2);
 144     errz = TMath::Sqrt(errz*errz+errmisalZ2);
 145   }
 146
 147   return retval;
 148
 149 }
 150 //--------------------------------------------------------------------------
 151 Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
 152                                    Float_t &erry,Float_t &errz)
 153 {
 154   //
 155   // Calculate cluster position error (just take error from AliITSRecPoint)
 156   //
 157   erry   = TMath::Sqrt(cl->GetSigmaY2());
 158   errz   = TMath::Sqrt(cl->GetSigmaZ2());
 159
 160   return 1;
 161 }
 162 //--------------------------------------------------------------------------
 163 Int_t AliITSClusterParam::GetErrorParamMI(Int_t layer,const AliITSRecPoint*cl,
 164                                           Float_t tgl,Float_t tgphitr,
 165                                           Float_t expQ,
 166                                           Float_t &erry,Float_t &errz)
 167 {
 168   //
 169   // Calculate cluster position error (original parametrization from
 170   // AliITStrackerMI)
 171   //
 172   Float_t nz,ny;
 173   GetNTeor(layer, cl,tgl,tgphitr,ny,nz);
 174   erry   = TMath::Sqrt(cl->GetSigmaY2());
 175   errz   = TMath::Sqrt(cl->GetSigmaZ2());
 176   //
 177   // PIXELS
 178   if (layer<2){
 179     if (TMath::Abs(ny-cl->GetNy())>0.6)  {
 180       if (ny<cl->GetNy()){
 181         erry*=0.4+TMath::Abs(ny-cl->GetNy());
 182         errz*=0.4+TMath::Abs(ny-cl->GetNy());
 183       }else{
 184         erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
 185         errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
 186       }
 187     }
 188     if (TMath::Abs(nz-cl->GetNz())>1.)  {
 189       erry*=TMath::Abs(nz-cl->GetNz());
 190       errz*=TMath::Abs(nz-cl->GetNz());
 191     }
 192     erry*=0.85;
 193     errz*=0.85;
 194     erry= TMath::Min(erry,float(0.0050));
 195     errz= TMath::Min(errz,float(0.0300));
 196     return 10;
 197   }
 198
 199   //STRIPS
 200   if (layer>3){
 201     //factor 1.8 appears in new simulation
 202     //
 203     Float_t scale=1.8;
 204     if (cl->GetNy()==100||cl->GetNz()==100){
 205       erry = 0.004*scale;
 206       errz = 0.2*scale;
 207       return 100;
 208     }
 209     if (cl->GetNy()+cl->GetNz()>12){
 210       erry = 0.06*scale;
 211       errz = 0.57*scale;
 212       return 100;
 213     }
 214     Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
 215     Float_t chargematch = TMath::Max(double(normq/expQ),2.);
 216     //
 217     if (cl->GetType()==1 || cl->GetType()==10 ){
 218       if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
 219         errz = 0.043*scale;
 220         erry = 0.00094*scale;
 221         return 101;
 222       }
 223       if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
 224         errz = 0.06*scale;
 225         erry =0.0013*scale;
 226         return 102;
 227       }
 228       erry = 0.0027*scale;
 229       errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
 230       return 103;
 231     }
 232     if (cl->GetType()==2 || cl->GetType()==11 ){
 233       erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
 234       errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
 235       return 104;
 236     }
 237
 238     if (cl->GetType()>100 ){
 239       if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
 240         errz = 0.05*scale;
 241         erry = 0.00096*scale;
 242         return 105;
 243       }
 244       if (cl->GetNy()+cl->GetNz()-nz-ny<1){
 245         errz = 0.10*scale;
 246         erry = 0.0025*scale;
 247         return 106;
 248       }
 249
 250       errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
 251       erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
 252       return 107;
 253     }
 254     Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
 255     if (diff<1) diff=1;
 256     if (diff>4) diff=4;
 257
 258     if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
 259       errz = 0.14*diff;
 260       erry = 0.003*diff;
 261       return 108;
 262     }
 263     erry = 0.04*diff;
 264     errz = 0.06*diff;
 265     return 109;
 266   }
 267   //DRIFTS
 268   Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
 269   Float_t chargematch = normq/expQ;
 270   chargematch/=2.4; // F. Prino Sept. 2007: SDD charge conversion keV->ADC
 271   Float_t factorz=1;
 272   Int_t   cnz = cl->GetNz()%10;
 273   //charge match
 274   if (cl->GetType()==1){
 275     if (chargematch<1.25){
 276       erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
 277     }
 278     else{
 279       erry = 0.003*chargematch;
 280       if (cl->GetNz()==3) erry*=1.5;
 281     }
 282     if (chargematch<1.0){
 283       errz =  0.0011*(1.+6./cl->GetQ());
 284     }
 285     else{
 286       errz = 0.002*(1+2*(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   if (cl->GetType()>1){
 294     if (chargematch<1){
 295       erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
 296       errz =  0.0016*(1.+6./cl->GetQ());
 297     }
 298     else{
 299       errz = 0.0014*(1+3*(chargematch-1.));
 300       erry = 0.003*(1+3*(chargematch-1.));
 301     }
 302     if (cnz>nz+0.6) {
 303       erry*=(cnz-nz+0.5);
 304       errz*=1.4*(cnz-nz+0.5);
 305     }
 306   }
 307
 308   if (TMath::Abs(cl->GetY())>2.5){
 309     factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
 310   }
 311   if (TMath::Abs(cl->GetY())<1){
 312     factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
 313   }
 314   factorz= TMath::Min(factorz,float(4.));
 315   errz*=factorz;
 316
 317   erry= TMath::Min(erry,float(0.05));
 318   errz= TMath::Min(errz,float(0.05));
 319   return 200;
 320 }
 321 //--------------------------------------------------------------------------
 322 Int_t AliITSClusterParam::GetErrorParamAngle(Int_t layer,
 323                                              const AliITSRecPoint *cl,
 324                                              Float_t tgl,Float_t tgphitr,
 325                                              Float_t &erry,Float_t &errz)
 326 {
 327   //
 328   // Calculate cluster position error (parametrization extracted from rp-hit
 329   // residuals, as a function of angle between track and det module plane.
 330   // Origin: M.Lunardon, S.Moretto)
 331   //
 332
 333   Double_t maxSigmaSPDx=100.;
 334   Double_t maxSigmaSPDz=400.;
 335   Double_t maxSigmaSDDx=100.;
 336   Double_t maxSigmaSDDz=400.;
 337   Double_t maxSigmaSSDx=100.;
 338   Double_t maxSigmaSSDz=1000.;
 339
 340   Double_t paramSPDx[3]={-6.417,0.18,11.14};
 341   Double_t paramSPDz[2]={118.,-0.155};
 342   Double_t paramSDDx[2]={30.93,0.059};
 343   Double_t paramSDDz[2]={33.09,0.011};
 344   Double_t paramSSDx[2]={18.64,-0.0046};
 345   Double_t paramSSDz[2]={784.4,-0.828};
 346   Double_t sigmax=1000.0,sigmaz=1000.0;
 347
 348   Int_t volId = (Int_t)cl->GetVolumeId();
 349   Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
 350   Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
 351
 352
 353   Double_t phitr = TMath::ATan(tgphitr);
 354   Double_t alpha = TMath::ATan2(tra[1],tra[0]);
 355   Double_t phiglob = alpha+phitr;
 356   Double_t p[3];
 357   p[0] = TMath::Cos(phiglob);
 358   p[1] = TMath::Sin(phiglob);
 359   p[2] = tgl;
 360   TVector3 pvec(p[0],p[1],p[2]);
 361   TVector3 normvec(rot[1],rot[4],rot[7]);
 362   Double_t angle = pvec.Angle(normvec);
 363   if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
 364   Double_t angleDeg = angle*180./TMath::Pi();
 365
 366   if(layer==0 || layer==1) { // SPD
 367
 368     sigmax = TMath::Exp(angleDeg*paramSPDx[1]+paramSPDx[0])+paramSPDx[2];
 369     sigmaz = paramSPDz[0]+paramSPDz[1]*angleDeg;
 370     if(sigmax > maxSigmaSPDx) sigmax = maxSigmaSPDx;
 371     if(sigmaz > maxSigmaSPDz) sigmax = maxSigmaSPDz;
 372
 373   } else if(layer==2 || layer==3) { // SDD
 374
 375     sigmax = angleDeg*paramSDDx[1]+paramSDDx[0];
 376     sigmaz = paramSDDz[0]+paramSDDz[1]*angleDeg;
 377     if(sigmax > maxSigmaSDDx) sigmax = maxSigmaSDDx;
 378     if(sigmaz > maxSigmaSDDz) sigmax = maxSigmaSDDz;
 379
 380   } else if(layer==4 || layer==5) { // SSD
 381
 382     sigmax = angleDeg*paramSSDx[1]+paramSSDx[0];
 383     sigmaz = paramSSDz[0]+paramSSDz[1]*angleDeg;
 384     if(sigmax > maxSigmaSSDx) sigmax = maxSigmaSSDx;
 385     if(sigmaz > maxSigmaSSDz) sigmax = maxSigmaSSDz;
 386
 387   }
 388
 389   // convert from micron to cm
 390   erry = 1.e-4*sigmax;
 391   errz = 1.e-4*sigmaz;
 392
 393   return 1;
 394 }
 395 //--------------------------------------------------------------------------
 396 void AliITSClusterParam::Print(Option_t* /*option*/) const {
 397   //
 398   // Print param Information
 399   //
 400
 401   //
 402   // Error parameterization
 403   //
 404   printf("NOT YET...\n");
 405   return;
 406 }
 407
 408
 409
 410
 411