X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCClusterParam.cxx;h=71c8f76fb68c01e15791f0b7269a31089ca5e5c9;hb=aa62360d3b63ed897ceeaa3af8f1b75cc71c9e3f;hp=80f6d05e98de844bbf2f79803565996c7f32d106;hpb=236a0d03e20ddc80ba1354bdf49fd818d4f1599b;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCClusterParam.cxx b/TPC/AliTPCClusterParam.cxx index 80f6d05e98d..71c8f76fb68 100644 --- a/TPC/AliTPCClusterParam.cxx +++ b/TPC/AliTPCClusterParam.cxx @@ -61,7 +61,7 @@ // // Example how to retrieve the paramterization: /* - AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT"); + AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); AliCDBManager::Instance()->SetRun(0) AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam(); @@ -91,9 +91,14 @@ AliTPCClusterParam::SetInstance(param); #include #include #include +#include #include #include #include +#include "AliTPCcalibDB.h" +#include "AliTPCParam.h" + +#include "AliMathBase.h" ClassImp(AliTPCClusterParam) @@ -144,25 +149,62 @@ AliTPCClusterParam* AliTPCClusterParam::Instance() AliTPCClusterParam::AliTPCClusterParam(): TObject(), fRatio(0), - fQNorm(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) + 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 @@ -170,6 +212,25 @@ AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& par 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; } @@ -180,7 +241,28 @@ 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]; + } } @@ -189,24 +271,24 @@ void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t // Fit z - angular dependence of resolution // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"Resol:AngleM:Zm"); - char varErr[100]; - sprintf(varErr,"Sigma:AngleS:Zs"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); - // - Int_t entries = tree->Draw(varVal,varCut); + TString varVal; + varVal="Resol:AngleM:Zm"; + TString varErr; + varErr="Sigma:AngleS:Zs"; + TString varCut; + varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type); + // + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[10000], py[10000], pz[10000]; Float_t ex[10000], ey[10000], ez[10000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint]= tree->GetV2()[ipoint]; ez[ipoint]= tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -242,24 +324,24 @@ void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float // Fit z - angular dependence of resolution // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"Resol:AngleM:Zm"); - char varErr[100]; - sprintf(varErr,"Sigma:AngleS:Zs"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); - // - Int_t entries = tree->Draw(varVal,varCut); + TString varVal; + varVal="Resol:AngleM:Zm"; + TString varErr; + varErr="Sigma:AngleS:Zs"; + TString varCut; + varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type); + // + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[10000], py[10000], pz[10000]; Float_t ex[10000], ey[10000], ez[10000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint]= tree->GetV2()[ipoint]; ez[ipoint]= tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -307,24 +389,24 @@ void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Flo // Fit z - angular dependence of resolution - pad length scaling // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length"); - char varErr[100]; - sprintf(varErr,"Sigma:AngleS:Zs"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&QMean<0",dim); - // - Int_t entries = tree->Draw(varVal,varCut); + TString varVal; + varVal="Resol:AngleM*sqrt(Length):Zm/Length"; + TString varErr; + varErr="Sigma:AngleS:Zs"; + TString varCut; + varCut=Form("Dim==%d&&QMean<0",dim); + // + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[10000], py[10000], pz[10000]; Float_t ex[10000], ey[10000], ez[10000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint]= tree->GetV2()[ipoint]; ez[ipoint]= tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -359,27 +441,27 @@ void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t // Fit z - angular dependence of resolution - Q scaling // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean"); + TString varVal; + varVal="Resol:AngleM/sqrt(QMean):Zm/QMean"; char varVal0[100]; sprintf(varVal0,"Resol:AngleM:Zm"); // - char varErr[100]; - sprintf(varErr,"Sigma:AngleS:Zs"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type); + TString varErr; + varErr="Sigma:AngleS:Zs"; + TString varCut; + varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type); // - Int_t entries = tree->Draw(varVal,varCut); + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000]; Float_t ex[20000], ey[20000], ez[20000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint]= tree->GetV2()[ipoint]; ez[ipoint]= tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -426,27 +508,27 @@ void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float // Fit z - angular dependence of resolution - Q scaling - parabolic correction // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean"); + TString varVal; + varVal="Resol:AngleM/sqrt(QMean):Zm/QMean"; char varVal0[100]; sprintf(varVal0,"Resol:AngleM:Zm"); // - char varErr[100]; - sprintf(varErr,"Sigma:AngleS:Zs"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type); + TString varErr; + varErr="Sigma:AngleS:Zs"; + TString varCut; + varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type); // - Int_t entries = tree->Draw(varVal,varCut); + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000]; Float_t ex[20000], ey[20000], ez[20000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint]= tree->GetV2()[ipoint]; ez[ipoint]= tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -506,24 +588,24 @@ void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *p // Fit z - angular dependence of resolution // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"RMSm:AngleM:Zm"); - char varErr[100]; - sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); - // - Int_t entries = tree->Draw(varVal,varCut); + TString varVal; + varVal="RMSm:AngleM:Zm"; + TString varErr; + varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs"; + TString varCut; + varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type); + // + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[10000], py[10000], pz[10000]; Float_t ex[10000], ey[10000], ez[10000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint]= tree->GetV2()[ipoint]; ez[ipoint]= tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -558,24 +640,24 @@ void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float // Fit z - angular dependence of resolution - pad length scaling // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"RMSm:AngleM*Length:Zm"); - char varErr[100]; - sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&QMean<0",dim); - // - Int_t entries = tree->Draw(varVal,varCut); + TString varVal; + varVal="RMSm:AngleM*Length:Zm"; + TString varErr; + varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad"; + TString varCut; + varCut=Form("Dim==%d&&QMean<0",dim); + // + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[10000], py[10000], pz[10000]; Float_t type[10000], ey[10000], ez[10000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint] = tree->GetV2()[ipoint]; ez[ipoint] = tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -613,27 +695,27 @@ void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *p // Fit z - angular dependence of resolution - Q scaling // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean"); + TString varVal; + varVal="RMSm:AngleM/sqrt(QMean):Zm/QMean"; char varVal0[100]; sprintf(varVal0,"RMSm:AngleM:Zm"); // - char varErr[100]; - sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs"); - char varCut[100]; - sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type); + TString varErr; + varErr="sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs"; + TString varCut; + varCut=Form("Dim==%d&&Pad==%d&&QMean>0",dim,type); // - Int_t entries = tree->Draw(varVal,varCut); + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000]; Float_t ex[20000], ey[20000], ez[20000]; // - tree->Draw(varErr,varCut); + tree->Draw(varErr.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; ey[ipoint]= tree->GetV2()[ipoint]; ez[ipoint]= tree->GetV1()[ipoint]; } - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV3()[ipoint]; py[ipoint]= tree->GetV2()[ipoint]; @@ -681,16 +763,16 @@ void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_ // Fit z - angular dependence of resolution - Q scaling // // Int_t dim=0, type=0; - char varVal[100]; - sprintf(varVal,"RMSs:RMSm"); + TString varVal; + varVal="RMSs:RMSm"; // - char varCut[100]; - sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); + TString varCut; + varCut=Form("Dim==%d&&Pad==%d&&QMean<0",dim,type); // - Int_t entries = tree->Draw(varVal,varCut); + Int_t entries = tree->Draw(varVal.Data(),varCut); Float_t px[20000], py[20000]; // - tree->Draw(varVal,varCut); + tree->Draw(varVal.Data(),varCut); for (Int_t ipoint=0; ipointGetV2()[ipoint]; py[ipoint]= tree->GetV1()[ipoint]; @@ -710,7 +792,7 @@ void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_ -Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){ +Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle) const { // // // @@ -723,7 +805,7 @@ Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t } -Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){ +Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle) const { // // // @@ -740,7 +822,7 @@ Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Floa -Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){ +Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle) const { // // // @@ -755,7 +837,7 @@ Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t return value; } -Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ +Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const { // // // @@ -771,7 +853,7 @@ Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t } -Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ +Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const { // // // @@ -790,7 +872,7 @@ Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Floa } -Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ +Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const { // // // @@ -811,7 +893,7 @@ Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z } -Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){ +Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle) const { // // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly // @@ -823,7 +905,7 @@ Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t a return value; } -Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){ +Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle) const { // // calculate mean RMS of cluster - z,angle - pad length scalling // @@ -842,7 +924,7 @@ Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t a return value; } -Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ +Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const { // // calculate mean RMS of cluster - z,angle, Q dependence // @@ -856,7 +938,7 @@ Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t a return value; } -Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ +Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean) const { // // calculates RMS of signal shape fluctuation // @@ -866,7 +948,7 @@ Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float return value; } -Float_t AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){ +Float_t AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM) const { // // calculates vallue - sigma distortion contribution // @@ -1248,13 +1330,14 @@ Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t t // type - 0 Qtot 1 Qmax // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm) // - //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000); + //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]+ + Float_t res = + no[0]+ no[1]*dr+ no[2]*ty+ no[3]*tz+ @@ -1270,7 +1353,34 @@ Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t t -void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){ +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; +} + + + +void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, const TVectorD *const norm){ // // set normalization // @@ -1282,32 +1392,406 @@ void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){ 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))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; +} + + -/* - gSystem->Load("libSTAT.so") - TStatToolkit toolkit; - Double_t chi2; - TVectorD fitParam; - TMatrixD covMatrix; - Int_t npoints; - - TString fstring=""; - // - fstring+="(Zm/250.)++" - fstring+="(Zm/250.)^2++" - fstring+="(Zm/250.)^3++" - fstring+="AngleM++" - fstring+="(AngleM)^2++" - fstring+="(AngleM)^3++" - fstring+="(Zm/250.)*AngleM++" - fstring+="(Zm/250.)^2*AngleM++" - fstring+="(Zm/250.)*AngleM^2++" - - TString *strRMSY0 = toolkit.FitPlane(treeRes,"RMSm^2",fstring->Data(), "Dim==0&&Pad==0&&QMean<0", chi2,npoints,fitParam,covMatrix); - - treeRes->SetAlias("RMSY0",strRMSY0->Data()); -*/