[u/mrichter/AliRoot.git] / ITS / AliITSClusterParam.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  ITS cluster error and shape parameterization                             //
//                                                                           //
//  andrea.dainese@lnl.infn.it                                               //
///////////////////////////////////////////////////////////////////////////////
//#include "TFile.h"
//#include "TTree.h"
//#include <TVectorF.h>
//#include <TLinearFitter.h>
//#include <TH1F.h>
//#include <TProfile2D.h>
#include "TMath.h"
#include "AliITSRecPoint.h"
#include "AliITSClusterParam.h"

ClassImp(AliITSClusterParam)


AliITSClusterParam* AliITSClusterParam::fgInstance = 0;


/*
  Example usage fitting parameterization (NOT YET...):
  TFile fres("resol.root");    //tree with resolution and shape 
  TTree * treeRes =(TTree*)fres.Get("Resol");
  
  AliITSClusterParam param;
  param.SetInstance(&param);
  param.FitResol(treeRes);
  param.FitRMS(treeRes);
  TFile fparam("ITSClusterParam.root","recreate");
  param.Write("Param");
  //
  //
  TFile fparam("ITSClusterParam.root");
  AliITSClusterParam *param2  =  (AliITSClusterParam *) fparam.Get("Param"); 
  param2->SetInstance(param2);
  param2->Test(treeRes);
  

  treeRes->Draw("(Resol-AliITSClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
*/

