fixed the tainted variables
[u/mrichter/AliRoot.git] / ITS / AliITSClusterParam.cxx
CommitLineData
572f41f9 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>
e50912db 30#include <TVector3.h>
572f41f9 31#include "TMath.h"
d9ead1a0 32#include "AliTracker.h"
e50912db 33#include "AliGeomManager.h"
572f41f9 34#include "AliITSRecPoint.h"
35#include "AliITSClusterParam.h"
e50912db 36#include "AliITSReconstructor.h"
37#include "AliExternalTrackParam.h"
572f41f9 38
39ClassImp(AliITSClusterParam)
40
41
42AliITSClusterParam* 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//-------------------------------------------------------------------------
68AliITSClusterParam* 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//-------------------------------------------------------------------------
80void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /*cl*/,
e50912db 81 Float_t tgl,Float_t tgphitr,
572f41f9 82 Float_t &ny,Float_t &nz)
83{
84 //
e50912db 85 // Get "mean shape" (original parametrization from AliITStrackerMI)
572f41f9 86 //
e50912db 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;
572f41f9 94 return;
95 }
e50912db 96 if (layer==1) {
97 ny = 1.+tgphitr*3.2;
98 nz = 1.+tgl*0.28;
572f41f9 99 return;
100 }
e50912db 101 // SSD
102 if (layer==4 || layer==5) {
103 ny = 2.02+tgphitr*1.95;
104 nz = 2.02+tgphitr*2.35;
572f41f9 105 return;
106 }
e50912db 107 // SDD
108 ny = 6.6-2.7*tgphitr;
109 nz = 2.8-3.11*tgphitr+0.45*tgl;
110 return;
111}
112//--------------------------------------------------------------------------
113Int_t AliITSClusterParam::GetError(Int_t layer,
114 const AliITSRecPoint *cl,
115 Float_t tgl,Float_t tgphitr,Float_t expQ,
d9ead1a0 116 Float_t &erry,Float_t &errz,Float_t &covyz,
8c139cf3 117 Bool_t addMisalErr)
e50912db 118{
119 //
120 // Calculate cluster position error
121 //
401eff16 122 Int_t retval=0;
d9ead1a0 123 covyz=0.;
e50912db 124 switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
125 case 0:
d9ead1a0 126 retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
e50912db 127 break;
128 case 1:
401eff16 129 retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
e50912db 130 break;
131 case 2:
401eff16 132 retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
e50912db 133 break;
134 default:
401eff16 135 retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
e50912db 136 break;
137 }
138
d9ead1a0 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
8c139cf3 147 if(addMisalErr) {
d9ead1a0 148 Double_t bz = (Double_t)AliTracker::GetBz();
8c139cf3 149 // add error due to misalignment (to be improved)
d9ead1a0 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);
f9119eb9 154 erry = TMath::Sqrt(erry*erry+errmisalY2);
155 errz = TMath::Sqrt(errz*errz+errmisalZ2);
8c139cf3 156 }
401eff16 157
158 return retval;
159
572f41f9 160}
161//--------------------------------------------------------------------------
e50912db 162Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
d9ead1a0 163 Float_t &erry,Float_t &errz,Float_t &covyz)
572f41f9 164{
e50912db 165 //
166 // Calculate cluster position error (just take error from AliITSRecPoint)
167 //
168 erry = TMath::Sqrt(cl->GetSigmaY2());
169 errz = TMath::Sqrt(cl->GetSigmaZ2());
d9ead1a0 170 covyz = cl->GetSigmaYZ();
e50912db 171
172 return 1;
173}
174//--------------------------------------------------------------------------
175Int_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)
572f41f9 183 //
184 Float_t nz,ny;
e50912db 185 GetNTeor(layer, cl,tgl,tgphitr,ny,nz);
572f41f9 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 }
e50912db 226 Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
572f41f9 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
e50912db 280 Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
572f41f9 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//--------------------------------------------------------------------------
e50912db 334Int_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();
242bd643 361 Double_t rotMA[9]; AliGeomManager::GetRotation(volId,rotMA); // misaligned rotation
362 Double_t rotOR[9]; AliGeomManager::GetOrigRotation(volId,rotOR); // original rotation
363 // difference in phi of original and misaligned sensors
364 double cross = rotOR[1]*rotMA[4]-rotOR[4]*rotMA[1];
365 cross /= TMath::Sqrt( (1.-rotOR[7]*rotOR[7]) * (1.-rotMA[7]*rotMA[7]) );
366 Double_t angleAzi = TMath::Abs(TMath::ATan(tgphitr) - TMath::ASin(cross) );
196b395e 367 Double_t anglePol = TMath::Abs(TMath::ATan(tgl));
242bd643 368
196b395e 369 if(angleAzi>0.5*TMath::Pi()) angleAzi = TMath::Pi()-angleAzi;
370 if(anglePol>0.5*TMath::Pi()) anglePol = TMath::Pi()-anglePol;
371 Double_t angleAziDeg = angleAzi*180./TMath::Pi();
372 Double_t anglePolDeg = anglePol*180./TMath::Pi();
e50912db 373
374 if(layer==0 || layer==1) { // SPD
375
196b395e 376 sigmax = TMath::Exp(angleAziDeg*paramSPDx[1]+paramSPDx[0])+paramSPDx[2];
377 sigmaz = paramSPDz[0]+paramSPDz[1]*anglePolDeg;
e50912db 378 if(sigmax > maxSigmaSPDx) sigmax = maxSigmaSPDx;
379 if(sigmaz > maxSigmaSPDz) sigmax = maxSigmaSPDz;
380
381 } else if(layer==2 || layer==3) { // SDD
382
196b395e 383 sigmax = angleAziDeg*paramSDDx[1]+paramSDDx[0];
384 sigmaz = paramSDDz[0]+paramSDDz[1]*anglePolDeg;
e50912db 385 if(sigmax > maxSigmaSDDx) sigmax = maxSigmaSDDx;
386 if(sigmaz > maxSigmaSDDz) sigmax = maxSigmaSDDz;
387
388 } else if(layer==4 || layer==5) { // SSD
389
196b395e 390 sigmax = angleAziDeg*paramSSDx[1]+paramSSDx[0];
391 sigmaz = paramSSDz[0]+paramSSDz[1]*anglePolDeg;
e50912db 392 if(sigmax > maxSigmaSSDx) sigmax = maxSigmaSSDx;
393 if(sigmaz > maxSigmaSSDz) sigmax = maxSigmaSSDz;
394
395 }
396
397 // convert from micron to cm
398 erry = 1.e-4*sigmax;
399 errz = 1.e-4*sigmaz;
400
401 return 1;
402}
403//--------------------------------------------------------------------------
572f41f9 404void AliITSClusterParam::Print(Option_t* /*option*/) const {
405 //
406 // Print param Information
407 //
408
409 //
410 // Error parameterization
411 //
412 printf("NOT YET...\n");
413 return;
414}
415
416
417
418
419