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