//-------------------------------------------------------------------------
AliITSClusterParam* AliITSClusterParam::Instance()
{
  //
  // Singleton implementation
  // Returns an instance of this class, it is created if neccessary
  //
  if (fgInstance == 0){
    fgInstance = new AliITSClusterParam();
  }
  return fgInstance;
}
//-------------------------------------------------------------------------
void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /*cl*/, 
				  Float_t theta,Float_t phi,
				  Float_t &ny,Float_t &nz)
{
  //
  //get "mean shape"
  //
  if (layer==0){
    ny = 1.+TMath::Abs(phi)*3.2;
    nz = 1.+TMath::Abs(theta)*0.34;
    return;
  }
  if (layer==1){
    ny = 1.+TMath::Abs(phi)*3.2;
    nz = 1.+TMath::Abs(theta)*0.28;
    return;
  }
  
  if (layer>3){
    ny = 2.02+TMath::Abs(phi)*1.95;
    nz = 2.02+TMath::Abs(phi)*2.35;
    return;
  }
  ny  = 6.6-2.7*TMath::Abs(phi);
  nz  = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
}
//--------------------------------------------------------------------------
Int_t AliITSClusterParam::GetError(Int_t layer,const AliITSRecPoint*cl,
				   Float_t theta,Float_t phi,Float_t expQ,
				   Float_t &erry,Float_t &errz)
{
  //calculate cluster position error
  //
  Float_t nz,ny;
  GetNTeor(layer, cl,theta,phi,ny,nz);  
  erry   = TMath::Sqrt(cl->GetSigmaY2()); 
  errz   = TMath::Sqrt(cl->GetSigmaZ2()); 
  //
  // PIXELS
  if (layer<2){
    if (TMath::Abs(ny-cl->GetNy())>0.6)  {
      if (ny<cl->GetNy()){
	erry*=0.4+TMath::Abs(ny-cl->GetNy());
	errz*=0.4+TMath::Abs(ny-cl->GetNy());
      }else{
	erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
	errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
      }
    }
    if (TMath::Abs(nz-cl->GetNz())>1.)  {
      erry*=TMath::Abs(nz-cl->GetNz());
      errz*=TMath::Abs(nz-cl->GetNz());	      
    }
    erry*=0.85;
    errz*=0.85;
    erry= TMath::Min(erry,float(0.0050));
    errz= TMath::Min(errz,float(0.0300));
    return 10;
  }

  //STRIPS
  if (layer>3){ 
    //factor 1.8 appears in new simulation
    //
    Float_t scale=1.8;
    if (cl->GetNy()==100||cl->GetNz()==100){
      erry = 0.004*scale;
      errz = 0.2*scale;
      return 100;
    }
    if (cl->GetNy()+cl->GetNz()>12){
      erry = 0.06*scale;
      errz = 0.57*scale;
      return 100;
    }
    Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
    Float_t chargematch = TMath::Max(double(normq/expQ),2.);
    //
    if (cl->GetType()==1 || cl->GetType()==10 ){     							       
      if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
	errz = 0.043*scale;
	erry = 0.00094*scale;
	return 101;
      }
      if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
	errz = 0.06*scale;
	erry =0.0013*scale;
	return 102;
      }
      erry = 0.0027*scale;
      errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
      return 103;
    }
    if (cl->GetType()==2 || cl->GetType()==11 ){ 
      erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
      errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
      return 104;
    }
    
    if (cl->GetType()>100 ){     							       
      if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
	errz = 0.05*scale;
	erry = 0.00096*scale;
	return 105;
      }
      if (cl->GetNy()+cl->GetNz()-nz-ny<1){
	errz = 0.10*scale;
	erry = 0.0025*scale;
	return 106;
      }

      errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
      erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
      return 107;
    }    
    Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
    if (diff<1) diff=1;
    if (diff>4) diff=4;
        
    if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
      errz = 0.14*diff;
      erry = 0.003*diff;
      return 108;
    }  
    erry = 0.04*diff;
    errz = 0.06*diff;
    return 109;
  }
  //DRIFTS
  Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
  Float_t chargematch = normq/expQ;
  chargematch/=2.4; // F. Prino Sept. 2007: SDD charge conversion keV->ADC
  Float_t factorz=1;
  Int_t   cnz = cl->GetNz()%10;
  //charge match
  if (cl->GetType()==1){
    if (chargematch<1.25){
      erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
    }
    else{
      erry = 0.003*chargematch;
      if (cl->GetNz()==3) erry*=1.5;
    }
    if (chargematch<1.0){
      errz =  0.0011*(1.+6./cl->GetQ());
    }
    else{
      errz = 0.002*(1+2*(chargematch-1.));
    }
    if (cnz>nz+0.6) {
      erry*=(cnz-nz+0.5);
      errz*=1.4*(cnz-nz+0.5);
    }
  }
  if (cl->GetType()>1){
    if (chargematch<1){
      erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
      errz =  0.0016*(1.+6./cl->GetQ());
    }
    else{
      errz = 0.0014*(1+3*(chargematch-1.));
      erry = 0.003*(1+3*(chargematch-1.));
    } 
    if (cnz>nz+0.6) {
      erry*=(cnz-nz+0.5);
      errz*=1.4*(cnz-nz+0.5);
    }
  }

  if (TMath::Abs(cl->GetY())>2.5){
    factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
  }
  if (TMath::Abs(cl->GetY())<1){
    factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
  }
  factorz= TMath::Min(factorz,float(4.));  
  errz*=factorz;

  erry= TMath::Min(erry,float(0.05));
  errz= TMath::Min(errz,float(0.05));  
  return 200;
}
//--------------------------------------------------------------------------
void AliITSClusterParam::Print(Option_t* /*option*/) const {
  //
  // Print param Information
  //

  //
  // Error parameterization
  //
  printf("NOT YET...\n");
  return;
}
Commit	Line	Data
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>
	30	#include "TMath.h"
	31	#include "AliITSRecPoint.h"
	32	#include "AliITSClusterParam.h"
	33
	34	ClassImp(AliITSClusterParam)
	35
	36
	37	AliITSClusterParam* AliITSClusterParam::fgInstance = 0;
	38
	39
	40	/*
	41	Example usage fitting parameterization (NOT YET...):
	42	TFile fres("resol.root"); //tree with resolution and shape
	43	TTree * treeRes =(TTree*)fres.Get("Resol");
	44
	45	AliITSClusterParam param;
	46	param.SetInstance(&param);
	47	param.FitResol(treeRes);
	48	param.FitRMS(treeRes);
	49	TFile fparam("ITSClusterParam.root","recreate");
	50	param.Write("Param");
	51	//
	52	//
	53	TFile fparam("ITSClusterParam.root");
	54	AliITSClusterParam param2 = (AliITSClusterParam ) fparam.Get("Param");
	55	param2->SetInstance(param2);
	56	param2->Test(treeRes);
	57
	58
	59	treeRes->Draw("(Resol-AliITSClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0")
	60	*/
	61
	62	//-------------------------------------------------------------------------
	63	AliITSClusterParam* AliITSClusterParam::Instance()
	64	{
65	//
66	// Singleton implementation
67	// Returns an instance of this class, it is created if neccessary
68	//
69	if (fgInstance == 0){
70	fgInstance = new AliITSClusterParam();
71	}
72	return fgInstance;
73	}
74	//-------------------------------------------------------------------------
75	void AliITSClusterParam::GetNTeor(Int_t layer,const AliITSRecPoint* /cl/,
76	Float_t theta,Float_t phi,
77	Float_t &ny,Float_t &nz)
78	{
79	//
80	//get "mean shape"
81	//
82	if (layer==0){
83	ny = 1.+TMath::Abs(phi)*3.2;
84	nz = 1.+TMath::Abs(theta)*0.34;
85	return;
86	}
87	if (layer==1){
88	ny = 1.+TMath::Abs(phi)*3.2;
89	nz = 1.+TMath::Abs(theta)*0.28;
90	return;
91	}
92
93	if (layer>3){
94	ny = 2.02+TMath::Abs(phi)*1.95;
95	nz = 2.02+TMath::Abs(phi)*2.35;
96	return;
97	}
98	ny = 6.6-2.7*TMath::Abs(phi);
99	nz = 2.8-3.11TMath::Abs(phi)+0.45TMath::Abs(theta);
100	}
101	//--------------------------------------------------------------------------
102	Int_t AliITSClusterParam::GetError(Int_t layer,const AliITSRecPoint*cl,
103	Float_t theta,Float_t phi,Float_t expQ,
104	Float_t &erry,Float_t &errz)
105	{
106	//calculate cluster position error
107	//
108	Float_t nz,ny;
109	GetNTeor(layer, cl,theta,phi,ny,nz);
110	erry = TMath::Sqrt(cl->GetSigmaY2());
111	errz = TMath::Sqrt(cl->GetSigmaZ2());
112	//
113	// PIXELS
114	if (layer<2){
115	if (TMath::Abs(ny-cl->GetNy())>0.6) {
116	if (ny<cl->GetNy()){
117	erry*=0.4+TMath::Abs(ny-cl->GetNy());
118	errz*=0.4+TMath::Abs(ny-cl->GetNy());
119	}else{
120	erry=0.7+0.5TMath::Abs(ny-cl->GetNy());
121	errz=0.7+0.5TMath::Abs(ny-cl->GetNy());
122	}
123	}
124	if (TMath::Abs(nz-cl->GetNz())>1.) {
125	erry*=TMath::Abs(nz-cl->GetNz());
126	errz*=TMath::Abs(nz-cl->GetNz());
127	}
128	erry*=0.85;
129	errz*=0.85;
130	erry= TMath::Min(erry,float(0.0050));
131	errz= TMath::Min(errz,float(0.0300));
132	return 10;
133	}
134
135	//STRIPS
136	if (layer>3){
137	//factor 1.8 appears in new simulation
138	//
139	Float_t scale=1.8;
140	if (cl->GetNy()==100\|\|cl->GetNz()==100){
141	erry = 0.004*scale;
142	errz = 0.2*scale;
143	return 100;
144	}
145	if (cl->GetNy()+cl->GetNz()>12){
146	erry = 0.06*scale;
147	errz = 0.57*scale;
148	return 100;
149	}
150	Float_t normq = cl->GetQ()/(TMath::Sqrt(1+thetatheta+phiphi));
151	Float_t chargematch = TMath::Max(double(normq/expQ),2.);
152	//
153	if (cl->GetType()==1 \|\| cl->GetType()==10 ){
154	if (chargematch<1.0 \|\| (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
155	errz = 0.043*scale;
156	erry = 0.00094*scale;
157	return 101;
158	}
159	if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
160	errz = 0.06*scale;
161	erry =0.0013*scale;
162	return 102;
163	}
164	erry = 0.0027*scale;
165	errz = TMath::Min(0.028(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)scale;
166	return 103;
167	}
168	if (cl->GetType()==2 \|\| cl->GetType()==11 ){
169	erry = TMath::Min(0.0010(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)scale;
170	errz = TMath::Min(0.025(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)scale;
171	return 104;
172	}
173
174	if (cl->GetType()>100 ){
175	if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
176	errz = 0.05*scale;
177	erry = 0.00096*scale;
178	return 105;
179	}
180	if (cl->GetNy()+cl->GetNz()-nz-ny<1){
181	errz = 0.10*scale;
182	erry = 0.0025*scale;
183	return 106;
184	}
185
186	errz = TMath::Min(0.05(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)scale;
187	erry = TMath::Min(0.003(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)scale;
188	return 107;
189	}
190	Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
191	if (diff<1) diff=1;
192	if (diff>4) diff=4;
193
194	if (cl->GetType()==5\|\|cl->GetType()==6\|\|cl->GetType()==7\|\|cl->GetType()==8){
195	errz = 0.14*diff;
196	erry = 0.003*diff;
197	return 108;
198	}
199	erry = 0.04*diff;
200	errz = 0.06*diff;
201	return 109;
202	}
203	//DRIFTS
204	Float_t normq = cl->GetQ()/(TMath::Sqrt(1+thetatheta+phiphi));
205	Float_t chargematch = normq/expQ;
206	chargematch/=2.4; // F. Prino Sept. 2007: SDD charge conversion keV->ADC
207	Float_t factorz=1;
208	Int_t cnz = cl->GetNz()%10;
209	//charge match
210	if (cl->GetType()==1){
211	if (chargematch<1.25){
212	erry = 0.0028*(1.+6./cl->GetQ()); // gold clusters
213	}
214	else{
215	erry = 0.003*chargematch;
216	if (cl->GetNz()==3) erry*=1.5;
217	}
218	if (chargematch<1.0){
219	errz = 0.0011*(1.+6./cl->GetQ());
220	}
221	else{
222	errz = 0.002(1+2(chargematch-1.));
223	}
224	if (cnz>nz+0.6) {
225	erry*=(cnz-nz+0.5);
226	errz=1.4(cnz-nz+0.5);
227	}
228	}
229	if (cl->GetType()>1){
230	if (chargematch<1){
231	erry = 0.00385*(1.+6./cl->GetQ()); // gold clusters
232	errz = 0.0016*(1.+6./cl->GetQ());
233	}
234	else{
235	errz = 0.0014(1+3(chargematch-1.));
236	erry = 0.003(1+3(chargematch-1.));
237	}
238	if (cnz>nz+0.6) {
239	erry*=(cnz-nz+0.5);
240	errz=1.4(cnz-nz+0.5);
241	}
242	}
243
244	if (TMath::Abs(cl->GetY())>2.5){
245	factorz=1+2(TMath::Abs(cl->GetY())-2.5);
246	}
247	if (TMath::Abs(cl->GetY())<1){
248	factorz=1.+0.5TMath::Abs(TMath::Abs(cl->GetY())-1.);
249	}
250	factorz= TMath::Min(factorz,float(4.));
251	errz*=factorz;
252
253	erry= TMath::Min(erry,float(0.05));
254	errz= TMath::Min(errz,float(0.05));
255	return 200;
256	}
257	//--------------------------------------------------------------------------
258	void AliITSClusterParam::Print(Option_t* /option/) const {
259	//
260	// Print param Information
261	//
262
263	//
264	// Error parameterization
265	//
266	printf("NOT YET...\n");
267	return;
268	}
269
270
271
272
273