///////////////////////////////////////////////////////////////////////////////
// //
-// TPC cluster error and shape parameterization //
-// //
-// //
+// TPC cluster error, shape and charge parameterization as function
+// of drift length, and inclination angle //
+//
+// Following notation is used in following
+// Int_t dim 0 - y direction
+// 1 - z direction
+//
+// Int_t type 0 - short pads
+// 1 - medium pads
+// 2 - long pads
+// Float_t z - drift length
+//
+// Float_t angle - tangent of inclination angle at given dimension
+//
+// Implemented parameterization
+//
+//
+// 1. Resolution as function of drift length and inclination angle
+// 1.a) GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle)
+// Simple error parameterization as derived from analytical formula
+// only linear term in drift length and angle^2
+// The formula is valid only with precission +-5%
+// Separate parameterization for differnt pad geometry
+// 1.b) GetError0Par
+// Parabolic term correction - better precision
+//
+// 1.c) GetError1 - JUST FOR Study
+// Similar to GetError1
+// The angular and diffusion effect is scaling with pad length
+// common parameterization for different pad length
+//
+// 2. Error parameterization using charge
+// 2.a) GetErrorQ
+// GetError0+
+// adding 1/Q component to diffusion and angluar part
+// 2.b) GetErrorQPar
+// GetError0Par+
+// adding 1/Q component to diffusion and angluar part
+// 2.c) GetErrorQParScaled - Just for study
+// One parameterization for all pad shapes
+// Smaller precission as previous one
+//
+//
+// Example how to retrieve the paramterization:
+/*
+ AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ AliCDBManager::Instance()->SetRun(0)
+ AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam();
+
+ //
+ //
+ AliTPCClusterParam::SetInstance(param);
+ TF1 f1("f1","AliTPCClusterParam::SGetError0Par(1,0,x,0)",0,250);
+*/
+
+// EXAMPLE hot to create parameterization
+/*
+// Note resol is the resolution tree created by AliTPCcalibTracks
+//
+AliTPCClusterParam *param = new AliTPCClusterParam;
+param->FitData(Resol);
+AliTPCClusterParam::SetInstance(param);
+
+*/
+
+//
+// //
///////////////////////////////////////////////////////////////////////////////
#include "AliTPCClusterParam.h"
#include "TMath.h"
#include <TVectorF.h>
#include <TLinearFitter.h>
#include <TH1F.h>
+#include <TH3F.h>
#include <TProfile2D.h>
#include <TVectorD.h>
#include <TObjArray.h>
+#include "AliTPCcalibDB.h"
+#include "AliTPCParam.h"
ClassImp(AliTPCClusterParam)
}
+AliTPCClusterParam::AliTPCClusterParam():
+ TObject(),
+ fRatio(0),
+ fQNorm(0),
+ fQNormCorr(0),
+ fQNormHis(0),
+ fQpadTnorm(0), // q pad normalization - Total charge
+ fQpadMnorm(0) // q pad normalization - Max charge
+ //
+{
+ //
+ // Default constructor
+ //
+ fPosQTnorm[0] = 0; fPosQTnorm[1] = 0; fPosQTnorm[2] = 0;
+ fPosQMnorm[0] = 0; fPosQMnorm[1] = 0; fPosQMnorm[2] = 0;
+ //
+ fPosYcor[0] = 0; fPosYcor[1] = 0; fPosYcor[2] = 0;
+ fPosZcor[0] = 0; fPosZcor[1] = 0; fPosZcor[2] = 0;
+}
+
+AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param):
+ TObject(param),
+ fRatio(0),
+ fQNorm(0),
+ fQNormCorr(0),
+ fQNormHis(0),
+ fQpadTnorm(new TVectorD(*(param.fQpadTnorm))), // q pad normalization - Total charge
+ fQpadMnorm(new TVectorD(*(param.fQpadMnorm))) // q pad normalization - Max charge
+
+{
+ //
+ // copy constructor
+ //
+ memcpy(this, ¶m,sizeof(AliTPCClusterParam));
+ if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
+ if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
+ //
+ if (param.fPosQTnorm[0]){
+ fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
+ fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
+ fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
+ //
+ fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
+ fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
+ fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
+ }
+ if (param.fPosYcor[0]){
+ fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
+ fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
+ fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
+ //
+ fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
+ fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
+ fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
+ }
+
+}
+
+
+AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){
+ //
+ // Assignment operator
+ //
+ if (this != ¶m) {
+ memcpy(this, ¶m,sizeof(AliTPCClusterParam));
+ if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone();
+ if (param.fQNormHis) fQNormHis = (TObjArray*) param.fQNormHis->Clone();
+ if (param.fPosQTnorm[0]){
+ fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0]));
+ fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1]));
+ fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2]));
+ //
+ fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0]));
+ fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1]));
+ fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2]));
+ }
+ if (param.fPosYcor[0]){
+ fPosYcor[0] = new TVectorD(*(param.fPosYcor[0]));
+ fPosYcor[1] = new TVectorD(*(param.fPosYcor[1]));
+ fPosYcor[2] = new TVectorD(*(param.fPosYcor[2]));
+ //
+ fPosZcor[0] = new TVectorD(*(param.fPosZcor[0]));
+ fPosZcor[1] = new TVectorD(*(param.fPosZcor[1]));
+ fPosZcor[2] = new TVectorD(*(param.fPosZcor[2]));
+ }
+ }
+ return *this;
+}
+
+
+AliTPCClusterParam::~AliTPCClusterParam(){
+ //
+ // destructor
+ //
+ if (fQNorm) fQNorm->Delete();
+ if (fQNormCorr) delete fQNormCorr;
+ if (fQNormHis) fQNormHis->Delete();
+ delete fQNorm;
+ delete fQNormHis;
+ if (fPosQTnorm[0]){
+ delete fPosQTnorm[0];
+ delete fPosQTnorm[1];
+ delete fPosQTnorm[2];
+ //
+ delete fPosQMnorm[0];
+ delete fPosQMnorm[1];
+ delete fPosQMnorm[2];
+ }
+ if (fPosYcor[0]){
+ delete fPosYcor[0];
+ delete fPosYcor[1];
+ delete fPosYcor[2];
+ //
+ delete fPosZcor[0];
+ delete fPosZcor[1];
+ delete fPosZcor[2];
+ }
+}
void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){
Float_t err = fRatio*px[ipoint];
Double_t x[4];
x[0] = px[ipoint];
- fitter.AddPoint(x,val,err);
+ if (err>0) fitter.AddPoint(x,val,err);
}
fitter.Eval();
param0[0]= fitter.GetParameter(0);
value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ;
}
}
- return TMath::Sqrt(value);
+ return TMath::Sqrt(TMath::Abs(value));
}
//
printf("\nResolution Scaled factors\n");
printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n");
- printf("Y\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[0][0])),TMath::Sqrt(fParamS1[0][1]),
- TMath::Sqrt(fParamS1[0][2]),TMath::Sqrt(fParamS1[0][3]));
+ printf("Y\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[0][0])),TMath::Sqrt(TMath::Abs(fParamS1[0][1])),
+ TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3])));
for (Int_t ipad=0; ipad<3; ipad++){
Float_t length=0.75;
if (ipad==1) length=1;
if (ipad==2) length=1.5;
printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])),
- TMath::Sqrt(fParamS0[0][ipad][1]*length),
- TMath::Sqrt(fParamS0[0][ipad][2]/length),
- TMath::Sqrt(fParamS0[0][ipad][3]));
+ TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)),
+ TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)),
+ TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3])));
}
for (Int_t ipad=0; ipad<3; ipad++){
Float_t length=0.75;
if (ipad==2) length=1.5;
printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])),
- TMath::Sqrt(fParamS0Par[0][ipad][1]*length),
- TMath::Sqrt(fParamS0Par[0][ipad][2]/length),
- TMath::Sqrt(fParamS0Par[0][ipad][6]));
+ TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)),
+ TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)),
+ TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6])));
}
printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]),
TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3]));
if (ipad==2) length=1.5;
printf("\t%d\t%f\t%f\t%f\t%f\n", ipad,
TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])),
- TMath::Sqrt(fParamS0[1][ipad][1]*length),
- TMath::Sqrt(fParamS0[1][ipad][2]/length),
- TMath::Sqrt(fParamS0[1][ipad][3]));
+ TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)),
+ TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)),
+ TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3])));
}
for (Int_t ipad=0; ipad<3; ipad++){
Float_t length=0.75;
if (ipad==1) length=1;
if (ipad==2) length=1.5;
printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad,
- TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][0])),
- TMath::Sqrt(fParamS0Par[1][ipad][1]*length),
- TMath::Sqrt(fParamS0Par[1][ipad][2]/length),
- TMath::Sqrt(fParamS0Par[1][ipad][6]));
+ TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))),
+ TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)),
+ TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)),
+ TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6])));
}
//
printf("\n");
printf("\nRMS Scaled factors\n");
printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n");
- printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),TMath::Sqrt(fParamRMS1[0][1]),
- TMath::Sqrt(fParamRMS1[0][2]),TMath::Sqrt(fParamRMS1[0][3]),TMath::Sqrt(fParamRMS1[0][4]));
+ printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n",
+ TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[0][4])));
for (Int_t ipad=0; ipad<3; ipad++){
Float_t length=0.75;
if (ipad==1) length=1;
printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
0.,
- TMath::Sqrt(fParamRMS0[0][ipad][1]),
- TMath::Sqrt(fParamRMS0[0][ipad][2]/(length*length)),
- TMath::Sqrt(fParamRMS0[0][ipad][3]));
+ TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
}else{
printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
0.,
TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])),
- TMath::Sqrt(fParamRMS0[0][ipad][1]),
- TMath::Sqrt(fParamRMS0[0][ipad][2]/(length*length)),
- TMath::Sqrt(fParamRMS0[0][ipad][3]));
+ TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3])));
}
}
printf("\n");
- printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),TMath::Sqrt(fParamRMS1[1][1]),
- TMath::Sqrt(fParamRMS1[1][2]),TMath::Sqrt(fParamRMS1[1][3]),TMath::Sqrt(fParamRMS1[1][4]));
+ printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n",
+ TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])),
+ TMath::Sqrt(TMath::Abs(fParamRMS1[1][4])));
for (Int_t ipad=0; ipad<3; ipad++){
Float_t length=0.75;
if (ipad==1) length=1;
printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
0.,
- TMath::Sqrt(fParamRMS0[1][ipad][1]),
- TMath::Sqrt(fParamRMS0[1][ipad][2]/(length*length)),
- TMath::Sqrt(fParamRMS0[1][ipad][3]));
+ TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
}else{
printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad,
0.,
TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])),
- TMath::Sqrt(fParamRMS0[1][ipad][1]),
- TMath::Sqrt(fParamRMS0[1][ipad][2]/(length*length)),
- TMath::Sqrt(fParamRMS0[1][ipad][3]));
+ TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))),
+ TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3])));
}
}
}
// type - 0 Qtot 1 Qmax
// ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
//
- //
- //formula= dr++tz++ty++dr*tz++dr*ty++ty*tz++dr**2++ty**2++tz**2
- if (!fQNorm) return 0;
+ //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++++dr*dr++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000);
+
+ if (fQNorm==0) return 0;
TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad);
if (!norm) return 0;
- TVectorD &no = *norm;
- Float_t res= no[0]+no[1]*dr+no[2]*tz+no[3]*ty+no[4]*dr*tz+no[5]*dr*ty+no[6]*ty*tz
- +no[7]*dr*dr+no[8]*ty*ty+no[9]*tz*tz;
+ TVectorD &no = *norm;
+ Float_t res =
+ no[0]+
+ no[1]*dr+
+ no[2]*ty+
+ no[3]*tz+
+ no[4]*dr*ty+
+ no[5]*dr*tz+
+ no[6]*ty*tz+
+ no[7]*dr*dr+
+ no[8]*ty*ty+
+ no[9]*tz*tz;
+ res/=no[0];
+ return res;
+}
+
+
+
+Float_t AliTPCClusterParam::QnormHis(Int_t ipad, Int_t itype, Float_t dr, Float_t p2, Float_t p3){
+ // get Q normalization
+ // type - 0 Qtot 1 Qmax
+ // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm)
+ //
+
+ if (fQNormHis==0) return 0;
+ TH3F * norm = (TH3F*)fQNormHis->At(4*itype+ipad);
+ if (!norm) return 1;
+ p2=TMath::Abs(p2);
+ dr=TMath::Min(dr,Float_t(norm->GetXaxis()->GetXmax()-norm->GetXaxis()->GetBinWidth(0)));
+ dr=TMath::Max(dr,Float_t(norm->GetXaxis()->GetXmin()+norm->GetXaxis()->GetBinWidth(0)));
+ //
+ p2=TMath::Min(p2,Float_t(norm->GetYaxis()->GetXmax()-norm->GetYaxis()->GetBinWidth(0)));
+ p2=TMath::Max(p2,Float_t(norm->GetYaxis()->GetXmin()+norm->GetYaxis()->GetBinWidth(0)));
+ //
+ p3=TMath::Min(p3,Float_t(norm->GetZaxis()->GetXmax()-norm->GetZaxis()->GetBinWidth(0)));
+ p3=TMath::Max(p3,Float_t(norm->GetZaxis()->GetXmin()+norm->GetZaxis()->GetBinWidth(0)));
+ //
+ Double_t res = norm->GetBinContent(norm->FindBin(dr,p2,p3));
+ if (res==0) res = norm->GetBinContent(norm->FindBin(0.5,0.5,0.5)); // This is just hack - to be fixed entries without
return res;
}
if (fQNorm==0) fQNorm = new TObjArray(6);
fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad);
}
+
+void AliTPCClusterParam::ResetQnormCorr(){
+ //
+ //
+ //
+ if (!fQNormCorr) fQNormCorr= new TMatrixD(12,6);
+ for (Int_t irow=0;irow<12; irow++)
+ for (Int_t icol=0;icol<6; icol++){
+ (*fQNormCorr)(irow,icol)=1.; // default - no correction
+ if (irow>5) (*fQNormCorr)(irow,icol)=0.; // default - no correction
+ }
+}
+
+void AliTPCClusterParam::SetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType, Float_t val){
+ //
+ // ipad - pad type
+ // itype - 0- qtot 1-qmax
+ // corrType - 0 - s0y corr - eff. PRF corr
+ // - 1 - s0z corr - eff. TRF corr
+ // - 2 - d0y - eff. diffusion correction y
+ // - 3 - d0z - eff. diffusion correction
+ // - 4 - eff length - eff.length - wire pitch + x diffsion
+ // - 5 - pad type normalization
+ if (!fQNormCorr) {
+ ResetQnormCorr();
+ }
+ //
+ // eff shap parameterization matrix
+ //
+ // rows
+ // itype*3+ipad - itype=0 qtot itype=1 qmax ipad=0
+ //
+ if (itype<2) (*fQNormCorr)(itype*3+ipad, corrType) *= val; // multiplicative correction
+ if (itype>=2) (*fQNormCorr)(itype*3+ipad, corrType)+= val; // additive correction
+}
+
+Double_t AliTPCClusterParam::GetQnormCorr(Int_t ipad, Int_t itype, Int_t corrType) const{
+ //
+ // see AliTPCClusterParam::SetQnormCorr
+ //
+ if (!fQNormCorr) return 0;
+ return (*fQNormCorr)(itype*3+ipad, corrType);
+}
+
+
+Float_t AliTPCClusterParam::QnormPos(Int_t ipad,Bool_t isMax, Float_t pad, Float_t time, Float_t z, Float_t sy2, Float_t sz2, Float_t qm, Float_t qt){
+ //
+ // Make Q normalization as function of following parameters
+ // Fit parameters to be used in corresponding correction function extracted in the AliTPCclaibTracksGain - Taylor expansion
+ // 1 - dp - relative pad position
+ // 2 - dt - relative time position
+ // 3 - di - drift length (norm to 1);
+ // 4 - dq0 - Tot/Max charge
+ // 5 - dq1 - Max/Tot charge
+ // 6 - sy - sigma y - shape
+ // 7 - sz - sigma z - shape
+ //
+
+ //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain -
+ // Following variable used - correspondance to the our variable conventions
+ //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)");
+ Double_t dp = ((pad-int(pad)-0.5)*2.);
+ //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)");
+ Double_t dt = ((time-int(time)-0.5)*2.);
+ //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))");
+ Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.);
+ //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))");
+ Double_t dq0 = 0.2*(qt+2.)/(qm+2.);
+ //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))");
+ Double_t dq1 = 5.*(qm+2.)/(qt+2.);
+ //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))");
+ Double_t sy = 0.32/TMath::Sqrt(0.01*0.01+sy2);
+ //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))");
+ Double_t sz = 0.32/TMath::Sqrt(0.01*0.01+sz2);
+ //
+ //
+ //
+ TVectorD * pvec = 0;
+ if (isMax){
+ pvec = fPosQMnorm[ipad];
+ }else{
+ pvec = fPosQTnorm[ipad];
+ }
+ TVectorD ¶m = *pvec;
+ //
+ // Eval part - in correspondance with fit part from debug streamer
+ //
+ Double_t result=param[0];
+ Int_t index =1;
+ //
+ result+=dp*param[index++]; //1
+ result+=dt*param[index++]; //2
+ result+=dp*dp*param[index++]; //3
+ result+=dt*dt*param[index++]; //4
+ result+=dt*dt*dt*param[index++]; //5
+ result+=dp*dt*param[index++]; //6
+ result+=dp*dt*dt*param[index++]; //7
+ result+=(dq0)*param[index++]; //8
+ result+=(dq1)*param[index++]; //9
+ //
+ //
+ result+=dp*dp*(di)*param[index++]; //10
+ result+=dt*dt*(di)*param[index++]; //11
+ result+=dp*dp*sy*param[index++]; //12
+ result+=dt*sz*param[index++]; //13
+ result+=dt*dt*sz*param[index++]; //14
+ result+=dt*dt*dt*sz*param[index++]; //15
+ //
+ result+=dp*dp*1*sy*sz*param[index++]; //16
+ result+=dt*sy*sz*param[index++]; //17
+ result+=dt*dt*sy*sz*param[index++]; //18
+ result+=dt*dt*dt*sy*sz*param[index++]; //19
+ //
+ result+=dp*dp*(dq0)*param[index++]; //20
+ result+=dt*1*(dq0)*param[index++]; //21
+ result+=dt*dt*(dq0)*param[index++]; //22
+ result+=dt*dt*dt*(dq0)*param[index++]; //23
+ //
+ result+=dp*dp*(dq1)*param[index++]; //24
+ result+=dt*(dq1)*param[index++]; //25
+ result+=dt*dt*(dq1)*param[index++]; //26
+ result+=dt*dt*dt*(dq1)*param[index++]; //27
+
+ if (result<0.75) result=0.75;
+ if (result>1.25) result=1.25;
+
+ return result;
+
+}
+
+
+
+
+
+Float_t AliTPCClusterParam::PosCorrection(Int_t type, Int_t ipad, Float_t pad, Float_t time, Float_t z, Float_t /*sy2*/, Float_t /*sz2*/, Float_t /*qm*/){
+
+ //
+ // Make postion correction
+ // type - 0 - y correction
+ // 1 - z correction
+ // ipad - 0, 1, 2 - short, medium long pads
+ // pad - float pad number
+ // time - float time bin number
+ // z - z of the cluster
+
+ //
+ //chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)");
+ //chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)");
+ //chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)");
+ //chainres->SetAlias("st","(sin(dt)-dt)");
+ //
+ //chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))");
+
+ //
+ // Derived variables
+ //
+ Double_t dp = (-1+(z>0)*2)*((pad-int(pad))-0.5);
+ Double_t dt = (-1+(z>0)*2)*((time-0.66-int(time-0.66))-0.5);
+ Double_t sp = (TMath::Sin(dp*TMath::Pi())-dp*TMath::Pi());
+ Double_t st = (TMath::Sin(dt)-dt);
+ //
+ Double_t di = TMath::Sqrt(TMath::Abs(1.-TMath::Abs(z/250.)));
+ //
+ //
+ //
+ TVectorD * pvec = 0;
+ if (type==0){
+ pvec = fPosYcor[ipad];
+ }else{
+ pvec = fPosZcor[ipad];
+ }
+ TVectorD ¶m = *pvec;
+ //
+ Double_t result=0;
+ Int_t index =1;
+
+ if (type==0){
+ // y corr
+ result+=(dp)*param[index++]; //1
+ result+=(dp)*di*param[index++]; //2
+ //
+ result+=(sp)*param[index++]; //3
+ result+=(sp)*di*param[index++]; //4
+ }
+ if (type==1){
+ result+=(dt)*param[index++]; //1
+ result+=(dt)*di*param[index++]; //2
+ //
+ result+=(st)*param[index++]; //3
+ result+=(st)*di*param[index++]; //4
+ }
+ if (TMath::Abs(result)>0.05) return 0;
+ return result;
+}
+
+
+
+Double_t AliTPCClusterParam::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
+ //
+ // 2 D gaus convoluted with angular effect
+ // See in mathematica:
+ //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-k1*xd)*(x1-k1*xd)/(2*s1*s1)]/(s0*s1),{xd,-1/2,1/2}]]
+ //
+ //TF1 f1("f1","AliTPCClusterParam::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
+ //TF2 f2("f2","AliTPCClusterParam::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
+ //
+ const Float_t kEpsilon = 0.0001;
+ if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
+ // small angular effect
+ Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
+ return val;
+ }
+ Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
+ Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
+ //
+ Double_t sigmaErf = 2*s0*s1*TMath::Sqrt(2*sigma2);
+ Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
+ Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
+ Double_t norm = 1./TMath::Sqrt(sigma2);
+ norm/=2.*TMath::Sqrt(2.*TMath::Pi());
+ Double_t val = norm*exp0*(erf0+erf1);
+ return val;
+}
+
+
+Double_t AliTPCClusterParam::GaussConvolutionTail(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
+ //
+ // 2 D gaus convoluted with angular effect and exponential tail in z-direction
+ // tail integrated numerically
+ // Integral normalized to one
+ // Mean at 0
+ //
+ // TF1 f1t("f1t","AliTPCClusterParam::GaussConvolutionTail(0,x,0,0,0.5,0.5,0.9)",-5,5)
+ Double_t sum =1, mean=0;
+ // the COG of exponent
+ for (Float_t iexp=0;iexp<5;iexp+=0.2){
+ mean+=iexp*TMath::Exp(-iexp/tau);
+ sum +=TMath::Exp(-iexp/tau);
+ }
+ mean/=sum;
+ //
+ sum = 1;
+ Double_t val = GaussConvolution(x0,x1+mean, k0, k1 , s0,s1);
+ for (Float_t iexp=0;iexp<5;iexp+=0.2){
+ val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*TMath::Exp(-iexp/tau);
+ sum+=TMath::Exp(-iexp/tau);
+ }
+ return val/sum;
+}
+
+Double_t AliTPCClusterParam::GaussConvolutionGamma4(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1, Double_t tau){
+ //
+ // 2 D gaus convoluted with angular effect and exponential tail in z-direction
+ // tail integrated numerically
+ // Integral normalized to one
+ // Mean at 0
+ //
+ // TF1 f1g4("f1g4","AliTPCClusterParam::GaussConvolutionGamma4(0,x,0,0,0.5,0.2,1.6)",-5,5)
+ // TF2 f2g4("f2g4","AliTPCClusterParam::GaussConvolutionGamma4(y,x,0,0,0.5,0.2,1.6)",-5,5,-5,5)
+ Double_t sum =0, mean=0;
+ // the COG of G4
+ for (Float_t iexp=0;iexp<5;iexp+=0.2){
+ Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
+ mean+=iexp*g4;
+ sum +=g4;
+ }
+ mean/=sum;
+ //
+ sum = 0;
+ Double_t val = 0;
+ for (Float_t iexp=0;iexp<5;iexp+=0.2){
+ Double_t g4 = TMath::Exp(-4.*iexp/tau)*TMath::Power(iexp/tau,4.);
+ val+=GaussConvolution(x0,x1+mean-iexp, k0, k1 , s0,s1)*g4;
+ sum+=g4;
+ }
+ return val/sum;
+}
+
+Double_t AliTPCClusterParam::QmaxCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0, Float_t effPad, Float_t effDiff){
+ //
+ //
+ // cpad - pad (y) coordinate
+ // ctime - time(z) coordinate
+ // ky - dy/dx
+ // kz - dz/dx
+ // rmsy0 - RF width in pad units
+ // rmsz0 - RF width in time bin units
+ // effLength - contibution of PRF and diffusion
+ // effDiff - overwrite diffusion
+
+ // Response function aproximated by convolution of gaussian with angular effect (integral=1)
+ //
+ // Gaus width sy and sz is determined by RF width and diffusion
+ // Integral of Q is equal 1
+ // Q max is calculated at position cpad, ctime
+ // Example function:
+ // TF1 f1("f1", "AliTPCClusterParam::QmaxCorrection(0,0.5,x,0,0,0.5,0.6)",0,1000)
+ //
+ AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
+ Double_t padLength= param->GetPadPitchLength(sector,row);
+ Double_t padWidth = param->GetPadPitchWidth(sector);
+ Double_t zwidth = param->GetZWidth();
+ Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
+
+ // diffusion in pad, time bin units
+ Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
+ Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
+ diffT*=effDiff; //
+ diffL*=effDiff; //
+ //
+ // transform angular effect to pad units
+ //
+ Double_t pky = ky*effLength/padWidth;
+ Double_t pkz = kz*effLength/zwidth;
+ // position in pad unit
+ Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
+ Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
+ //
+ //
+ Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
+ Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
+ //return GaussConvolutionGamma4(py,pz, pky,pkz,sy,sz,tau);
+ Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
+ return GaussConvolution(py,pz, pky,pkz,sy,sz)*length;
+}
+
+Double_t AliTPCClusterParam::QtotCorrection(Int_t sector, Int_t row, Float_t cpad, Float_t ctime, Float_t ky, Float_t kz, Float_t rmsy0, Float_t rmsz0, Float_t qtot, Float_t thr, Float_t effPad, Float_t effDiff){
+ //
+ //
+ // cpad - pad (y) coordinate
+ // ctime - time(z) coordinate
+ // ky - dy/dx
+ // kz - dz/dx
+ // rmsy0 - RF width in pad units
+ // rmsz0 - RF width in time bin units
+ // qtot - the sum of signal in cluster - without thr correction
+ // thr - threshold
+ // effLength - contibution of PRF and diffusion
+ // effDiff - overwrite diffusion
+
+ // Response function aproximated by convolution of gaussian with angular effect (integral=1)
+ //
+ // Gaus width sy and sz is determined by RF width and diffusion
+ // Integral of Q is equal 1
+ // Q max is calculated at position cpad, ctime
+ //
+ //
+ //
+ AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
+ Double_t padLength= param->GetPadPitchLength(sector,row);
+ Double_t padWidth = param->GetPadPitchWidth(sector);
+ Double_t zwidth = param->GetZWidth();
+ Double_t effLength= padLength+(param->GetWWPitch(0)+TMath::Sqrt(ctime*zwidth)*param->GetDiffT())*effPad;
+ //
+ // diffusion in pad units
+ Double_t diffT=TMath::Sqrt(ctime*zwidth)*param->GetDiffT()/padWidth;
+ Double_t diffL=TMath::Sqrt(ctime*zwidth)*param->GetDiffL()/zwidth;
+ diffT*=effDiff; //
+ diffL*=effDiff; //
+ //
+ // transform angular effect to pad units
+ Double_t pky = ky*effLength/padWidth;
+ Double_t pkz = kz*effLength/zwidth;
+ // position in pad unit
+ //
+ Double_t py = (cpad+0.5)-TMath::Nint(cpad+0.5);
+ Double_t pz = (ctime+0.5)-TMath::Nint(ctime+0.5);
+ //
+ Double_t sy = TMath::Sqrt(rmsy0*rmsy0+diffT*diffT);
+ Double_t sz = TMath::Sqrt(rmsz0*rmsz0+diffL*diffL);
+ //
+ //
+ //
+ Double_t sumAll=0,sumThr=0;
+ //
+ Double_t corr =1;
+ Double_t qnorm=qtot;
+ for (Float_t iy=-3;iy<=3;iy+=1.)
+ for (Float_t iz=-4;iz<=4;iz+=1.){
+ // Double_t val = GaussConvolutionGamma4(py-iy,pz-iz, pky,pkz, sy,sz,tau);
+ Double_t val = GaussConvolution(py-iy,pz-iz, pky,pkz, sy,sz);
+ Double_t qlocal =qnorm*val;
+ if (TMath::Abs(iy)<1.5&&TMath::Abs(iz)<1.5){
+ sumThr+=qlocal; // Virtual charge used in cluster finder
+ }
+ else{
+ if (qlocal>thr && TMath::Abs(iz)<2.5&&TMath::Abs(iy)<2.5) sumThr+=qlocal;
+ }
+ sumAll+=qlocal;
+ }
+ if (sumAll>0&&sumThr>0) {
+ corr=(sumThr)/sumAll;
+ }
+ //
+ Double_t length = padLength*TMath::Sqrt(1+ky*ky+kz*kz);
+ return corr*length;
+}
+
+
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff