remove props
[u/mrichter/AliRoot.git] / 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 #include "AliCheb3DCalc.h"
39
40 ClassImp(AliITSClusterParam)
41
42
43 AliITSClusterParam* AliITSClusterParam::fgInstance = 0;
44
45
46 /*
47   Example usage fitting parameterization (NOT YET...):
48   TFile fres("resol.root");    //tree with resolution and shape 
49   TTree * treeRes =(TTree*)fres.Get("Resol");
50   
51   AliITSClusterParam param;
52   param.SetInstance(&param);
53   param.FitResol(treeRes);
54   param.FitRMS(treeRes);
55   TFile fparam("ITSClusterParam.root","recreate");
56   param.Write("Param");
57   //
58   //
59   TFile fparam("ITSClusterParam.root");
60   AliITSClusterParam *param2  =  (AliITSClusterParam *) fparam.Get("Param"); 
61   param2->SetInstance(param2);
62   param2->Test(treeRes);
63   
64
65   treeRes->Draw("(Resol-AliITSClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
66 */
67
68 //-------------------------------------------------------------------------
69 AliITSClusterParam* AliITSClusterParam::Instance()
70 {
71   //
72   // Singleton implementation
73   // Returns an instance of this class, it is created if neccessary
74   //
75   if (fgInstance == 0){
76     fgInstance = new AliITSClusterParam();
77   }
78   return fgInstance;
79 }
80 //-------------------------------------------------------------------------
81 void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /*cl*/, 
82                                   Float_t tgl,Float_t tgphitr,
83                                   Float_t &ny,Float_t &nz)
84 {
85   //
86   // Get "mean shape" (original parametrization from AliITStrackerMI)
87   //
88   tgl = TMath::Abs(tgl);
89   tgphitr = TMath::Abs(tgphitr);
90
91   // SPD
92   if (layer==0) {
93     ny = 1.+tgphitr*3.2;
94     nz = 1.+tgl*0.34;
95     return;
96   }
97   if (layer==1) {
98     ny = 1.+tgphitr*3.2;
99     nz = 1.+tgl*0.28;
100     return;
101   }
102   // SSD
103   if (layer==4 || layer==5) {
104     ny = 2.02+tgphitr*1.95;
105     nz = 2.02+tgphitr*2.35;
106     return;
107   }
108   // SDD
109   ny  = 6.6-2.7*tgphitr;
110   nz  = 2.8-3.11*tgphitr+0.45*tgl;
111   return;
112 }
113 //--------------------------------------------------------------------------
114 Int_t AliITSClusterParam::GetError(Int_t layer,
115                                    const AliITSRecPoint *cl,
116                                    Float_t tgl,Float_t tgphitr,Float_t expQ,
117                                    Float_t &erry,Float_t &errz,Float_t &covyz,
118                                    Bool_t addMisalErr)
119 {
120   //
121   // Calculate cluster position error
122   //
123   Int_t retval=0;
124   covyz=0.;
125   switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
126   case 0: 
127     retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
128     break;
129   case 1: 
130     retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
131     break;
132   case 2: 
133     retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
134     break;
135   case 3: 
136     retval = GetErrorParamAngleOld(layer,cl,tgl,tgphitr,erry,errz);
137     break;
138   default: 
139     retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
140     break;
141   }
142
143   // for SSD use the error provided by the cluster finder 
144   // if single-sided clusters are enabled
145   if(layer>=4 && AliITSReconstructor::GetRecoParam()->GetUseBadChannelsInClusterFinderSSD()) { 
146     //printf("error 1 erry errz covyz %10.7f %10.7f %15.13f\n",erry,errz,covyz);
147     retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
148     //printf("type %d erry errz covyz %10.7f %10.7f %15.13f\n",cl->GetType(),erry,errz,covyz);
149   }
150   
151   if(addMisalErr) {
152     Double_t bz = (Double_t)AliTracker::GetBz();
153     // add error due to misalignment (to be improved)
154     Float_t errmisalY2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer,bz)
155       *AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorY(layer,bz);
156     Float_t errmisalZ2 = AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer,bz)
157       *AliITSReconstructor::GetRecoParam()->GetClusterMisalErrorZ(layer,bz);
158     erry = TMath::Sqrt(erry*erry+errmisalY2);
159     errz = TMath::Sqrt(errz*errz+errmisalZ2);
160   }
161
162   return retval;
163
164 }
165 //--------------------------------------------------------------------------
166 Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
167                                    Float_t &erry,Float_t &errz,Float_t &covyz)
168 {
169   //
170   // Calculate cluster position error (just take error from AliITSRecPoint)
171   //
172   erry   = TMath::Sqrt(cl->GetSigmaY2()); 
173   errz   = TMath::Sqrt(cl->GetSigmaZ2()); 
174   covyz  = cl->GetSigmaYZ();
175
176   return 1;
177 }
178 //--------------------------------------------------------------------------
179 Int_t AliITSClusterParam::GetErrorParamMI(Int_t layer,const AliITSRecPoint*cl,
180                                           Float_t tgl,Float_t tgphitr,
181                                           Float_t expQ,
182                                           Float_t &erry,Float_t &errz)
183 {
184   //
185   // Calculate cluster position error (original parametrization from 
186   // AliITStrackerMI)
187   //
188   Float_t nz,ny;
189   GetNTeor(layer, cl,tgl,tgphitr,ny,nz);  
190   erry   = TMath::Sqrt(cl->GetSigmaY2()); 
191   errz   = TMath::Sqrt(cl->GetSigmaZ2()); 
192   //
193   // PIXELS
194   if (layer<2){
195     if (TMath::Abs(ny-cl->GetNy())>0.6)  {
196       if (ny<cl->GetNy()){
197         erry*=0.4+TMath::Abs(ny-cl->GetNy());
198         errz*=0.4+TMath::Abs(ny-cl->GetNy());
199       }else{
200         erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
201         errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
202       }
203     }
204     if (TMath::Abs(nz-cl->GetNz())>1.)  {
205       erry*=TMath::Abs(nz-cl->GetNz());
206       errz*=TMath::Abs(nz-cl->GetNz());       
207     }
208     erry*=0.85;
209     errz*=0.85;
210     erry= TMath::Min(erry,float(0.0050));
211     errz= TMath::Min(errz,float(0.0300));
212     return 10;
213   }
214
215   //STRIPS
216   if (layer>3){ 
217     //factor 1.8 appears in new simulation
218     //
219     Float_t scale=1.8;
220     if (cl->GetNy()==100||cl->GetNz()==100){
221       erry = 0.004*scale;
222       errz = 0.2*scale;
223       return 100;
224     }
225     if (cl->GetNy()+cl->GetNz()>12){
226       erry = 0.06*scale;
227       errz = 0.57*scale;
228       return 100;
229     }
230     Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
231     Float_t chargematch = TMath::Max(double(normq/expQ),2.);
232     //
233     if (cl->GetType()==1 || cl->GetType()==10 ){                                                               
234       if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
235         errz = 0.043*scale;
236         erry = 0.00094*scale;
237         return 101;
238       }
239       if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
240         errz = 0.06*scale;
241         erry =0.0013*scale;
242         return 102;
243       }
244       erry = 0.0027*scale;
245       errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
246       return 103;
247     }
248     if (cl->GetType()==2 || cl->GetType()==11 ){ 
249       erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
250       errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
251       return 104;
252     }
253     
254     if (cl->GetType()>100 ){                                                                   
255       if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
256         errz = 0.05*scale;
257         erry = 0.00096*scale;
258         return 105;
259       }
260       if (cl->GetNy()+cl->GetNz()-nz-ny<1){
261         errz = 0.10*scale;
262         erry = 0.0025*scale;
263         return 106;
264       }
265
266       errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
267       erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
268       return 107;
269     }    
270     Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
271     if (diff<1) diff=1;
272     if (diff>4) diff=4;
273         
274     if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
275       errz = 0.14*diff;
276       erry = 0.003*diff;
277       return 108;
278     }  
279     erry = 0.04*diff;
280     errz = 0.06*diff;
281     return 109;
282   }
283   //DRIFTS
284   Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
285   Float_t chargematch = normq/expQ;
286   chargematch/=2.4; // F. Prino Sept. 2007: SDD charge conversion keV->ADC
287   Float_t factorz=1;
288   Int_t   cnz = cl->GetNz()%10;
289   //charge match
290   if (cl->GetType()==1){
291     if (chargematch<1.25){
292       erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
293     }
294     else{
295       erry = 0.003*chargematch;
296       if (cl->GetNz()==3) erry*=1.5;
297     }
298     if (chargematch<1.0){
299       errz =  0.0011*(1.+6./cl->GetQ());
300     }
301     else{
302       errz = 0.002*(1+2*(chargematch-1.));
303     }
304     if (cnz>nz+0.6) {
305       erry*=(cnz-nz+0.5);
306       errz*=1.4*(cnz-nz+0.5);
307     }
308   }
309   if (cl->GetType()>1){
310     if (chargematch<1){
311       erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
312       errz =  0.0016*(1.+6./cl->GetQ());
313     }
314     else{
315       errz = 0.0014*(1+3*(chargematch-1.));
316       erry = 0.003*(1+3*(chargematch-1.));
317     } 
318     if (cnz>nz+0.6) {
319       erry*=(cnz-nz+0.5);
320       errz*=1.4*(cnz-nz+0.5);
321     }
322   }
323
324   if (TMath::Abs(cl->GetY())>2.5){
325     factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
326   }
327   if (TMath::Abs(cl->GetY())<1){
328     factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
329   }
330   factorz= TMath::Min(factorz,float(4.));  
331   errz*=factorz;
332
333   erry= TMath::Min(erry,float(0.05));
334   errz= TMath::Min(errz,float(0.05));  
335   return 200;
336 }
337 //--------------------------------------------------------------------------
338 Int_t AliITSClusterParam::GetErrorParamAngle(Int_t layer,
339                                              const AliITSRecPoint *cl,
340                                              Float_t tgl,Float_t tgphitr,
341                                              Float_t &erry,Float_t &errz)
342 {
343   //
344   // Calculate cluster position error (parametrization extracted from rp-hit
345   // residuals, as a function of angle between track and det module plane.
346   // Origin: M.Lunardon, S.Moretto)
347   //
348   const int kNcfSPDResX = 21;
349   float kCfSPDResX[kNcfSPDResX] = {+1.1201e+01,+2.0903e+00,-2.2909e-01,-2.6413e-01,+4.2135e-01,-3.7190e-01,
350                             +4.2339e-01,+1.8679e-01,-5.1249e-01,+1.8421e-01,+4.8849e-02,-4.3127e-01,
351                             -1.1148e-01,+3.1984e-03,-2.5743e-01,-6.6408e-02,+3.0756e-01,+2.6809e-01,
352                             -5.0339e-03,-1.4964e-01,-1.1001e-01};
353   const float kSPDazMax=56.000000;
354   //
355   const int kNcfSPDMeanX = 16;
356   float kCfSPDMeanX[kNcfSPDMeanX] = {-1.2532e+00,-3.8185e-01,-8.9039e-01,+2.6648e+00,+7.0361e-01,+1.2298e+00,
357                                      +3.2871e-01,+7.8487e-02,-1.6792e-01,-1.3966e-01,-3.1670e-01,-2.1795e-01,
358                                      -1.9451e-01,-4.9347e-02,-1.9186e-01,-1.9195e-01};
359   //
360   const int kNcfSPDResZ = 5;
361   float kCfSPDResZ[kNcfSPDResZ] = {+9.2384e+01,+3.4352e-01,-2.7317e+01,-1.4642e-01,+2.0868e+00};
362   const float kSPDpolMin=34.358002, kSPDpolMax=145.000000;
363   //
364   Double_t maxSigmaSDDx=100.;
365   Double_t maxSigmaSDDz=400.;
366   Double_t maxSigmaSSDx=100.;
367   Double_t maxSigmaSSDz=1000.;
368   //  
369   Double_t paramSDDx[2]={30.93,0.059};
370   Double_t paramSDDz[2]={33.09,0.011};
371   Double_t paramSSDx[2]={18.64,-0.0046};
372   Double_t paramSSDz[2]={784.4,-0.828};
373   Double_t sigmax=1000.0,sigmaz=1000.0;
374   Double_t biasx = 0.0;
375
376   Int_t volId = (Int_t)cl->GetVolumeId();
377   Double_t rotMA[9]; AliGeomManager::GetRotation(volId,rotMA);      // misaligned rotation
378   Double_t rotOR[9]; AliGeomManager::GetOrigRotation(volId,rotOR);  // original rotation
379   // difference in phi of original and misaligned sensors
380   double cross = rotOR[1]*rotMA[4]-rotOR[4]*rotMA[1];
381   cross /= TMath::Sqrt( (1.-rotOR[7]*rotOR[7]) * (1.-rotMA[7]*rotMA[7]) );
382   Double_t angleAzi = TMath::Abs(TMath::ATan(tgphitr) - TMath::ASin(cross) );
383   Double_t anglePol = TMath::Abs(TMath::ATan(tgl));
384
385   if(angleAzi>0.5*TMath::Pi()) angleAzi = TMath::Pi()-angleAzi;
386   if(anglePol>0.5*TMath::Pi()) anglePol = TMath::Pi()-anglePol;
387   Double_t angleAziDeg = angleAzi*180./TMath::Pi();
388   Double_t anglePolDeg = anglePol*180./TMath::Pi();
389   
390   if(layer==0 || layer==1) { // SPD
391     //
392     float phiInt    = angleAziDeg/kSPDazMax; // mapped to -1:1
393     if (phiInt>1) phiInt = 1; else if (phiInt<-1) phiInt = -1;
394     float phiAbsInt = (TMath::Abs(angleAziDeg+angleAziDeg) - kSPDazMax)/kSPDazMax; // mapped to -1:1
395     if (phiAbsInt>1) phiAbsInt = 1; else if (phiAbsInt<-1) phiAbsInt = -1;
396     anglePolDeg += 90; // the parameterization was provided in polar angle (90 deg - normal to sensor)
397     float polInt   = (anglePolDeg+anglePolDeg - (kSPDpolMax+kSPDpolMin))/(kSPDpolMax-kSPDpolMin); // mapped to -1:1
398     if (polInt>1) polInt = 1; else if (polInt<-1) polInt = -1;
399     //
400     sigmax = AliCheb3DCalc::ChebEval1D(phiAbsInt, kCfSPDResX , kNcfSPDResX);
401     biasx  = AliCheb3DCalc::ChebEval1D(phiInt   , kCfSPDMeanX, kNcfSPDMeanX);
402     sigmaz = AliCheb3DCalc::ChebEval1D(polInt   , kCfSPDResZ , kNcfSPDResZ);
403     //
404     // for the moment for the SPD only, need to decide where to put it
405     biasx *= 1e-4;
406     
407   } else if(layer==2 || layer==3) { // SDD
408
409     sigmax = angleAziDeg*paramSDDx[1]+paramSDDx[0];
410     sigmaz = paramSDDz[0]+paramSDDz[1]*anglePolDeg;
411     if(sigmax > maxSigmaSDDx) sigmax = maxSigmaSDDx;
412     if(sigmaz > maxSigmaSDDz) sigmax = maxSigmaSDDz;
413     
414   } else if(layer==4 || layer==5) { // SSD
415
416     sigmax = angleAziDeg*paramSSDx[1]+paramSSDx[0];
417     sigmaz = paramSSDz[0]+paramSSDz[1]*anglePolDeg;
418     if(sigmax > maxSigmaSSDx) sigmax = maxSigmaSSDx;
419     if(sigmaz > maxSigmaSSDz) sigmax = maxSigmaSSDz;
420     
421   }
422
423   // convert from micron to cm
424   erry = 1.e-4*sigmax; 
425   errz = 1.e-4*sigmaz;
426   
427
428   return 1;
429 }
430 //--------------------------------------------------------------------------
431 Int_t AliITSClusterParam::GetErrorParamAngleOld(Int_t layer,
432                                                 const AliITSRecPoint *cl,
433                                                 Float_t tgl,Float_t tgphitr,
434                                                 Float_t &erry,Float_t &errz)
435 {
436   //
437   // Calculate cluster position error (parametrization extracted from rp-hit
438   // residuals, as a function of angle between track and det module plane.
439   // Origin: M.Lunardon, S.Moretto)
440   //
441
442   Double_t maxSigmaSPDx=100.;
443   Double_t maxSigmaSPDz=400.;
444   Double_t maxSigmaSDDx=100.;
445   Double_t maxSigmaSDDz=400.;
446   Double_t maxSigmaSSDx=100.;
447   Double_t maxSigmaSSDz=1000.;
448   
449   Double_t paramSPDx[3]={-6.417,0.18,11.14};
450   Double_t paramSPDz[2]={118.,-0.155};
451   Double_t paramSDDx[2]={30.93,0.059};
452   Double_t paramSDDz[2]={33.09,0.011};
453   Double_t paramSSDx[2]={18.64,-0.0046};
454   Double_t paramSSDz[2]={784.4,-0.828};
455   Double_t sigmax=1000.0,sigmaz=1000.0;
456   
457   Int_t volId = (Int_t)cl->GetVolumeId();
458   Double_t rotMA[9]; AliGeomManager::GetRotation(volId,rotMA);      // misaligned rotation
459   Double_t rotOR[9]; AliGeomManager::GetOrigRotation(volId,rotOR);  // original rotation
460   // difference in phi of original and misaligned sensors
461   double cross = rotOR[1]*rotMA[4]-rotOR[4]*rotMA[1];
462   cross /= TMath::Sqrt( (1.-rotOR[7]*rotOR[7]) * (1.-rotMA[7]*rotMA[7]) );
463   Double_t angleAzi = TMath::Abs(TMath::ATan(tgphitr) - TMath::ASin(cross) );
464   Double_t anglePol = TMath::Abs(TMath::ATan(tgl));
465
466   if(angleAzi>0.5*TMath::Pi()) angleAzi = TMath::Pi()-angleAzi;
467   if(anglePol>0.5*TMath::Pi()) anglePol = TMath::Pi()-anglePol;
468   Double_t angleAziDeg = angleAzi*180./TMath::Pi();
469   Double_t anglePolDeg = anglePol*180./TMath::Pi();
470   
471   if(layer==0 || layer==1) { // SPD
472
473     sigmax = TMath::Exp(angleAziDeg*paramSPDx[1]+paramSPDx[0])+paramSPDx[2];
474     sigmaz = paramSPDz[0]+paramSPDz[1]*anglePolDeg;
475     if(sigmax > maxSigmaSPDx) sigmax = maxSigmaSPDx;
476     if(sigmaz > maxSigmaSPDz) sigmax = maxSigmaSPDz;
477
478   } else if(layer==2 || layer==3) { // SDD
479
480     sigmax = angleAziDeg*paramSDDx[1]+paramSDDx[0];
481     sigmaz = paramSDDz[0]+paramSDDz[1]*anglePolDeg;
482     if(sigmax > maxSigmaSDDx) sigmax = maxSigmaSDDx;
483     if(sigmaz > maxSigmaSDDz) sigmax = maxSigmaSDDz;
484     
485   } else if(layer==4 || layer==5) { // SSD
486
487     sigmax = angleAziDeg*paramSSDx[1]+paramSSDx[0];
488     sigmaz = paramSSDz[0]+paramSSDz[1]*anglePolDeg;
489     if(sigmax > maxSigmaSSDx) sigmax = maxSigmaSSDx;
490     if(sigmaz > maxSigmaSSDz) sigmax = maxSigmaSSDz;
491     
492   }
493
494   // convert from micron to cm
495   erry = 1.e-4*sigmax; 
496   errz = 1.e-4*sigmaz;
497   
498   return 1;
499 }
500 //--------------------------------------------------------------------------
501 void AliITSClusterParam::Print(Option_t* /*option*/) const {
502   //
503   // Print param Information
504   //
505
506   //
507   // Error parameterization
508   //
509   printf("NOT YET...\n");
510   return;
511 }
512
513
514
515
516