]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSClusterParam.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
some histos added for TPC clusters
[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"
e50912db 32#include "AliGeomManager.h"
572f41f9 33#include "AliITSRecPoint.h"
34#include "AliITSClusterParam.h"
e50912db 35#include "AliITSReconstructor.h"
36#include "AliExternalTrackParam.h"
572f41f9 37
38ClassImp(AliITSClusterParam)
39
40
41AliITSClusterParam* 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//-------------------------------------------------------------------------
67AliITSClusterParam* 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//-------------------------------------------------------------------------
79void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /*cl*/,
e50912db 80 Float_t tgl,Float_t tgphitr,
572f41f9 81 Float_t &ny,Float_t &nz)
82{
83 //
e50912db 84 // Get "mean shape" (original parametrization from AliITStrackerMI)
572f41f9 85 //
e50912db 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;
572f41f9 93 return;
94 }
e50912db 95 if (layer==1) {
96 ny = 1.+tgphitr*3.2;
97 nz = 1.+tgl*0.28;
572f41f9 98 return;
99 }
e50912db 100 // SSD
101 if (layer==4 || layer==5) {
102 ny = 2.02+tgphitr*1.95;
103 nz = 2.02+tgphitr*2.35;
572f41f9 104 return;
105 }
e50912db 106 // SDD
107 ny = 6.6-2.7*tgphitr;
108 nz = 2.8-3.11*tgphitr+0.45*tgl;
109 return;
110}
111//--------------------------------------------------------------------------
112Int_t AliITSClusterParam::GetError(Int_t layer,
113 const AliITSRecPoint *cl,
114 Float_t tgl,Float_t tgphitr,Float_t expQ,
8c139cf3 115 Float_t &erry,Float_t &errz,
116 Bool_t addMisalErr)
e50912db 117{
118 //
119 // Calculate cluster position error
120 //
401eff16 121 Int_t retval=0;
e50912db 122 switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
123 case 0:
401eff16 124 retval = GetErrorOrigRecPoint(cl,erry,errz);
e50912db 125 break;
126 case 1:
401eff16 127 retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
e50912db 128 break;
129 case 2:
401eff16 130 retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
e50912db 131 break;
132 default:
401eff16 133 retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
e50912db 134 break;
135 }
136
8c139cf3 137 if(addMisalErr) {
138 // add error due to misalignment (to be improved)
f9119eb9 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);
8c139cf3 145 }
401eff16 146
147 return retval;
148
572f41f9 149}
150//--------------------------------------------------------------------------
e50912db 151Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
572f41f9 152 Float_t &erry,Float_t &errz)
153{
e50912db 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//--------------------------------------------------------------------------
163Int_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)
572f41f9 171 //
172 Float_t nz,ny;
e50912db 173 GetNTeor(layer, cl,tgl,tgphitr,ny,nz);
572f41f9 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 }
e50912db 214 Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
572f41f9 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
e50912db 268 Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
572f41f9 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//--------------------------------------------------------------------------
e50912db 322Int_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]);
c3a24b56 361 TVector3 normvec(rot[1],rot[4],rot[7]);
e50912db 362 Double_t angle = pvec.Angle(normvec);
363 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
0ffe9832 364 Double_t angleDeg = angle*180./TMath::Pi();
e50912db 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//--------------------------------------------------------------------------
572f41f9 396void 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