+ //
+ // linear fitters in planes
+ TLinearFitter fitterTC(2,"hyp2"); // fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
+ TLinearFitter fitterT2(4,"hyp4"); // fitting with tilting pads - kz not fixed
+ fitterTC.StoreData(kTRUE);
+ fitterT2.StoreData(kTRUE);
+ AliRieman rieman(1000); // rieman fitter
+ AliRieman rieman2(1000); // rieman fitter
+ //
+ // find the maximal and minimal layer for the planes
+ //
+ Int_t layers[6][2];
+ AliTRDpropagationLayer* reflayers[6];
+ for (Int_t i=0;i<6;i++){layers[i][0]=10000; layers[i][1]=0;}
+ for (Int_t ns=0;ns<maxSec;ns++){
+ for (Int_t ilayer=0;ilayer<fTrSec[ns]->GetNumberOfLayers();ilayer++){
+ AliTRDpropagationLayer& layer=*(fTrSec[ns]->GetLayer(ilayer));
+ if (layer==0) continue;
+ Int_t det = layer[0]->GetDetector();
+ Int_t plane = fGeom->GetPlane(det);
+ if (ilayer<layers[plane][0]) layers[plane][0] = ilayer;
+ if (ilayer>layers[plane][1]) layers[plane][1] = ilayer;
+ }
+ }
+ //
+ AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
+ Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
+ Double_t hL[6]; // tilting angle
+ Double_t xcl[6]; // x - position of reference cluster
+ Double_t ycl[6]; // y - position of reference cluster
+ Double_t zcl[6]; // z - position of reference cluster
+ AliTRDcluster *cl[6]={0,0,0,0,0,0}; // seeding clusters
+ Float_t padlength[6]={10,10,10,10,10,10}; //current pad-length
+ Double_t chi2R =0, chi2Z=0;
+ Double_t chi2RF =0, chi2ZF=0;
+ //
+ Int_t nclusters; // total number of clusters
+ for (Int_t i=0;i<6;i++) {hL[i]=h01; if (i%2==1) hL[i]*=-1.;}
+ //
+ //
+ // registered seed
+ AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
+ AliTRDseed *seed[kMaxSeed];
+ for (Int_t iseed=0;iseed<kMaxSeed;iseed++) seed[iseed]= &pseed[iseed*6];
+ AliTRDseed *cseed = seed[0];
+ //
+ Double_t seedquality[kMaxSeed];
+ Double_t seedquality2[kMaxSeed];
+ Double_t seedparams[kMaxSeed][7];
+ Int_t seedlayer[kMaxSeed];
+ Int_t registered =0;
+ Int_t sort[kMaxSeed];
+ //
+ // seeding part
+ //
+ for (Int_t ns = 0; ns<maxSec; ns++){ //loop over sectors
+ //for (Int_t ns = 0; ns<5; ns++){ //loop over sectors
+ registered = 0; // reset registerd seed counter
+ cseed = seed[registered];
+ Float_t iter=0;
+ for (Int_t sLayer=2; sLayer>=0;sLayer--){
+ //for (Int_t dseed=5;dseed<15; dseed+=3){ //loop over central seeding time bins
+ iter+=1.;
+ Int_t dseed = 5+Int_t(iter)*3;
+ // Initialize seeding layers
+ for (Int_t ilayer=0;ilayer<6;ilayer++){
+ reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
+ xcl[ilayer] = reflayers[ilayer]->GetX();
+ }
+ //
+ Double_t xref = (xcl[sLayer+1] + xcl[sLayer+2])*0.5;
+ AliTRDpropagationLayer& layer0=*reflayers[sLayer+0];
+ AliTRDpropagationLayer& layer1=*reflayers[sLayer+1];
+ AliTRDpropagationLayer& layer2=*reflayers[sLayer+2];
+ AliTRDpropagationLayer& layer3=*reflayers[sLayer+3];
+ //
+ Int_t maxn3 = layer3;
+ for (Int_t icl3=0;icl3<maxn3;icl3++){
+ AliTRDcluster *cl3 = layer3[icl3];
+ if (!cl3) continue;
+ padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2()*12.);
+ ycl[sLayer+3] = cl3->GetY();
+ zcl[sLayer+3] = cl3->GetZ();
+ Float_t yymin0 = ycl[sLayer+3] - 1- kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
+ Float_t yymax0 = ycl[sLayer+3] + 1+ kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
+ Int_t maxn0 = layer0; //
+ for (Int_t icl0=layer0.Find(yymin0);icl0<maxn0;icl0++){
+ AliTRDcluster *cl0 = layer0[icl0];
+ if (!cl0) continue;
+ if (cl3->IsUsed()&&cl0->IsUsed()) continue;
+ ycl[sLayer+0] = cl0->GetY();
+ zcl[sLayer+0] = cl0->GetZ();
+ if ( ycl[sLayer+0]>yymax0) break;
+ Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
+ if (TMath::Abs(tanphi)>kMaxPhi) continue;
+ Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
+ if (TMath::Abs(tantheta)>kMaxTheta) continue;
+ padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2()*12.);
+ //
+ // expected position in 1 layer
+ Double_t y1exp = ycl[sLayer+0]+(tanphi) *(xcl[sLayer+1]-xcl[sLayer+0]);
+ Double_t z1exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+1]-xcl[sLayer+0]);
+ Float_t yymin1 = y1exp - kRoad0y-tanphi;
+ Float_t yymax1 = y1exp + kRoad0y+tanphi;
+ Int_t maxn1 = layer1; //
+ //
+ for (Int_t icl1=layer1.Find(yymin1);icl1<maxn1;icl1++){
+ AliTRDcluster *cl1 = layer1[icl1];
+ if (!cl1) continue;
+ Int_t nusedCl = 0;
+ if (cl3->IsUsed()) nusedCl++;
+ if (cl0->IsUsed()) nusedCl++;
+ if (cl1->IsUsed()) nusedCl++;
+ if (nusedCl>1) continue;
+ ycl[sLayer+1] = cl1->GetY();
+ zcl[sLayer+1] = cl1->GetZ();
+ if ( ycl[sLayer+1]>yymax1) break;
+ if (TMath::Abs(ycl[sLayer+1]-y1exp)>kRoad0y+tanphi) continue;
+ if (TMath::Abs(zcl[sLayer+1]-z1exp)>kRoad0z) continue;
+ padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2()*12.);
+ //
+ Double_t y2exp = ycl[sLayer+0]+(tanphi) *(xcl[sLayer+2]-xcl[sLayer+0])+(ycl[sLayer+1]-y1exp);
+ Double_t z2exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+2]-xcl[sLayer+0]);
+ Int_t index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y, kRoad1z);
+ if (index2<=0) continue;
+ AliTRDcluster *cl2 = (AliTRDcluster*)GetCluster(index2);
+ padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2()*12.);
+ ycl[sLayer+2] = cl2->GetY();
+ zcl[sLayer+2] = cl2->GetZ();
+ if (TMath::Abs(cl2->GetZ()-z2exp)>kRoad0z) continue;
+ //
+ rieman.Reset();
+ rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
+ rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
+ rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
+ rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
+ rieman.Update();
+ //
+ // reset fitter
+ for (Int_t iLayer=0;iLayer<6;iLayer++){
+ cseed[iLayer].Reset();
+ }
+ chi2Z =0.; chi2R=0.;
+ for (Int_t iLayer=0;iLayer<4;iLayer++){
+ cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
+ chi2Z += (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer])*
+ (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer]);
+ cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
+ cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
+ chi2R += (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer])*
+ (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer]);
+ cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
+ }
+ if (TMath::Sqrt(chi2R)>1./iter) continue;
+ if (TMath::Sqrt(chi2Z)>7./iter) continue;
+ //
+ //
+ //
+ Float_t minmax[2]={-100,100};
+ for (Int_t iLayer=0;iLayer<4;iLayer++){
+ Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer]*0.5+1 -cseed[sLayer+iLayer].fZref[0];
+ if (max<minmax[1]) minmax[1]=max;
+ Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer]*0.5-1 -cseed[sLayer+iLayer].fZref[0];
+ if (min>minmax[0]) minmax[0]=min;
+ }
+ Bool_t isFake = kFALSE;
+ if (cl0->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
+ if (cl1->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
+ if (cl2->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
+ if (AliTRDReconstructor::StreamLevel()>0){
+ if ((!isFake) || (icl3%10)==0 ){ //debugging print
+ TTreeSRedirector& cstream = *fDebugStreamer;
+ cstream<<"Seeds0"<<
+ "isFake="<<isFake<<
+ "Cl0.="<<cl0<<
+ "Cl1.="<<cl1<<
+ "Cl2.="<<cl2<<
+ "Cl3.="<<cl3<<
+ "Xref="<<xref<<
+ "X0="<<xcl[sLayer+0]<<
+ "X1="<<xcl[sLayer+1]<<
+ "X2="<<xcl[sLayer+2]<<
+ "X3="<<xcl[sLayer+3]<<
+ "Y2exp="<<y2exp<<
+ "Z2exp="<<z2exp<<
+ "Chi2R="<<chi2R<<
+ "Chi2Z="<<chi2Z<<
+ "Seed0.="<<&cseed[sLayer+0]<<
+ "Seed1.="<<&cseed[sLayer+1]<<
+ "Seed2.="<<&cseed[sLayer+2]<<
+ "Seed3.="<<&cseed[sLayer+3]<<
+ "Zmin="<<minmax[0]<<
+ "Zmax="<<minmax[1]<<
+ "\n";
+ }
+ }
+
+ //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ //<<<<<<<<<<<<<<<<<< FIT SEEDING PART <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ cl[sLayer+0] = cl0;
+ cl[sLayer+1] = cl1;
+ cl[sLayer+2] = cl2;
+ cl[sLayer+3] = cl3;
+ Bool_t isOK=kTRUE;
+ for (Int_t jLayer=0;jLayer<4;jLayer++){
+ cseed[sLayer+jLayer].fTilt = hL[sLayer+jLayer];
+ cseed[sLayer+jLayer].fPadLength = padlength[sLayer+jLayer];
+ cseed[sLayer+jLayer].fX0 = xcl[sLayer+jLayer];
+ for (Int_t iter=0; iter<2; iter++){
+ //
+ // in iteration 0 we try only one pad-row
+ // if quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
+ //
+ AliTRDseed tseed = cseed[sLayer+jLayer];
+ Float_t roadz = padlength[sLayer+jLayer]*0.5;
+ if (iter>0) roadz = padlength[sLayer+jLayer];
+ //
+ Float_t quality =10000;
+ for (Int_t iTime=2;iTime<20;iTime++){
+ AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
+ Double_t dxlayer= layer.GetX()-xcl[sLayer+jLayer];
+ Double_t zexp = cl[sLayer+jLayer]->GetZ() ;
+ if (iter>0){
+ // try 2 pad-rows in second iteration
+ zexp = tseed.fZref[0]+ tseed.fZref[1]*dxlayer;
+ if (zexp>cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()+padlength[sLayer+jLayer]*0.5;
+ if (zexp<cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()-padlength[sLayer+jLayer]*0.5;
+ }
+ //
+ Double_t yexp = tseed.fYref[0]+
+ tseed.fYref[1]*dxlayer;
+ Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
+ if (index<=0) continue;
+ AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
+ //
+ tseed.fIndexes[iTime] = index;
+ tseed.fClusters[iTime] = cl; // register cluster
+ tseed.fX[iTime] = dxlayer; // register cluster
+ tseed.fY[iTime] = cl->GetY(); // register cluster
+ tseed.fZ[iTime] = cl->GetZ(); // register cluster
+ }
+ tseed.Update();
+ //count the number of clusters and distortions into quality
+ Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
+ Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
+ TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
+ 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
+ if (iter==0 && tseed.IsOK()) {
+ cseed[sLayer+jLayer] = tseed;
+ quality = tquality;
+ if (tquality<5) break;
+ }
+ if (tseed.IsOK() && tquality<quality)
+ cseed[sLayer+jLayer] = tseed;
+ }
+ if (!cseed[sLayer+jLayer].IsOK()){
+ isOK = kFALSE;
+ break;
+ }
+ cseed[sLayer+jLayer].CookLabels();
+ cseed[sLayer+jLayer].UpdateUsed();
+ nusedCl+= cseed[sLayer+jLayer].fNUsed;
+ if (nusedCl>25){
+ isOK = kFALSE;
+ break;
+ }
+ }
+ //
+ if (!isOK) continue;
+ nclusters=0;
+ for (Int_t iLayer=0;iLayer<4;iLayer++){
+ if (cseed[sLayer+iLayer].IsOK()){
+ nclusters+=cseed[sLayer+iLayer].fN2;
+ }
+ }
+ //
+ // iteration 0
+ rieman.Reset();
+ for (Int_t iLayer=0;iLayer<4;iLayer++){
+ rieman.AddPoint(xcl[sLayer+iLayer],cseed[sLayer+iLayer].fYfitR[0],
+ cseed[sLayer+iLayer].fZProb,1,10);
+ }
+ rieman.Update();
+ //
+ //
+ chi2R =0; chi2Z=0;
+ for (Int_t iLayer=0;iLayer<4;iLayer++){
+ cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
+ chi2R += (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0])*
+ (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0]);
+ cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
+ cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
+ chi2Z += (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz)*
+ (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz);
+ cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
+ }
+ Double_t curv = rieman.GetC();
+ //
+ // likelihoods
+ //
+ Double_t sumda =
+ TMath::Abs(cseed[sLayer+0].fYfitR[1]- cseed[sLayer+0].fYref[1])+
+ TMath::Abs(cseed[sLayer+1].fYfitR[1]- cseed[sLayer+1].fYref[1])+
+ TMath::Abs(cseed[sLayer+2].fYfitR[1]- cseed[sLayer+2].fYref[1])+
+ TMath::Abs(cseed[sLayer+3].fYfitR[1]- cseed[sLayer+3].fYref[1]);
+ Double_t likea = TMath::Exp(-sumda*10.6);
+ Double_t likechi2 = 0.0000000001;
+ if (chi2R<0.5) likechi2+=TMath::Exp(-TMath::Sqrt(chi2R)*7.73);
+ Double_t likechi2z = TMath::Exp(-chi2Z*0.088)/TMath::Exp(-chi2Z*0.019);
+ Double_t likeN = TMath::Exp(-(72-nclusters)*0.19);
+ Double_t like = likea*likechi2*likechi2z*likeN;
+ //
+ Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fYref[1]-130*curv)*1.9);
+ Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fZref[1]-
+ cseed[sLayer+0].fZref[0]/xcl[sLayer+0])*5.9);
+ Double_t likePrim = TMath::Max(likePrimY*likePrimZ,0.0005);
+
+ seedquality[registered] = like;
+ seedlayer[registered] = sLayer;
+ if (TMath::Log(0.000000000000001+like)<-15) continue;
+ AliTRDseed seedb[6];
+ for (Int_t iLayer=0;iLayer<6;iLayer++){
+ seedb[iLayer] = cseed[iLayer];
+ }
+ //
+ //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ //<<<<<<<<<<<<<<< FULL TRACK FIT PART <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+ //
+ Int_t nlayers = 0;
+ Int_t nusedf = 0;
+ Int_t findable = 0;
+ //
+ // add new layers - avoid long extrapolation
+ //
+ Int_t tLayer[2]={0,0};
+ if (sLayer==2) {tLayer[0]=1; tLayer[1]=0;}
+ if (sLayer==1) {tLayer[0]=5; tLayer[1]=0;}
+ if (sLayer==0) {tLayer[0]=4; tLayer[1]=5;}
+ //
+ for (Int_t iLayer=0;iLayer<2;iLayer++){
+ Int_t jLayer = tLayer[iLayer]; // set tracking layer
+ cseed[jLayer].Reset();
+ cseed[jLayer].fTilt = hL[jLayer];
+ cseed[jLayer].fPadLength = padlength[jLayer];
+ cseed[jLayer].fX0 = xcl[jLayer];
+ // get pad length and rough cluster
+ Int_t indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].fYref[0],
+ cseed[jLayer].fZref[0],kRoad2y,kRoad2z);
+ if (indexdummy<=0) continue;
+ AliTRDcluster *cldummy = (AliTRDcluster*)GetCluster(indexdummy);
+ padlength[jLayer] = TMath::Sqrt(cldummy->GetSigmaZ2()*12.);
+ }
+ AliTRDseed::FitRiemanTilt(cseed, kTRUE);
+ //
+ for (Int_t iLayer=0;iLayer<2;iLayer++){
+ Int_t jLayer = tLayer[iLayer]; // set tracking layer
+ if ( (jLayer==0) && !(cseed[1].IsOK())) continue; // break not allowed
+ if ( (jLayer==5) && !(cseed[4].IsOK())) continue; // break not allowed
+ Float_t zexp = cseed[jLayer].fZref[0];
+ Double_t zroad = padlength[jLayer]*0.5+1.;
+ //
+ //
+ for (Int_t iter=0;iter<2;iter++){
+ AliTRDseed tseed = cseed[jLayer];
+ Float_t quality = 10000;
+ for (Int_t iTime=2;iTime<20;iTime++){
+ AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
+ Double_t dxlayer = layer.GetX()-xcl[jLayer];
+ Double_t yexp = tseed.fYref[0]+tseed.fYref[1]*dxlayer;
+ Float_t yroad = kRoad1y;
+ Int_t index = layer.FindNearestCluster(yexp,zexp, yroad, zroad);
+ if (index<=0) continue;
+ AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
+ //
+ tseed.fIndexes[iTime] = index;
+ tseed.fClusters[iTime] = cl; // register cluster
+ tseed.fX[iTime] = dxlayer; // register cluster
+ tseed.fY[iTime] = cl->GetY(); // register cluster
+ tseed.fZ[iTime] = cl->GetZ(); // register cluster
+ }
+ tseed.Update();
+ if (tseed.IsOK()){
+ Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
+ Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
+ TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
+ 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
+ //
+ if (tquality<quality){
+ cseed[jLayer]=tseed;
+ quality = tquality;
+ }
+ }
+ zroad*=2.;
+ }
+ if ( cseed[jLayer].IsOK()){
+ cseed[jLayer].CookLabels();
+ cseed[jLayer].UpdateUsed();
+ nusedf+= cseed[jLayer].fNUsed;
+ AliTRDseed::FitRiemanTilt(cseed, kTRUE);
+ }
+ }
+ //
+ //
+ // make copy
+ AliTRDseed bseed[6];
+ for (Int_t jLayer=0;jLayer<6;jLayer++){
+ bseed[jLayer] = cseed[jLayer];
+ }
+ Float_t lastquality = 10000;
+ Float_t lastchi2 = 10000;
+ Float_t chi2 = 1000;
+
+ //
+ for (Int_t iter =0; iter<4;iter++){
+ //
+ // sort tracklets according "quality", try to "improve" 4 worst
+ //
+ Float_t sumquality = 0;
+ Float_t squality[6];
+ Int_t sortindexes[6];
+ for (Int_t jLayer=0;jLayer<6;jLayer++){
+ if (bseed[jLayer].IsOK()){
+ AliTRDseed &tseed = bseed[jLayer];
+ Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
+ Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
+ Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
+ TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
+ 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
+ squality[jLayer] = tquality;
+ }
+ else squality[jLayer]=-1;
+ sumquality +=squality[jLayer];
+ }
+
+ if (sumquality>=lastquality || chi2>lastchi2) break;
+ lastquality = sumquality;
+ lastchi2 = chi2;
+ if (iter>0){
+ for (Int_t jLayer=0;jLayer<6;jLayer++){
+ cseed[jLayer] = bseed[jLayer];
+ }
+ }
+ TMath::Sort(6,squality,sortindexes,kFALSE);
+ //
+ //
+ for (Int_t jLayer=5;jLayer>1;jLayer--){
+ Int_t bLayer = sortindexes[jLayer];
+ AliTRDseed tseed = bseed[bLayer];
+ for (Int_t iTime=2;iTime<20;iTime++){
+ AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
+ Double_t dxlayer= layer.GetX()-xcl[bLayer];
+ //
+ Double_t zexp = tseed.fZref[0];
+ Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
+ //
+ Float_t roadz = padlength[bLayer]+1;
+ if (TMath::Abs(tseed.fZProb-zexp)> padlength[bLayer]*0.5) {roadz = padlength[bLayer]*0.5;}
+ if (tseed.fZfit[1]*tseed.fZref[1]<0) {roadz = padlength[bLayer]*0.5;}
+ if (TMath::Abs(tseed.fZProb-zexp)<0.1*padlength[bLayer]) {
+ zexp = tseed.fZProb;
+ roadz = padlength[bLayer]*0.5;
+ }
+ //
+ Double_t yexp = tseed.fYref[0]+
+ tseed.fYref[1]*dxlayer-zcor;
+ Int_t index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
+ if (index<=0) continue;
+ AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
+ //
+ tseed.fIndexes[iTime] = index;
+ tseed.fClusters[iTime] = cl; // register cluster
+ tseed.fX[iTime] = dxlayer; // register cluster
+ tseed.fY[iTime] = cl->GetY(); // register cluster
+ tseed.fZ[iTime] = cl->GetZ(); // register cluster
+ }
+ tseed.Update();
+ if (tseed.IsOK()) {
+ Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
+ Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
+ //
+ Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
+ TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
+ 2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
+ //
+ if (tquality<squality[bLayer])
+ bseed[bLayer] = tseed;
+ }
+ }
+ chi2 = AliTRDseed::FitRiemanTilt(bseed, kTRUE);
+ }
+ //
+ //
+ //
+ nclusters = 0;
+ nlayers = 0;
+ findable = 0;
+ for (Int_t iLayer=0;iLayer<6;iLayer++) {
+ if (TMath::Abs(cseed[iLayer].fYref[0]/cseed[iLayer].fX0)<0.15)
+ findable++;
+ if (cseed[iLayer].IsOK()){
+ nclusters+=cseed[iLayer].fN2;
+ nlayers++;
+ }
+ }
+ if (nlayers<3) continue;
+ rieman.Reset();
+ for (Int_t iLayer=0;iLayer<6;iLayer++){
+ if (cseed[iLayer].IsOK()) rieman.AddPoint(xcl[iLayer],cseed[iLayer].fYfitR[0],
+ cseed[iLayer].fZProb,1,10);
+ }
+ rieman.Update();
+ //
+ chi2RF =0;
+ chi2ZF =0;
+ for (Int_t iLayer=0;iLayer<6;iLayer++){
+ if (cseed[iLayer].IsOK()){
+ cseed[iLayer].fYref[0] = rieman.GetYat(xcl[iLayer]);
+ chi2RF += (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0])*
+ (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0]);
+ cseed[iLayer].fYref[1] = rieman.GetDYat(xcl[iLayer]);
+ cseed[iLayer].fZref[0] = rieman.GetZat(xcl[iLayer]);
+ chi2ZF += (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz)*
+ (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz);
+ cseed[iLayer].fZref[1] = rieman.GetDZat(xcl[iLayer]);
+ }
+ }
+ chi2RF/=TMath::Max((nlayers-3.),1.);
+ chi2ZF/=TMath::Max((nlayers-3.),1.);
+ curv = rieman.GetC();
+
+ //
+
+ Double_t xref2 = (xcl[2]+xcl[3])*0.5; // middle of the chamber
+ Double_t dzmf = rieman.GetDZat(xref2);
+ Double_t zmf = rieman.GetZat(xref2);
+ //
+ // fit hyperplane
+ //
+ Int_t npointsT =0;
+ fitterTC.ClearPoints();
+ fitterT2.ClearPoints();
+ rieman2.Reset();
+ for (Int_t iLayer=0; iLayer<6;iLayer++){
+ if (!cseed[iLayer].IsOK()) continue;
+ 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];
+ // ExB correction to the correction
+ // 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*hL[iLayer]*uvt[1];
+ uvt[3] = 2.0*hL[iLayer]*x*uvt[1];
+ uvt[4] = 2.0*(y+hL[iLayer]*z)*uvt[1];
+ //
+ Double_t error = 2*0.2*uvt[1];
+ fitterT2.AddPoint(uvt,uvt[4],error);
+ //
+ // constrained rieman
+ //
+ z =cseed[iLayer].fZ[itime];
+ uvt[0] = 2.*x2*t; // u
+ uvt[1] = 2*hL[iLayer]*x2*uvt[1];
+ uvt[2] = 2*(y+hL[iLayer]*(z-GetZ()))*t;
+ fitterTC.AddPoint(uvt,uvt[2],error);
+ //
+ rieman2.AddPoint(x2,y,z,1,10);
+ npointsT++;
+ }
+ }
+ rieman2.Update();
+ fitterTC.Eval();
+ 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*(xcl[iLayer] - xref2);
+ if (TMath::Abs(cseed[iLayer].fZProb-zT2)>padlength[iLayer]*0.5+1)
+ acceptablez = kFALSE;
+ }
+ }
+ if (!acceptablez){
+ 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 chi2TC = fitterTC.GetChisquare()/Float_t(npointsT);
+ //
+ Double_t polz1c = fitterTC.GetParameter(2);
+ Double_t polz0c = polz1c*xref2;
+ //
+ Double_t aC = fitterTC.GetParameter(0);
+ Double_t bC = fitterTC.GetParameter(1);
+ Double_t cC = aC/TMath::Sqrt(bC*bC+1.); // curvature
+ //
+ Double_t aR = fitterT2.GetParameter(0);
+ Double_t bR = fitterT2.GetParameter(1);
+ Double_t dR = fitterT2.GetParameter(2);
+ Double_t cR = 1+bR*bR-dR*aR;
+ Double_t dca = 0.;
+ if (cR>0){
+ dca = -dR/(TMath::Sqrt(1+bR*bR-dR*aR)+TMath::Sqrt(1+bR*bR));
+ cR = aR/TMath::Sqrt(cR);
+ }
+ //
+ Double_t chi2ZT2=0, chi2ZTC=0;
+ for (Int_t iLayer=0; iLayer<6;iLayer++){
+ if (cseed[iLayer].IsOK()){
+ Double_t zT2 = rpolz0+rpolz1*(xcl[iLayer] - xref2);
+ Double_t zTC = polz0c+polz1c*(xcl[iLayer] - xref2);
+ chi2ZT2 += TMath::Abs(cseed[iLayer].fMeanz-zT2);
+ chi2ZTC += TMath::Abs(cseed[iLayer].fMeanz-zTC);
+ }
+ }
+ chi2ZT2/=TMath::Max((nlayers-3.),1.);
+ chi2ZTC/=TMath::Max((nlayers-3.),1.);
+ //
+ //
+ //
+ AliTRDseed::FitRiemanTilt(cseed, kTRUE);
+ Float_t sumdaf = 0;
+ for (Int_t iLayer=0;iLayer<6;iLayer++){
+ if (cseed[iLayer].IsOK())
+ sumdaf += TMath::Abs((cseed[iLayer].fYfit[1]-cseed[iLayer].fYref[1])/cseed[iLayer].fSigmaY2);
+ }
+ sumdaf /= Float_t (nlayers-2.);
+ //
+ // likelihoods for full track
+ //
+ Double_t likezf = TMath::Exp(-chi2ZF*0.14);
+ Double_t likechi2C = TMath::Exp(-chi2TC*0.677);
+ Double_t likechi2TR = TMath::Exp(-chi2TR*0.78);
+ Double_t likeaf = TMath::Exp(-sumdaf*3.23);
+ seedquality2[registered] = likezf*likechi2TR*likeaf;
+// Bool_t isGold = kFALSE;
+//
+// if (nlayers == 6 && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE; // gold
+// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
+// if (isGold &&nusedf<10){
+// for (Int_t jLayer=0;jLayer<6;jLayer++){
+// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
+// seed[index][jLayer].UseClusters(); //sign gold
+// }
+// }
+ //
+ //
+ //
+ Int_t index0=0;
+ if (!cseed[0].IsOK()){
+ index0 = 1;
+ if (!cseed[1].IsOK()) index0 = 2;
+ }
+ seedparams[registered][0] = cseed[index0].fX0;
+ seedparams[registered][1] = cseed[index0].fYref[0];
+ seedparams[registered][2] = cseed[index0].fZref[0];
+ seedparams[registered][5] = cR;
+ seedparams[registered][3] = cseed[index0].fX0*cR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
+ seedparams[registered][4] = cseed[index0].fZref[1]/
+ TMath::Sqrt(1+cseed[index0].fYref[1]*cseed[index0].fYref[1]);
+ seedparams[registered][6] = ns;
+ //
+ //
+ Int_t labels[12], outlab[24];
+ Int_t nlab=0;
+ for (Int_t iLayer=0;iLayer<6;iLayer++){
+ if (!cseed[iLayer].IsOK()) continue;
+ if (cseed[iLayer].fLabels[0]>=0) {
+ labels[nlab] = cseed[iLayer].fLabels[0];
+ nlab++;
+ }
+ if (cseed[iLayer].fLabels[1]>=0) {
+ labels[nlab] = cseed[iLayer].fLabels[1];
+ nlab++;
+ }
+ }
+ Freq(nlab,labels,outlab,kFALSE);
+ Int_t label = outlab[0];
+ Int_t frequency = outlab[1];
+ for (Int_t iLayer=0;iLayer<6;iLayer++){
+ cseed[iLayer].fFreq = frequency;
+ cseed[iLayer].fC = cR;
+ cseed[iLayer].fCC = cC;
+ cseed[iLayer].fChi2 = chi2TR;
+ cseed[iLayer].fChi2Z = chi2ZF;
+ }
+ //
+ if (1||(!isFake)){ //debugging print
+ Float_t zvertex = GetZ();
+ TTreeSRedirector& cstream = *fDebugStreamer;
+ if (AliTRDReconstructor::StreamLevel()>0)
+ cstream<<"Seeds1"<<
+ "isFake="<<isFake<<
+ "Vertex="<<zvertex<<
+ "Rieman2.="<<&rieman2<<
+ "Rieman.="<<&rieman<<
+ "Xref="<<xref<<
+ "X0="<<xcl[0]<<
+ "X1="<<xcl[1]<<
+ "X2="<<xcl[2]<<
+ "X3="<<xcl[3]<<
+ "X4="<<xcl[4]<<
+ "X5="<<xcl[5]<<
+ "Chi2R="<<chi2R<<
+ "Chi2Z="<<chi2Z<<
+ "Chi2RF="<<chi2RF<< //chi2 of trackletes on full track
+ "Chi2ZF="<<chi2ZF<< //chi2 z on tracklets on full track
+ "Chi2ZT2="<<chi2ZT2<< //chi2 z on tracklets on full track - rieman tilt
+ "Chi2ZTC="<<chi2ZTC<< //chi2 z on tracklets on full track - rieman tilt const
+ //
+ "Chi2TR="<<chi2TR<< //chi2 without vertex constrain
+ "Chi2TC="<<chi2TC<< //chi2 with vertex constrain
+ "C="<<curv<< // non constrained - no tilt correction
+ "DR="<<dR<< // DR parameter - tilt correction
+ "DCA="<<dca<< // DCA - tilt correction
+ "CR="<<cR<< // non constrained curvature - tilt correction
+ "CC="<<cC<< // constrained curvature
+ "Polz0="<<polz0c<<
+ "Polz1="<<polz1c<<
+ "RPolz0="<<rpolz0<<
+ "RPolz1="<<rpolz1<<
+ "Ncl="<<nclusters<<
+ "Nlayers="<<nlayers<<
+ "NUsedS="<<nusedCl<<
+ "NUsed="<<nusedf<<
+ "Findable="<<findable<<
+ "Like="<<like<<
+ "LikePrim="<<likePrim<<
+ "Likechi2C="<<likechi2C<<
+ "Likechi2TR="<<likechi2TR<<
+ "Likezf="<<likezf<<
+ "LikeF="<<seedquality2[registered]<<
+ "S0.="<<&cseed[0]<<
+ "S1.="<<&cseed[1]<<
+ "S2.="<<&cseed[2]<<
+ "S3.="<<&cseed[3]<<
+ "S4.="<<&cseed[4]<<
+ "S5.="<<&cseed[5]<<
+ "SB0.="<<&seedb[0]<<
+ "SB1.="<<&seedb[1]<<
+ "SB2.="<<&seedb[2]<<
+ "SB3.="<<&seedb[3]<<
+ "SB4.="<<&seedb[4]<<
+ "SB5.="<<&seedb[5]<<
+ "Label="<<label<<
+ "Freq="<<frequency<<
+ "sLayer="<<sLayer<<
+ "\n";
+ }
+ if (registered<kMaxSeed-1) {
+ registered++;
+ cseed = seed[registered];
+ }
+ }// end of loop over layer 1
+ } // end of loop over layer 0
+ } // end of loop over layer 3
+ } // end of loop over seeding time bins
+ //
+ // choos best
+ //
+ TMath::Sort(registered,seedquality2,sort,kTRUE);
+ Bool_t signedseed[kMaxSeed];
+ for (Int_t i=0;i<registered;i++){
+ signedseed[i]= kFALSE;
+ }
+ for (Int_t iter=0; iter<5; iter++){
+ for (Int_t iseed=0;iseed<registered;iseed++){
+ Int_t index = sort[iseed];
+ if (signedseed[index]) continue;
+ Int_t labelsall[1000];
+ Int_t nlabelsall=0;
+ Int_t naccepted=0;;
+ Int_t sLayer = seedlayer[index];
+ Int_t ncl = 0;
+ Int_t nused = 0;
+ Int_t nlayers =0;
+ Int_t findable = 0;
+ for (Int_t jLayer=0;jLayer<6;jLayer++){
+ if (TMath::Abs(seed[index][jLayer].fYref[0]/xcl[jLayer])<0.15)
+ findable++;
+ if (seed[index][jLayer].IsOK()){
+ seed[index][jLayer].UpdateUsed();
+ ncl +=seed[index][jLayer].fN2;
+ nused +=seed[index][jLayer].fNUsed;
+ nlayers++;
+ //cooking label
+ for (Int_t itime=0;itime<25;itime++){
+ if (seed[index][jLayer].fUsable[itime]){
+ naccepted++;
+ for (Int_t ilab=0;ilab<3;ilab++){
+ Int_t tindex = seed[index][jLayer].fClusters[itime]->GetLabel(ilab);
+ if (tindex>=0){
+ labelsall[nlabelsall] = tindex;
+ nlabelsall++;
+ }
+ }
+ }
+ }
+ }