From: hristov Date: Wed, 21 Jul 2004 21:47:49 +0000 (+0000) Subject: New version of the ITS tracking, presented by M.Ivanov during the off-line week in... X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=e43c066cbc76252a79f11f5bc1435a37cbc8bed3 New version of the ITS tracking, presented by M.Ivanov during the off-line week in June04. Now the stand-alone tracker inherits from AliITStrackerMI (M.Ivanov) --- diff --git a/ITS/AliITSReconstructor.cxx b/ITS/AliITSReconstructor.cxx index fbfd6d9126e..7f6be3193bf 100644 --- a/ITS/AliITSReconstructor.cxx +++ b/ITS/AliITSReconstructor.cxx @@ -26,6 +26,7 @@ #include "AliRunLoader.h" #include "AliRawReader.h" #include "AliITSclustererV2.h" +#include "AliITStrackerMI.h" #include "AliITStrackerSA.h" #include "AliITSVertexerIons.h" #include "AliITSVertexerFast.h" @@ -126,7 +127,9 @@ AliTracker* AliITSReconstructor::CreateTracker(AliRunLoader* runLoader) const // create a ITS tracker AliITSgeom* geom = GetITSgeom(runLoader); - if (!geom) return NULL; + if (!geom) return NULL; + TString selectedTracker = GetOption(); + if (selectedTracker.Contains("MI")) return new AliITStrackerMI(geom); return new AliITStrackerSA(geom); } diff --git a/ITS/AliITSclusterV2.h b/ITS/AliITSclusterV2.h index 26b1359e621..268bafc1d9e 100644 --- a/ITS/AliITSclusterV2.h +++ b/ITS/AliITSclusterV2.h @@ -16,14 +16,16 @@ //_____________________________________________________________________________ class AliITSclusterV2 : public AliCluster { public: - AliITSclusterV2() : AliCluster() {fQ=0; fLayer=0; fNz=0; fNy=0;} + AliITSclusterV2() : AliCluster() {fQ=0; fLayer=0; fNz=1; fNy=1; fType=0;fDeltaProb=0;} AliITSclusterV2(Int_t *lab,Float_t *hit, Int_t *info) : AliCluster(lab,hit) { fIndex=lab[3]; fQ=hit[4]; fNy = info[0]; fNz = info[1]; fLayer = info[2]; - + fChargeRatio = 0; + fType=0; + fDeltaProb=0.; } void Use() {fQ=-fQ;} void SetQ(Float_t q) {fQ=q;} @@ -31,21 +33,30 @@ public: void SetLayer(Int_t layer) {fLayer=layer;} void SetNz(Int_t nz) {fNz =nz;} void SetNy(Int_t ny){fNy=ny;} + void SetChargeRatio(Float_t ratio) { fChargeRatio = ratio;} + void SetType(Int_t type){ fType=type;} + void SetDeltaProbability(Float_t prob){fDeltaProb = prob;} + Int_t IsUsed() const {return (fQ<0) ? 1 : 0;} Float_t GetQ() const {return TMath::Abs(fQ);} Int_t GetDetectorIndex() const { return 0x3FF&fIndex; } Int_t GetLayer() const {return fLayer;} Int_t GetNz() const {return fNz;} Int_t GetNy() const {return fNy;} + Float_t GetChargeRatio() const {return fChargeRatio;} Int_t GetPindex() const { return 0xFFF00000&fIndex; } //SSD clusters only Int_t GetNindex() const { return 0xFFC00&fIndex; } //SSD clusters only - + Int_t GetType() const {return fType;} // type of the cluster + Float_t GetDeltaProbability() const{return fDeltaProb;} //probability to belong to the delta ray private: Int_t fIndex; // detector index Float_t fQ ; // Q of cluster (in ADC counts) Char_t fLayer; // layer number - Char_t fNz; //number of digits in Z direction - Char_t fNy; //number of digits in y direction + Short_t fNz; //number of digits in Z direction + Short_t fNy; //number of digits in y direction + Float_t fChargeRatio; //charge ratio + Int_t fType; //quality factor of the cluster + Float_t fDeltaProb; // probability to be deleta electron ClassDef(AliITSclusterV2,2) // ITS clusters }; diff --git a/ITS/AliITSclustererV2.cxx b/ITS/AliITSclustererV2.cxx index 3a72144719f..bba23509c4f 100644 --- a/ITS/AliITSclustererV2.cxx +++ b/ITS/AliITSclustererV2.cxx @@ -51,6 +51,7 @@ AliITSclustererV2::AliITSclustererV2(const AliITSgeom *geom) { fYshift[m] = x*ca + y*sa; fZshift[m] = (Double_t)z; fNdet[m] = (lad-1)*g->GetNdetectors(lay) + (det-1); + fNlayer[m] = lay-1; } fNModules = g->GetIndexMax(); @@ -283,22 +284,51 @@ static void CheckLabels2(Int_t lab[10]) { //------------------------------------------------------------ Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); + Int_t nlabels =0; + for (Int_t i=0;i<10;i++) if (lab[i]>=0) nlabels++; + for (Int_t i=0;i<10;i++){ Int_t label = lab[i]; if (label>=0 && labelGetMCApp()->Particle(label); - if (part->P() < 0.005) { - Int_t m=part->GetFirstMother(); - if (m<0) { - continue; - } - if (part->GetStatusCode()>0) { - continue; + if (part->P() < 0.02) { + Int_t m=part->GetFirstMother(); + if (m<0) { + continue; + } + if (part->GetStatusCode()>0) { + continue; + } + lab[i]=m; + } + else + if (part->P() < 0.12 && nlabels>3) { + lab[i]=-2; + nlabels--; + } + } + else{ + if ( (label>ntracks||label <0) && nlabels>3) { + lab[i]=-2; + nlabels--; + } + } + } + if (nlabels>3){ + for (Int_t i=0;i<10;i++){ + if (nlabels>3){ + Int_t label = lab[i]; + if (label>=0 && labelGetMCApp()->Particle(label); + if (part->P() < 0.1) { + lab[i]=-2; + nlabels--; + } } - lab[i]=m; } - } + } } + //compress labels -- if multi-times the same Int_t lab2[10]; for (Int_t i=0;i<10;i++) lab2[i]=-2; @@ -316,6 +346,8 @@ static void CheckLabels2(Int_t lab[10]) { } static void AddLabel(Int_t lab[10], Int_t label) { + Int_t ntracks = gAlice->GetMCApp()->GetNtrack(); + if (label>ntracks) return; for (Int_t i=0;i<10;i++){ if (label<0) break; if (lab[i]==label) break; @@ -455,21 +487,15 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { Info("FindClustersSPD","Too big cluster !"); continue; } - Int_t lab[4]; - lab[0]=-2; - lab[1]=-2; - lab[2]=-2; - lab[3]=fNdet[fI]; Int_t milab[10]; for (Int_t ilab=0;ilab<10;ilab++){ milab[ilab]=-2; } - d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); Int_t ymin=d->GetCoord2(),ymax=ymin; Int_t zmin=d->GetCoord1(),zmax=zmin; - Float_t y=0.,z=0.,q=0.; + for (Int_t l=0; lUncheckedAt(idx[l]); @@ -477,55 +503,65 @@ FindClustersSPD(const TClonesArray *digits, TClonesArray *clusters) { if (ymax < d->GetCoord2()) ymax=d->GetCoord2(); if (zmin > d->GetCoord1()) zmin=d->GetCoord1(); if (zmax < d->GetCoord1()) zmax=d->GetCoord1(); - - Int_t lab0=(d->GetTracks())[0]; - if (lab0>=0) { - if (lab[0]<0) { - lab[0]=lab0; - } else if (lab[1]<0) { - if (lab0!=lab[0]) lab[1]=lab0; - } else if (lab[2]<0) { - if (lab0!=lab[0]) - if (lab0!=lab[1]) lab[2]=lab0; - } - } - Float_t qq=d->GetSignal(); - y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; - // MI addition - find all labels - for (Int_t dlab=0;dlab<3;dlab++){ + // MI addition - find all labels in cluster + for (Int_t dlab=0;dlab<10;dlab++){ Int_t digitlab = (d->GetTracks())[dlab]; if (digitlab<0) continue; - for (Int_t index=0;index<10;index++){ - if (milab[index]<0) { - milab[index] = digitlab; - break; - } - if (milab[index]==digitlab) break; - } + AddLabel(milab,digitlab); } - } - y/=q; z/=q; - y-=fHwSPD; z-=fHlSPD; - - Float_t lp[5]; - lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0]; - lp[1]= -z+fZshift[fI]; - // Float_t factor=TMath::Max(double(ni-3.),1.5); - Float_t factor=1.5; - lp[2]= (fYpitchSPD*fYpitchSPD/12.)*factor; - lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.)*factor; - //lp[4]= q; - lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); - - CheckLabels(lab); - CheckLabels2(milab); + if (milab[9]>0) CheckLabels2(milab); + } CheckLabels2(milab); - milab[3]=fNdet[fI]; - d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); - Int_t info[3] = {ni,0,1}; - new (cl[n]) AliITSclusterV2(milab,lp,info); n++; + // + //Int_t idy = (fNlayer[fI]==0)? 2:3; + //for (Int_t iz=zmin; iz<=zmax;iz+=2) + //Int_t idy = (ymax-ymin)/4.; // max 2 clusters + Int_t idy = 0; // max 2 clusters + if (fNlayer[fI]==0 &&idy<3) idy=3; + if (fNlayer[fI]==1 &&idy<4) idy=4; + Int_t idz =3; + for (Int_t iz=zmin; iz<=zmax;iz+=idz) + for (Int_t iy=ymin; iy<=ymax;iy+=idy){ + // + Int_t ndigits =0; + Double_t y=0.,z=0.,q=0.; + for (Int_t l=0; lUncheckedAt(idx[l]); + if (zmax-zmin>=idz || ymax-ymin>=idy){ + if (TMath::Abs( d->GetCoord2()-iy)>0.75*idy) continue; + if (TMath::Abs( d->GetCoord1()-iz)>0.75*idz) continue; + } + ndigits++; + Double_t qq=d->GetSignal(); + y+=qq*fYSPD[d->GetCoord2()]; z+=qq*fZSPD[d->GetCoord1()]; q+=qq; + + } + if (ndigits==0) continue; + y/=q; z/=q; + y-=fHwSPD; z-=fHlSPD; + + Float_t lp[5]; + lp[0]=-(-y+fYshift[fI]); if (fI<=fLastSPD1) lp[0]=-lp[0]; + lp[1]= -z+fZshift[fI]; + // Float_t factor=TMath::Max(double(ni-3.),1.5); + Float_t factory=TMath::Max(ymax-ymin,1); + Float_t factorz=TMath::Max(zmax-zmin,1); + factory*= factory; + factorz*= factorz; + //lp[2]= (fYpitchSPD*fYpitchSPD/12.)*factory; + //lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.)*factorz; + lp[2]= (fYpitchSPD*fYpitchSPD/12.); + lp[3]= (fZ1pitchSPD*fZ1pitchSPD/12.); + //lp[4]= q; + lp[4]= (zmax-zmin+1)*100 + (ymax-ymin+1); + + milab[3]=fNdet[fI]; + d=(AliITSdigitSPD*)digits->UncheckedAt(idx[0]); + Int_t info[3] = {ymax-ymin+1,zmax-zmin+1,fNlayer[fI]}; + new (cl[n]) AliITSclusterV2(milab,lp,info); n++; + } } - + delete [] bins; } @@ -719,6 +755,7 @@ FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, FindPeaks(i, nzBins, bins[s], idx, msk, npeaks); if (npeaks>30) continue; + if (npeaks==0) continue; Int_t k,l; for (k=0; kGetPadPitchWidth(sec); - c.SetSigmaY2((s2 + 1./12.)*w*w); - if (s2 != 0.) { - c.SetSigmaY2(c.GetSigmaY2()*0.108); - if (secGetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); - } - - s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); - w=par->GetZWidth(); - c.SetSigmaZ2((s2 + 1./12.)*w*w); - if (s2 != 0.) { - c.SetSigmaZ2(c.GetSigmaZ2()*0.169); - if (secGetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); - } + Int_t maxi=0,mini=0,maxj=0,minj=0; + //AliBin *bmax=&bins[s][idx[k]]; + //Float_t max = TMath::Max(TMath::Abs(bmax->GetQ())/5.,3.); + Float_t max=3; + for (Int_t di=-2; di<=2;di++) + for (Int_t dj=-3;dj<=3;dj++){ + Int_t index = idx[k]+di+dj*nzBins; + if (index<0) continue; + if (index>=nMaxBin) continue; + AliBin *b=&bins[s][index]; + if (TMath::Abs(b->GetQ())>max){ + if (di>maxi) maxi=di; + if (dimaxj) maxj=dj; + if (djUncheckedAt(b->GetIndex()); + for (Int_t itrack=0;itrack<10;itrack++){ + Int_t track = (d->GetTracks())[itrack]; + if (track>=0) { + AddLabel(milab, track); + } + } + } + } + } + + /* + Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY(); + Float_t w=par->GetPadPitchWidth(sec); + c.SetSigmaY2(s2); + if (s2 != 0.) { + c.SetSigmaY2(c.GetSigmaY2()*0.108); + if (secGetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07); + } + s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ(); + w=par->GetZWidth(); + c.SetSigmaZ2(s2); + + if (s2 != 0.) { + c.SetSigmaZ2(c.GetSigmaZ2()*0.169); + if (secGetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77); + } */ c.SetSigmaY2(0.0030*0.0030); @@ -779,7 +843,12 @@ FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, Float_t y=c.GetY(),z=c.GetZ(), q=c.GetQ(); y/=q; z/=q; - + // + //Float_t s2 = c.GetSigmaY2()/c.GetQ() - y*y; + // c.SetSigmaY2(s2); + //s2 = c.GetSigmaZ2()/c.GetQ() - z*z; + //c.SetSigmaZ2(s2); + // y=(y-0.5)*fYpitchSDD; y-=fHwSDD; y-=fYoffSDD; //delay ? @@ -792,111 +861,29 @@ FindClustersSDD(AliBin* bins[2], Int_t nMaxBin, Int_t nzBins, z= -z+fZshift[fI]; c.SetY(y); c.SetZ(z); - + c.SetNy(maxj-minj+1); + c.SetNz(maxi-mini+1); + c.SetType(npeaks); c.SetQ(q/12.7); //to be consistent with the SSD charges if (c.GetQ() < 20.) continue; //noise cluster - - if (digits) { - AliBin *b=&bins[s][idx[k]]; - AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - Int_t l0=(d->GetTracks())[0]; - //if (l0<0) { - b=&bins[s][idx[k]-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - //} - //if (l0<0) { - b=&bins[s][idx[k]+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - // } - //if (l0<0) { - b=&bins[s][idx[k]-nzBins]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - //} - //if (l0<0) { - b=&bins[s][idx[k]+nzBins]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - //} - - //if (l0<0) { - b=&bins[s][idx[k]+nzBins+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - //} - //if (l0<0) { - b=&bins[s][idx[k]+nzBins-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - //} - //if (l0<0) { - b=&bins[s][idx[k]-nzBins+1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - //} - //if (l0<0) { - b=&bins[s][idx[k]-nzBins-1]; - if (b->GetQ()>0) { - d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); - l0=(d->GetTracks())[0]; - AddLabel(milab, (d->GetTracks())[0]); - AddLabel(milab, (d->GetTracks())[1]); - AddLabel(milab, (d->GetTracks())[2]); - } - //} - + + if (digits) { + // AliBin *b=&bins[s][idx[k]]; + // AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(b->GetIndex()); { - Int_t lab[3]; - lab[0]=(d->GetTracks())[0]; - lab[1]=(d->GetTracks())[1]; - lab[2]=(d->GetTracks())[2]; - CheckLabels(lab); + //Int_t lab[3]; + //lab[0]=(d->GetTracks())[0]; + //lab[1]=(d->GetTracks())[1]; + //lab[2]=(d->GetTracks())[2]; + //CheckLabels(lab); CheckLabels2(milab); c.SetLabel(milab[0],0); c.SetLabel(milab[1],1); c.SetLabel(milab[2],2); - c.SetLayer(3); + c.SetLayer(fNlayer[fI]); } } - new (cl[ncl]) AliITSclusterV2(c); ncl++; } } @@ -921,7 +908,7 @@ FindClustersSDD(const TClonesArray *digits, TClonesArray *clusters) { Int_t y=d->GetCoord2()+1; //y Int_t z=d->GetCoord1()+1; //z Int_t q=d->GetSignal(); - + if (q<3) continue; if (z <= fNzSDD) { bins[0][y*kNzBins+z].SetQ(q); @@ -987,6 +974,8 @@ void AliITSclustererV2::FindClustersSDD(AliITSRawStream* input, Info("FindClustersSDD", "found clusters in ITS SDD: %d", nClustersSDD); } + + void AliITSclustererV2:: FindClustersSSD(Ali1Dcluster* neg, Int_t nn, Ali1Dcluster* pos, Int_t np, @@ -995,90 +984,462 @@ FindClustersSSD(Ali1Dcluster* neg, Int_t nn, // Actual SSD cluster finder //------------------------------------------------------------ TClonesArray &cl=*clusters; - - Int_t lab[4]={-2,-2,-2,-2}; + // Float_t tanp=fTanP, tann=fTanN; if (fI>fLastSSD1) {tann=fTanP; tanp=fTanN;} - Int_t idet=fNdet[fI]; Int_t ncl=0; + // + Int_t negativepair[30000]; + Int_t cnegative[3000]; + Int_t cused1[3000]; + Int_t positivepair[30000]; + Int_t cpositive[3000]; + Int_t cused2[3000]; + for (Int_t i=0;i<3000;i++) {cnegative[i]=0; cused1[i]=0;} + for (Int_t i=0;i<3000;i++) {cpositive[i]=0; cused2[i]=0;} + Short_t pairs[1000][1000]; + memset(pairs,0,sizeof(Short_t)*1000000); + // + // find available pairs + // + for (Int_t i=0; iSetChargeRatio(ratio); + cl2->SetType(1); + pairs[ip][j]=1; + if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster + cl2->SetType(2); + pairs[ip][j]=2; + } + cused1[ip]++; + cused2[j]++; + } + } + + for (Int_t ip=0;ipSetChargeRatio(ratio); + cl2->SetType(5); + pairs[ip][in] = 5; + if ((pos[ip].GetNd()+neg[in].GetNd())>6){ //multi cluster + cl2->SetType(6); + pairs[ip][in] = 6; + } + } + // + // add second pair + + // if (!(cused1[ip2] || cused2[in])){ // + if (pairs[ip2][in]==100){ + Float_t yp=pos[ip2].GetY()*fYpitchSSD; + Float_t yn=neg[in].GetY()*fYpitchSSD; + Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); + Float_t yt=yn + tann*zt; + zt-=fHlSSD; yt-=fHwSSD; + ybest =yt; zbest=zt; + qbest =pos[ip2].GetQ(); + lp[0]=-(-ybest+fYshift[fI]); + lp[1]= -zbest+fZshift[fI]; + lp[2]=0.0025*0.0025; //SigmaY2 + lp[3]=0.110*0.110; //SigmaZ2 + + lp[4]=qbest; //Q + for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[ip2].GetLabel(ilab); + milab[ilab+3] = neg[in].GetLabel(ilab); + } + // + CheckLabels2(milab); + ratio = (pos[ip2].GetQ()-neg[in].GetQ())/(pos[ip2].GetQ()+neg[in].GetQ()); + milab[3]=(((ip2<<10) + in)<<10) + idet; // pos|neg|det + Int_t info[3] = {pos[ip2].GetNd(),neg[in].GetNd(),fNlayer[fI]}; + AliITSclusterV2 *cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + cl2->SetType(5); + pairs[ip2][in] =5; + if ((pos[ip2].GetNd()+neg[in].GetNd())>6){ //multi cluster + cl2->SetType(6); + pairs[ip2][in] =6; + } + } + cused1[ip]++; + cused1[ip2]++; + cused2[in]++; + } + } + } + + // + for (Int_t jn=0;jnSetChargeRatio(ratio); + cl2->SetType(7); + pairs[ip][jn] =7; + if ((pos[ip].GetNd()+neg[jn].GetNd())>6){ //multi cluster + cl2->SetType(8); + pairs[ip][jn]=8; + } + } + // + // add second pair + // if (!(cused1[ip]||cused2[jn2])){ + if (pairs[ip][jn2]==100){ + Float_t yn=neg[jn2].GetY()*fYpitchSSD; + Double_t yp=pos[ip].GetY()*fYpitchSSD; + Double_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); + Double_t yt=yn + tann*zt; + zt-=fHlSSD; yt-=fHwSSD; + ybest =yt; zbest=zt; + qbest =neg[jn2].GetQ(); + lp[0]=-(-ybest+fYshift[fI]); + lp[1]= -zbest+fZshift[fI]; + lp[2]=0.0025*0.0025; //SigmaY2 + lp[3]=0.110*0.110; //SigmaZ2 + + lp[4]=qbest; //Q + for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[ip].GetLabel(ilab); + milab[ilab+3] = neg[jn2].GetLabel(ilab); + } + // + CheckLabels2(milab); + ratio = (pos[ip].GetQ()-neg[jn2].GetQ())/(pos[ip].GetQ()+neg[jn2].GetQ()); + milab[3]=(((ip<<10) + jn2)<<10) + idet; // pos|neg|det + Int_t info[3] = {pos[ip].GetNd(),neg[jn2].GetNd(),fNlayer[fI]}; + AliITSclusterV2* cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + pairs[ip][jn2]=7; + cl2->SetType(7); + if ((pos[ip].GetNd()+neg[jn2].GetNd())>6){ //multi cluster + cl2->SetType(8); + pairs[ip][jn2]=8; + } + } + cused1[ip]++; + cused2[jn]++; + cused2[jn2]++; + } + } + } + + for (Int_t ip=0;ip1) continue; // more than one "proper" cluster for positive + // + count =0; + for (Int_t dj=0;dj1) continue; // more than one "proper" cluster for negative + + Int_t jp = 0; + + count =0; + for (Int_t dj=0;dj1) continue; + if (pairs[ip][j]<100) continue; + // + //almost gold clusters + Float_t yp=pos[ip].GetY()*fYpitchSSD; + Float_t yn=neg[j].GetY()*fYpitchSSD; + Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); + Float_t yt=yn + tann*zt; + zt-=fHlSSD; yt-=fHwSSD; + ybest=yt; zbest=zt; + qbest=0.5*(pos[ip].GetQ()+neg[j].GetQ()); + lp[0]=-(-ybest+fYshift[fI]); + lp[1]= -zbest+fZshift[fI]; + lp[2]=0.0025*0.0025; //SigmaY2 + lp[3]=0.110*0.110; //SigmaZ2 + lp[4]=qbest; //Q + for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[ip].GetLabel(ilab); + milab[ilab+3] = neg[j].GetLabel(ilab); + } + // + CheckLabels2(milab); + ratio = (pos[ip].GetQ()-neg[j].GetQ())/(pos[ip].GetQ()+neg[j].GetQ()); + milab[3]=(((ip<<10) + j)<<10) + idet; // pos|neg|det + Int_t info[3] = {pos[ip].GetNd(),neg[j].GetNd(),fNlayer[fI]}; + AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + cl2->SetType(10); + pairs[ip][j]=10; + if ((pos[ip].GetNd()+neg[j].GetNd())>6){ //multi cluster + cl2->SetType(11); + pairs[ip][j]=11; + } + cused1[ip]++; + cused2[j]++; + } + + } + + // for (Int_t i=0; i0 &&pairs[i][j]<100) continue; + ratio = (pos[i].GetQ()-neg[j].GetQ())/(pos[i].GetQ()+neg[j].GetQ()); Float_t yn=neg[j].GetY()*fYpitchSSD; Float_t zt=(2*fHlSSD*tanp + yp - yn)/(tann+tanp); Float_t yt=yn + tann*zt; zt-=fHlSSD; yt-=fHwSSD; if (TMath::Abs(yt) 4) { - lp[2]*=9; - lp[3]*=9; - } + lp[2]=0.0025*0.0025; //SigmaY2 + lp[3]=0.110*0.110; //SigmaZ2 + lp[4]=qbest; //Q - Int_t milab[10]; for (Int_t ilab=0;ilab<10;ilab++) milab[ilab]=-2; - milab[0]=pos[i].GetLabel(0); - milab[1]=neg[j].GetLabel(0); - milab[2]=pos[i].GetLabel(1); - milab[3]=neg[j].GetLabel(1); - milab[4]=pos[i].GetLabel(2); - milab[5]=neg[j].GetLabel(2); + for (Int_t ilab=0;ilab<3;ilab++){ + milab[ilab] = pos[i].GetLabel(ilab); + milab[ilab+3] = neg[j].GetLabel(ilab); + } // - CheckLabels(lab); CheckLabels2(milab); milab[3]=(((i<<10) + j)<<10) + idet; // pos|neg|det - Int_t info[3] = {0,0,5}; - new (cl[ncl]) AliITSclusterV2(milab,lp,info); ncl++; + Int_t info[3] = {pos[i].GetNd(),neg[j].GetNd(),fNlayer[fI]}; + AliITSclusterV2 * cl2 = new (cl[ncl]) AliITSclusterV2(milab,lp,info); + ncl++; + cl2->SetChargeRatio(ratio); + cl2->SetType(100+cpositive[j]+cnegative[i]); + //cl2->SetType(0); + /* + if (pairs[i][j]<100){ + printf("problem:- %d\n", pairs[i][j]); + } + if (cnegative[i]<2&&cpositive[j]<2){ + printf("problem:- %d\n", pairs[i][j]); + } + */ } } - /* - if (ybest<100) { - lab[3]=idet; - Float_t lp[5]; - lp[0]=-ybest-fYshift[fI]; - lp[1]= zbest+fZshift[fI]; - lp[2]=0.002*0.002; //SigmaY2 - lp[3]=0.080*0.080; //SigmaZ2 - lp[4]=qbest; //Q - // - new (cl[ncl]) AliITSclusterV2(lab,lp); ncl++; - } - */ } } + void AliITSclustererV2:: -FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { +FindClustersSSD(const TClonesArray *alldigits, TClonesArray *clusters) { //------------------------------------------------------------ // Actual SSD cluster finder //------------------------------------------------------------ - Int_t smax=digits->GetEntriesFast(); + Int_t smaxall=alldigits->GetEntriesFast(); + if (smaxall==0) return; + TObjArray *digits = new TObjArray; + for (Int_t i=0;iUncheckedAt(i); + if (d->GetSignal()<3) continue; + digits->AddLast(d); + } + Int_t smax = digits->GetEntriesFast(); if (smax==0) return; - + const Int_t MAX=1000; Int_t np=0, nn=0; Ali1Dcluster pos[MAX], neg[MAX]; Float_t y=0., q=0., qmax=0.; Int_t lab[4]={-2,-2,-2,-2}; - + AliITSdigitSSD *d=(AliITSdigitSSD*)digits->UncheckedAt(0); q += d->GetSignal(); y += d->GetCoord2()*d->GetSignal(); @@ -1089,28 +1450,60 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { Int_t *n=&nn; Ali1Dcluster *c=neg; Int_t nd=1; + Int_t milab[10]; + for (Int_t ilab=0;ilab<10;ilab++){ + milab[ilab]=-2; + } + milab[0]=d->GetTrack(0); milab[1]=d->GetTrack(1); milab[2]=d->GetTrack(2); + for (Int_t s=1; sUncheckedAt(s); + d=(AliITSdigitSSD*)digits->UncheckedAt(s); Int_t strip=d->GetCoord2(); if ((strip-curr) > 1 || flag!=d->GetCoord1()) { c[*n].SetY(y/q); c[*n].SetQ(q); c[*n].SetNd(nd); - c[*n].SetLabels(lab); + CheckLabels2(milab); + c[*n].SetLabels(milab); //Split suspiciously big cluster - if (nd>3) { - c[*n].SetY(y/q-0.5*nd); - c[*n].SetQ(0.5*q); - (*n)++; - if (*n==MAX) { - Error("FindClustersSSD","Too many 1D clusters !"); + /* + if (nd>10&&nd<16){ + c[*n].SetY(y/q-0.3*nd); + c[*n].SetQ(0.5*q); + (*n)++; + if (*n==MAX) { + Error("FindClustersSSD","Too many 1D clusters !"); return; - } - c[*n].SetY(y/q+0.5*nd); - c[*n].SetQ(0.5*q); - c[*n].SetNd(nd); - c[*n].SetLabels(lab); - } + } + c[*n].SetY(y/q-0.0*nd); + c[*n].SetQ(0.5*q); + c[*n].SetNd(nd); + (*n)++; + if (*n==MAX) { + Error("FindClustersSSD","Too many 1D clusters !"); + return; + } + // + c[*n].SetY(y/q+0.3*nd); + c[*n].SetQ(0.5*q); + c[*n].SetNd(nd); + c[*n].SetLabels(milab); + } + else{ + */ + if (nd>4&&nd<25) { + c[*n].SetY(y/q-0.25*nd); + c[*n].SetQ(0.5*q); + (*n)++; + if (*n==MAX) { + Error("FindClustersSSD","Too many 1D clusters !"); + return; + } + c[*n].SetY(y/q+0.25*nd); + c[*n].SetQ(0.5*q); + c[*n].SetNd(nd); + c[*n].SetLabels(milab); + } (*n)++; if (*n==MAX) { Error("FindClustersSSD","Too many 1D clusters !"); @@ -1119,6 +1512,11 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { y=q=qmax=0.; nd=0; lab[0]=lab[1]=lab[2]=-2; + // + for (Int_t ilab=0;ilab<10;ilab++){ + milab[ilab]=-2; + } + // if (flag!=d->GetCoord1()) { n=&np; c=pos; } } flag=d->GetCoord1(); @@ -1129,6 +1527,9 @@ FindClustersSSD(const TClonesArray *digits, TClonesArray *clusters) { qmax=d->GetSignal(); lab[0]=d->GetTrack(0); lab[1]=d->GetTrack(1); lab[2]=d->GetTrack(2); } + for (Int_t ilab=0;ilab<10;ilab++) { + if (d->GetTrack(ilab)>=0) AddLabel(milab, (d->GetTrack(ilab))); + } curr=strip; } c[*n].SetY(y/q); diff --git a/ITS/AliITSclustererV2.h b/ITS/AliITSclustererV2.h index f6e469bd091..af3b63aca2b 100644 --- a/ITS/AliITSclustererV2.h +++ b/ITS/AliITSclustererV2.h @@ -95,6 +95,7 @@ private: Float_t fYshift[2200]; //y-shifts of detector local coor. systems Float_t fZshift[2200]; //z-shifts of detector local coor. systems Int_t fNdet[2200]; //detector index + Int_t fNlayer[2200]; //detector layer Int_t fNModules; // total number of modules //SPD related values: diff --git a/ITS/AliITSrecoV2.h b/ITS/AliITSrecoV2.h index 714794db307..87d7ba04bf0 100644 --- a/ITS/AliITSrecoV2.h +++ b/ITS/AliITSrecoV2.h @@ -12,13 +12,16 @@ //namespace AliITSreco { const Int_t kMaxClusterPerLayer=7000*10; + const Int_t kMaxClusterPerLayer5=7000*10*2/5; + const Int_t kMaxClusterPerLayer10=7000*10*2/10; + const Int_t kMaxClusterPerLayer20=7000*10*2/20; const Int_t kMaxDetectorPerLayer=1000; const Int_t kLayersNotToSkip[]={0,0,0,0,0,0}; const Int_t kLastLayerToTrackTo=0; const Int_t kMaxLayer = 6; -const Double_t kMaxSnp = 0.6; +const Double_t kMaxSnp = 3.; const Double_t kSigmaY2[kMaxLayer]={ 1.44e-6, 1.44e-6, 1.444e-5, 1.444e-5, 4.0e-6, 4.0e-6 }; @@ -27,10 +30,18 @@ const Double_t kMaxSnp = 0.6; 1.44e-4, 1.44e-4, 7.84e-6, 7.84e-6, 0.006889, 0.006889 }; - const Double_t kChi2PerCluster=7.; - const Double_t kMaxChi2=25.; - const Double_t kMaxChi2In=16.; + const Double_t kChi2PerCluster=9.; +// const Double_t kMaxChi2PerCluster[5]={7.,5.,8.,8.,6.5}; + const Double_t kMaxChi2PerCluster[5]={11,12,12,5,12}; + const Double_t kMaxNormChi2NonC[6] = {7,8,8,11,14,25}; //max norm chi2 for non constrained tracks + const Double_t kMaxNormChi2C[6] = {11,13,15,18,30,35}; //max norm chi2 for constrained tracks + const Double_t kMaxChi2=35.; +// const Double_t kMaxChi2s[6]={40,40,40,40,40,40}; + const Double_t kMaxChi2s[6]={25,25,25,25,40,50}; + const Double_t kMaxChi2sR[6]={10,10,10,10,30,40}; + const Double_t kMaxChi2In=16.; + const Double_t kVertexCut=25; const Double_t kMaxRoad=6.0; const Double_t kXV=0.0e+0; diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx index 7b3ccd25ffa..814a66c9885 100644 --- a/ITS/AliITStrackV2.cxx +++ b/ITS/AliITStrackV2.cxx @@ -18,8 +18,6 @@ // // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch // dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch -// The class is used by AliITStrackerV2 and AliITStrackerSA classes -// (for the meaning of the track parametrization see AliITStrackV2.h) //------------------------------------------------------------------------- #include @@ -27,13 +25,14 @@ #include #include "AliCluster.h" +#include "AliTPCtrack.h" #include "AliESDtrack.h" #include "AliITStrackV2.h" -#define kWARN 5 - ClassImp(AliITStrackV2) +const Int_t kWARN=5; + //____________________________________________________________________________ AliITStrackV2::AliITStrackV2():AliKalmanTrack(), fX(0), @@ -61,15 +60,62 @@ AliITStrackV2::AliITStrackV2():AliKalmanTrack(), fC44(0), fNUsed(0), fNSkipped(0), - fReconstructed(kFALSE), + fNDeadZone(0), + fDeadZoneProbability(0), + fReconstructed(kFALSE), + fConstrain(kFALSE), fESDtrack(0) - { +{ + for(Int_t i=0; i ITS track //------------------------------------------------------------------ - for(Int_t i=0; i= TMath::Pi()) fAlpha -= 2*TMath::Pi(); + + //Conversion of the track parameters + Double_t x,p[5]; t.GetExternalParameters(x,p); + fX=x; x=GetConvConst(); + fP0=p[0]; + fP1=p[1]; + fP2=p[2]; + fP3=p[3]; + fP4=p[4]/x; + + //Conversion of the covariance matrix + Double_t c[15]; t.GetExternalCovariance(c); + + fC00=c[0 ]; + fC10=c[1 ]; fC11=c[2 ]; + fC20=c[3 ]; fC21=c[4 ]; fC22=c[5 ]; + fC30=c[6 ]; fC31=c[7 ]; fC32=c[8 ]; fC33=c[9 ]; + fC40=c[10]/x; fC41=c[11]/x; fC42=c[12]/x; fC43=c[13]/x; fC44=c[14]/x/x; + + for(Int_t i=0; i<6; i++) { fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000;} + for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0; } + fConstrain=kFALSE; + // + if (!Invariant()) throw "AliITStrackV2: conversion failed !\n"; + } //____________________________________________________________________________ @@ -82,6 +128,8 @@ AliKalmanTrack() { SetNumberOfClusters(t.GetITSclusters(fIndex)); SetLabel(t.GetLabel()); SetMass(t.GetMass()); + // + // fdEdx=t.GetITSsignal(); fAlpha = t.GetAlpha(); @@ -117,10 +165,19 @@ AliKalmanTrack() { fESDtrack=&t; fNUsed = 0; fReconstructed = kFALSE; - fNSkipped =0; - for(Int_t i=0; i<6; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0;; fChi2MIP[i]=0;} + fNSkipped =0; + fNDeadZone = 0; + fDeadZoneProbability = 0; + for(Int_t i=0; i<6; i++) {fClIndex[i]=-1; fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000;} + for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0;fChi2MIP[i]=0;} + fD[0]=0; fD[1]=0; + fExpQ=40; + fConstrain = kFALSE; + fdEdxMismatch=0; + fChi22 =0; + //if (!Invariant()) throw "AliITStrackV2: conversion failed !\n"; - SetFakeRatio(t.GetITSFakeRatio()); + } void AliITStrackV2::UpdateESDtrack(ULong_t flags) { @@ -150,14 +207,34 @@ AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) { Int_t n=GetNumberOfClusters(); for (Int_t i=0; iGet1Pt()); //Double_t c =TMath::Abs(Get1Pt()); - Double_t co=t->GetSigmaY2()*t->GetSigmaZ2()*TMath::Sqrt(TMath::Abs(fP4)); - Double_t c =GetSigmaY2()*GetSigmaZ2()*TMath::Sqrt(TMath::Abs(fP4)); + Double_t co=t->GetSigmaY2()*t->GetSigmaZ2()*(0.5+TMath::Sqrt(0.5*t->fD[0]*t->fD[0]+t->fD[1]*t->fD[1])); + Double_t c =GetSigmaY2()*GetSigmaZ2()*(0.5+TMath::Sqrt(0.5*fD[0]*fD[0]+fD[1]*fD[1])); if (c>co) return 1; else if (ckWARN) + Warning("GetPredictedChi2","Singular matrix (%d) !\n",n); + return 1e10; + } + Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01; + + Double_t dy=cy - fP0, dz=cz - fP1; + + return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; +} + //____________________________________________________________________________ Int_t AliITStrackV2::CorrectForMaterial(Double_t d, Double_t x0) { //------------------------------------------------------------------ //This function corrects the track parameters for crossed material //------------------------------------------------------------------ - Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()); - Double_t beta2=p2/(p2 + GetMass()*GetMass()); + // Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()); + Double_t p2=(1.+ fP3*fP3)/(Get1Pt()*Get1Pt()); + Double_t et = p2 + GetMass()*GetMass(); + Double_t beta2=p2/et; + et = sqrt(et); d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2)); + //d*=TMath::Sqrt(1.+ fP3*fP3 +fP2*fP2/(1.- fP2*fP2)); //Multiple scattering****************** if (d!=0) { @@ -275,10 +378,22 @@ Int_t AliITStrackV2::CorrectForMaterial(Double_t d, Double_t x0) { //Energy losses************************ if (x0!=0.) { d*=x0; - Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d; - if (beta2/(1-beta2)>3.5*3.5) + // Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d; + //Double_t dE=0; + Double_t dE = 0.265*0.153e-3*(39.2-55.6*beta2+28.7*beta2*beta2+27.41/beta2)*d; + /* + if (beta2/(1-beta2)>3.5*3.5){ dE=0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2)*d; - fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE); + } + else{ + dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d; + dE+=0.06e-3/(beta2*beta2)*d; + } + */ + fP4*=(1.- et/p2*dE); + Double_t delta44 = (dE*fP4*et/p2); + delta44*=delta44; + fC44+= delta44/400.; } if (!Invariant()) return 0; @@ -443,6 +558,84 @@ Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) { return 1; } +//____________________________________________________________________________ +Int_t AliITStrackV2::UpdateMI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz, Double_t chi2,UInt_t index) { + //------------------------------------------------------------------ + //This function updates track parameters + //------------------------------------------------------------------ + Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4; + Double_t c00=fC00; + Double_t c10=fC10, c11=fC11; + Double_t c20=fC20, c21=fC21, c22=fC22; + Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33; + Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44; + + + Double_t r00=cerry*cerry, r01=0., r11=cerrz*cerrz; + 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=cy - fP0, dz=cz - fP1; + Int_t layer = (index & 0xf0000000) >> 28; + fDy[layer] = dy; + fDz[layer] = dz; + fSigmaY[layer] = TMath::Sqrt(cerry*cerry+fC00); + fSigmaZ[layer] = TMath::Sqrt(cerrz*cerrz+fC11); + + Double_t sf=fP2 + k20*dy + k21*dz; + + fP0 += k00*dy + k01*dz; + fP1 += k10*dy + k11*dz; + fP2 = sf; + fP3 += k30*dy + k31*dz; + fP4 += k40*dy + k41*dz; + + Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40; + Double_t c12=fC21, c13=fC31, c14=fC41; + + fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11; + fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13; + fC40-=k00*c04+k01*c14; + + fC11-=k10*c01+k11*fC11; + fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13; + fC41-=k10*c04+k11*c14; + + fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13; + fC42-=k20*c04+k21*c14; + + fC33-=k30*c03+k31*c13; + fC43-=k30*c04+k31*c14; + + fC44-=k40*c04+k41*c14; + + if (!Invariant()) { + fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4; + fC00=c00; + fC10=c10; fC11=c11; + fC20=c20; fC21=c21; fC22=c22; + fC30=c30; fC31=c31; fC32=c32; fC33=c33; + fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44; + return 0; + } + + if (chi2<0) return 1; + Int_t n=GetNumberOfClusters(); + fIndex[n]=index; + SetNumberOfClusters(n+1); + SetChi2(GetChi2()+chi2); + + return 1; +} + Int_t AliITStrackV2::Invariant() const { //------------------------------------------------------------------ // This function is for debugging purpose only @@ -494,7 +687,7 @@ Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) { Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha); Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2); - TMatrixD *tT=0; + TMatrixD *T=0; // **** rotation ********************** { fAlpha = alp; @@ -502,25 +695,36 @@ Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) { fP0= -x*sa + p0*ca; fP2= sf*ca - cf*sa; - TMatrixD cC(5,5); - cC(0,0)=c00; - cC(1,0)=c10; cC(1,1)=c11; - cC(2,0)=c20; cC(2,1)=c21; cC(2,2)=c22; - cC(3,0)=c30; cC(3,1)=c31; cC(3,2)=c32; cC(3,3)=c33; - cC(4,0)=c40; cC(4,1)=c41; cC(4,2)=c42; cC(4,3)=c43; cC(4,4)=c44; - cC(0,1)=cC(1,0); - cC(0,2)=cC(2,0); cC(1,2)=cC(2,1); - cC(0,3)=cC(3,0); cC(1,3)=cC(3,1); cC(2,3)=cC(3,2); - cC(0,4)=cC(4,0); cC(1,4)=cC(4,1); cC(2,4)=cC(4,2); cC(3,4)=cC(4,3); - - TMatrixD mF(6,5); - mF(0,0)=sa; - mF(1,0)=ca; - mF(2,1)=mF(4,3)=mF(5,4)=1; - mF(3,2)=ca + sf/cf*sa; - - TMatrixD tmp(cC,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, mF)); - tT=new TMatrixD(mF,TMatrixD::kMult,tmp); + static TMatrixD C(5,5); + C(0,0)=c00; + C(1,0)=c10; C(1,1)=c11; + C(2,0)=c20; C(2,1)=c21; C(2,2)=c22; + C(3,0)=c30; C(3,1)=c31; C(3,2)=c32; C(3,3)=c33; + C(4,0)=c40; C(4,1)=c41; C(4,2)=c42; C(4,3)=c43; C(4,4)=c44; + C(0,1)=C(1,0); + C(0,2)=C(2,0); C(1,2)=C(2,1); + C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2); + C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3); + + + static TMatrixD F(6,5); + F(0,0)=sa; + F(1,0)=ca; + F(2,1)=F(4,3)=F(5,4)=1; + F(3,2)=ca + sf/cf*sa; + + //TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F)); + + static TMatrixD Ft(5,6); + Ft(0,0)=sa; + Ft(0,1)=ca; + Ft(1,2)=Ft(3,4)=Ft(4,5)=1; + Ft(2,3)=ca + sf/cf*sa; + + TMatrixD tmp(C,TMatrixD::kMult,Ft); + T=new TMatrixD(F,TMatrixD::kMult,tmp); + + } // **** translation ****************** @@ -540,27 +744,27 @@ Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) { fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3; fP2 += dx*fP4; - TMatrixD mF(5,6); - mF(0,1)=mF(1,2)=mF(2,3)=mF(3,4)=mF(4,5)=1; - mF(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2)); - mF(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2)); - mF(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1)); - mF(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1); - mF(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1)); - mF(2,5)=dx; - mF(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2)); - mF(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1))); - mF(2,0)=-fP4; - - TMatrixD tmp(*tT,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, mF)); - delete tT; - TMatrixD cC(mF,TMatrixD::kMult,tmp); - - fC00=cC(0,0); - fC10=cC(1,0); fC11=cC(1,1); - fC20=cC(2,0); fC21=cC(2,1); fC22=cC(2,2); - fC30=cC(3,0); fC31=cC(3,1); fC32=cC(3,2); fC33=cC(3,3); - fC40=cC(4,0); fC41=cC(4,1); fC42=cC(4,2); fC43=cC(4,3); fC44=cC(4,4); + static TMatrixD F(5,6); + F(0,1)=F(1,2)=F(2,3)=F(3,4)=F(4,5)=1; + F(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2)); + F(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2)); + F(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1)); + F(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1); + F(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1)); + F(2,5)=dx; + F(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2)); + F(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1))); + F(2,0)=-fP4; + + TMatrixD tmp(*T,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F)); + delete T; + TMatrixD C(F,TMatrixD::kMult,tmp); + + fC00=C(0,0); + fC10=C(1,0); fC11=C(1,1); + fC20=C(2,0); fC21=C(2,1); fC22=C(2,2); + fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3); + fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4); if (!Invariant()) { fAlpha=alpha; @@ -578,6 +782,30 @@ Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) { return 1; } + + +Int_t AliITStrackV2::GetProlongationFast(Double_t alp, Double_t xk,Double_t &y, Double_t &z) +{ + //----------------------------------------------------------------------------- + //get fast prolongation + //----------------------------------------------------------------------------- + Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha); + Double_t cf=TMath::Sqrt(1.- fP2*fP2); + // **** rotation ********************** + y= -fX*sa + fP0*ca; + // **** translation ****************** + Double_t dx = xk- fX*ca - fP0*sa; + Double_t f1=fP2*ca - cf*sa, f2=f1 + fP4*dx; + if (TMath::Abs(f2) >= 0.9999) { + return 0; + } + Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2); + y += dx*(f1+f2)/(r1+r2); + z = fP1+dx*(f1+f2)/(f1*r2 + f2*r1)*fP3; + return 1; +} + + Double_t AliITStrackV2::GetD(Double_t x, Double_t y) const { //------------------------------------------------------------------ // This function calculates the transverse impact parameter @@ -606,7 +834,6 @@ Double_t AliITStrackV2::GetZat(Double_t x) const { if (TMath::Abs(f2) >= 0.9999) { return 10000000; } - Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2); Double_t z = fP1 + dx*(f1+f2)/(f1*r2 + f2*r1)*fP3; return z; diff --git a/ITS/AliITStrackV2.h b/ITS/AliITStrackV2.h index a6ca38bd079..dd5615d5c1a 100644 --- a/ITS/AliITStrackV2.h +++ b/ITS/AliITStrackV2.h @@ -30,18 +30,24 @@ #include "AliITSrecoV2.h" class AliESDtrack; +class AliTPCtrack; //_____________________________________________________________________________ class AliITStrackV2 : public AliKalmanTrack { + friend class AliITStrackerV2; + friend class AliITStrackerMI; public: AliITStrackV2(); + AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *); AliITStrackV2(AliESDtrack& t,Bool_t c=kFALSE) throw (const Char_t *); AliITStrackV2(const AliITStrackV2& t); Int_t PropagateToVertex(Double_t d=0., Double_t x0=0.); Int_t Propagate(Double_t alpha, Double_t xr); + Int_t GetProlongationFast(Double_t alpha, Double_t xr,Double_t &y, Double_t &z); Int_t CorrectForMaterial(Double_t d, Double_t x0=21.82); Int_t PropagateTo(Double_t xr, Double_t d, Double_t x0=21.82); Int_t Update(const AliCluster* cl,Double_t chi2,UInt_t i); + Int_t UpdateMI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz, Double_t chi2,UInt_t i); Int_t Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]); void SetdEdx(Double_t dedx) {fdEdx=dedx;} void SetSampledEdx(Float_t q, Int_t i); @@ -51,15 +57,15 @@ public: void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); } void UpdateESDtrack(ULong_t flags); void SetConstrainedESDtrack(Double_t chi2); - void SetReconstructed(Bool_t sr=kTRUE){fReconstructed = sr;} + Bool_t GetReconstructed() const {return fReconstructed;} void SetChi2MIP(Int_t i,Float_t val){fChi2MIP[i]=val;} Float_t GetChi2MIP(Int_t i) const {return fChi2MIP[i];} void IncrementNSkipped(){fNSkipped++;} // increment by 1 the # of skipped cls - Int_t GetNSkipped() const {return fNSkipped;} + Float_t GetNSkipped() const {return fNSkipped;} void IncrementNUsed(){fNUsed++;} // increment by 1 the # of shared clusters - Int_t GetNUsed() const {return fNUsed;} + Float_t GetNUsed() const {return fNUsed;} AliESDtrack *GetESDtrack() const {return fESDtrack;} Double_t GetCov33() const {return fC33;} // cov. matrix el. 3,3 Double_t GetCov44() const {return fC44;} // cov. matrix el. 4,4 @@ -68,6 +74,8 @@ public: Float_t GetSigmaY(Int_t i) const {return fSigmaY[i];} Float_t GetSigmaZ(Int_t i) const {return fSigmaZ[i];} + + Int_t GetDetectorIndex() const {return GetLabel();} Double_t GetX() const {return fX;} Double_t GetAlpha()const {return fAlpha;} @@ -91,12 +99,9 @@ public: Int_t GetClusterIndex(Int_t i) const {return fIndex[i];} Int_t GetGlobalXYZat(Double_t r,Double_t &x,Double_t &y,Double_t &z) const; Double_t GetPredictedChi2(const AliCluster *cluster) const; + Double_t GetPredictedChi2MI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz) const; Int_t Invariant() const; - AliITStrackV2& operator=(const AliITStrackV2& /*t*/) { - Error("operator=","Assignment is not allowed !"); return *this; - } - protected: Double_t fX; // X-coordinate of this track (reference plane) Double_t fAlpha; // rotation angle @@ -114,18 +119,29 @@ protected: Double_t fC20, fC21, fC22; // of the Double_t fC30, fC31, fC32, fC33; // track Double_t fC40, fC41, fC42, fC43, fC44; // parameters - Int_t fNUsed; // number of shared clusters - Int_t fNSkipped; // number of skipped clusters + Float_t fNUsed; // number of shared clusters + Float_t fNSkipped; // number of skipped clusters + Float_t fNDeadZone; // number of clusters in dead zone + Float_t fDeadZoneProbability; // probability to cross dead zone Bool_t fReconstructed; // reconstructed - accepted flag - Float_t fChi2MIP[6]; // MIP chi squres + Float_t fChi2MIP[12]; // MIP chi squres UInt_t fIndex[kMaxLayer]; // indices of associated clusters Float_t fdEdxSample[4]; // array of dE/dx samples b.b. - Float_t fDy[6]; //dy in layer - Float_t fDz[6]; //dz in layer - Float_t fSigmaY[6]; //sigma y - Float_t fSigmaZ[6]; //sigma z + Float_t fDy[12]; //dy in layer + Float_t fDz[12]; //dz in layer + Float_t fSigmaY[12]; //sigma y + Float_t fSigmaZ[12]; //sigma z + Float_t fNy[6]; //expected size of cluster + Float_t fNz[6]; //expected size of cluster + Float_t fD[2]; //distance to the vertex + Float_t fNormQ[6]; // normalized Q + Float_t fExpQ; // expected Q + Float_t fNormChi2[6]; // normalized chi2 + Float_t fChi22; // chi22 + Float_t fdEdxMismatch; + Bool_t fConstrain; //indication of the vertex constrain + Int_t fClIndex[6]; //cluster Index AliESDtrack *fESDtrack; //! pointer to the connected ESD track - ClassDef(AliITStrackV2,3) //ITS reconstructed track }; diff --git a/ITS/AliITStrackerMI.cxx b/ITS/AliITStrackerMI.cxx new file mode 100644 index 00000000000..9eecd7f1c3b --- /dev/null +++ b/ITS/AliITStrackerMI.cxx @@ -0,0 +1,3322 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//------------------------------------------------------------------------- +// Implementation of the ITS tracker class +// It reads AliITSclusterV2 clusters and creates AliITStrackV2 tracks +// and fills with them the ESD +// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch +// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch +//------------------------------------------------------------------------- + +#include + +#include +#include +#include + +#include "AliITSgeom.h" +#include "AliITSRecPoint.h" +#include "AliTPCtrack.h" +#include "AliESD.h" +#include "AliITSclusterV2.h" +//#include "AliITStrackerV2.h" +#include "AliITStrackerMI.h" +#include "TMatrixD.h" +#include "AliHelix.h" +#include "AliV0vertex.h" + +ClassImp(AliITStrackerMI) +ClassImp(AliITSRecV0Info) + + + +AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[kMaxLayer]; // ITS layers + +AliITStrackerMI::AliITStrackerMI(const AliITSgeom *geom) : AliTracker() { + //-------------------------------------------------------------------- + //This is the AliITStrackerMI constructor + //-------------------------------------------------------------------- + AliITSgeom *g=(AliITSgeom*)geom; + + Float_t x,y,z; + Int_t i; + for (i=1; iGetNladders(i); + Int_t ndet=g->GetNdetectors(i); + + g->GetTrans(i,1,1,x,y,z); + Double_t r=TMath::Sqrt(x*x + y*y); + Double_t poff=TMath::ATan2(y,x); + Double_t zoff=z; + + g->GetTrans(i,1,2,x,y,z); + r += TMath::Sqrt(x*x + y*y); + g->GetTrans(i,2,1,x,y,z); + r += TMath::Sqrt(x*x + y*y); + g->GetTrans(i,2,2,x,y,z); + r += TMath::Sqrt(x*x + y*y); + r*=0.25; + + new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet); + + for (Int_t j=1; jGetTrans(i,j,k,x,y,zshift); + Double_t rot[9]; g->GetRotMatrix(i,j,k,rot); + + Double_t phi=TMath::ATan2(rot[1],rot[0])+TMath::Pi(); + phi+=TMath::Pi()/2; + if (i==1) phi+=TMath::Pi(); + Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi); + Double_t r=x*cp+y*sp; + + AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1); + new(&det) AliITSdetector(r,phi); + } + } + + } + + fI=kMaxLayer; + + fPass=0; + fConstraint[0]=1; fConstraint[1]=0; + + Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV}; + SetVertex(xyz,ers); + + for (Int_t i=0; iGetBranch("Clusters"); + if (!branch) { + Error("LoadClusters"," can't get the branch !\n"); + return 1; + } + + TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy; + branch->SetAddress(&clusters); + + Int_t j=0; + Int_t detector=0; + for (Int_t i=0; iGetEvent(j)) continue; + Int_t ncl=clusters->GetEntriesFast(); + SignDeltas(clusters,GetZ()); + while (ncl--) { + AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl); + detector = c->GetDetectorIndex(); + fgLayers[i].InsertCluster(new AliITSclusterV2(*c)); + } + clusters->Delete(); + //add dead zone virtual "cluster" + if (i<2){ + for (Float_t ydead = 0; ydead < 1.31 ; ydead+=(i+1.)*0.018){ + Int_t lab[4] = {0,0,0,detector}; + Int_t info[3] = {0,0,0}; + Float_t hit[5]={0,0,0.004/12.,0.001/12.,0}; + if (i==0) hit[0] =ydead-0.4; + if (i==1) hit[0]=ydead-3.75; + hit[1] =-0.04; + if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.) + fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); + hit[1]=-7.05; + if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.) + fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); + hit[1]=-7.15; + if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.) + fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); + hit[1] =0.06; + if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.) + fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); + hit[1]=7.05; + if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.) + fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); + hit[1]=7.25; + if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.) + fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info)); + } + } + + } + // + fgLayers[i].ResetRoad(); //road defined by the cluster density + fgLayers[i].SortClusters(); + } + + return 0; +} + +void AliITStrackerMI::UnloadClusters() { + //-------------------------------------------------------------------- + //This function unloads ITS clusters + //-------------------------------------------------------------------- + for (Int_t i=0; iGetX() > riw) { + if (!t->PropagateTo(riw,diw,x0iw)) return 1; + if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr); + if (TMath::Abs(t->GetZ())CorrectForMaterial(dm); + if (!t->PropagateTo(rcd,dcd,x0cd)) return 1; + //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z); + //if (TMath::Abs(y)PropagateTo(rr,dr,x0r); + if (!t->PropagateTo(rs,ds)) return 1; + } else if (t->GetX() < rs) { + if (!t->PropagateTo(rs,-ds)) return 1; + //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z); + //if (TMath::Abs(y)PropagateTo(rr,-dr,x0r); + if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1; + if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1; + } else { + ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !"); + return 1; + } + + return 0; +} + +Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) { + //-------------------------------------------------------------------- + // This functions reconstructs ITS tracks + // The clusters must be already loaded ! + //-------------------------------------------------------------------- + TObjArray itsTracks(15000); + + {/* Read ESD tracks */ + Int_t nentr=event->GetNumberOfTracks(); + Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr); + while (nentr--) { + AliESDtrack *esd=event->GetTrack(nentr); + + if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue; + if (esd->GetStatus()&AliESDtrack::kTPCout) continue; + if (esd->GetStatus()&AliESDtrack::kITSin) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("Clusters2Tracks",msg); + delete t; + continue; + } + t->fD[0] = t->GetD(GetX(),GetY()); + t->fD[1] = t->GetZat(GetX())-GetZ(); + Double_t vdist = TMath::Sqrt(t->fD[0]*t->fD[0]+t->fD[1]*t->fD[1]); + if (t->GetMass()<0.13) t->SetMass(0.13957); // MI look to the esd - mass hypothesys !!!!!!!!!!! + // write expected q + t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.); + + if (TMath::Abs(t->fD[0])>10) { + delete t; + continue; + } + + if (TMath::Abs(vdist)>20) { + delete t; + continue; + } + if (TMath::Abs(1/t->Get1Pt())<0.120) { + delete t; + continue; + } + + if (CorrectForDeadZoneMaterial(t)!=0) { + Warning("Clusters2Tracks", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + t->fReconstructed = kFALSE; + itsTracks.AddLast(t); + } + } /* End Read ESD tracks */ + + itsTracks.Sort(); + Int_t nentr=itsTracks.GetEntriesFast(); + fTrackHypothesys.Expand(nentr); + MakeCoeficients(nentr); + Int_t ntrk=0; + for (fPass=0; fPass<2; fPass++) { + Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue; + for (Int_t i=0; ifReconstructed&&(t->fNUsed<1.5)) continue; //this track was already "succesfully" reconstructed + if ( (TMath::Abs(t->GetD(GetX(),GetY())) >3.) && fConstraint[fPass]) continue; + if ( (TMath::Abs(t->GetZat(GetX())-GetZ())>3.) && fConstraint[fPass]) continue; + + Int_t tpcLabel=t->GetLabel(); //save the TPC track label + fI = 6; + ResetTrackToFollow(*t); + ResetBestTrack(); + + FollowProlongationTree(t,i); + + + SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change + // + AliITStrackV2 * besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15); + if (!besttrack) continue; + besttrack->SetLabel(tpcLabel); + // besttrack->CookdEdx(); + CookdEdx(besttrack); + besttrack->fFakeRatio=1.; + CookLabel(besttrack,0.); //For comparison only + // besttrack->UpdateESDtrack(AliESDtrack::kITSin); + // + + UpdateESDtrack(besttrack,AliESDtrack::kITSin); + + if ( besttrack->GetNumberOfClusters()<5 && fConstraint[fPass]) { + continue; + } + if (besttrack->fChi2MIP[0]+besttrack->fNUsed>3.5) continue; + if ( (TMath::Abs(besttrack->fD[0]*besttrack->fD[0]+besttrack->fD[1]*besttrack->fD[1])>0.1) && fConstraint[fPass]) continue; + //delete itsTracks.RemoveAt(i); + t->fReconstructed = kTRUE; + ntrk++; + } + GetBestHypothesysMIP(itsTracks); + } + + //GetBestHypothesysMIP(itsTracks); + //FindV0(event); + + itsTracks.Delete(); + // + Int_t entries = fTrackHypothesys.GetEntriesFast(); + for (Int_t ientry=0;ientryDelete(); + delete fTrackHypothesys.RemoveAt(ientry); + } + + fTrackHypothesys.Delete(); + Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk); + + return 0; +} + + + +Int_t AliITStrackerMI::Clusters2Tracks(TTree *tpcTree, TTree *itsTree) { + //-------------------------------------------------------------------- + // This functions reconstructs ITS tracks + // The clusters must be already loaded ! + //-------------------------------------------------------------------- + Int_t nentr=0; TObjArray itsTracks(15000); + + Warning("Clusters2Tracks(TTree *, TTree *)", + "Will be removed soon ! Use Clusters2Tracks(AliESD *) instead."); + + {/* Read TPC tracks */ + AliTPCtrack *itrack=new AliTPCtrack; + TBranch *branch=tpcTree->GetBranch("tracks"); + if (!branch) { + Error("Clusters2Tracks","Can't get the branch !"); + return 1; + } + tpcTree->SetBranchAddress("tracks",&itrack); + nentr=(Int_t)tpcTree->GetEntries(); + + Info("Clusters2Tracks","Number of TPC tracks: %d\n",nentr); + + for (Int_t i=0; iGetEvent(i); + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*itrack); + } catch (const Char_t *msg) { + Warning("Clusters2Tracks",msg); + delete t; + continue; + } + if (TMath::Abs(t->GetD())>4) { + delete t; + continue; + } + + if (CorrectForDeadZoneMaterial(t)!=0) { + Warning("Clusters2Tracks", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + + itsTracks.AddLast(t); + } + delete itrack; + } + itsTracks.Sort(); + nentr=itsTracks.GetEntriesFast(); + + + AliITStrackV2 *otrack=&fBestTrack; + TBranch *branch=itsTree->GetBranch("tracks"); + if (!branch) itsTree->Branch("tracks","AliITStrackV2",&otrack,32000,3); + else branch->SetAddress(&otrack); + + for (fPass=0; fPass<2; fPass++) { + Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue; + for (Int_t i=0; iGetLabel(); //save the TPC track label + + ResetTrackToFollow(*t); + ResetBestTrack(); + /* + for (FollowProlongation(); fIFill(); + //UseClusters(&fBestTrack); + delete itsTracks.RemoveAt(i); + } + } + + nentr=(Int_t)itsTree->GetEntries(); + Info("Clusters2Tracks","Number of prolonged tracks: %d\n",nentr); + + itsTracks.Delete(); + + return 0; +} + +Int_t AliITStrackerMI::PropagateBack(AliESD *event) { + //-------------------------------------------------------------------- + // This functions propagates reconstructed ITS tracks back + // The clusters must be loaded ! + //-------------------------------------------------------------------- + Int_t nentr=event->GetNumberOfTracks(); + Info("PropagateBack", "Number of ESD tracks: %d\n", nentr); + + Int_t ntrk=0; + for (Int_t i=0; iGetTrack(i); + + if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue; + if (esd->GetStatus()&AliESDtrack::kITSout) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("PropagateBack",msg); + delete t; + continue; + } + t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.); + + ResetTrackToFollow(*t); + + // propagete to vertex [SR, GSI 17.02.2003] + // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov + if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) { + if (fTrackToFollow.PropagateToVertex()) { + fTrackToFollow.StartTimeIntegral(); + } + fTrackToFollow.PropagateTo(3.,-0.0028,65.19); + } + + fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters(); + if (RefitAt(49.,&fTrackToFollow,t)) { + if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) { + Warning("PropagateBack", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + fTrackToFollow.SetLabel(t->GetLabel()); + //fTrackToFollow.CookdEdx(); + CookLabel(&fTrackToFollow,0.); //For comparison only + fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout); + //UseClusters(&fTrackToFollow); + ntrk++; + } + delete t; + } + + Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk); + + return 0; +} + +Int_t AliITStrackerMI::RefitInward(AliESD *event) { + //-------------------------------------------------------------------- + // This functions refits ITS tracks using the + // "inward propagated" TPC tracks + // The clusters must be loaded ! + //-------------------------------------------------------------------- + Int_t nentr=event->GetNumberOfTracks(); + Info("RefitInward", "Number of ESD tracks: %d\n", nentr); + + Int_t ntrk=0; + for (Int_t i=0; iGetTrack(i); + + if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue; + if (esd->GetStatus()&AliESDtrack::kITSrefit) continue; + if (esd->GetStatus()&AliESDtrack::kTPCout) + if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue; + + AliITStrackV2 *t=0; + try { + t=new AliITStrackV2(*esd); + } catch (const Char_t *msg) { + Warning("RefitInward",msg); + delete t; + continue; + } + t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.); + if (CorrectForDeadZoneMaterial(t)!=0) { + Warning("RefitInward", + "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } + + ResetTrackToFollow(*t); + fTrackToFollow.ResetClusters(); + + if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) + fTrackToFollow.ResetCovariance(); + + //Refitting... + if (RefitAt(3.7, &fTrackToFollow, t)) { + fTrackToFollow.SetLabel(t->GetLabel()); + // fTrackToFollow.CookdEdx(); + CookdEdx(&fTrackToFollow); + + CookLabel(&fTrackToFollow,0.0); //For comparison only + + if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe + Double_t a=fTrackToFollow.GetAlpha(); + Double_t cs=TMath::Cos(a),sn=TMath::Sin(a); + Double_t xv= GetX()*cs + GetY()*sn; + Double_t yv=-GetX()*sn + GetY()*cs; + + Double_t c=fTrackToFollow.GetC(), snp=fTrackToFollow.GetSnp(); + Double_t x=fTrackToFollow.GetX(), y=fTrackToFollow.GetY(); + Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp)); + Double_t fv=TMath::ATan(tgfv); + + cs=TMath::Cos(fv); sn=TMath::Sin(fv); + x = xv*cs + yv*sn; + yv=-xv*sn + yv*cs; xv=x; + + if (fTrackToFollow.Propagate(fv+a,xv)) { + fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit); + //UseClusters(&fTrackToFollow); + { + AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ()); + c.SetSigmaY2(GetSigmaY()*GetSigmaY()); + c.SetSigmaZ2(GetSigmaZ()*GetSigmaZ()); + Double_t chi2=fTrackToFollow.GetPredictedChi2(&c); + //Double_t chi2=GetPredictedChi2MI(&fTrackToFollow,&c,fI); + if (chi2> 28; + Int_t c=(index & 0x0fffffff) >> 00; + return fgLayers[l].GetCluster(c); +} + + +void AliITStrackerMI::FollowProlongationTree(AliITStrackV2 * otrack, Int_t esdindex) +{ + //-------------------------------------------------------------------- + // Follow prolongation tree + //-------------------------------------------------------------------- + + //setup tree of the prolongations + // + static AliITStrackV2 tracks[7][100]; + AliITStrackV2 *currenttrack; + static AliITStrackV2 currenttrack1; + static AliITStrackV2 currenttrack2; + static AliITStrackV2 backuptrack; + Int_t ntracks[7]; + Int_t nindexes[7][100]; + Float_t normalizedchi2[100]; + for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0; + otrack->fNSkipped=0; + new (&(tracks[6][0])) AliITStrackV2(*otrack); + ntracks[6]=1; + nindexes[6][0]=0; + // + // + // follow prolongations + for (Int_t ilayer=5;ilayer>=0;ilayer--){ + // + AliITSlayer &layer=fgLayers[ilayer]; + Double_t r=layer.GetR(); + ntracks[ilayer]=0; + // + // + Int_t nskipped=0; + Float_t nused =0; + for (Int_t itrack =0;itrack=100) break; + if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNSkipped>0) nskipped++; + if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNUsed>2.) nused++; + if (ntracks[ilayer]>15+ilayer){ + if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNSkipped>0 && nskipped>4+ilayer) continue; + if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNUsed>2. && nused>3) continue; + } + + new(¤ttrack1) AliITStrackV2(tracks[ilayer+1][nindexes[ilayer+1][itrack]]); + if (ilayer==3 || ilayer==1) { + Double_t rs=0.5*(fgLayers[ilayer+1].GetR() + r); + Double_t d=0.0034, x0=38.6; + if (ilayer==1) {rs=9.; d=0.0097; x0=42;} + if (!currenttrack1.PropagateTo(rs,d,x0)) { + continue; + } + } + // + //find intersection with layer + Double_t x,y,z; + if (!currenttrack1.GetGlobalXYZat(r,x,y,z)) { + continue; + } + Double_t phi=TMath::ATan2(y,x); + Int_t idet=layer.FindDetectorIndex(phi,z); + if (idet<0) { + continue; + } + //propagate to the intersection + const AliITSdetector &det=layer.GetDetector(idet); + phi=det.GetPhi(); + new(¤ttrack2) AliITStrackV2(currenttrack1); + if (!currenttrack1.Propagate(phi,det.GetR())) { + continue; + } + currenttrack2.Propagate(phi,det.GetR()); // + currenttrack1.SetDetectorIndex(idet); + currenttrack2.SetDetectorIndex(idet); + + // + // + Double_t dz=7.5*TMath::Sqrt(currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]); + Double_t dy=7.5*TMath::Sqrt(currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]); + // + Bool_t isBoundary=kFALSE; + if (currenttrack1.GetY()-dy< det.GetYmin()+0.2) isBoundary = kTRUE; + if (currenttrack1.GetY()+dy> det.GetYmax()-0.2) isBoundary = kTRUE; + if (currenttrack1.GetZ()-dz< det.GetZmin()+0.2) isBoundary = kTRUE; + if (currenttrack1.GetZ()+dz> det.GetZmax()-0.2) isBoundary = kTRUE; + + if (isBoundary){ // track at boundary between detectors + Float_t maxtgl = TMath::Abs(currenttrack1.GetTgl()); + if (maxtgl>1) maxtgl=1; + dz = TMath::Sqrt(dz*dz+0.25*maxtgl*maxtgl); + // + Float_t maxsnp = TMath::Abs(currenttrack1.GetSnp()); + if (maxsnp>0.95) continue; + //if (maxsnp>0.5) maxsnp=0.5; + dy=TMath::Sqrt(dy*dy+0.25*maxsnp*maxsnp); + } + + Double_t zmin=currenttrack1.GetZ() - dz; + Double_t zmax=currenttrack1.GetZ() + dz; + Double_t ymin=currenttrack1.GetY() + r*phi - dy; + Double_t ymax=currenttrack1.GetY() + r*phi + dy; + layer.SelectClusters(zmin,zmax,ymin,ymax); + // + // loop over all possible prolongations + // + Double_t msz=1./((currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer])); + Double_t msy=1./((currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer])); + if (fConstraint[fPass]){ + msy/=60; msz/=60.; + } + else{ + msy/=50; msz/=50.; + } + // + const AliITSclusterV2 *c=0; Int_t ci=-1; + Double_t chi2=12345.; + Int_t deadzone=0; + currenttrack = ¤ttrack1; + while ((c=layer.GetNextCluster(ci))!=0) { + if (ntracks[ilayer]>95) break; //space for skipped clusters + Bool_t change =kFALSE; + if (c->GetQ()==0 && (deadzone==1)) continue; + Int_t idet=c->GetDetectorIndex(); + if (currenttrack->GetDetectorIndex()!=idet) { + const AliITSdetector &det=layer.GetDetector(idet); + Double_t y,z; + if (!currenttrack2.GetProlongationFast(det.GetPhi(),det.GetR(),y,z)) continue; + Float_t pz = (z - c->GetZ()) , py=(y - c->GetY()); + if (pz*pz*msz+py*py*msy>1.) continue; + // + new (&backuptrack) AliITStrackV2(currenttrack2); + change = kTRUE; + currenttrack =¤ttrack2; + if (!currenttrack->Propagate(det.GetPhi(),det.GetR())) { + new (currenttrack) AliITStrackV2(backuptrack); + change = kFALSE; + continue; + } + currenttrack->SetDetectorIndex(idet); + } + else{ + Float_t pz = (currenttrack->GetZ() - c->GetZ()) , py=(currenttrack->GetY() - c->GetY()); + if (pz*pz*msz+py*py*msy>1.) continue; + } + + chi2=GetPredictedChi2MI(currenttrack,c,ilayer); + if (chi2GetQ()==0) deadzone=1; // take dead zone only once + if (ntracks[ilayer]>=100) continue; + AliITStrackV2 * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackV2(*currenttrack); + updatetrack->fClIndex[ilayer]=0; + if (change){ + new (¤ttrack2) AliITStrackV2(backuptrack); + } + if (c->GetQ()!=0){ + if (!UpdateMI(updatetrack,c,chi2,(ilayer<<28)+ci)) continue; + updatetrack->SetSampledEdx(c->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b. + } + else { + updatetrack->fNDeadZone++; + updatetrack->fDeadZoneProbability=GetDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())); + } + if (c->IsUsed()){ + updatetrack->fNUsed++; + } + Double_t x0; + Double_t d=layer.GetThickness(updatetrack->GetY(),updatetrack->GetZ(),x0); + updatetrack->CorrectForMaterial(d,x0); + if (fConstraint[fPass]) { + updatetrack->fConstrain = fConstraint[fPass]; + fI = ilayer; + Double_t d=GetEffectiveThickness(0,0); //Think of this !!!! + Double_t xyz[]={GetX(),GetY(),GetZ()}; + Double_t ptfactor = 1; + Double_t ers[]={GetSigmaX()*ptfactor,GetSigmaY()*ptfactor,GetSigmaZ()}; + Bool_t isPrim = kTRUE; + if (ilayer<4){ + updatetrack->fD[0] = updatetrack->GetD(GetX(),GetY()); + updatetrack->fD[1] = updatetrack->GetZat(GetX())-GetZ(); + if ( TMath::Abs(updatetrack->fD[0]/(1.+ilayer))>0.4 || TMath::Abs(updatetrack->fD[1]/(1.+ilayer))>0.4) isPrim=kFALSE; + } + if (isPrim) updatetrack->Improve(d,xyz,ers); + } //apply vertex constrain + ntracks[ilayer]++; + } // create new hypothesy + } // loop over possible cluster prolongation + // if (fConstraint[fPass]&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0){ + if (itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){ + AliITStrackV2* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackV2(currenttrack1); + vtrack->fClIndex[ilayer]=0; + fI = ilayer; + Double_t d=GetEffectiveThickness(0,0); //Think of this !!!! + Double_t xyz[]={GetX(),GetY(),GetZ()}; + Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()}; + vtrack->Improve(d,xyz,ers); + vtrack->fNSkipped++; + ntracks[ilayer]++; + } + + } //loop over track candidates + // + // + Int_t accepted=0; + + Int_t golds=0; + for (Int_t itrack=0;itrack4) accepted++; + else{ + if ( fConstraint[fPass] && normalizedchi2[itrack]90) ntracks[ilayer]=90; + } //loop over layers + //printf("%d\t%d\t%d\t%d\t%d\t%d\n",ntracks[0],ntracks[1],ntracks[2],ntracks[3],ntracks[4],ntracks[5]); + + for (Int_t i=0;i7.)continue; + AddTrackHypothesys(new AliITStrackV2(track), esdindex); + } + for (Int_t i=0;i7.)continue; + if (fConstraint[fPass]) track.fNSkipped+=1; + if (!fConstraint[fPass]) { + track.fD[0] = track.GetD(GetX(),GetY()); + track.fNSkipped+=4./(4.+8.*abs(track.fD[0])); + if (track.fN+track.fNDeadZone+track.fNSkipped>6) { + track.fNSkipped = 6-track.fN+track.fNDeadZone; + } + } + AddTrackHypothesys(new AliITStrackV2(track), esdindex); + } + //} + + if (!fConstraint[fPass]){ + for (Int_t i=0;i7.)continue; + if (fConstraint[fPass]) track.fNSkipped+=2; + if (!fConstraint[fPass]){ + track.fD[0] = track.GetD(GetX(),GetY()); + track.fNSkipped+= 7./(7.+8.*abs(track.fD[0])); + if (track.fN+track.fNDeadZone+track.fNSkipped>6) { + track.fNSkipped = 6-track.fN+track.fNDeadZone; + } + } + AddTrackHypothesys(new AliITStrackV2(track), esdindex); + } + } +} + + +AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const +{ + //-------------------------------------------------------------------- + // + // + return fgLayers[layer]; +} +AliITStrackerMI::AliITSlayer::AliITSlayer() { + //-------------------------------------------------------------------- + //default AliITSlayer constructor + //-------------------------------------------------------------------- + fN=0; + fDetectors=0; + fSkip = 0; + fCurrentSlice=-1; + for (Int_t i=0; iIsUsed()) cl->Use(); + } + +} + +void AliITStrackerMI::AliITSlayer::ResetRoad() { + //-------------------------------------------------------------------- + // This function calculates the road defined by the cluster density + //-------------------------------------------------------------------- + Int_t n=0; + for (Int_t i=0; iGetZ())1) fRoad=2*fR*TMath::Sqrt(3.14/n); + if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n); +} + +Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSclusterV2 *c) { + //-------------------------------------------------------------------- + //This function adds a cluster to this layer + //-------------------------------------------------------------------- + if (fN==kMaxClusterPerLayer) { + ::Error("InsertCluster","Too many clusters !\n"); + return 1; + } + fCurrentSlice=-1; + if (fN==0) {fClusters[fN++]=c; return 0;} + Int_t i=FindClusterIndex(c->GetZ()); + memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*)); + memmove(fY+i+1 ,fY+i,(fN-i)*sizeof(Float_t)); + memmove(fZ+i+1 ,fZ+i,(fN-i)*sizeof(Float_t)); + fClusters[i]=c; fN++; + // + AliITSdetector &det=GetDetector(c->GetDetectorIndex()); + Double_t y=fR*det.GetPhi() + c->GetY(); + if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi(); + fY[i] = y; + fZ[i] = c->GetZ(); + if (c->GetY()GetY()); + if (c->GetY()>det.GetYmax()) det.SetYmax(c->GetY()); + if (c->GetZ()GetZ()); + if (c->GetZ()>det.GetZmax()) det.SetZmax(c->GetZ()); + + return 0; +} + +void AliITStrackerMI::AliITSlayer::SortClusters() +{ + // + //sort clusters + // + fYB[0]=10000000; + fYB[1]=-10000000; + for (Int_t i=0;ifYB[1]) fYB[1]=fY[i]; + fClusterIndex[i] = i; + } + // + // fill slices + fDy5 = (fYB[1]-fYB[0])/5.; + fDy10 = (fYB[1]-fYB[0])/10.; + fDy20 = (fYB[1]-fYB[0])/20.; + for (Int_t i=0;i<6;i++) fN5[i] =0; + for (Int_t i=0;i<11;i++) fN10[i]=0; + for (Int_t i=0;i<21;i++) fN20[i]=0; + // + for (Int_t i=0;i<6;i++) {fBy5[i][0] = fYB[0]+(i-0.75)*fDy5; fBy5[i][1] = fYB[0]+(i+0.75)*fDy5;} + for (Int_t i=0;i<11;i++) {fBy10[i][0] = fYB[0]+(i-0.75)*fDy10; fBy10[i][1] = fYB[0]+(i+0.75)*fDy10;} + for (Int_t i=0;i<21;i++) {fBy20[i][0] = fYB[0]+(i-0.75)*fDy20; fBy20[i][1] = fYB[0]+(i+0.75)*fDy20;} + // + // + for (Int_t i=0;ifYB[1]+fDy5) continue; + // + // slice 10 + Float_t fslice = TMath::Nint(10.*(curY-fYB[0])/(fYB[1]-fYB[0])); + Float_t ymiddle = fYB[0]+fslice*fDy10; + for (Int_t di =-1;di<=1;di++){ + if (TMath::Abs(curY-(ymiddle+(float)di*fDy10))<0.75*fDy10){ + // + Int_t slice = int(fslice+21.0001)-21+di; + if (slice<0) continue; + if (slice>10) continue; + fClusters10[slice][fN10[slice]] = fClusters[i]; + fY10[slice][fN10[slice]] = curY; + fZ10[slice][fN10[slice]] = fZ[i]; + fClusterIndex10[slice][fN10[slice]]=i; + fN10[slice]++; + } + } + // + // slice 5 + fslice = TMath::Nint(5.*(curY-fYB[0])/(fYB[1]-fYB[0])); + ymiddle = fYB[0]+fslice*fDy5; + for (Int_t di =-1;di<=1;di++){ + if (TMath::Abs(curY-(ymiddle+(float)di*fDy5))<0.75*fDy5){ + // + Int_t slice = int(fslice+21.0001)-21+di; + if (slice<0) continue; + if (slice>5) continue; + fClusters5[slice][fN5[slice]] = fClusters[i]; + fY5[slice][fN5[slice]] = curY; + fZ5[slice][fN5[slice]] = fZ[i]; + fClusterIndex5[slice][fN5[slice]]=i; + fN5[slice]++; + } + } + // + // slice 20 + fslice = TMath::Nint(20.*(curY-fYB[0])/(fYB[1]-fYB[0])); + ymiddle = fYB[0]+fslice*fDy20; + for (Int_t di =-1;di<=1;di++){ + if (TMath::Abs(curY-(ymiddle+(float)di*fDy20))<0.75*fDy20){ + // + Int_t slice = int(fslice+21.0001)-21+di; + if (slice<0) continue; + if (slice>20) continue; + fClusters20[slice][fN20[slice]] = fClusters[i]; + fY20[slice][fN20[slice]] = curY; + fZ20[slice][fN20[slice]] = fZ[i]; + fClusterIndex20[slice][fN20[slice]]=i; + fN20[slice]++; + } + } + } +} + + + + +Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const { + //-------------------------------------------------------------------- + // This function returns the index of the nearest cluster + //-------------------------------------------------------------------- + Int_t ncl=0; + const Float_t *zcl; + if (fCurrentSlice<0) { + ncl = fN; + zcl = fZ; + } + else{ + ncl = fNcs; + zcl = fZcs;; + } + + if (ncl==0) return 0; + Int_t b=0, e=ncl-1, m=(b+e)/2; + for (; b fClusters[m]->GetZ()) b=m+1; + if (z > zcl[m]) b=m+1; + else e=m; + } + return m; +} +/* +void AliITStrackerMI::AliITSlayer:: +SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) { + //-------------------------------------------------------------------- + // This function sets the "window" + //-------------------------------------------------------------------- + fI=FindClusterIndex(zmin); fZmax=zmax; + fImax = TMath::Min(FindClusterIndex(zmax)+1,fN); + Double_t circle=2*TMath::Pi()*fR; + if (ymax>circle) { ymax-=circle; ymin-=circle; } + fYmin=ymin; fYmax=ymax; + fSkip = 0; +} +*/ + + +void AliITStrackerMI::AliITSlayer:: +SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) { + //-------------------------------------------------------------------- + // This function sets the "window" + //-------------------------------------------------------------------- + + Double_t circle=2*TMath::Pi()*fR; + fYmin = ymin; fYmax =ymax; + Float_t ymiddle = (fYmin+fYmax)*0.5; + if (ymiddlefYB[1]) {fYmin-=circle; fYmax-=circle;ymiddle-=circle;} + } + // + fCurrentSlice =-1; + // defualt take all + fClustersCs = fClusters; + fClusterIndexCs = fClusterIndex; + fYcs = fY; + fZcs = fZ; + fNcs = fN; + // + //is in 20 slice? + if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){ + Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20); + if (slice<0) slice=0; + if (slice>20) slice=20; + Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax10) slice=10; + Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax5) slice=5; + Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmaxGetZ() > fZmax) break; + // if (c->IsUsed()) continue; + const AliITSdetector &det=GetDetector(c->GetDetectorIndex()); + Double_t y=fR*det.GetPhi() + c->GetY(); + + if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi(); + if (y>1.*fR*TMath::Pi() && fYmaxfYmax) continue; + cluster=c; ci=i; + fI=i+1; + break; + } + return cluster; +} +*/ + +const AliITSclusterV2 *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci){ + //-------------------------------------------------------------------- + // This function returns clusters within the "window" + //-------------------------------------------------------------------- + + if (fCurrentSlice<0){ + Double_t rpi2 = 2.*fR*TMath::Pi(); + for (Int_t i=fI; ifYmax) continue; + if (fClusters[i]->GetQ()==0&&fSkip==2) continue; + ci=i; + fI=i+1; + return fClusters[i]; + } + } + else{ + for (Int_t i=fI; ifYmax) continue; + if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue; + ci=fClusterIndexCs[i]; + fI=i+1; + return fClustersCs[i]; + } + } + return 0; +} + + + +Int_t AliITStrackerMI::AliITSlayer:: +FindDetectorIndex(Double_t phi, Double_t z) const { + //-------------------------------------------------------------------- + //This function finds the detector crossed by the track + //-------------------------------------------------------------------- + Double_t dphi=-(phi-fPhiOffset); + if (dphi < 0) dphi += 2*TMath::Pi(); + else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi(); + Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5); + if (np>=fNladders) np-=fNladders; + if (np<0) np+=fNladders; + + Double_t dz=fZOffset-z; + Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5); + if (nz>=fNdetectors) return -1; + if (nz<0) return -1; + + return np*fNdetectors + nz; +} + +Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0) +const { + //-------------------------------------------------------------------- + //This function returns the layer thickness at this point (units X0) + //-------------------------------------------------------------------- + Double_t d=0.0085; + x0=21.82; + if (433.40) d+=dd; + if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);} + if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);} + for (Int_t i=0; i<12; i++) { + if (TMath::Abs(z-3.9*(i+0.5))<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=0.0034; + break; + } + if (TMath::Abs(z+3.9*(i+0.5))<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=0.0034; + break; + } + if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;} + if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;} + } + } else + if (373.40) d+=dd; + if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);} + if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);} + for (Int_t i=0; i<11; i++) { + if (TMath::Abs(z-3.9*i)<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=dd; + break; + } + if (TMath::Abs(z+3.9*i)<0.15) { + if (TMath::Abs(y-0.00)>3.40) d+=dd; + d+=dd; + break; + } + if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;} + if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;} + } + } else + if (133.30) d+=dd; + + if (TMath::Abs(y-1.80)<0.55) { + d+=0.016; + for (Int_t j=0; j<20; j++) { + if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;} + if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;} + } + } + if (TMath::Abs(y+1.80)<0.55) { + d+=0.016; + for (Int_t j=0; j<20; j++) { + if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;} + if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;} + } + } + + for (Int_t i=0; i<4; i++) { + if (TMath::Abs(z-7.3*i)<0.60) { + d+=dd; + if (TMath::Abs(y-0.00)>3.30) d+=dd; + break; + } + if (TMath::Abs(z+7.3*i)<0.60) { + d+=dd; + if (TMath::Abs(y-0.00)>3.30) d+=dd; + break; + } + } + } else + if (60.5) d+=dd; + //if (TMath::Abs(y-3.08)>0.45) d+=dd; + if (TMath::Abs(y-3.03)<0.10) {d+=0.014;} + } else + if (30.6) d+=dd; + //if (TMath::Abs(y+0.21)>0.45) d+=dd; + if (TMath::Abs(y+0.10)<0.10) {d+=0.014;} + } + + return d; +} + +Double_t AliITStrackerMI::GetEffectiveThickness(Double_t y,Double_t z) const +{ + //-------------------------------------------------------------------- + //Returns the thickness between the current layer and the vertex (units X0) + //-------------------------------------------------------------------- + Double_t d=0.0028*3*3; //beam pipe + Double_t x0=0; + + Double_t xn=fgLayers[fI].GetR(); + for (Int_t i=0; i1) { + Double_t xi=9.; + d+=0.0097*xi*xi; + } + + if (fI>3) { + Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR()); + d+=0.0034*xi*xi; + } + + return d/(xn*xn); +} + +Int_t AliITStrackerMI::AliITSlayer::InRoad() const { + //-------------------------------------------------------------------- + // This function returns number of clusters within the "window" + //-------------------------------------------------------------------- + Int_t ncl=0; + for (Int_t i=fI; iGetZ() > fZmax) break; + if (c->IsUsed()) continue; + const AliITSdetector &det=GetDetector(c->GetDetectorIndex()); + Double_t y=fR*det.GetPhi() + c->GetY(); + + if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi(); + if (y>1.*fR*TMath::Pi() && fYmaxfYmax) continue; + ncl++; + } + return ncl; +} + +Bool_t +AliITStrackerMI::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) { + //-------------------------------------------------------------------- + // This function refits the track "t" at the position "x" using + // the clusters from "c" + //-------------------------------------------------------------------- + Int_t index[kMaxLayer]; + Int_t k; + for (k=0; kGetNumberOfClusters(); + for (k=0; kGetClusterIndex(k),nl=(idx&0xf0000000)>>28; + index[nl]=idx; + } + + Int_t from, to, step; + if (xx > t->GetX()) { + from=0; to=kMaxLayer; + step=+1; + } else { + from=kMaxLayer-1; to=-1; + step=-1; + } + + for (Int_t i=from; i != to; i += step) { + AliITSlayer &layer=fgLayers[i]; + Double_t r=layer.GetR(); + + { + Double_t hI=i-0.5*step; + if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) { + Double_t rs=0.5*(fgLayers[i-step].GetR() + r); + Double_t d=0.0034, x0=38.6; + if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;} + if (!t->PropagateTo(rs,-step*d,x0)) { + return kFALSE; + } + } + } + + // remember old position [SR, GSI 18.02.2003] + Double_t oldX=0., oldY=0., oldZ=0.; + if (t->IsStartedTimeIntegral() && step==1) { + t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ); + } + // + + Double_t x,y,z; + if (!t->GetGlobalXYZat(r,x,y,z)) { + return kFALSE; + } + Double_t phi=TMath::ATan2(y,x); + Int_t idet=layer.FindDetectorIndex(phi,z); + if (idet<0) { + return kFALSE; + } + const AliITSdetector &det=layer.GetDetector(idet); + phi=det.GetPhi(); + if (!t->Propagate(phi,det.GetR())) { + return kFALSE; + } + t->SetDetectorIndex(idet); + + const AliITSclusterV2 *cl=0; + Double_t maxchi2=1000.*kMaxChi2; + + Int_t idx=index[i]; + if (idx>0) { + const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx); + if (idet != c->GetDetectorIndex()) { + idet=c->GetDetectorIndex(); + const AliITSdetector &det=layer.GetDetector(idet); + if (!t->Propagate(det.GetPhi(),det.GetR())) { + return kFALSE; + } + t->SetDetectorIndex(idet); + } + //Double_t chi2=t->GetPredictedChi2(c); + Int_t layer = (idx & 0xf0000000) >> 28;; + Double_t chi2=GetPredictedChi2MI(t,c,layer); + if (chi2GetNumberOfClusters()>2) { + Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]); + Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]); + Double_t zmin=t->GetZ() - dz; + Double_t zmax=t->GetZ() + dz; + Double_t ymin=t->GetY() + phi*r - dy; + Double_t ymax=t->GetY() + phi*r + dy; + layer.SelectClusters(zmin,zmax,ymin,ymax); + + const AliITSclusterV2 *c=0; Int_t ci=-1; + while ((c=layer.GetNextCluster(ci))!=0) { + if (idet != c->GetDetectorIndex()) continue; + Double_t chi2=t->GetPredictedChi2(c); + if (chi2Update(cl,maxchi2,idx)) { + if (!UpdateMI(t,cl,maxchi2,idx)) { + return kFALSE; + } + t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1); + } + + { + Double_t x0; + Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0); + t->CorrectForMaterial(-step*d,x0); + } + + // track time update [SR, GSI 17.02.2003] + if (t->IsStartedTimeIntegral() && step==1) { + Double_t newX, newY, newZ; + t->GetGlobalXYZat(t->GetX(),newX,newY,newZ); + Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) + + (oldZ-newZ)*(oldZ-newZ); + t->AddTimeStep(TMath::Sqrt(dL2)); + } + // + + } + + if (!t->PropagateTo(xx,0.,0.)) return kFALSE; + return kTRUE; +} + + +Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackV2 * track, Int_t mode) +{ + // + // calculate normalized chi2 + // return NormalizedChi2(track,0); + Float_t chi2 = 0; + Float_t sum=0; + Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack); + // track->fdEdxMismatch=0; + Float_t dedxmismatch =0; + Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack); + if (mode<100){ + for (Int_t i = 0;i<6;i++){ + if (track->fClIndex[i]>0){ + Float_t cerry, cerrz; + if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];} + else + { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];} + cerry*=cerry; + cerrz*=cerrz; + Float_t cchi2 = (track->fDy[i]*track->fDy[i]/cerry)+(track->fDz[i]*track->fDz[i]/cerrz); + if (i>1){ + Float_t ratio = track->fNormQ[i]/track->fExpQ; + if (ratio<0.5) { + cchi2+=(0.5-ratio)*10.; + //track->fdEdxMismatch+=(0.5-ratio)*10.; + dedxmismatch+=(0.5-ratio)*10.; + } + } + if (i<2 ||i>3){ + AliITSclusterV2 * cl = (AliITSclusterV2*)GetCluster( track->fClIndex[i]); + Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i]; + if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.); + if (i<2) chi2+=2*cl->GetDeltaProbability(); + } + chi2+=cchi2; + sum++; + } + } + if (TMath::Abs(track->fdEdxMismatch-dedxmismatch)>0.0001){ + track->fdEdxMismatch = dedxmismatch; + } + } + else{ + for (Int_t i = 0;i<4;i++){ + if (track->fClIndex[i]>0){ + Float_t cerry, cerrz; + if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];} + else { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];} + cerry*=cerry; + cerrz*=cerrz; + chi2+= (track->fDy[i]*track->fDy[i]/cerry); + chi2+= (track->fDz[i]*track->fDz[i]/cerrz); + sum++; + } + } + for (Int_t i = 4;i<6;i++){ + if (track->fClIndex[i]>0){ + Float_t cerry, cerrz; + if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];} + else { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];} + cerry*=cerry; + cerrz*=cerrz; + Float_t cerryb, cerrzb; + if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];} + else { cerryb= track->fSigmaY[i+6]; cerrzb = track->fSigmaZ[i+6];} + cerryb*=cerryb; + cerrzb*=cerrzb; + chi2+= TMath::Min((track->fDy[i+6]*track->fDy[i+6]/cerryb),track->fDy[i]*track->fDy[i]/cerry); + chi2+= TMath::Min((track->fDz[i+6]*track->fDz[i+6]/cerrzb),track->fDz[i]*track->fDz[i]/cerrz); + sum++; + } + } + } + if (track->fESDtrack->GetTPCsignal()>85){ + Float_t ratio = track->fdEdx/track->fESDtrack->GetTPCsignal(); + if (ratio<0.5) { + chi2+=(0.5-ratio)*5.; + } + if (ratio>2){ + chi2+=(ratio-2.0)*3; + } + } + // + Double_t match = TMath::Sqrt(track->fChi22); + if (track->fConstrain) match/=track->GetNumberOfClusters(); + if (!track->fConstrain) match/=track->GetNumberOfClusters()-2.; + if (match<0) match=0; + Float_t deadzonefactor = (track->fNDeadZone>0) ? 3*(1.1-track->fDeadZoneProbability):0.; + Double_t normchi2 = 2*track->fNSkipped+match+deadzonefactor+(1+(2*track->fNSkipped+deadzonefactor)/track->GetNumberOfClusters())* + (chi2)/TMath::Max(double(sum-track->fNSkipped), + 1./(1.+track->fNSkipped)); + + return normchi2; +} + + +Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackV2 * track1, AliITStrackV2 * track2) +{ + // + // return matching chi2 between two tracks + AliITStrackV2 track3(*track2); + track3.Propagate(track1->GetAlpha(),track1->GetX()); + TMatrixD vec(5,1); + vec(0,0)=track1->fP0-track3.fP0; + vec(1,0)=track1->fP1-track3.fP1; + vec(2,0)=track1->fP2-track3.fP2; + vec(3,0)=track1->fP3-track3.fP3; + vec(4,0)=track1->fP4-track3.fP4; + // + TMatrixD cov(5,5); + cov(0,0) = track1->fC00+track3.fC00; + cov(1,1) = track1->fC11+track3.fC11; + cov(2,2) = track1->fC22+track3.fC22; + cov(3,3) = track1->fC33+track3.fC33; + cov(4,4) = track1->fC44+track3.fC44; + + cov(0,1)=cov(1,0) = track1->fC10+track3.fC10; + cov(0,2)=cov(2,0) = track1->fC20+track3.fC20; + cov(0,3)=cov(3,0) = track1->fC30+track3.fC30; + cov(0,4)=cov(4,0) = track1->fC40+track3.fC40; + // + cov(1,2)=cov(2,1) = track1->fC21+track3.fC21; + cov(1,3)=cov(3,1) = track1->fC31+track3.fC31; + cov(1,4)=cov(4,1) = track1->fC41+track3.fC41; + // + cov(2,3)=cov(3,2) = track1->fC32+track3.fC32; + cov(2,4)=cov(4,2) = track1->fC42+track3.fC42; + // + cov(3,4)=cov(4,3) = track1->fC43+track3.fC43; + + cov.Invert(); + TMatrixD vec2(cov,TMatrixD::kMult,vec); + TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec); + return chi2(0,0); +} + +Double_t AliITStrackerMI::GetDeadZoneProbability(Double_t zpos, Double_t zerr) +{ + // + // return probability that given point - characterized by z position and error is in dead zone + // + Double_t probability =0; + Double_t absz = TMath::Abs(zpos); + Double_t nearestz = (absz<2)? 0.:7.1; + if (TMath::Abs(absz-nearestz)>0.25+3*zerr) return 0; + Double_t zmin=0, zmax=0; + if (zpos<-6.){ + zmin = -7.25; zmax = -6.95; + } + if (zpos>6){ + zmin = 7.0; zmax =7.3; + } + if (absz<2){ + zmin = -0.75; zmax = 1.5; + } + probability = (TMath::Erf((zpos-zmin)/zerr) - TMath::Erf((zpos-zmax)/zerr))*0.5; + return probability; +} + + +Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackV2 * track, Float_t fac) +{ + // + // calculate normalized chi2 + Float_t chi2[6]; + Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack); + Float_t ncl = 0; + for (Int_t i = 0;i<6;i++){ + if (TMath::Abs(track->fDy[i])>0){ + chi2[i]= (track->fDy[i]/erry[i])*(track->fDy[i]/erry[i]); + chi2[i]+= (track->fDz[i]/errz[i])*(track->fDz[i]/errz[i]); + ncl++; + } + else{chi2[i]=10000;} + } + Int_t index[6]; + TMath::Sort(6,chi2,index,kFALSE); + Float_t max = float(ncl)*fac-1.; + Float_t sumchi=0, sumweight=0; + for (Int_t i=0;ifNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000; + Int_t npoints = 0; + Double_t res =0; + for (Int_t i=0;i<6;i++){ + if ( (backtrack->fSigmaY[i]<0.000000001) || (forwardtrack->fSigmaY[i]<0.000000001)) continue; + Double_t sy1 = forwardtrack->fSigmaY[i]; + Double_t sz1 = forwardtrack->fSigmaZ[i]; + Double_t sy2 = backtrack->fSigmaY[i]; + Double_t sz2 = backtrack->fSigmaZ[i]; + if (i<2){ sy2=1000.;sz2=1000;} + // + Double_t dy0 = (forwardtrack->fDy[i]/(sy1*sy1) +backtrack->fDy[i]/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2)); + Double_t dz0 = (forwardtrack->fDz[i]/(sz1*sz1) +backtrack->fDz[i]/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2)); + // + Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2))); + Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2))); + // + res+= nz0*nz0+ny0*ny0; + npoints++; + } + if (npoints>1) return + TMath::Max(TMath::Abs(0.3*forwardtrack->Get1Pt())-0.5,0.)+ + //2*forwardtrack->fNUsed+ + res/TMath::Max(double(npoints-forwardtrack->fNSkipped), + 1./(1.+forwardtrack->fNSkipped)); + return 1000; +} + + + + + +Float_t *AliITStrackerMI::GetWeight(Int_t index) { + //-------------------------------------------------------------------- + // Return pointer to a given cluster + //-------------------------------------------------------------------- + Int_t l=(index & 0xf0000000) >> 28; + Int_t c=(index & 0x0fffffff) >> 00; + return fgLayers[l].GetWeight(c); +} + +void AliITStrackerMI::RegisterClusterTracks(AliITStrackV2* track,Int_t id) +{ + //--------------------------------------------- + // register track to the list + for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ + Int_t index = track->GetClusterIndex(icluster); + Int_t l=(index & 0xf0000000) >> 28; + Int_t c=(index & 0x0fffffff) >> 00; + if (c>fgLayers[l].fN) continue; + for (Int_t itrack=0;itrack<4;itrack++){ + if (fgLayers[l].fClusterTracks[itrack][c]<0){ + fgLayers[l].fClusterTracks[itrack][c]=id; + break; + } + } + } +} +void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackV2* track, Int_t id) +{ + //--------------------------------------------- + // unregister track from the list + for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ + Int_t index = track->GetClusterIndex(icluster); + Int_t l=(index & 0xf0000000) >> 28; + Int_t c=(index & 0x0fffffff) >> 00; + if (c>fgLayers[l].fN) continue; + for (Int_t itrack=0;itrack<4;itrack++){ + if (fgLayers[l].fClusterTracks[itrack][c]==id){ + fgLayers[l].fClusterTracks[itrack][c]=-1; + } + } + } +} +Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackV2* track,Int_t id, Int_t list[6], AliITSclusterV2 *clist[6]) +{ + //------------------------------------------------------------- + //get number of shared clusters + //------------------------------------------------------------- + Float_t shared=0; + for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;} + // mean number of clusters + Float_t *ny = GetNy(id), *nz = GetNz(id); + + + for (Int_t icluster=0;iclusterGetNumberOfClusters();icluster++){ + Int_t index = track->GetClusterIndex(icluster); + Int_t l=(index & 0xf0000000) >> 28; + Int_t c=(index & 0x0fffffff) >> 00; + if (c>fgLayers[l].fN) continue; + if (ny[l]==0){ + printf("problem\n"); + } + AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(index); + Float_t weight=1; + // + Float_t deltan = 0; + if (l>3&&cl->GetNy()+cl->GetNz()>6) continue; + if (l>2&&track->fNormQ[l]/track->fExpQ>3.5) continue; + if (l<2 || l>3){ + deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]); + } + else{ + deltan = (cl->GetNz()-nz[l]); + } + if (deltan>2.0) continue; // extended - highly probable shared cluster + weight = 2./TMath::Max(3.+deltan,2.); + // + for (Int_t itrack=0;itrack<4;itrack++){ + if (fgLayers[l].fClusterTracks[itrack][c]>=0 && fgLayers[l].fClusterTracks[itrack][c]!=id){ + list[l]=index; + clist[l] = (AliITSclusterV2*)GetCluster(index); + shared+=weight; + break; + } + } + } + track->fNUsed=shared; + return shared; +} + +Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackV2 *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6]) +{ + // + // find first shared track + // + // mean number of clusters + Float_t *ny = GetNy(trackID), *nz = GetNz(trackID); + // + for (Int_t i=0;i<6;i++) overlist[i]=-1; + Int_t sharedtrack=100000; + Int_t tracks[24],trackindex=0; + for (Int_t i=0;i<24;i++) {tracks[i]=-1;} + // + for (Int_t icluster=0;icluster<6;icluster++){ + if (clusterlist[icluster]<0) continue; + Int_t index = clusterlist[icluster]; + Int_t l=(index & 0xf0000000) >> 28; + Int_t c=(index & 0x0fffffff) >> 00; + if (ny[l]==0){ + printf("problem\n"); + } + if (c>fgLayers[l].fN) continue; + //if (l>3) continue; + AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(index); + // + Float_t deltan = 0; + if (l>3&&cl->GetNy()+cl->GetNz()>6) continue; + if (l>2&&track->fNormQ[l]/track->fExpQ>3.5) continue; + if (l<2 || l>3){ + deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]); + } + else{ + deltan = (cl->GetNz()-nz[l]); + } + if (deltan>2.0) continue; // extended - highly probable shared cluster + // + for (Int_t itrack=3;itrack>=0;itrack--){ + if (fgLayers[l].fClusterTracks[itrack][c]<0) continue; + if (fgLayers[l].fClusterTracks[itrack][c]!=trackID){ + tracks[trackindex] = fgLayers[l].fClusterTracks[itrack][c]; + trackindex++; + } + } + } + if (trackindex==0) return -1; + if (trackindex==1){ + sharedtrack = tracks[0]; + }else{ + if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]); + else{ + // + Int_t track[24], cluster[24]; + for (Int_t i=0;imax) { + sharedtrack=track[index]; + max=cluster[index]; + } + } + } + } + + if (sharedtrack>=100000) return -1; + // + // make list of overlaps + shared =0; + for (Int_t icluster=0;icluster<6;icluster++){ + if (clusterlist[icluster]<0) continue; + Int_t index = clusterlist[icluster]; + Int_t l=(index & 0xf0000000) >> 28; + Int_t c=(index & 0x0fffffff) >> 00; + if (c>fgLayers[l].fN) continue; + AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(index); + if (l==0 || l==1){ + if (cl->GetNy()>2) continue; + if (cl->GetNz()>2) continue; + } + if (l==4 || l==5){ + if (cl->GetNy()>3) continue; + if (cl->GetNz()>3) continue; + } + // + for (Int_t itrack=3;itrack>=0;itrack--){ + if (fgLayers[l].fClusterTracks[itrack][c]<0) continue; + if (fgLayers[l].fClusterTracks[itrack][c]==sharedtrack){ + overlist[l]=index; + shared++; + } + } + } + return sharedtrack; +} + + +AliITStrackV2 * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){ + // + // try to find track hypothesys without conflicts + // with minimal chi2; + TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1); + Int_t entries1 = arr1->GetEntriesFast(); + TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2); + if (!arr2) return (AliITStrackV2*) arr1->UncheckedAt(0); + Int_t entries2 = arr2->GetEntriesFast(); + if (entries2<=0) return (AliITStrackV2*) arr1->UncheckedAt(0); + // + AliITStrackV2 * track10=(AliITStrackV2*) arr1->UncheckedAt(0); + AliITStrackV2 * track20=(AliITStrackV2*) arr2->UncheckedAt(0); + if (TMath::Abs(1./track10->Get1Pt())>0.5+TMath::Abs(1/track20->Get1Pt())) return track10; + + for (Int_t itrack=0;itrackUncheckedAt(itrack); + UnRegisterClusterTracks(track,trackID1); + } + // + for (Int_t itrack=0;itrackUncheckedAt(itrack); + UnRegisterClusterTracks(track,trackID2); + } + Int_t index1=0; + Int_t index2=0; + Float_t maxconflicts=6; + Double_t maxchi2 =1000.; + // + // get weight of hypothesys - using best hypothesy + Double_t w1,w2; + + Int_t list1[6],list2[6]; + AliITSclusterV2 *clist1[6], *clist2[6] ; + RegisterClusterTracks(track10,trackID1); + RegisterClusterTracks(track20,trackID2); + Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1); + Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2); + UnRegisterClusterTracks(track10,trackID1); + UnRegisterClusterTracks(track20,trackID2); + // + // normalized chi2 + Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0; + Float_t nerry[6],nerrz[6]; + Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1); + Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2); + for (Int_t i=0;i<6;i++){ + if ( (erry1[i]>0) && (erry2[i]>0)) { + nerry[i] = TMath::Min(erry1[i],erry2[i]); + nerrz[i] = TMath::Min(errz1[i],errz2[i]); + }else{ + nerry[i] = TMath::Max(erry1[i],erry2[i]); + nerrz[i] = TMath::Max(errz1[i],errz2[i]); + } + if (TMath::Abs(track10->fDy[i])>0.000000000000001){ + chi21 += track10->fDy[i]*track10->fDy[i]/(nerry[i]*nerry[i]); + chi21 += track10->fDz[i]*track10->fDz[i]/(nerrz[i]*nerrz[i]); + ncl1++; + } + if (TMath::Abs(track20->fDy[i])>0.000000000000001){ + chi22 += track20->fDy[i]*track20->fDy[i]/(nerry[i]*nerry[i]); + chi22 += track20->fDz[i]*track20->fDz[i]/(nerrz[i]*nerrz[i]); + ncl2++; + } + } + chi21/=ncl1; + chi22/=ncl2; + // + // + Float_t d1 = TMath::Sqrt(track10->fD[0]*track10->fD[0]+track10->fD[1]*track10->fD[1])+0.1; + Float_t d2 = TMath::Sqrt(track20->fD[0]*track20->fD[0]+track20->fD[1]*track20->fD[1])+0.1; + Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2()); + Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2()); + // + w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+ + +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.) + +1.*TMath::Abs(1./track10->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt())) + ); + w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+ + s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.) + +1.*TMath::Abs(1./track20->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt())) + ); + + Double_t sumw = w1+w2; + w1/=sumw; + w2/=sumw; + if (w1fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.; + //Float_t maxconflicts0 = w1*conflict1+w2*conflict2; + // + // get pair of "best" hypothesys + // + Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1); + Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2); + + for (Int_t itrack1=0;itrack1UncheckedAt(itrack1); + //if (track1->fFakeRatio>0) continue; + RegisterClusterTracks(track1,trackID1); + for (Int_t itrack2=0;itrack2UncheckedAt(itrack2); + + // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0]; + //if (track2->fFakeRatio>0) continue; + Float_t nskipped=0; + RegisterClusterTracks(track2,trackID2); + Int_t list1[6],list2[6]; + AliITSclusterV2 *clist1[6], *clist2[6] ; + Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1); + Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2); + UnRegisterClusterTracks(track2,trackID2); + // + if (track1->fConstrain) nskipped+=w1*track1->fNSkipped; + if (track2->fConstrain) nskipped+=w2*track2->fNSkipped; + if (nskipped>0.5) continue; + // + //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue; + if (conflict1+13){ + deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i])); + } + else{ + deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i])); + } + c1 = 2./TMath::Max(3.+deltan,2.); + c2 = 2./TMath::Max(3.+deltan,2.); + } + else{ + if (clist1[i]){ + Float_t deltan = 0; + if (i<2 || i>3){ + deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]); + } + else{ + deltan = (clist1[i]->GetNz()-nz1[i]); + } + c1 = 2./TMath::Max(3.+deltan,2.); + c2 = 0; + } + + if (clist2[i]){ + Float_t deltan = 0; + if (i<2 || i>3){ + deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]); + } + else{ + deltan = (clist2[i]->GetNz()-nz2[i]); + } + c2 = 2./TMath::Max(3.+deltan,2.); + c1 = 0; + } + } + // + Double_t chi21=0,chi22=0; + if (TMath::Abs(track1->fDy[i])>0.) { + chi21 = (track1->fDy[i]/track1->fSigmaY[i])*(track1->fDy[i]/track1->fSigmaY[i])+ + (track1->fDz[i]/track1->fSigmaZ[i])*(track1->fDz[i]/track1->fSigmaZ[i]); + //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+ + // (track1->fDz[i]*track1->fDz[i])/(nerrz[i]*nerrz[i]); + }else{ + if (TMath::Abs(track1->fSigmaY[i]>0.)) c1=1; + } + // + if (TMath::Abs(track2->fDy[i])>0.) { + chi22 = (track2->fDy[i]/track2->fSigmaY[i])*(track2->fDy[i]/track2->fSigmaY[i])+ + (track2->fDz[i]/track2->fSigmaZ[i])*(track2->fDz[i]/track2->fSigmaZ[i]); + //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+ + // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]); + } + else{ + if (TMath::Abs(track2->fSigmaY[i]>0.)) c2=1; + } + sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2); + if (chi21>0) sum+=w1; + if (chi22>0) sum+=w2; + conflict+=(c1+c2); + } + Double_t norm = sum-w1*track1->fNSkipped-w2*track2->fNSkipped; + if (norm<0) norm =1/(w1*track1->fNSkipped+w2*track2->fNSkipped); + Double_t normchi2 = 2*conflict+sumchi2/sum; + if ( normchi2 fFakeRatio*track10->GetNumberOfClusters(); + AliITStrackV2* track1=(AliITStrackV2*) arr1->UncheckedAt(index1); + track1->fChi2MIP[5] = maxconflicts; + track1->fChi2MIP[6] = maxchi2; + track1->fChi2MIP[7] = 0.01+orig-(track1->fFakeRatio*track1->GetNumberOfClusters()); + // track1->UpdateESDtrack(AliESDtrack::kITSin); + track1->fChi2MIP[8] = index1; + fBestTrackIndex[trackID1] =index1; + UpdateESDtrack(track1, AliESDtrack::kITSin); + } + else if (track10->fChi2MIP[0]fChi2MIP[5] = maxconflicts; + track10->fChi2MIP[6] = maxchi2; + // track10->UpdateESDtrack(AliESDtrack::kITSin); + UpdateESDtrack(track10,AliESDtrack::kITSin); + } + + for (Int_t itrack=0;itrackUncheckedAt(itrack); + UnRegisterClusterTracks(track,trackID1); + } + // + for (Int_t itrack=0;itrackUncheckedAt(itrack); + UnRegisterClusterTracks(track,trackID2); + } + + if (track10->fConstrain&&track10->fChi2MIP[0]fChi2MIP[1]fChi2MIP[2]fChi2MIP[3]fChi2MIP[0]fChi2MIP[1]fChi2MIP[2]fChi2MIP[3]fConstrain&&track20->fChi2MIP[0]fChi2MIP[1]fChi2MIP[2]fChi2MIP[3]fChi2MIP[0]fChi2MIP[1]fChi2MIP[2]fChi2MIP[3]GetClusterIndex(0)); + //if (c->GetQ()>2) c->Use(); + if (c->GetSigmaZ2()>0.1) c->Use(); + c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(1)); + //if (c->GetQ()>2) c->Use(); + if (c->GetSigmaZ2()>0.1) c->Use(); + +} + + +void AliITStrackerMI::AddTrackHypothesys(AliITStrackV2 * track, Int_t esdindex) +{ + //------------------------------------------------------------------ + // add track to the list of hypothesys + //------------------------------------------------------------------ + + if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10); + // + TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); + if (!array) { + array = new TObjArray(10); + fTrackHypothesys.AddAt(array,esdindex); + } + array->AddLast(track); +} + +void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode) +{ + //------------------------------------------------------------------- + // compress array of track hypothesys + // keep only maxsize best hypothesys + //------------------------------------------------------------------- + if (esdindex>fTrackHypothesys.GetEntriesFast()) return; + if (! (fTrackHypothesys.At(esdindex)) ) return; + TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); + Int_t entries = array->GetEntriesFast(); + // + //- find preliminary besttrack as a reference + Float_t minchi2=10000; + Int_t maxn=0; + AliITStrackV2 * besttrack=0; + for (Int_t itrack=0;itrackGetEntriesFast();itrack++){ + AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack); + if (!track) continue; + Float_t chi2 = NormalizedChi2(track,0); + // + Int_t tpcLabel=track->fESDtrack->GetTPCLabel(); + track->SetLabel(tpcLabel); + CookdEdx(track); + track->fFakeRatio=1.; + CookLabel(track,0.); //For comparison only + // + //if (chi2fFakeRatio==0){ + if (chi2GetNumberOfClusters()GetNumberOfClusters(); + if (chi2RemoveAt(itrack); + } + } + if (!besttrack) return; + // + // + //take errors of best track as a reference + Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex); + Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex); + for (Int_t i=0;i<6;i++) { + if (besttrack->fClIndex[i]>0){ + erry[i] = besttrack->fSigmaY[i]; erry[i+6] = besttrack->fSigmaY[i+6]; + errz[i] = besttrack->fSigmaZ[i]; errz[i+6] = besttrack->fSigmaZ[i+6]; + ny[i] = besttrack->fNy[i]; + nz[i] = besttrack->fNz[i]; + } + } + // + // calculate normalized chi2 + // + Float_t * chi2 = new Float_t[entries]; + Int_t * index = new Int_t[entries]; + for (Int_t i=0;iAt(itrack); + if (track){ + track->fChi2MIP[0] = GetNormalizedChi2(track, mode); + if (track->fChi2MIP[0]fChi2MIP[0]; + else + delete array->RemoveAt(itrack); + } + } + // + TMath::Sort(entries,chi2,index,kFALSE); + besttrack = (AliITStrackV2*)array->At(index[0]); + if (besttrack&&besttrack->fChi2MIP[0]fClIndex[i]>0){ + erry[i] = besttrack->fSigmaY[i]; erry[i+6] = besttrack->fSigmaY[i+6]; + errz[i] = besttrack->fSigmaZ[i]; erry[i+6] = besttrack->fSigmaY[i+6]; + ny[i] = besttrack->fNy[i]; + nz[i] = besttrack->fNz[i]; + } + } + } + // + // calculate one more time with updated normalized errors + for (Int_t i=0;iAt(itrack); + if (track){ + track->fChi2MIP[0] = GetNormalizedChi2(track,mode); + if (track->fChi2MIP[0]fChi2MIP[0]-0*(track->GetNumberOfClusters()+track->fNDeadZone); + else + delete array->RemoveAt(itrack); + } + } + entries = array->GetEntriesFast(); + // + if (entries>0){ + TObjArray * newarray = new TObjArray(); + TMath::Sort(entries,chi2,index,kFALSE); + besttrack = (AliITStrackV2*)array->At(index[0]); + if (besttrack){ + // + for (Int_t i=0;i<6;i++){ + if (besttrack->fNz[i]>0&&besttrack->fNy[i]>0){ + erry[i] = besttrack->fSigmaY[i]; erry[i+6] = besttrack->fSigmaY[i+6]; + errz[i] = besttrack->fSigmaZ[i]; errz[i+6] = besttrack->fSigmaZ[i+6]; + ny[i] = besttrack->fNy[i]; + nz[i] = besttrack->fNz[i]; + } + } + besttrack->fChi2MIP[0] = GetNormalizedChi2(besttrack,mode); + Float_t minchi2 = TMath::Min(besttrack->fChi2MIP[0]+5.+besttrack->fNUsed, double(kMaxChi2PerCluster[0])); + Float_t minn = besttrack->GetNumberOfClusters()-3; + Int_t accepted=0; + for (Int_t i=0;iAt(index[i]); + if (!track) continue; + if (accepted>maxcut) break; + track->fChi2MIP[0] = GetNormalizedChi2(track,mode); + if (track->GetNumberOfClusters()<6 && (track->fChi2MIP[0]+track->fNUsed>minchi2)){ + delete array->RemoveAt(index[i]); + continue; + } + if (track->fChi2MIP[0]+track->fNUsedGetNumberOfClusters()>=minn){ + accepted++; + // + newarray->AddLast(array->RemoveAt(index[i])); + for (Int_t i=0;i<6;i++){ + if (nz[i]==0){ + erry[i] = track->fSigmaY[i]; erry[i+6] = track->fSigmaY[i+6]; + errz[i] = track->fSigmaZ[i]; errz[i] = track->fSigmaZ[i+6]; + ny[i] = track->fNy[i]; + nz[i] = track->fNz[i]; + } + } + } + else{ + delete array->RemoveAt(index[i]); + } + } + array->Delete(); + delete fTrackHypothesys.RemoveAt(esdindex); + fTrackHypothesys.AddAt(newarray,esdindex); + } + else{ + array->Delete(); + delete fTrackHypothesys.RemoveAt(esdindex); + } + } + delete [] chi2; + delete [] index; +} + + + +AliITStrackV2 * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackV2 * original, Int_t checkmax) +{ + //------------------------------------------------------------- + // try to find best hypothesy + // currently - minimal chi2 of track+backpropagated track+matching to the tpc track + //------------------------------------------------------------- + if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0; + TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); + if (!array) return 0; + Int_t entries = array->GetEntriesFast(); + if (!entries) return 0; + Float_t minchi2 = 100000; + AliITStrackV2 * besttrack=0; + // + AliITStrackV2 * backtrack = new AliITStrackV2(*original); + AliITStrackV2 * forwardtrack = new AliITStrackV2(*original); + // + for (Int_t i=0;iAt(i); + if (!track) continue; + track->fChi2MIP[1] = 1000000; + track->fChi2MIP[2] = 1000000; + track->fChi2MIP[3] = 1000000; + // + // backtrack + backtrack = new(backtrack) AliITStrackV2(*track); + backtrack->ResetCovariance(); + backtrack->ResetCovariance(); + backtrack->ResetClusters(); + Double_t x = original->GetX(); + if (!RefitAt(x,backtrack,track)) continue; + track->fChi2MIP[1] = NormalizedChi2(backtrack,0); + //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];} + if (track->fChi2MIP[1]>kMaxChi2PerCluster[1]*6.) continue; + track->fChi22 = GetMatchingChi2(backtrack,original); + + if ((track->fConstrain) && track->fChi22>90.) continue; + if ((!track->fConstrain) && track->fChi22>30.) continue; + if ( track->fChi22/track->GetNumberOfClusters()>11.) continue; + + + if (!(track->fConstrain)&&track->fChi2MIP[1]>kMaxChi2PerCluster[1]) continue; + Bool_t isOK=kTRUE; + /* + for (Int_t i=0;i<6;i++){ + if (track->fClIndex[i]>0){ + Double_t dy1 = (track->fDy[i]/track->fSigmaY[i]); + Double_t dz1 = (track->fDz[i]/track->fSigmaZ[i]); + Double_t dy2 = (backtrack->fDy[i]/backtrack->fSigmaY[i]); + Double_t dz2 = (backtrack->fDz[i]/backtrack->fSigmaZ[i]); + if (TMath::Min(dy1*dy1+dz1*dz1,dy2*dy2+dz2*dz2)> kMaxChi2sR[i]) isOK =kFALSE; + track->fDy[i+6] = backtrack->fDy[i]; + track->fDz[i+6] = backtrack->fDz[i]; + track->fSigmaY[i+6] = backtrack->fSigmaY[i]; + track->fSigmaZ[i+6] = backtrack->fSigmaZ[i]; + } + else{ + if (i==5){ + if (track->fClIndex[i-1]>0){ + Double_t dy2 = (backtrack->fDy[i-1]/backtrack->fSigmaY[i-1]); + Double_t dz2 = (backtrack->fDz[i-1]/backtrack->fSigmaZ[i-1]); + if (dy2*dy2+dz2*dz2> kMaxChi2sR[i]) isOK =kFALSE; + } + else isOK = kFALSE; + } + } + } + */ + if(!isOK) continue; + // + //forward track - without constraint + forwardtrack = new(forwardtrack) AliITStrackV2(*original); + forwardtrack->ResetClusters(); + x = track->GetX(); + if (!RefitAt(x,forwardtrack,track)) continue; + track->fChi2MIP[2] = NormalizedChi2(forwardtrack,0); + if (track->fChi2MIP[2]>kMaxChi2PerCluster[2]*6.0) continue; + if (!(track->fConstrain)&&track->fChi2MIP[2]>kMaxChi2PerCluster[2]) continue; + + track->fD[0] = forwardtrack->GetD(GetX(),GetY()); + track->fD[1] = forwardtrack->GetZat(GetX())-GetZ(); + forwardtrack->fD[0] = track->fD[0]; + forwardtrack->fD[1] = track->fD[1]; + { + Int_t list[6]; + AliITSclusterV2* clist[6]; + track->fChi2MIP[4] = GetNumberOfSharedClusters(track,esdindex,list,clist); + if ( (!track->fConstrain) && track->fChi2MIP[4]>1.0) continue; + } + + track->fChi2MIP[3] = GetInterpolatedChi2(forwardtrack,backtrack); + if ( (track->fChi2MIP[3]>6.*kMaxChi2PerCluster[3])) continue; + if ( (!track->fConstrain) && (track->fChi2MIP[3]>2*kMaxChi2PerCluster[3])) { + track->fChi2MIP[3]=1000; + continue; + } + Double_t chi2 = track->fChi2MIP[0]+track->fNUsed; + // + for (Int_t ichi=0;ichi<5;ichi++){ + forwardtrack->fChi2MIP[ichi] = track->fChi2MIP[ichi]; + } + if (chi2 < minchi2){ + besttrack = new AliITStrackV2(*forwardtrack); + besttrack->SetLabel(track->GetLabel()); + besttrack->fFakeRatio = track->fFakeRatio; + minchi2 = chi2; + original->fD[0] = forwardtrack->GetD(GetX(),GetY()); + original->fD[1] = forwardtrack->GetZat(GetX())-GetZ(); + } + } + delete backtrack; + delete forwardtrack; + Int_t accepted=0; + for (Int_t i=0;iAt(i); + if (!track) continue; + if (accepted>checkmax || track->fChi2MIP[3]>kMaxChi2PerCluster[3]*6. || + (track->GetNumberOfClusters()GetNumberOfClusters()-1.)|| + track->fChi2MIP[0]>besttrack->fChi2MIP[0]+2.*besttrack->fNUsed+3.){ + delete array->RemoveAt(i); + continue; + } + else{ + accepted++; + } + } + // + array->Compress(); + SortTrackHypothesys(esdindex,checkmax,1); + array = (TObjArray*) fTrackHypothesys.At(esdindex); + besttrack = (AliITStrackV2*)array->At(0); + if (!besttrack) return 0; + besttrack->fChi2MIP[8]=0; + fBestTrackIndex[esdindex]=0; + entries = array->GetEntriesFast(); + AliITStrackV2 *longtrack =0; + minchi2 =1000; + Float_t minn=besttrack->GetNumberOfClusters()+besttrack->fNDeadZone; + for (Int_t itrack=entries-1;itrack>0;itrack--){ + AliITStrackV2 * track = (AliITStrackV2*)array->At(itrack); + if (!track->fConstrain) continue; + if (track->GetNumberOfClusters()+track->fNDeadZonefChi2MIP[0]-besttrack->fChi2MIP[0]>0.0) continue; + if (track->fChi2MIP[0]>4.) continue; + minn = track->GetNumberOfClusters()+track->fNDeadZone; + longtrack =track; + } + //if (longtrack) besttrack=longtrack; + + Int_t list[6]; + AliITSclusterV2 * clist[6]; + Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist); + if (besttrack->fConstrain&&besttrack->fChi2MIP[0]fChi2MIP[1]fChi2MIP[2]fChi2MIP[3]fChi2MIP[0]fChi2MIP[1]fChi2MIP[2]fChi2MIP[3]0.0){ + Int_t nshared; + Int_t overlist[6]; + Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist); + if (sharedtrack>=0){ + // + besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5); + if (besttrack){ + shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist); + } + else return 0; + } + } + + if (shared>2.5) return 0; + if (shared>1.0) return besttrack; + // + // don't sign clusters if not gold track + Double_t deltad = besttrack->GetD(GetX(),GetY()); + Double_t deltaz = besttrack->GetZat(GetX()) - GetZ(); + Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz); + if (deltaprim>0.1 && (fConstraint[fPass])) return besttrack; + if (TMath::Abs(deltad)>0.1) return besttrack; + if (TMath::Abs(deltaz)>0.1) return besttrack; + if (besttrack->fChi2MIP[0]>4.) return besttrack; + if (besttrack->fChi2MIP[1]>4.) return besttrack; + if (besttrack->fChi2MIP[2]>4.) return besttrack; + if (besttrack->fChi2MIP[3]>4.) return besttrack; + if (fConstraint[fPass]){ + // + // sign clusters + // + Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex); + for (Int_t i=0;i<6;i++){ + Int_t index = besttrack->fClIndex[i]; + if (index<=0) continue; + Int_t ilayer = (index & 0xf0000000) >> 28; + if (besttrack->fSigmaY[ilayer]<0.00000000001) continue; + AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index); + if (!c) continue; + if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue; + if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track + if ( c->GetNz()> nz[i]+0.7) continue; //shared track + if ( ilayer>2&& besttrack->fNormQ[ilayer]/besttrack->fExpQ>1.5) continue; + //if ( c->GetNy()> ny[i]+0.7) continue; //shared track + + Bool_t cansign = kTRUE; + for (Int_t itrack=0;itrackAt(i); + if (!track) continue; + if (track->fChi2MIP[0]>besttrack->fChi2MIP[0]+2.*shared+1.) break; + if ( (track->fClIndex[ilayer]>0) && (track->fClIndex[ilayer]!=besttrack->fClIndex[ilayer])){ + cansign = kFALSE; + break; + } + } + if (cansign){ + if (TMath::Abs(besttrack->fDy[ilayer]/besttrack->fSigmaY[ilayer])>3.) continue; + if (TMath::Abs(besttrack->fDz[ilayer]/besttrack->fSigmaZ[ilayer])>3.) continue; + if (!c->IsUsed()) c->Use(); + } + } + } + return besttrack; +} + + + +void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks) +{ + // + // get "best" hypothesys + //for (Int_t ilayer=0;ilayer<6;ilayer++) fgLayers[ilayer].ResetWeights(); + + + Int_t nentries = itsTracks.GetEntriesFast(); + for (Int_t i=0;iGetEntriesFast()<=0) continue; + // + AliITStrackV2* longtrack=0; + Float_t minn=0; + Float_t maxchi2=1000; + for (Int_t j=0;jGetEntriesFast();j++){ + AliITStrackV2* track = (AliITStrackV2*)array->At(j); + if (!track) continue; + if (track->GetNumberOfClusters()+track->fNDeadZoneGetNumberOfClusters()+track->fNDeadZone>minn) maxchi2 = track->fChi2MIP[0]; + if (track->fChi2MIP[0]>maxchi2) continue; + minn= track->GetNumberOfClusters()+track->fNDeadZone; + maxchi2 = track->fChi2MIP[0]; + longtrack=track; + break; + } + AliITStrackV2 * besttrack = (AliITStrackV2*)array->At(0); + if (!longtrack) {longtrack = besttrack;} + else besttrack= longtrack; + if (besttrack){ + Int_t list[6]; + AliITSclusterV2 * clist[6]; + Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist); + // + track->fNUsed = shared; + track->fNSkipped = besttrack->fNSkipped; + track->fChi2MIP[0] = besttrack->fChi2MIP[0]; + if (shared>0){ + Int_t nshared; + Int_t overlist[6]; + // + Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist); + //if (sharedtrack==-1) sharedtrack=0; + if (sharedtrack>=0){ + besttrack = GetBest2Tracks(i,sharedtrack,10,5.5); + } + } + if (besttrack&&fConstraint[fPass]) + UpdateESDtrack(besttrack,AliESDtrack::kITSin); + //if (besttrack&&besttrack->fConstrain) + // UpdateESDtrack(besttrack,AliESDtrack::kITSin); + if (besttrack->fChi2MIP[0]+besttrack->fNUsed>1.5){ + if ( (TMath::Abs(besttrack->fD[0])>0.1) && fConstraint[fPass]) { + track->fReconstructed= kFALSE; + } + if ( (TMath::Abs(besttrack->fD[1])>0.1) && fConstraint[fPass]){ + track->fReconstructed= kFALSE; + } + } + + } + } +} + + +void AliITStrackerMI::CookLabel(AliITStrackV2 *track,Float_t wrong) const { + //-------------------------------------------------------------------- + //This function "cooks" a track label. If label<0, this track is fake. + //-------------------------------------------------------------------- + Int_t tpcLabel=-1; + + if ( track->fESDtrack) tpcLabel = TMath::Abs(track->fESDtrack->GetTPCLabel()); + + track->fChi2MIP[9]=0; + Int_t nwrong=0; + for (Int_t i=0;iGetNumberOfClusters();i++){ + Int_t cindex = track->GetClusterIndex(i); + Int_t l=(cindex & 0xf0000000) >> 28; + AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(cindex); + Int_t isWrong=1; + for (Int_t ind=0;ind<3;ind++){ + if (tpcLabel>0) + if (cl->GetLabel(ind)==tpcLabel) isWrong=0; + } + track->fChi2MIP[9]+=isWrong*(2<fFakeRatio = double(nwrong)/double(track->GetNumberOfClusters()); + if (tpcLabel>0){ + if (track->fFakeRatio>wrong) track->fLab = -tpcLabel; + else + track->fLab = tpcLabel; + } + +} + + + +void AliITStrackerMI::CookdEdx(AliITStrackV2* track) +{ + // + // + // Int_t list[6]; + //AliITSclusterV2 * clist[6]; + // Int_t shared = GetNumberOfSharedClusters(track,index,list,clist); + Float_t dedx[4]; + Int_t accepted=0; + track->fChi2MIP[9]=0; + for (Int_t i=0;iGetNumberOfClusters();i++){ + Int_t cindex = track->GetClusterIndex(i); + Int_t l=(cindex & 0xf0000000) >> 28; + AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(cindex); + Int_t lab = TMath::Abs(track->fESDtrack->GetTPCLabel()); + Int_t isWrong=1; + for (Int_t ind=0;ind<3;ind++){ + if (cl->GetLabel(ind)==lab) isWrong=0; + } + track->fChi2MIP[9]+=isWrong*(2<3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue; //shared track + //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue; + //if (l<4&& !(cl->GetType()==1)) continue; + dedx[accepted]= track->fdEdxSample[i]; + //dedx[accepted]= track->fNormQ[l]; + accepted++; + } + if (accepted<1) { + track->SetdEdx(0); + return; + } + Int_t indexes[4]; + TMath::Sort(accepted,dedx,indexes,kFALSE); + Double_t low=0.; + Double_t up=0.51; + Double_t nl=low*accepted, nu =up*accepted; + Float_t sumamp = 0; + Float_t sumweight =0; + for (Int_t i=0; inu-1) weight = TMath::Max(nu-i,0.); + sumamp+= dedx[indexes[i]]*weight; + sumweight+=weight; + } + track->SetdEdx(sumamp/sumweight); +} + + +void AliITStrackerMI::MakeCoeficients(Int_t ntracks){ + // + // + fCoeficients = new Float_t[ntracks*48]; + for (Int_t i=0;iGetTgl(); + Float_t phi = track->GetSnp(); + phi = TMath::Sqrt(phi*phi/(1.-phi*phi)); + GetError(layer,cluster,theta,phi,track->fExpQ,erry,errz); + Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz); + Float_t ny,nz; + GetNTeor(layer,cluster, theta,phi,ny,nz); + Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny; + if (delta>1){ + chi2+=0.5*TMath::Min(delta/2,2.); + chi2+=2.*cluster->GetDeltaProbability(); + } + // + track->fNy[layer] =ny; + track->fNz[layer] =nz; + track->fSigmaY[layer] = erry; + track->fSigmaZ[layer] = errz; + //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi); + track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt((1.+ track->fP3*track->fP3)/(1.- track->fP2*track->fP2)); + return chi2; + +} + +Int_t AliITStrackerMI::UpdateMI(AliITStrackV2* track, const AliITSclusterV2* cl,Double_t chi2,Int_t index) +{ + // + // + // + Int_t layer = (index & 0xf0000000) >> 28; + track->fClIndex[layer] = index; + if ( (layer>1) &&track->fNormQ[layer]/track->fExpQ<0.5 ) { + chi2+= (0.5-track->fNormQ[layer]/track->fExpQ)*10.; + track->fdEdxMismatch+=(0.5-track->fNormQ[layer]/track->fExpQ)*10.; + } + return track->UpdateMI(cl->GetY(),cl->GetZ(),track->fSigmaY[layer],track->fSigmaZ[layer],chi2,index); +} + +void AliITStrackerMI::GetNTeor(Int_t layer, const AliITSclusterV2* /*cl*/, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz) +{ + // + //get "mean shape" + // + if (layer==0){ + ny = 1.+TMath::Abs(phi)*3.2; + nz = 1.+TMath::Abs(theta)*0.34; + return; + } + if (layer==1){ + ny = 1.+TMath::Abs(phi)*3.2; + nz = 1.+TMath::Abs(theta)*0.28; + return; + } + + if (layer>3){ + ny = 2.02+TMath::Abs(phi)*1.95; + nz = 2.02+TMath::Abs(phi)*2.35; + return; + } + ny = 6.6-2.7*abs(phi); + nz = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta); +} + + + +Int_t AliITStrackerMI::GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi,Float_t expQ, Float_t &erry, Float_t &errz) +{ + //calculate cluster position error + // + Float_t nz,ny; + GetNTeor(layer, cl,theta,phi,ny,nz); + erry = TMath::Sqrt(cl->GetSigmaY2()); + errz = TMath::Sqrt(cl->GetSigmaZ2()); + // + // PIXELS + if (layer<2){ + + if (TMath::Abs(ny-cl->GetNy())>0.6) { + if (nyGetNy()){ + erry*=0.4+TMath::Abs(ny-cl->GetNy()); + errz*=0.4+TMath::Abs(ny-cl->GetNy()); + }else{ + erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy()); + errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy()); + } + } + if (TMath::Abs(nz-cl->GetNz())>1.) { + erry*=TMath::Abs(nz-cl->GetNz()); + errz*=TMath::Abs(nz-cl->GetNz()); + } + erry*=0.85; + errz*=0.85; + erry= TMath::Min(erry,float(0.005)); + errz= TMath::Min(errz,float(0.03)); + return 10; + } + +//STRIPS + if (layer>3){ + if (cl->GetNy()==100||cl->GetNz()==100){ + erry = 0.004; + errz = 0.2; + return 100; + } + if (cl->GetNy()+cl->GetNz()>12){ + erry = 0.06; + errz = 0.57; + return 100; + } + Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi)); + Float_t chargematch = TMath::Max(double(normq/expQ),2.); + // + if (cl->GetType()==1 || cl->GetType()==10 ){ + if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()GetNy()+cl->GetNz()GetNy()+cl->GetNz()-nz+ny),0.15); + return 103; + } + if (cl->GetType()==2 || cl->GetType()==11 ){ + erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05); + errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5); + return 104; + } + + if (cl->GetType()>100 ){ + if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){ + errz = 0.05; + erry = 0.00096; + return 105; + } + if (cl->GetNy()+cl->GetNz()-nz-ny<1){ + errz = 0.10; + erry = 0.0025; + return 106; + } + + errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4); + erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05); + return 107; + } + Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz; + if (diff<1) diff=1; + if (diff>4) diff=4; + + if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){ + errz = 0.14*diff; + erry = 0.003*diff; + return 108; + } + erry = 0.04*diff; + errz = 0.06*diff; + return 109; + } + //DRIFTS + Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi)); + Float_t chargematch = normq/expQ; + Float_t factorz=1; + Int_t cnz = cl->GetNz()%10; + //charge match + if (cl->GetType()==1){ + if (chargematch<1.25){ + erry = 0.0028*(1.+6./cl->GetQ()); // gold clusters + } + else{ + erry = 0.003*chargematch; + if (cl->GetNz()==3) erry*=1.5; + } + if (chargematch<1.0){ + errz = 0.0011*(1.+6./cl->GetQ()); + } + else{ + errz = 0.002*(1+2*(chargematch-1.)); + } + if (cnz>nz+0.6) { + erry*=(cnz-nz+0.5); + errz*=1.4*(cnz-nz+0.5); + } + } + if (cl->GetType()>1){ + if (chargematch<1){ + erry = 0.00385*(1.+6./cl->GetQ()); // gold clusters + errz = 0.0016*(1.+6./cl->GetQ()); + } + else{ + errz = 0.0014*(1+3*(chargematch-1.)); + erry = 0.003*(1+3*(chargematch-1.)); + } + if (cnz>nz+0.6) { + erry*=(cnz-nz+0.5); + errz*=1.4*(cnz-nz+0.5); + } + } + + if (TMath::Abs(cl->GetY())>2.5){ + factorz*=1+2*(TMath::Abs(cl->GetY())-2.5); + } + if (TMath::Abs(cl->GetY())<1){ + factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.); + } + factorz= TMath::Min(factorz,float(4.)); + errz*=factorz; + + erry= TMath::Min(erry,float(0.05)); + errz= TMath::Min(errz,float(0.05)); + return 200; +} + + + +void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz) +{ + // + // + Int_t entries = ClusterArray->GetEntriesFast(); + if (entries<4) return; + AliITSclusterV2* cluster = (AliITSclusterV2*)ClusterArray->At(0); + Int_t layer = cluster->GetLayer(); + if (layer>1) return; + Int_t index[10000]; + Int_t ncandidates=0; + Float_t r = (layer>0)? 7:4; + // + for (Int_t i=0;iAt(i); + Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r); + if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue; + index[ncandidates] = i; //candidate to belong to delta electron track + ncandidates++; + if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) { + cl0->SetDeltaProbability(1); + } + } + // + // + // + for (Int_t i=0;iAt(index[i]); + if (cl0->GetDeltaProbability()>0.8) continue; + // + Int_t ncl = 0; + Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw; + sumy=sumz=sumy2=sumyz=sumw=0.0; + for (Int_t j=0;jAt(index[j]); + // + Float_t dz = cl0->GetZ()-cl1->GetZ(); + Float_t dy = cl0->GetY()-cl1->GetY(); + if (TMath::Sqrt(dz*dz+dy*dy)<0.2){ + Float_t weight = cl1->GetNy()+cl1->GetNz()-2; + y[ncl] = cl1->GetY(); + z[ncl] = cl1->GetZ(); + sumy+= y[ncl]*weight; + sumz+= z[ncl]*weight; + sumy2+=y[ncl]*y[ncl]*weight; + sumyz+=y[ncl]*z[ncl]*weight; + sumw+=weight; + ncl++; + } + } + if (ncl<4) continue; + Float_t det = sumw*sumy2 - sumy*sumy; + Float_t delta=1000; + if (TMath::Abs(det)>0.01){ + Float_t z0 = (sumy2*sumz - sumy*sumyz)/det; + Float_t k = (sumyz*sumw - sumy*sumz)/det; + delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY())); + } + else{ + Float_t z0 = sumyz/sumy; + delta = TMath::Abs(cl0->GetZ()-z0); + } + if ( delta<0.05) { + cl0->SetDeltaProbability(1-20.*delta); + } + } +} + + +void AliITStrackerMI::UpdateESDtrack(AliITStrackV2* track, ULong_t flags) +{ + // + // + track->UpdateESDtrack(flags); + track->fESDtrack->SetITStrack(new AliITStrackV2(*track)); +} + + + + +void AliITStrackerMI::FindV0(AliESD */*event*/) +{ + // + // fast V0 finder + // + AliHelix helixes[30000]; + TObjArray trackarray(30000); + Float_t dist[30000]; + AliITSRecV0Info vertex; + // + Int_t entries = fTrackHypothesys.GetEntriesFast(); + for (Int_t i=0;iAt(fBestTrackIndex[i]); + if (track){ + dist[i] = TMath::Sqrt(track->fD[0]*track->fD[0]+track->fD[1]*track->fD[1]); + trackarray.AddAt(track,i); + new (&helixes[i]) AliHelix(*track); + } + } + for (Int_t itrack0=0;itrack0fP4*track0->fP4>0) continue; //the same sign + AliHelix *h1 = &helixes[itrack0]; + AliHelix *h2 = &helixes[itrack1]; + + Double_t distance = TestV0(h1,h2,&vertex); + if (distance>0.4) continue; + if (vertex.fRr<0.3) continue; + if (vertex.fRr>27) continue; + Float_t v[3]={GetX(),GetY(),GetZ()}; + vertex.Update(v,track0->fMass, track1->fMass); + + if ((TMath::Abs(vertex.fInvMass-0.4976)<0.05 || TMath::Abs(vertex.fInvMass-0.0)<0.1|| TMath::Abs(vertex.fInvMass-1.1)<0.1)) { + if (vertex.fPointAngle<0.85) continue; + } + else{ + if (vertex.fPointAngle<0.98) continue; + } + if (TMath::Abs(TMath::Abs(track0->fLab)-TMath::Abs(track1->fLab))<2){ + Float_t mindist = FindBestPair(itrack0,itrack1,&vertex); + if (mindist>1) FindBestPair(itrack0,itrack1,&vertex); + } + } + } +} + +Double_t AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliITSRecV0Info *vertex) +{ + // + // try to find best pair from the tree of track hyp. + // + TObjArray * array0 = (TObjArray*)fTrackHypothesys.At(esdtrack0); + Int_t entries0 = array0->GetEntriesFast(); + TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(esdtrack1); + Int_t entries1 = array1->GetEntriesFast(); + // + // + AliITSRecV0Info v0; + AliITSRecV0Info vbest; + Double_t criticalradius = vertex->fRr; + Double_t mindist =1000; + Int_t bestpair[2]; + // + for (Int_t itrack0=0;itrack0At(itrack0); + if (!track0) continue; + if (track0->fXfX>criticalradius+5) continue; + for (Int_t itrack1=0;itrack1At(itrack1); + if (!track1) continue; + if (track1->fXfX>criticalradius+5) continue; + + AliHelix h0(*track0); + AliHelix h1(*track1); + Double_t distance = TestV0(&h0,&h1,&v0); + if (v0.fRr>30) continue; + if (v0.fRr<0.2) continue; + Float_t v[3]={GetX(),GetY(),GetZ()}; + v0.Update(v,track0->fMass, track1->fMass); + if (distanceDump(); + } + return mindist; +} + + +void AliITSRecV0Info::Update(Float_t vertex[3], Float_t mass1, Float_t mass2) +{ + Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]}; + Float_t p[3] = {fPdr[0]+fPm[0], fPdr[1]+fPm[1],fPdr[2]+fPm[2]}; + + Float_t vnorm2 = v[0]*v[0]+v[1]*v[1]; + Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2); + vnorm2 = TMath::Sqrt(vnorm2); + Float_t pnorm2 = p[0]*p[0]+p[1]*p[1]; + Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2); + pnorm2 = TMath::Sqrt(pnorm2); + + fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2); + fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3); + fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3); + + Float_t e1 = TMath::Sqrt(mass1*mass1+ + fPdr[0]*fPdr[0]+ + fPdr[1]*fPdr[1]+ + fPdr[2]*fPdr[2]); + Float_t e2 = TMath::Sqrt(mass2*mass2+ + fPm[0]*fPm[0]+ + fPm[1]*fPm[1]+ + fPm[2]*fPm[2]); + + fInvMass = + (fPm[0]+fPdr[0])*(fPm[0]+fPdr[0])+ + (fPm[1]+fPdr[1])*(fPm[1]+fPdr[1])+ + (fPm[2]+fPdr[2])*(fPm[2]+fPdr[2]); + + fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass); + + +} + + + +Double_t AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliITSRecV0Info *vertex) +{ + // + // test the helixes for the distnce calculate vertex characteristic + // + Float_t distance1,distance2; + AliHelix & dhelix1 = *helix1; + Double_t pp[3],xx[3]; + dhelix1.GetMomentum(0,pp,0); + dhelix1.Evaluate(0,xx); + AliHelix &mhelix = *helix2; + // + //find intersection linear + // + Double_t phase[2][2],radius[2]; + Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius); + Double_t delta1=10000,delta2=10000; + + if (points>0){ + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + distance1 = TMath::Min(delta1,delta2); + vertex->fDist1 = TMath::Sqrt(distance1); + // + //find intersection parabolic + // + points = dhelix1.GetRPHIintersections(mhelix, phase, radius); + delta1=10000,delta2=10000; + + if (points>0){ + dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); + } + if (points==2){ + dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); + } + + distance2 = TMath::Min(delta1,delta2); + vertex->fDist2 = TMath::Sqrt(distance2); + if (distance2<100){ + if (delta1fXr); + dhelix1.GetMomentum(phase[0][0],vertex->fPdr); + mhelix.GetMomentum(phase[0][1],vertex->fPm); + dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],vertex->fAngle); + vertex->fRr = TMath::Sqrt(radius[0]); + } + else{ + dhelix1.Evaluate(phase[1][0],vertex->fXr); + dhelix1.GetMomentum(phase[1][0],vertex->fPdr); + mhelix.GetMomentum(phase[1][1],vertex->fPm); + dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],vertex->fAngle); + vertex->fRr = TMath::Sqrt(radius[1]); + } + } + // + // + return vertex->fDist2; +} + diff --git a/ITS/AliITStrackerMI.h b/ITS/AliITStrackerMI.h new file mode 100644 index 00000000000..43b39b8e51e --- /dev/null +++ b/ITS/AliITStrackerMI.h @@ -0,0 +1,321 @@ +#ifndef ALIITSTRACKERMI_H +#define ALIITSTRACKERMI_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------------------------- +// ITS tracker +// reads AliITSclusterMI clusters and creates AliITStrackV2 tracks +// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch +//------------------------------------------------------------------------- + +#include + +#include "AliTracker.h" +#include "AliITSrecoV2.h" +#include "AliITStrackV2.h" +#include "AliITSclusterV2.h" + +class AliITSclusterV2; +class AliESD; +class AliITSgeom; +class TTree; +class AliHelix; +class AliV0vertex; + + + +class AliITSRecV0Info: public TObject { +public: + void Update(Float_t vertex[3], Float_t mass1, Float_t mass2); + Double_t fDist1; //info about closest distance according closest MC - linear DCA + Double_t fDist2; //info about closest distance parabolic DCA + Double_t fInvMass; //reconstructed invariant mass - + // + Double_t fPdr[3]; //momentum at vertex daughter - according approx at DCA + Double_t fXr[3]; //rec. position according helix + // + Double_t fPm[3]; //momentum at the vertex mother + Double_t fAngle[3]; //three angles + Double_t fRr; // rec position of the vertex + Int_t fLab[2]; //MC label of the partecle + Float_t fPointAngleFi; //point angle fi + Float_t fPointAngleTh; //point angle theta + Float_t fPointAngle; //point angle full + ClassDef(AliITSRecV0Info,1) // container for +}; + + + + +//------------------------------------------------------------------------- +class AliITStrackerMI : public AliTracker { +public: + AliITStrackerMI():AliTracker(){} + AliITStrackerMI(const AliITSgeom *geom); + AliCluster *GetCluster(Int_t index) const; + AliITSclusterV2 *GetClusterLayer(Int_t layn, Int_t ncl) const + {return fgLayers[layn].GetCluster(ncl);} + Int_t GetNumberOfClustersLayer(Int_t layn) const + {return fgLayers[layn].GetNumberOfClusters();} + Int_t LoadClusters(TTree *cf); + void UnloadClusters(); + Int_t Clusters2Tracks(TTree *in, TTree *out); + Int_t Clusters2Tracks(AliESD *event); + Int_t PropagateBack(AliESD *event); + Int_t RefitInward(AliESD *event); + Bool_t RefitAt(Double_t x, AliITStrackV2 *seed, const AliITStrackV2 *t); + void SetupFirstPass(Int_t *flags, Double_t *cuts=0); + void SetupSecondPass(Int_t *flags, Double_t *cuts=0); + + void SetLastLayerToTrackTo(Int_t l=0) {fLastLayerToTrackTo=l;} + void SetLayersNotToSkip(Int_t *l); + void UseClusters(const AliKalmanTrack *t, Int_t from=0) const; + void GetNTeor(Int_t layer, const AliITSclusterV2* cl, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz); + Int_t GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi, Float_t expQ, Float_t &erry, Float_t &errz); + Double_t GetPredictedChi2MI(AliITStrackV2* track, const AliITSclusterV2 *cluster,Int_t layer); + Int_t UpdateMI(AliITStrackV2* track, const AliITSclusterV2* cl,Double_t chi2,Int_t layer); + class AliITSdetector { + public: + AliITSdetector(){} + AliITSdetector(Double_t r,Double_t phi) {fR=r; fPhi=phi; fSinPhi = TMath::Sin(phi); fCosPhi = TMath::Cos(phi); + fYmin=10000;fYmax=-1000; fZmin=10000;fZmax=-1000;} + inline void GetGlobalXYZ( const AliITSclusterV2 *cl, Double_t xyz[3]) const; + Double_t GetR() const {return fR;} + Double_t GetPhi() const {return fPhi;} + Double_t GetYmin() const {return fYmin;} + Double_t GetYmax() const {return fYmax;} + Double_t GetZmin() const {return fZmin;} + Double_t GetZmax() const {return fZmax;} + void SetYmin(Double_t min) {fYmin = min;} + void SetYmax(Double_t max) {fYmax = max;} + void SetZmin(Double_t min) {fZmin = min;} + void SetZmax(Double_t max) {fZmax = max;} + private: + Double_t fR; // polar coordinates + Double_t fPhi; // of this detector + Double_t fSinPhi; // sin of phi; + Double_t fCosPhi; // cos of phi + Double_t fYmin; // local y minimal + Double_t fYmax; // local max y + Double_t fZmin; // local z min + Double_t fZmax; // local z max + }; + + class AliITSlayer { + friend class AliITStrackerMI; + public: + AliITSlayer(); + AliITSlayer(Double_t r, Double_t p, Double_t z, Int_t nl, Int_t nd); + ~AliITSlayer(); + Int_t InsertCluster(AliITSclusterV2 *c); + void SortClusters(); + void ResetClusters(); + void ResetWeights(); + void SelectClusters(Double_t zmi,Double_t zma,Double_t ymi,Double_t yma); + const AliITSclusterV2 *GetNextCluster(Int_t &ci); + void ResetRoad(); + Double_t GetRoad() const {return fRoad;} + Double_t GetR() const {return fR;} + Int_t FindClusterIndex(Float_t z) const; + AliITSclusterV2 *GetCluster(Int_t i) const {return iGetLabel(); + if (tpcLabel<0) return; + AliTracker::CookLabel(t,wrong); + if (tpcLabel!=TMath::Abs(t->GetLabel())){ + t->SetFakeRatio(1.); + } + if (tpcLabel !=t->GetLabel()) { + t->SetLabel(-tpcLabel); + } +} + +inline Double_t AliITStrackerMI::NormalizedChi2(AliITStrackV2 * track, Int_t layer) +{ + //-------------------------------------------------------------------- + //get normalize chi2 + //-------------------------------------------------------------------- + track->fNormChi2[layer] = 2.*track->fNSkipped+0.25*track->fNDeadZone+track->fdEdxMismatch+track->GetChi2()/ + //track->fNormChi2[layer] = 2.*track->fNSkipped+0.25*track->fNDeadZone+track->fdEdxMismatch+track->fChi22/ + TMath::Max(double(track->GetNumberOfClusters()-track->fNSkipped), + 1./(1.+track->fNSkipped)); + return track->fNormChi2[layer]; +} +inline void AliITStrackerMI::AliITSdetector::GetGlobalXYZ(const AliITSclusterV2 *cl, Double_t xyz[3]) const +{ + // + // get cluster coordinates in global cooordinate + // + xyz[2] = cl->GetZ(); + xyz[0] = fR*fCosPhi - cl->GetY()*fSinPhi; + xyz[1] = fR*fSinPhi + cl->GetY()*fCosPhi; +} +#endif diff --git a/ITS/AliITStrackerSA.cxx b/ITS/AliITStrackerSA.cxx index f292403e6cf..ee8f1af5522 100644 --- a/ITS/AliITStrackerSA.cxx +++ b/ITS/AliITStrackerSA.cxx @@ -41,13 +41,13 @@ ClassImp(AliITStrackerSA) //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA():AliITStrackerV2(){ +AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(){ // Default constructor Init(); } //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom):AliITStrackerV2(geom) +AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom):AliITStrackerMI(geom) { // Standard constructor (Vertex is known and passed to this obj.) Init(); @@ -57,7 +57,7 @@ AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom):AliITStrackerV2(geom) } //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliESDVertex *vert):AliITStrackerV2(geom) +AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliESDVertex *vert):AliITStrackerMI(geom) { // Standard constructor (Vertex is known and passed to this obj.) Init(); @@ -68,7 +68,7 @@ AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliESDVertex *vert):AliITStra //______________________________________________________________________ AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA &trkr) : - AliITStrackerV2(trkr) { + AliITStrackerMI(trkr) { // Copy constructor // Copies are not allowed. The method is protected to avoid misuse. Error("AliITStrackerSA","Copy constructor not allowed\n"); @@ -84,7 +84,7 @@ AliITStrackerSA& AliITStrackerSA::operator=(const } //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliITSVertexer *vertexer):AliITStrackerV2(geom) +AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliITSVertexer *vertexer):AliITStrackerMI(geom) { // Standard constructor (Vertex is unknown - vertexer is passed to this obj) Init(); @@ -94,7 +94,7 @@ AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliITSVertexer *vertexer):Ali } //____________________________________________________________________________ -AliITStrackerSA::AliITStrackerSA(AliITStrackerSA& tracker):AliITStrackerV2(){ +AliITStrackerSA::AliITStrackerSA(AliITStrackerSA& tracker):AliITStrackerMI(){ // Copy constructor fPhiEstimate = tracker.fPhiEstimate; for(Int_t i=0;i<2;i++){ @@ -815,7 +815,7 @@ AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVert 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() + //fit with Kalman filter using AliITStrackerMI::RefitAt() AliITStrackV2* ot = new AliITStrackV2(*trac); @@ -902,6 +902,7 @@ AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVert Int_t numberofpoints; if(fSixPoints) numberofpoints=6; else numberofpoints=5; + CookLabel(otrack,0.); //MI change - to see fake ratio Int_t label = Label(cl0->GetLabel(0),cl1->GetLabel(0), cl2->GetLabel(0),cl3->GetLabel(0), cl4->GetLabel(0),labl[0], diff --git a/ITS/AliITStrackerSA.h b/ITS/AliITStrackerSA.h index 41480f58da8..5ed8eb51c54 100644 --- a/ITS/AliITStrackerSA.h +++ b/ITS/AliITStrackerSA.h @@ -3,7 +3,7 @@ -#include "AliITStrackerV2.h" +#include "AliITStrackerMI.h" /* Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ @@ -20,7 +20,7 @@ class AliESDVertex; class AliITSVertexer; class TTree; -class AliITStrackerSA : public AliITStrackerV2 { +class AliITStrackerSA : public AliITStrackerMI { public: @@ -31,14 +31,14 @@ class AliITStrackerSA : public AliITStrackerV2 { AliITStrackerSA(AliITSgeom *geom,AliITSVertexer *vertexer); AliITStrackerSA(AliITStrackerSA& tracker); virtual ~AliITStrackerSA(); - virtual Int_t Clusters2Tracks(AliESD *event){Int_t rc = AliITStrackerV2::Clusters2Tracks(event); if(!rc) rc=FindTracks(event); return rc;} + virtual Int_t Clusters2Tracks(AliESD *event){Int_t rc = AliITStrackerMI::Clusters2Tracks(event); if(!rc) rc=FindTracks(event); return rc;} Int_t FindTracks(AliESD* event); void FindTracks(TTree *out,Int_t evnumber=0); AliITStrackV2* FitTrack(AliITStrackSA* tr,Double_t* primaryVertex, Double_t *errorprimvert); AliITStrackV2* FindTrackLowChiSquare(TObjArray* tracklist, Int_t dim) const; - Int_t LoadClusters(TTree *cf) {Int_t rc=AliITStrackerV2::LoadClusters(cf); SetClusterTree(cf);SetSixPoints(kTRUE); return rc;} + Int_t LoadClusters(TTree *cf) {Int_t rc=AliITStrackerMI::LoadClusters(cf); SetClusterTree(cf);SetSixPoints(kTRUE); return rc;} void SetVertex(AliESDVertex *vtx){fVert = vtx;} void SetClusterTree(TTree * itscl){fITSclusters = itscl;} void SetSixPoints(Bool_t sp = kTRUE){fSixPoints = sp;} diff --git a/ITS/ITSLinkDef.h b/ITS/ITSLinkDef.h index 7d7f4c173bc..b9969935513 100644 --- a/ITS/ITSLinkDef.h +++ b/ITS/ITSLinkDef.h @@ -139,6 +139,8 @@ #pragma link C++ class AliITSclusterV2+; #pragma link C++ class AliITStrackV2+; #pragma link C++ class AliITStrackerV2+; +#pragma link C++ class AliITStrackerMI+; +#pragma link C++ class AliITSRecV0Info+; #pragma link C++ class AliV0vertex+; #pragma link C++ class AliV0vertexer+; #pragma link C++ class AliCascadeVertex+; diff --git a/ITS/libITS.pkg b/ITS/libITS.pkg index e5f8866c702..a78c700aff3 100644 --- a/ITS/libITS.pkg +++ b/ITS/libITS.pkg @@ -63,6 +63,7 @@ SRCS = AliITS.cxx \ AliITSclusterV2.cxx \ AliITStrackV2.cxx \ AliITStrackerV2.cxx \ + AliITStrackerMI.cxx \ AliITSVertexer.cxx \ AliITSVertexerIons.cxx \ AliITSVertexerPPZ.cxx \