+ for (Int_t it=0;it<=t1-t0; it++){
+ x[it]=0;
+ yt[it]=0;
+ zt[it]=0;
+ clusters[it+t0]=-2;
+ zmean[it]=0;
+ nmean[it]=0;
+ //
+ for (Int_t ih=0;ih<10;ih++){
+ indexes[ih][it]=-2; //reset indexes1
+ cl[ih][it]=0;
+ dz[ih][it]=-100;
+ dy[ih][it]=-100;
+ best[ih][it]=0;
+ }
+ }
+ //
+ Double_t x0 = track->GetX();
+ Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
+ Int_t nall=0;
+ Int_t nfound=0;
+ Double_t h01 =0;
+ Int_t plane =-1;
+ Float_t padlength=0;
+ AliTRDtrack track2(*track);
+ Float_t snpy = track->GetSnp();
+ Float_t tany = TMath::Sqrt(snpy*snpy/(1.-snpy*snpy));
+ if (snpy<0) tany*=-1;
+ //
+ Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
+ Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
+ Double_t road = 15.*sqrt(track->GetSigmaY2() + sy2);
+ if (road>6.) road=6.;
+
+ //
+ for (Int_t it=0;it<t1-t0;it++){
+ Double_t maxChi2[2]={fgkMaxChi2,fgkMaxChi2};
+ AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it+t0));
+ if (timeBin==0) continue; // no indexes1
+ Int_t maxn = timeBin;
+ x[it] = timeBin.GetX();
+ track2.PropagateTo(x[it]);
+ yt[it] = track2.GetY();
+ zt[it] = track2.GetZ();
+
+ Double_t y=yt[it],z=zt[it];
+ Double_t chi2 =1000000;
+ nall++;
+ //
+ // find 2 nearest cluster at given time bin
+ //
+ //
+ for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
+ AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
+ h01 = GetTiltFactor(c);
+ if (plane<0){
+ Int_t det = c->GetDetector();
+ plane = fGeom->GetPlane(det);
+ padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
+ }
+ // if (c->GetLocalTimeBin()==0) continue;
+ if (c->GetY() > y+road) break;
+ if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue;
+
+ Double_t dist = TMath::Abs(c->GetZ()-z);
+ if (dist> (0.5*padlength+6.*sigmaz)) continue; // 6 sigma boundary cut
+ Double_t cost = 0;
+ //
+ if (dist> (0.5*padlength-sigmaz)){ // sigma boundary cost function
+ cost = (dist-0.5*padlength)/(2.*sigmaz);
+ if (cost>-1) cost= (cost+1.)*(cost+1.);
+ else cost=0;
+ }
+ // Int_t label = TMath::Abs(track->GetLabel());
+ // if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
+ chi2=track2.GetPredictedChi2(c,h01)+cost;
+ //
+ clfound++;
+ if (chi2 > maxChi2[1]) continue;
+
+ for (Int_t ih=2;ih<9; ih++){ //store the clusters in the road
+ if (cl[ih][it]==0){
+ cl[ih][it] = c;
+ indexes[ih][it] =timeBin.GetIndex(i); // index - 9 - reserved for outliers
+ break;
+ }
+ }
+ //
+ if (chi2 <maxChi2[0]){
+ maxChi2[1] = maxChi2[0];
+ maxChi2[0] = chi2;
+ indexes[1][it] = indexes[0][it];
+ cl[1][it] = cl[0][it];
+ indexes[0][it] = timeBin.GetIndex(i);
+ cl[0][it] = c;
+ continue;
+ }
+ maxChi2[1]=chi2;
+ cl[1][it] = c;
+ indexes[1][it] =timeBin.GetIndex(i);
+ }
+ if (cl[0][it]){
+ nfound++;
+ xmean += x[it];
+ }
+ }
+ //
+ if (nfound<4) return 0;
+ xmean /=Float_t(nfound); // middle x
+ track2.PropagateTo(xmean); // propagate track to the center
+ //
+ // choose one of the variants
+ //
+ Int_t changes[10];
+ Float_t sumz = 0;
+ Float_t sum = 0;
+ Double_t sumdy = 0;
+ Double_t sumdy2 = 0;
+ Double_t sumx = 0;
+ Double_t sumxy = 0;
+ Double_t sumx2 = 0;
+ Double_t mpads = 0;
+ //
+ Int_t ngood[10];
+ Int_t nbad[10];
+ //
+ Double_t meanz[10];
+ Double_t moffset[10]; // mean offset
+ Double_t mean[10]; // mean value
+ Double_t angle[10]; // angle
+ //
+ Double_t smoffset[10]; // sigma of mean offset
+ Double_t smean[10]; // sigma of mean value
+ Double_t sangle[10]; // sigma of angle
+ Double_t smeanangle[10]; // correlation
+ //
+ Double_t sigmas[10];
+ Double_t tchi2s[10]; // chi2s for tracklet
+ //
+ // calculate zmean
+ //
+ for (Int_t it=0;it<t1-t0;it++){
+ if (!cl[0][it]) continue;
+ for (Int_t dt=-3;dt<=3;dt++){
+ if (it+dt<0) continue;
+ if (it+dt>t1-t0) continue;
+ if (!cl[0][it+dt]) continue;
+ zmean[it]+=cl[0][it+dt]->GetZ();
+ nmean[it]+=1.;
+ }
+ zmean[it]/=nmean[it];
+ }
+ //
+ for (Int_t it=0; it<t1-t0;it++){
+ best[0][it]=0;
+ for (Int_t ih=0;ih<10;ih++){
+ dz[ih][it]=-100;
+ dy[ih][it]=-100;
+ if (!cl[ih][it]) continue;
+ Double_t xcluster = cl[ih][it]->GetX();
+ Double_t ytrack,ztrack;
+ track2.GetProlongation(xcluster, ytrack, ztrack );
+ dz[ih][it] = cl[ih][it]->GetZ()- ztrack; // calculate distance from track in z
+ dy[ih][it] = cl[ih][it]->GetY()+ dz[ih][it]*h01 -ytrack; // in y
+ }
+ // minimize changes
+ if (!cl[0][it]) continue;
+ if (TMath::Abs(cl[0][it]->GetZ()-zmean[it])> padlength*0.8 &&cl[1][it])
+ if (TMath::Abs(cl[1][it]->GetZ()-zmean[it])< padlength*0.5){
+ best[0][it]=1;
+ }
+ }
+ //
+ // iterative choosing of "best path"
+ //
+ //
+ Int_t label = TMath::Abs(track->GetLabel());
+ Int_t bestiter=0;
+ //
+ for (Int_t iter=0;iter<9;iter++){
+ //
+ changes[iter]= 0;
+ sumz = 0; sum=0; sumdy=0;sumdy2=0;sumx=0;sumx2=0;sumxy=0;mpads=0; ngood[iter]=0; nbad[iter]=0;
+ // linear fit
+ for (Int_t it=0;it<t1-t0;it++){
+ if (!cl[best[iter][it]][it]) continue;
+ //calculates pad-row changes
+ Double_t zbefore= cl[best[iter][it]][it]->GetZ();
+ Double_t zafter = cl[best[iter][it]][it]->GetZ();
+ for (Int_t itd = it-1; itd>=0;itd--) {
+ if (cl[best[iter][itd]][itd]) {
+ zbefore= cl[best[iter][itd]][itd]->GetZ();
+ break;
+ }
+ }
+ for (Int_t itd = it+1; itd<t1-t0;itd++) {
+ if (cl[best[iter][itd]][itd]) {
+ zafter= cl[best[iter][itd]][itd]->GetZ();
+ break;
+ }
+ }
+ if (TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore)>0.1&&TMath::Abs(cl[best[iter][it]][it]->GetZ()-zafter)>0.1) changes[iter]++;
+ //
+ Double_t dx = x[it]-xmean; // distance to reference x
+ sumz += cl[best[iter][it]][it]->GetZ();
+ sum++;
+ sumdy += dy[best[iter][it]][it];
+ sumdy2+= dy[best[iter][it]][it]*dy[best[iter][it]][it];
+ sumx += dx;
+ sumx2 += dx*dx;
+ sumxy += dx*dy[best[iter][it]][it];
+ mpads += cl[best[iter][it]][it]->GetNPads();
+ if (cl[best[iter][it]][it]->GetLabel(0)==label || cl[best[iter][it]][it]->GetLabel(1)==label||cl[best[iter][it]][it]->GetLabel(2)==label){
+ ngood[iter]++;
+ }
+ else{
+ nbad[iter]++;
+ }
+ }
+ //
+ // calculates line parameters
+ //
+ Double_t det = sum*sumx2-sumx*sumx;
+ angle[iter] = (sum*sumxy-sumx*sumdy)/det;
+ mean[iter] = (sumx2*sumdy-sumx*sumxy)/det;
+ meanz[iter] = sumz/sum;
+ moffset[iter] = sumdy/sum;
+ mpads /= sum; // mean number of pads
+ //
+ //
+ Double_t sigma2 = 0; // normalized residuals - for line fit
+ Double_t sigma1 = 0; // normalized residuals - constant fit
+ //
+ for (Int_t it=0;it<t1-t0;it++){
+ if (!cl[best[iter][it]][it]) continue;
+ Double_t dx = x[it]-xmean;
+ Double_t ytr = mean[iter]+angle[iter]*dx;
+ sigma2 += (dy[best[iter][it]][it]-ytr)*(dy[best[iter][it]][it]-ytr);
+ sigma1 += (dy[best[iter][it]][it]-moffset[iter])*(dy[best[iter][it]][it]-moffset[iter]);
+ sum++;
+ }
+ sigma2 /=(sum-2); // normalized residuals
+ sigma1 /=(sum-1); // normalized residuals
+ //
+ smean[iter] = sigma2*(sumx2/det); // estimated error2 of mean
+ sangle[iter] = sigma2*(sum/det); // estimated error2 of angle
+ smeanangle[iter] = sigma2*(-sumx/det); // correlation
+ //
+ //
+ sigmas[iter] = TMath::Sqrt(sigma1); //
+ smoffset[iter]= (sigma1/sum)+0.01*0.01; // sigma of mean offset + unisochronity sigma
+ //
+ // iterative choosing of "better path"
+ //
+ for (Int_t it=0;it<t1-t0;it++){
+ if (!cl[best[iter][it]][it]) continue;
+ //
+ Double_t sigmatr2 = smoffset[iter]+0.5*tany*tany; //add unisochronity + angular effect contribution
+ Double_t sweight = 1./sigmatr2+1./track->GetSigmaY2();
+ Double_t weighty = (moffset[iter]/sigmatr2)/sweight; // weighted mean
+ Double_t sigmacl = TMath::Sqrt(sigma1*sigma1+track->GetSigmaY2()); //
+ Double_t mindist=100000;
+ Int_t ihbest=0;
+ for (Int_t ih=0;ih<10;ih++){
+ if (!cl[ih][it]) break;
+ Double_t dist2 = (dy[ih][it]-weighty)/sigmacl;
+ dist2*=dist2; //chi2 distance
+ if (dist2<mindist){
+ mindist = dist2;
+ ihbest =ih;
+ }
+ }
+ best[iter+1][it]=ihbest;
+ }
+ //
+ // update best hypothesy if better chi2 according tracklet position and angle
+ //
+ Double_t sy2 = smean[iter] + track->GetSigmaY2();
+ Double_t sa2 = sangle[iter] + track->fCee;
+ Double_t say = track->fCey;
+ // Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2;
+ // Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2;
+
+ Double_t detchi = sy2*sa2-say*say;
+ Double_t invers[3] = {sa2/detchi, sy2/detchi, -say/detchi}; //inverse value of covariance matrix
+
+ Double_t chi20 = mean[bestiter]*mean[bestiter]*invers[0]+angle[bestiter]*angle[bestiter]*invers[1]+
+ 2.*mean[bestiter]*angle[bestiter]*invers[2];
+ Double_t chi21 = mean[iter]*mean[iter]*invers[0]+angle[iter]*angle[iter]*invers[1]+
+ 2*mean[iter]*angle[iter]*invers[2];
+ tchi2s[iter] =chi21;
+ //
+ if (changes[iter]<=changes[bestiter] && chi21<chi20) {
+ bestiter =iter;
+ }
+ }
+ //
+ //set clusters
+ //
+ Double_t sigma2 = sigmas[0]; // choose as sigma from 0 iteration
+ Short_t maxpos = -1;
+ Float_t maxcharge = 0;
+ Short_t maxpos4 = -1;
+ Float_t maxcharge4 = 0;
+ Short_t maxpos5 = -1;
+ Float_t maxcharge5 = 0;
+
+ //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
+ //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept
+
+ Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0));
+ Double_t expectederr = sigma2*sigma2+0.01*0.01;
+ if (mpads>3.5) expectederr += (mpads-3.5)*0.04;
+ if (changes[bestiter]>1) expectederr+= changes[bestiter]*0.01;
+ expectederr+=(0.03*(tany-exB)*(tany-exB))*15;
+ // if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
+ //expectederr+=10000;
+ for (Int_t it=0;it<t1-t0;it++){
+ if (!cl[best[bestiter][it]][it]) continue;
+ cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // set cluster error
+ if (!cl[best[bestiter][it]][it]->IsUsed()){
+ cl[best[bestiter][it]][it]->SetY( cl[best[bestiter][it]][it]->GetY());
+ // cl[best[bestiter][it]][it]->Use();
+ }
+ //
+ // time bins with maximal charge
+ if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
+ maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
+ maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
+ }
+
+ if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
+ if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
+ maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
+ maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
+ }
+ }
+ if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
+ if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
+ maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
+ maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
+ }
+ }
+ //
+ // time bins with maximal charge
+ if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
+ maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
+ maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
+ }
+
+ if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
+ if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
+ maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
+ maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
+ }
+ }
+ if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
+ if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
+ maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
+ maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
+ }
+ }
+ clusters[it+t0] = indexes[best[bestiter][it]][it];
+ //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 && cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0] = indexes[best[bestiter][it]][it]; //Test
+ }
+ //
+ // set tracklet parameters
+ //
+ Double_t trackleterr2 = smoffset[bestiter]+0.01*0.01;
+ if (mpads>3.5) trackleterr2 += (mpads-3.5)*0.04;
+ trackleterr2+= changes[bestiter]*0.01;
+ trackleterr2*= TMath::Max(14.-nfound,1.);
+ trackleterr2+= 0.2*(tany-exB)*(tany-exB);
+ //
+ tracklet.Set(xmean, track2.GetY()+moffset[bestiter], meanz[bestiter], track2.GetAlpha(), trackleterr2); //set tracklet parameters
+ tracklet.SetTilt(h01);
+ tracklet.SetP0(mean[bestiter]);
+ tracklet.SetP1(angle[bestiter]);
+ tracklet.SetN(nfound);
+ tracklet.SetNCross(changes[bestiter]);
+ tracklet.SetPlane(plane);
+ tracklet.SetSigma2(expectederr);
+ tracklet.SetChi2(tchi2s[bestiter]);
+ tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
+ track->fTracklets[plane] = tracklet;
+ track->fNWrong+=nbad[0];
+ //
+ // Debuging part
+ //
+ TClonesArray array0("AliTRDcluster");
+ TClonesArray array1("AliTRDcluster");
+ array0.ExpandCreateFast(t1-t0+1);
+ array1.ExpandCreateFast(t1-t0+1);
+ TTreeSRedirector& cstream = *fDebugStreamer;
+ AliTRDcluster dummy;
+ Double_t dy0[100];
+ Double_t dyb[100];
+
+ for (Int_t it=0;it<t1-t0;it++){
+ dy0[it] = dy[0][it];
+ dyb[it] = dy[best[bestiter][it]][it];
+ if(cl[0][it]) {
+ new(array0[it]) AliTRDcluster(*cl[0][it]);
+ }
+ else{
+ new(array0[it]) AliTRDcluster(dummy);
+ }
+ if(cl[best[bestiter][it]][it]) {
+ new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
+ }
+ else{
+ new(array1[it]) AliTRDcluster(dummy);
+ }
+ }
+ TGraph graph0(t1-t0,x,dy0);
+ TGraph graph1(t1-t0,x,dyb);
+ TGraph graphy(t1-t0,x,yt);
+ TGraph graphz(t1-t0,x,zt);
+ //
+ //
+ cstream<<"tracklet"<<
+ "track.="<<track<< // track parameters
+ "tany="<<tany<< // tangent of the local track angle
+ "xmean="<<xmean<< // xmean - reference x of tracklet
+ "tilt="<<h01<< // tilt angle
+ "nall="<<nall<< // number of foundable clusters
+ "nfound="<<nfound<< // number of found clusters
+ "clfound="<<clfound<< // total number of found clusters in road
+ "mpads="<<mpads<< // mean number of pads per cluster
+ "plane="<<plane<< // plane number
+ "road="<<road<< // the width of the used road
+ "graph0.="<<&graph0<< // x - y = dy for closest cluster
+ "graph1.="<<&graph1<< // x - y = dy for second closest cluster
+ "graphy.="<<&graphy<< // y position of the track
+ "graphz.="<<&graphz<< // z position of the track
+ // "fCl.="<<&array0<< // closest cluster
+ //"fCl2.="<<&array1<< // second closest cluster
+ "maxpos="<<maxpos<< // maximal charge postion
+ "maxcharge="<<maxcharge<< // maximal charge
+ "maxpos4="<<maxpos4<< // maximal charge postion - after bin 4
+ "maxcharge4="<<maxcharge4<< // maximal charge - after bin 4
+ "maxpos5="<<maxpos5<< // maximal charge postion - after bin 5
+ "maxcharge5="<<maxcharge5<< // maximal charge - after bin 5
+ //
+ "bestiter="<<bestiter<< // best iteration number
+ "tracklet.="<<&tracklet<< // corrspond to the best iteration
+ "tchi20="<<tchi2s[0]<< // chi2 of cluster in the 0 iteration
+ "tchi2b="<<tchi2s[bestiter]<< // chi2 of cluster in the best iteration
+ "sigmas0="<<sigmas[0]<< // residuals sigma
+ "sigmasb="<<sigmas[bestiter]<< // residulas sigma
+ //
+ "ngood0="<<ngood[0]<< // number of good clusters in 0 iteration
+ "nbad0="<<nbad[0]<< // number of bad clusters in 0 iteration
+ "ngoodb="<<ngood[bestiter]<< // in best iteration
+ "nbadb="<<nbad[bestiter]<< // in best iteration
+ //
+ "changes0="<<changes[0]<< // changes of pardrows in iteration number 0
+ "changesb="<<changes[bestiter]<< // changes of pardrows in best iteration
+ //
+ "moffset0="<<moffset[0]<< // offset fixing angle in iter=0
+ "smoffset0="<<smoffset[0]<< // sigma of offset fixing angle in iter=0
+ "moffsetb="<<moffset[bestiter]<< // offset fixing angle in iter=best
+ "smoffsetb="<<smoffset[bestiter]<< // sigma of offset fixing angle in iter=best
+ //
+ "mean0="<<mean[0]<< // mean dy in iter=0;
+ "smean0="<<smean[0]<< // sigma of mean dy in iter=0
+ "meanb="<<mean[bestiter]<< // mean dy in iter=best
+ "smeanb="<<smean[bestiter]<< // sigma of mean dy in iter=best
+ //
+ "angle0="<<angle[0]<< // angle deviation in the iteration number 0
+ "sangle0="<<sangle[0]<< // sigma of angular deviation in iteration number 0
+ "angleb="<<angle[bestiter]<< // angle deviation in the best iteration
+ "sangleb="<<sangle[bestiter]<< // sigma of angle deviation in the best iteration
+ //
+ "expectederr="<<expectederr<< // expected error of cluster position
+ "\n";
+ //
+ //
+ return nfound;
+}
+
+
+Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down)
+{
+ //
+ // Sort eleements according occurancy
+ // The size of output array has is 2*n
+ //
+ Int_t * sindexS = new Int_t[n]; // temp array for sorting
+ Int_t * sindexF = new Int_t[2*n];
+ for (Int_t i=0;i<n;i++) sindexF[i]=0;
+ //
+ TMath::Sort(n,inlist, sindexS, down);
+ Int_t last = inlist[sindexS[0]];
+ Int_t val = last;
+ sindexF[0] = 1;
+ sindexF[0+n] = last;
+ Int_t countPos = 0;
+ //
+ // find frequency
+ for(Int_t i=1;i<n; i++){
+ val = inlist[sindexS[i]];
+ if (last == val) sindexF[countPos]++;
+ else{
+ countPos++;
+ sindexF[countPos+n] = val;
+ sindexF[countPos]++;
+ last =val;
+ }
+ }
+ if (last==val) countPos++;
+ // sort according frequency
+ TMath::Sort(countPos, sindexF, sindexS, kTRUE);
+ for (Int_t i=0;i<countPos;i++){
+ outlist[2*i ] = sindexF[sindexS[i]+n];
+ outlist[2*i+1] = sindexF[sindexS[i]];
+ }
+ delete [] sindexS;
+ delete [] sindexF;
+
+ return countPos;
+}
+
+AliTRDtrack * AliTRDtracker::RegisterSeed(AliTRDseed * seeds, Double_t * params)
+{
+ //
+ //
+ //
+ Double_t alpha=AliTRDgeometry::GetAlpha();
+ Double_t shift=AliTRDgeometry::GetAlpha()/2.;
+ Double_t c[15];
+ c[0] = 0.2;
+ c[1] = 0 ; c[2] = 2;
+ c[3] = 0 ; c[4] = 0; c[5] = 0.02;
+ c[6] = 0 ; c[7] = 0; c[8] = 0; c[9] = 0.1;
+ c[10] = 0 ; c[11] = 0; c[12] = 0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
+ //
+ Int_t index =0;
+ AliTRDcluster *cl =0;
+ for (Int_t ilayer=0;ilayer<6;ilayer++){
+ if (seeds[ilayer].isOK()){
+ for (Int_t itime=22;itime>0;itime--){
+ if (seeds[ilayer].fIndexes[itime]>0){
+ index = seeds[ilayer].fIndexes[itime];
+ cl = seeds[ilayer].fClusters[itime];
+ break;
+ }
+ }
+ }
+ if (index>0) break;
+ }
+ if (cl==0) return 0;
+ AliTRDtrack * track = new AliTRDtrack(cl,index,¶ms[1],c, params[0],params[6]*alpha+shift);
+ track->PropagateTo(params[0]-5.);
+ track->ResetCovariance(1);
+ //
+ Int_t rc=FollowBackProlongation(*track);
+ if (rc<30) {
+ delete track;
+ track =0;
+ }else{
+ track->CookdEdx();
+ CookdEdxTimBin(*track);
+ CookLabel(track, 0.9);
+ }
+ return track;
+}
+
+
+
+
+
+
+AliTRDseed::AliTRDseed()
+{
+ //
+ //
+ fTilt =0; // tilting angle
+ fPadLength = 0; // pad length
+ fX0 = 0; // x0 position
+ for (Int_t i=0;i<25;i++){
+ fX[i]=0; // !x position
+ fY[i]=0; // !y position
+ fZ[i]=0; // !z position
+ fIndexes[i]=0; // !indexes
+ fClusters[i]=0; // !clusters
+ }
+ for (Int_t i=0;i<2;i++){
+ fYref[i]=0; // reference y
+ fZref[i]=0; // reference z
+ fYfit[i]=0; // y fit position +derivation
+ fYfitR[i]=0; // y fit position +derivation
+ fZfit[i]=0; // z fit position
+ fZfitR[i]=0; // z fit position
+ fLabels[i]=0; // labels
+ }
+ fSigmaY = 0;
+ fSigmaY2 = 0;
+ fMeanz=0; // mean vaue of z
+ fZProb=0; // max probbable z
+ fMPads=0;
+ //
+ fN=0; // number of associated clusters
+ fN2=0; // number of not crossed
+ fNUsed=0; // number of used clusters
+ fNChange=0; // change z counter
+}
+
+void AliTRDseed::Reset(){
+ //
+ // reset seed
+ //
+ for (Int_t i=0;i<25;i++){
+ fX[i]=0; // !x position
+ fY[i]=0; // !y position
+ fZ[i]=0; // !z position
+ fIndexes[i]=0; // !indexes
+ fClusters[i]=0; // !clusters
+ fUsable[i] = kFALSE;
+ }
+ for (Int_t i=0;i<2;i++){
+ fYref[i]=0; // reference y
+ fZref[i]=0; // reference z
+ fYfit[i]=0; // y fit position +derivation
+ fYfitR[i]=0; // y fit position +derivation
+ fZfit[i]=0; // z fit position
+ fZfitR[i]=0; // z fit position
+ fLabels[i]=-1; // labels
+ }
+ fSigmaY =0; //"robust" sigma in y
+ fSigmaY2=0; //"robust" sigma in y
+ fMeanz =0; // mean vaue of z
+ fZProb =0; // max probbable z
+ fMPads =0;
+ //
+ fN=0; // number of associated clusters
+ fN2=0; // number of not crossed
+ fNUsed=0; // number of used clusters
+ fNChange=0; // change z counter
+}
+
+void AliTRDseed::CookLabels(){
+ //
+ // cook 2 labels for seed
+ //
+ Int_t labels[200];
+ Int_t out[200];
+ Int_t nlab =0;
+ for (Int_t i=0;i<25;i++){
+ if (!fClusters[i]) continue;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ if (fClusters[i]->GetLabel(ilab)>=0){
+ labels[nlab] = fClusters[i]->GetLabel(ilab);
+ nlab++;
+ }
+ }
+ }
+ Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
+ fLabels[0] = out[0];
+ if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
+}
+
+void AliTRDseed::UseClusters()
+{
+ //
+ // use clusters
+ //
+ for (Int_t i=0;i<25;i++){
+ if (!fClusters[i]) continue;
+ if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
+ }
+}
+
+
+void AliTRDseed::Update(){
+ //
+ //
+ //
+ const Float_t ratio = 0.8;
+ const Int_t kClmin = 6;
+ const Float_t kmaxtan = 2;
+ if (TMath::Abs(fYref[1])>kmaxtan) return; // too much inclined track
+ //
+ Float_t sigmaexp = 0.05+TMath::Abs(fYref[1]*0.25); // expected r.m.s in y direction
+ Float_t ycrosscor = fPadLength*fTilt*0.5; // y correction for crossing
+ fNChange =0;
+ //
+ Double_t sumw, sumwx,sumwx2;
+ Double_t sumwy, sumwxy, sumwz,sumwxz;
+ Int_t zints[25]; // histograming of the z coordinate - get 1 and second max probable coodinates in z
+ Int_t zouts[50]; //
+ Float_t allowedz[25]; // allowed z for given time bin
+ Float_t yres[25]; // residuals from reference
+ Float_t anglecor = fTilt*fZref[1]; //correction to the angle
+ //
+ //
+ fN=0; fN2 =0;
+ for (Int_t i=0;i<25;i++){
+ yres[i] =10000;
+ if (!fClusters[i]) continue;
+ yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
+ zints[fN] = Int_t(fZ[i]);
+ fN++;
+ }
+ if (fN<kClmin) return;
+ Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
+ fZProb = zouts[0];
+ if (nz<=1) zouts[3]=0;
+ if (zouts[1]+zouts[3]<kClmin) return;
+ //
+ if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0; // z distance bigger than pad - length
+ //
+ Int_t breaktime = -1;
+ Bool_t mbefore = kFALSE;
+ Int_t cumul[25][2];
+ Int_t counts[2]={0,0};
+ //
+ if (zouts[3]>=3){
+ //
+ // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
+ //
+ fNChange=1;
+ for (Int_t i=0;i<25;i++){
+ cumul[i][0] = counts[0];
+ cumul[i][1] = counts[1];
+ if (TMath::Abs(fZ[i]-zouts[0])<2) counts[0]++;
+ if (TMath::Abs(fZ[i]-zouts[2])<2) counts[1]++;
+ }
+ Int_t maxcount = 0;
+ for (Int_t i=0;i<24;i++) {
+ Int_t after = cumul[24][0]-cumul[i][0];
+ Int_t before = cumul[i][1];
+ if (after+before>maxcount) {
+ maxcount=after+before;
+ breaktime=i;
+ mbefore=kFALSE;
+ }
+ after = cumul[24][1]-cumul[i][1];
+ before = cumul[i][0];
+ if (after+before>maxcount) {
+ maxcount=after+before;
+ breaktime=i;
+ mbefore=kTRUE;
+ }
+ }
+ breaktime-=1;
+ }
+ for (Int_t i=0;i<25;i++){
+ if (i>breaktime) allowedz[i] = mbefore ? zouts[2]:zouts[0];
+ if (i<=breaktime) allowedz[i] = (!mbefore) ? zouts[2]:zouts[0];
+ }
+ if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] && fZref[1]>0)){
+ //
+ // tracklet z-direction not in correspondance with track z direction
+ //
+ fNChange =0;
+ for (Int_t i=0;i<25;i++){
+ allowedz[i] = zouts[0]; //only longest taken
+ }
+ }
+ //
+ if (fNChange>0){
+ //
+ // cross pad -row tracklet - take the step change into account
+ //
+ for (Int_t i=0;i<25;i++){
+ if (!fClusters[i]) continue;
+ if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
+ yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
+ if (TMath::Abs(fZ[i]-fZProb)>2){
+ if (fZ[i]>fZProb) yres[i]+=fTilt*fPadLength;
+ if (fZ[i]<fZProb) yres[i]-=fTilt*fPadLength;
+ }
+ }
+ }
+ //
+ Double_t yres2[25];
+ Double_t mean,sigma;
+ for (Int_t i=0;i<25;i++){
+ if (!fClusters[i]) continue;
+ if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
+ yres2[fN2] = yres[i];
+ fN2++;
+ }
+ if (fN2<kClmin){
+ fN2 = 0;
+ return;
+ }
+ EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*ratio-2));
+ if (sigma<sigmaexp*0.8) sigma=sigmaexp;
+ fSigmaY = sigma;
+ //
+ //
+ // reset sums
+ sumw=0; sumwx=0; sumwx2=0;
+ sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
+ fN2 =0;
+ fMeanz =0;
+ fMPads =0;
+ //
+ for (Int_t i=0;i<25;i++){
+ fUsable[i]=kFALSE;
+ if (!fClusters[i]) continue;
+ if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
+ if (TMath::Abs(yres[i]-mean)>4.*sigma) continue;
+ fUsable[i] = kTRUE;
+ fN2++;
+ fMPads+=fClusters[i]->GetNPads();
+ Float_t weight =1;
+ if (fClusters[i]->GetNPads()>4) weight=0.5;
+ if (fClusters[i]->GetNPads()>5) weight=0.2;
+ //
+ Double_t x = fX[i];
+ sumw+=weight; sumwx+=x*weight; sumwx2+=x*x*weight;
+ sumwy+=weight*yres[i]; sumwxy+=weight*(yres[i])*x;
+ sumwz+=weight*fZ[i]; sumwxz+=weight*fZ[i]*x;
+ }
+ if (fN2<kClmin){
+ fN2 = 0;
+ return;
+ }
+ fMeanz = sumwz/sumw;
+ Float_t correction =0;
+ if (fNChange>0){
+ // tracklet on boundary
+ if (fMeanz<fZProb) correction = ycrosscor;
+ if (fMeanz>fZProb) correction = -ycrosscor;
+ }
+ Double_t det = sumw*sumwx2-sumwx*sumwx;
+ fYfitR[0] = (sumwx2*sumwy-sumwx*sumwxy)/det;
+ fYfitR[1] = (sumw*sumwxy-sumwx*sumwy)/det;
+ //
+ fSigmaY2 =0;
+ for (Int_t i=0;i<25;i++){
+ if (!fUsable[i]) continue;
+ Float_t delta = yres[i]-fYfitR[0]-fYfitR[1]*fX[i];
+ fSigmaY2+=delta*delta;
+ }
+ fSigmaY2 = TMath::Sqrt(fSigmaY2/Float_t(fN2-2));
+ //
+ fZfitR[0] = (sumwx2*sumwz-sumwx*sumwxz)/det;
+ fZfitR[1] = (sumw*sumwxz-sumwx*sumwz)/det;
+ fZfit[0] = (sumwx2*sumwz-sumwx*sumwxz)/det;
+ fZfit[1] = (sumw*sumwxz-sumwx*sumwz)/det;
+ fYfitR[0] += fYref[0]+correction;
+ fYfitR[1] += fYref[1];
+ fYfit[0] = fYfitR[0];
+ fYfit[1] = fYfitR[1];
+ //
+ //
+ UpdateUsed();
+}
+
+
+
+
+
+
+void AliTRDseed::UpdateUsed(){
+ //
+ fNUsed =0;
+ for (Int_t i=0;i<25;i++){
+ if (!fClusters[i]) continue;
+ if ((fClusters[i]->IsUsed())) fNUsed++;
+ }
+}
+
+
+void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh)
+{
+ //
+ // robust estimator in 1D case MI version
+ //
+ //for the univariate case
+ //estimates of location and scatter are returned in mean and sigma parameters
+ //the algorithm works on the same principle as in multivariate case -
+ //it finds a subset of size hh with smallest sigma, and then returns mean and
+ //sigma of this subset
+
+ if (hh==0)
+ hh=(nvectors+2)/2;
+ Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
+ Int_t *index=new Int_t[nvectors];
+ TMath::Sort(nvectors, data, index, kFALSE);
+ //
+ Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
+ Double_t factor = faclts[nquant-1];
+ //
+ //
+ Double_t sumx =0;
+ Double_t sumx2 =0;
+ Int_t bestindex = -1;
+ Double_t bestmean = 0;
+ Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma
+ for (Int_t i=0; i<hh; i++){
+ sumx += data[index[i]];
+ sumx2 += data[index[i]]*data[index[i]];
+ }
+ //
+ Double_t norm = 1./Double_t(hh);
+ Double_t norm2 = 1./Double_t(hh-1);
+ for (Int_t i=hh; i<nvectors; i++){
+ Double_t cmean = sumx*norm;
+ Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
+ if (csigma<bestsigma){
+ bestmean = cmean;
+ bestsigma = csigma;
+ bestindex = i-hh;
+ }
+ //
+ //
+ sumx += data[index[i]]-data[index[i-hh]];
+ sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
+ }
+
+ Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
+ mean = bestmean;
+ sigma = bstd;
+ delete [] index;
+}
+
+
+Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror){
+ //
+ //
+ //
+ TLinearFitter fitterT2(4,"hyp4"); // fitting with tilting pads - kz not fixed
+ fitterT2.StoreData(kTRUE);
+ Float_t xref2 = (cseed[2].fX0+cseed[3].fX0)*0.5; // reference x0 for z
+ //
+ Int_t npointsT =0;
+ fitterT2.ClearPoints();
+ for (Int_t iLayer=0; iLayer<6;iLayer++){
+ if (!cseed[iLayer].isOK()) continue;
+ Double_t tilt = cseed[iLayer].fTilt;
+
+ for (Int_t itime=0;itime<25;itime++){
+ if (!cseed[iLayer].fUsable[itime]) continue;
+ Double_t x = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2; // x relative to the midle chamber
+ Double_t y = cseed[iLayer].fY[itime];
+ Double_t z = cseed[iLayer].fZ[itime];
+ // tilted rieman
+ //
+ Double_t uvt[6];
+ Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0; // global x
+ Double_t t = 1./(x2*x2+y*y);
+ uvt[1] = t; // t
+ uvt[0] = 2.*x2*uvt[1]; // u
+ uvt[2] = 2.0*tilt*uvt[1];
+ uvt[3] = 2.0*tilt*x*uvt[1];
+ uvt[4] = 2.0*(y+tilt*z)*uvt[1];
+ //
+ Double_t error = 2*uvt[1];
+ if (terror) error*=cseed[iLayer].fSigmaY;
+ else {error *=0.2;} //default error
+ fitterT2.AddPoint(uvt,uvt[4],error);
+ npointsT++;
+ }
+ }
+ fitterT2.Eval();
+ Double_t rpolz0 = fitterT2.GetParameter(3);
+ Double_t rpolz1 = fitterT2.GetParameter(4);
+ //
+ // linear fitter - not possible to make boundaries
+ // non accept non possible z and dzdx combination
+ //
+ Bool_t acceptablez =kTRUE;
+ for (Int_t iLayer=0; iLayer<6;iLayer++){
+ if (cseed[iLayer].isOK()){
+ Double_t zT2 = rpolz0+rpolz1*(cseed[iLayer].fX0 - xref2);
+ if (TMath::Abs(cseed[iLayer].fZProb-zT2)>cseed[iLayer].fPadLength*0.5+1)
+ acceptablez = kFALSE;
+ }
+ }
+ if (!acceptablez){
+ Double_t zmf = cseed[2].fZref[0]+cseed[2].fZref[1]*(xref2-cseed[2].fX0);
+ Double_t dzmf = (cseed[2].fZref[1]+ cseed[3].fZref[1])*0.5;
+ fitterT2.FixParameter(3,zmf);
+ fitterT2.FixParameter(4,dzmf);
+ fitterT2.Eval();
+ fitterT2.ReleaseParameter(3);
+ fitterT2.ReleaseParameter(4);
+ rpolz0 = fitterT2.GetParameter(3);
+ rpolz1 = fitterT2.GetParameter(4);
+ }
+ //
+ Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
+ Double_t params[3];
+ params[0] = fitterT2.GetParameter(0);
+ params[1] = fitterT2.GetParameter(1);
+ params[2] = fitterT2.GetParameter(2);
+ Double_t CR = 1+params[1]*params[1]-params[2]*params[0];
+ for (Int_t iLayer = 0; iLayer<6;iLayer++){
+ Double_t x = cseed[iLayer].fX0;
+ Double_t y=0,dy=0, z=0, dz=0;
+ // y
+ Double_t res2 = (x*params[0]+params[1]);
+ res2*=res2;
+ res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
+ if (res2>=0){
+ res2 = TMath::Sqrt(res2);
+ y = (1-res2)/params[0];
+ }
+ //dy
+ Double_t x0 = -params[1]/params[0];
+ if (-params[2]*params[0]+params[1]*params[1]+1>0){
+ Double_t Rm1 = params[0]/TMath::Sqrt(-params[2]*params[0]+params[1]*params[1]+1);
+ if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)>0){
+ Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
+ if (params[0]<0) res*=-1.;
+ dy = res;
+ }
+ }
+ z = rpolz0+rpolz1*(x-xref2);
+ dz = rpolz1;
+ cseed[iLayer].fYref[0] = y;
+ cseed[iLayer].fYref[1] = dy;
+ cseed[iLayer].fZref[0] = z;
+ cseed[iLayer].fZref[1] = dz;
+ cseed[iLayer].fC = CR;
+ //
+ }
+ return chi2TR;
+}