X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCseed.cxx;h=a8d4909f0b6daecc86ce253da36c1a9222846f4a;hp=a1219154ec6b7d09b97f35724721bc00ff95a4ab;hb=2efc897ff2dff3be721a1d887b5ab5d2b3010ffe;hpb=4cc7ba77206b30206572605955cdfd126932032b diff --git a/TPC/AliTPCseed.cxx b/TPC/AliTPCseed.cxx index a1219154ec6..a8d4909f0b6 100644 --- a/TPC/AliTPCseed.cxx +++ b/TPC/AliTPCseed.cxx @@ -22,13 +22,19 @@ // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch //----------------------------------------------------------------- #include "TClonesArray.h" +#include "TGraphErrors.h" #include "AliTPCseed.h" #include "AliTPCReconstructor.h" #include "AliTPCClusterParam.h" #include "AliTPCCalPad.h" #include "AliTPCCalROC.h" #include "AliTPCcalibDB.h" - +#include "AliTPCParam.h" +#include "AliMathBase.h" +#include "AliTPCTransform.h" +#include "AliSplineFit.h" +#include "AliCDBManager.h" +#include "AliTPCcalibDButil.h" ClassImp(AliTPCseed) @@ -441,8 +447,8 @@ Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const } // Double_t y1=fP0, z1=fP1; - Double_t c1=GetSnp(), r1=sqrt(1.- c1*c1); - Double_t c2=c1 + GetC()*dx, r2=sqrt(1.- c2*c2); + Double_t c1=GetSnp(), r1=sqrt((1.-c1)*(1.+c1)); + Double_t c2=c1 + GetC()*dx, r2=sqrt((1.-c2)*(1.+c2)); y = GetY(); z = GetZ(); @@ -479,6 +485,15 @@ Double_t AliTPCseed::GetPredictedChi2(const AliCluster *c) const //----------------------------------------------------------------- Double_t p[2]={c->GetY(), c->GetZ()}; Double_t cov[3]={fErrorY2, 0., fErrorZ2}; + + Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX(); + if (TMath::Abs(dx)>0){ + Float_t ty = TMath::Tan(TMath::ASin(GetSnp())); + Float_t dy = dx*ty; + Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl(); + p[0] = c->GetY()-dy; + p[1] = c->GetZ()-dz; + } return AliExternalTrackParam::GetPredictedChi2(p,cov); } @@ -528,6 +543,15 @@ Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_t index) Int_t idx=GetClusterIndex(n); // save the current cluster index AliCluster cl(*c); cl.SetSigmaY2(fErrorY2); cl.SetSigmaZ2(fErrorZ2); + Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX(); + if (TMath::Abs(dx)>0){ + Float_t ty = TMath::Tan(TMath::ASin(GetSnp())); + Float_t dy = dx*ty; + Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl(); + cl.SetY(c->GetY()-dy); + cl.SetZ(c->GetZ()-dz); + } + if (!AliTPCtrack::Update(&cl,chisq,index)) return kFALSE; if (fCMeanSigmaY2p30<0){ @@ -552,252 +576,272 @@ Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_t index) //_____________________________________________________________________________ -Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) { +Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t /* onlyused */) { //----------------------------------------------------------------- // This funtion calculates dE/dX within the "low" and "up" cuts. //----------------------------------------------------------------- - - Float_t amp[200]; - Float_t angular[200]; - Float_t weight[200]; - Int_t index[200]; - //Int_t nc = 0; - Float_t meanlog = 100.; + // CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal) + AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); - Float_t mean[4] = {0,0,0,0}; - Float_t sigma[4] = {1000,1000,1000,1000}; - Int_t nc[4] = {0,0,0,0}; - Float_t norm[4] = {1000,1000,1000,1000}; + Int_t row0 = param->GetNRowLow(); + Int_t row1 = row0+param->GetNRowUp1(); + Int_t row2 = row1+param->GetNRowUp2(); + const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam(); + Int_t useTot = 0; + if (recoParam) useTot = (recoParam->GetUseTotCharge())? 0:1; + // + // // + fDEDX[0] = CookdEdxAnalytical(low,up,useTot ,i1 ,i2, 0); + fDEDX[1] = CookdEdxAnalytical(low,up,useTot ,0 ,row0, 0); + fDEDX[2] = CookdEdxAnalytical(low,up,useTot ,row0,row1, 0); + fDEDX[3] = CookdEdxAnalytical(low,up,useTot ,row1,row2, 0); // - fNShared =0; + fSDEDX[0] = CookdEdxAnalytical(low,up,useTot ,i1 ,i2, 1); + fSDEDX[1] = CookdEdxAnalytical(low,up,useTot ,0 ,row0, 1); + fSDEDX[2] = CookdEdxAnalytical(low,up,useTot ,row0,row1, 1); + fSDEDX[3] = CookdEdxAnalytical(low,up,useTot ,row1,row2, 1); + // + fNCDEDX[0] = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,i1 ,i2, 2)); + fNCDEDX[1] = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,0 ,row0, 2)); + fNCDEDX[2] = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,row0,row1, 2)); + fNCDEDX[3] = TMath::Nint(CookdEdxAnalytical(low,up,useTot ,row1,row2, 2)); - Float_t gainGG = 1; - if (AliTPCcalibDB::Instance()->GetParameters()){ - gainGG= 20000./AliTPCcalibDB::Instance()->GetParameters()->GetGasGain(); //relative gas gain - } + SetdEdx(fDEDX[0]); + return fDEDX[0]; +// return CookdEdxNorm(low,up,0,i1,i2,1,0,2); - for (Int_t of =0; of<4; of++){ - for (Int_t i=of+i1;iIsUsed(10))) continue; - if (cl->IsUsed(11)) { - fNShared++; - continue; - } - Int_t type = cl->GetType(); - //if (point->fIsShared){ - // fNShared++; - // continue; - //} - //if (pointm) - // if (pointm->fIsShared) continue; - //if (pointp) - // if (pointp->fIsShared) continue; - - if (type<0) continue; - //if (type>10) continue; - //if (point->GetErrY()==0) continue; - //if (point->GetErrZ()==0) continue; - - //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY(); - //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ(); - //if ((ddy*ddy+ddz*ddz)>10) continue; - - - // if (point->GetCPoint().GetMax()<5) continue; - if (cl->GetMax()<5) continue; - Float_t angley = point->GetAngleY(); - Float_t anglez = point->GetAngleZ(); - - Float_t rsigmay2 = point->GetSigmaY(); - Float_t rsigmaz2 = point->GetSigmaZ(); - /* - Float_t ns = 1.; - if (pointm){ - rsigmay += pointm->GetTPoint().GetSigmaY(); - rsigmaz += pointm->GetTPoint().GetSigmaZ(); - ns+=1.; - } - if (pointp){ - rsigmay += pointp->GetTPoint().GetSigmaY(); - rsigmaz += pointp->GetTPoint().GetSigmaZ(); - ns+=1.; - } - rsigmay/=ns; - rsigmaz/=ns; - */ - - Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2); - - Float_t ampc = 0; // normalization to the number of electrons - if (i>64){ - // ampc = 1.*point->GetCPoint().GetMax(); - ampc = 1.*cl->GetMax(); - //ampc = 1.*point->GetCPoint().GetQ(); - // AliTPCClusterPoint & p = point->GetCPoint(); - // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5); - // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; - //Float_t dz = - // TMath::Abs( Int_t(iz) - iz + 0.5); - //ampc *= 1.15*(1-0.3*dy); - //ampc *= 1.15*(1-0.3*dz); - // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ())); - //ampc *=zfactor; - } - else{ - //ampc = 1.0*point->GetCPoint().GetMax(); - ampc = 1.0*cl->GetMax(); - //ampc = 1.0*point->GetCPoint().GetQ(); - //AliTPCClusterPoint & p = point->GetCPoint(); - // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5); - //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; - //Float_t dz = - // TMath::Abs( Int_t(iz) - iz + 0.5); - - //ampc *= 1.15*(1-0.3*dy); - //ampc *= 1.15*(1-0.3*dz); - // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ())); - //ampc *=zfactor; - } - ampc *= 2.0; // put mean value to channel 50 - //ampc *= 0.58; // put mean value to channel 50 - Float_t w = 1.; - // if (type>0) w = 1./(type/2.-0.5); - // Float_t z = TMath::Abs(cl->GetZ()); - if (i<64) { - ampc /= 0.6; - //ampc /= (1+0.0008*z); - } else - if (i>128){ - ampc /=1.5; - //ampc /= (1+0.0008*z); - }else{ - //ampc /= (1+0.0008*z); - } +// Float_t amp[200]; +// Float_t angular[200]; +// Float_t weight[200]; +// Int_t index[200]; +// //Int_t nc = 0; +// Float_t meanlog = 100.; + +// Float_t mean[4] = {0,0,0,0}; +// Float_t sigma[4] = {1000,1000,1000,1000}; +// Int_t nc[4] = {0,0,0,0}; +// Float_t norm[4] = {1000,1000,1000,1000}; +// // +// // +// fNShared =0; + +// Float_t gainGG = 1; +// if (AliTPCcalibDB::Instance()->GetParameters()){ +// gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000.; //relative gas gain +// } + + +// for (Int_t of =0; of<4; of++){ +// for (Int_t i=of+i1;iIsUsed(10))) continue; +// if (cl->IsUsed(11)) { +// fNShared++; +// continue; +// } +// Int_t type = cl->GetType(); +// //if (point->fIsShared){ +// // fNShared++; +// // continue; +// //} +// //if (pointm) +// // if (pointm->fIsShared) continue; +// //if (pointp) +// // if (pointp->fIsShared) continue; + +// if (type<0) continue; +// //if (type>10) continue; +// //if (point->GetErrY()==0) continue; +// //if (point->GetErrZ()==0) continue; + +// //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY(); +// //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ(); +// //if ((ddy*ddy+ddz*ddz)>10) continue; + + +// // if (point->GetCPoint().GetMax()<5) continue; +// if (cl->GetMax()<5) continue; +// Float_t angley = point->GetAngleY(); +// Float_t anglez = point->GetAngleZ(); + +// Float_t rsigmay2 = point->GetSigmaY(); +// Float_t rsigmaz2 = point->GetSigmaZ(); +// /* +// Float_t ns = 1.; +// if (pointm){ +// rsigmay += pointm->GetTPoint().GetSigmaY(); +// rsigmaz += pointm->GetTPoint().GetSigmaZ(); +// ns+=1.; +// } +// if (pointp){ +// rsigmay += pointp->GetTPoint().GetSigmaY(); +// rsigmaz += pointp->GetTPoint().GetSigmaZ(); +// ns+=1.; +// } +// rsigmay/=ns; +// rsigmaz/=ns; +// */ + +// Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2); + +// Float_t ampc = 0; // normalization to the number of electrons +// if (i>64){ +// // ampc = 1.*point->GetCPoint().GetMax(); +// ampc = 1.*cl->GetMax(); +// //ampc = 1.*point->GetCPoint().GetQ(); +// // AliTPCClusterPoint & p = point->GetCPoint(); +// // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5); +// // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; +// //Float_t dz = +// // TMath::Abs( Int_t(iz) - iz + 0.5); +// //ampc *= 1.15*(1-0.3*dy); +// //ampc *= 1.15*(1-0.3*dz); +// // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ())); +// //ampc *=zfactor; +// } +// else{ +// //ampc = 1.0*point->GetCPoint().GetMax(); +// ampc = 1.0*cl->GetMax(); +// //ampc = 1.0*point->GetCPoint().GetQ(); +// //AliTPCClusterPoint & p = point->GetCPoint(); +// // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5); +// //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; +// //Float_t dz = +// // TMath::Abs( Int_t(iz) - iz + 0.5); + +// //ampc *= 1.15*(1-0.3*dy); +// //ampc *= 1.15*(1-0.3*dz); +// // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ())); +// //ampc *=zfactor; + +// } +// ampc *= 2.0; // put mean value to channel 50 +// //ampc *= 0.58; // put mean value to channel 50 +// Float_t w = 1.; +// // if (type>0) w = 1./(type/2.-0.5); +// // Float_t z = TMath::Abs(cl->GetZ()); +// if (i<64) { +// ampc /= 0.6; +// //ampc /= (1+0.0008*z); +// } else +// if (i>128){ +// ampc /=1.5; +// //ampc /= (1+0.0008*z); +// }else{ +// //ampc /= (1+0.0008*z); +// } - if (type<0) { //amp at the border - lower weight - // w*= 2.; +// if (type<0) { //amp at the border - lower weight +// // w*= 2.; - continue; - } - if (rsigma>1.5) ampc/=1.3; // if big backround - amp[nc[of]] /=gainGG; - amp[nc[of]] = ampc; - angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez); - weight[nc[of]] = w; - nc[of]++; - } +// continue; +// } +// if (rsigma>1.5) ampc/=1.3; // if big backround +// amp[nc[of]] = ampc; +// amp[nc[of]] /=gainGG; +// angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez); +// weight[nc[of]] = w; +// nc[of]++; +// } - TMath::Sort(nc[of],amp,index,kFALSE); - Float_t sumamp=0; - Float_t sumamp2=0; - Float_t sumw=0; - //meanlog = amp[index[Int_t(nc[of]*0.33)]]; - meanlog = 50; - for (Int_t i=int(nc[of]*low+0.5);i0.1) - sigma[of] = TMath::Sqrt(sigma[of]); - else - sigma[of] = 1000; +// TMath::Sort(nc[of],amp,index,kFALSE); +// Float_t sumamp=0; +// Float_t sumamp2=0; +// Float_t sumw=0; +// //meanlog = amp[index[Int_t(nc[of]*0.33)]]; +// meanlog = 50; +// for (Int_t i=int(nc[of]*low+0.5);i0.1) +// sigma[of] = TMath::Sqrt(sigma[of]); +// else +// sigma[of] = 1000; - mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; - //mean *=(1-0.02*(sigma/(mean*0.17)-1.)); - //mean *=(1-0.1*(norm-1.)); - } - } - - Float_t dedx =0; - fSdEdx =0; - fMAngular =0; - // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1)); - // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1)); +// mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; +// //mean *=(1-0.02*(sigma/(mean*0.17)-1.)); +// //mean *=(1-0.1*(norm-1.)); +// } +// } + +// Float_t dedx =0; +// fSdEdx =0; +// fMAngular =0; +// // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1)); +// // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1)); - // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ - // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1]))); - - Int_t norm2 = 0; - Int_t norm3 = 0; - for (Int_t i =0;i<4;i++){ - if (nc[i]>2&&nc[i]<1000){ - dedx += mean[i] *nc[i]; - fSdEdx += sigma[i]*(nc[i]-2); - fMAngular += norm[i] *nc[i]; - norm2 += nc[i]; - norm3 += nc[i]-2; - } - fDEDX[i] = mean[i]; - fSDEDX[i] = sigma[i]; - fNCDEDX[i]= nc[i]; - } - - if (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; - } - else{ - SetdEdx(0); - return 0; - } - // Float_t dedx1 =dedx; - /* - dedx =0; - for (Int_t i =0;i<4;i++){ - if (nc[i]>2&&nc[i]<1000){ - mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.)); - dedx += mean[i] *nc[i]; - } - fDEDX[i] = mean[i]; - } - dedx /= norm2; - */ +// // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ +// // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1]))); + +// Int_t norm2 = 0; +// Int_t norm3 = 0; +// for (Int_t i =0;i<4;i++){ +// if (nc[i]>2&&nc[i]<1000){ +// dedx += mean[i] *nc[i]; +// fSdEdx += sigma[i]*(nc[i]-2); +// fMAngular += norm[i] *nc[i]; +// norm2 += nc[i]; +// norm3 += nc[i]-2; +// } +// fDEDX[i] = mean[i]; +// fSDEDX[i] = sigma[i]; +// fNCDEDX[i]= nc[i]; +// } + +// if (norm3>0){ +// dedx /=norm2; +// fSdEdx /=norm3; +// fMAngular/=norm2; +// } +// else{ +// SetdEdx(0); +// return 0; +// } +// // Float_t dedx1 =dedx; +// /* +// dedx =0; +// for (Int_t i =0;i<4;i++){ +// if (nc[i]>2&&nc[i]<1000){ +// mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.)); +// dedx += mean[i] *nc[i]; +// } +// fDEDX[i] = mean[i]; +// } +// dedx /= norm2; +// */ - SetdEdx(dedx); - return dedx; -} -Double_t AliTPCseed::Bethe(Double_t bg){ - // - // This is the Bethe-Bloch function normalised to 1 at the minimum - // - Double_t bg2=bg*bg; - Double_t bethe; - if (bg<3.5e1) - bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2)); - else // Density effect ( approximately :) - bethe=1.15*(1.+ bg2)/bg2*(log(3.5*5940*bg) - bg2/(1.+ bg2)); - return bethe/11.091; +// SetdEdx(dedx); +// return dedx; } void AliTPCseed::CookPID() @@ -815,7 +859,7 @@ void AliTPCseed::CookPID() Double_t mass=AliPID::ParticleMass(j); Double_t mom=GetP(); Double_t dedx=fdEdx/fMIP; - Double_t bethe=Bethe(mom/mass); + Double_t bethe=AliMathBase::BetheBlochAleph(mom/mass); Double_t sigma=fRes*bethe; if (sigma>0.001){ if (TMath::Abs(dedx-bethe) > fRange*sigma) { @@ -836,187 +880,15 @@ void AliTPCseed::CookPID() } } -/* -void AliTPCseed::CookdEdx2(Double_t low, Double_t up) { - //----------------------------------------------------------------- - // This funtion calculates dE/dX within the "low" and "up" cuts. - //----------------------------------------------------------------- - - Float_t amp[200]; - Float_t angular[200]; - Float_t weight[200]; - Int_t index[200]; - Bool_t inlimit[200]; - for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE; - for (Int_t i=0;i<200;i++) amp[i]=10000; - for (Int_t i=0;i<200;i++) angular[i]= 1;; - - - // - Float_t meanlog = 100.; - Int_t indexde[4]={0,64,128,160}; - - Float_t amean =0; - Float_t asigma =0; - Float_t anc =0; - Float_t anorm =0; - - Float_t mean[4] = {0,0,0,0}; - Float_t sigma[4] = {1000,1000,1000,1000}; - Int_t nc[4] = {0,0,0,0}; - Float_t norm[4] = {1000,1000,1000,1000}; - // - // - fNShared =0; - - // for (Int_t of =0; of<3; of++){ - // for (Int_t i=indexde[of];ifIsShared){ - fNShared++; - continue; - } - Int_t type = point->GetCPoint().GetType(); - if (type<0) continue; - if (point->GetCPoint().GetMax()<5) continue; - Float_t angley = point->GetTPoint().GetAngleY(); - Float_t anglez = point->GetTPoint().GetAngleZ(); - Float_t rsigmay = point->GetCPoint().GetSigmaY(); - Float_t rsigmaz = point->GetCPoint().GetSigmaZ(); - Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz); - - Float_t ampc = 0; // normalization to the number of electrons - if (i>64){ - ampc = point->GetCPoint().GetMax(); - } - else{ - ampc = point->GetCPoint().GetMax(); - } - ampc *= 2.0; // put mean value to channel 50 - // ampc *= 0.565; // put mean value to channel 50 - - Float_t w = 1.; - Float_t z = TMath::Abs(point->GetCPoint().GetZ()); - if (i<64) { - ampc /= 0.63; - } else - if (i>128){ - ampc /=1.51; - } - if (type<0) { //amp at the border - lower weight - continue; - } - if (rsigma>1.5) ampc/=1.3; // if big backround - angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez); - amp[i] = ampc/angular[i]; - weight[i] = w; - anc++; - } - - TMath::Sort(159,amp,index,kFALSE); - for (Int_t i=int(anc*low+0.5);i0.1) - sigma[of] = TMath::Sqrt(sigma[of]); - else - sigma[of] = 1000; - mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; - } - } - - Float_t dedx =0; - fSdEdx =0; - fMAngular =0; - // - Int_t norm2 = 0; - Int_t norm3 = 0; - Float_t www[3] = {12.,14.,17.}; - //Float_t www[3] = {1.,1.,1.}; - - for (Int_t i =0;i<3;i++){ - if (nc[i]>2&&nc[i]<1000){ - dedx += mean[i] *nc[i]*www[i]/sigma[i]; - fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i]; - fMAngular += norm[i] *nc[i]; - norm2 += nc[i]*www[i]/sigma[i]; - norm3 += (nc[i]-2)*www[i]/sigma[i]; - } - fDEDX[i] = mean[i]; - fSDEDX[i] = sigma[i]; - fNCDEDX[i]= nc[i]; - } - - if (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; - } - else{ - SetdEdx(0); - return; - } - // Float_t dedx1 =dedx; - - dedx =0; - Float_t norm4 = 0; - for (Int_t i =0;i<3;i++){ - if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){ - //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.)); - dedx += mean[i] *(nc[i])/(sigma[i]); - norm4 += (nc[i])/(sigma[i]); - } - fDEDX[i] = mean[i]; - } - if (norm4>0) dedx /= norm4; - - - - SetdEdx(dedx); - - //mi deDX - -} -*/ Double_t AliTPCseed::GetYat(Double_t xk) const { //----------------------------------------------------------------- // This function calculates the Y-coordinate of a track at the plane x=xk. //----------------------------------------------------------------- if (TMath::Abs(GetSnp())>AliTPCReconstructor::GetMaxSnpTrack()) return 0.; //patch 01 jan 06 - Double_t c1=GetSnp(), r1=TMath::Sqrt(1.- c1*c1); + Double_t c1=GetSnp(), r1=TMath::Sqrt((1.-c1)*(1.+c1)); Double_t c2=c1+GetC()*(xk-GetX()); if (TMath::Abs(c2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0; - Double_t r2=TMath::Sqrt(1.- c2*c2); + Double_t r2=TMath::Sqrt((1.-c2)*(1.+c2)); return GetY() + (xk-GetX())*(c1+c2)/(r1+r2); } @@ -1041,7 +913,7 @@ Bool_t AliTPCseed::GetSharedMapBit(int ibit) -Float_t AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, AliTPCCalPad * gainMap, Bool_t posNorm, Bool_t padNorm){ +Float_t AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Int_t posNorm, Int_t padNorm, Int_t returnVal){ // // calculates dedx using the cluster @@ -1049,20 +921,39 @@ Float_t AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i // type - 1 - max charge or 0- total charge in cluster // //2- max no corr 3- total+ correction // i1-i2 - the pad-row range used for calculation - // + // shapeNorm - kTRUE -taken from OCDB + // + // posNorm - usage of pos normalization + // padNorm - pad type normalization + // returnVal - 0 return mean + // - 1 return RMS + // - 2 return number of clusters + // // normalization parametrization taken from AliTPCClusterParam // - AliTPCClusterParam * parcl = AliTPCClusterParam::Instance(); - if (!parcl) parcl = AliTPCcalibDB::Instance()->GetClusterParam(); - if (!parcl) return 0; + AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam(); + AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters(); + if (!parcl) return 0; + if (!param) return 0; + Int_t row0 = param->GetNRowLow(); + Int_t row1 = row0+param->GetNRowUp1(); + Float_t amp[160]; Int_t indexes[160]; Int_t ncl=0; // // - Float_t gainGG = 1; + Float_t gainGG = 1; // gas gain factor -always enabled + Float_t gainPad = 1; // gain map - used always + Float_t corrShape = 1; // correction due angular effect, diffusion and electron attachment + Float_t corrPos = 1; // local position correction - if posNorm enabled + Float_t corrPadType = 1; // pad type correction - if padNorm enabled + Float_t corrNorm = 1; // normalization factor - set Q to channel 50 + // + // + // if (AliTPCcalibDB::Instance()->GetParameters()){ - gainGG= 20000./AliTPCcalibDB::Instance()->GetParameters()->GetGasGain(); //relative gas gain + gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain } const Float_t ktany = TMath::Tan(TMath::DegToRad()*10); @@ -1074,59 +965,96 @@ Float_t AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i if (!cluster) continue; if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ(); - if (!gainMap) gainMap = AliTPCcalibDB::Instance()->GetDedxGainFactor(); - if (gainMap) { + Int_t ipad= 0; + if (irow>=row0) ipad=1; + if (irow>=row1) ipad=2; + // + // + // + AliTPCCalPad * gainMap = AliTPCcalibDB::Instance()->GetDedxGainFactor(); + if (gainMap) { + // + // Get gainPad - pad by pad calibration + // Float_t factor = 1; AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector()); - if (irow < 63) { // IROC - factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()))*1.55; + if (irow < row0) { // IROC + factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad())); } else { // OROC - factor = roc->GetValue(irow - 63, TMath::Nint(cluster->GetPad())); + factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad())); + } + if (factor>0.5) gainPad=factor; + } + // + //do position and angular normalization + // + if (shapeNorm){ + if (type<=1){ + // + AliTPCTrackerPoint * point = GetTrackPoint(irow); + Float_t ty = TMath::Abs(point->GetAngleY()); + Float_t tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty)); + + Float_t dr = (250.-TMath::Abs(cluster->GetZ()))/250.; + corrShape = parcl->Qnorm(ipad,type,dr,ty,tz); } - if (factor>0.5) charge/=factor; } - //do normalization - Float_t corr=1; - Int_t ipad= 0; - if (irow>62) ipad=1; - if (irow>127) ipad=2; - if (type<=1){ - // + if (posNorm>0){ + // + // Do position normalization - relative distance to + // center of pad- time bin + // Work in progress + // corrPos = parcl->QnormPos(ipad,type, cluster->GetPad(), + // cluster->GetTimeBin(), cluster->GetZ(), + // cluster->GetSigmaY2(),cluster->GetSigmaZ2(), + // cluster->GetMax(),cluster->GetQ()); + // scaled response function + Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector()); + Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth(); + // + AliTPCTrackerPoint * point = GetTrackPoint(irow); Float_t ty = TMath::Abs(point->GetAngleY()); - Float_t tz = TMath::Abs(point->GetAngleZ()); + Float_t tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty)); - Float_t dr = (250.-TMath::Abs(cluster->GetZ()))/250.; - corr = parcl->Qnorm(ipad,type,dr,ty,tz); - } - amp[ncl]=charge/corr; - amp[ncl]/gainGG; - if (posNorm){ - // - // - // - corr = parcl->QnormPos(ipad,type, cluster->GetPad(),cluster->GetTimeBin(), cluster->GetZ(), - cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax(),cluster->GetQ()); - amp[ncl]/=corr; + if (type==1) corrPos = + parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), + cluster->GetTimeBin(),ty,tz,yres0,zres0,0.4); + if (type==0) corrPos = + parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), + cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,0.4); + if (posNorm==3){ + Float_t dr = (250.-TMath::Abs(cluster->GetZ()))/250.; + Double_t signtgl = (cluster->GetZ()*point->GetAngleZ()>0)? 1:-1; + Double_t p2 = TMath::Abs(TMath::Sin(TMath::ATan(ty))); + Float_t corrHis = parcl->QnormHis(ipad,type,dr,p2,TMath::Abs(point->GetAngleZ())*signtgl); + if (corrHis>0) corrPos*=corrHis; + } + } + if (padNorm==1){ + //taken from OCDB + if (type==0 && parcl->QpadTnorm()) corrPadType = (*parcl->QpadTnorm())[ipad]; + if (type==1 && parcl->QpadMnorm()) corrPadType = (*parcl->QpadMnorm())[ipad]; - amp[ncl] *= 2.0; // put mean value to channel 50 - if (padNorm){ - corr=1; - if (type==0 && parcl->fQpadTnorm) corr = (*parcl->fQpadTnorm)[ipad]; - if (type==1 && parcl->fQpadTnorm) corr = (*parcl->fQpadMnorm)[ipad]; - amp[ncl]/=corr; } - - // if (ipad==0) { -// amp[ncl] /= 0.65; // this we will take form OCDB -// } else -// if (ipad==2){ -// amp[ncl] /=1.57; -// }else{ -// } + if (padNorm==2){ + corrPadType =param->GetPadPitchLength(cluster->GetDetector(),cluster->GetRow()); + //use hardwired - temp fix + if (type==0) corrNorm=3.; + if (type==1) corrNorm=1.; + } + // + amp[ncl]=charge; + amp[ncl]/=gainGG; + amp[ncl]/=gainPad; + amp[ncl]/=corrShape; + amp[ncl]/=corrPadType; + amp[ncl]/=corrPos; + amp[ncl]/=corrNorm; + // ncl++; } @@ -1136,39 +1064,254 @@ Float_t AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i if (ncl<10) return 0; Float_t suma=0; + Float_t suma2=0; Float_t sumn=0; Int_t icl0=TMath::Nint(ncl*low); Int_t icl1=TMath::Nint(ncl*up); for (Int_t icl=icl0; iclGetTransform(); + const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam(); + if (trans && recoParam->GetUseGainCorrectionTime()>0) { + runNumber = trans->GetCurrentRunNumber(); + //AliTPCcalibDB::Instance()->SetRun(runNumber); + TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber); + if (timeGainSplines) { + UInt_t time = trans->GetCurrentTimeStamp(); + AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0); + AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1); + if (fitMIP) { + corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/ + } else { + if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/ + } + } + } + mean /= corrTimeGain; + rms /= corrTimeGain; + // + if (returnVal==1) return rms; + if (returnVal==2) return ncl; + return mean; } -Double_t AliTPCseed::BetheMass(Double_t mass){ +Float_t AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal, Int_t rowThres){ + + // + // calculates dedx using the cluster + // low - up specify trunc mean range - default form 0-0.7 + // type - 1 - max charge or 0- total charge in cluster + // //2- max no corr 3- total+ correction + // i1-i2 - the pad-row range used for calculation + // + // posNorm - usage of pos normalization + // returnVal - 0 return mean + // - 1 return RMS + // - 2 return number of clusters // - // return bethe-bloch + // rowThres - number of rows before and after given pad row to check for clusters below threshold + // + // normalization parametrization taken from AliTPCClusterParam // - Float_t bg= P()/mass; - const Double_t kp1=0.76176e-1; - const Double_t kp2=10.632; - const Double_t kp3=0.13279e-4; - const Double_t kp4=1.8631; - const Double_t kp5=1.9479; + AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam(); + AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters(); + if (!parcl) return 0; + if (!param) return 0; + Int_t row0 = param->GetNRowLow(); + Int_t row1 = row0+param->GetNRowUp1(); - Double_t dbg = (Double_t) bg; + Float_t amp[160]; + Int_t indexes[160]; + Int_t ncl=0; + Int_t nclBelowThr = 0; // counts number of clusters below threshold + // + // + Float_t gainGG = 1; // gas gain factor -always enabled + Float_t gainPad = 1; // gain map - used always + Float_t corrPos = 1; // local position correction - if posNorm enabled + // + // + // + if (AliTPCcalibDB::Instance()->GetParameters()){ + gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain + } + // + // extract time-dependent correction for pressure and temperature variations + // + UInt_t runNumber = 1; + Float_t corrTimeGain = 1; + TObjArray * timeGainSplines = 0x0; + // + AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform(); + const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam(); + if (trans) { + runNumber = trans->GetCurrentRunNumber(); + //AliTPCcalibDB::Instance()->SetRun(runNumber); + timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber); + if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) { + UInt_t time = trans->GetCurrentTimeStamp(); + AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0); + AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1); + if (fitMIP) { + corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time); /*fitMIP->Eval(time);*/ + } else { + if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time); /*fitFPcosmic->Eval(time); */ + } + } + } + + const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam + const Float_t ktany = TMath::Tan(TMath::DegToRad()*10); + const Float_t kedgey =3.; + // + // + for (Int_t irow=i1; irowGetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster + // + AliTPCTrackerPoint * point = GetTrackPoint(irow); + if (point==0) continue; + Float_t rsigmay = TMath::Sqrt(point->GetSigmaY()); + if (rsigmay > kClusterShapeCut) continue; + // + if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb + // + Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ(); + Int_t ipad= 0; + if (irow>=row0) ipad=1; + if (irow>=row1) ipad=2; + // + // + // + AliTPCCalPad * gainMap = AliTPCcalibDB::Instance()->GetDedxGainFactor(); + if (gainMap) { + // + // Get gainPad - pad by pad calibration + // + Float_t factor = 1; + AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector()); + if (irow < row0) { // IROC + factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad())); + } else { // OROC + factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad())); + } + if (factor>0.3) gainPad=factor; + } + // + // Do position normalization - relative distance to + // center of pad- time bin + + Float_t ty = TMath::Abs(point->GetAngleY()); + Float_t tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty)); + Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector()); + Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth(); + + yres0 *=parcl->GetQnormCorr(ipad, type,0); + zres0 *=parcl->GetQnormCorr(ipad, type,1); + Float_t effLength=parcl->GetQnormCorr(ipad, type,4)*0.5; + Float_t effDiff =(parcl->GetQnormCorr(ipad, type,2)+parcl->GetQnormCorr(ipad, type,3))*0.5; + // + if (type==1) { + corrPos = parcl->GetQnormCorr(ipad, type,5)* + parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), + cluster->GetTimeBin(),ty,tz,yres0,zres0,effLength,effDiff); + Float_t drm = 0.5-TMath::Abs(cluster->GetZ()/250.); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz); + // + } + if (type==0) { + corrPos = parcl->GetQnormCorr(ipad, type,5)* + parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), + cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,effLength,effDiff); + + Float_t drm = 0.5-TMath::Abs(cluster->GetZ()/250.); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz); + // + } + // + // pad region equalization outside of cluster param + // + Float_t gainEqualPadRegion = 1; + if (timeGainSplines) { //1 - max charge or 0- total charge in cluster + TGraphErrors * grPadEqual = 0x0; + if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL"); + if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL"); + if (grPadEqual) gainEqualPadRegion = grPadEqual->Eval(ipad); + } + // + amp[ncl]=charge; + amp[ncl]/=gainGG; + amp[ncl]/=gainPad; + amp[ncl]/=corrPos; + amp[ncl]/=gainEqualPadRegion; + // + ncl++; + } - Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg); + if (type>3) return ncl; + TMath::Sort(ncl,amp, indexes, kFALSE); + // + if (ncl<10) return 0; + // + Double_t * ampWithBelow = new Double_t[ncl + nclBelowThr]; + for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) { + if (iCl < nclBelowThr) { + ampWithBelow[iCl] = amp[indexes[0]]; + } else { + ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]]; + } + } + //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]); + // + Float_t suma=0; + Float_t suma2=0; + Float_t sumn=0; + Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low); + Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up); + for (Int_t icl=icl0; icl0)? means/meanc:0; return mean; } + + + +Int_t AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, AliExternalTrackParam * parout){ + // + // Refit the track + // return value - number of used clusters + // + // + const Int_t kMinNcl =10; + AliTPCseed *track=new AliTPCseed(*seed); + Int_t sector=-1; + // reset covariance + // + Double_t covar[15]; + for (Int_t i=0;i<15;i++) covar[i]=0; + covar[0]=10.*10.; + covar[2]=10.*10.; + covar[5]=10.*10./(64.*64.); + covar[9]=10.*10./(64.*64.); + covar[14]=1*1; + // + + Float_t xmin=1000, xmax=-10000; + Int_t imin=158, imax=0; + for (Int_t i=0;i<160;i++) { + AliTPCclusterMI *c=track->GetClusterPointer(i); + if (!c) continue; + if (sector<0) sector = c->GetDetector(); + if (c->GetX()GetX(); + if (c->GetX()>xmax) xmax=c->GetX(); + if (iimax) imax=i; + } + if(imax-iminRotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) { + delete track; + return 0; + } + // + // + // fit from inner to outer row + // + AliExternalTrackParam paramIn; + AliExternalTrackParam paramOut; + Bool_t isOK=kTRUE; + Int_t ncl=0; + // + // + // + for (Int_t i=imin; i<=imax; i++){ + AliTPCclusterMI *c=track->GetClusterPointer(i); + if (!c) continue; + // if (RejectCluster(c,track)) continue; + sector = (c->GetDetector()%18); + if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) { + //continue; + } + Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; + Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation + if (!track->PropagateTo(r[0])) { + isOK=kFALSE; + } + if ( !((static_cast(track)->Update(&r[1],cov)))) isOK=kFALSE; + } + if (!isOK) { delete track; return 0;} + track->AddCovariance(covar); + // + // + // + for (Int_t i=imax; i>=imin; i--){ + AliTPCclusterMI *c=track->GetClusterPointer(i); + if (!c) continue; + //if (RejectCluster(c,track)) continue; + sector = (c->GetDetector()%18); + if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) { + //continue; + } + Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; + Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation + if (!track->PropagateTo(r[0])) { + isOK=kFALSE; + } + if ( !((static_cast(track)->Update(&r[1],cov)))) isOK=kFALSE; + } + //if (!isOK) { delete track; return 0;} + paramIn = *track; + track->AddCovariance(covar); + // + // + for (Int_t i=imin; i<=imax; i++){ + AliTPCclusterMI *c=track->GetClusterPointer(i); + if (!c) continue; + sector = (c->GetDetector()%18); + if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) { + //continue; + } + ncl++; + //if (RejectCluster(c,track)) continue; + Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; + Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation + if (!track->PropagateTo(r[0])) { + isOK=kFALSE; + } + if ( !((static_cast(track)->Update(&r[1],cov)))) isOK=kFALSE; + } + //if (!isOK) { delete track; return 0;} + paramOut=*track; + // + // + // + if (parin) (*parin)=paramIn; + if (parout) (*parout)=paramOut; + delete track; + return ncl; +} + + + +Bool_t AliTPCseed::RefitTrack(AliTPCseed* /*seed*/, Bool_t /*out*/){ + // + // + // + return kFALSE; +} + + + + + + +void AliTPCseed::GetError(AliTPCclusterMI* cluster, AliExternalTrackParam * param, + Double_t& erry, Double_t &errz) +{ + // + // Get cluster error at given position + // + AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam(); + Double_t tany,tanz; + Double_t snp1=param->GetSnp(); + tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Double_t tgl1=param->GetTgl(); + tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Int_t padSize = 0; // short pads + if (cluster->GetDetector() >= 36) { + padSize = 1; // medium pads + if (cluster->GetRow() > 63) padSize = 2; // long pads + } + + erry = clusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany) ); + errz = clusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ); +} + + +void AliTPCseed::GetShape(AliTPCclusterMI* cluster, AliExternalTrackParam * param, + Double_t& rmsy, Double_t &rmsz) +{ + // + // Get cluster error at given position + // + AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam(); + Double_t tany,tanz; + Double_t snp1=param->GetSnp(); + tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Double_t tgl1=param->GetTgl(); + tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Int_t padSize = 0; // short pads + if (cluster->GetDetector() >= 36) { + padSize = 1; // medium pads + if (cluster->GetRow() > 63) padSize = 2; // long pads + } + + rmsy = clusterParam->GetRMSQ( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany), TMath::Abs(cluster->GetMax()) ); + rmsz = clusterParam->GetRMSQ( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ,TMath::Abs(cluster->GetMax())); +} + + + +Double_t AliTPCseed::GetQCorrGeom(Float_t ty, Float_t tz){ + //Geoetrical + //ty - tangent in local y direction + //tz - tangent + // + Float_t norm=TMath::Sqrt(1+ty*ty+tz*tz); + return norm; +} + +Double_t AliTPCseed::GetQCorrShape(Int_t ipad, Int_t type,Float_t z, Float_t ty, Float_t tz, Float_t /*q*/, Float_t /*thr*/){ + // + // Q normalization + // + // return value = Q Normalization factor + // Normalization - 1 - shape factor part for full drift + // 1 - electron attachment for 0 drift + + // Input parameters: + // + // ipad - 0 short pad + // 1 medium pad + // 2 long pad + // + // type - 0 qmax + // - 1 qtot + // + //z - z position (-250,250 cm) + //ty - tangent in local y direction + //tz - tangent + // + + AliTPCClusterParam * paramCl = AliTPCcalibDB::Instance()->GetClusterParam(); + AliTPCParam * paramTPC = AliTPCcalibDB::Instance()->GetParameters(); + + if (!paramCl) return 1; + // + Double_t dr = 250.-TMath::Abs(z); + Double_t sy = paramCl->GetRMS0( 0,ipad, dr, TMath::Abs(ty)); + Double_t sy0= paramCl->GetRMS0(0,ipad, 250, 0); + Double_t sz = paramCl->GetRMS0( 1,ipad, dr, TMath::Abs(tz)); + Double_t sz0= paramCl->GetRMS0(1,ipad, 250, 0); + + Double_t sfactorMax = TMath::Sqrt(sy0*sz0/(sy*sz)); + + + Double_t dt = 1000000*(dr/paramTPC->GetDriftV()); //time in microsecond + Double_t attProb = TMath::Exp(-paramTPC->GetAttCoef()*paramTPC->GetOxyCont()*dt); + // + // + if (type==0) return sfactorMax*attProb; + + return attProb; + + +} +