1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // ITS cluster error and shape parameterization //
22 // andrea.dainese@lnl.infn.it //
23 ///////////////////////////////////////////////////////////////////////////////
26 //#include <TVectorF.h>
27 //#include <TLinearFitter.h>
29 //#include <TProfile2D.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"
39 ClassImp(AliITSClusterParam)
42 AliITSClusterParam* AliITSClusterParam::fgInstance = 0;
46 Example usage fitting parameterization (NOT YET...):
47 TFile fres("resol.root"); //tree with resolution and shape
48 TTree * treeRes =(TTree*)fres.Get("Resol");
50 AliITSClusterParam param;
51 param.SetInstance(¶m);
52 param.FitResol(treeRes);
53 param.FitRMS(treeRes);
54 TFile fparam("ITSClusterParam.root","recreate");
58 TFile fparam("ITSClusterParam.root");
59 AliITSClusterParam *param2 = (AliITSClusterParam *) fparam.Get("Param");
60 param2->SetInstance(param2);
61 param2->Test(treeRes);
64 treeRes->Draw("(Resol-AliITSClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
67 //-------------------------------------------------------------------------
68 AliITSClusterParam* AliITSClusterParam::Instance()
71 // Singleton implementation
72 // Returns an instance of this class, it is created if neccessary
75 fgInstance = new AliITSClusterParam();
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)
85 // Get "mean shape" (original parametrization from AliITStrackerMI)
87 tgl = TMath::Abs(tgl);
88 tgphitr = TMath::Abs(tgphitr);
102 if (layer==4 || layer==5) {
103 ny = 2.02+tgphitr*1.95;
104 nz = 2.02+tgphitr*2.35;
108 ny = 6.6-2.7*tgphitr;
109 nz = 2.8-3.11*tgphitr+0.45*tgl;
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,
120 // Calculate cluster position error
124 switch(AliITSReconstructor::GetRecoParam()->GetClusterErrorsParam()) {
126 retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
129 retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
132 retval = GetErrorParamAngle(layer,cl,tgl,tgphitr,erry,errz);
135 retval = GetErrorParamMI(layer,cl,tgl,tgphitr,expQ,erry,errz);
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);
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);
161 //--------------------------------------------------------------------------
162 Int_t AliITSClusterParam::GetErrorOrigRecPoint(const AliITSRecPoint*cl,
163 Float_t &erry,Float_t &errz,Float_t &covyz)
166 // Calculate cluster position error (just take error from AliITSRecPoint)
168 erry = TMath::Sqrt(cl->GetSigmaY2());
169 errz = TMath::Sqrt(cl->GetSigmaZ2());
170 covyz = cl->GetSigmaYZ();
174 //--------------------------------------------------------------------------
175 Int_t AliITSClusterParam::GetErrorParamMI(Int_t layer,const AliITSRecPoint*cl,
176 Float_t tgl,Float_t tgphitr,
178 Float_t &erry,Float_t &errz)
181 // Calculate cluster position error (original parametrization from
185 GetNTeor(layer, cl,tgl,tgphitr,ny,nz);
186 erry = TMath::Sqrt(cl->GetSigmaY2());
187 errz = TMath::Sqrt(cl->GetSigmaZ2());
191 if (TMath::Abs(ny-cl->GetNy())>0.6) {
193 erry*=0.4+TMath::Abs(ny-cl->GetNy());
194 errz*=0.4+TMath::Abs(ny-cl->GetNy());
196 erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
197 errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
200 if (TMath::Abs(nz-cl->GetNz())>1.) {
201 erry*=TMath::Abs(nz-cl->GetNz());
202 errz*=TMath::Abs(nz-cl->GetNz());
206 erry= TMath::Min(erry,float(0.0050));
207 errz= TMath::Min(errz,float(0.0300));
213 //factor 1.8 appears in new simulation
216 if (cl->GetNy()==100||cl->GetNz()==100){
221 if (cl->GetNy()+cl->GetNz()>12){
226 Float_t normq = cl->GetQ()/(TMath::Sqrt(1+tgl*tgl+tgphitr*tgphitr));
227 Float_t chargematch = TMath::Max(double(normq/expQ),2.);
229 if (cl->GetType()==1 || cl->GetType()==10 ){
230 if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
232 erry = 0.00094*scale;
235 if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
241 errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
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;
250 if (cl->GetType()>100 ){
251 if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
253 erry = 0.00096*scale;
256 if (cl->GetNy()+cl->GetNz()-nz-ny<1){
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;
266 Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
270 if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
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
284 Int_t cnz = cl->GetNz()%10;
286 if (cl->GetType()==1){
287 if (chargematch<1.25){
288 erry = 0.0028*(1.+6./cl->GetQ()); // gold clusters
291 erry = 0.003*chargematch;
292 if (cl->GetNz()==3) erry*=1.5;
294 if (chargematch<1.0){
295 errz = 0.0011*(1.+6./cl->GetQ());
298 errz = 0.002*(1+2*(chargematch-1.));
302 errz*=1.4*(cnz-nz+0.5);
305 if (cl->GetType()>1){
307 erry = 0.00385*(1.+6./cl->GetQ()); // gold clusters
308 errz = 0.0016*(1.+6./cl->GetQ());
311 errz = 0.0014*(1+3*(chargematch-1.));
312 erry = 0.003*(1+3*(chargematch-1.));
316 errz*=1.4*(cnz-nz+0.5);
320 if (TMath::Abs(cl->GetY())>2.5){
321 factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
323 if (TMath::Abs(cl->GetY())<1){
324 factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
326 factorz= TMath::Min(factorz,float(4.));
329 erry= TMath::Min(erry,float(0.05));
330 errz= TMath::Min(errz,float(0.05));
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)
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)
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.;
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;
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);
365 Double_t phitr = TMath::ATan(tgphitr);
366 Double_t alpha = TMath::ATan2(tra[1],tra[0]);
367 Double_t phiglob = alpha+phitr;
369 p[0] = TMath::Cos(phiglob);
370 p[1] = TMath::Sin(phiglob);
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();
378 if(layer==0 || layer==1) { // SPD
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;
385 } else if(layer==2 || layer==3) { // SDD
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;
392 } else if(layer==4 || layer==5) { // SSD
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;
401 // convert from micron to cm
407 //--------------------------------------------------------------------------
408 void AliITSClusterParam::Print(Option_t* /*option*/) const {
410 // Print param Information
414 // Error parameterization
416 printf("NOT YET...\n");