Modifications to use 1D SDD clusters for tracking (A. Dainese)
[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
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