X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCtrackerMI.cxx;h=ce187049a08e3fdf6f045e6513002858e976640b;hb=4f543b2feaf125022621d2f7dc16a14fc3a48350;hp=ab30679a7e7d700011560330a8669ec6510052ba;hpb=6bdc18d6e05aba04ff7afd0bab3452f33c9e687b;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCtrackerMI.cxx b/TPC/AliTPCtrackerMI.cxx index ab30679a7e7..ce187049a08 100644 --- a/TPC/AliTPCtrackerMI.cxx +++ b/TPC/AliTPCtrackerMI.cxx @@ -19,42 +19,42 @@ // // Origin: Marian Ivanov Marian.Ivanov@cern.ch // -// AliTPC parallel tracker - -// How to use? - -// run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded -// run AliTPCFindTracksMI.C macro - to find tracks -// tracks are written to AliTPCtracks.root file -// for comparison also seeds are written to the same file - to special branch +// AliTPC parallel tracker //------------------------------------------------------- /* $Id$ */ - -#include +#include "Riostream.h" +#include #include +#include #include -#include - -#include "Riostream.h" +#include "AliLog.h" -#include "AliTPCclusterMI.h" #include "AliComplexCluster.h" -#include "AliTPCParam.h" -#include "AliTPCClustersRow.h" -#include "AliTPCpolyTrack.h" -#include "TStopwatch.h" #include "AliESD.h" +#include "AliESDkink.h" #include "AliHelix.h" -// #include "AliRunLoader.h" -// -#include "AliTPCreco.h" +#include "AliTPCClustersRow.h" +#include "AliTPCParam.h" +#include "AliTPCReconstructor.h" +#include "AliTPCclusterMI.h" +#include "AliTPCpolyTrack.h" +#include "AliTPCreco.h" +#include "AliTPCseed.h" #include "AliTPCtrackerMI.h" +#include "TStopwatch.h" +#include "AliTPCReconstructor.h" +#include "AliESDkink.h" +#include "AliPID.h" +#include "TTreeStream.h" +#include "AliAlignObj.h" +#include "AliTrackPointArray.h" +// - -ClassImp(AliTPCseed) ClassImp(AliTPCtrackerMI) @@ -67,7 +67,7 @@ public: }; Double_t AliTPCFastMath::fgFastAsin[20000]; -AliTPCFastMath gAliTPCFastMath; +AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT AliTPCFastMath::AliTPCFastMath(){ // @@ -124,13 +124,13 @@ Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){ track->fClusterPointer[track->fRow] = c; // - Float_t angle2 = track->GetSnp()*track->GetSnp(); - angle2 = TMath::Sqrt(angle2/(1-angle2)); + Double_t angle2 = track->GetSnp()*track->GetSnp(); // //SET NEW Track Point // - // if (debug) + if (angle2<1) //PH sometimes angle2 is very big. To be investigated... { + angle2 = TMath::Sqrt(angle2/(1-angle2)); AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->fRow)); // point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2); @@ -257,8 +257,24 @@ AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2) fNewIO =0; fDebug =0; fEvent =0; + fDebugStreamer = new TTreeSRedirector("TPCdebug.root"); +} +//________________________________________________________________________ +AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): + AliTracker(t), + fkNIS(t.fkNIS), + fkNOS(t.fkNOS) +{ + //------------------------------------ + // dummy copy constructor + //------------------------------------------------------------------ +} +AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){ + //------------------------------ + // dummy + //-------------------------------------------------------------- + return *this; } - //_____________________________________________________________________________ AliTPCtrackerMI::~AliTPCtrackerMI() { //------------------------------------------------------------------ @@ -270,15 +286,16 @@ AliTPCtrackerMI::~AliTPCtrackerMI() { fSeeds->Delete(); delete fSeeds; } + if (fDebugStreamer) delete fDebugStreamer; } void AliTPCtrackerMI::SetIO() { // fNewIO = kTRUE; - fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::fgkDefaultEventFolderName); + fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName()); - fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::fgkDefaultEventFolderName); + fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName()); if (fOutput){ AliTPCtrack *iotrack= new AliTPCtrack; fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100); @@ -344,18 +361,113 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) // write tracks to the event // store index of the track Int_t nseed=arr->GetEntriesFast(); + //FindKinks(arr,fEvent); for (Int_t i=0; iUncheckedAt(i); if (!pt) continue; - pt->PropagateTo(fParam->GetInnerRadiusLow()); - if (pt->GetNumberOfClusters()>70) { + pt->UpdatePoints(); + // pt->PropagateTo(fParam->GetInnerRadiusLow()); + if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks + pt->PropagateTo(fParam->GetInnerRadiusLow()); + } + + if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){ + AliESDtrack iotrack; + iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); + iotrack.SetTPCPoints(pt->GetPoints()); + iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + // iotrack.SetTPCpid(pt->fTPCr); + //iotrack.SetTPCindex(i); + fEvent->AddTrack(&iotrack); + continue; + } + + if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) { + AliESDtrack iotrack; + iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); + iotrack.SetTPCPoints(pt->GetPoints()); + //iotrack.SetTPCindex(i); + iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + // iotrack.SetTPCpid(pt->fTPCr); + fEvent->AddTrack(&iotrack); + continue; + } + // + // short tracks - maybe decays + + if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.70) { + Int_t found,foundable,shared; + pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE); + if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2)){ + AliESDtrack iotrack; + iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); + //iotrack.SetTPCindex(i); + iotrack.SetTPCPoints(pt->GetPoints()); + iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + //iotrack.SetTPCpid(pt->fTPCr); + fEvent->AddTrack(&iotrack); + continue; + } + } + + if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.8) { + Int_t found,foundable,shared; + pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE); + if (found<20) continue; + if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue; + // + AliESDtrack iotrack; + iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); + iotrack.SetTPCPoints(pt->GetPoints()); + iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + //iotrack.SetTPCpid(pt->fTPCr); + //iotrack.SetTPCindex(i); + fEvent->AddTrack(&iotrack); + continue; + } + // short tracks - secondaties + // + if ( (pt->GetNumberOfClusters()>30) ) { + Int_t found,foundable,shared; + pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE); + if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){ + AliESDtrack iotrack; + iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); + iotrack.SetTPCPoints(pt->GetPoints()); + iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + //iotrack.SetTPCpid(pt->fTPCr); + //iotrack.SetTPCindex(i); + fEvent->AddTrack(&iotrack); + continue; + } + } + + if ( (pt->GetNumberOfClusters()>15)) { + Int_t found,foundable,shared; + pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE); + if (found<15) continue; + if (foundable<=0) continue; + if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue; + if (float(found)/float(foundable)<0.8) continue; + // AliESDtrack iotrack; iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); + iotrack.SetTPCPoints(pt->GetPoints()); + iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + // iotrack.SetTPCpid(pt->fTPCr); //iotrack.SetTPCindex(i); fEvent->AddTrack(&iotrack); - } + continue; + } } } + printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()); } void AliTPCtrackerMI::WriteTracks(TTree * tree) @@ -530,7 +642,7 @@ Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ())); Int_t ctype = cl->GetType(); Float_t padlength= GetPadPitchLength(seed->fRow); - Float_t angle2 = seed->GetSnp()*seed->GetSnp(); + Double_t angle2 = seed->GetSnp()*seed->GetSnp(); angle2 = angle2/(1-angle2); // //cluster "quality" @@ -665,7 +777,7 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ Int_t ctype = cl->GetType(); Float_t padlength= GetPadPitchLength(seed->fRow); // - Float_t angle2 = seed->GetSnp()*seed->GetSnp(); + Double_t angle2 = seed->GetSnp()*seed->GetSnp(); // if (angle2<0.6) angle2 = 0.6; angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2)); // @@ -800,160 +912,6 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ -void AliTPCseed::Reset(Bool_t all) -{ - // - // - SetNumberOfClusters(0); - fNFoundable = 0; - SetChi2(0); - ResetCovariance(); - /* - if (fTrackPoints){ - for (Int_t i=0;i<8;i++){ - delete [] fTrackPoints[i]; - } - delete fTrackPoints; - fTrackPoints =0; - } - */ - - if (all){ - for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); - for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; - } - -} - - -void AliTPCseed::Modify(Double_t factor) -{ - - //------------------------------------------------------------------ - //This function makes a track forget its history :) - //------------------------------------------------------------------ - if (factor<=0) { - ResetCovariance(); - return; - } - fC00*=factor; - fC10*=0; fC11*=factor; - fC20*=0; fC21*=0; fC22*=factor; - fC30*=0; fC31*=0; fC32*=0; fC33*=factor; - fC40*=0; fC41*=0; fC42*=0; fC43*=0; fC44*=factor; - SetNumberOfClusters(0); - fNFoundable =0; - SetChi2(0); - fRemoval = 0; - fCurrentSigmaY2 = 0.000005; - fCurrentSigmaZ2 = 0.000005; - fNoCluster = 0; - //fFirstPoint = 160; - //fLastPoint = 0; -} - - - - -Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const -{ - //----------------------------------------------------------------- - // This function find proloncation of a track to a reference plane x=xk. - // doesn't change internal state of the track - //----------------------------------------------------------------- - - Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1; - - if (TMath::Abs(fP4*xk - fP2) >= 0.999) { - return 0; - } - - // Double_t y1=fP0, z1=fP1; - Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1); - Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2); - - y = fP0; - z = fP1; - //y += dx*(c1+c2)/(r1+r2); - //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3; - - Double_t dy = dx*(c1+c2)/(r1+r2); - Double_t dz = 0; - // - Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1); - /* - if (TMath::Abs(delta)>0.0001){ - dz = fP3*TMath::ASin(delta)/fP4; - }else{ - dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1); - } - */ - dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4; - // - y+=dy; - z+=dz; - - - return 1; -} - - -//_____________________________________________________________________________ -Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const -{ - //----------------------------------------------------------------- - // This function calculates a predicted chi2 increment. - //----------------------------------------------------------------- - //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2(); - Double_t r00=fErrorY2, r01=0., r11=fErrorZ2; - r00+=fC00; r01+=fC10; r11+=fC11; - - Double_t det=r00*r11 - r01*r01; - if (TMath::Abs(det) < 1.e-10) { - Int_t n=GetNumberOfClusters(); - if (n>4) cerr<GetY() - fP0, dz=c->GetZ() - fP1; - - return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; -} - - -//_________________________________________________________________________________________ - - -Int_t AliTPCseed::Compare(const TObject *o) const { - //----------------------------------------------------------------- - // This function compares tracks according to the sector - for given sector according z - //----------------------------------------------------------------- - AliTPCseed *t=(AliTPCseed*)o; - - if (fSort == 0){ - if (t->fRelativeSector>fRelativeSector) return -1; - if (t->fRelativeSectorGetZ(); - Double_t z1 = GetZ(); - if (z2>z1) return 1; - if (z2fC44)/(TMath::Abs(t->GetC())+0.0066); - if (t->fBConstrain) f2=1.2; - - Float_t f1 =1; - f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066); - - if (fBConstrain) f1=1.2; - - if (t->GetNumberOfClusters()*f2 GetY() - fP0, dz=c->GetZ() - fP1; - Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz; - if (TMath::Abs(cur*fX-eta) >= 0.9) { - return 0; - } - - fP0 += k00*dy + k01*dz; - fP1 += k10*dy + k11*dz; - fP2 = eta; - fP3 += k30*dy + k31*dz; - fP4 = cur; - - Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40; - Double_t c12=fC21, c13=fC31, c14=fC41; - - fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11; - fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13; - fC40-=k00*c04+k01*c14; - - fC11-=k10*c01+k11*fC11; - fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13; - fC41-=k10*c04+k11*c14; - - fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13; - fC42-=k20*c04+k21*c14; - - fC33-=k30*c03+k31*c13; - fC43-=k40*c03+k41*c13; - - fC44-=k40*c04+k41*c14; - - Int_t n=GetNumberOfClusters(); - // fIndex[n]=index; - SetNumberOfClusters(n+1); - SetChi2(GetChi2()+chisq); - - return 1; -} - - - //_____________________________________________________________________________ Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1, Double_t x2,Double_t y2, @@ -1229,6 +1128,9 @@ Int_t AliTPCtrackerMI::LoadClusters() // Int_t sec,row; fParam->AdjustSectorRow(clrow->GetID(),sec,row); + for (Int_t icl=0; iclGetArray()->GetEntriesFast(); icl++){ + Transform((AliCluster*)(clrow->GetArray()->At(icl))); + } // AliTPCRow * tpcrow=0; Int_t left=0; @@ -1291,6 +1193,21 @@ void AliTPCtrackerMI::UnloadClusters() return ; } +void AliTPCtrackerMI::Transform(AliCluster * cluster){ + // + // + // + //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices(); + TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector()); + //TGeoHMatrix mat; + Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()}; + Double_t posC[3]; + //mat.LocalToMaster(pos,posC); + mat->LocalToMaster(pos,posC); + cluster->SetX(posC[0]); + cluster->SetY(posC[1]); + cluster->SetZ(posC[2]); +} //_____________________________________________________________________________ Int_t AliTPCtrackerMI::LoadOuterSectors() { @@ -1396,29 +1313,46 @@ AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const { //-------------------------------------------------------------------- // Return pointer to a given cluster //-------------------------------------------------------------------- + if (index<0) return 0; // no cluster Int_t sec=(index&0xff000000)>>24; Int_t row=(index&0x00ff0000)>>16; Int_t ncl=(index&0x00007fff)>>00; const AliTPCRow * tpcrow=0; AliTPCclusterMI * clrow =0; + + if (sec<0 || sec>=fkNIS*4) { + AliWarning(Form("Wrong sector %d",sec)); + return 0x0; + } + if (secfN1<=ncl) return 0; clrow = tpcrow->fClusters1; - else + } + else { + if (tpcrow->fN2<=ncl) return 0; clrow = tpcrow->fClusters2; + } } - else{ + else { tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]); - if (sec-2*fkNISfN1<=ncl) return 0; clrow = tpcrow->fClusters1; - else + } + else { + if (tpcrow->fN2<=ncl) return 0; clrow = tpcrow->fClusters2; + } } - if (tpcrow==0) return 0; - if (tpcrow->GetN()<=ncl) return 0; - // return (AliTPCclusterMI*)(*tpcrow)[ncl]; + return &(clrow[ncl]); } @@ -1431,16 +1365,64 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { //----------------------------------------------------------------- // Double_t x= GetXrow(nr), ymax=GetMaxY(nr); - - // if (t.GetRadius()>x+10 ) return 0; - // t.PropagateTo(x+0.02); - //t.PropagateTo(x+0.01); - if (!t.PropagateTo(x)) { - t.fRemoval = 10; + AliTPCclusterMI *cl=0; + Int_t tpcindex= t.GetClusterIndex2(nr); + // + // update current shape info every 5 pad-row + // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){ + GetShape(&t,nr); + //} + // + if (fIteration>0 && tpcindex>=-1){ //if we have already clusters + // + if (tpcindex==-1) return 0; //track in dead zone + if (tpcindex>0){ // + cl = t.fClusterPointer[nr]; + if ( (cl==0) ) cl = GetClusterMI(tpcindex); + t.fCurrentClusterIndex1 = tpcindex; + } + if (cl){ + Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector + Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift(); + // + if (angle<-TMath::Pi()) angle += 2*TMath::Pi(); + if (angle>=TMath::Pi()) angle -= 2*TMath::Pi(); + + if (TMath::Abs(angle-t.GetAlpha())>0.001){ + Double_t rotation = angle-t.GetAlpha(); + t.fRelativeSector= relativesector; + if (!t.Rotate(rotation)) return 0; + } + if (!t.PropagateTo(x)) return 0; + // + t.fCurrentCluster = cl; + t.fRow = nr; + Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.); + if ((tpcindex&0x8000)==0) accept =0; + if (accept<3) { + //if founded cluster is acceptible + if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty + t.fErrorY2 += 0.03; + t.fErrorZ2 += 0.03; + t.fErrorY2 *= 3; + t.fErrorZ2 *= 3; + } + t.fNFoundable++; + UpdateTrack(&t,accept); + return 1; + } + } + } + if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle + if (fIteration>1){ + // not look for new cluster during refitting + t.fNFoundable++; return 0; } // - Double_t y=t.GetY(), z=t.GetZ(); + UInt_t index=0; + // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06 + Double_t y=t.GetYat(x); if (TMath::Abs(y)>ymax){ if (y > ymax) { t.fRelativeSector= (t.fRelativeSector+1) % fN; @@ -1451,33 +1433,26 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { if (!t.Rotate(-fSectors->GetAlpha())) return 0; } - //if (!t.PropagateTo(x)) { - // return 0; - //} - return 1; + //return 1; } // - // update current shape info every 3 pad-row - if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){ - //t.fCurrentSigmaY = GetSigmaY(&t); - //t.fCurrentSigmaZ = GetSigmaZ(&t); - GetShape(&t,nr); + if (!t.PropagateTo(x)) { + if (fIteration==0) t.fRemoval = 10; + return 0; } - // - AliTPCclusterMI *cl=0; - UInt_t index=0; - - - //Int_t nr2 = nr; - if (t.GetClusterIndex2(nr)>0){ - // - //cl = GetClusterMI(t.GetClusterIndex2(nr)); - index = t.GetClusterIndex2(nr); - cl = t.fClusterPointer[nr]; - if ( (cl==0) && (index>0)) cl = GetClusterMI(index); - t.fCurrentClusterIndex1 = index; + y=t.GetY(); + Double_t z=t.GetZ(); + // + if (!IsActive(t.fRelativeSector,nr)) { + t.fInDead = kTRUE; + t.SetClusterIndex2(nr,-1); + return 0; } + //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha())); + Bool_t isActive = IsActive(t.fRelativeSector,nr); + Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.fRelativeSector][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.fRelativeSector][nr].GetN()>0; + if (!isActive || !isActive2) return 0; const AliTPCRow &krow=GetRow(t.fRelativeSector,nr); if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0; Double_t roady =1.; @@ -1490,34 +1465,28 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { } else { - if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++; + if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)GetZLength() && (TMath::Abs(t.GetSnp())0) accept =0; - if (accept<3) { - //if founded cluster is acceptible - UpdateTrack(&t,accept); - return 1; - } - } - if (krow) { // cl = krow.FindNearest2(y+10.,z,roady,roadz,index); cl = krow.FindNearest2(y,z,roady,roadz,index); if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index); } - // t.fNoCluster++; - if (cl) { t.fCurrentCluster = cl; t.fRow = nr; + if (fIteration==2&&cl->IsUsed(10)) return 0; Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.); + if (fIteration==2&&cl->IsUsed(11)) { + t.fErrorY2 += 0.03; + t.fErrorZ2 += 0.03; + t.fErrorY2 *= 3; + t.fErrorZ2 *= 3; + } /* if (t.fCurrentCluster->IsUsed(10)){ // @@ -1594,7 +1563,7 @@ Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) { } else { - if (TMath::Abs(z)>(1.05*x+10)) t.SetClusterIndex2(row,-1); + if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1); } //calculate @@ -1618,6 +1587,53 @@ Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) { +//_________________________________________________________________________ +Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const +{ + // Get track space point by index + // return false in case the cluster doesn't exist + AliTPCclusterMI *cl = GetClusterMI(index); + if (!cl) return kFALSE; + Int_t sector = (index&0xff000000)>>24; + // Int_t row = (index&0x00ff0000)>>16; + Float_t xyz[3]; + // xyz[0] = fParam->GetPadRowRadii(sector,row); + xyz[0] = cl->GetX(); + xyz[1] = cl->GetY(); + xyz[2] = cl->GetZ(); + Float_t sin,cos; + fParam->AdjustCosSin(sector,cos,sin); + Float_t x = cos*xyz[0]-sin*xyz[1]; + Float_t y = cos*xyz[1]+sin*xyz[0]; + Float_t cov[6]; + Float_t sigmaY2 = 0.027*cl->GetSigmaY2(); + if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07; + Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2(); + if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77; + cov[0] = sin*sin*sigmaY2; + cov[1] = -sin*cos*sigmaY2; + cov[2] = 0.; + cov[3] = cos*cos*sigmaY2; + cov[4] = 0.; + cov[5] = sigmaZ2; + p.SetXYZ(x,y,xyz[2],cov); + AliAlignObj::ELayerID iLayer; + Int_t idet; + if (sector < fParam->GetNInnerSector()) { + iLayer = AliAlignObj::kTPC1; + idet = sector; + } + else { + iLayer = AliAlignObj::kTPC2; + idet = sector - fParam->GetNInnerSector(); + } + UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,idet); + p.SetVolumeID(volid); + return kTRUE; +} + + + Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) { //----------------------------------------------------------------- // This function tries to find a track prolongation to next pad row @@ -1630,14 +1646,14 @@ Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) { Double_t ymax= GetMaxY(nr); if (row < nr) return 1; // don't prolongate if not information until now - - if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) { - t.fRemoval =10; - return 0; // not prolongate strongly inclined tracks - } - if (TMath::Abs(t.GetSnp())>0.95) { - t.fRemoval =10; - return 0; // not prolongate strongly inclined tracks - } +// if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) { +// t.fRemoval =10; +// return 0; // not prolongate strongly inclined tracks +// } +// if (TMath::Abs(t.GetSnp())>0.95) { +// t.fRemoval =10; +// return 0; // not prolongate strongly inclined tracks +// }// patch 28 fev 06 Double_t x= GetXrow(nr); Double_t y,z; @@ -1667,7 +1683,13 @@ Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) { //y = t.GetY(); } // - + if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; + if (!IsActive(t.fRelativeSector,nr)) { + t.fInDead = kTRUE; + t.SetClusterIndex2(nr,-1); + return 0; + } + //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha())); AliTPCRow &krow=GetRow(t.fRelativeSector,nr); if (TMath::Abs(TMath::Abs(y)-ymax)GetAlpha()+0.0001)%fN; Int_t first = GetRowNumber(xt)-1; - for (Int_t nr= first; nr>=rf; nr-=step) { + for (Int_t nr= first; nr>=rf; nr-=step) { + // update kink info + if (t.GetKinkIndexes()[0]>0){ + for (Int_t i=0;i<3;i++){ + Int_t index = t.GetKinkIndexes()[i]; + if (index==0) break; + if (index<0) continue; + // + AliESDkink * kink = fEvent->GetKink(index-1); + if (!kink){ + printf("PROBLEM\n"); + } + else{ + Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2))); + if (kinkrow==nr){ + AliExternalTrackParam paramd(t); + kink->SetDaughter(paramd); + kink->SetStatus(2,5); + kink->Update(); + } + } + } + } + + if (nr==80) t.UpdateReference(); if (nrGetNRows()) fSectors = fInnerSec; else fSectors = fOuterSec; if (FollowToNext(t,nr)==0) - if (!t.IsActive()) return 0; + if (!t.IsActive()) + return 0; } return 1; @@ -1815,19 +1865,41 @@ Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) { //----------------------------------------------------------------- // This function tries to find a track prolongation. //----------------------------------------------------------------- - // Double_t xt=t.GetX(); // + Double_t xt=t.GetX(); Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift(); if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); if (alpha < 0. ) alpha += 2.*TMath::Pi(); t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN; - Int_t first = 0; - first = t.fFirstPoint+3; + Int_t first = t.fFirstPoint; + if (first0.95)) break;//patch 28 fev 06 + if (t.GetKinkIndexes()[0]<0){ + for (Int_t i=0;i<3;i++){ + Int_t index = t.GetKinkIndexes()[i]; + if (index==0) break; + if (index>0) continue; + index = TMath::Abs(index); + AliESDkink * kink = fEvent->GetKink(index-1); + if (!kink){ + printf("PROBLEM\n"); + } + else{ + Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2))); + if (kinkrow==nr){ + AliExternalTrackParam paramm(t); + kink->SetMother(paramm); + kink->SetStatus(2,1); + kink->Update(); + } + } + } + } + // if (nrGetNRows()) fSectors = fInnerSec; else @@ -2157,7 +2229,96 @@ void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1, Float_t fact } } -void AliTPCtrackerMI::UnsignClusters() + +void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal) +{ + + //Loop over all tracks and remove "overlaps" + // + // + UnsignClusters(); + // + Int_t nseed = arr->GetEntriesFast(); + Float_t * quality = new Float_t[nseed]; + Int_t * indexes = new Int_t[nseed]; + Int_t good =0; + // + // + for (Int_t i=0; iUncheckedAt(i); + if (!pt){ + quality[i]=-1; + continue; + } + pt->UpdatePoints(); //select first last max dens points + Float_t * points = pt->GetPoints(); + if (points[3]<0.8) quality[i] =-1; + // + quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters(); + } + TMath::Sort(nseed,quality,indexes); + // + // + for (Int_t itrack=0; itrackUncheckedAt(trackindex); + if (quality[trackindex]<0){ + if (pt) { + delete arr->RemoveAt(trackindex); + } + else{ + arr->RemoveAt(trackindex); + } + continue; + } + // + Int_t first = Int_t(pt->GetPoints()[0]); + Int_t last = Int_t(pt->GetPoints()[2]); + Double_t factor = (pt->fBConstrain) ? factor1: factor2; + // + Int_t found,foundable,shared; + pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE); + Float_t sharedfactor = Float_t(shared+1)/Float_t(found+1); + Bool_t itsgold =kFALSE; + if (pt->fEsd){ + Int_t dummy[12]; + if (pt->fEsd->GetITSclusters(dummy)>4) itsgold= kTRUE; + } + if (!itsgold){ + // + if (Float_t(shared+1)/Float_t(found+1)>factor){ + if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks + delete arr->RemoveAt(trackindex); + continue; + } + if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks + delete arr->RemoveAt(trackindex); + continue; + } + } + + good++; + if (sharedfactor>0.4) continue; + if (pt->GetKinkIndexes()[0]>0) continue; + for (Int_t i=first; iGetClusterIndex2(i); + // if (index<0 || index&0x8000 ) continue; + if (index<0 || index&0x8000 ) continue; + AliTPCclusterMI *c= pt->fClusterPointer[i]; + if (!c) continue; + c->Use(10); + } + } + fNtracks = good; + if (fDebug>0){ + Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks); + } + delete []quality; + delete []indexes; +} + +void AliTPCtrackerMI::UnsignClusters() { // // loop over all clusters and unsign them @@ -2415,20 +2576,50 @@ Int_t AliTPCtrackerMI::RefitInward(AliESD *event) fIteration=2; //PrepareForProlongation(fSeeds,1); PropagateForward2(fSeeds); + + Int_t ntracks=0; Int_t nseed = fSeeds->GetEntriesFast(); for (Int_t i=0;iUncheckedAt(i); + if (!seed) continue; + if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks + seed->PropagateTo(fParam->GetInnerRadiusLow()); + seed->UpdatePoints(); AliESDtrack *esd=event->GetTrack(i); seed->CookdEdx(0.02,0.6); CookLabel(seed,0.1); //For comparison only - if (seed->GetNumberOfClusters()>20){ - esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit); + // + if (AliTPCReconstructor::StreamLevel()>0 && seed!=0&&esd!=0) { + TTreeSRedirector &cstream = *fDebugStreamer; + cstream<<"Crefit"<< + "Esd.="<GetNumberOfClusters()>15){ + esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit); + esd->SetTPCPoints(seed->GetPoints()); + esd->SetTPCPointsF(seed->fNFoundable); + Int_t ndedx = seed->fNCDEDX[0]+seed->fNCDEDX[1]+seed->fNCDEDX[2]+seed->fNCDEDX[3]; + Float_t sdedx = (seed->fSDEDX[0]+seed->fSDEDX[1]+seed->fSDEDX[2]+seed->fSDEDX[3])*0.25; + Float_t dedx = seed->GetdEdx(); + esd->SetTPCsignal(dedx, sdedx, ndedx); + // + // add seed to the esd track in Calib level + // + if (AliTPCReconstructor::StreamLevel()>0){ + AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE); + esd->AddCalibObject(seedCopy); + } + ntracks++; } else{ //printf("problem\n"); } } + //FindKinks(fSeeds,event); + Info("RefitInward","Number of refitted tracks %d",ntracks); fEvent =0; //WriteTracks(); return 0; @@ -2443,16 +2634,38 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESD *event) fEvent = event; fIteration = 1; - ReadSeeds(event,0); - PropagateBack(fSeeds); + ReadSeeds(event,1); + PropagateBack(fSeeds); + RemoveUsed2(fSeeds,0.4,0.4,20); + // Int_t nseed = fSeeds->GetEntriesFast(); + Int_t ntracks=0; for (Int_t i=0;iUncheckedAt(i); + if (!seed) continue; + if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks + seed->UpdatePoints(); AliESDtrack *esd=event->GetTrack(i); seed->CookdEdx(0.02,0.6); CookLabel(seed,0.1); //For comparison only - esd->UpdateTrackParams(seed,AliESDtrack::kTPCout); + if (seed->GetNumberOfClusters()>15){ + esd->UpdateTrackParams(seed,AliESDtrack::kTPCout); + esd->SetTPCPoints(seed->GetPoints()); + esd->SetTPCPointsF(seed->fNFoundable); + Int_t ndedx = seed->fNCDEDX[0]+seed->fNCDEDX[1]+seed->fNCDEDX[2]+seed->fNCDEDX[3]; + Float_t sdedx = (seed->fSDEDX[0]+seed->fSDEDX[1]+seed->fSDEDX[2]+seed->fSDEDX[3])*0.25; + Float_t dedx = seed->GetdEdx(); + esd->SetTPCsignal(dedx, sdedx, ndedx); + ntracks++; + Int_t eventnumber = event->GetEventNumber();// patch 28 fev 06 + (*fDebugStreamer)<<"Cback"<< + "Tr0.="<GetTrack(i); - ULong_t status=esd->GetStatus(); + ULong_t status=esd->GetStatus(); + if (!(status&AliESDtrack::kTPCin)) continue; AliTPCtrack t(*esd); - AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha()); - if ((status==AliESDtrack::kTPCin)&&(direction==1)) seed->ResetCovariance(); + t.SetNumberOfClusters(0); + // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha()); + AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/); + for (Int_t ikink=0;ikink<3;ikink++) { + Int_t index = esd->GetKinkIndex(ikink); + seed->GetKinkIndexes()[ikink] = index; + if (index==0) continue; + index = TMath::Abs(index); + AliESDkink * kink = fEvent->GetKink(index-1); + if (kink&&esd->GetKinkIndex(ikink)<0){ + if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2); + if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0); + } + if (kink&&esd->GetKinkIndex(ikink)>0){ + if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6); + if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4); + } + + } + if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(); if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(); + if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) { + fSeeds->AddAt(0,i); + delete seed; + continue; + } + if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) { + Double_t par0[5],par1[5],alpha,x; + esd->GetInnerExternalParameters(alpha,x,par0); + esd->GetExternalParameters(x,par1); + Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4])); + Double_t delta2 = TMath::Abs(par0[3]-par1[3]); + Double_t trdchi2=0; + if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls(); + //reset covariance if suspicious + if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.) + seed->ResetCovariance(); + } // // // rotate to the local coordinate system - - fSectors=fInnerSec; fN=fkNIS; - + // + fSectors=fInnerSec; fN=fkNIS; Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift(); if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); if (alpha < 0. ) alpha += 2.*TMath::Pi(); Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN; alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift(); + if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi(); + if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi(); alpha-=seed->GetAlpha(); - if (!seed->Rotate(alpha)) continue; + if (!seed->Rotate(alpha)) { + delete seed; + continue; + } seed->fEsd = esd; - // - //seed->PropagateTo(fSectors->GetX(0)); - // - // Int_t index = esd->GetTPCindex(); - //AliTPCseed * seed2= (AliTPCseed*)fSeeds->At(index); - //if (direction==2){ - // AliTPCseed * seed2 = ReSeed(seed,0.,0.5,1.); - // if (seed2) { - // delete seed; - // seed = seed2; - // } - //} - - fSeeds->AddLast(seed); + // sign clusters + if (esd->GetKinkIndex(0)<=0){ + for (Int_t irow=0;irow<160;irow++){ + Int_t index = seed->GetClusterIndex2(irow); + if (index>0){ + // + AliTPCclusterMI * cl = GetClusterMI(index); + seed->fClusterPointer[irow] = cl; + if (cl){ + if ((index & 0x8000)==0){ + cl->Use(10); // accepted cluster + }else{ + cl->Use(6); // close cluster not accepted + } + }else{ + Info("ReadSeeds","Not found cluster"); + } + } + } + } + fSeeds->AddAt(seed,i); } } @@ -3371,15 +3631,14 @@ AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point Int_t pp2=0; Double_t x0[3],x1[3],x2[3]; - x0[0]=-1; - x0[0]=-1; - x0[0]=-1; + for (Int_t i=0;i<3;i++){ + x0[i]=-1; + x1[i]=-1; + x2[i]=-1; + } // find track position at given ratio of the length - Int_t sec0, sec1, sec2; - sec0=0; - sec1=0; - sec2=0; + Int_t sec0=0, sec1=0, sec2=0; Int_t index=-1; Int_t clindex; for (Int_t i=0;i<160;i++){ @@ -3606,282 +3865,1516 @@ AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, F return seed; } -Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed, Float_t th) + +AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward) { // // - // - for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0; - // - if (TMath::Abs(seed->GetC())>0.01) return 0; - // - - Float_t x[160], y[160], erry[160], z[160], errz[160]; - Int_t sec[160]; - Float_t xt[160], yt[160], zt[160]; - Int_t i1 = 200; - Int_t i2 = 0; - Int_t secm = -1; - Int_t padm = -1; - Int_t middle = seed->GetNumberOfClusters()/2; - // - // - // find central sector, get local cooordinates - Int_t count = 0; - for (Int_t i=seed->fFirstPoint;i<=seed->fLastPoint;i++) { - sec[i]= seed->GetClusterSector(i)%18; - x[i] = GetXrow(i); - if (sec[i]>=0) { - AliTPCclusterMI * cl = seed->fClusterPointer[i]; - // if (cl==0) cl = GetClusterMI(seed->GetClusterIndex2(i)); - if (cl==0) { - sec[i] = -1; - continue; + //reseed using founded clusters + // + Double_t xyz[3][3]; + Int_t row[3]={0,0,0}; + Int_t sec[3]={0,0,0}; + // + // forward direction + if (forward){ + for (Int_t irow=r0;irow<160;irow++){ + if (track->GetClusterIndex(irow)>0){ + row[0] = irow; + break; } - // - // - if (i>i2) i2 = i; //last point with cluster - if (i2GetY(); - z[i] = cl->GetZ(); - AliTPCTrackerPoint * point = seed->GetTrackPoint(i); - xt[i] = x[i]; - yt[i] = point->GetY(); - zt[i] = point->GetZ(); - - if (point->GetX()>0){ - erry[i] = point->GetErrY(); - errz[i] = point->GetErrZ(); + } + for (Int_t irow=160;irow>r0;irow--){ + if (track->GetClusterIndex(irow)>0){ + row[2] = irow; + break; } - - count++; - if (countrow[0];irow--){ + if (track->GetClusterIndex(irow)>0){ + row[1] = irow; + break; } } - } - // - // rotate position to global coordinate system connected to sector at last the point - // - for (Int_t i=i1;i<=i2;i++){ - // - if (sec[i]<0) continue; - Double_t alpha = (sec[i2]-sec[i])*fSectors->GetAlpha(); - Double_t cs = TMath::Cos(alpha); - Double_t sn = TMath::Sin(alpha); - Float_t xx2= x[i]*cs+y[i]*sn; - Float_t yy2= -x[i]*sn+y[i]*cs; - x[i] = xx2; - y[i] = yy2; // - xx2= xt[i]*cs+yt[i]*sn; - yy2= -xt[i]*sn+yt[i]*cs; - xt[i] = xx2; - yt[i] = yy2; - - } - //get "state" vector - Double_t xh[5],xm = x[padm]; - xh[0]=yt[i2]; - xh[1]=zt[i2]; - xh[4]=F1(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]); - xh[2]=F2(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]); - xh[3]=F3n(xt[i2],yt[i2],xt[i1],yt[i1],zt[i2],zt[i1],xh[4]); - // - // - for (Int_t i=i1;i<=i2;i++){ - Double_t yy,zz; - if (sec[i]<0) continue; - GetProlongation(x[i2], x[i],xh,yy,zz); - if (TMath::Abs(y[i]-yy)>4||TMath::Abs(z[i]-zz)>4){ - //Double_t xxh[5]; - //xxh[4]=F1old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]); - //xxh[2]=F2old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]); - Error("AliTPCtrackerMI::CheckKinkPoint","problem\n"); - } - y[i] = y[i] - yy; - z[i] = z[i] - zz; } - Float_t dyup[160],dydown[160], dzup[160], dzdown[160]; - Float_t yup[160], ydown[160], zup[160], zdown[160]; - - AliTPCpolyTrack ptrack1,ptrack2; - // - // derivation up - for (Int_t i=i1;i<=i2;i++){ - AliTPCclusterMI * cl = seed->fClusterPointer[i]; - if (!cl) continue; - if (cl->GetType()<0) continue; - if (cl->GetType()>10) continue; - - if (sec[i]>=0){ - ptrack1.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1); - } - if (ptrack1.GetN()>4.){ - ptrack1.UpdateParameters(); - Double_t ddy,ddz; - ptrack1.GetFitDerivation(x[i]-xm,ddy,ddz); - Double_t yy,zz; - ptrack1.GetFitPoint(x[i]-xm,yy,zz); - - dyup[i] = ddy; - dzup[i] = ddz; - yup[i] = yy; - zup[i] = zz; - - } - else{ - dyup[i]=0.; //not enough points + if (!forward){ + for (Int_t irow=0;irowGetClusterIndex(irow)>0){ + row[0] = irow; + break; + } } + for (Int_t irow=r0;irow>0;irow--){ + if (track->GetClusterIndex(irow)>0){ + row[2] = irow; + break; + } + } + for (Int_t irow=row[2]-15;irow>row[0];irow--){ + if (track->GetClusterIndex(irow)>0){ + row[1] = irow; + break; + } + } } // - // derivation down - for (Int_t i=i2;i>=i1;i--){ - AliTPCclusterMI * cl = seed->fClusterPointer[i]; - if (!cl) continue; - if (cl->GetType()<0) continue; - if (cl->GetType()>10) continue; - if (sec[i]>=0){ - ptrack2.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1); + if ((row[2]-row[0])<20) return 0; + if (row[1]==0) return 0; + // + // + //Get cluster and sector position + for (Int_t ipoint=0;ipoint<3;ipoint++){ + Int_t clindex = track->GetClusterIndex2(row[ipoint]); + AliTPCclusterMI * cl = GetClusterMI(clindex); + if (cl==0) { + //Error("Bug\n"); + // AliTPCclusterMI * cl = GetClusterMI(clindex); + return 0; } - if (ptrack2.GetN()>4){ - ptrack2.UpdateParameters(); - Double_t ddy,ddz; - ptrack2.GetFitDerivation(x[i]-xm,ddy,ddz); - Double_t yy,zz; - ptrack2.GetFitPoint(x[i]-xm,yy,zz); - - dydown[i] = ddy; - dzdown[i] = ddz; - ydown[i] = yy; - zdown[i] = zz; + sec[ipoint] = ((clindex&0xff000000)>>24)%18; + xyz[ipoint][0] = GetXrow(row[ipoint]); + AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]); + if (point&&ipoint<2){ + // + xyz[ipoint][1] = point->GetY(); + xyz[ipoint][2] = point->GetZ(); } else{ - dydown[i]=0.; //not enough points + xyz[ipoint][1] = cl->GetY(); + xyz[ipoint][2] = cl->GetZ(); } } // // - // find maximal difference of the derivation - for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0; - - - for (Int_t i=i1+10;i th){ - seed->fKinkPoint[0] = i; - seed->fKinkPoint[1] = ddy; - seed->fKinkPoint[2] = ddz; - th = ddy+ddz; - } - } + // + // + // Calculate seed state vector and covariance matrix - if (fTreeDebug){ - // - //write information to the debug tree - TBranch * br = fTreeDebug->GetBranch("debug"); - TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2"); - arr->ExpandCreateFast(i2-i1); - br->SetAddress(&arr); - // - AliTPCclusterMI cldummy; - cldummy.SetQ(0); - AliTPCTrackPoint2 pdummy; - pdummy.GetTPoint().fIsShared = 10; - // - Double_t alpha = sec[i2]*fSectors->GetAlpha(); - Double_t cs = TMath::Cos(alpha); - Double_t sn = TMath::Sin(alpha); - - for (Int_t i=i1;iUncheckedAt(i-i1); - //cluster info - AliTPCclusterMI * cl0 = seed->fClusterPointer[i]; - // - AliTPCTrackerPoint * point = seed->GetTrackPoint(i); - - if (cl0){ - Double_t x = GetXrow(i); - trpoint->GetTPoint() = *point; - trpoint->GetCPoint() = *cl0; - trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ())); - trpoint->fID = seed->GetUniqueID(); - trpoint->fLab = seed->GetLabel(); - // - trpoint->fGX = cs *x + sn*point->GetY(); - trpoint->fGY = -sn *x + cs*point->GetY() ; - trpoint->fGZ = point->GetZ(); - // - trpoint->fDY = y[i]; - trpoint->fDZ = z[i]; - // - trpoint->fDYU = dyup[i]; - trpoint->fDZU = dzup[i]; - // - trpoint->fDYD = dydown[i]; - trpoint->fDZD = dzdown[i]; - // - if (TMath::Abs(dyup[i])>0.00000000001 &&TMath::Abs(dydown[i])>0.00000000001){ - trpoint->fDDY = dydown[i]-dyup[i]; - trpoint->fDDZ = dzdown[i]-dzup[i]; - }else{ - trpoint->fDDY = 0.; - trpoint->fDDZ = 0.; - } - } - else{ - *trpoint = pdummy; - trpoint->GetCPoint()= cldummy; - trpoint->fID = -1; - } - // - } - fTreeDebug->Fill(); - } + Double_t alpha, cs,sn, xx2,yy2; + // + alpha = (sec[1]-sec[2])*fSectors->GetAlpha(); + cs = TMath::Cos(alpha); + sn = TMath::Sin(alpha); + xx2= xyz[1][0]*cs-xyz[1][1]*sn; + yy2= xyz[1][0]*sn+xyz[1][1]*cs; + xyz[1][0] = xx2; + xyz[1][1] = yy2; + // + alpha = (sec[0]-sec[2])*fSectors->GetAlpha(); + cs = TMath::Cos(alpha); + sn = TMath::Sin(alpha); + xx2= xyz[0][0]*cs-xyz[0][1]*sn; + yy2= xyz[0][0]*sn+xyz[0][1]*cs; + xyz[0][0] = xx2; + xyz[0][1] = yy2; + // + // + // + Double_t x[5],c[15]; + // + x[0]=xyz[2][1]; + x[1]=xyz[2][2]; + x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]); + x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]); + x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]); + // + Double_t sy =0.1, sz =0.1; + // + Double_t sy1=0.2, sz1=0.2; + Double_t sy2=0.2, sz2=0.2; + Double_t sy3=0.2; + // + Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy; + Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy; + Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy; + Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy; + Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy; + Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy; + // + Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy; + Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz; + Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy; + Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz; - return 0; + c[0]=sy1; + c[1]=0.; c[2]=sz1; + c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23; + c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22; + c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34; + c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23; + c[13]=f30*sy1*f40+f32*sy2*f42; + c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43; -} - - - + // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]); + AliTPCseed *seed=new AliTPCseed(0, x, c, xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift()); + seed->fLastPoint = row[2]; + seed->fFirstPoint = row[2]; + for (Int_t i=row[0];ifIndex[i] = track->fIndex[i]; + } + return seed; +} -AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t) +void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) { // - // reseed - refit - track + // find kinks // - Int_t first = 0; - // Int_t last = fSectors->GetNRows()-1; // - if (fSectors == fOuterSec){ - first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows()); - //last = + + TObjArray *kinks= new TObjArray(10000); + // TObjArray *v0s= new TObjArray(10000); + Int_t nentries = array->GetEntriesFast(); + AliHelix *helixes = new AliHelix[nentries]; + Int_t *sign = new Int_t[nentries]; + Int_t *nclusters = new Int_t[nentries]; + Float_t *alpha = new Float_t[nentries]; + AliESDkink * kink = new AliESDkink(); + Int_t * usage = new Int_t[nentries]; + Float_t *zm = new Float_t[nentries]; + Float_t *z0 = new Float_t[nentries]; + Float_t *fim = new Float_t[nentries]; + Float_t *shared = new Float_t[nentries]; + Bool_t *circular = new Bool_t[nentries]; + Float_t *dca = new Float_t[nentries]; + //const AliESDVertex * primvertex = esd->GetVertex(); + // + // nentries = array->GetEntriesFast(); + // + + // + // + for (Int_t i=0;iAt(i); + if (!track) continue; + track->fCircular =0; + shared[i] = kFALSE; + track->UpdatePoints(); + if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){ + } + nclusters[i]=track->GetNumberOfClusters(); + alpha[i] = track->GetAlpha(); + new (&helixes[i]) AliHelix(*track); + Double_t xyz[3]; + helixes[i].Evaluate(0,xyz); + sign[i] = (track->GetC()>0) ? -1:1; + Double_t x,y,z; + x=160; + if (track->GetProlongation(x,y,z)){ + zm[i] = z; + fim[i] = alpha[i]+TMath::ATan2(y,x); + } + else{ + zm[i] = track->GetZ(); + fim[i] = alpha[i]; + } + z0[i]=1000; + circular[i]= kFALSE; + if (track->GetProlongation(0,y,z)) z0[i] = z; + dca[i] = track->GetD(0,0); } - else - first = t->fFirstPoint; // - AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9); - FollowBackProlongation(*t,fSectors->GetNRows()-1); - t->Reset(kFALSE); - FollowProlongation(*t,first); - return seed; -} - - - - - - - -//_____________________________________________________________________________ + // + TStopwatch timer; + timer.Start(); + Int_t ncandidates =0; + Int_t nall =0; + Int_t ntracks=0; + Double_t phase[2][2],radius[2]; + + // + // Find circling track + TTreeSRedirector &cstream = *fDebugStreamer; + // + for (Int_t i0=0;i0At(i0); + if (!track0) continue; + if (track0->fN<40) continue; + if (TMath::Abs(1./track0->fP4)>200) continue; + for (Int_t i1=i0+1;i1At(i1); + if (!track1) continue; + if (track1->fN<40) continue; + if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue; + if (track0->fBConstrain&&track1->fBConstrain) continue; + if (TMath::Abs(1./track1->fP4)>200) continue; + if (track1->fP4*track0->fP4>0) continue; + if (track1->fP3*track0->fP3>0) continue; + if (max(TMath::Abs(1./track0->fP4),TMath::Abs(1./track1->fP4))>190) continue; + if (track0->fBConstrain&&TMath::Abs(track1->fP4)fP4)) continue; //returning - lower momenta + if (track1->fBConstrain&&TMath::Abs(track0->fP4)fP4)) continue; //returning - lower momenta + // + Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1])); + if (mindcar<5) continue; + Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ())); + if (mindcaz<5) continue; + if (mindcar+mindcaz<20) continue; + // + // + Float_t xc0 = helixes[i0].GetHelix(6); + Float_t yc0 = helixes[i0].GetHelix(7); + Float_t r0 = helixes[i0].GetHelix(8); + Float_t xc1 = helixes[i1].GetHelix(6); + Float_t yc1 = helixes[i1].GetHelix(7); + Float_t r1 = helixes[i1].GetHelix(8); + + Float_t rmean = (r0+r1)*0.5; + Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0)); + //if (delta>30) continue; + if (delta>rmean*0.25) continue; + if (TMath::Abs(r0-r1)/rmean>0.3) continue; + // + Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10); + if (npoints==0) continue; + helixes[i0].GetClosestPhases(helixes[i1], phase); + // + Double_t xyz0[3]; + Double_t xyz1[3]; + Double_t hangles[3]; + helixes[i0].Evaluate(phase[0][0],xyz0); + helixes[i1].Evaluate(phase[0][1],xyz1); + + helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles); + Double_t deltah[2],deltabest; + if (hangles[2]<2.8) continue; + /* + cstream<<"C"<fLab<fLab<< + track0->fP3<fP3<< + track0->fP4<fP4<< + delta<0){ + Int_t ibest=0; + helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2); + if (npoints==2){ + helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2); + if (deltah[1]6) continue; + if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue; + Bool_t sign =kFALSE; + if (hangles[2]>3.06) sign =kTRUE; + // + if (sign){ + circular[i0] = kTRUE; + circular[i1] = kTRUE; + if (TMath::Abs(track0->fP4)fP4)){ + track0->fCircular += 1; + track1->fCircular += 2; + } + else{ + track1->fCircular += 1; + track0->fCircular += 2; + } + } + if (sign&&AliTPCReconstructor::StreamLevel()>1){ + //debug stream + cstream<<"Curling"<< + "lab0="<fLab<< + "lab1="<fLab<< + "Tr0.="<At(i); + ntracks++; + // + Double_t cradius0 = 40*40; + Double_t cradius1 = 270*270; + Double_t cdist1=8.; + Double_t cdist2=8.; + Double_t cdist3=0.55; + for (Int_t j =i+1;j200) continue; + if ( (nclusters[i]+nclusters[j])<80) continue; + if ( TMath::Abs(zm[i]-zm[j])>60.) continue; + if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue; + //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2]; + Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20); + if (npoints<1) continue; + // cuts on radius + if (npoints==1){ + if (radius[0]cradius1) continue; + } + else{ + if ( (radius[0]cradius1) && (radius[1]cradius1) ) continue; + } + // + Double_t delta1=10000,delta2=10000; + // cuts on the intersection radius + helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1); + if (radius[0]<20&&delta1<1) continue; //intersection at vertex + if (radius[0]<10&&delta1<3) continue; //intersection at vertex + if (npoints==2){ + helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2); + if (radius[1]<20&&delta2<1) continue; //intersection at vertex + if (radius[1]<10&&delta2<3) continue; //intersection at vertex + } + // + Double_t distance1 = TMath::Min(delta1,delta2); + if (distance1>cdist1) continue; // cut on DCA linear approximation + // + npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20); + helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1); + if (radius[0]<20&&delta1<1) continue; //intersection at vertex + if (radius[0]<10&&delta1<3) continue; //intersection at vertex + // + if (npoints==2){ + helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2); + if (radius[1]<20&&delta2<1) continue; //intersection at vertex + if (radius[1]<10&&delta2<3) continue; //intersection at vertex + } + distance1 = TMath::Min(delta1,delta2); + Float_t rkink =0; + if (delta1cdist2) continue; + // + // + AliTPCseed * track1 = (AliTPCseed*)array->At(j); + // + // + Int_t row0 = GetRowNumber(rkink); + if (row0<10) continue; + if (row0>150) continue; + // + // + Float_t dens00=-1,dens01=-1; + Float_t dens10=-1,dens11=-1; + // + Int_t found,foundable,shared; + track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE); + if (foundable>5) dens00 = Float_t(found)/Float_t(foundable); + track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE); + if (foundable>5) dens01 = Float_t(found)/Float_t(foundable); + // + track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE); + if (foundable>10) dens10 = Float_t(found)/Float_t(foundable); + track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE); + if (foundable>10) dens11 = Float_t(found)/Float_t(foundable); + // + if (dens00dens10 && dens01>dens11) continue; + if (TMath::Max(dens00,dens10)<0.1) continue; + if (TMath::Max(dens01,dens11)<0.3) continue; + // + if (TMath::Min(dens00,dens10)>0.6) continue; + if (TMath::Min(dens01,dens11)>0.6) continue; + + // + AliTPCseed * ktrack0, *ktrack1; + if (dens00>dens10){ + ktrack0 = track0; + ktrack1 = track1; + } + else{ + ktrack0 = track1; + ktrack1 = track0; + } + if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle + AliExternalTrackParam paramm(*ktrack0); + AliExternalTrackParam paramd(*ktrack1); + if (row0>60&&ktrack1->GetReference().GetX()>90.) new (¶md) AliExternalTrackParam(ktrack1->GetReference()); + // + // + kink->SetMother(paramm); + kink->SetDaughter(paramd); + kink->Update(); + + Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]}; + Int_t index[4]; + fParam->Transform0to1(x,index); + fParam->Transform1to2(x,index); + row0 = GetRowNumber(x[0]); + + if (kink->GetR()<100) continue; + if (kink->GetR()>240) continue; + if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume + if (kink->GetDistance()>cdist3) continue; + Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate + if (dird<0) continue; + + Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate + if (dirm<0) continue; + Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]); + if (mpt<0.2) continue; + + if (mpt<1){ + //for high momenta momentum not defined well in first iteration + Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->P(); + if (qt>0.35) continue; + } + + kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0); + kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1); + if (dens00>dens10){ + kink->SetTPCDensity(dens00,0,0); + kink->SetTPCDensity(dens01,0,1); + kink->SetTPCDensity(dens10,1,0); + kink->SetTPCDensity(dens11,1,1); + kink->SetIndex(i,0); + kink->SetIndex(j,1); + } + else{ + kink->SetTPCDensity(dens10,0,0); + kink->SetTPCDensity(dens11,0,1); + kink->SetTPCDensity(dens00,1,0); + kink->SetTPCDensity(dens01,1,1); + kink->SetIndex(j,0); + kink->SetIndex(i,1); + } + + if (mpt<1||kink->GetAngle(2)>0.1){ + // angle and densities not defined yet + if (kink->GetTPCDensityFactor()<0.8) continue; + if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue; + if (kink->GetAngle(2)*ktrack0->P()<0.003) continue; //too small angle + if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue; + if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue; + + Float_t criticalangle = track0->fC22+track0->fC33; + criticalangle+= track1->fC22+track1->fC33; + criticalangle= 3*TMath::Sqrt(criticalangle); + if (criticalangle>0.02) criticalangle=0.02; + if (kink->GetAngle(2)GetAngle(2))); // overlap region defined + Float_t shapesum =0; + Float_t sum = 0; + for ( Int_t row = row0-drow; row155) continue; + if (ktrack0->fClusterPointer[row]){ + AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row); + shapesum+=point->GetSigmaY()+point->GetSigmaZ(); + sum++; + } + if (ktrack1->fClusterPointer[row]){ + AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row); + shapesum+=point->GetSigmaY()+point->GetSigmaZ(); + sum++; + } + } + if (sum<4){ + kink->SetShapeFactor(-1.); + } + else{ + kink->SetShapeFactor(shapesum/sum); + } + // esd->AddKink(kink); + kinks->AddLast(kink); + kink = new AliESDkink; + ncandidates++; + } + } + // + // sort the kinks according quality - and refit them towards vertex + // + Int_t nkinks = kinks->GetEntriesFast(); + Float_t *quality = new Float_t[nkinks]; + Int_t *indexes = new Int_t[nkinks]; + AliTPCseed *mothers = new AliTPCseed[nkinks]; + AliTPCseed *daughters = new AliTPCseed[nkinks]; + // + // + for (Int_t i=0;iAt(i); + // + // refit kinks towards vertex + // + Int_t index0 = kink->GetIndex(0); + Int_t index1 = kink->GetIndex(1); + AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0); + AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1); + // + Int_t sumn=ktrack0->fN+ktrack1->fN; + // + // Refit Kink under if too small angle + // + if (kink->GetAngle(2)<0.05){ + kink->SetTPCRow0(GetRowNumber(kink->GetR())); + Int_t row0 = kink->GetTPCRow0(); + Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); + // + // + Int_t last = row0-drow; + if (last<40) last=40; + if (lastfFirstPoint+25) last = ktrack0->fFirstPoint+25; + AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE); + // + // + Int_t first = row0+drow; + if (first>130) first=130; + if (first>ktrack1->fLastPoint-25) first = TMath::Max(ktrack1->fLastPoint-25,30); + AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE); + // + if (seed0 && seed1){ + kink->SetStatus(1,8); + if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9); + row0 = GetRowNumber(kink->GetR()); + sumn = seed0->fN+seed1->fN; + new (&mothers[i]) AliTPCseed(*seed0); + new (&daughters[i]) AliTPCseed(*seed1); + } + else{ + delete kinks->RemoveAt(i); + if (seed0) delete seed0; + if (seed1) delete seed1; + continue; + } + if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) { + delete kinks->RemoveAt(i); + if (seed0) delete seed0; + if (seed1) delete seed1; + continue; + } + // + delete seed0; + delete seed1; + } + // + if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win + } + TMath::Sort(nkinks,quality,indexes,kFALSE); + // + //remove double find kinks + // + for (Int_t ikink0=1;ikink0At(indexes[ikink0]); + if (!kink0) continue; + // + for (Int_t ikink1=0;ikink1At(indexes[ikink1]); + if (!kink1) continue; + // if not close kink continue + if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue; + if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue; + if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue; + // + AliTPCseed &mother0 = mothers[indexes[ikink0]]; + AliTPCseed &daughter0 = daughters[indexes[ikink0]]; + AliTPCseed &mother1 = mothers[indexes[ikink1]]; + AliTPCseed &daughter1 = daughters[indexes[ikink1]]; + Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2; + // + Int_t same = 0; + Int_t both = 0; + Int_t samem = 0; + Int_t bothm = 0; + Int_t samed = 0; + Int_t bothd = 0; + // + for (Int_t i=0;i0 && mother1.fIndex[i]>0){ + both++; + bothm++; + if (mother0.fIndex[i]==mother1.fIndex[i]){ + same++; + samem++; + } + } + } + + for (Int_t i=row0;i<158;i++){ + if (daughter0.fIndex[i]>0 && daughter0.fIndex[i]>0){ + both++; + bothd++; + if (mother0.fIndex[i]==mother1.fIndex[i]){ + same++; + samed++; + } + } + } + Float_t ratio = Float_t(same+1)/Float_t(both+1); + Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1); + Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1); + if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) { + Int_t sum0 = mother0.fN+daughter0.fN; + Int_t sum1 = mother1.fN+daughter1.fN; + if (sum1>sum0){ + shared[kink0->GetIndex(0)]= kTRUE; + shared[kink0->GetIndex(1)]= kTRUE; + delete kinks->RemoveAt(indexes[ikink0]); + } + else{ + shared[kink1->GetIndex(0)]= kTRUE; + shared[kink1->GetIndex(1)]= kTRUE; + delete kinks->RemoveAt(indexes[ikink1]); + } + } + } + } + + + for (Int_t i=0;iAt(indexes[i]); + if (!kink) continue; + kink->SetTPCRow0(GetRowNumber(kink->GetR())); + Int_t index0 = kink->GetIndex(0); + Int_t index1 = kink->GetIndex(1); + if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue; + kink->SetMultiple(usage[index0],0); + kink->SetMultiple(usage[index1],1); + if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue; + if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue; + if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue; + if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue; + + AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0); + AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1); + if (!ktrack0 || !ktrack1) continue; + Int_t index = esd->AddKink(kink); + // + // + if ( ktrack0->fKinkIndexes[0]==0 && ktrack1->fKinkIndexes[0]==0) { //best kink + if (mothers[indexes[i]].fN>20 && daughters[indexes[i]].fN>20 && (mothers[indexes[i]].fN+daughters[indexes[i]].fN)>100){ + new (ktrack0) AliTPCseed(mothers[indexes[i]]); + new (ktrack1) AliTPCseed(daughters[indexes[i]]); + } + } + // + ktrack0->fKinkIndexes[usage[index0]] = -(index+1); + ktrack1->fKinkIndexes[usage[index1]] = (index+1); + usage[index0]++; + usage[index1]++; + } + // + // Remove tracks corresponding to shared kink's + // + for (Int_t i=0;iAt(i); + if (!track0) continue; + if (track0->fKinkIndexes[0]!=0) continue; + if (shared[i]) delete array->RemoveAt(i); + } + + // + // + RemoveUsed2(array,0.5,0.4,30); + UnsignClusters(); + for (Int_t i=0;iAt(i); + if (!track0) continue; + track0->CookdEdx(0.02,0.6); + track0->CookPID(); + } + // + for (Int_t i=0;iAt(i); + if (!track0) continue; + if (track0->Pt()<1.4) continue; + //remove double high momenta tracks - overlapped with kink candidates + Int_t shared=0; + Int_t all =0; + for (Int_t icl=track0->fFirstPoint;iclfLastPoint; icl++){ + if (track0->fClusterPointer[icl]!=0){ + all++; + if (track0->fClusterPointer[icl]->IsUsed(10)) shared++; + } + } + if (Float_t(shared+1)/Float_t(nall+1)>0.5) { + delete array->RemoveAt(i); + continue; + } + // + if (track0->fKinkIndexes[0]!=0) continue; + if (track0->GetNumberOfClusters()<80) continue; + + AliTPCseed *pmother = new AliTPCseed(); + AliTPCseed *pdaughter = new AliTPCseed(); + AliESDkink *pkink = new AliESDkink; + + AliTPCseed & mother = *pmother; + AliTPCseed & daughter = *pdaughter; + AliESDkink & kink = *pkink; + if (CheckKinkPoint(track0,mother,daughter, kink)){ + if (mother.fN<30||daughter.fN<20) { + delete pmother; + delete pdaughter; + delete pkink; + continue; //too short tracks + } + if (mother.Pt()<1.4) { + delete pmother; + delete pdaughter; + delete pkink; + continue; + } + Int_t row0= kink.GetTPCRow0(); + if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) { + delete pmother; + delete pdaughter; + delete pkink; + continue; + } + // + Int_t index = esd->AddKink(&kink); + mother.fKinkIndexes[0] = -(index+1); + daughter.fKinkIndexes[0] = index+1; + if (mother.fN>50) { + delete array->RemoveAt(i); + array->AddAt(new AliTPCseed(mother),i); + } + else{ + array->AddLast(new AliTPCseed(mother)); + } + array->AddLast(new AliTPCseed(daughter)); + for (Int_t icl=0;iclUse(20); + } + // + for (Int_t icl=row0;icl<158;icl++) { + if (daughter.fClusterPointer[icl]) daughter.fClusterPointer[icl]->Use(20); + } + // + } + delete pmother; + delete pdaughter; + delete pkink; + } + + delete [] daughters; + delete [] mothers; + // + // + delete [] dca; + delete []circular; + delete []shared; + delete []quality; + delete []indexes; + // + delete kink; + delete[]fim; + delete[] zm; + delete[] z0; + delete [] usage; + delete[] alpha; + delete[] nclusters; + delete[] sign; + delete[] helixes; + kinks->Delete(); + delete kinks; + + printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall); + timer.Print(); +} + +void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESD *esd) +{ + // + // find V0s + // + // + TObjArray *tpcv0s = new TObjArray(100000); + Int_t nentries = array->GetEntriesFast(); + AliHelix *helixes = new AliHelix[nentries]; + Int_t *sign = new Int_t[nentries]; + Float_t *alpha = new Float_t[nentries]; + Float_t *z0 = new Float_t[nentries]; + Float_t *dca = new Float_t[nentries]; + Float_t *sdcar = new Float_t[nentries]; + Float_t *cdcar = new Float_t[nentries]; + Float_t *pulldcar = new Float_t[nentries]; + Float_t *pulldcaz = new Float_t[nentries]; + Float_t *pulldca = new Float_t[nentries]; + Bool_t *isPrim = new Bool_t[nentries]; + const AliESDVertex * primvertex = esd->GetVertex(); + Double_t zvertex = primvertex->GetZv(); + // + // nentries = array->GetEntriesFast(); + // + for (Int_t i=0;iAt(i); + if (!track) continue; + track->GetV0Indexes()[0] = 0; //rest v0 indexes + track->GetV0Indexes()[1] = 0; //rest v0 indexes + track->GetV0Indexes()[2] = 0; //rest v0 indexes + // + alpha[i] = track->GetAlpha(); + new (&helixes[i]) AliHelix(*track); + Double_t xyz[3]; + helixes[i].Evaluate(0,xyz); + sign[i] = (track->GetC()>0) ? -1:1; + Double_t x,y,z; + x=160; + z0[i]=1000; + if (track->GetProlongation(0,y,z)) z0[i] = z; + dca[i] = track->GetD(0,0); + // + // dca error parrameterezation + pulls + // + sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->fP4)*(100*track->fP4)); + if (TMath::Abs(track->fP3)>1) sdcar[i]*=2.5; + cdcar[i] = TMath::Exp((TMath::Abs(track->fP4)-0.0106)*525.3); + pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i]; + pulldcaz[i] = (z0[i]-zvertex)/sdcar[i]; + pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]); + if (track->fTPCr[1]+track->fTPCr[2]+track->fTPCr[3]>0.5) { + if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut + } + if (track->fTPCr[4]>0.5) { + if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut + } + if (track->fTPCr[0]>0.4) { + isPrim[i]=kFALSE; //electron no sigma cut + } + } + // + // + TStopwatch timer; + timer.Start(); + Int_t ncandidates =0; + Int_t nall =0; + Int_t ntracks=0; + Double_t phase[2][2],radius[2]; + // + // Finf V0s loop + // + // + // // + TTreeSRedirector &cstream = *fDebugStreamer; + Float_t fprimvertex[3]={GetX(),GetY(),GetZ()}; + AliESDV0MI vertex; + Double_t cradius0 = 10*10; + Double_t cradius1 = 200*200; + Double_t cdist1=3.; + Double_t cdist2=4.; + Double_t cpointAngle = 0.95; + // + Double_t delta[2]={10000,10000}; + for (Int_t i =0;iAt(i); + if (!track0) continue; + if (AliTPCReconstructor::StreamLevel()>0){ + cstream<<"Tracks"<< + "Tr0.="<fP4<0) continue; + if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink + // + if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue; + ntracks++; + // debug output + + + for (Int_t j =0;jAt(j); + if (!track1) continue; + if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink + if (sign[j]*sign[i]>0) continue; + if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue; + if (track0->fCircular+track1->fCircular>1) continue; //circling -returning track + nall++; + // + // DCA to prim vertex cut + // + // + delta[0]=10000; + delta[1]=10000; + + Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2); + if (npoints<1) continue; + Int_t iclosest=0; + // cuts on radius + if (npoints==1){ + if (radius[0]cradius1) continue; + helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]); + if (delta[0]>cdist1) continue; + } + else{ + if (TMath::Max(radius[0],radius[1])cradius1) continue; + helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]); + helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]); + if (delta[1]cdist1) continue; + } + helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]); + if (radius[iclosest]cradius1 || delta[iclosest]>cdist1) continue; + // + Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex); + if (pointAnglefTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate + Double_t pointAngle2 = vertex.GetPointAngle(); + //continue; + if (vertex.GetPointAngle()2&&(!isGamma)) continue; // point angle cut + Float_t sigmae = 0.15*0.15; + if (vertex.GetRr()<80) + sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.); + sigmae+= TMath::Sqrt(sigmae); + if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue; + Float_t densb0=0,densb1=0,densa0=0,densa1=0; + Int_t row0 = GetRowNumber(vertex.GetRr()); + if (row0>15){ + if (vertex.GetDist2()>0.2) continue; + densb0 = track0->Density2(0,row0-5); + densb1 = track1->Density2(0,row0-5); + if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex + densa0 = track0->Density2(row0+5,row0+40); + densa1 = track1->Density2(row0+5,row0+40); + if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex + } + else{ + densa0 = track0->Density2(0,40); //cluster density + densa1 = track1->Density2(0,40); //cluster density + if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue; + } + vertex.SetLab(0,track0->GetLabel()); + vertex.SetLab(1,track1->GetLabel()); + vertex.SetChi2After((densa0+densa1)*0.5); + vertex.SetChi2Before((densb0+densb1)*0.5); + vertex.SetIndex(0,i); + vertex.SetIndex(1,j); + vertex.SetStatus(1); // TPC v0 candidate + vertex.SetRp(track0->fTPCr); + vertex.SetRm(track1->fTPCr); + tpcv0s->AddLast(new AliESDV0MI(vertex)); + ncandidates++; + { + Int_t eventNr = esd->GetEventNumber(); + Double_t radiusm= (delta[0]0) + cstream<<"V0"<< + "Event="<At(i); + quality[i] = 1./(1.00001-v0->GetPointAngle()); //base point angle + // quality[i] /= (0.5+v0->GetDist2()); + // quality[i] *= v0->GetChi2After(); //density factor + Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull + Int_t index0 = v0->GetIndex(0); + Int_t index1 = v0->GetIndex(1); + AliTPCseed * track0 = (AliTPCseed*)array->At(index0); + AliTPCseed * track1 = (AliTPCseed*)array->At(index1); + if (track0->fTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate + if (track0->fTPCr[4]>0.9||track1->fTPCr[4]>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate + } + + TMath::Sort(ncandidates,quality,indexes,kTRUE); + // + // + for (Int_t i=0;iAt(indexes[i]); + if (!v0) continue; + Int_t index0 = v0->GetIndex(0); + Int_t index1 = v0->GetIndex(1); + AliTPCseed * track0 = (AliTPCseed*)array->At(index0); + AliTPCseed * track1 = (AliTPCseed*)array->At(index1); + if (!track0||!track1) { + printf("Bug\n"); + continue; + } + Bool_t accept =kTRUE; //default accept + Int_t *v0indexes0 = track0->GetV0Indexes(); + Int_t *v0indexes1 = track1->GetV0Indexes(); + // + Int_t order0 = (v0indexes0[0]!=0) ? 1:0; + Int_t order1 = (v0indexes1[0]!=0) ? 1:0; + if (v0indexes0[1]!=0) order0 =2; + if (v0indexes1[1]!=0) order1 =2; + // + if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;} + if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;} + // + AliESDV0MI * v02 = v0; + if (accept){ + v0->SetOrder(0,order0); + v0->SetOrder(1,order1); + v0->SetOrder(1,order0+order1); + Int_t index = esd->AddV0MI(v0); + v02 = esd->GetV0MI(index); + v0indexes0[order0]=index; + v0indexes1[order1]=index; + naccepted++; + } + { + Int_t eventNr = esd->GetEventNumber(); + if (AliTPCReconstructor::StreamLevel()>0) + cstream<<"V02"<< + "Event="<240.) continue; + if (TMath::Abs(kinks[irow].GetR())<100.) continue; + // + Float_t normdist = TMath::Abs(param0[irow].fX-kinks[irow].GetR())*(0.1+kink.GetDistance()); + normdist/= (param0[irow].fN+param1[irow].fN+40.); + if (normdist < mindist){ + mindist = normdist; + index = irow; + } + } + // + if (index==-1) return 0; + // + // + param0[index].Reset(kTRUE); + FollowProlongation(param0[index],0); + // + new (&mother) AliTPCseed(param0[index]); + new (&daughter) AliTPCseed(param1[index]); // daughter in vertex + // + kink.SetMother(mother); + kink.SetDaughter(daughter); + kink.Update(); + kink.SetTPCRow0(GetRowNumber(kink.GetR())); + kink.SetTPCncls(param0[index].fN,0); + kink.SetTPCncls(param1[index].fN,1); + kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0); + kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1); + mother.SetLabel(kink.GetLabel(0)); + daughter.SetLabel(kink.GetLabel(1)); + + return 1; +} + + +void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){ + // + // update Kink quality information for mother after back propagation + // + if (seed->GetKinkIndex(0)>=0) return; + for (Int_t ikink=0;ikink<3;ikink++){ + Int_t index = seed->GetKinkIndex(ikink); + if (index>=0) break; + index = TMath::Abs(index)-1; + AliESDkink * kink = fEvent->GetKink(index); + //kink->fTPCdensity2[0][0]=-1; + //kink->fTPCdensity2[0][1]=-1; + kink->SetTPCDensity2(-1,0,0); + kink->SetTPCDensity2(1,0,1); + // + Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2))); + if (row0<15) row0=15; + // + Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2))); + if (row1>145) row1=145; + // + Int_t found,foundable,shared; + seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE); + if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),0,0); + seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE); + if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),0,1); + } + +} + +void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){ + // + // update Kink quality information for daughter after refit + // + if (seed->GetKinkIndex(0)<=0) return; + for (Int_t ikink=0;ikink<3;ikink++){ + Int_t index = seed->GetKinkIndex(ikink); + if (index<=0) break; + index = TMath::Abs(index)-1; + AliESDkink * kink = fEvent->GetKink(index); + kink->SetTPCDensity2(-1,1,0); + kink->SetTPCDensity2(-1,1,1); + // + Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2))); + if (row0<15) row0=15; + // + Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2))); + if (row1>145) row1=145; + // + Int_t found,foundable,shared; + seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE); + if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),1,0); + seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE); + if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),1,1); + } + +} + + +Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink) +{ + // + // check kink point for given track + // if return value=0 kink point not found + // otherwise seed0 correspond to mother particle + // seed1 correspond to daughter particle + // kink parameter of kink point + + Int_t middlerow = (seed->fFirstPoint+seed->fLastPoint)/2; + Int_t first = seed->fFirstPoint; + Int_t last = seed->fLastPoint; + if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows + + + AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber + if (!seed1) return 0; + FollowProlongation(*seed1,seed->fLastPoint-20); + seed1->Reset(kTRUE); + FollowProlongation(*seed1,158); + seed1->Reset(kTRUE); + last = seed1->fLastPoint; + // + AliTPCseed *seed0 = new AliTPCseed(*seed); + seed0->Reset(kFALSE); + seed0->Reset(); + // + AliTPCseed param0[20]; // parameters along the track + AliTPCseed param1[20]; // parameters along the track + AliESDkink kinks[20]; // corresponding kink parameters + Int_t rows[20]; + for (Int_t irow=0; irow<20;irow++){ + rows[irow] = first +((last-first)*irow)/19; + } + // store parameters along the track + // + for (Int_t irow=0;irow<20;irow++){ + FollowBackProlongation(*seed0, rows[irow]); + FollowProlongation(*seed1,rows[19-irow]); + new(¶m0[irow]) AliTPCseed(*seed0); + new(¶m1[19-irow]) AliTPCseed(*seed1); + } + // + // define kinks + for (Int_t irow=0; irow<19;irow++){ + kinks[irow].SetMother(param0[irow]); + kinks[irow].SetDaughter(param1[irow]); + kinks[irow].Update(); + } + // + // choose kink with biggest change of angle + Int_t index =-1; + Double_t maxchange= 0; + for (Int_t irow=1;irow<19;irow++){ + if (TMath::Abs(kinks[irow].GetR())>240.) continue; + if (TMath::Abs(kinks[irow].GetR())<110.) continue; + Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].fX)); + if ( quality > maxchange){ + maxchange = quality; + index = irow; + // + } + } + delete seed0; + delete seed1; + if (index<0) return 0; + // + Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate + seed0 = new AliTPCseed(param0[index]); + seed1 = new AliTPCseed(param1[index]); + seed0->Reset(kFALSE); + seed1->Reset(kFALSE); + seed0->ResetCovariance(); + seed1->ResetCovariance(); + FollowProlongation(*seed0,0); + FollowBackProlongation(*seed1,158); + new (&mother) AliTPCseed(*seed0); // backup mother at position 0 + seed0->Reset(kFALSE); + seed1->Reset(kFALSE); + seed0->ResetCovariance(); + seed1->ResetCovariance(); + // + first = TMath::Max(row0-20,0); + last = TMath::Min(row0+20,158); + // + for (Int_t irow=0; irow<20;irow++){ + rows[irow] = first +((last-first)*irow)/19; + } + // store parameters along the track + // + for (Int_t irow=0;irow<20;irow++){ + FollowBackProlongation(*seed0, rows[irow]); + FollowProlongation(*seed1,rows[19-irow]); + new(¶m0[irow]) AliTPCseed(*seed0); + new(¶m1[19-irow]) AliTPCseed(*seed1); + } + // + // define kinks + for (Int_t irow=0; irow<19;irow++){ + kinks[irow].SetMother(param0[irow]); + kinks[irow].SetDaughter(param1[irow]); + // param0[irow].Dump(); + //param1[irow].Dump(); + kinks[irow].Update(); + } + // + // choose kink with biggest change of angle + index =-1; + maxchange= 0; + for (Int_t irow=0;irow<20;irow++){ + if (TMath::Abs(kinks[irow].GetR())>250.) continue; + if (TMath::Abs(kinks[irow].GetR())<90.) continue; + Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].fX)); + if ( quality > maxchange){ + maxchange = quality; + index = irow; + // + } + } + // + // + if (index==-1 || param0[index].fN+param1[index].fN<100){ + delete seed0; + delete seed1; + return 0; + } + // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33); + + kink.SetMother(param0[index]); + kink.SetDaughter(param1[index]); + kink.Update(); + row0 = GetRowNumber(kink.GetR()); + kink.SetTPCRow0(row0); + kink.SetLabel(CookLabel(seed0,0.5,0,row0),0); + kink.SetLabel(CookLabel(seed1,0.5,row0,158),1); + kink.SetIndex(-10,0); + kink.SetIndex(int(param0[index].fN+param1[index].fN),1); + kink.SetTPCncls(param0[index].fN,0); + kink.SetTPCncls(param1[index].fN,1); + // + // + // new (&mother) AliTPCseed(param0[index]); + new (&daughter) AliTPCseed(param1[index]); + daughter.SetLabel(kink.GetLabel(1)); + param0[index].Reset(kTRUE); + FollowProlongation(param0[index],0); + new (&mother) AliTPCseed(param0[index]); + mother.SetLabel(kink.GetLabel(0)); + delete seed0; + delete seed1; + // + return 1; +} + + + + +AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t) +{ + // + // reseed - refit - track + // + Int_t first = 0; + // Int_t last = fSectors->GetNRows()-1; + // + if (fSectors == fOuterSec){ + first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows()); + //last = + } + else + first = t->fFirstPoint; + // + AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9); + FollowBackProlongation(*t,fSectors->GetNRows()-1); + t->Reset(kFALSE); + FollowProlongation(*t,first); + return seed; +} + + + + + + + +//_____________________________________________________________________________ Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) { //----------------------------------------------------------------- // This function reades track seeds. @@ -3909,7 +5402,7 @@ Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) { Int_t n=(Int_t)seedTree->GetEntries(); for (Int_t i=0; iGetEvent(i); - fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha())); + fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/)); } delete seed; @@ -3921,6 +5414,7 @@ Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) { Int_t AliTPCtrackerMI::Clusters2Tracks (AliESD *esd) { // + if (fSeeds) DeleteSeeds(); fEvent = esd; Clusters2Tracks(); if (!fSeeds) return 1; @@ -3952,7 +5446,8 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { if (nc<20) { delete fSeeds->RemoveAt(i); continue; - } + } + CookLabel(pt,0.1); if (pt->fRemoval==10) { if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7) pt->Desactivate(10); // make track again active @@ -3962,11 +5457,50 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { } } } - RemoveDouble(fSeeds,0.2,0.6,11); - RemoveUsed(fSeeds,0.5,0.5,6); - // + RemoveUsed2(fSeeds,0.85,0.85,0); + if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent); + RemoveUsed2(fSeeds,0.5,0.4,20); + // // +// // refit short tracks +// // Int_t nseed=fSeeds->GetEntriesFast(); +// for (Int_t i=0; iUncheckedAt(i), &t=*pt; +// if (!pt) continue; +// Int_t nc=t.GetNumberOfClusters(); +// if (nc<15) { +// delete fSeeds->RemoveAt(i); +// continue; +// } +// if (pt->GetKinkIndexes()[0]!=0) continue; // ignore kink candidates +// if (nc>100) continue; // hopefully, enough for ITS +// AliTPCseed *seed2 = new AliTPCseed(*pt); +// //seed2->Reset(kFALSE); +// //pt->ResetCovariance(); +// seed2->Modify(1); +// FollowBackProlongation(*seed2,158); +// //seed2->Reset(kFALSE); +// seed2->Modify(10); +// FollowProlongation(*seed2,0); +// TTreeSRedirector &cstream = *fDebugStreamer; +// cstream<<"Crefit"<< +// "Tr0.="<fN>pt->fN){ +// delete fSeeds->RemoveAt(i); +// fSeeds->AddAt(seed2,i); +// }else{ +// delete seed2; +// } +// } +// RemoveUsed2(fSeeds,0.6,0.6,50); + +// FindV0s(fSeeds,fEvent); + //RemoveDouble(fSeeds,0.2,0.6,11); + + // Int_t found = 0; for (Int_t i=0; iUncheckedAt(i), &t=*pt; @@ -3979,7 +5513,8 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { CookLabel(pt,0.1); //For comparison only //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){ if ((pt->IsActive() || (pt->fRemoval==10) )){ - cerr<0) cerr<fLab2 = i; } else @@ -4122,6 +5657,7 @@ TObjArray * AliTPCtrackerMI::Tracking() { // // + if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial(); TStopwatch timer; timer.Start(); Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows(); @@ -4245,12 +5781,12 @@ TObjArray * AliTPCtrackerMI::Tracking() cuts[0]=0.0080; // find secondaries - for (Int_t delta = 30; delta<70; delta+=10){ + for (Int_t delta = 30; delta<90; delta+=10){ // cuts[0] = 0.3; - cuts[1] = 1.5; + cuts[1] = 3.5; cuts[2] = 3.; - cuts[3] = 1.5; + cuts[3] = 3.5; arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1); SumTracks(seeds,arr); SignClusters(seeds,fnumber,fdensity); @@ -4272,6 +5808,49 @@ TObjArray * AliTPCtrackerMI::Tracking() } +TObjArray * AliTPCtrackerMI::TrackingSpecial() +{ + // + // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle + // no primary vertex seeding tried + // + TStopwatch timer; + timer.Start(); + Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows(); + + TObjArray * seeds = new TObjArray; + TObjArray * arr=0; + + Int_t gap = 15; + Float_t cuts[4]; + Float_t fnumber = 3.0; + Float_t fdensity = 3.0; + + // find secondaries + cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature + cuts[1] = 3.5; // max tan(phi) angle for seeding + cuts[2] = 3.; // not used (cut on z primary vertex) + cuts[3] = 3.5; // max tan(theta) angle for seeding + + for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){ + // + arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + } + + if (fDebug>0){ + Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast()); + timer.Print(); + timer.Start(); + } + + return seeds; + // + +} + + void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2) const { // @@ -4281,7 +5860,18 @@ void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2) const for (Int_t i=0;iUncheckedAt(i); if (pt){ - + // + // remove tracks with too big curvature + // + if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){ + delete arr2->RemoveAt(i); + continue; + } + // REMOVE VERY SHORT TRACKS + if (pt->GetNumberOfClusters()<20){ + delete arr2->RemoveAt(i); + continue; + }// patch 28 fev06 // NORMAL ACTIVE TRACK if (pt->IsActive()){ arr1->AddLast(arr2->RemoveAt(i)); @@ -4292,11 +5882,7 @@ void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2) const delete arr2->RemoveAt(i); continue; } - // REMOVE VERY SHORT TRACKS - if (pt->GetNumberOfClusters()<20){ - delete arr2->RemoveAt(i); - continue; - } + // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7) arr1->AddLast(arr2->RemoveAt(i)); @@ -4339,6 +5925,7 @@ void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rla for (Int_t i=0; iUncheckedAt(i), &t=*pt; if (!pt) continue; + if (nr==80) pt->UpdateReference(); if (!pt->IsActive()) continue; // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())fRelativeSector>17) { @@ -4423,21 +6010,24 @@ Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) Int_t nseed= arr->GetEntriesFast(); for (Int_t i=0;iUncheckedAt(i); - if (pt) { - AliTPCseed *pt2 = new AliTPCseed(*pt); + if (pt&& pt->GetKinkIndex(0)<=0) { + //AliTPCseed *pt2 = new AliTPCseed(*pt); fSectors = fInnerSec; //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1); //fSectors = fOuterSec; - FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); - - if (fDebug>1 && pt->GetNumberOfClusters()<20 && pt->GetLabel()>0 ){ - Error("PropagateBack","Not prolonged track %d",pt->GetLabel()); - fSectors = fInnerSec; - //FollowBackProlongation(*pt2,fInnerSec->GetNRows()-1); - //fSectors = fOuterSec; - FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); - } - } + FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); + //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){ + // Error("PropagateBack","Not prolonged track %d",pt->GetLabel()); + // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); + //} + } + if (pt&& pt->GetKinkIndex(0)>0) { + AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1); + pt->fFirstPoint = kink->GetTPCRow0(); + fSectors = fInnerSec; + FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); + } + } return 0; } @@ -4449,16 +6039,12 @@ Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr) // make forward propagation // Int_t nseed= arr->GetEntriesFast(); + // for (Int_t i=0;iUncheckedAt(i); if (pt) { - AliTPCseed *pt2 = new AliTPCseed(*pt); FollowProlongation(*pt,0); - if (pt->GetNumberOfClusters()<35 && pt->GetLabel()>0 ){ - FollowProlongation(*pt2,0); - } - delete pt2; - } + } } return 0; } @@ -4468,7 +6054,7 @@ Int_t AliTPCtrackerMI::PropagateForward() { // // propagate track forward - UnsignClusters(); + //UnsignClusters(); Int_t nseed = fSeeds->GetEntriesFast(); for (Int_t i=0;iUncheckedAt(i); @@ -4535,7 +6121,7 @@ void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row) Float_t padlength = GetPadPitchLength(row); // Float_t sresy = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3; - Float_t angulary = seed->GetSnp(); + Double_t angulary = seed->GetSnp(); angulary = angulary*angulary/(1-angulary*angulary); seed->fCurrentSigmaY2 = sd2+padlength*padlength*angulary/12.+sresy*sresy; // @@ -4550,41 +6136,123 @@ void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row) if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){ printf("problem\n"); } - */ -} + */ +} + + +Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed) +{ + // + // + Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); + Float_t padlength = fParam->GetPadPitchLength(seed->fSector); + Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3; + Float_t angular = seed->GetSnp(); + angular = angular*angular/(1-angular*angular); + // angular*=angular; + //angular = TMath::Sqrt(angular/(1-angular)); + Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres); + return res; +} +Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed) +{ + // + // + Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); + Float_t padlength = fParam->GetPadPitchLength(seed->fSector); + Float_t sres = fParam->GetZSigma(); + Float_t angular = seed->GetTgl(); + Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres); + return res; +} + + +//__________________________________________________________________________ +void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const { + //-------------------------------------------------------------------- + //This function "cooks" a track label. If label<0, this track is fake. + //-------------------------------------------------------------------- + AliTPCseed * t = (AliTPCseed*)tk; + Int_t noc=t->GetNumberOfClusters(); + if (noc<10){ + //printf("\nnot founded prolongation\n\n\n"); + //t->Dump(); + return ; + } + Int_t lb[160]; + Int_t mx[160]; + AliTPCclusterMI *clusters[160]; + // + for (Int_t i=0;i<160;i++) { + clusters[i]=0; + lb[i]=mx[i]=0; + } + + Int_t i; + Int_t current=0; + for (i=0; i<160 && currentGetClusterIndex2(i); + if (index<=0) continue; + if (index&0x8000) continue; + // + //clusters[current]=GetClusterMI(index); + if (t->fClusterPointer[i]){ + clusters[current]=t->fClusterPointer[i]; + current++; + } + } + noc = current; + + Int_t lab=123456789; + for (i=0; iGetLabel(0)); + Int_t j; + for (j=0; jmax) {max=mx[i]; lab=lb[i];} + + for (i=0; iGetLabel(1)) == lab || + TMath::Abs(c->GetLabel(2)) == lab ) max++; + } -Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed) -{ - // - // - Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); - Float_t padlength = fParam->GetPadPitchLength(seed->fSector); - Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3; - Float_t angular = seed->GetSnp(); - angular = angular*angular/(1-angular*angular); - // angular*=angular; - //angular = TMath::Sqrt(angular/(1-angular)); - Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres); - return res; -} -Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed) -{ - // - // - Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); - Float_t padlength = fParam->GetPadPitchLength(seed->fSector); - Float_t sres = fParam->GetZSigma(); - Float_t angular = seed->GetTgl(); - Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres); - return res; -} + if ((1.- Float_t(max)/noc) > wrong) lab=-lab; + + else { + Int_t tail=Int_t(0.10*noc); + max=0; + Int_t ind=0; + for (i=1; i<=160&&indGetLabel(0)) || + lab == TMath::Abs(c->GetLabel(1)) || + lab == TMath::Abs(c->GetLabel(2))) max++; + ind++; + } + if (max < Int_t(0.5*tail)) lab=-lab; + } + t->SetLabel(lab); + // delete[] lb; + //delete[] mx; + //delete[] clusters; +} //__________________________________________________________________________ -void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const { +Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const { //-------------------------------------------------------------------- //This function "cooks" a track label. If label<0, this track is fake. //-------------------------------------------------------------------- @@ -4592,7 +6260,7 @@ void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const { if (noc<10){ //printf("\nnot founded prolongation\n\n\n"); //t->Dump(); - return ; + return -1; } Int_t lb[160]; Int_t mx[160]; @@ -4606,7 +6274,8 @@ void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const { Int_t i; Int_t current=0; for (i=0; i<160 && currentlast) continue; Int_t index=t->GetClusterIndex2(i); if (index<=0) continue; if (index&0x8000) continue; @@ -4618,7 +6287,7 @@ void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const { } } noc = current; - + if (noc<5) return -1; Int_t lab=123456789; for (i=0; iSetLabel(lab); - + // t->SetLabel(lab); + return lab; // delete[] lb; //delete[] mx; //delete[] clusters; @@ -4689,6 +6358,18 @@ Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const } } +Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const +{ + //return pad row number for given x vector + Float_t phi = TMath::ATan2(x[1],x[0]); + if(phi<0) phi=2.*TMath::Pi()+phi; + // Get the local angle in the sector philoc + const Float_t kRaddeg = 180/3.14159265358979312; + Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg; + Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle); + return GetRowNumber(localx); +} + //_________________________________________________________________________ void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) { //----------------------------------------------------------------------- @@ -4815,10 +6496,16 @@ AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t // Return the index of the nearest cluster in z y //----------------------------------------------------------------------- Float_t maxdistance = roady*roady + roadz*roadz; - Int_t iz1 = TMath::Max(fFastCluster[Int_t(z-roadz+254.5)]-1,0); - Int_t iz2 = TMath::Min(fFastCluster[Int_t(z+roadz+255.5)]+1,fN); - AliTPCclusterMI *cl =0; + + //PH Check boundaries. 510 is the size of fFastCluster + Int_t iz1 = Int_t(z-roadz+254.5); + if (iz1<0 || iz1>=510) return cl; + iz1 = TMath::Max(fFastCluster[iz1]-1,0); + Int_t iz2 = Int_t(z+roadz+255.5); + if (iz2<0 || iz2>=510) return cl; + iz2 = TMath::Min(fFastCluster[iz2]+1,fN); + //FindNearest3(y,z,roady,roadz,index); // for (Int_t i=Find(z-roadz); i0) { - SetClusterIndex2(i,index); - } - else{ - SetClusterIndex2(i,-3); - } - } - fFirstPoint =0; - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; - -} - -AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){ - // - //copy constructor - fRow=0; - for (Int_t i=0;i<160;i++) { - fClusterPointer[i] = 0; - Int_t index = t.GetClusterIndex(i); - SetClusterIndex2(i,index); - } - - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared =0; - // fTrackPoints =0; - fRemoval =0; - fSort = 0; - fFirstPoint =0; - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; -} - -AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15], - Double_t xr, Double_t alpha): - AliTPCtrack(index, xx, cc, xr, alpha) { - // - // - //constructor - fRow =0; - for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); - for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared = 0; - // fTrackPoints =0; - fRemoval =0; - fSort =0; - fFirstPoint =0; - // fHelixIn = new TClonesArray("AliHelix",0); - //fHelixOut = new TClonesArray("AliHelix",0); - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; - -} - -AliTPCseed::~AliTPCseed(){ - // - // destructor - if (fPoints) delete fPoints; - fPoints =0; - if (fEPoints) delete fEPoints; - fEPoints = 0; - fNoCluster =0; -} - -AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i) -{ - // - // - return &fTrackPoints[i]; -} - -void AliTPCseed::RebuildSeed() -{ - // - // rebuild seed to be ready for storing - AliTPCclusterMI cldummy; - cldummy.SetQ(0); - AliTPCTrackPoint pdummy; - pdummy.GetTPoint().fIsShared = 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; - } - - } - -} - - -Double_t AliTPCseed::GetDensityFirst(Int_t n) -{ - // - // - // return cluster for n rows bellow first point - Int_t nfoundable = 1; - Int_t nfound = 1; - for (Int_t i=fLastPoint-1;i>0&&nfoundable0) nfound++; - } - if (nfoundableIsUsed(10)) { - shared++; - continue; - } - if (!plus2) continue; //take also neighborhoud - // - if ( (i>0) && fClusterPointer[i-1]){ - if (fClusterPointer[i-1]->IsUsed(10)) { - shared++; - continue; - } - } - if ( fClusterPointer[i+1]){ - if (fClusterPointer[i+1]->IsUsed(10)) { - shared++; - continue; - } - } - - } - //if (shared>found){ - //Error("AliTPCseed::GetClusterStatistic","problem\n"); - //} -} - -//_____________________________________________________________________________ -void 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.; - - 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 = (1.05-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.; - - 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]++; - } - - 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)); - - - // 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; - } - // 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); - - //mi deDX - - - - //Very rough PID - Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt())); - - if (p<0.6) { - if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;} - if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;} - SetMass(0.93827); return; - } - - if (p<1.2) { - if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;} - SetMass(0.93827); return; - } - - SetMass(0.13957); return; - -} - - - -/* - - - -void AliTPCseed::CookdEdx2(Double_t low, Double_t up) { - //----------------------------------------------------------------- - // This funtion calculates dE/dX within the "low" and "up" cuts. - //----------------------------------------------------------------- - - Float_t amp[200]; - Float_t angular[200]; - Float_t weight[200]; - Int_t index[200]; - Bool_t inlimit[200]; - for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE; - for (Int_t i=0;i<200;i++) amp[i]=10000; - for (Int_t i=0;i<200;i++) angular[i]= 1;; - - - // - Float_t meanlog = 100.; - Int_t indexde[4]={0,64,128,160}; - - Float_t amean =0; - Float_t asigma =0; - Float_t anc =0; - Float_t anorm =0; - - Float_t mean[4] = {0,0,0,0}; - Float_t sigma[4] = {1000,1000,1000,1000}; - Int_t nc[4] = {0,0,0,0}; - Float_t norm[4] = {1000,1000,1000,1000}; - // - // - fNShared =0; - - // for (Int_t of =0; of<3; of++){ - // for (Int_t i=indexde[of];ifIsShared){ - fNShared++; - continue; - } - Int_t type = point->GetCPoint().GetType(); - if (type<0) continue; - if (point->GetCPoint().GetMax()<5) continue; - Float_t angley = point->GetTPoint().GetAngleY(); - Float_t anglez = point->GetTPoint().GetAngleZ(); - Float_t rsigmay = point->GetCPoint().GetSigmaY(); - Float_t rsigmaz = point->GetCPoint().GetSigmaZ(); - Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz); - - Float_t ampc = 0; // normalization to the number of electrons - if (i>64){ - ampc = point->GetCPoint().GetMax(); - } - else{ - ampc = point->GetCPoint().GetMax(); - } - ampc *= 2.0; // put mean value to channel 50 - // ampc *= 0.565; // put mean value to channel 50 - - Float_t w = 1.; - Float_t z = TMath::Abs(point->GetCPoint().GetZ()); - if (i<64) { - ampc /= 0.63; - } else - if (i>128){ - ampc /=1.51; - } - if (type<0) { //amp at the border - lower weight - continue; - } - if (rsigma>1.5) ampc/=1.3; // if big backround - angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez); - amp[i] = ampc/angular[i]; - weight[i] = w; - anc++; - } - - TMath::Sort(159,amp,index,kFALSE); - for (Int_t i=int(anc*low+0.5);i0.1) - sigma[of] = TMath::Sqrt(sigma[of]); - else - sigma[of] = 1000; - mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; - } - } - - Float_t dedx =0; - fSdEdx =0; - fMAngular =0; - // - Int_t norm2 = 0; - Int_t norm3 = 0; - Float_t www[3] = {12.,14.,17.}; - //Float_t www[3] = {1.,1.,1.}; - - for (Int_t i =0;i<3;i++){ - if (nc[i]>2&&nc[i]<1000){ - dedx += mean[i] *nc[i]*www[i]/sigma[i]; - fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i]; - fMAngular += norm[i] *nc[i]; - norm2 += nc[i]*www[i]/sigma[i]; - norm3 += (nc[i]-2)*www[i]/sigma[i]; - } - fDEDX[i] = mean[i]; - fSDEDX[i] = sigma[i]; - fNCDEDX[i]= nc[i]; - } - - if (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; - } - else{ - SetdEdx(0); - return; - } - // Float_t dedx1 =dedx; - - dedx =0; - Float_t norm4 = 0; - for (Int_t i =0;i<3;i++){ - if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){ - //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.)); - dedx += mean[i] *(nc[i])/(sigma[i]); - norm4 += (nc[i])/(sigma[i]); - } - fDEDX[i] = mean[i]; - } - if (norm4>0) dedx /= norm4; - - - SetdEdx(dedx); - - //mi deDX +// AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i) +// { +// // +// // +// return &fTrackPoints[i]; +// } -} -*/