X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCseed.cxx;h=55d89d5df819919997470f8c7765abd3099c7ecb;hb=4ab8985c1f7fbe9b9642049b922752f866ac00a2;hp=77b7b39fa9ae59c9bb541c188479657120459584;hpb=6c23ffeda23cbfdabd4e870e0591de923b08a9a0;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCseed.cxx b/TPC/AliTPCseed.cxx index 77b7b39fa9a..55d89d5df81 100644 --- a/TPC/AliTPCseed.cxx +++ b/TPC/AliTPCseed.cxx @@ -24,6 +24,17 @@ #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) @@ -33,13 +44,16 @@ AliTPCseed::AliTPCseed(): AliTPCtrack(), fEsd(0x0), fClusterOwner(kFALSE), - fPoints(0x0), - fEPoints(0x0), 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% + // fErrorY2(1e10), fErrorZ2(1e10), fCurrentCluster(0x0), @@ -68,21 +82,26 @@ AliTPCseed::AliTPCseed(): fNCDEDX[i] = 0; } for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1; - for (Int_t i=0;i<160;i++) fClusterMap[i]=kFALSE; - for (Int_t i=0;i<160;i++) fSharedMap[i]=kFALSE; + // for (Int_t i=0;i<160;i++) fClusterMap[i]=kFALSE; + //for (Int_t i=0;i<160;i++) fSharedMap[i]=kFALSE; + fClusterMap.ResetAllBits(kFALSE); + fSharedMap.ResetAllBits(kFALSE); + } AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner): AliTPCtrack(s), fEsd(0x0), fClusterOwner(clusterOwner), - fPoints(0x0), - fEPoints(0x0), 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), @@ -129,13 +148,15 @@ AliTPCseed::AliTPCseed(const AliTPCtrack &t): AliTPCtrack(t), fEsd(0x0), fClusterOwner(kFALSE), - fPoints(0x0), - fEPoints(0x0), 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), @@ -174,8 +195,12 @@ AliTPCseed::AliTPCseed(const AliTPCtrack &t): fNCDEDX[i] = 0; } for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1; - for (Int_t i=0;i<160;i++) fClusterMap[i]=kFALSE; - for (Int_t i=0;i<160;i++) fSharedMap[i]=kFALSE; + + //for (Int_t i=0;i<160;i++) fClusterMap[i]=kFALSE; + //for (Int_t i=0;i<160;i++) fSharedMap[i]=kFALSE; + fClusterMap.ResetAllBits(kFALSE); + fSharedMap.ResetAllBits(kFALSE); + } AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5], @@ -183,13 +208,15 @@ AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5], AliTPCtrack(xr, alpha, xx, cc, index), fEsd(0x0), fClusterOwner(kFALSE), - fPoints(0x0), - fEPoints(0x0), 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), @@ -225,26 +252,58 @@ AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5], AliTPCseed::~AliTPCseed(){ // // destructor - if (fPoints) delete fPoints; - fPoints =0; - if (fEPoints) delete fEPoints; - fEPoints = 0; fNoCluster =0; if (fClusterOwner){ for (Int_t icluster=0; icluster<160; icluster++){ delete fClusterPointer[icluster]; } } - for (Int_t i=0;i<160;i++) fClusterMap[i]=kFALSE; - for (Int_t i=0;i<160;i++) fSharedMap[i]=kFALSE; + } //_________________________________________________ -AliTPCseed & AliTPCseed::operator =(const AliTPCseed & param) +AliTPCseed & AliTPCseed::operator=(const AliTPCseed ¶m) { // - // assignment operator - dummy + // assignment operator // - fRow=param.fRow; + 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;iUncheckedAt(i); - if (cl0){ - trpoint->GetTPoint() = *(GetTrackPoint(i)); - trpoint->GetCPoint() = *cl0; - trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ())); - } - else{ - *trpoint = pdummy; - trpoint->GetCPoint()= cldummy; - } - - } - -} Double_t AliTPCseed::GetDensityFirst(Int_t n) @@ -309,6 +344,7 @@ void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_ for (Int_t i=first;iGetY(), 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); } @@ -488,267 +533,314 @@ Int_t AliTPCseed::Compare(const TObject *o) const { //_____________________________________________________________________________ -Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_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 p[2]={c->GetY(), c->GetZ()}; - Double_t cov[3]={fErrorY2, 0., fErrorZ2}; + 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); + } - if (!AliExternalTrackParam::Update(p,cov)) return kFALSE; + 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); + 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. //----------------------------------------------------------------- + // 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)); - 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.; + 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 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() @@ -766,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) { @@ -787,203 +879,659 @@ 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); +} - 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;; - +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 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);iQpadTnorm()) corrPadType = (*parcl->QpadTnorm())[ipad]; + if (type==1 && parcl->QpadMnorm()) corrPadType = (*parcl->QpadMnorm())[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(); + 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 + } } - 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; +} - if (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; +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 (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(); + 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 + } } - 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; } -*/ -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 c2=c1+GetC()*(xk-GetX()); - if (TMath::Abs(c2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0; - Double_t r2=TMath::Sqrt(1.- c2*c2); - return GetY() + (xk-GetX())*(c1+c2)/(r1+r2); + + + +Bool_t AliTPCseed::RefitTrack(AliTPCseed* /*seed*/, Bool_t /*out*/){ + // + // + // + return kFALSE; } -void AliTPCseed::SetClusterMapBit(int ibit, Bool_t state) + + + + + +void AliTPCseed::GetError(AliTPCclusterMI* cluster, AliExternalTrackParam * param, + Double_t& erry, Double_t &errz) { - fClusterMap[ibit] = state; + // + // 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) ); } -Bool_t AliTPCseed::GetClusterMapBit(int ibit) + + +void AliTPCseed::GetShape(AliTPCclusterMI* cluster, AliExternalTrackParam * param, + Double_t& rmsy, Double_t &rmsz) { - return fClusterMap[ibit]; + // + // 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())); } -void AliTPCseed::SetSharedMapBit(int ibit, Bool_t state) -{ - fSharedMap[ibit] = state; + + + +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; } -Bool_t AliTPCseed::GetSharedMapBit(int ibit) -{ - return fSharedMap[ibit]; + +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; + + } +