X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCseed.cxx;h=61543888d4d178147f783c1983e6739061601726;hb=8f1b4e07c12a9ba77f9ff0bce824dd2a6975253d;hp=b7eeb526f2ddac6b2ff236779d6757f96db4e8e9;hpb=465c437ba393c3f5760ddcf54dd3ee3941ce7996;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCseed.cxx b/TPC/AliTPCseed.cxx index b7eeb526f2d..61543888d4d 100644 --- a/TPC/AliTPCseed.cxx +++ b/TPC/AliTPCseed.cxx @@ -25,6 +25,16 @@ #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" + + ClassImp(AliTPCseed) @@ -39,6 +49,11 @@ AliTPCseed::AliTPCseed(): fRelativeSector(-1), fCurrentSigmaY2(1e10), fCurrentSigmaZ2(1e10), + fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30% + fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30% + fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2% + fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2% + // fErrorY2(1e10), fErrorZ2(1e10), fCurrentCluster(0x0), @@ -81,8 +96,12 @@ AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner): fRow(0), fSector(-1), fRelativeSector(-1), - fCurrentSigmaY2(1e10), - fCurrentSigmaZ2(1e10), + fCurrentSigmaY2(-1), + fCurrentSigmaZ2(-1), + fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30% + fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30% + fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2% + fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2% fErrorY2(1e10), fErrorZ2(1e10), fCurrentCluster(0x0), @@ -132,8 +151,12 @@ AliTPCseed::AliTPCseed(const AliTPCtrack &t): fRow(0), fSector(-1), fRelativeSector(-1), - fCurrentSigmaY2(1e10), - fCurrentSigmaZ2(1e10), + fCurrentSigmaY2(-1), + fCurrentSigmaZ2(-1), + fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30% + fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30% + fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2% + fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2% fErrorY2(1e10), fErrorZ2(1e10), fCurrentCluster(0x0), @@ -188,8 +211,12 @@ AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5], fRow(0), fSector(-1), fRelativeSector(-1), - fCurrentSigmaY2(1e10), - fCurrentSigmaZ2(1e10), + fCurrentSigmaY2(-1), + fCurrentSigmaZ2(-1), + fCMeanSigmaY2p30(-1.), //! current mean sigma Y2 - mean30% + fCMeanSigmaZ2p30(-1.), //! current mean sigma Z2 - mean30% + fCMeanSigmaY2p30R(-1.), //! current mean sigma Y2 - mean2% + fCMeanSigmaZ2p30R(-1.), //! current mean sigma Z2 - mean2% fErrorY2(1e10), fErrorZ2(1e10), fCurrentCluster(0x0), @@ -419,8 +446,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(); @@ -457,6 +484,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); } @@ -506,7 +542,31 @@ 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){ + fCMeanSigmaY2p30= c->GetSigmaY2(); //! current mean sigma Y2 - mean30% + fCMeanSigmaZ2p30= c->GetSigmaZ2(); //! current mean sigma Z2 - mean30% + fCMeanSigmaY2p30R = 1; //! current mean sigma Y2 - mean5% + fCMeanSigmaZ2p30R = 1; //! current mean sigma Z2 - mean5% + } + // + fCMeanSigmaY2p30= 0.70*fCMeanSigmaY2p30 +0.30*c->GetSigmaY2(); + fCMeanSigmaZ2p30= 0.70*fCMeanSigmaZ2p30 +0.30*c->GetSigmaZ2(); + if (fCurrentSigmaY2>0){ + fCMeanSigmaY2p30R = 0.7*fCMeanSigmaY2p30R +0.3*c->GetSigmaY2()/fCurrentSigmaY2; + fCMeanSigmaZ2p30R = 0.7*fCMeanSigmaZ2p30R +0.3*c->GetSigmaZ2()/fCurrentSigmaZ2; + } + SetClusterIndex(n,idx); // restore the current cluster index return kTRUE; @@ -515,245 +575,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. //----------------------------------------------------------------- + // CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal) + AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); + + 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); + // + 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)); + + SetdEdx(fDEDX[0]); + return fDEDX[0]; + +// return CookdEdxNorm(low,up,0,i1,i2,1,0,2); + - 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 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; - - 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 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]] = 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() @@ -771,7 +858,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) { @@ -792,187 +879,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); } @@ -997,7 +912,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){ +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 @@ -1005,48 +920,140 @@ 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) 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; // 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= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain + } + const Float_t ktany = TMath::Tan(TMath::DegToRad()*10); - const Float_t kedgey =4.; + const Float_t kedgey =3.; + // // for (Int_t irow=i1; irowGetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ(); - //do normalization - Float_t corr=1; Int_t ipad= 0; - if (irow>62) ipad=1; - if (irow>127) ipad=2; - if (type<=1){ - // + 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.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 (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); + 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; + } + } - amp[ncl]=charge/corr; - - amp[ncl] *= 2.0; // put mean value to channel 50 - if (ipad==0) { - amp[ncl] /= 0.65; // this we will take form OCDB - } else - if (ipad==2){ - amp[ncl] /=1.57; - }else{ - } + + if (padNorm==1){ + //taken from OCDB + if (type==0 && parcl->fQpadTnorm) corrPadType = (*parcl->fQpadTnorm)[ipad]; + if (type==1 && parcl->fQpadTnorm) corrPadType = (*parcl->fQpadMnorm)[ipad]; + + } + 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++; } @@ -1056,41 +1063,210 @@ 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 = fitMIP->Eval(time); + } else { + if (fitFPcosmic) corrTimeGain = fitFPcosmic->Eval(time); // This value describes the ratio FP-to-MIP, hardwired for the moment + } + } + } + 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){ + + // + // 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 + // + // normalization parametrization taken from AliTPCClusterParam + // + 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; // - // return bethe-bloch // - 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; + 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 - Double_t dbg = (Double_t) bg; + // + // + // + if (AliTPCcalibDB::Instance()->GetParameters()){ + gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain + } + + 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 + 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.5) gainPad=factor; + } + + // + // Do position normalization - relative distance to + // center of pad- time bin + + 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 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); + // + } - Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg); + // + amp[ncl]=charge; + amp[ncl]/=gainGG; + amp[ncl]/=gainPad; + amp[ncl]/=corrPos; + // + ncl++; + } - Double_t aa = TMath::Power(beta,kp4); - Double_t bb = TMath::Power(1./dbg,kp5); + if (type>3) return ncl; + TMath::Sort(ncl,amp, indexes, kFALSE); - bb=TMath::Log(kp3+bb); + if (ncl<10) return 0; - return ((Float_t)((kp2-aa-bb)*kp1/aa)); + 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 = fitMIP->Eval(time); + } else { + if (fitFPcosmic) corrTimeGain = fitFPcosmic->Eval(time); // This value describes the ratio FP-to-MIP, hardwired for the moment + } + } + } + mean /= corrTimeGain; + rms /= corrTimeGain; + // + if (returnVal==1) return rms; + if (returnVal==2) return ncl; + return mean; } + + Float_t AliTPCseed::CookShape(Int_t type){ // // @@ -1118,3 +1294,244 @@ Float_t AliTPCseed::CookShape(Int_t type){ Float_t mean = (meanc>0)? 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; + 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; + + +} +