X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCtrackerMI.cxx;h=d62e100f2abcd3f179f637664779a974cfc0856a;hb=95fda180b699f8781d9ba179e79f1d3c3716c139;hp=86e022cfa7b59c1700c31b8a3f4848172f5bab59;hpb=f8aae3779b58ff58faf59c4c1104d068068b06f6;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCtrackerMI.cxx b/TPC/AliTPCtrackerMI.cxx index 86e022cfa7b..d62e100f2ab 100644 --- a/TPC/AliTPCtrackerMI.cxx +++ b/TPC/AliTPCtrackerMI.cxx @@ -13,215 +13,316 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ - -/* - 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 -*/ //------------------------------------------------------- // Implementation of the TPC tracker // // Origin: Marian Ivanov Marian.Ivanov@cern.ch // +// AliTPC parallel tracker //------------------------------------------------------- -#include + +/* $Id$ */ + +#include "Riostream.h" +#include #include +#include #include -#include "Riostream.h" +#include "AliLog.h" -#include "AliTPCtrackerMI.h" -#include "AliTPCclusterMI.h" -#include "AliTPCParam.h" -#include "AliTPCClustersRow.h" #include "AliComplexCluster.h" -#include "AliTPCpolyTrack.h" +#include "AliESD.h" +#include "AliKink.h" +#include "AliV0.h" +#include "AliHelix.h" #include "AliRunLoader.h" +#include "AliTPCClustersRow.h" +#include "AliTPCParam.h" +#include "AliTPCReconstructor.h" +#include "AliTPCpolyTrack.h" +#include "AliTPCreco.h" +#include "AliTPCseed.h" +#include "AliTPCtrackerMI.h" #include "TStopwatch.h" +#include "AliTPCReconstructor.h" +#include "AliPID.h" +#include "TTreeStream.h" +#include "AliAlignObj.h" +#include "AliTrackPointArray.h" +// -ClassImp(AliTPCseed) -ClassImp(AliTPCKalmanSegment) - - - -//_____________________________________________________________________________ - -AliTPCKalmanSegment::AliTPCKalmanSegment(){ - // - // - fX=fAlpha=fChi2=0; - for (Int_t i=0;i<5;i++) fState[i] = 0.; - for (Int_t i=0;i<15;i++) fCovariance[i] = 0.; - fNCFoundable = 0; - fNC = 0; - // fN = 0; -} +ClassImp(AliTPCtrackerMI) -//_____________________________________________________________________________ -void AliTPCKalmanSegment::Init(AliTPCseed* seed) -{ - // in initialization - // initial entrance integral chi2, fNCFoundable and fNC stored - fNCFoundable = seed->fNFoundable; - fNC = seed->GetNumberOfClusters(); - fChi2 = seed->GetChi2(); -} +class AliTPCFastMath { +public: + AliTPCFastMath(); + static Double_t FastAsin(Double_t x); + private: + static Double_t fgFastAsin[20000]; //lookup table for fast asin computation +}; +Double_t AliTPCFastMath::fgFastAsin[20000]; +AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT -void AliTPCKalmanSegment::Finish(AliTPCseed* seed) -{ - // - // in finish state vector stored and chi2 and fNC... calculated - Double_t x; - Double_t state[5]; - Double_t cov[15]; - seed->GetExternalParameters(x,state); - seed->GetExternalCovariance(cov); - //float precision for tracklet - for (Int_t i=0;i<5;i++) fState[i] = state[i]; - for (Int_t i=0;i<15;i++) fCovariance[i] = cov[i]; - // - // in current seed integral track characteristic - // for tracklet differenciation between beginning (Init state) and Finish state - fNCFoundable = seed->fNFoundable - fNCFoundable; - fNC = seed->GetNumberOfClusters() - fNC; - fChi2 = seed->GetChi2()-fChi2; +AliTPCFastMath::AliTPCFastMath(){ // + // initialized lookup table; + for (Int_t i=0;i<10000;i++){ + fgFastAsin[2*i] = TMath::ASin(i/10000.); + fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]); + } } -void AliTPCKalmanSegment::GetState(Double_t &x, Double_t & alpha, Double_t state[5]) -{ +Double_t AliTPCFastMath::FastAsin(Double_t x){ // - x = fX; - alpha = fAlpha; - for (Int_t i=0;i<5;i++) state[i] = fState[i]; -} - -void AliTPCKalmanSegment::GetCovariance(Double_t covariance[15]) -{ - // - for (Int_t i=0;i<5;i++) covariance[i] = fCovariance[i]; + // return asin using lookup table + if (x>0){ + Int_t index = int(x*10000); + return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]; + } + x*=-1; + Int_t index = int(x*10000); + return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]); } - -void AliTPCKalmanSegment::GetStatistic(Int_t & nclusters, Int_t & nfoundable, Float_t & chi2) +//__________________________________________________________________ +AliTPCtrackerMI::AliTPCtrackerMI() + :AliTracker(), + fkNIS(0), + fInnerSec(0), + fkNOS(0), + fOuterSec(0), + fN(0), + fSectors(0), + fInput(0), + fOutput(0), + fSeedTree(0), + fTreeDebug(0), + fEvent(0), + fDebug(0), + fNewIO(kFALSE), + fNtracks(0), + fSeeds(0), + fIteration(0), + fParam(0), + fDebugStreamer(0) { // + // default constructor // - nclusters = fNC; - nfoundable = fNCFoundable; - chi2 = fChi2; -} - - - -AliTPCclusterTracks::AliTPCclusterTracks(){ - // class for storing overlaping info - fTrackIndex[0]=-1; - fTrackIndex[1]=-1; - fTrackIndex[2]=-1; - fDistance[0]=1000; - fDistance[1]=1000; - fDistance[2]=1000; } +//_____________________________________________________________________ +Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){ + // + //update track information using current cluster - track->fCurrentCluster + AliTPCclusterMI* c =track->GetCurrentCluster(); + if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters - - - -Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){ + UInt_t i = track->GetCurrentClusterIndex1(); Int_t sec=(i&0xff000000)>>24; - Int_t row = (i&0x00ff0000)>>16; - track->fRow=(i&0x00ff0000)>>16; - track->fSector = sec; + //Int_t row = (i&0x00ff0000)>>16; + track->SetRow((i&0x00ff0000)>>16); + track->SetSector(sec); // Int_t index = i&0xFFFF; - if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow(); - track->fClusterIndex[track->fRow] = i; - track->fFirstPoint = row; - if ( track->fLastPointfLastPoint =row; + if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow()); + track->SetClusterIndex2(track->GetRow(), i); + //track->fFirstPoint = row; + //if ( track->fLastPointfLastPoint =row; + // if (track->fRow<0 || track->fRow>160) { + // printf("problem\n"); + //} + if (track->GetFirstPoint()>track->GetRow()) + track->SetFirstPoint(track->GetRow()); + if (track->GetLastPoint()GetRow()) + track->SetLastPoint(track->GetRow()); + + + track->SetClusterPointer(track->GetRow(),c); // - AliTPCTrackPoint *trpoint =track->GetTrackPoint(track->fRow); - 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 (c!=0){ - //if we have a cluster - trpoint->GetCPoint().SetY(c->GetY()); - trpoint->GetCPoint().SetZ(c->GetZ()); - // - trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY)); - trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ)); + if (angle2<1) //PH sometimes angle2 is very big. To be investigated... + { + angle2 = TMath::Sqrt(angle2/(1-angle2)); + AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow())); // - trpoint->GetCPoint().SetType(c->GetType()); - trpoint->GetCPoint().SetQ(c->GetQ()); - trpoint->GetCPoint().SetMax(c->GetMax()); - // - trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2)); - trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2)); + point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2()); + point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2()); + point.SetErrY(sqrt(track->GetErrorY2())); + point.SetErrZ(sqrt(track->GetErrorZ2())); // + point.SetX(track->GetX()); + point.SetY(track->GetY()); + point.SetZ(track->GetZ()); + point.SetAngleY(angle2); + point.SetAngleZ(track->GetTgl()); + if (point.IsShared()){ + track->SetErrorY2(track->GetErrorY2()*4); + track->SetErrorZ2(track->GetErrorZ2()*4); + } + } + + Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster()); + // + track->SetErrorY2(track->GetErrorY2()*1.3); + track->SetErrorY2(track->GetErrorY2()+0.01); + track->SetErrorZ2(track->GetErrorZ2()*1.3); + track->SetErrorZ2(track->GetErrorZ2()+0.005); + //} + if (accept>0) return 0; + if (track->GetNumberOfClusters()%20==0){ + // if (track->fHelixIn){ + // TClonesArray & larr = *(track->fHelixIn); + // Int_t ihelix = larr.GetEntriesFast(); + // new(larr[ihelix]) AliHelix(*track) ; + //} } - trpoint->GetTPoint().SetX(track->GetX()); - trpoint->GetTPoint().SetY(track->GetY()); - trpoint->GetTPoint().SetZ(track->GetZ()); + track->SetNoCluster(0); + return track->Update(c,chi2,i); +} + + + +Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor, + Float_t cory, Float_t corz) +{ // - trpoint->GetTPoint().SetAngleY(angle2); - trpoint->GetTPoint().SetAngleZ(track->GetTgl()); + // decide according desired precision to accept given + // cluster for tracking + Double_t sy2=ErrY2(seed,cluster)*cory; + Double_t sz2=ErrZ2(seed,cluster)*corz; + //sy2=ErrY2(seed,cluster)*cory; + //sz2=ErrZ2(seed,cluster)*cory; - + Double_t sdistancey2 = sy2+seed->GetSigmaY2(); + Double_t sdistancez2 = sz2+seed->GetSigmaZ2(); + + Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())* + (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2; + Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())* + (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2; + + Double_t rdistance2 = rdistancey2+rdistancez2; + //Int_t accept =0; + + if (rdistance2>16) return 3; + + + if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0) + return 2; //suspisiouce - will be changed + + if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0) + // strict cut on overlaped cluster + return 2; //suspisiouce - will be changed - if (chi2>10){ - // printf("suspicious chi2 %f\n",chi2); + if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor ) + && cluster->GetType()<0){ + seed->SetNFoundable(seed->GetNFoundable()-1); + return 2; } - // if (track->fIsSeeding){ - track->fErrorY2 *= 1.2; - track->fErrorY2 += 0.0064; - track->fErrorZ2 *= 1.2; - track->fErrorY2 += 0.005; - - //} + return 0; +} + + - return track->Update(c,chi2,i); -} //_____________________________________________________________________________ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par): -AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2) +AliTracker(), + fkNIS(par->GetNInnerSector()/2), + fInnerSec(0), + fkNOS(par->GetNOuterSector()/2), + fOuterSec(0), + fN(0), + fSectors(0), + fInput(0), + fOutput(0), + fSeedTree(0), + fTreeDebug(0), + fEvent(0), + fDebug(0), + fNewIO(0), + fNtracks(0), + fSeeds(0), + fIteration(0), + fParam(0), + fDebugStreamer(0) { //--------------------------------------------------------------------- // The main TPC tracker constructor //--------------------------------------------------------------------- fInnerSec=new AliTPCSector[fkNIS]; fOuterSec=new AliTPCSector[fkNOS]; - + Int_t i; for (i=0; iGetNRowLow(); + Int_t nrowup = par->GetNRowUp(); + + + for (Int_t i=0;iGetPadRowRadiiLow(i); + fPadLength[i]= par->GetPadPitchLength(0,i); + fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle()); + } - fClustersArray.Setup(par); - fClustersArray.SetClusterType("AliTPCclusterMI"); + + for (Int_t i=0;iGetPadRowRadiiUp(i); + fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i); + fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle()); + } - fSeeds=0; - fNtracks = 0; - fParam = par; + fDebugStreamer = new TTreeSRedirector("TPCdebug.root"); +} +//________________________________________________________________________ +AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): + AliTracker(t), + fkNIS(t.fkNIS), + fInnerSec(0), + fkNOS(t.fkNOS), + fOuterSec(0), + fN(0), + fSectors(0), + fInput(0), + fOutput(0), + fSeedTree(0), + fTreeDebug(0), + fEvent(0), + fDebug(0), + fNewIO(kFALSE), + fNtracks(0), + fSeeds(0), + fIteration(0), + fParam(0), + fDebugStreamer(0) +{ + //------------------------------------ + // dummy copy constructor + //------------------------------------------------------------------ + fOutput=t.fOutput; +} +AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){ + //------------------------------ + // dummy + //-------------------------------------------------------------- + return *this; } - //_____________________________________________________________________________ AliTPCtrackerMI::~AliTPCtrackerMI() { //------------------------------------------------------------------ @@ -233,228 +334,632 @@ AliTPCtrackerMI::~AliTPCtrackerMI() { fSeeds->Delete(); delete fSeeds; } + if (fDebugStreamer) delete fDebugStreamer; } - -Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ - // +void AliTPCtrackerMI::SetIO() +{ // - Float_t snoise2; - Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ())); + fNewIO = kTRUE; + fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName()); + + fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName()); + if (fOutput){ + AliTPCtrack *iotrack= new AliTPCtrack; + fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100); + delete iotrack; + } +} - //cluster "quality" - Float_t rsigmay = 1; - Int_t ctype = 0; - //standard if we don't have cluster - take MIP - const Float_t chmip = 50.; - Float_t amp = chmip/0.3; - Float_t nel; - Float_t nprim; - if (cl){ - amp = cl->GetQ(); - rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY); - ctype = cl->GetType(); +void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event) +{ + + // set input + fNewIO = kFALSE; + fInput = 0; + fOutput = 0; + fSeedTree = 0; + fTreeDebug =0; + fInput = input; + if (input==0){ + return; + } + //set output + fOutput = output; + if (output){ + AliTPCtrack *iotrack= new AliTPCtrack; + // iotrack->fHelixIn = new TClonesArray("AliHelix"); + //iotrack->fHelixOut = new TClonesArray("AliHelix"); + fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100); + delete iotrack; + } + if (output && (fDebug&2)){ + //write the full seed information if specified in debug mode + // + fSeedTree = new TTree("Seeds","Seeds"); + AliTPCseed * vseed = new AliTPCseed; + // + TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160); + arrtr->ExpandCreateFast(160); + TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160); + // + vseed->SetPoints(arrtr); + vseed->SetEPoints(arre); + // vseed->fClusterPoints = arrcl; + fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99); + delete arrtr; + delete arre; + fTreeDebug = new TTree("trackDebug","trackDebug"); + TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0); + fTreeDebug->Branch("debug",&arrd,32000,99); } - - Float_t landau=2 ; //landau fluctuation part - Float_t gg=2; // gg fluctuation part - Float_t padlength= fSectors->GetPadPitchLength(seed->GetX()); + //set ESD event + fEvent = event; +} - if (fSectors==fInnerSec){ - snoise2 = 0.0004; - nel = 0.268*amp; - nprim = 0.155*amp; - gg = (2+0.0002*amp)/nel; - landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim; - if (landau>1) landau=1; +void AliTPCtrackerMI::FillESD(TObjArray* arr) +{ + // + // + //fill esds using updated tracks + if (fEvent){ + // 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->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->GetNFoundable()))>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->GetNFoundable()))>0.70) { + Int_t found,foundable,shared; + pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE); + if ( (found>20) && (pt->GetNShared()/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->GetNFoundable()))>0.8) { + Int_t found,foundable,shared; + pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE); + if (found<20) continue; + if (pt->GetNShared()/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->GetNShared()/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->GetNShared()/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; + } + } } - else { - snoise2 = 0.0004; - nel = 0.3*amp; - nprim = 0.133*amp; - gg = (2+0.0002*amp)/nel; - landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim; - if (landau>1) landau=1; + printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()); +} + +void AliTPCtrackerMI::WriteTracks(TTree * tree) +{ + // + // write tracks from seed array to selected tree + // + fOutput = tree; + if (fOutput){ + AliTPCtrack *iotrack= new AliTPCtrack; + fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100); } + WriteTracks(); +} +void AliTPCtrackerMI::WriteTracks() +{ + // + // write tracks to the given output tree - + // output specified with SetIO routine + if (!fSeeds) return; + if (!fOutput){ + SetIO(); + } - Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z; - Float_t angle2 = seed->GetSnp()*seed->GetSnp(); - angle2 = angle2/(1-angle2); - Float_t angular = landau*angle2*padlength*padlength/12.; - Float_t res = sdiff + angular; - + if (fOutput){ + AliTPCtrack *iotrack= 0; + Int_t nseed=fSeeds->GetEntriesFast(); + //for (Int_t i=0; iUncheckedAt(i); + // if (iotrack) break; + //} + //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100); + TBranch * br = fOutput->GetBranch("tracks"); + br->SetAddress(&iotrack); + // + for (Int_t i=0; iUncheckedAt(i); + if (!pt) continue; + AliTPCtrack * track = new AliTPCtrack(*pt); + iotrack = track; + pt->SetLab2(i); + // br->SetAddress(&iotrack); + fOutput->Fill(); + delete track; + iotrack =0; + } + //fOutput->GetDirectory()->cd(); + //fOutput->Write(); + } + // delete iotrack; + // + if (fSeedTree){ + //write the full seed information if specified in debug mode + + AliTPCseed * vseed = new AliTPCseed; + // + TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160); + arrtr->ExpandCreateFast(160); + //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160); + //arrcl->ExpandCreateFast(160); + TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160); + // + vseed->SetPoints(arrtr); + vseed->SetEPoints(arre); + // vseed->fClusterPoints = arrcl; + //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99); + TBranch * brseed = fSeedTree->GetBranch("seeds"); + + Int_t nseed=fSeeds->GetEntriesFast(); + + for (Int_t i=0; iUncheckedAt(i); + if (!pt) continue; + pt->SetPoints(arrtr); + // pt->fClusterPoints = arrcl; + pt->SetEPoints(arre); + pt->RebuildSeed(); + vseed = pt; + brseed->SetAddress(&vseed); + fSeedTree->Fill(); + pt->SetPoints(0); + pt->SetEPoints(0); + // pt->fClusterPoints = 0; + } + fSeedTree->Write(); + if (fTreeDebug) fTreeDebug->Write(); + } + +} - if ((ctype==0) && (fSectors ==fOuterSec)) - res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2)); - if ((ctype==0) && (fSectors ==fInnerSec)) - res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2)); - if ((ctype>0)) - res*= TMath::Power((rsigmay+0.5),1.5)+0.0064; - - if (ctype<0) +Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ + // + // + //seed->SetErrorY2(0.1); + //return 0.1; + //calculate look-up table at the beginning + static Bool_t ginit = kFALSE; + static Float_t gnoise1,gnoise2,gnoise3; + static Float_t ggg1[10000]; + static Float_t ggg2[10000]; + static Float_t ggg3[10000]; + static Float_t glandau1[10000]; + static Float_t glandau2[10000]; + static Float_t glandau3[10000]; + // + static Float_t gcor01[500]; + static Float_t gcor02[500]; + static Float_t gcorp[500]; + // + + // + if (ginit==kFALSE){ + for (Int_t i=1;i<500;i++){ + Float_t rsigma = float(i)/100.; + gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6); + gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6); + gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2); + } + + // + for (Int_t i=3;i<10000;i++){ + // + // + // inner sector + Float_t amp = float(i); + Float_t padlength =0.75; + gnoise1 = 0.0004/padlength; + Float_t nel = 0.268*amp; + Float_t nprim = 0.155*amp; + ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel; + glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (glandau1[i]>1) glandau1[i]=1; + glandau1[i]*=padlength*padlength/12.; + // + // outer short + padlength =1.; + gnoise2 = 0.0004/padlength; + nel = 0.3*amp; + nprim = 0.133*amp; + ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; + glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (glandau2[i]>1) glandau2[i]=1; + glandau2[i]*=padlength*padlength/12.; + // + // + // outer long + padlength =1.5; + gnoise3 = 0.0004/padlength; + nel = 0.3*amp; + nprim = 0.133*amp; + ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; + glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (glandau3[i]>1) glandau3[i]=1; + glandau3[i]*=padlength*padlength/12.; + // + } + ginit = kTRUE; + } + // + // + // + Int_t amp = int(TMath::Abs(cl->GetQ())); + if (amp>9999) { + seed->SetErrorY2(1.); + return 1.; + } + Float_t snoise2; + Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ())); + Int_t ctype = cl->GetType(); + Float_t padlength= GetPadPitchLength(seed->GetRow()); + Double_t angle2 = seed->GetSnp()*seed->GetSnp(); + angle2 = angle2/(1-angle2); + // + //cluster "quality" + Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2())); + Float_t res; + // + if (fSectors==fInnerSec){ + snoise2 = gnoise1; + res = ggg1[amp]*z+glandau1[amp]*angle2; + if (ctype==0) res *= gcor01[rsigmay]; + if ((ctype>0)){ + res+=0.002; + res*= gcorp[rsigmay]; + } + } + else { + if (padlength<1.1){ + snoise2 = gnoise2; + res = ggg2[amp]*z+glandau2[amp]*angle2; + if (ctype==0) res *= gcor02[rsigmay]; + if ((ctype>0)){ + res+=0.002; + res*= gcorp[rsigmay]; + } + } + else{ + snoise2 = gnoise3; + res = ggg3[amp]*z+glandau3[amp]*angle2; + if (ctype==0) res *= gcor02[rsigmay]; + if ((ctype>0)){ + res+=0.002; + res*= gcorp[rsigmay]; + } + } + } + + if (ctype<0){ + res+=0.005; res*=2.4; // overestimate error 2 times - + } res+= snoise2; if (res<2*snoise2) res = 2*snoise2; - + seed->SetErrorY2(res); return res; -} +} + Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ // // - Float_t snoise2; - Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ())); - //signal quality - Float_t rsigmaz=1; - Int_t ctype =0; - - const Float_t chmip = 50.; - Float_t amp = chmip/0.3; - Float_t nel; - Float_t nprim; - if (cl){ - amp = cl->GetQ(); - rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ); - ctype = cl->GetType(); - } + //seed->SetErrorY2(0.1); + //return 0.1; + //calculate look-up table at the beginning + static Bool_t ginit = kFALSE; + static Float_t gnoise1,gnoise2,gnoise3; + static Float_t ggg1[10000]; + static Float_t ggg2[10000]; + static Float_t ggg3[10000]; + static Float_t glandau1[10000]; + static Float_t glandau2[10000]; + static Float_t glandau3[10000]; + // + static Float_t gcor01[1000]; + static Float_t gcor02[1000]; + static Float_t gcorp[1000]; + // // - Float_t landau=2 ; //landau fluctuation part - Float_t gg=2; // gg fluctuation part - Float_t padlength= fSectors->GetPadPitchLength(seed->GetX()); + if (ginit==kFALSE){ + for (Int_t i=1;i<1000;i++){ + Float_t rsigma = float(i)/100.; + gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6); + gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6); + gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2); + } - if (fSectors==fInnerSec){ - snoise2 = 0.0004; - nel = 0.268*amp; - nprim = 0.155*amp; - gg = (2+0.0002*amp)/nel; - landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim; - if (landau>1) landau=1; - } - else { - snoise2 = 0.0004; - nel = 0.3*amp; - nprim = 0.133*amp; - gg = (2+0.0002*amp)/nel; - landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim; - if (landau>1) landau=1; + // + for (Int_t i=3;i<10000;i++){ + // + // + // inner sector + Float_t amp = float(i); + Float_t padlength =0.75; + gnoise1 = 0.0004/padlength; + Float_t nel = 0.268*amp; + Float_t nprim = 0.155*amp; + ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel; + glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (glandau1[i]>1) glandau1[i]=1; + glandau1[i]*=padlength*padlength/12.; + // + // outer short + padlength =1.; + gnoise2 = 0.0004/padlength; + nel = 0.3*amp; + nprim = 0.133*amp; + ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; + glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (glandau2[i]>1) glandau2[i]=1; + glandau2[i]*=padlength*padlength/12.; + // + // + // outer long + padlength =1.5; + gnoise3 = 0.0004/padlength; + nel = 0.3*amp; + nprim = 0.133*amp; + ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; + glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (glandau3[i]>1) glandau3[i]=1; + glandau3[i]*=padlength*padlength/12.; + // + } + ginit = kTRUE; } - Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z; - - Float_t angle = seed->GetTgl(); - Float_t angular = landau*angle*angle*padlength*padlength/12.; - Float_t res = sdiff + angular; - - if ((ctype==0) && (fSectors ==fOuterSec)) - res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2)); + // + // + // + Int_t amp = int(TMath::Abs(cl->GetQ())); + if (amp>9999) { + seed->SetErrorY2(1.); + return 1.; + } + Float_t snoise2; + Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ())); + Int_t ctype = cl->GetType(); + Float_t padlength= GetPadPitchLength(seed->GetRow()); + // + Double_t angle2 = seed->GetSnp()*seed->GetSnp(); + // if (angle2<0.6) angle2 = 0.6; + angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2)); + // + //cluster "quality" + Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2())); + Float_t res; + // + if (fSectors==fInnerSec){ + snoise2 = gnoise1; + res = ggg1[amp]*z+glandau1[amp]*angle2; + if (ctype==0) res *= gcor01[rsigmaz]; + if ((ctype>0)){ + res+=0.002; + res*= gcorp[rsigmaz]; + } + } + else { + if (padlength<1.1){ + snoise2 = gnoise2; + res = ggg2[amp]*z+glandau2[amp]*angle2; + if (ctype==0) res *= gcor02[rsigmaz]; + if ((ctype>0)){ + res+=0.002; + res*= gcorp[rsigmaz]; + } + } + else{ + snoise2 = gnoise3; + res = ggg3[amp]*z+glandau3[amp]*angle2; + if (ctype==0) res *= gcor02[rsigmaz]; + if ((ctype>0)){ + res+=0.002; + res*= gcorp[rsigmaz]; + } + } + } - if ((ctype==0) && (fSectors ==fInnerSec)) - res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2)); - if ((ctype>0)) - res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064; //0.31+0.147*ctype; - if (ctype<0) + if (ctype<0){ + res+=0.002; res*=1.3; - if ((ctype<0) &&<70) + } + if ((ctype<0) &&<70){ + res+=0.002; res*=1.3; - + } res += snoise2; if (res<2*snoise2) res = 2*snoise2; - + if (res>3) res =3; seed->SetErrorZ2(res); return res; } - - -void AliTPCseed::Reset() -{ +/* +Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ // - //PH SetN(0); - fNFoundable = 0; - ResetCovariance(); - SetChi2(0); - for (Int_t i=0;i<200;i++) fClusterIndex[i]=-1; -} - + // + //seed->SetErrorZ2(0.1); + //return 0.1; -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; - // 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); + Float_t snoise2; + Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ())); + // + Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2); + Int_t ctype = cl->GetType(); + Float_t amp = TMath::Abs(cl->GetQ()); - y = fP0; - z = fP1; - y += dx*(c1+c2)/(r1+r2); - z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3; - return 0; -} + Float_t nel; + Float_t nprim; + // + Float_t landau=2 ; //landau fluctuation part + Float_t gg=2; // gg fluctuation part + Float_t padlength= GetPadPitchLength(seed->GetX()); + + if (fSectors==fInnerSec){ + snoise2 = 0.0004/padlength; + nel = 0.268*amp; + nprim = 0.155*amp; + gg = (2+0.001*nel/(padlength*padlength))/nel; + landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (landau>1) landau=1; + } + else { + snoise2 = 0.0004/padlength; + nel = 0.3*amp; + nprim = 0.133*amp; + gg = (2+0.0008*nel/(padlength*padlength))/nel; + landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; + if (landau>1) landau=1; + } + Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z; + // + Float_t angle2 = seed->GetSnp()*seed->GetSnp(); + angle2 = TMath::Sqrt((1-angle2)); + if (angle2<0.6) angle2 = 0.6; + //angle2 = 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; + Float_t angle = seed->GetTgl()/angle2; + Float_t angular = landau*angle*angle*padlength*padlength/12.; + Float_t res = sdiff + angular; - 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; + if ((ctype==0) && (fSectors ==fOuterSec)) + res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2)); + + if ((ctype==0) && (fSectors ==fInnerSec)) + res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2)); - return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; -} + if ((ctype>0)){ + res+=0.005; + res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype; + } + if (ctype<0){ + res+=0.002; + res*=1.3; + } + if ((ctype<0) &&<70){ + res+=0.002; + res*=1.3; + } + res += snoise2; + if (res<2*snoise2) + res = 2*snoise2; + seed->SetErrorZ2(res); + return res; +} +*/ -//_________________________________________________________________________________________ -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 (t->fRelativeSector>fRelativeSector) return -1; - if (t->fRelativeSectorGetZ(); - Double_t z1 = GetZ(); - if (z2>z1) return 1; - if (z2GetX(); Float_t y = seed->GetY(); Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha()); + if (y > ymax) { - seed->fRelativeSector= (seed->fRelativeSector+1) % fN; + seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN); if (!seed->Rotate(fSectors->GetAlpha())) return; } else if (y <-ymax) { - seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN; + seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN); if (!seed->Rotate(-fSectors->GetAlpha())) return; } @@ -476,91 +982,85 @@ void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed) - //_____________________________________________________________________________ -Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) { +Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1, + Double_t x2,Double_t y2, + Double_t x3,Double_t y3) +{ //----------------------------------------------------------------- - // This function associates a cluster with this track. + // Initial approximation of the track curvature //----------------------------------------------------------------- - // Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2(); - //Double_t r00=sigmay2, r01=0., r11=sigmaz2; - Double_t r00=fErrorY2, r01=0., r11=fErrorZ2; - - r00+=fC00; r01+=fC10; r11+=fC11; - Double_t det=r00*r11 - r01*r01; - Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; - - Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11; - Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11; - Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11; - Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11; - Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11; - - Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1; - Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz; - if (TMath::Abs(cur*fX-eta) >= 0.9) { - // Int_t n=GetNumberOfClusters(); - //if (n>4) cerr<GetEntries()); - for (Int_t i=0; i1) return 0; + // Double_t angle2 = TMath::ASin(d*c*0.5); + // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5); + Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5); + + angle2 = (z1-z2)*c/(angle2*2.); + return angle2; +} + +Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) +{//----------------------------------------------------------------- + // This function find proloncation of a track to a reference plane x=x2. + //----------------------------------------------------------------- + + Double_t dx=x2-x1; + + if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) { + return kFALSE; } - LoadOuterSectors(); - LoadInnerSectors(); + Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1); + Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2); + y = x[0]; + z = x[1]; + + Double_t dy = dx*(c1+c2)/(r1+r2); + Double_t dz = 0; + // + Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1); + + if (TMath::Abs(delta)>0.01){ + dz = x[3]*TMath::ASin(delta)/x[4]; + }else{ + dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4]; + } + + //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4]; - return 0; + y+=dy; + z+=dz; + + return kTRUE; } -void AliTPCtrackerMI::UnloadClusters() +Int_t AliTPCtrackerMI::LoadClusters (TTree *tree) +{ + // + // + fInput = tree; + return LoadClusters(); +} + +Int_t AliTPCtrackerMI::LoadClusters() { // // load clusters to the memory - Int_t j=Int_t(fClustersArray.GetTree()->GetEntries()); + AliTPCClustersRow *clrow= new AliTPCClustersRow; + clrow->SetClass("AliTPCclusterMI"); + clrow->SetArray(0); + clrow->GetArray()->ExpandCreateFast(10000); + // + // TTree * tree = fClustersArray.GetTree(); + + TTree * tree = fInput; + TBranch * br = tree->GetBranch("Segment"); + br->SetAddress(&clrow); + // + Int_t j=Int_t(tree->GetEntries()); for (Int_t i=0; iGetEntry(i); + // + Int_t sec,row; + fParam->AdjustSectorRow(clrow->GetID(),sec,row); + for (Int_t icl=0; iclGetArray()->GetEntriesFast(); icl++){ + Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl))); + } + // + AliTPCRow * tpcrow=0; + Int_t left=0; + if (secSetN1(clrow->GetArray()->GetEntriesFast()); + tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]); + for (Int_t i=0;iGetN1();i++) + tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i))); + } + if (left ==1){ + tpcrow->SetN2(clrow->GetArray()->GetEntriesFast()); + tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]); + for (Int_t i=0;iGetN2();i++) + tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i))); + } } + // + delete clrow; + LoadOuterSectors(); + LoadInnerSectors(); + return 0; } +void AliTPCtrackerMI::UnloadClusters() +{ + // + // unload clusters from the memory + // + Int_t nrows = fOuterSec->GetNRows(); + for (Int_t sec = 0;secfClusters1) delete []tpcrow->fClusters1; + // if (tpcrow->fClusters2) delete []tpcrow->fClusters2; + //} + tpcrow->ResetClusters(); + } + // + nrows = fInnerSec->GetNRows(); + for (Int_t sec = 0;secfClusters1) delete []tpcrow->fClusters1; + //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; + //} + tpcrow->ResetClusters(); + } + + return ; +} + +void AliTPCtrackerMI::Transform(AliTPCclusterMI * 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]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; + if (mat) mat->LocalToMaster(pos,posC); + else{ + // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX + } + //cluster->SetX(posC[0]); + //cluster->SetY(posC[1]); + //cluster->SetZ(posC[2]); +} //_____________________________________________________________________________ -void AliTPCtrackerMI::LoadOuterSectors() { +Int_t AliTPCtrackerMI::LoadOuterSectors() { //----------------------------------------------------------------- // This function fills outer TPC sectors with clusters. //----------------------------------------------------------------- - UInt_t index; - //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries()); - Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal(); - for (Int_t i=0; i(fClustersArray.At(i)); - if (!s) continue; - Int_t sec,row; - AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam(); - par->AdjustSectorRow(s->GetID(),sec,row); - if (secGetArray()->GetEntriesFast(); - while (ncl--) { - AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl]; - index=(((sec<<8)+row)<<16)+ncl; - fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index); - } - } + Int_t nrows = fOuterSec->GetNRows(); + UInt_t index=0; + for (Int_t sec = 0;secGetN1(); + while (ncl--) { + AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl)); + index=(((sec2<<8)+row)<<16)+ncl; + tpcrow->InsertCluster(c,index); + } + //right + ncl = tpcrow->GetN2(); + while (ncl--) { + AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl)); + index=((((sec2+fkNOS)<<8)+row)<<16)+ncl; + tpcrow->InsertCluster(c,index); + } + // + // write indexes for fast acces + // + for (Int_t i=0;i<510;i++) + tpcrow->SetFastCluster(i,-1); + for (Int_t i=0;iGetN();i++){ + Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.); + tpcrow->SetFastCluster(zi,i); // write index + } + Int_t last = 0; + for (Int_t i=0;i<510;i++){ + if (tpcrow->GetFastCluster(i)<0) + tpcrow->SetFastCluster(i,last); + else + last = tpcrow->GetFastCluster(i); + } + } fN=fkNOS; fSectors=fOuterSec; + return 0; } //_____________________________________________________________________________ -void AliTPCtrackerMI::LoadInnerSectors() { +Int_t AliTPCtrackerMI::LoadInnerSectors() { //----------------------------------------------------------------- // This function fills inner TPC sectors with clusters. //----------------------------------------------------------------- - UInt_t index; - //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries()); - Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal(); - for (Int_t i=0; i(fClustersArray.At(i)); - if (!s) continue; - Int_t sec,row; - AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam(); - par->AdjustSectorRow(s->GetID(),sec,row); - if (sec>=fkNIS*2) continue; - AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row); - Int_t ncl=clrow->GetArray()->GetEntriesFast(); - while (ncl--) { - AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl]; - index=(((sec<<8)+row)<<16)+ncl; - fInnerSec[sec%fkNIS][row].InsertCluster(c,index); - } - } + Int_t nrows = fInnerSec->GetNRows(); + UInt_t index=0; + for (Int_t sec = 0;secGetN1(); + while (ncl--) { + AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl)); + index=(((sec<<8)+row)<<16)+ncl; + tpcrow->InsertCluster(c,index); + } + //right + ncl = tpcrow->GetN2(); + while (ncl--) { + AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl)); + index=((((sec+fkNIS)<<8)+row)<<16)+ncl; + tpcrow->InsertCluster(c,index); + } + // + // write indexes for fast acces + // + for (Int_t i=0;i<510;i++) + tpcrow->SetFastCluster(i,-1); + for (Int_t i=0;iGetN();i++){ + Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.); + tpcrow->SetFastCluster(zi,i); // write index + } + Int_t last = 0; + for (Int_t i=0;i<510;i++){ + if (tpcrow->GetFastCluster(i)<0) + tpcrow->SetFastCluster(i,last); + else + last = tpcrow->GetFastCluster(i); + } + } + fN=fkNIS; fSectors=fInnerSec; + return 0; } -Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { - //----------------------------------------------------------------- - // This function tries to find a track prolongation to next pad row - //----------------------------------------------------------------- - // Double_t xt=t.GetX(); - // Int_t row = fSectors->GetRowNumber(xt)-1; - // if (row < nr) return 1; // don't prolongate if not information until now - - // - Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr); - // if (t.GetRadius()>x+10 ) return 0; - if (!t.PropagateTo(x)) { - t.fStopped = kTRUE; - return 0; + +//_________________________________________________________________________ +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; } - // update current - t.fCurrentSigmaY = GetSigmaY(&t); - t.fCurrentSigmaZ = GetSigmaZ(&t); - // - AliTPCclusterMI *cl=0; - UInt_t index=0; - const AliTPCRow &krow=fSectors[t.fRelativeSector][nr]; - Double_t sy2=ErrY2(&t)*2; - Double_t sz2=ErrZ2(&t)*2; + if (secGetN1()<=ncl) return 0; + clrow = tpcrow->GetClusters1(); + } + else { + if (tpcrow->GetN2()<=ncl) return 0; + clrow = tpcrow->GetClusters2(); + } + } + else { + tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]); + if (tpcrow==0) return 0; - if (TMath::Abs(TMath::Abs(y)-ymax)GetNRowLow(); - t.fClusterIndex[row] = -1; - return 0; - } - else + if (sec-2*fkNISGetN1()<=ncl) return 0; + clrow = tpcrow->GetClusters1(); + } + else { + if (tpcrow->GetN2()<=ncl) return 0; + clrow = tpcrow->GetClusters2(); + } + } + + return &(clrow[ncl]); + +} + + + +Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { + //----------------------------------------------------------------- + // This function tries to find a track prolongation to next pad row + //----------------------------------------------------------------- + // + Double_t x= GetXrow(nr), ymax=GetMaxY(nr); + 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.GetClusterPointer(nr); + if ( (cl==0) ) cl = GetClusterMI(tpcindex); + t.SetCurrentClusterIndex1(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.SetRelativeSector(relativesector); + if (!t.Rotate(rotation)) return 0; + } + if (!t.PropagateTo(x)) return 0; + // + t.SetCurrentCluster(cl); + t.SetRow(nr); + Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),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.SetErrorY2(t.GetErrorY2()+0.03); + t.SetErrorZ2(t.GetErrorZ2()+0.03); + t.SetErrorY2(t.GetErrorY2()*3); + t.SetErrorZ2(t.GetErrorZ2()*3); + } + t.SetNFoundable(t.GetNFoundable()+1); + 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.SetNFoundable(t.GetNFoundable()+1); + return 0; + } + // + 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.SetRelativeSector((t.GetRelativeSector()+1) % fN); + if (!t.Rotate(fSectors->GetAlpha())) + return 0; + } else if (y <-ymax) { + t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN); + if (!t.Rotate(-fSectors->GetAlpha())) + return 0; + } + //return 1; + } + // + if (!t.PropagateTo(x)) { + if (fIteration==0) t.SetRemoval(10); + return 0; + } + y=t.GetY(); + Double_t z=t.GetZ(); + // + + if (!IsActive(t.GetRelativeSector(),nr)) { + t.SetInDead(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.GetRelativeSector(),nr); + Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0; + + if (!isActive || !isActive2) return 0; + + const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr); + if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0; + Double_t roady =1.; + Double_t roadz = 1.; + // + if (TMath::Abs(TMath::Abs(y)-ymax)GetZLength() && (TMath::Abs(t.GetSnp())GetZ() > z+roadz) break; - if ( (c->GetY()-y) > roady ) continue; - Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y); - if (maxdistance>distance) { - maxdistance = distance; - cl=c; - // index=krow.GetIndex(i); - index =i; - } - } - } + // cl = krow.FindNearest2(y+10.,z,roady,roadz,index); + cl = krow.FindNearest2(y,z,roady,roadz,index); + if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index)); + } if (cl) { - // Double_t sy2= ErrY2(&t,cl); - // Double_t sz2= ErrZ2(&t,cl); - // Double_t chi2= t.GetPredictedChi2(cl); - // UpdateTrack(&t,cl,chi2,index); - - t.fCurrentCluster = cl; - t.fCurrentClusterIndex1 = krow.GetIndex(index); - t.fCurrentClusterIndex2 = index; - Double_t sy2=ErrY2(&t,t.fCurrentCluster); - Double_t sz2=ErrZ2(&t,t.fCurrentCluster); + t.SetCurrentCluster(cl); + t.SetRow(nr); + if (fIteration==2&&cl->IsUsed(10)) return 0; + Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.); + if (fIteration==2&&cl->IsUsed(11)) { + t.SetErrorY2(t.GetErrorY2()+0.03); + t.SetErrorZ2(t.GetErrorZ2()+0.03); + t.SetErrorY2(t.GetErrorY2()*3); + t.SetErrorZ2(t.GetErrorZ2()*3); + } + /* + if (t.fCurrentCluster->IsUsed(10)){ + // + // - Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2()); - Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2()); + t.fNShared++; + if (t.fNShared>0.7*t.GetNumberOfClusters()) { + t.fRemoval =10; + return 0; + } + } + */ + if (accept<3) UpdateTrack(&t,accept); - Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY()); - Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ()); + } else { + if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10); - Double_t rdistance = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2)); - - - // printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance); - if ( (rdistancey>1) || (rdistancez>1)) return 0; - if (rdistance>4) return 0; - - if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0) - return 0; //suspisiouce - will be changed - - if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0) - // strict cut on overlaped cluster - return 0; //suspisiouce - will be changed - - if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 ) - && t.fCurrentCluster->GetType()<0) - return 0; - - // t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t)); - UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1); + } + return 1; +} - } else { +Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) { + //----------------------------------------------------------------- + // This function tries to find a track prolongation to next pad row + //----------------------------------------------------------------- + // + Double_t x= GetXrow(nr), ymax=GetMaxY(nr); + Double_t y,z; + if (!t.GetProlongation(x,y,z)) { + t.SetRemoval(10); + return 0; + } + // + // + if (TMath::Abs(y)>ymax){ + if (y > ymax) { - t.fRelativeSector= (t.fRelativeSector+1) % fN; + t.SetRelativeSector((t.GetRelativeSector()+1) % fN); if (!t.Rotate(fSectors->GetAlpha())) return 0; } else if (y <-ymax) { - t.fRelativeSector= (t.fRelativeSector-1+fN) % fN; + t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN); if (!t.Rotate(-fSectors->GetAlpha())) return 0; + } + if (!t.PropagateTo(x)) { + return 0; + } + t.GetProlongation(x,y,z); + } + // + // update current shape info every 3 pad-row + if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){ + // t.fCurrentSigmaY = GetSigmaY(&t); + //t.fCurrentSigmaZ = GetSigmaZ(&t); + GetShape(&t,nr); + } + // + AliTPCclusterMI *cl=0; + UInt_t index=0; + + + //Int_t nr2 = nr; + const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr); + if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0; + Double_t roady =1.; + Double_t roadz = 1.; + // + Int_t row = nr; + if (TMath::Abs(TMath::Abs(y)-ymax)(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1); } + //calculate + + if ((cl==0)&&(krow)) { + // cl = krow.FindNearest2(y+10,z,roady,roadz,index); + cl = krow.FindNearest2(y,z,roady,roadz,index); + + if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index)); + } + + if (cl) { + t.SetCurrentCluster(cl); + // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.); + //if (accept<3){ + t.SetClusterIndex2(row,index); + t.SetClusterPointer(row, cl); + //} } return 1; } -Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex, 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 //----------------------------------------------------------------- - t.fCurrentCluster = 0; - t.fCurrentClusterIndex1 = 0; - t.fCurrentClusterIndex2 = 0; + t.SetCurrentCluster(0); + t.SetCurrentClusterIndex1(0); Double_t xt=t.GetX(); - Int_t row = fSectors->GetRowNumber(xt)-1; + Int_t row = GetRowNumber(xt)-1; + Double_t ymax= GetMaxY(nr); + if (row < nr) return 1; // don't prolongate if not information until now - - Double_t x=fSectors->GetX(nr); - // if (t.fStopped) return 0; - // if (t.GetRadius()>x+10 ) return 0; +// 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; + //t.PropagateTo(x+0.02); + //t.PropagateTo(x+0.01); if (!t.PropagateTo(x)){ - t.fStopped =kTRUE; return 0; } - // update current - t.fCurrentSigmaY = GetSigmaY(&t); - t.fCurrentSigmaZ = GetSigmaZ(&t); - - AliTPCclusterMI *cl=0; - UInt_t index=0; - AliTPCRow &krow=fSectors[t.fRelativeSector][nr]; - // - Double_t y=t.GetY(), z=t.GetZ(); - Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY); - Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ); // + y=t.GetY(); + z=t.GetZ(); - Float_t maxdistance = 1000000; - if (krow) { - for (Int_t i=krow.Find(z-roadz); iGetZ() > z+roadz) break; - if (TMath::Abs(c->GetY()-y)>roady) continue; - - //krow.UpdateClusterTrack(i,trindex,&t); - - Float_t dy2 = (c->GetY()- t.GetY()); - dy2*=dy2; - Float_t dz2 = (c->GetZ()- t.GetZ()); - dz2*=dz2; - // - Float_t distance = dy2+dz2; - // - if (distance > maxdistance) continue; - maxdistance = distance; - cl=c; - index=i; + if (TMath::Abs(y)>ymax){ + if (y > ymax) { + t.SetRelativeSector((t.GetRelativeSector()+1) % fN); + if (!t.Rotate(fSectors->GetAlpha())) + return 0; + } else if (y <-ymax) { + t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN); + if (!t.Rotate(-fSectors->GetAlpha())) + return 0; } + // if (!t.PropagateTo(x)){ + // return 0; + //} + return 1; + //y = t.GetY(); } - t.fCurrentCluster = cl; - t.fCurrentClusterIndex1 = krow.GetIndex(index); - t.fCurrentClusterIndex2 = index; - return 1; -} + // + if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; + if (!IsActive(t.GetRelativeSector(),nr)) { + t.SetInDead(kTRUE); + t.SetClusterIndex2(nr,-1); + return 0; + } + //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha())); -Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) { - //----------------------------------------------------------------- - // This function tries to find a track prolongation to next pad row - //----------------------------------------------------------------- - AliTPCseed & t = *((AliTPCseed*)(fSeeds->At(trindex))); - AliTPCRow &krow=fSectors[t.fRelativeSector][nr]; - // Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt(); - Double_t y=t.GetY(); - Double_t ymax=fSectors->GetMaxY(nr); - - if (TMath::Abs(TMath::Abs(y)-ymax)GetNRowLow(); - t.fClusterIndex[row] = -1; + AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr); + + if (TMath::Abs(TMath::Abs(y)-ymax)GetPadPitchLength(); - // AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1); - Double_t sy2=ErrY2(&t,t.fCurrentCluster); - Double_t sz2=ErrZ2(&t,t.fCurrentCluster); + // update current + if ( (nr%6==0) || t.GetNumberOfClusters()<2){ + // t.fCurrentSigmaY = GetSigmaY(&t); + //t.fCurrentSigmaZ = GetSigmaZ(&t); + GetShape(&t,nr); + } + + AliTPCclusterMI *cl=0; + Int_t index=0; + // + Double_t roady = 1.; + Double_t roadz = 1.; + // + if (!cl){ + index = t.GetClusterIndex2(nr); + if ( (index>0) && (index&0x8000)==0){ + cl = t.GetClusterPointer(nr); + if ( (cl==0) && (index>0)) cl = GetClusterMI(index); + t.SetCurrentClusterIndex1(index); + if (cl) { + t.SetCurrentCluster(cl); + return 1; + } + } + } - Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2()); - Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2()); + // if (index<0) return 0; + UInt_t uindex = TMath::Abs(index); - Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY()); - Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ()); - - Double_t rdistance = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2)); + if (krow) { + //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex); + cl = krow.FindNearest2(y,z,roady,roadz,uindex); + } + if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex)); + t.SetCurrentCluster(cl); - // printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance); - if ( (rdistancey>1) || (rdistancez>1)) return 0; - if (rdistance>4) return 0; + return 1; +} - if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0) - return 0; //suspisiouce - will be changed - if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0) - // strict cut on overlaped cluster - return 0; //suspisiouce - will be changed +Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) { + //----------------------------------------------------------------- + // This function tries to find a track prolongation to next pad row + //----------------------------------------------------------------- - if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 ) - && t.fCurrentCluster->GetType()<0) - return 0; + //update error according neighborhoud - // t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t)); - UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1); - - } else { - if (y > ymax) { - t.fRelativeSector= (t.fRelativeSector+1) % fN; - if (!t.Rotate(fSectors->GetAlpha())) - return 0; - } else if (y <-ymax) { - t.fRelativeSector= (t.fRelativeSector-1+fN) % fN; - if (!t.Rotate(-fSectors->GetAlpha())) + if (t.GetCurrentCluster()) { + t.SetRow(nr); + Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.); + + if (t.GetCurrentCluster()->IsUsed(10)){ + // + // + // t.fErrorZ2*=2; + // t.fErrorY2*=2; + t.SetNShared(t.GetNShared()+1); + if (t.GetNShared()>0.7*t.GetNumberOfClusters()) { + t.SetRemoval(10); return 0; + } + } + if (fIteration>0) accept = 0; + if (accept<3) UpdateTrack(&t,accept); + + } else { + if (fIteration==0){ + if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10); + if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10); + + if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10); } } return 1; @@ -923,129 +1838,60 @@ Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) { - -/* -Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step) -{ +//_____________________________________________________________________________ +Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) { //----------------------------------------------------------------- - // fast prolongation mathod - - // don't update track only after step clusters + // This function tries to find a track prolongation. //----------------------------------------------------------------- Double_t xt=t.GetX(); // - Double_t alpha=t.GetAlpha(); - alpha =- fSectors->GetAlphaShift(); + 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())%fN; - Int_t row0 = fSectors->GetRowNumber(xt); - Double_t x = fSectors->GetX(row0); - Double_t ymax = fSectors->GetMaxY(row0); - // - Double_t sy2=ErrY2(&t)*2; - Double_t sz2=ErrZ2(&t)*2; - Double_t roady =3.*sqrt(t.GetSigmaY2() + sy2); - Double_t roadz = 3 *sqrt(t.GetSigmaZ2() + sz2); - Float_t maxdistance = roady*roady + roadz*roadz; - t.fCurrentSigmaY = GetSigmaY(&t); - t.fCurrentSigmaZ = GetSigmaZ(&t); // - Int_t nclusters = 0; - Double_t y; - Double_t z; - Double_t yy[200]; //track prolongation - Double_t zz[200]; - Double_t cy[200]; // founded cluster position - Double_t cz[200]; - Double_t sy[200]; // founded cluster error - Double_t sz[200]; - Bool_t hitted[200]; // indication of cluster presence - // - - // - for (Int_t drow = step; drow>=0; drow--) { - Int_t row = row0-drow; - if (row<0) break; - Double_t x = fSectors->GetX(row); - Double_t ymax = fSectors->GetMaxY(row); - t.GetProlongation(x,y,z); - yy[drow] =y; - zz[drow] =z; - const AliTPCRow &krow=fSectors[t.fRelativeSector][row]; - if (TMath::Abs(TMath::Abs(y)-ymax)GetAlpha()+0.0001)%fN); - //find nearest cluster - AliTPCclusterMI *cl= 0; - if (krow) { - for (Int_t i=krow.Find(z-roadz); iGetZ() > z+roadz) break; - if ( (c->GetY()-y) > roady ) continue; - Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y); - if (maxdistance>distance) { - maxdistance = distance; - cl=c; - // index=krow.GetIndex(i); - } - } - } // end of seearch - //update cluster information - if (cl){ - cy[drow] = cl->GetY(); - cz[drow] = cl->GetZ(); - sy[drow] = ErrY2(&t,cl); - sz[drow] = ErrZ2(&t,cl); - hitted[drow] = kTRUE; - nclusters++; + Int_t first = GetRowNumber(xt)-1; + 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; + // + AliKink * kink = (AliKink*)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 - hitted[drow] = kFALSE; - } - //if we have information - update track - if (nclusters>0){ - Float_t sumyw0 = 0; - Float_t sumzw0 = 0; - Float_t sumyw = 0; - Float_t sumzw = 0; - Float_t sumrow = 0; - for (Int_t i=0;iGetX(mrow); - t.PropagateTo(x); - AliTPCclusterMI cvirtual; - cvirtual.SetZ(dz+t.GetZ()); - cvirtual.SetY(dy+t.GetY()); - t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters))); - t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters))); - Float_t chi2 = t.GetPredictedChi2(&cvirtual); - t.Update(&cvirtual,chi2,0); - Int_t ncl = t.GetNumberOfClusters(); - ncl = ncl-1+nclusters; - t.SetN(ncl); - } - return 1; -} -*/ + fSectors = fOuterSec; + if (FollowToNext(t,nr)==0) + if (!t.IsActive()) + return 0; + + } + return 1; +} + //_____________________________________________________________________________ -Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) { +Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) { //----------------------------------------------------------------- // This function tries to find a track prolongation. //----------------------------------------------------------------- @@ -1054,32 +1900,66 @@ Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) { 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())%fN; + t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN); + + for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) { + + if (FollowToNextFast(t,nr)==0) + if (!t.IsActive()) return 0; - for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) { - - if (FollowToNext(t,nr)==0) { - } } return 1; } -//_____________________________________________________________________________ + + + 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())%fN; + t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN); - for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) { - if (t.GetSnp()<0.8) - FollowToNext(t,nr); + Int_t first = t.GetFirstPoint(); + 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); + AliKink * kink = (AliKink*)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 + fSectors = fOuterSec; + + FollowToNext(t,nr); } return 1; } @@ -1098,43 +1978,34 @@ Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t & // Float_t dz2 =(s1->GetZ() - s2->GetZ()); dz2*=dz2; - /* - Float_t x = s1->GetX(); - Float_t x2 = s2->GetX(); - - Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha()); - */ + Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY())); - //if (TMath::Abs(dy2)>2*ymax-3) - // dy2-=2*ymax; dy2*=dy2; Float_t distance = TMath::Sqrt(dz2+dy2); if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics - Int_t offset =0; - if (fSectors==fOuterSec) offset = fParam->GetNRowLow(); - Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint); - Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint); - lastpoint +=offset; - firstpoint+=offset; + // Int_t offset =0; + Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint()); + Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint()); if (lastpoint>160) lastpoint =160; if (firstpoint<0) firstpoint = 0; - if (firstpointlastpoint) { + firstpoint =lastpoint; + // lastpoint =160; } - for (Int_t i=firstpoint;ifClusterIndex[i]>0) sum1++; - if (s2->fClusterIndex[i]>0) sum2++; - if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) { + for (Int_t i=firstpoint-1;iGetClusterIndex2(i)>0) sum1++; + if (s2->GetClusterIndex2(i)>0) sum2++; + if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) { sum++; } } - + if (sum<5) return 0; + Float_t summin = TMath::Min(sum1+1,sum2+1); Float_t ratio = (sum+1)/Float_t(summin); return ratio; @@ -1144,132 +2015,210 @@ void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2) { // // - if (s1->fSector!=s2->fSector) return; - // + if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return; + if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return; + Float_t dz2 =(s1->GetZ() - s2->GetZ()); dz2*=dz2; Float_t dy2 =(s1->GetY() - s2->GetY()); - dy2*=dy2; - Float_t distance = TMath::Sqrt(dz2+dy2); - if (distance>15.) return ; // if there are far away - not overlap - to reduce combinatorics - //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint; - // TClonesArray &pointarray1 = *(s1->fPoints); - //TClonesArray &pointarray2 = *(s2->fPoints); - // - for (Int_t i=0;i<160;i++){ - if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) { - // AliTPCTrackPoint *p1 = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i)); - //AliTPCTrackPoint *p2 = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i)); - AliTPCTrackPoint *p1 = s1->GetTrackPoint(i); - AliTPCTrackPoint *p2 = s2->GetTrackPoint(i);; - p1->fIsShared = kTRUE; - p2->fIsShared = kTRUE; - } - } -} - - - - -void AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared){ - + Float_t distance = dz2+dy2; + if (distance>325.) return ; // if there are far away - not overlap - to reduce combinatorics - // - // remove overlap - used removal factor - removal index stored in the track - for (Int_t i=0; iGetEntriesFast(); i++) { - AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i); - if (pt) RotateToLocal(pt); - } - arr->Sort(); // sorting according z - arr->Expand(arr->GetEntries()); - Int_t nseed=arr->GetEntriesFast(); - // printf("seeds \t%p \t%d\n",arr, nseed); - // arr->Expand(arr->GetEntries()); //remove 0 pointers - nseed = arr->GetEntriesFast(); - Int_t removed = 0; + Int_t sumshared=0; + // + Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint()); + Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint()); + // + if (firstpoint>=lastpoint-5) return;; + + for (Int_t i=firstpoint;iGetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) { + if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) { + sumshared++; + } + } + if (sumshared>4){ + // sign clusters + // + for (Int_t i=firstpoint;iGetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) { + if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) { + AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i); + AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);; + if (s1->IsActive()&&s2->IsActive()){ + p1->SetShared(kTRUE); + p2->SetShared(kTRUE); + } + } + } + } + // + if (sumshared>10){ + for (Int_t i=0;i<4;i++){ + if (s1->GetOverlapLabel(3*i)==0){ + s1->SetOverlapLabel(3*i, s2->GetLabel()); + s1->SetOverlapLabel(3*i+1,sumshared); + s1->SetOverlapLabel(3*i+2,s2->GetUniqueID()); + break; + } + } + for (Int_t i=0;i<4;i++){ + if (s2->GetOverlapLabel(3*i)==0){ + s2->SetOverlapLabel(3*i, s1->GetLabel()); + s2->SetOverlapLabel(3*i+1,sumshared); + s2->SetOverlapLabel(3*i+2,s1->GetUniqueID()); + break; + } + } + } + +} + +void AliTPCtrackerMI::SignShared(TObjArray * arr) +{ + // + //sort trackss according sectors + // + for (Int_t i=0; iGetEntriesFast(); i++) { + AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i); + if (!pt) continue; + //if (pt) RotateToLocal(pt); + pt->SetSort(0); + } + arr->UnSort(); + arr->Sort(); // sorting according z + arr->Expand(arr->GetEntries()); + // + // + Int_t nseed=arr->GetEntriesFast(); for (Int_t i=0; iUncheckedAt(i); - if (!pt) { - continue; + if (!pt) continue; + for (Int_t j=0;j<=12;j++){ + pt->SetOverlapLabel(j,0); } - if (!(pt->IsActive())) continue; + } + for (Int_t i=0; iUncheckedAt(i); + if (!pt) continue; + if (pt->GetRemoval()>10) continue; for (Int_t j=i+1; jUncheckedAt(j); - // - if (!pt2) continue; - if (!(pt2->IsActive())) continue; - if (TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break; - if (TMath::Abs(pt2->GetZ()-pt->GetZ())<4){ - Int_t sum1,sum2; - Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2); - //if (sum1==0) { - // pt->Desactivate(removalindex); // arr->RemoveAt(i); - // break; - //} - if (ratio>factor){ - // if (pt->GetChi2()GetChi2()) pt2->Desactivate(removalindex); // arr->RemoveAt(j); - Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1); - Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2); - removed++; - if (TMath::Abs(ratio3)>0.025){ // if much more points - if (sum1>sum2) pt2->Desactivate(removalindex); - else { - pt->Desactivate(removalindex); // arr->RemoveAt(i); - break; - } - } - else{ //decide on mean chi2 - if (ratio2<1) - pt2->Desactivate(removalindex); - else { - pt->Desactivate(removalindex); // arr->RemoveAt(i); - break; - } - } - - } // if suspicious ratio + // if (pt2){ + if (pt2->GetRemoval()<=10) { + if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break; + SignShared(pt,pt2); } - else - break; + } + } +} + +void AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex) +{ + // + //sort trackss according sectors + // + if (fDebug&1) { + Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries()); + } + // + for (Int_t i=0; iGetEntriesFast(); i++) { + AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i); + if (!pt) continue; + pt->SetSort(0); + } + arr->UnSort(); + arr->Sort(); // sorting according z + arr->Expand(arr->GetEntries()); + // + //reset overlap labels + // + Int_t nseed=arr->GetEntriesFast(); + for (Int_t i=0; iUncheckedAt(i); + if (!pt) continue; + pt->SetUniqueID(i); + for (Int_t j=0;j<=12;j++){ + pt->SetOverlapLabel(j,0); } } - // printf("removed\t%d\n",removed); - Int_t good =0; + // + //sign shared tracks for (Int_t i=0; iUncheckedAt(i); - if (!pt) break; - if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) { - //desactivate tracks with small number of points - // printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable)); - pt->Desactivate(10); //desactivate - small muber of points + if (!pt) continue; + if (pt->GetRemoval()>10) continue; + Float_t deltac = pt->GetC()*0.1; + for (Int_t j=i+1; jUncheckedAt(j); + // if (pt2){ + if (pt2->GetRemoval()<=10) { + if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break; + if (TMath::Abs(pt->GetC() -pt2->GetC())>deltac) continue; + if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue; + // + SignShared(pt,pt2); + } } - if (!(pt->IsActive())) continue; - good++; } - - - if (shared) - for (Int_t i=0; iUncheckedAt(i); - if (!pt) continue; - if (!(pt->IsActive())) continue; - for (Int_t j=i+1; jUncheckedAt(j); - if ((pt2) && pt2->IsActive()) { - if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break; - SignShared(pt,pt2); + // + // remove highly shared tracks + for (Int_t i=0; iUncheckedAt(i); + if (!pt) continue; + if (pt->GetRemoval()>10) continue; + // + Int_t sumshared =0; + for (Int_t j=0;j<4;j++){ + sumshared = pt->GetOverlapLabel(j*3+1); + } + Float_t factor = factor1; + if (pt->GetRemoval()>0) factor = factor2; + if (sumshared/pt->GetNumberOfClusters()>factor){ + for (Int_t j=0;j<4;j++){ + if (pt->GetOverlapLabel(3*j)==0) continue; + if (pt->GetOverlapLabel(3*j+1)<5) continue; + if (pt->GetRemoval()==removalindex) continue; + AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->GetOverlapLabel(3*j+2)); + if (!pt2) continue; + if (pt2->GetSigma2C()GetSigma2C()){ + // pt->fRemoval = removalindex; + delete arr->RemoveAt(i); + break; } - } + } } - fNtracks = good; - printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks); + } + arr->Compress(); + if (fDebug&1) { + Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries()); + } +} + + + + +void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const +{ + // + //sort tracks in array according mode criteria + Int_t nseed = arr->GetEntriesFast(); + for (Int_t i=0; iUncheckedAt(i); + if (!pt) { + continue; + } + pt->SetSort(mode); + } + arr->UnSort(); + arr->Sort(); } -void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor, Int_t removalindex) +void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex) { //Loop over all tracks and remove "overlaps" @@ -1277,426 +2226,3232 @@ void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor, Int_t removali // Int_t nseed = arr->GetEntriesFast(); Int_t good =0; + for (Int_t i=0; iUncheckedAt(i); if (!pt) { - continue; + delete arr->RemoveAt(i); } - if (!(pt->IsActive())) continue; - Int_t noc=pt->GetNumberOfClusters(); - Int_t shared =0; - for (Int_t i=0; iGetClusterIndex(i); - AliTPCclusterMI *c=(AliTPCclusterMI*)GetClusterMI(index); - if (!c) continue; - if (c->IsUsed()) shared++; + else{ + pt->SetSort(1); + pt->SetBSigned(kFALSE); } - if ((Float_t(shared)/Float_t(noc))>factor) + } + arr->Compress(); + nseed = arr->GetEntriesFast(); + arr->UnSort(); + arr->Sort(); + // + //unsign used + UnsignClusters(); + // + for (Int_t i=0; iUncheckedAt(i); + if (!pt) { + continue; + } + Int_t found,foundable,shared; + if (pt->IsActive()) + pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE); + else + pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE); + // + Double_t factor = factor2; + if (pt->GetBConstrain()) factor = factor1; + + if ((Float_t(shared)/Float_t(found))>factor){ pt->Desactivate(removalindex); - else{ - good++; - for (Int_t i=0; iGetClusterIndex(i); - AliTPCclusterMI *c=(AliTPCclusterMI*)GetClusterMI(index); - if (!c) continue; - c->Use(); - } + continue; + } + + good++; + for (Int_t i=0; i<160; i++) { + Int_t index=pt->GetClusterIndex2(i); + if (index<0 || index&0x8000 ) continue; + AliTPCclusterMI *c= pt->GetClusterPointer(i); + if (!c) continue; + // if (!c->IsUsed(10)) c->Use(10); + //if (pt->IsActive()) + c->Use(10); + //else + // c->Use(5); } + } fNtracks = good; - printf("\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks); - + if (fDebug>0){ + Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks); + } } -void AliTPCtrackerMI::MakeSeedsAll() +void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal) { - if (fSeeds == 0) fSeeds = new TObjArray; - TObjArray * arr; - for (Int_t sec=0;secGetEntriesFast(); - for (Int_t i=0;iAddLast(arr->RemoveAt(i)); - } - // fSeeds = MakeSeedsSectors(0,fkNOS); -} -TObjArray * AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2) -{ + //Loop over all tracks and remove "overlaps" // - // loop over all sectors and make seed - //find track seeds - TStopwatch timer; - Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows(); - Int_t nrows=nlow+nup; - Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap); - // if (fSeeds==0) fSeeds = new TObjArray; - TObjArray * arr = new TObjArray; + // + 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->GetBConstrain()) ? 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->GetESD()){ + Int_t dummy[12]; + if (pt->GetESD()->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; + } + } - for (Int_t sec=sec1; sec0.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->GetClusterPointer(i); + if (!c) continue; + c->Use(10); + } } - gap = Int_t(0.2* nrows); - for (Int_t sec=sec1; sec0){ + 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 + // - Int_t nseed=arr->GetEntriesFast(); - Int_t i; - - gap=Int_t(0.3*nrows); - // continue seeds - for (i=0; iUncheckedAt(i), &t=*pt; - if (!pt) continue; - if (FollowProlongation(t,nup-gap)) { - pt->fIsSeeding =kFALSE; - continue; - } - delete arr->RemoveAt(i); + for (Int_t sec=0;secGetNRows();row++){ + AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1(); + for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++) + // if (cl[icl].IsUsed(10)) + cl[icl].Use(-1); + cl = fInnerSec[sec][row].GetClusters2(); + for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++) + //if (cl[icl].IsUsed(10)) + cl[icl].Use(-1); } + } - - // - //remove seeds which overlaps - RemoveOverlap(arr,0.4,1); - //delete seeds - which were sign - nseed=arr->GetEntriesFast(); - for (i=0; iUncheckedAt(i); - //, &t=*pt; - if (!pt) continue; - if ((pt->IsActive()) && pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) { - //pt->Reset(); - //FollowBackProlongation(*pt,nup-1); - //if ( pt->GetNumberOfClusters() < pt->fNFoundable*0.5 || pt->GetNumberOfClusters()<10 ) - //delete arr->RemoveAt(i); - //else - // pt->Reset(); - continue; + for (Int_t sec=0;secGetNRows();row++){ + AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1(); + for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++) + //if (cl[icl].IsUsed(10)) + cl[icl].Use(-1); + cl = fOuterSec[sec][row].GetClusters2(); + for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++) + //if (cl[icl].IsUsed(10)) + cl[icl].Use(-1); } - delete arr->RemoveAt(i); - } - //RemoveOverlap(arr,0.6,1); - return arr; + } + } -//_____________________________________________________________________________ -void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) { - //----------------------------------------------------------------- - // This function creates track seeds. - //----------------------------------------------------------------- - // if (fSeeds==0) fSeeds=new TObjArray(15000); - - Double_t x[5], c[15]; - - Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift(); - Double_t cs=cos(alpha), sn=sin(alpha); - - Double_t x1 =fOuterSec->GetX(i1); - Double_t xx2=fOuterSec->GetX(i2); +void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity) +{ // - // for (Int_t ns=0; nsGetY(), z1=kr1[is]->GetZ(); - Double_t x3=GetX(), y3=GetY(), z3=GetZ(); - - Float_t anglez = (z1-z3)/(x1-x3); - Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z - - for (Int_t js=0; js < nl+nm+nu; js++) { - const AliTPCclusterMI *kcl; - Double_t x2, y2, z2; - if (js=nl) continue; - } - kcl=kr21[js]; - z2=kcl->GetZ(); - if ((extraz-z2)>10) continue; - if ((extraz-z2)<-10) { - js = nl-1; - continue; - } - y2=kcl->GetY(); - x2= xx2*cs+y2*sn; - y2=-xx2*sn+y2*cs; - } else - if (js=nl+nm) continue; - } - kcl=kr22[js-nl]; - z2=kcl->GetZ(); - if ((extraz-z2)>10) continue; - if ((extraz-z2)<-10) { - js = nl+nm-1; - continue; - } - x2=xx2; y2=kcl->GetY(); - } else { - //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2]; - if (js==nl+nm) { - js = nl+nm+kr23.Find(extraz-15.); - if (js>=nl+nm+nu) break; - } - kcl=kr23[js-nl-nm]; - z2=kcl->GetZ(); - if ((extraz-z2)>10) continue; - if ((extraz-z2)<-10) { - break; - } - y2=kcl->GetY(); - x2=xx2*cs-y2*sn; - y2=xx2*sn+y2*cs; - } - - Double_t zz=z1 - anglez*(x1-x2); - if (TMath::Abs(zz-z2)>10.) continue; - - Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1); - if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;} - - x[0]=y1; - x[1]=z1; - x[4]=f1(x1,y1,x2,y2,x3,y3); - if (TMath::Abs(x[4]) >= 0.0066) continue; - x[2]=f2(x1,y1,x2,y2,x3,y3); - //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue; - x[3]=f3(x1,y1,x2,y2,z1,z2); - if (TMath::Abs(x[3]) > 1.2) continue; - Double_t a=asin(x[2]); - Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2])); - if (TMath::Abs(zv-z3)>10.) continue; - - Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4; - Double_t sy2=kcl->GetSigmaY2()*2, sz2=kcl->GetSigmaZ2()*4; - //Double_t sy3=400*3./12., sy=0.1, sz=0.1; - Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1; - //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1; + //sign clusters to be "used" + // + // snumber and sdensity sign number of sigmas - bellow mean value to be accepted + // loop over "primaries" + + Float_t sumdens=0; + Float_t sumdens2=0; + Float_t sumn =0; + Float_t sumn2 =0; + Float_t sumchi =0; + Float_t sumchi2 =0; - Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy; - Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy; - Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy; - Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; - Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; - Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; - Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy; - Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz; - Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy; - Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz; + Float_t sum =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; + TStopwatch timer; + timer.Start(); - UInt_t index=kr1.GetIndex(is); - AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift); - track->fIsSeeding = kTRUE; - //Int_t rc=FollowProlongation(*track, i2); - Int_t delta4 = Int_t((i2-i1)/4.); - - FollowProlongation(*track, i1-delta4); - if (track->GetNumberOfClusters() < track->fNFoundable/2.) { - delete track; - continue; - } - FollowProlongation(*track, i1-2*delta4); - if (track->GetNumberOfClusters() < track->fNFoundable/2.) { - delete track; - continue; - } - FollowProlongation(*track, i1-3*delta4); - if (track->GetNumberOfClusters() < track->fNFoundable/2.) { - delete track; - continue; - } - FollowProlongation(*track, i2); - //Int_t rc = 1; - - track->fLastPoint = i1; // first cluster in track position - if (track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track; - else arr->AddLast(track); - } + Int_t nseed = arr->GetEntriesFast(); + for (Int_t i=0; iUncheckedAt(i); + if (!pt) { + continue; + } + if (!(pt->IsActive())) continue; + Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable()); + if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){ + sumdens += dens; + sumdens2+= dens*dens; + sumn += pt->GetNumberOfClusters(); + sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters(); + Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters(); + if (chi2>5) chi2=5; + sumchi +=chi2; + sumchi2 +=chi2*chi2; + sum++; } } -} + Float_t mdensity = 0.9; + Float_t meann = 130; + Float_t meanchi = 1; + Float_t sdensity = 0.1; + Float_t smeann = 10; + Float_t smeanchi =0.4; + -//_____________________________________________________________________________ -void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) { - //----------------------------------------------------------------- - // This function creates track seeds - without vertex constraint - //----------------------------------------------------------------- + if (sum>20){ + mdensity = sumdens/sum; + meann = sumn/sum; + meanchi = sumchi/sum; + // + sdensity = sumdens2/sum-mdensity*mdensity; + sdensity = TMath::Sqrt(sdensity); + // + smeann = sumn2/sum-meann*meann; + smeann = TMath::Sqrt(smeann); + // + smeanchi = sumchi2/sum - meanchi*meanchi; + smeanchi = TMath::Sqrt(smeanchi); + } - Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift(); - // Double_t cs=cos(alpha), sn=sin(alpha); - Int_t row0 = (i1+i2)/2; - Int_t drow = (i1-i2)/2; - const AliTPCRow& kr0=fSectors[sec][row0]; - const AliTPCRow& krm=fSectors[sec][row0-1]; - const AliTPCRow& krp=fSectors[sec][row0+1]; - AliTPCRow * kr=0; - AliTPCpolyTrack polytrack; - Int_t nclusters=fSectors[sec][row0]; + //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC + // + for (Int_t i=0; iUncheckedAt(i); + if (!pt) { + continue; + } + if (pt->GetBSigned()) continue; + if (pt->GetBConstrain()) continue; + //if (!(pt->IsActive())) continue; + /* + Int_t found,foundable,shared; + pt->GetClusterStatistic(0,160,found, foundable,shared); + if (shared/float(found)>0.3) { + if (shared/float(found)>0.9 ){ + //delete arr->RemoveAt(i); + } + continue; + } + */ + Bool_t isok =kFALSE; + if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60) + isok = kTRUE; + if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7)) + isok =kTRUE; + if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1) + isok =kTRUE; + if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.)) + isok =kTRUE; + + if (isok) + for (Int_t i=0; i<160; i++) { + Int_t index=pt->GetClusterIndex2(i); + if (index<0) continue; + AliTPCclusterMI *c= pt->GetClusterPointer(i); + if (!c) continue; + //if (!(c->IsUsed(10))) c->Use(); + c->Use(10); + } + } + + + // + Double_t maxchi = meanchi+2.*smeanchi; + + for (Int_t i=0; iUncheckedAt(i); + if (!pt) { + continue; + } + //if (!(pt->IsActive())) continue; + if (pt->GetBSigned()) continue; + Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters(); + if (chi>maxchi) continue; + + Float_t bfactor=1; + Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable()); + + //sign only tracks with enoug big density at the beginning + + if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65); + if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens)) + minn=0; + + if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chiGetNumberOfClusters(); + pt->SetBSigned(kTRUE); + for (Int_t i=0; i<160; i++) { + + Int_t index=pt->GetClusterIndex2(i); + if (index<0) continue; + AliTPCclusterMI *c= pt->GetClusterPointer(i); + if (!c) continue; + // if (!(c->IsUsed(10))) c->Use(); + c->Use(10); + } + } + } + // gLastCheck = nseed; + // arr->Compress(); + if (fDebug>0){ + timer.Print(); + } +} + + +void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const +{ + // stop not active tracks + // take th1 as threshold for number of founded to number of foundable on last 10 active rows + // take th2 as threshold for number of founded to number of foundable on last 20 active rows + Int_t nseed = arr->GetEntriesFast(); + // + for (Int_t i=0; iUncheckedAt(i); + if (!pt) { + continue; + } + if (!(pt->IsActive())) continue; + StopNotActive(pt,row0,th0, th1,th2); + } +} + + + +void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1, + Float_t th2) const +{ + // stop not active tracks + // take th1 as threshold for number of founded to number of foundable on last 10 active rows + // take th2 as threshold for number of founded to number of foundable on last 20 active rows + Int_t sumgood1 = 0; + Int_t sumgood2 = 0; + Int_t foundable = 0; + Int_t maxindex = seed->GetLastPoint(); //last foundable row + if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) { + seed->Desactivate(10) ; + return; + } + + for (Int_t i=row0; iGetClusterIndex2(i); + if (index!=-1) foundable++; + //if (!c) continue; + if (foundable<=30) sumgood1++; + if (foundable<=50) { + sumgood2++; + } + else{ + break; + } + } + if (foundable>=30.){ + if (sumgood1<(th1*30.)) seed->Desactivate(10); + } + if (foundable>=50) + if (sumgood2<(th2*50.)) seed->Desactivate(10); +} + + +Int_t AliTPCtrackerMI::RefitInward(AliESD *event) +{ + // + // back propagation of ESD tracks + // + //return 0; + fEvent = event; + ReadSeeds(event,2); + 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 (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->GetNFoundable()); + Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3); + Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(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; +} + + +Int_t AliTPCtrackerMI::PropagateBack(AliESD *event) +{ + // + // back propagation of ESD tracks + // + + fEvent = event; + fIteration = 1; + 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 + if (seed->GetNumberOfClusters()>15){ + esd->UpdateTrackParams(seed,AliESDtrack::kTPCout); + esd->SetTPCPoints(seed->GetPoints()); + esd->SetTPCPointsF(seed->GetNFoundable()); + Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3); + Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25; + Float_t dedx = seed->GetdEdx(); + esd->SetTPCsignal(dedx, sdedx, ndedx); + ntracks++; + Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06 + // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number + (*fDebugStreamer)<<"Cback"<< + "Tr0.="<GetEntriesFast(); + for (Int_t i=0;iAt(i); + if (seed) delete fSeeds->RemoveAt(i); + } + delete fSeeds; + + fSeeds =0; +} + +void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction) +{ + // + //read seeds from the event + + Int_t nentr=event->GetNumberOfTracks(); + if (fDebug>0){ + Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr); + } + if (fSeeds) + DeleteSeeds(); + if (!fSeeds){ + fSeeds = new TObjArray(nentr); + } + UnsignClusters(); + // Int_t ntrk=0; + for (Int_t i=0; iGetTrack(i); + ULong_t status=esd->GetStatus(); + if (!(status&AliESDtrack::kTPCin)) continue; + AliTPCtrack t(*esd); + 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(10.); + if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.); + 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(10.); + } + + // + // + // rotate to the local coordinate system + // + 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)) { + delete seed; + continue; + } + seed->SetESD(esd); + // 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->SetClusterPointer(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); + } +} + + + +//_____________________________________________________________________________ +void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4], + Float_t deltay, Int_t ddsec) { + //----------------------------------------------------------------- + // This function creates track seeds. + // SEEDING WITH VERTEX CONSTRAIN + //----------------------------------------------------------------- + // cuts[0] - fP4 cut + // cuts[1] - tan(phi) cut + // cuts[2] - zvertex cut + // cuts[3] - fP3 cut + Int_t nin0 = 0; + Int_t nin1 = 0; + Int_t nin2 = 0; + Int_t nin = 0; + Int_t nout1 = 0; + Int_t nout2 = 0; + + Double_t x[5], c[15]; + // Int_t di = i1-i2; + // + AliTPCseed * seed = new AliTPCseed; + Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift(); + Double_t cs=cos(alpha), sn=sin(alpha); + // + // Double_t x1 =fOuterSec->GetX(i1); + //Double_t xx2=fOuterSec->GetX(i2); + + Double_t x1 =GetXrow(i1); + Double_t xx2=GetXrow(i2); + + Double_t x3=GetX(), y3=GetY(), z3=GetZ(); + + Int_t imiddle = (i2+i1)/2; //middle pad row index + Double_t xm = GetXrow(imiddle); // radius of middle pad-row + const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row + // + Int_t ns =sec; + + const AliTPCRow& kr1=GetRow(ns,i1); + Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5; + Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5; + + // + // change cut on curvature if it can't reach this layer + // maximal curvature set to reach it + Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3)); + if (dvertexmax*0.5*cuts[0]>0.85){ + cuts[0] = 0.85/(dvertexmax*0.5+1.); + } + Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut + + // Int_t ddsec = 1; + if (deltay>0) ddsec = 0; + // loop over clusters + for (Int_t is=0; is < kr1; is++) { + // + if (kr1[is]->IsUsed(10)) continue; + Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ(); + //if (TMath::Abs(y1)>ymax) continue; + + if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge + + // find possible directions + Float_t anglez = (z1-z3)/(x1-x3); + Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z + // + // + //find rotation angles relative to line given by vertex and point 1 + Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3); + Double_t dvertex = TMath::Sqrt(dvertex2); + Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3)); + Double_t cs13 = cos(-angle13), sn13 = sin(-angle13); + + // + // loop over 2 sectors + Int_t dsec1=-ddsec; + Int_t dsec2= ddsec; + if (y1<0) dsec2= 0; + if (y1>0) dsec1= 0; + + Double_t dddz1=0; // direction of delta inclination in z axis + Double_t dddz2=0; + if ( (z1-z3)>0) + dddz1 =1; + else + dddz2 =1; + // + for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){ + Int_t sec2 = sec + dsec; + // + // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2]; + //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle]; + AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2); + AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle); + Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0); + Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2); + + // rotation angles to p1-p3 + Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex; + Double_t x2, y2, z2; + // + // Double_t dymax = maxangle*TMath::Abs(x1-xx2); + + // + Double_t dxx0 = (xx2-x3)*cs13r; + Double_t dyy0 = (xx2-x3)*sn13r; + for (Int_t js=index1; js < index2; js++) { + const AliTPCclusterMI *kcl = kr2[js]; + if (kcl->IsUsed(10)) continue; + // + //calcutate parameters + // + Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r; + // stright track + if (TMath::Abs(yy0)<0.000001) continue; + Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r; + Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0; + Double_t r02 = (0.25+y0*y0)*dvertex2; + //curvature (radius) cut + if (r020) c0*=-1.; + + + //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5); + //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5); + Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5); + Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5); + // + // + Double_t z0 = kcl->GetZ(); + Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0; + if (TMath::Abs(zzzz2-z0)>0.5) continue; + nin1++; + // + Double_t dip = (z1-z0)*c0/dfi1; + Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0; + // + y2 = kcl->GetY(); + if (dsec==0){ + x2 = xx2; + z2 = kcl->GetZ(); + } + else + { + // rotation + z2 = kcl->GetZ(); + x2= xx2*cs-y2*sn*dsec; + y2=+xx2*sn*dsec+y2*cs; + } + + x[0] = y1; + x[1] = z1; + x[2] = x0; + x[3] = dip; + x[4] = c0; + // + // + // do we have cluster at the middle ? + Double_t ym,zm; + GetProlongation(x1,xm,x,ym,zm); + UInt_t dummy; + AliTPCclusterMI * cm=0; + if (TMath::Abs(ym)-ymaxm<0){ + cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy); + if ((!cm) || (cm->IsUsed(10))) { + continue; + } + } + else{ + // rotate y1 to system 0 + // get state vector in rotated system + Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0; + Double_t xr2 = x0*cs+yr1*sn*dsec; + Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0}; + // + GetProlongation(xx2,xm,xr,ym,zm); + if (TMath::Abs(ym)-ymaxm<0){ + cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy); + if ((!cm) || (cm->IsUsed(10))) { + continue; + } + } + } + + + Double_t dym = 0; + Double_t dzm = 0; + if (cm){ + dym = ym - cm->GetY(); + dzm = zm - cm->GetZ(); + } + nin2++; + + + // + // + Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.; + Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.; + //Double_t sy3=400*3./12., sy=0.1, sz=0.1; + Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1; + //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1; + + Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy; + Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy; + Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy; + Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; + Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; + Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; + + Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy; + Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz; + Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy; + Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz; + + 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; + + // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue; + + UInt_t index=kr1.GetIndex(is); + AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index); + + track->SetIsSeeding(kTRUE); + track->SetSeed1(i1); + track->SetSeed2(i2); + track->SetSeedType(3); + + + //if (dsec==0) { + FollowProlongation(*track, (i1+i2)/2,1); + Int_t foundable,found,shared; + track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE); + if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){ + seed->Reset(); + seed->~AliTPCseed(); + continue; + } + //} + + nin++; + FollowProlongation(*track, i2,1); + + + //Int_t rc = 1; + track->SetBConstrain(1); + // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position + track->SetLastPoint(i1); // first cluster in track position + track->SetFirstPoint(track->GetLastPoint()); + + if (track->GetNumberOfClusters()<(i1-i2)*0.5 || + track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || + track->GetNShared()>0.4*track->GetNumberOfClusters() ) { + seed->Reset(); + seed->~AliTPCseed(); + continue; + } + nout1++; + // Z VERTEX CONDITION + Double_t zv, bz=GetBz(); + if ( !track->GetZAt(0.,bz,zv) ) continue; + if (TMath::Abs(zv-z3)>cuts[2]) { + FollowProlongation(*track, TMath::Max(i2-20,0)); + if ( !track->GetZAt(0.,bz,zv) ) continue; + if (TMath::Abs(zv-z3)>cuts[2]){ + FollowProlongation(*track, TMath::Max(i2-40,0)); + if ( !track->GetZAt(0.,bz,zv) ) continue; + if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){ + // make seed without constrain + AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.); + FollowProlongation(*track2, i2,1); + track2->SetBConstrain(kFALSE); + track2->SetSeedType(1); + arr->AddLast(track2); + seed->Reset(); + seed->~AliTPCseed(); + continue; + } + else{ + seed->Reset(); + seed->~AliTPCseed(); + continue; + + } + } + } + + track->SetSeedType(0); + arr->AddLast(track); + seed = new AliTPCseed; + nout2++; + // don't consider other combinations + if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8) + break; + } + } + } + if (fDebug>3){ + Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2); + } + delete seed; +} + + +void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4], + Float_t deltay) { + + + + //----------------------------------------------------------------- + // This function creates track seeds. + //----------------------------------------------------------------- + // cuts[0] - fP4 cut + // cuts[1] - tan(phi) cut + // cuts[2] - zvertex cut + // cuts[3] - fP3 cut + + + Int_t nin0 = 0; + Int_t nin1 = 0; + Int_t nin2 = 0; + Int_t nin = 0; + Int_t nout1 = 0; + Int_t nout2 = 0; + Int_t nout3 =0; + Double_t x[5], c[15]; + // + // make temporary seed + AliTPCseed * seed = new AliTPCseed; + Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift(); + // Double_t cs=cos(alpha), sn=sin(alpha); + // + // + + // first 3 padrows + Double_t x1 = GetXrow(i1-1); + const AliTPCRow& kr1=GetRow(sec,i1-1); + Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5; + // + Double_t x1p = GetXrow(i1); + const AliTPCRow& kr1p=GetRow(sec,i1); + // + Double_t x1m = GetXrow(i1-2); + const AliTPCRow& kr1m=GetRow(sec,i1-2); + + // + //last 3 padrow for seeding + AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7); + Double_t x3 = GetXrow(i1-7); + // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5; + // + AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6); + Double_t x3p = GetXrow(i1-6); + // + AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8); + Double_t x3m = GetXrow(i1-8); + + // + // + // middle padrow + Int_t im = i1-4; //middle pad row index + Double_t xm = GetXrow(im); // radius of middle pad-row + const AliTPCRow& krm=GetRow(sec,im); //middle pad -row + // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5; + // + // + Double_t deltax = x1-x3; + Double_t dymax = deltax*cuts[1]; + Double_t dzmax = deltax*cuts[3]; + // + // loop over clusters + for (Int_t is=0; is < kr1; is++) { + // + if (kr1[is]->IsUsed(10)) continue; + Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ(); + // + if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge + // + Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0); + Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3); + // + Double_t y3, z3; + // + // + UInt_t index; + for (Int_t js=index1; js < index2; js++) { + const AliTPCclusterMI *kcl = kr3[js]; + if (kcl->IsUsed(10)) continue; + y3 = kcl->GetY(); + // apply angular cuts + if (TMath::Abs(y1-y3)>dymax) continue; + x3 = x3; + z3 = kcl->GetZ(); + if (TMath::Abs(z1-z3)>dzmax) continue; + // + Double_t angley = (y1-y3)/(x1-x3); + Double_t anglez = (z1-z3)/(x1-x3); + // + Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5; + Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5; + // + Double_t yyym = angley*(xm-x1)+y1; + Double_t zzzm = anglez*(xm-x1)+z1; + + const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index); + if (!kcm) continue; + if (kcm->IsUsed(10)) continue; + + erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5; + errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5; + // + // + // + Int_t used =0; + Int_t found =0; + // + // look around first + const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1, + anglez*(x1m-x1)+z1, + erry,errz,index); + // + if (kc1m){ + found++; + if (kc1m->IsUsed(10)) used++; + } + const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1, + anglez*(x1p-x1)+z1, + erry,errz,index); + // + if (kc1p){ + found++; + if (kc1p->IsUsed(10)) used++; + } + if (used>1) continue; + if (found<1) continue; + + // + // look around last + const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3, + anglez*(x3m-x3)+z3, + erry,errz,index); + // + if (kc3m){ + found++; + if (kc3m->IsUsed(10)) used++; + } + else + continue; + const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3, + anglez*(x3p-x3)+z3, + erry,errz,index); + // + if (kc3p){ + found++; + if (kc3p->IsUsed(10)) used++; + } + else + continue; + if (used>1) continue; + if (found<3) continue; + // + Double_t x2,y2,z2; + x2 = xm; + y2 = kcm->GetY(); + z2 = kcm->GetZ(); + // + + x[0]=y1; + x[1]=z1; + x[4]=F1(x1,y1,x2,y2,x3,y3); + //if (TMath::Abs(x[4]) >= cuts[0]) continue; + nin0++; + // + x[2]=F2(x1,y1,x2,y2,x3,y3); + nin1++; + // + x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]); + //if (TMath::Abs(x[3]) > cuts[3]) continue; + nin2++; + // + // + Double_t sy1=0.1, sz1=0.1; + Double_t sy2=0.1, sz2=0.1; + Double_t sy3=0.1, sy=0.1, sz=0.1; + + Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy; + Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy; + Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy; + Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; + Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; + Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; + + Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy; + Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz; + Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy; + Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz; + + 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; + + // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue; + + UInt_t index=kr1.GetIndex(is); + AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index); + + track->SetIsSeeding(kTRUE); + + nin++; + FollowProlongation(*track, i1-7,1); + if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 || + track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){ + seed->Reset(); + seed->~AliTPCseed(); + continue; + } + nout1++; + nout2++; + //Int_t rc = 1; + FollowProlongation(*track, i2,1); + track->SetBConstrain(0); + track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position + track->SetFirstPoint(track->GetLastPoint()); + + if (track->GetNumberOfClusters()<(i1-i2)*0.5 || + track->GetNumberOfClusters()GetNFoundable()*0.7 || + track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) { + seed->Reset(); + seed->~AliTPCseed(); + continue; + } + + { + FollowProlongation(*track, TMath::Max(i2-10,0),1); + AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9); + FollowProlongation(*track2, i2,1); + track2->SetBConstrain(kFALSE); + track2->SetSeedType(4); + arr->AddLast(track2); + seed->Reset(); + seed->~AliTPCseed(); + } + + + //arr->AddLast(track); + //seed = new AliTPCseed; + nout3++; + } + } + + if (fDebug>3){ + Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3); + } + delete seed; +} + + +//_____________________________________________________________________________ +void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/, + Float_t deltay, Bool_t /*bconstrain*/) { + //----------------------------------------------------------------- + // This function creates track seeds - without vertex constraint + //----------------------------------------------------------------- + // cuts[0] - fP4 cut - not applied + // cuts[1] - tan(phi) cut + // cuts[2] - zvertex cut - not applied + // cuts[3] - fP3 cut + Int_t nin0=0; + Int_t nin1=0; + Int_t nin2=0; + Int_t nin3=0; + // Int_t nin4=0; + //Int_t nin5=0; + + + + Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift(); + // Double_t cs=cos(alpha), sn=sin(alpha); + Int_t row0 = (i1+i2)/2; + Int_t drow = (i1-i2)/2; + const AliTPCRow& kr0=fSectors[sec][row0]; + AliTPCRow * kr=0; + + AliTPCpolyTrack polytrack; + Int_t nclusters=fSectors[sec][row0]; + AliTPCseed * seed = new AliTPCseed; + + Int_t sumused=0; + Int_t cused=0; + Int_t cnused=0; + for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters + Int_t nfound =0; + Int_t nfoundable =0; + for (Int_t iter =1; iter<2; iter++){ //iterations + const AliTPCRow& krm=fSectors[sec][row0-iter]; + const AliTPCRow& krp=fSectors[sec][row0+iter]; + const AliTPCclusterMI * cl= kr0[is]; + + if (cl->IsUsed(10)) { + cused++; + } + else{ + cnused++; + } + Double_t x = kr0.GetX(); + // Initialization of the polytrack + nfound =0; + nfoundable =0; + polytrack.Reset(); + // + Double_t y0= cl->GetY(); + Double_t z0= cl->GetZ(); + Float_t erry = 0; + Float_t errz = 0; + + Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5; + if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge + + erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6; + errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6; + polytrack.AddPoint(x,y0,z0,erry, errz); + + sumused=0; + if (cl->IsUsed(10)) sumused++; + + + Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter; + Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter; + // + x = krm.GetX(); + AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz); + if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) { + erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3; + errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3; + if (cl1->IsUsed(10)) sumused++; + polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz); + } + // + x = krp.GetX(); + AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz); + if (cl2) { + erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3; + errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3; + if (cl2->IsUsed(10)) sumused++; + polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz); + } + // + if (sumused>0) continue; + nin0++; + polytrack.UpdateParameters(); + // follow polytrack + roadz = 1.2; + roady = 1.2; + // + Double_t yn,zn; + nfoundable = polytrack.GetN(); + nfound = nfoundable; + // + for (Int_t ddrow = iter+1; ddrowGetX(); + Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5; + polytrack.GetFitPoint(xn,yn,zn); + if (TMath::Abs(yn)>ymax) continue; + nfoundable++; + AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz); + if (cln) { + Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ())); + if (distGetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow)); + errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow)); + if (cln->IsUsed(10)) { + // printf("used\n"); + sumused++; + erry*=2; + errz*=2; + } + */ + erry=0.1; + errz=0.1; + polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz); + nfound++; + } + } + } + if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break; + polytrack.UpdateParameters(); + } + } + if ( (sumused>3) || (sumused>0.5*nfound)) { + //printf("sumused %d\n",sumused); + continue; + } + nin1++; + Double_t dy,dz; + polytrack.GetFitDerivation(kr0.GetX(),dy,dz); + AliTPCpolyTrack track2; + + polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3); + if (track2.GetN()<0.5*nfoundable) continue; + nin2++; + + if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) { + // + // test seed with and without constrain + for (Int_t constrain=0; constrain<=0;constrain++){ + // add polytrack candidate + + Double_t x[5], c[15]; + Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3; + track2.GetBoundaries(x3,x1); + x2 = (x1+x3)/2.; + track2.GetFitPoint(x1,y1,z1); + track2.GetFitPoint(x2,y2,z2); + track2.GetFitPoint(x3,y3,z3); + // + //is track pointing to the vertex ? + Double_t x0,y0,z0; + x0=0; + polytrack.GetFitPoint(x0,y0,z0); + + if (constrain) { + x2 = x3; + y2 = y3; + z2 = z3; + + x3 = 0; + y3 = 0; + z3 = 0; + } + x[0]=y1; + x[1]=z1; + x[4]=F1(x1,y1,x2,y2,x3,y3); + + // if (TMath::Abs(x[4]) >= cuts[0]) continue; // + x[2]=F2(x1,y1,x2,y2,x3,y3); + + //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue; + //x[3]=F3(x1,y1,x2,y2,z1,z2); + x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]); + //if (TMath::Abs(x[3]) > cuts[3]) continue; + + + Double_t sy =0.1, sz =0.1; + Double_t sy1=0.02, sz1=0.02; + Double_t sy2=0.02, sz2=0.02; + Double_t sy3=0.02; + + if (constrain){ + sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1; + } + + Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy; + Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy; + Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy; + Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; + Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; + Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; + + Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy; + Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz; + Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy; + Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz; + + + 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(x1); + Int_t row1 = GetRowNumber(x1); + + UInt_t index=0; + //kr0.GetIndex(is); + AliTPCseed *track=new (seed) AliTPCseed(x1,sec*alpha+shift,x,c,index); + track->SetIsSeeding(kTRUE); + Int_t rc=FollowProlongation(*track, i2); + if (constrain) track->SetBConstrain(1); + else + track->SetBConstrain(0); + track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position + track->SetFirstPoint(track->GetLastPoint()); + + if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 || + track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || + track->GetNShared()>0.4*track->GetNumberOfClusters()) { + //delete track; + seed->Reset(); + seed->~AliTPCseed(); + } + else { + arr->AddLast(track); + seed = new AliTPCseed; + } + nin3++; + } + } // if accepted seed + } + if (fDebug>3){ + Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3); + } + delete seed; +} + + +AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2) +{ + // + // + //reseed using track points + Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0 + Int_t p1 = int(r1*track->GetNumberOfClusters()); + Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point + Int_t pp2=0; + Double_t x0[3],x1[3],x2[3]; + 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=0, sec1=0, sec2=0; + Int_t index=-1; + Int_t clindex; + for (Int_t i=0;i<160;i++){ + if (track->GetClusterPointer(i)){ + index++; + AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i); + if ( (indexGetX()>1){ + clindex = track->GetClusterIndex2(i); + if (clindex>0){ + x0[0] = trpoint->GetX(); + x0[1] = trpoint->GetY(); + x0[2] = trpoint->GetZ(); + sec0 = ((clindex&0xff000000)>>24)%18; + } + } + } + + if ( (indexGetX()>1)){ + clindex = track->GetClusterIndex2(i); + if (clindex>0){ + x1[0] = trpoint->GetX(); + x1[1] = trpoint->GetY(); + x1[2] = trpoint->GetZ(); + sec1 = ((clindex&0xff000000)>>24)%18; + } + } + if ( (indexGetX()>1)){ + clindex = track->GetClusterIndex2(i); + if (clindex>0){ + x2[0] = trpoint->GetX(); + x2[1] = trpoint->GetY(); + x2[2] = trpoint->GetZ(); + sec2 = ((clindex&0xff000000)>>24)%18; + pp2 = i; + } + } + } + } + + Double_t alpha, cs,sn, xx2,yy2; + // + alpha = (sec1-sec2)*fSectors->GetAlpha(); + cs = TMath::Cos(alpha); + sn = TMath::Sin(alpha); + xx2= x1[0]*cs-x1[1]*sn; + yy2= x1[0]*sn+x1[1]*cs; + x1[0] = xx2; + x1[1] = yy2; + // + alpha = (sec0-sec2)*fSectors->GetAlpha(); + cs = TMath::Cos(alpha); + sn = TMath::Sin(alpha); + xx2= x0[0]*cs-x0[1]*sn; + yy2= x0[0]*sn+x0[1]*cs; + x0[0] = xx2; + x0[1] = yy2; + // + // + // + Double_t x[5],c[15]; + // + x[0]=x2[1]; + x[1]=x2[2]; + x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]); + // if (x[4]>1) return 0; + x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]); + x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]); + //if (TMath::Abs(x[3]) > 2.2) return 0; + //if (TMath::Abs(x[2]) > 1.99) return 0; + // + Double_t sy =0.1, sz =0.1; + // + Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2(); + Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2(); + Double_t sy3=0.01+track->GetSigmaY2(); + // + Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy; + Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy; + Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy; + Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy; + Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy; + Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy; + // + Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy; + Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz; + Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy; + Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz; + + + 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(x2[0]); + AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + // Double_t y0,z0,y1,z1, y2,z2; + //seed->GetProlongation(x0[0],y0,z0); + // seed->GetProlongation(x1[0],y1,z1); + //seed->GetProlongation(x2[0],y2,z2); + // seed =0; + seed->SetLastPoint(pp2); + seed->SetFirstPoint(pp2); + + + return seed; +} + + +AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2) +{ + // + // + //reseed using founded clusters + // + // Find the number of clusters + Int_t nclusters = 0; + for (Int_t irow=0;irow<160;irow++){ + if (track->GetClusterIndex(irow)>0) nclusters++; + } + // + Int_t ipos[3]; + ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster + ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); // + ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point + // + // + Double_t xyz[3][3]; + Int_t row[3],sec[3]={0,0,0}; + // + // find track row position at given ratio of the length + Int_t index=-1; + for (Int_t irow=0;irow<160;irow++){ + if (track->GetClusterIndex2(irow)<0) continue; + index++; + for (Int_t ipoint=0;ipoint<3;ipoint++){ + if (index<=ipos[ipoint]) row[ipoint] = irow; + } + } + // + //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; + } + sec[ipoint] = ((clindex&0xff000000)>>24)%18; + xyz[ipoint][0] = GetXrow(row[ipoint]); + xyz[ipoint][1] = cl->GetY(); + xyz[ipoint][2] = cl->GetZ(); + } + // + // + // Calculate seed state vector and covariance matrix + + 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; + + + 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(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + seed->SetLastPoint(row[2]); + seed->SetFirstPoint(row[2]); + return seed; +} + + +AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward) +{ + // + // + //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; + } + } + for (Int_t irow=160;irow>r0;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; + } + } + // + } + 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; + } + } + } + // + 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; + } + 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{ + xyz[ipoint][1] = cl->GetY(); + xyz[ipoint][2] = cl->GetZ(); + } + } + // + // + // + // + // Calculate seed state vector and covariance matrix + + 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; + + + 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(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + seed->SetLastPoint(row[2]); + seed->SetFirstPoint(row[2]); + for (Int_t i=row[0];iSetClusterIndex(i, track->GetClusterIndex(i)); + } + + return seed; +} + +void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) +{ + // + // find kinks + // + // + + 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]; + AliKink *kink = new AliKink(); + 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->SetCircular(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); + } + // + // + 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->GetNumberOfClusters()<40) continue; + if (TMath::Abs(1./track0->GetC())>200) continue; + for (Int_t i1=i0+1;i1At(i1); + if (!track1) continue; + if (track1->GetNumberOfClusters()<40) continue; + if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue; + if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; + if (TMath::Abs(1./track1->GetC())>200) continue; + if (track1->Get1Pt()*track0->Get1Pt()>0) continue; + if (track1->GetTgl()*track0->GetTgl()>0) continue; + if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue; + if (track0->GetBConstrain()&&TMath::Abs(track1->Get1Pt())Get1Pt())) continue; //returning - lower momenta + if (track1->GetBConstrain()&&TMath::Abs(track0->Get1Pt())Get1Pt())) 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->Get1Pt())Get1Pt())){ + track0->SetCircular(track0->GetCircular()+1); + track1->SetCircular(track1->GetCircular()+2); + } + else{ + track1->SetCircular(track1->GetCircular()+1); + track0->SetCircular(track0->GetCircular()+2); + } + } + if (sign&&AliTPCReconstructor::StreamLevel()>1){ + //debug stream + Int_t lab0=track0->GetLabel(); + Int_t lab1=track1->GetLabel(); + cstream<<"Curling"<< + "lab0="<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->GetP(); + 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); + } - for (Int_t is=0; is < nclusters; is++) { - const AliTPCclusterMI * cl= kr0[is]; - Double_t x = kr0.GetX(); - - // Initialization of the polytrack - polytrack.Reset(); - - Double_t y0= cl->GetY(); - Double_t z0= cl->GetZ(); - polytrack.AddPoint(x,y0,z0); - Float_t roady = 5*TMath::Sqrt(cl->GetSigmaY2()+0.2); - Float_t roadz = 5*TMath::Sqrt(cl->GetSigmaZ2()+0.2); - // - x = krm.GetX(); - cl = krm.FindNearest(y0,z0,roady,roadz); - if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),roady,roadz); - // - x = krp.GetX(); - cl = krp.FindNearest(y0,z0,roady,roadz); - if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),cl->GetSigmaY2()+0.05,cl->GetSigmaZ2()+0.05); - // - polytrack.UpdateParameters(); - // follow polytrack - roadz = 0.6; - roady = 0.6; - // - Double_t yn,zn; - Int_t nfoundable = polytrack.GetN(); - Int_t nfound = nfoundable; - for (Int_t ddrow = 2; ddrowGetX(); - Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone; - polytrack.GetFitPoint(xn,yn,zn); - if (TMath::Abs(yn)>ymax) continue; - nfoundable++; - AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz); - if (cln) { - polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),cln->GetSigmaY2()+0.05,cln->GetSigmaZ2()+0.05); - nfound++; + 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->GetP()<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->GetSigmaSnp2()+track0->GetSigmaTgl2(); + criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2(); + 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->GetClusterPointer(row)){ + AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row); + shapesum+=point->GetSigmaY()+point->GetSigmaZ(); + sum++; } + if (ktrack1->GetClusterPointer(row)){ + AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row); + shapesum+=point->GetSigmaY()+point->GetSigmaZ(); + sum++; + } } - polytrack.UpdateParameters(); - if (nfound<0.45*nfoundable) break; - } - if ((nfound>0.5*nfoundable) &&( nfoundable>0.4*(i1-i2))) { - // add polytrack candidate - Double_t x[5], c[15]; - Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3; - polytrack.GetBoundaries(x3,x1); - x2 = (x1+x3)/2.; - polytrack.GetFitPoint(x1,y1,z1); - polytrack.GetFitPoint(x2,y2,z2); - polytrack.GetFitPoint(x3,y3,z3); - // - //is track pointing to the vertex ? - Double_t x0,y0,z0; - x0=0; - polytrack.GetFitPoint(x0,y0,z0); - if ( (TMath::Abs(z0-GetZ())<10) && (TMath::Abs(y0-GetY())<5)){ //if yes apply vertex constraint - // x3 = 0; - //y3 = GetY(); - //z3 = GetZ(); + if (sum<4){ + kink->SetShapeFactor(-1.); + } + else{ + kink->SetShapeFactor(shapesum/sum); + } + // esd->AddKink(kink); + kinks->AddLast(kink); + kink = new AliKink; + 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->GetNumberOfClusters()+ktrack1->GetNumberOfClusters(); + // + // 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 (lastGetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25; + AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE); + // + // + Int_t first = row0+drow; + if (first>130) first=130; + if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-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->GetNumberOfClusters()+seed1->GetNumberOfClusters(); + 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.GetClusterIndex(i)>0){ + both++; + bothm++; + if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){ + same++; + samem++; + } + } } - x[0]=y1; - x[1]=z1; - x[4]=f1(x1,y1,x2,y2,x3,y3); - if (TMath::Abs(x[4]) >= 0.05) continue; //MI change - x[2]=f2(x1,y1,x2,y2,x3,y3); - //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue; - x[3]=f3(x1,y1,x2,y2,z1,z2); - if (TMath::Abs(x[3]) > 1.2) continue; - if (TMath::Abs(x[2]) > 0.99) continue; - // Double_t a=asin(x[2]); - + for (Int_t i=row0;i<158;i++){ + if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ + both++; + bothd++; + if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(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.GetNumberOfClusters()+daughter0.GetNumberOfClusters(); + Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters(); + 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]); + } + } + } + } - Double_t sy=1.5, sz=1.5; - Double_t sy1=1.5, sz1=1.5; - Double_t sy2=1.3, sz2=1.3; - Double_t sy3=1.5; - //sz3=1.5; - //Double_t sy3=400*3./12., sy=0.1, sz=0.1; - // Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1; - //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1; - - Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy; - Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy; - Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy; - Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy; - Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy; - Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy; - Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy; - Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz; - Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy; - Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz; - - 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; - - UInt_t index=0; - //kr0.GetIndex(is); - AliTPCseed *track=new AliTPCseed(index, x, c, x1, sec*alpha+shift); - track->fStopped =kFALSE; - track->fIsSeeding = kTRUE; - Int_t rc=FollowProlongation(*track, i2); - track->fLastPoint = i1; // first cluster in track position - if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track; - else arr->AddLast(track); + + 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->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink + if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){ + new (ktrack0) AliTPCseed(mothers[indexes[i]]); + new (ktrack1) AliTPCseed(daughters[indexes[i]]); + } } + // + ktrack0->SetKinkIndex(usage[index0],-(index+1)); + ktrack1->SetKinkIndex(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->GetKinkIndex(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->GetPt()<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->GetFirstPoint();iclGetLastPoint(); icl++){ + if (track0->GetClusterPointer(icl)!=0){ + all++; + if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++; + } + } + if (Float_t(shared+1)/Float_t(nall+1)>0.5) { + delete array->RemoveAt(i); + continue; + } + // + if (track0->GetKinkIndex(0)!=0) continue; + if (track0->GetNumberOfClusters()<80) continue; + + AliTPCseed *pmother = new AliTPCseed(); + AliTPCseed *pdaughter = new AliTPCseed(); + AliKink *pkink = new AliKink; + + AliTPCseed & mother = *pmother; + AliTPCseed & daughter = *pdaughter; + AliKink & kink = *pkink; + if (CheckKinkPoint(track0,mother,daughter, kink)){ + if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) { + delete pmother; + delete pdaughter; + delete pkink; + continue; //too short tracks + } + if (mother.GetPt()<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.SetKinkIndex(0,-(index+1)); + daughter.SetKinkIndex(0,index+1); + if (mother.GetNumberOfClusters()>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.GetClusterPointer(icl)) daughter.GetClusterPointer(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->GetC())*(100*track->GetC())); + if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5; + cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-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->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) { + if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut + } + if (track->TPCrPID(4)>0.5) { + if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut + } + if (track->TPCrPID(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()}; + AliV0 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.="<Get1Pt()<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->GetCircular()+track1->GetCircular()>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 (pointAngleTPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate + Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle(); + //continue; + if (vertex.GetV0CosineOfPointingAngle()2&&(!isGamma)) continue; // point angle cut + if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // 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); + //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue; + if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue; + Float_t densb0=0,densb1=0,densa0=0,densa1=0; + Int_t row0 = GetRowNumber(vertex.GetRr()); + if (row0>15){ + //Bo: if (vertex.GetDist2()>0.2) continue; + if (vertex.GetDcaV0Daughters()>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; + } +//Bo: vertex.SetLab(0,track0->GetLabel()); +//Bo: 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); +//Bo: vertex.SetStatus(1); // TPC v0 candidate + vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate +//Bo: vertex.SetRp(track0->TPCrPIDs()); +//Bo: vertex.SetRm(track1->TPCrPIDs()); + tpcv0s->AddLast(new AliESDv0(vertex)); + ncandidates++; + { + Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number + Double_t radiusm= (delta[0]0) { + Int_t lab0=track0->GetLabel(); + Int_t lab1=track1->GetLabel(); + Char_t c0=track0->GetCircular(); + Char_t c1=track1->GetCircular(); + cstream<<"V0"<< + "Event="<At(i); + quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle + // quality[i] /= (0.5+v0->GetDist2()); + // quality[i] *= v0->GetChi2After(); //density factor + + Int_t index0 = v0->GetIndex(0); + Int_t index1 = v0->GetIndex(1); + //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull + Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo: + + + + AliTPCseed * track0 = (AliTPCseed*)array->At(index0); + AliTPCseed * track1 = (AliTPCseed*)array->At(index1); + if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate + if (track0->TPCrPID(4)>0.9||track1->TPCrPID(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;} + // + AliESDv0 * v02 = v0; + if (accept){ + //Bo: v0->SetOrder(0,order0); + //Bo: v0->SetOrder(1,order1); + //Bo: v0->SetOrder(1,order0+order1); + v0->SetOnFlyStatus(kTRUE); + Int_t index = esd->AddV0(v0); + v02 = esd->GetV0(index); + v0indexes0[order0]=index; + v0indexes1[order1]=index; + naccepted++; + } + { + Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number + if (AliTPCReconstructor::StreamLevel()>0) { + Int_t lab0=track0->GetLabel(); + Int_t lab1=track1->GetLabel(); + cstream<<"V02"<< + "Event="<240.) continue; + if (TMath::Abs(kinks[irow].GetR())<100.) continue; + // + Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance()); + normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+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].GetNumberOfClusters(),0); + kink.SetTPCncls(param1[index].GetNumberOfClusters(),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 &knk) +{ + // + // 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 + AliKink &kink=(AliKink &)knk; + + Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2; + Int_t first = seed->GetFirstPoint(); + Int_t last = seed->GetLastPoint(); + 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->GetLastPoint()-20); + seed1->Reset(kTRUE); + FollowProlongation(*seed1,158); + seed1->Reset(kTRUE); + last = seed1->GetLastPoint(); + // + AliTPCseed *seed0 = new AliTPCseed(*seed); + seed0->Reset(kFALSE); + seed0->Reset(); + // + AliTPCseed param0[20]; // parameters along the track + AliTPCseed param1[20]; // parameters along the track + AliKink 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].GetX())); + 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(10.); + seed1->ResetCovariance(10.); + FollowProlongation(*seed0,0); + FollowBackProlongation(*seed1,158); + new (&mother) AliTPCseed(*seed0); // backup mother at position 0 + seed0->Reset(kFALSE); + seed1->Reset(kFALSE); + seed0->ResetCovariance(10.); + seed1->ResetCovariance(10.); + // + 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].GetX())); + if ( quality > maxchange){ + maxchange = quality; + index = irow; + // + } + } + // + // + if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<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].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1); + kink.SetTPCncls(param0[index].GetNumberOfClusters(),0); + kink.SetTPCncls(param1[index].GetNumberOfClusters(),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->GetFirstPoint()-fInnerSec->GetNRows()); + //last = + } + else + first = t->GetFirstPoint(); + // + 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) { //----------------------------------------------------------------- @@ -1725,7 +5480,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; @@ -1734,188 +5489,742 @@ Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) { return 0; } +Int_t AliTPCtrackerMI::Clusters2Tracks (AliESD *esd) +{ + // + if (fSeeds) DeleteSeeds(); + fEvent = esd; + Clusters2Tracks(); + if (!fSeeds) return 1; + FillESD(fSeeds); + return 0; + // +} + + //_____________________________________________________________________________ Int_t AliTPCtrackerMI::Clusters2Tracks() { //----------------------------------------------------------------- // This is a track finder. //----------------------------------------------------------------- - TTree* clustersTree = AliRunLoader::GetTreeR("TPC", kFALSE); - if (!clustersTree) { - Error("Clusters2Tracks", "no clusters found"); - return 1; + TDirectory *savedir=gDirectory; + TStopwatch timer; + + fIteration = 0; + fSeeds = Tracking(); + + if (fDebug>0){ + Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start(); } - fClustersArray.ConnectTree(clustersTree); + //activate again some tracks + for (Int_t i=0; iGetEntriesFast(); i++) { + AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt; + if (!pt) continue; + Int_t nc=t.GetNumberOfClusters(); + if (nc<20) { + delete fSeeds->RemoveAt(i); + continue; + } + CookLabel(pt,0.1); + if (pt->GetRemoval()==10) { + if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7) + pt->Desactivate(10); // make track again active + else{ + pt->Desactivate(20); + delete fSeeds->RemoveAt(i); + } + } + } + // + 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); - TTree* tracksTree = AliRunLoader::GetTreeT("TPC", kTRUE); - TTree& tracktree = *tracksTree; -// TTree seedtree("Seeds","Seeds"); - AliTPCtrack *iotrack=0; - AliTPCseed *ioseed=0; - tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0); - TStopwatch timer; + // + Int_t found = 0; + 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; + } + 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->GetRemoval()==10) )){ + found++; + if (fDebug>0) cerr<SetLab2(i); + } + else + delete fSeeds->RemoveAt(i); + } + + + //RemoveOverlap(fSeeds,0.99,7,kTRUE); + SignShared(fSeeds); + //RemoveUsed(fSeeds,0.9,0.9,6); + // + nseed=fSeeds->GetEntriesFast(); + found = 0; + 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; + } + t.SetUniqueID(i); + t.CookdEdx(0.02,0.6); + // CheckKinkPoint(&t,0.05); + //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){ + if ((pt->IsActive() || (pt->GetRemoval()==10) )){ + found++; + if (fDebug>0){ + cerr<SetLab2(i); + } + else + delete fSeeds->RemoveAt(i); + //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1); + //if (seed1){ + // FollowProlongation(*seed1,0); + // Int_t n = seed1->GetNumberOfClusters(); + // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC()); + // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters()); + // + //} + //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05); + + } + + SortTracks(fSeeds, 1); + + /* + fIteration = 1; + PrepareForBackProlongation(fSeeds,5.); + PropagateBack(fSeeds); + printf("Time for back propagation: \t");timer.Print();timer.Start(); + + fIteration = 2; - printf("Loading clusters \n"); - LoadClusters(); - printf("Time for loading clusters: \t");timer.Print();timer.Start(); + PrepareForProlongation(fSeeds,5.); + PropagateForward2(fSeeds); + + printf("Time for FORWARD propagation: \t");timer.Print();timer.Start(); + // RemoveUsed(fSeeds,0.7,0.7,6); + //RemoveOverlap(fSeeds,0.9,7,kTRUE); + + nseed=fSeeds->GetEntriesFast(); + found = 0; + 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; + } + t.CookdEdx(0.02,0.6); + // 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<RemoveAt(i); + pt->fLab2 = i; + } + */ + + // fNTracks = found; + if (fDebug>0){ + Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start(); + } + // + // cerr<<"Number of found tracks : "<<"\t"<cd(); + // UnloadClusters(); + // + return 0; +} + +void AliTPCtrackerMI::Tracking(TObjArray * arr) +{ + // + // tracking of the seeds + // fSectors = fOuterSec; - fN=fkNOS; + ParallelTracking(arr,150,63); + fSectors = fOuterSec; + ParallelTracking(arr,63,0); +} + +TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec) +{ + // + // + //tracking routine + TObjArray * arr = new TObjArray; + // + fSectors = fOuterSec; + TStopwatch timer; + timer.Start(); + for (Int_t sec=0;sec0){ + Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast()); + timer.Print(); + timer.Start(); + } + Tracking(arr); + if (fDebug>0){ + timer.Print(); + } + + return arr; +} + +TObjArray * AliTPCtrackerMI::Tracking() +{ + // + // + if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial(); + TStopwatch timer; + timer.Start(); + Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows(); + + TObjArray * seeds = new TObjArray; + TObjArray * arr=0; + + Int_t gap =20; + Float_t cuts[4]; + cuts[0] = 0.002; + cuts[1] = 1.5; + cuts[2] = 3.; + cuts[3] = 3.; + Float_t fnumber = 3.0; + Float_t fdensity = 3.0; + // + //find primaries + cuts[0]=0.0066; + for (Int_t delta = 0; delta<18; delta+=6){ + // + cuts[0]=0.0070; + cuts[1] = 1.5; + arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + for (Int_t i=2;i<6;i+=2){ + // seed high pt tracks + cuts[0]=0.0022; + cuts[1]=0.3; + arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + } + } + fnumber = 4; + fdensity = 4.; + // RemoveUsed(seeds,0.9,0.9,1); + // UnsignClusters(); + // SignClusters(seeds,fnumber,fdensity); + + //find primaries + cuts[0]=0.0077; + for (Int_t delta = 20; delta<120; delta+=10){ + // + // seed high pt tracks + cuts[0]=0.0060; + cuts[1]=0.3; + cuts[2]=6.; + arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + + cuts[0]=0.003; + cuts[1]=0.3; + cuts[2]=6.; + arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + } + + cuts[0] = 0.01; + cuts[1] = 2.0; + cuts[2] = 3.; + cuts[3] = 2.0; + fnumber = 2.; + fdensity = 2.; - //find track seeds - MakeSeedsAll(); - printf("Time for seeding: \t"); timer.Print();timer.Start(); - Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows(); - Int_t nrows=nlow+nup; + if (fDebug>0){ + Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast()); + timer.Print(); + timer.Start(); + } + // RemoveUsed(seeds,0.75,0.75,1); + //UnsignClusters(); + //SignClusters(seeds,fnumber,fdensity); - Int_t gap=Int_t(0.3*nrows); - Int_t i; - //RemoveOverlap(fSeeds,0.6,2); - Int_t nseed=fSeeds->GetEntriesFast(); - // outer sectors parallel tracking - ParallelTracking(fSectors->GetNRows()-gap-1,0); - printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start(); + // find secondaries - RemoveOverlap(fSeeds, 0.4,3); - printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start(); - //parallel tracking - fSectors = fInnerSec; - fN=fkNIS; + cuts[0] = 0.3; + cuts[1] = 1.5; + cuts[2] = 3.; + cuts[3] = 1.5; - ParallelTracking(fSectors->GetNRows()-1,0); - /* - ParallelTracking(fSectors->GetNRows()-1,2*fSectors->GetNRows()/3); - RemoveOverlap(fSeeds,0.4,5,kTRUE); - ParallelTracking(2*fSectors->GetNRows()/3-1,fSectors->GetNRows()/3); - RemoveOverlap(fSeeds,0.4,5,kTRUE); - ParallelTracking(fSectors->GetNRows()/3-1,0); - */ - printf("Number of tracks after inner tracking %d\n",fNtracks); - printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start(); - // - for (Int_t i=0;iGetEntriesFast();i++){ - AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i); - if (!pt) continue; - if (!pt->IsActive()) continue; - pt->PropagateTo(90.); + arr = Tracking(4,nup-1,nup-1-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + arr = Tracking(4,nup-2,nup-2-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + arr = Tracking(4,nup-3,nup-3-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + + + for (Int_t delta = 3; delta<30; delta+=5){ + // + cuts[0] = 0.3; + cuts[1] = 1.5; + cuts[2] = 3.; + cuts[3] = 1.5; + arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // } - RemoveOverlap(fSeeds,0.4,5,kTRUE); // remove overlap - shared points signed - RemoveUsed(fSeeds,0.4,6); - printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start(); + fnumber = 1; + fdensity = 1; // - // + // change cuts + fnumber = 2.; + fdensity = 2.; + cuts[0]=0.0080; + + // find secondaries + for (Int_t delta = 30; delta<90; delta+=10){ + // + cuts[0] = 0.3; + cuts[1] = 3.5; + cuts[2] = 3.; + cuts[3] = 3.5; + arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 ); + 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; + // + +} + + +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; - ioseed = (AliTPCseed*)(fSeeds->UncheckedAt(0)); - AliTPCseed * vseed = new AliTPCseed; - vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1); - vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1); - vseed->fPoints->ExpandCreateFast(2); + Int_t gap = 15; + Float_t cuts[4]; + Float_t fnumber = 3.0; + Float_t fdensity = 3.0; - //TBranch * seedbranch = -// seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99); - //delete vseed; - nseed=fSeeds->GetEntriesFast(); + // 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 - Int_t found = 0; - for (i=0; iUncheckedAt(i), &t=*pt; - if (!pt) continue; - Int_t nc=t.GetNumberOfClusters(); - if (nc<20) continue; - t.CookdEdx(0.02,0.6); - CookLabel(pt,0.1); //For comparison only - // if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){ - if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (t.fNFoundable>Int_t(0.3*nrows)))){ - iotrack=pt; - tracktree.Fill(); -// cerr<RebuildSeed(); - seedbranch->SetAddress(&pt); + 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; + // - seedtree.Fill(); - for (Int_t j=0;j<160;j++){ - delete pt->fPoints->RemoveAt(j); +} + + +void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2) const +{ + // + //sum tracks to common container + //remove suspicious tracks + Int_t nseed = arr2->GetEntriesFast(); + 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; } - delete pt->fPoints; - pt->fPoints =0; - */ - delete fSeeds->RemoveAt(i); + // 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)); + continue; + } + //remove not usable tracks + if (pt->GetRemoval()!=10){ + 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)); + else{ + delete arr2->RemoveAt(i); + } + } + } + delete arr2; +} + + + +void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast) +{ + // + // try to track in parralel + + Int_t nseed=arr->GetEntriesFast(); + //prepare seeds for tracking + for (Int_t i=0; iUncheckedAt(i), &t=*pt; + if (!pt) continue; + if (!t.IsActive()) continue; + // follow prolongation to the first layer + if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) ) + FollowProlongation(t, rfirst+1); + } + + + // + for (Int_t nr=rfirst; nr>=rlast; nr--){ + if (nrGetNRows()) + fSectors = fInnerSec; + else + fSectors = fOuterSec; + // make indexes with the cluster tracks for given + + // find nearest cluster + 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())GetRelativeSector()>17) { + continue; + } + UpdateClusters(t,nr); + } + // prolonagate to the nearest cluster - if founded + for (Int_t i=0; iUncheckedAt(i); + if (!pt) continue; + if (!pt->IsActive()) continue; + // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())GetRelativeSector()>17) { + continue; + } + FollowToNextCluster(*pt,nr); + } + } +} + +void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const +{ + // + // + // if we use TPC track itself we have to "update" covariance + // + Int_t nseed= arr->GetEntriesFast(); + for (Int_t i=0;iUncheckedAt(i); + if (pt) { + pt->Modify(fac); + // + //rotate to current local system at first accepted point + Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint()); + Int_t sec = (index&0xff000000)>>24; + sec = sec%18; + Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift(); + if (angle1>TMath::Pi()) + angle1-=2.*TMath::Pi(); + Float_t angle2 = pt->GetAlpha(); + + if (TMath::Abs(angle1-angle2)>0.001){ + pt->Rotate(angle1-angle2); + //angle2 = pt->GetAlpha(); + //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha(); + //if (pt->GetAlpha()<0) + // pt->fRelativeSector+=18; + //sec = pt->fRelativeSector; + } + + } + + } + + +} +void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const +{ + // + // + // if we use TPC track itself we have to "update" covariance + // + Int_t nseed= arr->GetEntriesFast(); + for (Int_t i=0;iUncheckedAt(i); + if (pt) { + pt->Modify(fac); + pt->SetFirstPoint(pt->GetLastPoint()); + } + } - // fNTracks = found; - printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start(); - UnloadClusters(); - printf("Time for unloading cluster: \t"); timer.Print();timer.Start(); -// seedtree.Write(); - cerr<<"Number of found tracks : "<<"\t"<WriteTracks("OVERWRITE"); +} + +Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) +{ + // + // make back propagation + // + Int_t nseed= arr->GetEntriesFast(); + for (Int_t i=0;iUncheckedAt(i); + 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 (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->SetFirstPoint(kink->GetTPCRow0()); + fSectors = fInnerSec; + FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); + } + + } + return 0; +} + +Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr) +{ + // + // make forward propagation + // + Int_t nseed= arr->GetEntriesFast(); + // + for (Int_t i=0;iUncheckedAt(i); + if (pt) { + FollowProlongation(*pt,0); + } + } return 0; } -void AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast) +Int_t AliTPCtrackerMI::PropagateForward() { // - // try to track in parralel - - Int_t nseed=fSeeds->GetEntriesFast(); - //prepare seeds for tracking - for (Int_t i=0; iUncheckedAt(i), &t=*pt; - if (!pt) continue; - if (!t.IsActive()) continue; - // follow prolongation to the first layer - if ( (fSectors ==fInnerSec) || (t.fFirstPoint>rfirst+1)) - FollowProlongation(t, rfirst+1); + // propagate track forward + //UnsignClusters(); + Int_t nseed = fSeeds->GetEntriesFast(); + for (Int_t i=0;iUncheckedAt(i); + if (pt){ + AliTPCseed &t = *pt; + Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift(); + if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); + if (alpha < 0. ) alpha += 2.*TMath::Pi(); + t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN); + } } + + fSectors = fOuterSec; + ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows()); + fSectors = fInnerSec; + ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0); + //WriteTracks(); + return 1; +} - // - for (Int_t nr=rfirst; nr>=rlast; nr--){ - // make indexes with the cluster tracks for given - // for (Int_t i = 0;iUncheckedAt(i), &t=*pt; - if (!pt) continue; - if (!pt->IsActive()) continue; - if ( (fSectors ==fOuterSec) && pt->fFirstPointfRelativeSector>17) { - continue; - } - UpdateClusters(t,i,nr); - } - // prolonagate to the nearest cluster - if founded - for (Int_t i=0; iUncheckedAt(i); - if (!pt) continue; - if (!pt->IsActive()) continue; - if ((fSectors ==fOuterSec) &&pt->fFirstPointfRelativeSector>17) { - continue; - } - FollowToNextCluster(i,nr); - } - // for (Int_t i= 0;iGetNRows()) + r1 = row1; + else + r1 = fSectors->GetNRows()-1; + + if (row0GetNRows()&& r1>0 ) + FollowBackProlongation(*pt,r1); + if (row1<=fSectors->GetNRows()) + return 0; + // + r1 = row1 - fSectors->GetNRows(); + if (r1<=0) return 0; + if (r1>=fOuterSec->GetNRows()) return 0; + fSectors = fOuterSec; + return FollowBackProlongation(*pt,r1); + } + return 0; +} + + + + +void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row) +{ + // + // + Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); + // Float_t padlength = fParam->GetPadPitchLength(seed->fSector); + Float_t padlength = GetPadPitchLength(row); + // + Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3; + Double_t angulary = seed->GetSnp(); + angulary = angulary*angulary/(1-angulary*angulary); + seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy); + // + Float_t sresz = fParam->GetZSigma(); + Float_t angularz = seed->GetTgl(); + seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz); + /* + Float_t wy = GetSigmaY(seed); + Float_t wz = GetSigmaZ(seed); + wy*=wy; + wz*=wz; + 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 padlength = fParam->GetPadPitchLength(seed->GetSector()); + Float_t sres = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3; Float_t angular = seed->GetSnp(); angular = angular*angular/(1-angular*angular); // angular*=angular; @@ -1928,7 +6237,7 @@ 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 padlength = fParam->GetPadPitchLength(seed->GetSector()); Float_t sres = fParam->GetZSigma(); Float_t angular = seed->GetTgl(); Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres); @@ -1936,42 +6245,131 @@ Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed) } - -//_________________________________________________________________________ -AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const { +//__________________________________________________________________________ +void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const { //-------------------------------------------------------------------- - // Return pointer to a given cluster + //This function "cooks" a track label. If label<0, this track is fake. //-------------------------------------------------------------------- - Int_t sec=(index&0xff000000)>>24; - Int_t row=(index&0x00ff0000)>>16; - Int_t ncl=(index&0x0000ffff)>>00; + 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->GetClusterPointer(i)){ + clusters[current]=t->GetClusterPointer(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++; + } + + 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; + } - AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row); - if (!clrow) return 0; - return (AliTPCclusterMI*)(*clrow)[ncl]; + t->SetLabel(lab); + + // delete[] lb; + //delete[] mx; + //delete[] clusters; } + //__________________________________________________________________________ -void AliTPCtrackerMI::CookLabel(AliKalmanTrack *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. //-------------------------------------------------------------------- Int_t noc=t->GetNumberOfClusters(); - Int_t *lb=new Int_t[noc]; - Int_t *mx=new Int_t[noc]; - AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc]; + if (noc<10){ + //printf("\nnot founded prolongation\n\n\n"); + //t->Dump(); + return -1; + } + 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; - for (i=0; iGetClusterIndex(i); - clusters[i]=GetClusterMI(index); + 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; + // + //clusters[current]=GetClusterMI(index); + if (t->GetClusterPointer(i)){ + clusters[current]=t->GetClusterPointer(i); + current++; + } } - + noc = current; + if (noc<5) return -1; Int_t lab=123456789; for (i=0; iGetLabel(0)); Int_t j; for (j=0; jGetLabel(1)) == lab || TMath::Abs(c->GetLabel(2)) == lab ) max++; } @@ -1994,21 +6392,60 @@ void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const { else { Int_t tail=Int_t(0.10*noc); max=0; - for (i=1; i<=tail; i++) { - AliTPCclusterMI *c=clusters[noc-i]; - if (!clusters[i]) continue; + 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); + // t->SetLabel(lab); + return lab; + // delete[] lb; + //delete[] mx; + //delete[] clusters; +} + + +Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const +{ + //return pad row number for this x + Double_t r; + if (fN < 64){ + r=fRow[fN-1].GetX(); + if (x > r) return fN; + r=fRow[0].GetX(); + if (x < r) return -1; + return Int_t((x-r)/fPadPitchLength + 0.5);} + else{ + r=fRow[fN-1].GetX(); + if (x > r) return fN; + r=fRow[0].GetX(); + if (x < r) return -1; + Double_t r1=fRow[64].GetX(); + if(xGetPadRowRadiiLow(i)); - fRow[i].fDeadZone =1.5; //1.5 cm of dead zone + fRow[i].SetDeadZone(1.5); //1.5 cm of dead zone } } else { fAlpha=par->GetOuterAngle(); @@ -2039,70 +6476,32 @@ void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) { fRow=new AliTPCRow[fN]; for (Int_t i=0; iGetPadRowRadiiUp(i)); - fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone + fRow[i].SetDeadZone(1.5); // 1.5 cm of dead zone } } } - -AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){ +AliTPCtrackerMI::AliTPCRow::AliTPCRow(): + fDeadZone(0.), + fClusters1(0), + fN1(0), + fClusters2(0), + fN2(0), + fN(0), + fX(0.) +{ + // + // default constructor // - if (fClusterTracks) delete [] fClusterTracks; - fClusterTracks = 0; } -void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){ - //create cluster tracks - if (fN>0) - fClusterTracks = new AliTPCclusterTracks[fN]; -} +AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){ + // -void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){ - if (fClusterTracks) delete[] fClusterTracks; - fClusterTracks =0; } -void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){ - // - // - // update information of the cluster tracks - if track is nearer then other tracks to the - // given track - const AliTPCclusterMI * cl = (*this)[clindex]; - AliTPCclusterTracks * cltracks = GetClusterTracks(clindex); - // find the distance of the cluster to the track - Float_t dy2 = (cl->GetY()- seed->GetY()); - dy2*=dy2; - Float_t dz2 = (cl->GetZ()- seed->GetZ()); - dz2*=dz2; - // - Float_t distance = TMath::Sqrt(dy2+dz2); - if (distance> 3) - return; // MI - to be changed - AliTPCtrackerParam - - if ( distance < cltracks->fDistance[0]){ - cltracks->fDistance[2] =cltracks->fDistance[1]; - cltracks->fDistance[1] =cltracks->fDistance[0]; - cltracks->fDistance[0] =distance; - cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1]; - cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0]; - cltracks->fTrackIndex[0] =trindex; - } - else - if ( distance < cltracks->fDistance[1]){ - cltracks->fDistance[2] =cltracks->fDistance[1]; - cltracks->fDistance[1] =distance; - cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1]; - cltracks->fTrackIndex[1] =trindex; - } else - if (distance < cltracks->fDistance[2]){ - cltracks->fDistance[2] =distance; - cltracks->fTrackIndex[2] =trindex; - } -} - - //_________________________________________________________________________ void AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) { @@ -2119,6 +6518,21 @@ AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index fIndex[i]=index; fClusters[i]=c; fN++; } +void AliTPCtrackerMI::AliTPCRow::ResetClusters() { + // + // reset clusters + fN = 0; + fN1 = 0; + fN2 = 0; + //delete[] fClusterArray; + if (fClusters1) delete []fClusters1; + if (fClusters2) delete []fClusters2; + //fClusterArray=0; + fClusters1 = 0; + fClusters2 = 0; +} + + //___________________________________________________________________ Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const { //----------------------------------------------------------------------- @@ -2158,299 +6572,76 @@ AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z return cl; } - - - - -AliTPCseed::AliTPCseed():AliTPCtrack(){ - // - fRow=0; - fRemoval =0; - memset(fClusterIndex,0,sizeof(Int_t)*200); - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared =0; - fTrackPoints =0; - fRemoval = 0; -} - -AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){ - fPoints = 0; - fEPoints = 0; - fNShared =0; - fTrackPoints =0; - fRemoval =0; -} - -AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){ - fRow=0; - memset(fClusterIndex,0,sizeof(Int_t)*200); - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared =0; - fTrackPoints =0; - fRemoval =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) { - // - // - fRow =0; - memset(fClusterIndex,0,sizeof(Int_t)*200); - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared = 0; - fTrackPoints =0; - fRemoval =0; -} - -AliTPCseed::~AliTPCseed(){ - if (fPoints) delete fPoints; - fPoints =0; - fEPoints = 0; - if (fTrackPoints){ - for (Int_t i=0;i<8;i++){ - delete [] fTrackPoints[i]; - } - delete fTrackPoints; - fTrackPoints =0; - } - -} - -AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i) +AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const { - // - // - if (!fTrackPoints) { - fTrackPoints = new AliTPCTrackPoint*[8]; - for ( Int_t i=0;i<8;i++) - fTrackPoints[i]=0; - } - Int_t index1 = i/20; - if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20]; - return &(fTrackPoints[index1][i%20]); -} + //----------------------------------------------------------------------- + // Return the index of the nearest cluster in z y + //----------------------------------------------------------------------- + Float_t maxdistance = roady*roady + roadz*roadz; + AliTPCclusterMI *cl =0; -void AliTPCseed::RebuildSeed() -{ - // - // rebuild seed to be ready for storing - fPoints = new TClonesArray("AliTPCTrackPoint",160); - fPoints->ExpandCreateFast(160); - fEPoints = new TClonesArray("AliTPCExactPoint",1); - for (Int_t i=0;i<160;i++){ - AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i); - *trpoint = *(GetTrackPoint(i)); + //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(GetFastCluster(iz1)-1,0); + Int_t iz2 = Int_t(z+roadz+255.5); + if (iz2<0 || iz2>=510) return cl; + iz2 = TMath::Min(GetFastCluster(iz2)+1,fN); + + //FindNearest3(y,z,roady,roadz,index); + // for (Int_t i=Find(z-roadz); iGetZ() > z+roadz) break; + if ( c->GetY()-y > roady ) continue; + if ( y-c->GetY() > roady ) continue; + Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y); + if (maxdistance>distance) { + maxdistance = distance; + cl=c; + index =i; + //roady = TMath::Sqrt(maxdistance); + } } - + return cl; } -//_____________________________________________________________________________ -void AliTPCseed::CookdEdx(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]; - //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;i<160;i+=4) - { - //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i); - AliTPCTrackPoint * point = GetTrackPoint(i); - if (point==0) continue; - if (point->fIsShared){ - fNShared++; - continue; - } - if (point->GetCPoint().GetMax()<5) continue; - Float_t angley = point->GetTPoint().GetAngleY(); - Float_t anglez = point->GetTPoint().GetAngleZ(); - Int_t type = point->GetCPoint().GetType(); - 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 = 1.*point->GetCPoint().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*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(point->GetCPoint().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]++; +AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const +{ + //----------------------------------------------------------------------- + // Return the index of the nearest cluster in z y + //----------------------------------------------------------------------- + Float_t maxdistance = roady*roady + roadz*roadz; + // Int_t iz = Int_t(z+255.); + AliTPCclusterMI *cl =0; + for (Int_t i=Find(z-roadz); iGetZ() > z+roadz) break; + if ( c->GetY()-y > roady ) continue; + if ( y-c->GetY() > roady ) continue; + Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y); + if (maxdistance>distance) { + maxdistance = distance; + cl=c; + index =i; + //roady = TMath::Sqrt(maxdistance); } - - 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 = 200; - 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 - - + return cl; +} - //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; -} +// AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i) +// { +// // +// // +// return &fTrackPoints[i]; +// }