X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=ITS%2FAliITStrackerSA.cxx;h=ea8abebb330c94bfd3515b6d7877e8c36e4c35eb;hp=f04b682d5619cdf0659d51a984681b4e71d817b6;hb=f1d424e80153f84428dd5c169fc57ad7437d3d1c;hpb=c152044c832d37acc5f95d057932613f1decbd77 diff --git a/ITS/AliITStrackerSA.cxx b/ITS/AliITStrackerSA.cxx index f04b682d561..ea8abebb330 100644 --- a/ITS/AliITStrackerSA.cxx +++ b/ITS/AliITStrackerSA.cxx @@ -1,4 +1,3 @@ - /************************************************************************** * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. * * * @@ -14,134 +13,146 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -//////////////////////////////////////////////////// -// Stand alone tracker class // -// Origin: Elisabetta Crescio // -// e-mail: crescio@to.infn.it // -// tracks are saved as AliITStrackV2 objects // -//////////////////////////////////////////////////// +/* $Id$ */ + +/////////////////////////////////////////////////////////// +// Stand alone ITS tracker class // +// Origin: Elisabetta Crescio - crescio@to.infn.it // +// Updated: Francesco Prino - prino@to.infn.it // +/////////////////////////////////////////////////////////// #include -#include "TArrayI.h" + +#include #include -#include #include #include -#include "AliRun.h" -#include "AliITSclusterTable.h" -#include "AliITSclusterV2.h" -#include "AliITSgeom.h" -#include "AliITSRiemannFit.h" -#include "AliITStrackerSA.h" -#include "AliITStrackSA.h" -#include "AliITSVertexer.h" + +#include "AliESDEvent.h" #include "AliESDVertex.h" -#include "AliESD.h" #include "AliESDtrack.h" +#include "AliITSVertexer.h" +#include "AliITSclusterTable.h" +#include "AliITSRecPoint.h" +#include "AliITSgeomTGeo.h" +#include "AliITStrackSA.h" +#include "AliITStrackerSA.h" +#include "AliITSReconstructor.h" +#include "AliLog.h" +#include "AliRun.h" ClassImp(AliITStrackerSA) //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA():AliITStrackerV2(){ +AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(), +fPhiEstimate(0), +fITSStandAlone(0), +fLambdac(0), +fPhic(0), +fCoef1(0), +fCoef2(0), +fCoef3(0), +fNloop(0), +fPhiWin(0), +fLambdaWin(0), +fListOfTracks(0), +fListOfSATracks(0), +fITSclusters(0), +fInwardFlag(0), +fOuterStartLayer(0), +fInnerStartLayer(5), +fMinNPoints(0), +fMinQ(0.), +fCluCoord(0){ // Default constructor Init(); + } //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom):AliITStrackerV2(geom) +AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0), +fPhiEstimate(0), +fITSStandAlone(0), +fLambdac(0), +fPhic(0), +fCoef1(0), +fCoef2(0), +fCoef3(0), +fNloop(0), +fPhiWin(0), +fLambdaWin(0), +fListOfTracks(0), +fListOfSATracks(0), +fITSclusters(0), +fInwardFlag(0), +fOuterStartLayer(0), +fInnerStartLayer(5), +fMinNPoints(0), +fMinQ(0.), +fCluCoord(0) { // Standard constructor (Vertex is known and passed to this obj.) - Init(); - fVert = 0; - fGeom = geom; -} + if (geom) { + AliWarning("\"geom\" is actually a dummy argument !"); + } -//____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliESDVertex *vert):AliITStrackerV2(geom) -{ - // Standard constructor (Vertex is known and passed to this obj.) Init(); - fVert = vert; - fGeom = geom; -} - -//______________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA &trkr) : - AliITStrackerV2(trkr) { - // Copy constructor - // Copies are not allowed. The method is protected to avoid misuse. - Error("AliITStrackerSA","Copy constructor not allowed\n"); -} - -//______________________________________________________________________ -AliITStrackerSA& AliITStrackerSA::operator=(const - AliITStrackerSA& /* trkr */){ - // Assignment operator - // Assignment is not allowed. The method is protected to avoid misuse. - Error("= operator","Assignment operator not allowed\n"); - return *this; + } //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliITSVertexer *vertexer):AliITStrackerV2(geom) -{ - // Standard constructor (Vertex is unknown - vertexer is passed to this obj) - Init(); - fVertexer = vertexer; - fGeom = geom; +AliITStrackerSA::~AliITStrackerSA(){ + // destructor + + if(fPhiWin)delete []fPhiWin; + if(fLambdaWin)delete []fLambdaWin; + fListOfTracks->Delete(); + delete fListOfTracks; + fListOfSATracks->Delete(); + delete fListOfSATracks; + if(fCluCoord){ + for(Int_t i=0;iDelete(); + delete fCluCoord[i]; + } + } + delete [] fCluCoord; + } } //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITStrackerSA& tracker):AliITStrackerV2(){ - // Copy constructor - fPhiEstimate = tracker.fPhiEstimate; - for(Int_t i=0;i<2;i++){ - fPoint1[i]=tracker.fPoint1[i]; - fPoint2[i]=tracker.fPoint2[i]; - fPoint3[i]=tracker.fPoint3[i]; - fPointc[i]=tracker.fPointc[i]; - } - fLambdac = tracker.fLambdac; - fPhic = tracker.fPhic; - fCoef1 = tracker.fCoef1; - fCoef2 = tracker.fCoef2; - fCoef3 = tracker.fCoef3; - fNloop = tracker.fNloop; - fPhiWin = tracker.fPhiWin; - fLambdaWin = tracker.fLambdaWin; - if(tracker.fVertexer && tracker.fVert){ - fVert = new AliESDVertex(*tracker.fVert); +Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){ +// This method is used to find and fit the tracks. By default the corresponding +// method in the parent class is invoked. In this way a combined tracking +// TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking +// is done in the ITS only. In the standard reconstruction chain this option +// can be set via AliReconstruction::SetOption("ITS","onlyITS") + Int_t rc=0; + + if(!fITSStandAlone){ + rc=AliITStrackerMI::Clusters2Tracks(event); } else { - fVert = tracker.fVert; + AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n"); } - fVertexer = tracker.fVertexer; - fGeom = tracker.fGeom; - fTable = tracker.fTable; - fListOfTracks = tracker.fListOfTracks; -} - -//____________________________________________________________________________ -AliITStrackerSA::~AliITStrackerSA(){ - // destructor - // if fVertexer is not null, the AliESDVertex obj. is owned by this class - // and is deleted here - if(fVertexer){ - if(fVert)delete fVert; + if(!rc){ + rc=FindTracks(event,kFALSE); + Int_t nSPDcontr=0; + const AliESDVertex *spdv = event->GetPrimaryVertexSPD(); + if(spdv) nSPDcontr = spdv->GetNContributors(); + if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE && + nSPDcontr<=AliITSReconstructor::GetRecoParam()->GetMaxSPDcontrForSAToUseAllClusters()) { + rc=FindTracks(event,kTRUE); + } } - fVert = 0; - fVertexer = 0; - fGeom = 0; - if(fPhiWin)delete []fPhiWin; - if(fLambdaWin)delete []fLambdaWin; - fTable =0; - fListOfTracks->Delete(); - } + return rc; +} //____________________________________________________________________________ void AliITStrackerSA::Init(){ // Reset all data members fPhiEstimate=0; - for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;} + for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;} fLambdac=0; fPhic=0; fCoef1=0; @@ -149,20 +160,30 @@ void AliITStrackerSA::Init(){ fCoef3=0; fPointc[0]=0; fPointc[1]=0; - fVert = 0; - fVertexer = 0; - fGeom = 0; - SetWindowSizes(); - fTable = 0; + Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA(); + if(nLoops==32){ + SetFixedWindowSizes(); + }else{ + Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA(); + Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA(); + Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA(); + Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA(); + SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax); + } + fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge(); fITSclusters = 0; - SetSixPoints(); - fListOfTracks=new TObjArray(0,0); -} + SetOuterStartLayer(1); + SetSAFlag(kFALSE); + fListOfTracks=new TClonesArray("AliITStrackMI",100); + fListOfSATracks=new TClonesArray("AliITStrackSA",100); + fCluCoord = 0; + fMinNPoints = 3; + } //_______________________________________________________________________ void AliITStrackerSA::ResetForFinding(){ // Reset data members used in all loops during track finding fPhiEstimate=0; - for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;} + for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;} fLambdac=0; fPhic=0; fCoef1=0; @@ -170,922 +191,578 @@ void AliITStrackerSA::ResetForFinding(){ fCoef3=0; fPointc[0]=0; fPointc[1]=0; - fListOfTracks->Delete(); + fListOfTracks->Clear(); + fListOfSATracks->Clear(); } -//____________________________________________________________________________ -void AliITStrackerSA::FindTracks(TTree *out,Int_t evnumber){ - - /************************************************************************** - * This function finds primary tracks. - * * - * * - * Example: to execute function with only the ITS (no combined tracking * - * with TPC+ITS) and requiring 5/6 points to define a good track * - * call SetSixPoinbts(kFALSE) in advance and then * - * use: FindTracks(treein,treeout,evnumber) * - * to execute combined tracking, before using FindTracks, use * - * UseFoundTracksV2 * - *************************************************************************/ - if(!fITSclusters){ - Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n"); - exit(1); - } - //Get primary vertex - if(fVertexer){ - if(fVert)delete fVert; - fVert = fVertexer->FindVertexForCurrentEvent(evnumber); - } - else { - gAlice->GetEvent(evnumber); - if(!fVert){ - Fatal("FindTracks","Vertex is missing\n"); - return; - } - } - Double_t primaryVertex[3]; - Double_t errorsprimvert[3]; - fVert->GetXYZ(primaryVertex); - fVert->GetSigmaXYZ(errorsprimvert); - if(errorsprimvert[0]==0 || errorsprimvert[1]==0){ - Warning("FindTracks","Set errors on vertex positions x and y at 0.0001"); - errorsprimvert[0]=0.0001; - errorsprimvert[1]=0.0001; - } - fVert->PrintStatus(); + +//______________________________________________________________________ +Int_t AliITStrackerSA::FindTracks(AliESDEvent* event, Bool_t useAllClusters){ - //Fill array with cluster indices for each module - if(!fTable){ - fTable = new AliITSclusterTable(fGeom,this,primaryVertex); - fTable->FillArray(fITSclusters); - fTable->FillArrayCoorAngles(); - } +// Track finder using the ESD object - -//Fill tree for found tracks - AliITStrackV2* outrack=0; - TBranch* branch=out->Branch("tracks","AliITStrackV2",&outrack,32000,0); - if (!branch) out->Branch("tracks","AliITStrackV2",&outrack,32000,3); - else branch->SetAddress(&outrack); + AliDebug(2,Form(" field is %f",event->GetMagneticField())); + AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5))); - - Int_t * firstmod = new Int_t[fGeom->GetNlayers()]; - for(Int_t i=0;iGetNlayers();i++){ - firstmod[i]=fGeom->GetModuleIndex(i+1,1,1); + if(!fITSclusters){ + Fatal("FindTracks","ITS cluster tree is not accessed!!!\n Please use method SetClusterTree to pass the pointer to the tree\n"); + return -1; } - // firstmod [i] number of the first module in the ITS layer i. - - AliITSlayer &layer=fgLayers[0]; // first layer - Int_t ntrack=0; - Int_t dim=layer.GetNumberOfClusters(); - //loop on the different windows - for(Int_t nloop=0;nloopIsUsed()==1) continue; - fPhic = fTable->GetPhiCluster(0,ncl); - fLambdac = fTable->GetLambdaCluster(0,ncl); - fPhiEstimate = fPhic; - AliITStrackSA* trs = new AliITStrackSA(); - fPoint1[0]=primaryVertex[0]; - fPoint1[1]=primaryVertex[1]; - fPoint2[0]=fTable->GetXCluster(0,ncl); - fPoint2[1]=fTable->GetYCluster(0,ncl); - - Int_t * nn = new Int_t[fGeom->GetNlayers()];//counter for clusters on each layer - for(Int_t i=0;iGetNlayers();i++){ nn[i]=0;} - nn[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - nn[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - - if(nn[1]>0){ - pflag=1; - fPoint3[0] = fPointc[0]; - fPoint3[1] = fPointc[1]; - } - nn[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - if(nn[1]==0 && nn[2]==0) pflag=0; - if(nn[2]!=0 && nn[1]!=0){ pflag=1; UpdatePoints();} - if(nn[2]!=0 && nn[1]==0){ - pflag=1; - fPoint3[0]=fPointc[0]; - fPoint3[1]=fPointc[1]; - } - - nn[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - pflag=1; - if(nn[3]!=0) UpdatePoints(); - nn[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - pflag=1; - if(nn[4]!=0) UpdatePoints(); - nn[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - - - Int_t layOK=0; - Int_t numberofpoints; - if(fSixPoints) numberofpoints=6; //check of the candidate track - else numberofpoints=5; //if track is good (with the required number - for(Int_t nnp=0;nnpGetNlayers();nnp++){ //of points) it is written on file - if(nn[nnp]!=0) layOK+=1; + //Reads event and mark clusters of traks already found, with flag kITSin + Int_t nentr=event->GetNumberOfTracks(); + if(!useAllClusters) { + while (nentr--) { + AliESDtrack *track=event->GetTrack(nentr); + if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){ + Int_t idx[12]; + Int_t ncl = track->GetITSclusters(idx); + for(Int_t k=0;kSetBit(kSAflag); + } } - if(layOK>=numberofpoints){ - AliITStrackV2* tr2 = FitTrack(trs,primaryVertex,errorsprimvert); - if(tr2==0){ - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==1) kl->Use(); - } - continue; - } - outrack=tr2; - out->Fill(); - ntrack++; - Int_t nct = tr2->GetNumberOfClusters(); - - while(nct--){ - Int_t index = tr2->GetClusterIndex(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==0) kl->Use(); - - } - } - else{ - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==1) kl->Use(); - } + } + }else{ + while (nentr--) { + AliESDtrack *track=event->GetTrack(nentr); + if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){ + Int_t idx[12]; + Int_t ncl = track->GetITSclusters(idx); + for(Int_t k=0;kResetBit(kSAflag); + } } - delete trs; - - }//end loop on clusters of layer1 - - }//end loop2 - - //if 5/6 points are required, second loop starting - //from second layer, to find tracks with point of - //layer 1 missing - - if(!fSixPoints){ - // counter for clusters on each layer - Int_t * nn = new Int_t[fGeom->GetNlayers()-1]; - for(Int_t nloop=0;nloopIsUsed()==1) continue; - fPhic = fTable->GetPhiCluster(1,ncl2); - fLambdac = fTable->GetLambdaCluster(1,ncl2); - fPhiEstimate = fPhic; - AliITStrackSA* trs = new AliITStrackSA(); - - fPoint1[0]=primaryVertex[0]; - fPoint1[1]=primaryVertex[1]; - fPoint2[0]=fTable->GetXCluster(1,ncl2);; - fPoint2[1]=fTable->GetYCluster(1,ncl2);; - - for(Int_t kk=0;kkGetNlayers()-1;kk++)nn[kk] = 0; - for(Int_t kk=0;kkGetNlayers()-1;kk++){ - nn[kk] = SearchClusters(kk+1,fPhiWin[nloop],fLambdaWin[nloop], - trs,primaryVertex[2],pflag,fTable); - if(nn[kk]==0)break; - if(kk>0){ - UpdatePoints(); - pflag = 1; - } - } - Int_t fl=0; - for(Int_t nnp=0;nnpGetNlayers()-1;nnp++){ - if(nn[nnp]!=0) fl+=1; - } - if(fl>=5){ // 5/6 - AliITStrackV2* tr2 = FitTrack(trs,primaryVertex,errorsprimvert); - if(tr2==0){ - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==1) kl->Use(); - } - continue; - } - outrack=tr2; - out->Fill(); - Int_t nct = tr2->GetNumberOfClusters(); - while(nct--){ - Int_t index = tr2->GetClusterIndex(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl==0) continue; - if(kl->IsUsed()==0) kl->Use(); - } - } - else{ - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl==0) continue; - if(kl->IsUsed()==1) kl->Use(); - } - } - delete trs; - }//end loop on clusters of layer2 } - delete [] nn; - } // if(!fSixPoints.... - - delete [] firstmod; - delete fTable; - fTable = 0; -} - - -//______________________________________________________________________ -Int_t AliITStrackerSA::FindTracks(AliESD* event){ - - // Track finder using the ESD object - - - if(!fITSclusters){ - Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n"); - return -1; } - //Get primary vertex - Double_t errorsprimvert[3]; Double_t primaryVertex[3]; event->GetVertex()->GetXYZ(primaryVertex); - event->GetVertex()->GetSigmaXYZ(errorsprimvert); - - if(errorsprimvert[0]==0 || errorsprimvert[1]==0){ - // Warning("FindTracks","Set errors on vertex positions x and y at 0.005"); - errorsprimvert[0]=0.005; - errorsprimvert[1]=0.005; + //Creates TClonesArray with clusters for each layer. The clusters already used + //by AliITStrackerMI are not considered + Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0}; + Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0}; + if (fCluCoord == 0) { + fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers]; + for(Int_t i=0;iTestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol. + if(cls->GetQ()==0) continue; //fake clusters dead zones + if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge + nclusters[i]++; + } + } + dmar[i]=0; + if(!fCluCoord[i]){ + fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]); + }else{ + fCluCoord[i]->Delete(); + fCluCoord[i]->Expand(nclusters[i]); + } } - //Fill array with cluster indices for each module - if(!fTable){ - fTable = new AliITSclusterTable(fGeom,this,primaryVertex); - fTable->FillArray(fITSclusters); - fTable->FillArrayCoorAngles(); - } - - Int_t * firstmod = new Int_t[fGeom->GetNlayers()]; - for(Int_t i=0;iGetNlayers();i++){ - firstmod[i]=fGeom->GetModuleIndex(i+1,1,1); + for(Int_t ilay=0;ilayTestBit(kSAflag)==kTRUE) continue; + if(cls->GetQ()==0) continue; + if(ilay>1 && cls->GetQ()<=fMinQ) continue; + Double_t phi=0;Double_t lambda=0; + Double_t x=0;Double_t y=0;Double_t z=0; + Double_t sx=0;Double_t sy=0;Double_t sz=0; + GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex); + GetCoorErrors(cls,sx,sy,sz); + new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli); + dmar[ilay]++; + } + } + fCluCoord[ilay]->Sort(); } - // firstmod [i] number of the first module in the ITS layer i. - - - AliITSlayer &layer=fgLayers[0]; // first layer + + // track counter Int_t ntrack=0; - Int_t dim=layer.GetNumberOfClusters(); - //loop on the different windows - for(Int_t nloop=0;nloopIsUsed()==1) continue; - fPhic = fTable->GetPhiCluster(0,ncl); - fLambdac = fTable->GetLambdaCluster(0,ncl); - fPhiEstimate = fPhic; - AliITStrackSA* trs = new AliITStrackSA(); - fPoint1[0]=primaryVertex[0]; - fPoint1[1]=primaryVertex[1]; - fPoint2[0]=fTable->GetXCluster(0,ncl); - fPoint2[1]=fTable->GetYCluster(0,ncl); - Int_t * nn = new Int_t[fGeom->GetNlayers()];//counter for clusters on each layer - for(Int_t i=0;iGetNlayers();i++){ nn[i]=0;} - nn[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - - nn[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - if(nn[1]>0){ - pflag=1; - fPoint3[0] = fPointc[0]; - fPoint3[1] = fPointc[1]; - } - nn[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - if(nn[1]==0 && nn[2]==0) pflag=0; - if(nn[2]!=0 && nn[1]!=0){ pflag=1; UpdatePoints();} - if(nn[2]!=0 && nn[1]==0){ - pflag=1; - fPoint3[0]=fPointc[0]; - fPoint3[1]=fPointc[1]; - } + static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer + Int_t startLayForSeed=0; + Int_t lastLayForSeed=fOuterStartLayer; + Int_t nSeedSteps=lastLayForSeed-startLayForSeed; + Int_t seedStep=1; + if(fInwardFlag){ + startLayForSeed=AliITSgeomTGeo::GetNLayers()-1; + lastLayForSeed=fInnerStartLayer; + nSeedSteps=startLayForSeed-lastLayForSeed; + seedStep=-1; + } - nn[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - pflag=1; - if(nn[3]!=0) UpdatePoints(); - nn[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - pflag=1; - if(nn[4]!=0) UpdatePoints(); - nn[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag,fTable); - - - Int_t layOK=0; - Int_t numberofpoints; - if(fSixPoints) numberofpoints=6; //check of the candidate track - else numberofpoints=5; //if track is good (with the required number - for(Int_t nnp=0;nnpGetNlayers();nnp++){ //of points) it is written on file - if(nn[nnp]!=0) layOK+=1; - } - if(layOK>=numberofpoints){ - AliITStrackV2* tr2 = FitTrack(trs,primaryVertex,errorsprimvert); - if(tr2==0){ - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==1) kl->Use(); - } - continue; - } - - AliESDtrack outtrack; - outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin); - event->AddTrack(&outtrack); - ntrack++; - Int_t nct = tr2->GetNumberOfClusters(); - while(nct--){ - Int_t index = tr2->GetClusterIndex(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==0) kl->Use(); + // loop on minimum number of points + for(Int_t iMinNPoints=AliITSgeomTGeo::GetNLayers(); iMinNPoints>=fMinNPoints; iMinNPoints--) { + + // loop on starting layer for track finding + for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) { + Int_t theLay=startLayForSeed+iSeedLay*seedStep; + if(ForceSkippingOfLayer(theLay)) continue; + Int_t minNPoints=iMinNPoints-theLay; + if(fInwardFlag) minNPoints=iMinNPoints-(AliITSgeomTGeo::GetNLayers()-1-theLay); + for(Int_t i=theLay+1;iGetEntries(); + while(nclTheLay--){ + ResetForFinding(); + Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex); + if(!useRP) continue; + AliITStrackSA trs; + + Int_t pflag=0; + Int_t kk; + for(kk=0;kk=0; nnp--){ + if(nClusLay[nnp]!=0) layOK+=1; + } + } + if(layOK>=minNPoints){ + AliDebug(2,Form("---NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5])); + AliITStrackV2* tr2 = 0; + tr2 = FitTrack(&trs,primaryVertex); + if(!tr2){ + continue; + } + AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters())); + + StoreTrack(tr2,event,useAllClusters); + ntrack++; + + } - } - } - else{ - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==1) kl->Use(); - } - } - delete trs; - delete[] nn; - - }//end loop on clusters of layer1 - - }//end loop2 - - //if 5/6 points are required, second loop starting - //from second layer, to find tracks with point of - //layer 1 missing - - if(!fSixPoints){ - // counter for clusters on each layer - Int_t * nn = new Int_t[fGeom->GetNlayers()-1]; - for(Int_t nloop=0;nloopIsUsed()==1) continue; - fPhic = fTable->GetPhiCluster(1,ncl2); - fLambdac = fTable->GetLambdaCluster(1,ncl2); - fPhiEstimate = fPhic; - AliITStrackSA* trs = new AliITStrackSA(); - fPoint1[0]=primaryVertex[0]; - fPoint1[1]=primaryVertex[1]; - fPoint2[0]=fTable->GetXCluster(1,ncl2); - fPoint2[1]=fTable->GetYCluster(1,ncl2); - - for(Int_t kk=0;kkGetNlayers()-1;kk++)nn[kk] = 0; - for(Int_t kk=0;kkGetNlayers()-1;kk++){ - nn[kk] = SearchClusters(kk+1,fPhiWin[nloop],fLambdaWin[nloop], - trs,primaryVertex[2],pflag,fTable); - if(nn[kk]==0)break; - if(kk>0){ - UpdatePoints(); - pflag = 1; - } - } - Int_t fl=0; - for(Int_t nnp=0;nnpGetNlayers()-1;nnp++){ - if(nn[nnp]!=0) fl+=1; - } - if(fl>=5){ // 5/6 - AliITStrackV2* tr2 = FitTrack(trs,primaryVertex,errorsprimvert); - if(tr2==0){ + }//end loop on clusters of theLay + } //end loop on window sizes + } //end loop on theLay + }//end loop on min points + + // search for 1-point tracks in SPD, only for cosmics + // (A.Dainese 21.03.08) + if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() && + TMath::Abs(event->GetMagneticField())<0.01) { + Int_t outerLayer=1; // only SPD + for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) { + // counter for clusters on each layer + + for(Int_t nloop=0;nloopGetEntries(); + while(nclInnLay--){ //loop starting from layer innLay + ResetForFinding(); + Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex); + if(!useRP) continue; + AliITStrackSA trs; + + Int_t pflag=0; + Int_t kk; + for(kk=0;kkGetNumberOfClusters())); + + StoreTrack(tr2,event,useAllClusters); ntrack++; - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl->IsUsed()==1) kl->Use(); - } - - continue; - } - - AliESDtrack outtrack; - outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin); - event->AddTrack(&outtrack); - Int_t nct = tr2->GetNumberOfClusters(); - while(nct--){ - Int_t index = tr2->GetClusterIndex(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl==0) continue; - if(kl->IsUsed()==0) kl->Use(); - } - } - else{ - Int_t nct = trs->GetNumberOfClustersSA(); - while(nct--){ - Int_t index = trs->GetClusterIndexSA(nct); - AliITSclusterV2* kl = (AliITSclusterV2*)GetCluster(index); - if(kl==0) continue; - if(kl->IsUsed()==1) kl->Use(); - } - } - delete trs; - }//end loop on clusters of layer2 - } - delete [] nn; - } //end opt="5/6" - - delete [] firstmod; - delete fTable; - fTable = 0; - Info("FindTracks","Number of found tracks: %d",event->GetNumberOfTracks()); - + + } + + }//end loop on clusters of innLay + } //end loop on window sizes + + } //end loop on innLay + } // end search 1-point tracks + + if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks())); + ResetForFinding(); return 0; } - - + //________________________________________________________________________ -AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Double_t *errorsprimvert){ - //fit of the found track +AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) { + //fit of the found track (most general case, also <6 points, layers missing) + // A.Dainese 16.11.07 - Int_t * firstmod = new Int_t[fGeom->GetNlayers()]; - for(Int_t i=0;iGetNlayers();i++){ - firstmod[i]=fGeom->GetModuleIndex(i+1,1,1); - } - AliITStrackV2* otrack2; - Int_t nclusters = tr->GetNumberOfClustersSA(); - TObjArray** listlayer = new TObjArray*[fGeom->GetNlayers()]; - for(Int_t i=0;iGetNlayers();i++){ - listlayer[i] = new TObjArray(0,0); + const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters; + + static Int_t firstmod[AliITSgeomTGeo::kNLayers]; + static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu]; + static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu]; + static Int_t end[AliITSgeomTGeo::kNLayers]; + static Int_t indices[AliITSgeomTGeo::kNLayers]; + + static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu]; + + for(Int_t i=0;iGetNlayers()]; - for(Int_t i=0;iGetNlayers();i++)nnn[i]=0; - for(Int_t ncl=0;nclGetClusterIndexSA(ncl); - AliITSclusterV2* cl = (AliITSclusterV2*)GetCluster(index); - if(cl->IsUsed()==1) cl->Use(); + Int_t nclusters = tr->GetNumberOfClustersSA(); + for(Int_t ncl=0;nclGetClusterIndexSA(ncl); + AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index); Int_t lay = (index & 0xf0000000) >> 28; - if(lay==0) { listlayer[0]->AddLast(cl); clind0[nnn[0]]=index;nnn[0]++;} - if(lay==1) { listlayer[1]->AddLast(cl); clind1[nnn[1]]=index;nnn[1]++;} - if(lay==2) { listlayer[2]->AddLast(cl); clind2[nnn[2]]=index;nnn[2]++;} - if(lay==3) { listlayer[3]->AddLast(cl); clind3[nnn[3]]=index;nnn[3]++;} - if(lay==4) { listlayer[4]->AddLast(cl); clind4[nnn[4]]=index;nnn[4]++;} - if(lay==5) { listlayer[5]->AddLast(cl); clind5[nnn[5]]=index;nnn[5]++;} + Int_t nInLay=end[lay]; + listlayer[lay][nInLay]=cl; + clind[lay][nInLay]=index; + end[lay]++; } - delete [] nnn; - - Int_t * end = new Int_t[fGeom->GetNlayers()]; - for(Int_t i=0;iGetNlayers();i++){ - if(listlayer[i]->GetEntries()==0) end[i]=1; - else end[i]=listlayer[i]->GetEntries(); + for(Int_t nlay=0;nlayGetNumberOfMarked(nlay);ncl++){ + Int_t mark = tr->GetClusterMark(nlay,ncl); + clmark[nlay][ncl]=mark; + } } - for(Int_t l1=0;l1At(l1); - TVector3** recp = new TVector3*[3]; - TVector3** errs = new TVector3*[3]; - recp[0] = new TVector3(primaryVertex[0],primaryVertex[1],primaryVertex[2]); - errs[0] = new TVector3(errorsprimvert[0],errorsprimvert[1],errorsprimvert[2]); - Double_t x1,y1,z1,sx1,sy1,sz1; - Double_t x2,y2,z2,sx2,sy2,sz2; - AliITSclusterV2* p1=0; - AliITSclusterV2* p2=0; - Int_t index1=clind0[l1]; - Int_t index2=0; - for(Int_t l2=0;l2At(l2); - index2=clind1[l2]; - for(Int_t l3=0;l3At(l3); - - if(cl0==0 && cl1!=0) { - p2 = cl2;index1=clind2[l3]; - p1=cl1; - - } - if(cl0!=0 && cl1==0){ - p1=cl0; - p2=cl2;index2=clind2[l3]; - } - if(cl0!=0 && cl1!=0){ - p1=cl0; - p2=cl1; - } - Int_t lay1=(index1 & 0xf0000000) >> 28; - Int_t cln1=(index1 & 0x0fffffff) >> 00; - Int_t lay2=(index2 & 0xf0000000) >> 28; - Int_t cln2=(index2 & 0x0fffffff) >> 00; - x1 = fTable->GetXCluster(lay1,cln1); - x2 = fTable->GetXCluster(lay2,cln2); - y1 = fTable->GetYCluster(lay1,cln1); - y2 = fTable->GetYCluster(lay2,cln2); - z1 = fTable->GetZCluster(lay1,cln1); - z2 = fTable->GetZCluster(lay2,cln2); - sx1 = fTable->GetXClusterError(lay1,cln1); - sx2 = fTable->GetXClusterError(lay2,cln2); - sy1 = fTable->GetYClusterError(lay1,cln1); - sy2 = fTable->GetYClusterError(lay2,cln2); - sz1 = fTable->GetZClusterError(lay1,cln1); - sz2 = fTable->GetZClusterError(lay2,cln2); - Double_t phi1 = fTable->GetPhiCluster(lay1,cln1); - Int_t module1 = p1->GetDetectorIndex()+firstmod[0]; - recp[1] = new TVector3(x1,y1,z1); - errs[1] = new TVector3(sx1,sy1,sz1); - recp[2] = new TVector3(x2,y2,z2); - errs[2] = new TVector3(sx2,sy2,sz2); - - //fit on the Riemann sphere - Float_t seed1,seed2,seed3; - AliITSRiemannFit fit; - Int_t rf = fit.FitHelix(3,recp,errs,seed1,seed2,seed3); //this gives phi,tgl,curvature to start Kalman Filter - if(rf==0) continue; - Double_t phi=seed1; - Double_t tgl=seed2; - - if(phi1>0){ - if(seed1>-TMath::Pi() && seed1<-0.5*TMath::Pi()){ - phi=seed1+1.5*TMath::Pi(); - tgl=seed2; - } - if(seed1>-0.5*TMath::Pi() && seed1<0.5*TMath::Pi()){ - phi=seed1+0.5*TMath::Pi(); - tgl=(-1)*seed2; - } - if(seed1>0.5*TMath::Pi() && seed1-TMath::Pi() && seed1<-0.5*TMath::Pi()){ - phi=seed1+0.5*TMath::Pi(); - tgl=(-1)*seed2; - } - if(seed1>-0.5*TMath::Pi() && seed1<0.5*TMath::Pi()){ - phi=seed1-0.5*TMath::Pi(); - tgl=seed2; - } - if(seed1>0.5*TMath::Pi() && seed1GetModuleId(module1,layer,ladder,detector); - Float_t yclu1 = p1->GetY(); - Float_t zclu1 = p1->GetZ(); - Double_t cv=Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2); - - for(Int_t l4=0;l4At(l4); - for(Int_t l5=0;l5At(l5); - for(Int_t l6=0;l6At(l6); - AliITStrackSA* trac = new AliITStrackSA(layer,ladder,detector,yclu1,zclu1,phi,tgl,cv,1); - - if(cl5!=0) trac->AddClusterV2(5,(clind5[l6] & 0x0fffffff)>>0); - if(cl4!=0) trac->AddClusterV2(4,(clind4[l5] & 0x0fffffff)>>0); - if(cl3!=0) trac->AddClusterV2(3,(clind3[l4] & 0x0fffffff)>>0); - if(cl2!=0) trac->AddClusterV2(2,(clind2[l3] & 0x0fffffff)>>0); - if(cl1!=0) trac->AddClusterV2(1,(clind1[l2] & 0x0fffffff)>>0); - if(cl0!=0) trac->AddClusterV2(0,(clind0[l1] & 0x0fffffff)>>0); - - //fit with Kalman filter using AliITStrackerV2::RefitAt() - - AliITStrackV2* ot = new AliITStrackV2(*trac); - - ot->ResetCovariance(); - ot->ResetClusters(); - - if(RefitAt(49.,ot,trac)){ //fit from layer 1 to layer 6 - otrack2 = new AliITStrackV2(*ot); - otrack2->ResetCovariance(); - otrack2->ResetClusters(); - //fit from layer 6 to layer 1 - if(RefitAt(3.7,otrack2,ot)) fListOfTracks->AddLast(otrack2); + Int_t firstLay=-1,secondLay=-1; + for(Int_t i=0;iGetDetectorIndex()+firstmod[firstLay]; + Int_t layer,ladder,detector; + AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector); + Double_t yclu1 = p1->GetY(); + Double_t zclu1 = p1->GetZ(); + + Double_t x1,y1,z1; + Double_t x2,y2,z2; + Double_t cv=0,tgl2=0,phi2=0; + AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1); + x1 = arr1->GetX(); + y1 = arr1->GetY(); + z1 = arr1->GetZ(); + + if(secondLay>0) { + Int_t iSecondLay=indices[secondLay]; + Int_t mrk2=clmark[secondLay][iSecondLay]; + AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2); + x2 = arr2->GetX(); + y2 = arr2->GetY(); + z2 = arr2->GetZ(); + cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2); + tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); + phi2 = TMath::ATan2((y2-y1),(x2-x1)); + } else { // special case of 1-point tracks, only for cosmics (B=0) + x2 = primaryVertex[0]; + y2 = primaryVertex[1]; + z2 = primaryVertex[2]; + cv = 0; + tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); + phi2 = TMath::ATan2((y1-y2),(x1-x2)); + } + + // create track and attach it the RecPoints + AliITStrackSA trac(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1); + for(Int_t iLay=5; iLay>=0; iLay--){ + Int_t iInLay=indices[iLay]; + AliITSRecPoint* cl=(AliITSRecPoint*)listlayer[iLay][iInLay]; + if(cl!=0){ + trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0); + trac.AddClusterMark(iLay,clmark[iLay][iInLay]); + } + } + + //fit with Kalman filter using AliITStrackerMI::RefitAt() + AliITStrackSA ot(trac); + + ot.ResetCovariance(10.); + ot.ResetClusters(); + + // Propagate inside the innermost layer with a cluster + if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) { + + if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&ot,&trac)){ //fit from layer 1 to layer 6 + AliITStrackMI otrack2(ot); + otrack2.ResetCovariance(10.); + otrack2.ResetClusters(); + //fit from layer 6 to layer 1 + if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),&otrack2,&ot)) { + new(arrMI[nFoundTracks]) AliITStrackMI(otrack2); + new(arrSA[nFoundTracks]) AliITStrackSA(trac); + ++nFoundTracks; + } - } - - delete ot; - delete trac; + } + } }//end loop layer 6 }//end loop layer 5 - }//end loop layer 4 - - for(Int_t i=1;i<3;i++){ - delete recp[i]; - delete errs[i]; - } + }//end loop layer 4 }//end loop layer 3 }//end loop layer 2 - delete recp[0]; - delete errs[0]; - delete[] recp; - delete[] errs; }//end loop layer 1 - delete [] end; - Int_t dim=fListOfTracks->GetEntries(); - if(dim==0){ - for(Int_t i=0;iGetNlayers();i++){ - delete listlayer[i]; - } - delete listlayer; - return 0; - } - AliITStrackV2* otrack =(AliITStrackV2*)FindTrackLowChiSquare(fListOfTracks,dim); + if(fListOfTracks->GetEntries()==0) return 0; + + Int_t lowchi2 = FindTrackLowChiSquare(); + AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2); + AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2); + if(otrack==0) return 0; - Int_t * indexc = new Int_t[fGeom->GetNlayers()]; - for(Int_t i=0;iGetNlayers();i++) indexc[i]=0; - for(Int_t nind=0;nindGetNumberOfClusters();nind++){ - indexc[nind] = otrack->GetClusterIndex(nind); - } - AliITSclusterV2* cl0 = (AliITSclusterV2*)GetCluster(indexc[0]); - AliITSclusterV2* cl1 = (AliITSclusterV2*)GetCluster(indexc[1]); - AliITSclusterV2* cl2 = (AliITSclusterV2*)GetCluster(indexc[2]); - AliITSclusterV2* cl3 = (AliITSclusterV2*)GetCluster(indexc[3]); - AliITSclusterV2* cl4 = (AliITSclusterV2*)GetCluster(indexc[4]); - Int_t labl[3]={-1,-1,-1}; - if(otrack->GetNumberOfClusters()==fGeom->GetNlayers()){ - AliITSclusterV2* cl5 = (AliITSclusterV2*)GetCluster(indexc[5]); - labl[0]=cl5->GetLabel(0); - labl[1]=cl5->GetLabel(1); - labl[2]=cl5->GetLabel(2); - } - delete [] indexc; - if(otrack->GetNumberOfClusters()==(fGeom->GetNlayers()-1)){ - labl[0]=-1; - labl[1]=-1; - labl[2]=-1; - } - Int_t numberofpoints; - if(fSixPoints) numberofpoints=6; - else numberofpoints=5; - Int_t label = Label(cl0->GetLabel(0),cl1->GetLabel(0), - cl2->GetLabel(0),cl3->GetLabel(0), - cl4->GetLabel(0),labl[0], - cl0->GetLabel(1),cl1->GetLabel(1), - cl2->GetLabel(1),cl3->GetLabel(1), - cl4->GetLabel(1),labl[1], - cl0->GetLabel(2),cl1->GetLabel(2), - cl2->GetLabel(2),cl3->GetLabel(2), - cl4->GetLabel(2),labl[2],numberofpoints); - + + CookLabel(otrack,0.); //MI change - to see fake ratio + Int_t label=FindLabel(otrack); otrack->SetLabel(label); - for(Int_t i=0;iGetNlayers();i++){ - delete listlayer[i]; + Double_t low=0.; + Double_t up=0.51; + otrack->CookdEdx(low,up); + + //remove clusters of found track + for(Int_t nlay=0;nlayGetNumberOfMarked(nlay);cln++){ + Int_t index = trsa->GetClusterMark(nlay,cln); + RemoveClusterCoord(nlay,index); + } } - delete listlayer; - delete [] firstmod; + return otrack; } -//_______________________________________________________________________ -void AliITStrackerSA::UseFoundTracksV2(Int_t evnum,TTree* treev2){ - // Marks as used clusters belonging to tracks found with V2 TPC+ITS tracking - //(or AliITStrackV2 tracks found with function FindTracks of this class) - - - //Get primary vertex - if(fVertexer){ - if(fVert)delete fVert; - fVert = fVertexer->FindVertexForCurrentEvent(evnum); +//_______________________________________________________ +void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) const +{ + // + // Add new track to the ESD + // + AliESDtrack outtrack; + outtrack.UpdateTrackParams(t,AliESDtrack::kITSin); + if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA); + for(Int_t i=0;i<12;i++) { + outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i)); } - else { - gAlice->GetEvent(evnum); - if(!fVert){ - Fatal("FindTracks","Vertex is missing\n"); - return; + Double_t sdedx[4]={0.,0.,0.,0.}; + for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i); + outtrack.SetITSdEdxSamples(sdedx); + + + if(AliITSReconstructor::GetRecoParam()->GetSAUsedEdxInfo()){ + Double_t mom=t->P(); + Double_t ppid[AliPID::kSPECIES]; + for(Int_t isp=0;ispGetdEdx(); + Int_t ide=fITSPid->GetParticleIdFromdEdxVsP(mom,truncmean,kTRUE); + if(ide==AliPID::kProton){ + ppid[AliPID::kProton]=1.; + ppid[AliPID::kPion]=0.; + } + else if(ide==AliPID::kKaon){ + ppid[AliPID::kKaon]=1.; + ppid[AliPID::kPion]=0.; + } } + outtrack.SetITSpid(ppid); + outtrack.SetESDpid(ppid); } - Double_t primaryVertex[3]; - fVert->GetXYZ(primaryVertex); - - if(!fTable){ - fTable = new AliITSclusterTable(fGeom,this,primaryVertex); - fTable->FillArray(fITSclusters); - fTable->FillArrayCoorAngles(); - } + event->AddTrack(&outtrack); - TBranch* bra = (TBranch*)treev2->GetBranch("tracks"); - if(!bra) Warning("UseFoundTracksV2","No branch for track tree"); - AliITStrackV2* ttrrt = new AliITStrackV2; - bra->SetAddress(&ttrrt); - - for(Int_t nj=0;njGetEntries();nj++){ - treev2->GetEvent(nj); - Int_t ncl = ttrrt->GetNumberOfClusters(); - for(Int_t k=0;kGetClusterIndex(k); - AliITSclusterV2* clui = (AliITSclusterV2*)GetCluster(index); - if(clui->IsUsed()==0) clui->Use(); - } - } - delete ttrrt; - + return; } -//_______________________________________________________________________ -void AliITStrackerSA::UseFoundTracksV2(AliESD *event){ - // Marks as used clusters belonging to tracks found with V2 TPC+ITS tracking - //Get primary vertex - - Double_t primaryVertex[3]; - event->GetVertex()->GetXYZ(primaryVertex); - - if(!fTable){ - fTable = new AliITSclusterTable(fGeom,this,primaryVertex); - fTable->FillArray(fITSclusters); - fTable->FillArrayCoorAngles(); - } +//_______________________________________________________ +Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){ + //function used to to find the clusters associated to the track - Int_t ntracks = event->GetNumberOfTracks(); - while (ntracks--) { - AliESDtrack *esd=event->GetTrack(ntracks); - if ((esd->GetStatus()& - AliESDtrack::kITSin|AliESDtrack::kTPCin)==0) continue; - UInt_t idx[6]; - Int_t ncl = esd->GetITSclusters(idx); - for(Int_t clu=0; cluIsUsed()==0) cl->Use(); - } - } - - Info("UseFoundTracksV2","Clusters of tracks prolonged from TPC deleted"); - + if(ForceSkippingOfLayer(layer)) return 0; -} -//_______________________________________________________ -Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t zvertex,Int_t pflag,AliITSclusterTable* table){ - //function used to to find the clusters associated to the track Int_t nc=0; AliITSlayer &lay = fgLayers[layer]; - Int_t * firstmod = new Int_t[fGeom->GetNlayers()]; - for(Int_t i=0;iGetNlayers();i++){ - firstmod[i]=fGeom->GetModuleIndex(i+1,1,1); - } - if(pflag==1){ - - Float_t cx1,cx2,cy1,cy2; + Double_t r=lay.GetR(); + if(pflag==1){ + Double_t cx1,cx2,cy1,cy2; FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3); - Int_t fun = FindIntersection(fCoef1,fCoef2,fCoef3,-(lay.GetR()*lay.GetR()),cx1,cy1,cx2,cy2); - if(fun==0) { - delete[] firstmod; - return 0; - } - Double_t fi1 =TMath::ATan2(cy1,cx1); - Double_t fi2 =TMath::ATan2(cy2,cx2); - fPhiEstimate = ChoosePoint(fi1,fi2,fPhic); + if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0) + return 0; + Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]); + Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]); + fPhiEstimate=ChoosePoint(fi1,fi2,fPhic); } - Double_t zed = TMath::Tan(fLambdac)*lay.GetR()+zvertex; - Double_t zed1 = TMath::Tan(fLambdac+lambdawindow)*lay.GetR()+zvertex; - Double_t zed2 = TMath::Tan(fLambdac-lambdawindow)*lay.GetR()+zvertex; - - Double_t fi = fPhiEstimate; - Int_t nmod = lay.FindDetectorIndex(fi,zed); - if (nmod < 0) { - delete[] firstmod; - return 0; - } - nmod += firstmod[layer]; - - Int_t nm[8]={0,0,0,0,0,0,0,0}; - nm[0] = lay.FindDetectorIndex(fi+phiwindow,zed); - nm[1] = lay.FindDetectorIndex(fi-phiwindow,zed); - nm[2] = lay.FindDetectorIndex(fi,zed1); - nm[3] = lay.FindDetectorIndex(fi,zed2); - nm[4] = lay.FindDetectorIndex(fi+phiwindow,zed1); - nm[5] = lay.FindDetectorIndex(fi-phiwindow,zed1); - nm[6] = lay.FindDetectorIndex(fi+phiwindow,zed2); - nm[7] = lay.FindDetectorIndex(fi-phiwindow,zed2); - - - Int_t nn=0; - TArrayI* array =(TArrayI*)table->GetListOfClusters(nmod); - TArrayI* listc = new TArrayI(array->GetSize()); - for(Int_t i=0;iGetSize();i++){ - Int_t in=(Int_t)array->At(i); - listc->AddAt(in,nn); - nn++; - } - - Int_t k=0; - Int_t val; - while(k<8){ - for(Int_t h=k+1;h<8;h++){ - if(nm[k]>nm[h]){ - val=nm[k]; - nm[k]=nm[h]; - nm[h]=val; - } - - } - k++; - } - - Int_t value=-5; - - for(Int_t ii=0;ii<8;ii++){ - if(nm[ii]!=value && nm[ii]!=nmod && nm[ii]>=0){ - TArrayI* ar =(TArrayI*)table->GetListOfClusters(nm[ii]+firstmod[layer]); - listc->Set(listc->GetSize()+ar->GetSize()); - for(Int_t j=0;jGetSize();j++){ - Int_t in=(Int_t)ar->At(j); - listc->AddAt(in,nn); - nn++; - value=nm[ii]; - } - } - } - - - for(Int_t i=0;iGetSize();i++){ - Int_t index = (Int_t)listc->At(i); - AliITSclusterV2* cllay = lay.GetCluster(index); - if(cllay==0) continue; - if(cllay->IsUsed()==1) continue; - Double_t phi = fTable->GetPhiCluster(layer,index); - Double_t lambda= fTable->GetLambdaCluster(layer,index); - if(TMath::Abs(fLambdac-lambda)GetNumberOfClustersSA()==15){ - delete[] firstmod; - delete listc; - return 0; - } - trs->AddClusterSA(layer,index); - cllay->Use(); - fPhiEstimate=phi; - fPointc[0]=fTable->GetXCluster(layer,index); - fPointc[1]=fTable->GetYCluster(layer,index); - } + Double_t phiExpect=fPhiEstimate; + Double_t lamExpect=fLambdac; + + Int_t ncl = fCluCoord[layer]->GetEntriesFast(); + Int_t startcl=FindIndex(layer,lamExpect-lambdawindow*1.02); + Int_t endcl=FindIndex(layer,lamExpect+lambdawindow*1.02)+1; + if(endcl>=ncl) endcl=ncl-1; + + for (Int_t index=startcl; index<=endcl; index++) { + //for (Int_t index=0; indexGetPhi(); + Double_t deltaPhi = phi-phiExpect; + if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi(); + else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi(); + if (TMath::Abs(deltaPhi)>phiwindow) continue; + + Double_t lambda = arr->GetLambda(); + if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue; + + if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0; + if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0; + Int_t orind = arr->GetOrInd(); + trs->AddClusterSA(layer,orind); + trs->AddClusterMark(layer,index); + nc++; + fLambdac=lambda; + fPhiEstimate=phi; + + fPointc[0]=arr->GetX(); + fPointc[1]=arr->GetY(); + } - delete listc; - delete [] firstmod; return nc; - } - +//________________________________________________________________ +Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){ + // Sets the first point (seed) for tracking + + AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu); + fPhic = arr->GetPhi(); + fLambdac = arr->GetLambda(); + fPhiEstimate = fPhic; + fPoint1[0]=primaryVertex[0]; + fPoint1[1]=primaryVertex[1]; + fPoint2[0]=arr->GetX(); + fPoint2[1]=arr->GetY(); + return kTRUE; +} //________________________________________________________________ void AliITStrackerSA::UpdatePoints(){ //update of points for the estimation of the curvature - fPoint1[0]=fPoint2[0]; - fPoint1[1]=fPoint2[1]; fPoint2[0]=fPoint3[0]; fPoint2[1]=fPoint3[1]; fPoint3[0]=fPointc[0]; @@ -1094,36 +771,42 @@ void AliITStrackerSA::UpdatePoints(){ } - - //___________________________________________________________________ -Int_t AliITStrackerSA::FindEquation(Float_t x1, Float_t y1, Float_t x2, Float_t y2, Float_t x3, Float_t y3,Float_t& a, Float_t& b, Float_t& c){ +Int_t AliITStrackerSA::FindEquation(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3,Double_t& a, Double_t& b, Double_t& c){ //given (x,y) of three recpoints (in global coordinates) //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c - - Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1); + double dx31=x3-x1,dy31=y3-y1,dx21=x2-x1,dy21=y2-y1; + Double_t den = dx31*dy21-dx21*dy31; if(den==0) return 0; - a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den; - b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1); + else den = 1./den; + // + double r31 = -dx31*(x1+x3) - dy31*(y1+y3); + double r21 = -dx21*(x1+x2) - dy21*(y1+y2); + // + double da = r31*dy21 - r21*dy31; + double db = r21*dx31 - r31*dx21; + a = da*den; + b = db*den; c = -x1*x1-y1*y1-a*x1-b*y1; return 1; } //__________________________________________________________________________ - Int_t AliITStrackerSA::FindIntersection(Float_t a1, Float_t b1, Float_t c1, Float_t c2,Float_t& x1,Float_t& y1, Float_t& x2, Float_t& y2){ + Int_t AliITStrackerSA::FindIntersection(Double_t a1, Double_t b1, Double_t c1, Double_t c2,Double_t& x1,Double_t& y1, Double_t& x2, Double_t& y2){ //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer //c2 is -rlayer*rlayer if(a1==0) return 0; - Float_t m = c2-c1; - Float_t aA = (b1*b1)/(a1*a1)+1; - Float_t bB = (-2*m*b1/(a1*a1)); - Float_t cC = c2+(m*m)/(a1*a1); - if((bB*bB-4*aA*cC)<0) return 0; + Double_t m = c2-c1; + Double_t aA = (b1*b1)/(a1*a1)+1; + Double_t bB = (-2*m*b1/(a1*a1)); + Double_t cC = c2+(m*m)/(a1*a1); + Double_t dD = bB*bB-4*aA*cC; + if(dD<0) return 0; - y1 = (-bB+TMath::Sqrt(bB*bB-4*aA*cC))/(2*aA); - y2 = (-bB-TMath::Sqrt(bB*bB-4*aA*cC))/(2*aA); + y1 = (-bB+TMath::Sqrt(dD))/(2*aA); + y2 = (-bB-TMath::Sqrt(dD))/(2*aA); x1 = (c2-c1-b1*y1)/a1; x2 = (c2-c1-b1*y2)/a1; @@ -1153,6 +836,8 @@ x2,Double_t y2,Double_t x3,Double_t y3){ return 1/rad; } + + //____________________________________________________________________ Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){ @@ -1169,113 +854,99 @@ Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){ //_________________________________________________________________ -AliITStrackV2* AliITStrackerSA::FindTrackLowChiSquare(TObjArray* tracklist, Int_t dim) const { - // returns track with lowes chi square - if(dim==1){ - AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(0); - return trk; - } - if(dim==0) return 0; - Double_t * chi2 = new Double_t[dim]; - Int_t * index = new Int_t[dim]; - for(Int_t i=0;iAt(i); - chi2[i]=trk->GetChi2(); - index[i]=i; - } - - Int_t w=0;Double_t value; - Int_t lp; - while(wGetEntries(); + if(dim<=1) return 0; + AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0); + Double_t minChi2=trk->GetChi2(); + Int_t index=0; + for(Int_t i=1;iAt(i); + Double_t chi2=trk->GetChi2(); + if(chi2At(index[dim-1]); - delete [] chi2; - delete [] index; - return trk; - + return index; } //__________________________________________________________ -Int_t AliITStrackerSA::FindLabel(Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5, Int_t l6){ - - //function used to determine the track label - - Int_t lb[6] = {l1,l2,l3,l4,l5,l6}; - Int_t aa[6]={1,1,1,1,1,1}; - Int_t ff=0; - Int_t ll=0; - Int_t k=0;Int_t w=0;Int_t num=6; - if(lb[5]==-1) num=5; - - while(kGetNumberOfClusters(); i++) { + Int_t indexc = track->GetClusterIndex(i); + AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc); + Int_t iLayer=cl->GetLayer(); + for(Int_t k=0;k<3;k++){ + labl[iLayer][k]=cl->GetLabel(k); + if(labl[iLayer][k]<0) iNotLabel++; + } + } + if(iNotLabel==3*track->GetNumberOfClusters()) return -2; + + for(Int_t j1=0;j1=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){ + cnts[j2][k2]++; + cnts[j1][k1]=0; + } + } + } } - w++; } - if(num==6) return lb[5]; - else return lb[4]; -} - -//_____________________________________________________________________________ -Int_t AliITStrackerSA::Label(Int_t gl1, Int_t gl2, Int_t gl3, Int_t gl4, Int_t gl5, Int_t gl6,Int_t gl7, Int_t gl8, Int_t gl9, Int_t gl10,Int_t gl11, -Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t numberofpoints){ - - - //function used to assign label to the found track. If track is fake, the label is negative - Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6}; - Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12}; - Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18}; - Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]); - Int_t lflag=0;Int_t num=6; - if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5; - for(Int_t i=0;icntMax && labl[j][k]>=0){ + cntMax=cnts[j][k]; + label=labl[j][k]; + } + } } - if(lflag>=numberofpoints) return ll; - else return -ll; - + Int_t lflag=0; + for(Int_t i=0;iGetNumberOfClusters()) label = -label; + return label; } - //_____________________________________________________________________________ -void AliITStrackerSA::SetWindowSizes(Int_t n, Double_t *phi, Double_t *lam){ +void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Double_t phimin, Double_t phimax, Double_t lambdamin, Double_t lambdamax){ + // Set sizes of the phi and lambda windows used for track finding + fNloop = n; + if(fPhiWin) delete [] fPhiWin; + if(fLambdaWin) delete [] fLambdaWin; + fPhiWin = new Double_t[fNloop]; + fLambdaWin = new Double_t[fNloop]; + Double_t stepPhi=(phimax-phimin)/(Double_t)(fNloop-1); + Double_t stepLambda=(lambdamax-lambdamin)/(Double_t)(fNloop-1); + for(Int_t k=0;kGetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese) + for(Int_t k=0;kGetRotMatrix(module,rot); + Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det); + Double_t tx,ty,tz; fGeom->GetTrans(lay,lad,det,tx,ty,tz); + + Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi(); + Double_t phi1=TMath::Pi()/2+alpha; + if (lay==1) phi1+=TMath::Pi(); + + Double_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1); + Double_t r=tx*cp+ty*sp; + + xyz= r*cp - cl->GetY()*sp; + y= r*sp + cl->GetY()*cp; + z=cl->GetZ(); +*/ + Float_t xyz[3]; + cl->GetGlobalXYZ(xyz); + x=xyz[0]; + y=xyz[1]; + z=xyz[2]; + + phi=TMath::ATan2(y-vertex[1],x-vertex[0]); + lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1]))); +} +//________________________________________________________________________ +void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Double_t &sx,Double_t &sy, Double_t &sz){ + + //returns sigmax, y, z of cluster in global coordinates +/* + Double_t rot[9]; fGeom->GetRotMatrix(module,rot); + Int_t lay,lad,det; + AliITSgeomTGeo::GetModuleId(module,lay,lad,det); + + Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi(); + Double_t phi=TMath::Pi()/2+alpha; + if (lay==1) phi+=TMath::Pi(); + + Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi); +*/ + Float_t covm[6]; + cl->GetGlobalCov(covm); + sx=TMath::Sqrt(covm[0]); + sy=TMath::Sqrt(covm[3]); + sz=TMath::Sqrt(covm[5]); +/* + sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2()); + sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2()); + sz = TMath::Sqrt(cl->GetSigmaZ2()); +*/ +} + +//________________________________________________________________________ +Int_t AliITStrackerSA::FindIndex(Int_t lay, Double_t lamVal) const { + // Find the cluster at limit of lambda window + + Int_t base = 0; + Int_t last = fCluCoord[lay]->GetEntriesFast()-1; + if(last<0) return 0; + Int_t position; + Double_t lamfirst=((AliITSclusterTable*)fCluCoord[lay]->At(base))->GetLambda(); + if(lamfirst>lamVal) return base; + Double_t lamlast=((AliITSclusterTable*)fCluCoord[lay]->At(last))->GetLambda(); + if(lamlast<=lamVal) return last; + while (last >= base) { + position = (base+last) / 2; + Double_t a=((AliITSclusterTable*)fCluCoord[lay]->At(position))->GetLambda()-lamVal; + Double_t b=((AliITSclusterTable*)fCluCoord[lay]->At(position+1))->GetLambda()-lamVal; + if(a*b<=0) return position; + if(a>0) last = position; + else base = position; + } + return 0; +}