X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCseed.cxx;h=3109109ff7a61b71f0e8f0b2923b6c8fc148ffe7;hb=f8b468657705805f0c80bc532ebeedc2823f7103;hp=c54f3e9daa854498f90aae207581d52195053796;hpb=3f82c4f228f0e53fc1421725b00f791cae0f5929;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCseed.cxx b/TPC/AliTPCseed.cxx index c54f3e9daa8..3109109ff7a 100644 --- a/TPC/AliTPCseed.cxx +++ b/TPC/AliTPCseed.cxx @@ -24,58 +24,160 @@ #include "TClonesArray.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" + + ClassImp(AliTPCseed) -AliTPCseed::AliTPCseed():AliTPCtrack(){ +AliTPCseed::AliTPCseed(): + AliTPCtrack(), + fEsd(0x0), + fClusterOwner(kFALSE), + fRow(0), + fSector(-1), + 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% // - fRow=0; - fRemoval =0; - for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); + fErrorY2(1e10), + fErrorZ2(1e10), + fCurrentCluster(0x0), + fCurrentClusterIndex1(-1), + fInDead(kFALSE), + fIsSeeding(kFALSE), + fNoCluster(0), + fSort(0), + fBSigned(kFALSE), + fSeedType(0), + fSeed1(-1), + fSeed2(-1), + fMAngular(0), + fCircular(0), + fClusterMap(159), + fSharedMap(159) +{ + // + for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3); for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0; - for (Int_t i=0;i<5;i++) fTPCr[i]=0.2; - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared =0; - fRemoval = 0; - fSort =0; - fFirstPoint =0; - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; - fCurrentSigmaY2=0; - fCurrentSigmaZ2=0; - fEsd =0; - fCircular = 0; // not curling track + for (Int_t i=0;iUncheckedAt(i); - if (cl0){ - trpoint->GetTPoint() = *(GetTrackPoint(i)); - trpoint->GetCPoint() = *cl0; - trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ())); - } - else{ - *trpoint = pdummy; - trpoint->GetCPoint()= cldummy; + // assignment operator + // + if(this!=¶m){ + AliTPCtrack::operator=(param); + fEsd =param.fEsd; + for(Int_t i = 0;i<160;++i)fClusterPointer[i] = param.fClusterPointer[i]; // this is not allocated by AliTPCSeed + fClusterOwner = param.fClusterOwner; + // leave out fPoint, they are also not copied in the copy ctor... + // but deleted in the dtor... strange... + fRow = param.fRow; + fSector = param.fSector; + fRelativeSector = param.fRelativeSector; + fCurrentSigmaY2 = param.fCurrentSigmaY2; + fCurrentSigmaZ2 = param.fCurrentSigmaZ2; + fErrorY2 = param.fErrorY2; + fErrorZ2 = param.fErrorZ2; + fCurrentCluster = param.fCurrentCluster; // this is not allocated by AliTPCSeed + fCurrentClusterIndex1 = param.fCurrentClusterIndex1; + fInDead = param.fInDead; + fIsSeeding = param.fIsSeeding; + fNoCluster = param.fNoCluster; + fSort = param.fSort; + fBSigned = param.fBSigned; + for(Int_t i = 0;i<4;++i){ + fDEDX[i] = param.fDEDX[i]; + fSDEDX[i] = param.fSDEDX[i]; + fNCDEDX[i] = param.fNCDEDX[i]; } + for(Int_t i = 0;i= AliTPCReconstructor::GetMaxSnpTrack()) { + if (TMath::Abs(GetSnp()+GetC()*dx) >= AliTPCReconstructor::GetMaxSnpTrack()) { return 0; } // Double_t y1=fP0, z1=fP1; - Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1); - Double_t c2=fP4*x2 - fP2, 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 = fP0; - z = fP1; + y = GetY(); + z = GetZ(); //y += dx*(c1+c2)/(r1+r2); //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3; Double_t dy = dx*(c1+c2)/(r1+r2); Double_t dz = 0; // - Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1); + Double_t delta = GetC()*dx*(c1+c2)/(c1*r2 + c2*r1); /* if (TMath::Abs(delta)>0.0001){ dz = fP3*TMath::ASin(delta)/fP4; @@ -326,7 +466,7 @@ Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const } */ // dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4; - dz = fP3*TMath::ASin(delta)/fP4; + dz = GetTgl()*TMath::ASin(delta)/GetC(); // y+=dy; z+=dz; @@ -342,24 +482,20 @@ Double_t AliTPCseed::GetPredictedChi2(const AliCluster *c) const //----------------------------------------------------------------- // This function calculates a predicted chi2 increment. //----------------------------------------------------------------- - //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2(); - Double_t r00=fErrorY2, r01=0., r11=fErrorZ2; - r00+=fC00; r01+=fC10; r11+=fC11; - - Double_t det=r00*r11 - r01*r01; - if (TMath::Abs(det) < 1.e-10) { - //Int_t n=GetNumberOfClusters(); - //if (n>4) cerr<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; } - Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01; - - Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1; - - return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; + return AliExternalTrackParam::GetPredictedChi2(p,cov); } - //_________________________________________________________________________________________ @@ -380,11 +516,11 @@ Int_t AliTPCseed::Compare(const TObject *o) const { } else { Float_t f2 =1; - f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066); + f2 = 1-20*TMath::Sqrt(t->GetSigma1Pt2())/(t->OneOverPt()+0.0066); if (t->fBConstrain) f2=1.2; Float_t f1 =1; - f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066); + f1 = 1-20*TMath::Sqrt(GetSigma1Pt2())/(OneOverPt()+0.0066); if (fBConstrain) f1=1.2; @@ -397,304 +533,313 @@ Int_t AliTPCseed::Compare(const TObject *o) const { //_____________________________________________________________________________ -Int_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, UInt_t /*index*/) { +Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_t index) +{ //----------------------------------------------------------------- // This function associates a cluster with this track. //----------------------------------------------------------------- - Double_t r00=fErrorY2, r01=0., r11=fErrorZ2; - - r00+=fC00; r01+=fC10; r11+=fC11; - Double_t det=r00*r11 - r01*r01; - Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; - - Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11; - Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11; - Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11; - Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11; - Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11; - - Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1; - Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz; - if (TMath::Abs(cur*fX-eta) >= AliTPCReconstructor::GetMaxSnpTrack()) { - return 0; + Int_t n=GetNumberOfClusters(); + 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); } - fP0 += k00*dy + k01*dz; - fP1 += k10*dy + k11*dz; - fP2 = eta; - fP3 += k30*dy + k31*dz; - fP4 = cur; - - Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40; - Double_t c12=fC21, c13=fC31, c14=fC41; - - fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11; - fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13; - fC40-=k00*c04+k01*c14; - - fC11-=k10*c01+k11*fC11; - fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13; - fC41-=k10*c04+k11*c14; - - fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13; - fC42-=k20*c04+k21*c14; - - fC33-=k30*c03+k31*c13; - fC43-=k40*c03+k41*c13; - - fC44-=k40*c04+k41*c14; + 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; + } - Int_t n=GetNumberOfClusters(); - // fIndex[n]=index; - SetNumberOfClusters(n+1); - SetChi2(GetChi2()+chisq); - return 1; + SetClusterIndex(n,idx); // restore the current cluster index + return kTRUE; } //_____________________________________________________________________________ -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. //----------------------------------------------------------------- + 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] = CookdEdxNorm(low,up,useTot ,i1 ,i2, kTRUE,kFALSE,2,0); + fDEDX[1] = CookdEdxNorm(low,up,useTot ,0 ,row0,kTRUE,kFALSE,2,0); + fDEDX[2] = CookdEdxNorm(low,up,useTot ,row0,row1,kTRUE,kFALSE,2,0); + fDEDX[3] = CookdEdxNorm(low,up,useTot ,row1,row2,kTRUE,kFALSE,2,0); + // + fSDEDX[0] = CookdEdxNorm(low,up,useTot ,i1 ,i2, kTRUE,kFALSE,2,1); + fSDEDX[1] = CookdEdxNorm(low,up,useTot ,0 ,row0,kTRUE,kFALSE,2,1); + fSDEDX[2] = CookdEdxNorm(low,up,useTot ,row0,row1,kTRUE,kFALSE,2,1); + fSDEDX[3] = CookdEdxNorm(low,up,useTot ,row1,row2,kTRUE,kFALSE,2,1); + // + fNCDEDX[0] = TMath::Nint(CookdEdxNorm(low,up,useTot ,i1 ,i2, kTRUE,kFALSE,2,2)); + fNCDEDX[1] = TMath::Nint(CookdEdxNorm(low,up,useTot ,0 ,row0,kTRUE,kFALSE,2,2)); + fNCDEDX[2] = TMath::Nint(CookdEdxNorm(low,up,useTot ,row0,row1,kTRUE,kFALSE,2,2)); + fNCDEDX[3] = TMath::Nint(CookdEdxNorm(low,up,useTot ,row1,row2,kTRUE,kFALSE,2,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; - // TClonesArray & arr = *fPoints; - 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() @@ -710,9 +855,9 @@ void AliTPCseed::CookPID() Double_t sumr =0; for (Int_t j=0; j0.001){ if (TMath::Abs(dedx-bethe) > fRange*sigma) { @@ -733,175 +878,657 @@ void AliTPCseed::CookPID() } } -/* -void AliTPCseed::CookdEdx2(Double_t low, Double_t up) { - //----------------------------------------------------------------- - // This funtion calculates dE/dX within the "low" and "up" cuts. - //----------------------------------------------------------------- +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)*(1.+c1)); + Double_t c2=c1+GetC()*(xk-GetX()); + if (TMath::Abs(c2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0; + Double_t r2=TMath::Sqrt((1.-c2)*(1.+c2)); + return GetY() + (xk-GetX())*(c1+c2)/(r1+r2); +} + +void AliTPCseed::SetClusterMapBit(int ibit, Bool_t state) +{ + fClusterMap[ibit] = state; +} +Bool_t AliTPCseed::GetClusterMapBit(int ibit) +{ + return fClusterMap[ibit]; +} +void AliTPCseed::SetSharedMapBit(int ibit, Bool_t state) +{ + fSharedMap[ibit] = state; +} +Bool_t AliTPCseed::GetSharedMapBit(int ibit) +{ + return fSharedMap[ibit]; +} + + + - 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 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 + // 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 + // 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 = 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 =3.; + // // - 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++; + 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 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()*TMath::Sqrt(1+ty*ty)); + + 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; + } + } - TMath::Sort(159,amp,index,kFALSE); - for (Int_t i=int(anc*low+0.5);ifQpadTnorm) 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++; } + + if (type>3) return ncl; + TMath::Sort(ncl,amp, indexes, kFALSE); + + if (ncl<10) return 0; - // meanlog = amp[index[Int_t(anc*0.3)]]; - meanlog =10000.; - for (Int_t of =0; of<3; of++){ - Float_t sumamp=0; - Float_t sumamp2=0; - Float_t sumw=0; - for (Int_t i=indexde[of];i0.1) - sigma[of] = TMath::Sqrt(sigma[of]); - else - sigma[of] = 1000; - mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; - } + 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; icl2&&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]; + Float_t mean =suma/sumn; + Float_t rms =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean)); + // + // do time-dependent correction for pressure and temperature variations + UInt_t runNumber = 1; + Float_t corrTimeGain = 1; + AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform(); + if (trans) { + 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 + } } - fDEDX[i] = mean[i]; - fSDEDX[i] = sigma[i]; - fNCDEDX[i]= nc[i]; } + mean /= corrTimeGain; + rms /= corrTimeGain; + // + if (returnVal==1) return rms; + if (returnVal==2) return ncl; + return mean; +} + +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; + // + // + 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 (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; + // + // + // + if (AliTPCcalibDB::Instance()->GetParameters()){ + gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain } - else{ - SetdEdx(0); - return; + + 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); + // + } + + // + amp[ncl]=charge; + amp[ncl]/=gainGG; + amp[ncl]/=gainPad; + amp[ncl]/=corrPos; + // + ncl++; } - // Float_t dedx1 =dedx; + + if (type>3) return ncl; + TMath::Sort(ncl,amp, indexes, kFALSE); + + if (ncl<10) return 0; - 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]); + 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(); + if (trans) { + 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 + } } - fDEDX[i] = mean[i]; } - if (norm4>0) dedx /= norm4; - + mean /= corrTimeGain; + rms /= corrTimeGain; + // + if (returnVal==1) return rms; + if (returnVal==2) return ncl; + return mean; +} - - SetdEdx(dedx); + + + +Float_t AliTPCseed::CookShape(Int_t type){ + // + // + // + //----------------------------------------------------------------- + // This funtion calculates dE/dX within the "low" and "up" cuts. + //----------------------------------------------------------------- + Float_t means=0; + Float_t meanc=0; + for (Int_t i =0; i<160;i++) { + AliTPCTrackerPoint * point = GetTrackPoint(i); + if (point==0) continue; + + AliTPCclusterMI * cl = fClusterPointer[i]; + if (cl==0) continue; - //mi deDX + Float_t rsigmay = TMath::Sqrt(point->GetSigmaY()); + Float_t rsigmaz = TMath::Sqrt(point->GetSigmaZ()); + Float_t rsigma = (rsigmay+rsigmaz)*0.5; + if (type==0) means+=rsigma; + if (type==1) means+=rsigmay; + if (type==2) means+=rsigmaz; + meanc++; + } + 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; + + +} +