X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCseed.cxx;h=9bbe38f9caaa4ea78484c961d672c6e42ab2aea7;hb=5e053cadc0ff2b182fc581c52af2b97020b17e4b;hp=bdfaa40f72ce95c335d014e2c4e68a420adac0e3;hpb=b96715746b109eb6dd6231a91c1abc8bac34ef89;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCseed.cxx b/TPC/AliTPCseed.cxx index bdfaa40f72c..9bbe38f9caa 100644 --- a/TPC/AliTPCseed.cxx +++ b/TPC/AliTPCseed.cxx @@ -17,13 +17,26 @@ //----------------------------------------------------------------- +// // Implementation of the TPC seed class // This class is used by the AliTPCtrackerMI class // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch //----------------------------------------------------------------- #include "TClonesArray.h" +#include "TGraphErrors.h" #include "AliTPCseed.h" #include "AliTPCReconstructor.h" +#include "AliTPCClusterParam.h" +#include "AliTPCCalPad.h" +#include "AliTPCCalROC.h" +#include "AliTPCcalibDB.h" +#include "AliTPCParam.h" +#include "AliMathBase.h" +#include "AliTPCTransform.h" +#include "AliSplineFit.h" +#include "AliCDBManager.h" +#include "AliTPCcalibDButil.h" + ClassImp(AliTPCseed) @@ -33,13 +46,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), @@ -53,7 +69,8 @@ AliTPCseed::AliTPCseed(): fSeed1(-1), fSeed2(-1), fMAngular(0), - fCircular(0) + fCircular(0), + fPoolID(-1) { // for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3); @@ -64,7 +81,9 @@ AliTPCseed::AliTPCseed(): fDEDX[i] = 0.; fSDEDX[i] = 1e10; fNCDEDX[i] = 0; + fNCDEDXInclThres[i] = 0; } + fDEDX[4] = 0; for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1; } @@ -72,13 +91,15 @@ 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), @@ -92,7 +113,8 @@ AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner): fSeed1(-1), fSeed2(-1), fMAngular(0), - fCircular(0) + fCircular(0), + fPoolID(-1) { //--------------------- // dummy copy constructor @@ -113,8 +135,11 @@ AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner): fDEDX[i] = s.fDEDX[i]; fSDEDX[i] = s.fSDEDX[i]; fNCDEDX[i] = s.fNCDEDX[i]; + fNCDEDXInclThres[i] = s.fNCDEDXInclThres[i]; } + fDEDX[4] = s.fDEDX[4]; for (Int_t i=0;i<12;i++) fOverlapLabels[i] = s.fOverlapLabels[i]; + } @@ -122,13 +147,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), @@ -142,7 +169,8 @@ AliTPCseed::AliTPCseed(const AliTPCtrack &t): fSeed1(-1), fSeed2(-1), fMAngular(0), - fCircular(0) + fCircular(0), + fPoolID(-1) { // // Constructor from AliTPCtrack @@ -163,7 +191,9 @@ AliTPCseed::AliTPCseed(const AliTPCtrack &t): fDEDX[i] = 0.; fSDEDX[i] = 1e10; fNCDEDX[i] = 0; + fNCDEDXInclThres[i] = 0; } + fDEDX[4] = 0; for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1; } @@ -172,13 +202,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), @@ -192,7 +224,8 @@ AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5], fSeed1(-1), fSeed2(-1), fMAngular(0), - fCircular(0) + fCircular(0), + fPoolID(-1) { // // Constructor @@ -205,58 +238,84 @@ AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5], fDEDX[i] = 0.; fSDEDX[i] = 1e10; fNCDEDX[i] = 0; + fNCDEDXInclThres[i] = 0; } + fDEDX[4] = 0; for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1; } 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]; } } -} -AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i) -{ - // - // - return &fTrackPoints[i]; } - -void AliTPCseed::RebuildSeed() +//_________________________________________________ +AliTPCseed & AliTPCseed::operator=(const AliTPCseed ¶m) { // - // rebuild seed to be ready for storing - AliTPCclusterMI cldummy; - cldummy.SetQ(0); - AliTPCTrackPoint pdummy; - pdummy.GetTPoint().SetShared(10); - for (Int_t i=0;i<160;i++){ - AliTPCclusterMI * cl0 = fClusterPointer[i]; - AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(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 + // don't touch pool ID + // + if(this!=¶m){ + AliTPCtrack::operator=(param); + fEsd =param.fEsd; + fClusterOwner = param.fClusterOwner; + if (!fClusterOwner) for(Int_t i = 0;i<160;++i)fClusterPointer[i] = param.fClusterPointer[i]; + else for(Int_t i = 0;i<160;++i) { + delete fClusterPointer[i]; + fClusterPointer[i] = new AliTPCclusterMI(*(param.fClusterPointer[i])); + } + // 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]; + fNCDEDXInclThres[i] = param.fNCDEDXInclThres[i]; } + fDEDX[4] = param.fDEDX[4]; + for(Int_t i = 0;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); } @@ -447,11 +517,11 @@ Int_t AliTPCseed::Compare(const TObject *o) const { } else { Float_t f2 =1; - f2 = 1-20*TMath::Sqrt(t->GetSigma1Pt2())/(TMath::Abs(t->Get1Pt())+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(GetSigma1Pt2())/(TMath::Abs(Get1Pt())+0.0066); + f1 = 1-20*TMath::Sqrt(GetSigma1Pt2())/(OneOverPt()+0.0066); if (fBConstrain) f1=1.2; @@ -464,267 +534,320 @@ 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. //----------------------------------------------------------------- - - 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.; + // CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal) + AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); - Float_t mean[4] = {0,0,0,0}; - Float_t sigma[4] = {1000,1000,1000,1000}; - Int_t nc[4] = {0,0,0,0}; - Float_t norm[4] = {1000,1000,1000,1000}; - // - // - 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; + 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); + fDEDX[4] = CookdEdxAnalytical(low,up,useTot ,row0,row2, 0); // full OROC truncated mean + // + 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(GetTPCClustInfo(2, 1, i1 , i2)); + fNCDEDX[1] = TMath::Nint(GetTPCClustInfo(2, 1, 0 , row0)); + fNCDEDX[2] = TMath::Nint(GetTPCClustInfo(2, 1, row0, row1)); + fNCDEDX[3] = TMath::Nint(GetTPCClustInfo(2, 1, row1, row2)); + // + fNCDEDXInclThres[0] = TMath::Nint(GetTPCClustInfo(2, 2, i1 , i2)); + fNCDEDXInclThres[1] = TMath::Nint(GetTPCClustInfo(2, 2, 0 , row0)); + fNCDEDXInclThres[2] = TMath::Nint(GetTPCClustInfo(2, 2, row0, row1)); + fNCDEDXInclThres[3] = TMath::Nint(GetTPCClustInfo(2, 2, row1, row2)); + // + SetdEdx(fDEDX[0]); + return fDEDX[0]; - } - 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); - } +// 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; + +// 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() @@ -742,7 +865,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) { @@ -763,187 +886,845 @@ 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;; - - // - 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 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 + } - 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}; + const Float_t ktany = TMath::Tan(TMath::DegToRad()*10); + const Float_t kedgey =3.; // // - fNShared =0; + 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; + } - // 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);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];iGetTransform(); + const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam(); + if (trans && recoParam->GetUseGainCorrectionTime()>0) { + runNumber = trans->GetCurrentRunNumber(); + //AliTPCcalibDB::Instance()->SetRun(runNumber); + TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber); + if (timeGainSplines) { + UInt_t time = trans->GetCurrentTimeStamp(); + AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0); + AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1); + if (fitMIP) { + corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/ + } else { + if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/ + } + } + } + mean /= corrTimeGain; + rms /= corrTimeGain; + // + if (returnVal==1) return rms; + if (returnVal==2) return ncl; + return mean; +} + +Float_t AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal, Int_t rowThres, Int_t mode){ + + // + // 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 + // - 3 ratio + // - 4 mean upper half + // - 5 mean - lower half + // - 6 third moment + // mode - 0 - linear + // - 1 - logatithmic + // rowThres - number of rows before and after given pad row to check for clusters below threshold + // + // 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; + Int_t nclBelowThr = 0; // counts number of clusters below threshold + // + // + Float_t gainGG = 1; // gas gain factor -always enabled + Float_t gainPad = 1; // gain map - used always + Float_t corrPos = 1; // local position correction - if posNorm enabled + // + // + // + if (AliTPCcalibDB::Instance()->GetParameters()){ + gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000; //relative gas gain + } + // + // extract time-dependent correction for pressure and temperature variations + // + UInt_t runNumber = 1; + Float_t corrTimeGain = 1; + TObjArray * timeGainSplines = 0x0; + TGraphErrors * grPadEqual = 0x0; + TGraphErrors* grChamberGain[3]={0x0,0x0,0x0}; + // + AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform(); + const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam(); + // + if (recoParam->GetNeighborRowsDedx() == 0) rowThres = 0; + // + if (trans) { + runNumber = trans->GetCurrentRunNumber(); + //AliTPCcalibDB::Instance()->SetRun(runNumber); + timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber); + if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) { + UInt_t time = trans->GetCurrentTimeStamp(); + AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0); + AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1); + if (fitMIP) { + corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time); /*fitMIP->Eval(time);*/ + } else { + if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time); /*fitFPcosmic->Eval(time); */ + } // - sumw += weight[i]; - sumamp += weight[i]*ampl; - sumamp2 += weight[i]*ampl*ampl; - norm[of] += angular[i]*weight[i]; - nc[of]++; + if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL"); + if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL"); + const char* names[3]={"SHORT","MEDIUM","LONG"}; + for (Int_t iPadRegion=0; iPadRegion<3; ++iPadRegion) + grChamberGain[iPadRegion]=(TGraphErrors*)timeGainSplines->FindObject(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[iPadRegion])); } - if (sumw<1){ - SetdEdx(0); - } - else { - norm[of] /= sumw; - mean[of] = sumamp/sumw; - sigma[of] = sumamp2/sumw-mean[of]*mean[of]; - if (sigma[of]>0.1) - sigma[of] = TMath::Sqrt(sigma[of]); - else - sigma[of] = 1000; - mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; - } } + + const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam + const Float_t ktany = TMath::Tan(TMath::DegToRad()*10); + const Float_t kedgey =3.; + // + // + for (Int_t irow=i1; irow 1 && irow < 157) { + Bool_t isClBefore = kFALSE; + Bool_t isClAfter = kFALSE; + for(Int_t ithres = 1; ithres <= rowThres; ithres++) { + AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres); + if (clusterBefore) isClBefore = kTRUE; + AliTPCclusterMI * clusterAfter = GetClusterPointer(irow + ithres); + if (clusterAfter) isClAfter = kTRUE; + } + if (isClBefore && isClAfter) nclBelowThr++; + } + if (!cluster) continue; + // + // + if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster + // + AliTPCTrackerPoint * point = GetTrackPoint(irow); + if (point==0) continue; + Float_t rsigmay = TMath::Sqrt(point->GetSigmaY()); + if (rsigmay > kClusterShapeCut) continue; + // + if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb + // + Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ(); + Int_t ipad= 0; + if (irow>=row0) ipad=1; + if (irow>=row1) ipad=2; + // + // + // + AliTPCCalPad * gainMap = AliTPCcalibDB::Instance()->GetDedxGainFactor(); + if (gainMap) { + // + // Get gainPad - pad by pad calibration + // + Float_t factor = 1; + AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector()); + if (irow < row0) { // IROC + factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad())); + } else { // OROC + factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad())); + } + if (factor>0.3) gainPad=factor; + } + // + // Do position normalization - relative distance to + // center of pad- time bin - Float_t 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]; + 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); + // } - fDEDX[i] = mean[i]; - fSDEDX[i] = sigma[i]; - fNCDEDX[i]= nc[i]; + if (type==0) { + corrPos = parcl->GetQnormCorr(ipad, type,5)* + parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), + cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,effLength,effDiff); + + Float_t drm = 0.5-TMath::Abs(cluster->GetZ()/250.); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty); + corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz); + // + } + // + // pad region equalization outside of cluster param + // + Float_t gainEqualPadRegion = 1; + if (grPadEqual && recoParam->GetUseGainCorrectionTime()>0) gainEqualPadRegion = grPadEqual->Eval(ipad); + // + // chamber-by-chamber equalization outside gain map + // + Float_t gainChamber = 1; + if (grChamberGain[ipad] && recoParam->GetUseGainCorrectionTime()>0) gainChamber = grChamberGain[ipad]->Eval(cluster->GetDetector()); + // + amp[ncl]=charge; + amp[ncl]/=gainGG; + amp[ncl]/=gainPad; + amp[ncl]/=corrPos; + amp[ncl]/=gainEqualPadRegion; + amp[ncl]/=gainChamber; + // + ncl++; } - if (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; + if (type==2) return ncl; + TMath::Sort(ncl,amp, indexes, kFALSE); + // + if (ncl<10) return 0; + // + Double_t * ampWithBelow = new Double_t[ncl + nclBelowThr]; + for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) { + if (iCl < nclBelowThr) { + ampWithBelow[iCl] = amp[indexes[0]]; + } else { + ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]]; + } } - else{ - SetdEdx(0); - return; + //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]); + // + Float_t suma=0; + Float_t suma2=0; + Float_t suma3=0; + Float_t sumaS=0; + Float_t sumn=0; + // upper,and lower part statistic + Float_t sumL=0, sumL2=0, sumLN=0; + Float_t sumD=0, sumD2=0, sumDN=0; + + Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low); + Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up); + Int_t iclm=TMath::Nint((ncl + nclBelowThr)*(low +(up+low)*0.5)); + // + for (Int_t icl=icl0; icliclm){ + sumL+=camp; + sumL2+=camp*camp; + sumLN++; + } + if (icl<=iclm){ + sumD+=camp; + sumD2+=camp*camp; + sumDN++; + } } - // Float_t dedx1 =dedx; + // + Float_t mean = 0; + Float_t meanL = 0; + Float_t meanD = 0; // lower half mean + if (sumn > 1e-30) mean =suma/sumn; + if (sumLN > 1e-30) meanL =sumL/sumLN; + if (sumDN > 1e-30) meanD =(sumD/sumDN); + /* + Float_t mean =suma/sumn; + Float_t meanL = sumL/sumLN; + Float_t meanD =(sumD/sumDN); // lower half mean + */ + + Float_t rms = 0; + Float_t mean2=0; + Float_t mean3=0; + Float_t meanS=0; + + if(sumn>0){ + rms = TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean)); + mean2=suma2/sumn; + mean3=suma3/sumn; + meanS=sumaS/sumn; + } + + if (mean2>0) mean2=TMath::Power(TMath::Abs(mean2),1./2.); + if (mean3>0) mean3=TMath::Power(TMath::Abs(mean3),1./3.); + if (meanS>0) meanS=TMath::Power(TMath::Abs(meanS),3.); + // + if (mode==1) mean=TMath::Exp(mean); + if (mode==1) meanL=TMath::Exp(meanL); // upper truncation + if (mode==1) meanD=TMath::Exp(meanD); // lower truncation + // + delete [] ampWithBelow; - 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]); + + // + if (returnVal==1) return rms; + if (returnVal==2) return ncl; + if (returnVal==3) return Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl); + if (returnVal==4) return meanL; + if (returnVal==5) return meanD; + if (returnVal==6) return mean2; + if (returnVal==7) return mean3; + if (returnVal==8) return meanS; + return mean; +} + + + + +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; + + 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 || (track->GetClusterIndex(i) & 0x8000)) 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 || (track->GetClusterIndex(i) & 0x8000)) continue; + // if (RejectCluster(c,track)) continue; + sector = (c->GetDetector()%18); + if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) { + //continue; } - fDEDX[i] = mean[i]; + 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 (norm4>0) dedx /= norm4; - + if (!isOK) { delete track; return 0;} + track->AddCovariance(covar); + // + // + // + for (Int_t i=imax; i>=imin; i--){ + AliTPCclusterMI *c=track->GetClusterPointer(i); + if (!c || (track->GetClusterIndex(i) & 0x8000)) 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 || (track->GetClusterIndex(i) & 0x8000)) continue; + sector = (c->GetDetector()%18); + if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) { + //continue; + } + ncl++; + //if (RejectCluster(c,track)) continue; + Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()}; + Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation + if (!track->PropagateTo(r[0])) { + isOK=kFALSE; + } + if ( !((static_cast(track)->Update(&r[1],cov)))) isOK=kFALSE; + } + //if (!isOK) { delete track; return 0;} + paramOut=*track; + // + // + // + if (parin) (*parin)=paramIn; + if (parout) (*parout)=paramOut; + delete track; + return ncl; +} + + +Bool_t AliTPCseed::RefitTrack(AliTPCseed* /*seed*/, Bool_t /*out*/){ + // + // + // + return kFALSE; +} + + + + + + +void AliTPCseed::GetError(AliTPCclusterMI* cluster, AliExternalTrackParam * param, + Double_t& erry, Double_t &errz) +{ + // + // Get cluster error at given position + // + AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam(); + Double_t tany,tanz; + Double_t snp1=param->GetSnp(); + tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Double_t tgl1=param->GetTgl(); + tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Int_t padSize = 0; // short pads + if (cluster->GetDetector() >= 36) { + padSize = 1; // medium pads + if (cluster->GetRow() > 63) padSize = 2; // long pads + } + + erry = clusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany) ); + errz = clusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ); +} + + +void AliTPCseed::GetShape(AliTPCclusterMI* cluster, AliExternalTrackParam * param, + Double_t& rmsy, Double_t &rmsz) +{ + // + // Get cluster error at given position + // + AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam(); + Double_t tany,tanz; + Double_t snp1=param->GetSnp(); + tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Double_t tgl1=param->GetTgl(); + tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1)); + // + Int_t padSize = 0; // short pads + if (cluster->GetDetector() >= 36) { + padSize = 1; // medium pads + if (cluster->GetRow() > 63) padSize = 2; // long pads + } + + rmsy = clusterParam->GetRMSQ( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany), TMath::Abs(cluster->GetMax()) ); + rmsz = clusterParam->GetRMSQ( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ,TMath::Abs(cluster->GetMax())); +} + + + +Double_t AliTPCseed::GetQCorrGeom(Float_t ty, Float_t tz){ + //Geoetrical + //ty - tangent in local y direction + //tz - tangent + // + Float_t norm=TMath::Sqrt(1+ty*ty+tz*tz); + return norm; +} + +Double_t AliTPCseed::GetQCorrShape(Int_t ipad, Int_t type,Float_t z, Float_t ty, Float_t tz, Float_t /*q*/, Float_t /*thr*/){ + // + // Q normalization + // + // return value = Q Normalization factor + // Normalization - 1 - shape factor part for full drift + // 1 - electron attachment for 0 drift + + // Input parameters: + // + // ipad - 0 short pad + // 1 medium pad + // 2 long pad + // + // type - 0 qmax + // - 1 qtot + // + //z - z position (-250,250 cm) + //ty - tangent in local y direction + //tz - tangent + // + + AliTPCClusterParam * paramCl = AliTPCcalibDB::Instance()->GetClusterParam(); + AliTPCParam * paramTPC = AliTPCcalibDB::Instance()->GetParameters(); + + if (!paramCl) return 1; + // + Double_t dr = 250.-TMath::Abs(z); + Double_t sy = paramCl->GetRMS0( 0,ipad, dr, TMath::Abs(ty)); + Double_t sy0= paramCl->GetRMS0(0,ipad, 250, 0); + Double_t sz = paramCl->GetRMS0( 1,ipad, dr, TMath::Abs(tz)); + Double_t sz0= paramCl->GetRMS0(1,ipad, 250, 0); + + Double_t sfactorMax = TMath::Sqrt(sy0*sz0/(sy*sz)); + + + Double_t dt = 1000000*(dr/paramTPC->GetDriftV()); //time in microsecond + Double_t attProb = TMath::Exp(-paramTPC->GetAttCoef()*paramTPC->GetOxyCont()*dt); + // + // + if (type==0) return sfactorMax*attProb; + + return attProb; + + +} + + +//_______________________________________________________________________ +Float_t AliTPCseed::GetTPCClustInfo(Int_t nNeighbours, Int_t type, Int_t row0, Int_t row1) +{ + // + // TPC cluster information + // type 0: get fraction of found/findable clusters with neighbourhood definition + // 1: found clusters + // 2: findable (number of clusters above and below threshold) + // + // definition of findable clusters: + // a cluster is defined as findable if there is another cluster + // within +- nNeighbours pad rows. The idea is to overcome threshold + // effects with a very simple algorithm. + // + + const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam + const Float_t ktany = TMath::Tan(TMath::DegToRad()*10); + const Float_t kedgey =3.; - SetdEdx(dedx); - - //mi deDX + Float_t ncl = 0; + Float_t nclBelowThr = 0; // counts number of clusters below threshold + + for (Int_t irow=row0; irow 1 && irow < 157) { + Bool_t isClBefore = kFALSE; + Bool_t isClAfter = kFALSE; + for(Int_t ithres = 1; ithres <= nNeighbours; ithres++) { + AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres); + if (clusterBefore) isClBefore = kTRUE; + AliTPCclusterMI * clusterAfter = GetClusterPointer(irow + ithres); + if (clusterAfter) isClAfter = kTRUE; + } + if (isClBefore && isClAfter) nclBelowThr++; + } + if (!cluster) continue; + // + // + if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster + // + AliTPCTrackerPoint * point = GetTrackPoint(irow); + if (point==0) continue; + Float_t rsigmay = TMath::Sqrt(point->GetSigmaY()); + if (rsigmay > kClusterShapeCut) continue; + // + if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb + ncl++; + } + if(ncl<10) + return 0; + if(type==0) + if(nclBelowThr+ncl>0) + return ncl/(nclBelowThr+ncl); + if(type==1) + return ncl; + if(type==2) + return ncl+nclBelowThr; + return 0; } -*/ -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); +//_______________________________________________________________________ +Int_t AliTPCseed::GetNumberOfClustersIndices() { + Int_t ncls = 0; + for (int i=0; i < 160; i++) { + if ((fIndex[i] & 0x8000) == 0) + ncls++; + } + return ncls; } +//_______________________________________________________________________ +void AliTPCseed::Clear(Option_t*) +{ + // formally seed may allocate memory for clusters (althought this should not happen for + // the seeds in the pool). Hence we need this method for fwd. compatibility + if (fClusterOwner) for (int i=160;i--;) {delete fClusterPointer[i]; fClusterPointer[i] = 0;} +}