+ //
+ TStopwatch timer;
+ timer.Start();
+ Int_t ncandidates =0;
+ Int_t nall =0;
+ Int_t ntracks=0;
+ Double_t phase[2][2],radius[2];
+
+ //
+ // Find circling track
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ //
+ for (Int_t i0=0;i0<nentries;i0++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
+ if (!track0) continue;
+ if (track0->fN<40) continue;
+ if (TMath::Abs(1./track0->fP4)>200) continue;
+ for (Int_t i1=i0+1;i1<nentries;i1++){
+ AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
+ if (!track1) continue;
+ if (track1->fN<40) continue;
+ if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue;
+ if (track0->fBConstrain&&track1->fBConstrain) continue;
+ if (TMath::Abs(1./track1->fP4)>200) continue;
+ if (track1->fP4*track0->fP4>0) continue;
+ if (track1->fP3*track0->fP3>0) continue;
+ if (max(TMath::Abs(1./track0->fP4),TMath::Abs(1./track1->fP4))>190) continue;
+ if (track0->fBConstrain&&TMath::Abs(track1->fP4)<TMath::Abs(track0->fP4)) continue; //returning - lower momenta
+ if (track1->fBConstrain&&TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)) continue; //returning - lower momenta
+ //
+ Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
+ if (mindcar<5) continue;
+ Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
+ if (mindcaz<5) continue;
+ if (mindcar+mindcaz<20) continue;
+ //
+ //
+ Float_t xc0 = helixes[i0].GetHelix(6);
+ Float_t yc0 = helixes[i0].GetHelix(7);
+ Float_t r0 = helixes[i0].GetHelix(8);
+ Float_t xc1 = helixes[i1].GetHelix(6);
+ Float_t yc1 = helixes[i1].GetHelix(7);
+ Float_t r1 = helixes[i1].GetHelix(8);
+
+ Float_t rmean = (r0+r1)*0.5;
+ Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
+ //if (delta>30) continue;
+ if (delta>rmean*0.25) continue;
+ if (TMath::Abs(r0-r1)/rmean>0.3) continue;
+ //
+ Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
+ if (npoints==0) continue;
+ helixes[i0].GetClosestPhases(helixes[i1], phase);
+ //
+ Double_t xyz0[3];
+ Double_t xyz1[3];
+ Double_t hangles[3];
+ helixes[i0].Evaluate(phase[0][0],xyz0);
+ helixes[i1].Evaluate(phase[0][1],xyz1);
+
+ helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
+ Double_t deltah[2],deltabest;
+ if (hangles[2]<2.8) continue;
+ /*
+ cstream<<"C"<<track0->fLab<<track1->fLab<<
+ track0->fP3<<track1->fP3<<
+ track0->fP4<<track1->fP4<<
+ delta<<rmean<<npoints<<
+ hangles[0]<<hangles[2]<<
+ xyz0[2]<<xyz1[2]<<radius[0]<<"\n";
+ */
+ if (npoints>0){
+ Int_t ibest=0;
+ helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
+ if (npoints==2){
+ helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
+ if (deltah[1]<deltah[0]) ibest=1;
+ }
+ deltabest = TMath::Sqrt(deltah[ibest]);
+ helixes[i0].Evaluate(phase[ibest][0],xyz0);
+ helixes[i1].Evaluate(phase[ibest][1],xyz1);
+ helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
+ Double_t radiusbest = TMath::Sqrt(radius[ibest]);
+ //
+ if (deltabest>6) continue;
+ if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
+ Bool_t sign =kFALSE;
+ if (hangles[2]>3.06) sign =kTRUE;
+ //
+ if (sign){
+ circular[i0] = kTRUE;
+ circular[i1] = kTRUE;
+ if (TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)){
+ track0->fCircular += 1;
+ track1->fCircular += 2;
+ }
+ else{
+ track1->fCircular += 1;
+ track0->fCircular += 2;
+ }
+ }
+ if (sign&&AliTPCReconstructor::StreamLevel()>1){
+ //debug stream
+ cstream<<"Curling"<<
+ "lab0="<<track0->fLab<<
+ "lab1="<<track1->fLab<<
+ "Tr0.="<<track0<<
+ "Tr1.="<<track1<<
+ "dca0="<<dca[i0]<<
+ "dca1="<<dca[i1]<<
+ "mindcar="<<mindcar<<
+ "mindcaz="<<mindcaz<<
+ "delta="<<delta<<
+ "rmean="<<rmean<<
+ "npoints="<<npoints<<
+ "hangles0="<<hangles[0]<<
+ "hangles2="<<hangles[2]<<
+ "xyz0="<<xyz0[2]<<
+ "xyzz1="<<xyz1[2]<<
+ "z0="<<z0[i0]<<
+ "z1="<<z0[i1]<<
+ "radius="<<radiusbest<<
+ "deltabest="<<deltabest<<
+ "phase0="<<phase[ibest][0]<<
+ "phase1="<<phase[ibest][1]<<
+ "\n";
+ }
+ }
+ }
+ }
+ //
+ // Finf kinks loop
+ //
+ //
+ for (Int_t i =0;i<nentries;i++){
+ if (sign[i]==0) continue;
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+ ntracks++;
+ //
+ Double_t cradius0 = 40*40;
+ Double_t cradius1 = 270*270;
+ Double_t cdist1=8.;
+ Double_t cdist2=8.;
+ Double_t cdist3=0.55;
+ for (Int_t j =i+1;j<nentries;j++){
+ nall++;
+ if (sign[j]*sign[i]<1) continue;
+ if ( (nclusters[i]+nclusters[j])>200) continue;
+ if ( (nclusters[i]+nclusters[j])<80) continue;
+ if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
+ if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
+ //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
+ Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
+ if (npoints<1) continue;
+ // cuts on radius
+ if (npoints==1){
+ if (radius[0]<cradius0||radius[0]>cradius1) continue;
+ }
+ else{
+ if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
+ }
+ //
+ Double_t delta1=10000,delta2=10000;
+ // cuts on the intersection radius
+ helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
+ if (radius[0]<20&&delta1<1) continue; //intersection at vertex
+ if (radius[0]<10&&delta1<3) continue; //intersection at vertex
+ if (npoints==2){
+ helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
+ if (radius[1]<20&&delta2<1) continue; //intersection at vertex
+ if (radius[1]<10&&delta2<3) continue; //intersection at vertex
+ }
+ //
+ Double_t distance1 = TMath::Min(delta1,delta2);
+ if (distance1>cdist1) continue; // cut on DCA linear approximation
+ //
+ npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
+ helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
+ if (radius[0]<20&&delta1<1) continue; //intersection at vertex
+ if (radius[0]<10&&delta1<3) continue; //intersection at vertex
+ //
+ if (npoints==2){
+ helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
+ if (radius[1]<20&&delta2<1) continue; //intersection at vertex
+ if (radius[1]<10&&delta2<3) continue; //intersection at vertex
+ }
+ distance1 = TMath::Min(delta1,delta2);
+ Float_t rkink =0;
+ if (delta1<delta2){
+ rkink = TMath::Sqrt(radius[0]);
+ }
+ else{
+ rkink = TMath::Sqrt(radius[1]);
+ }
+ if (distance1>cdist2) continue;
+ //
+ //
+ AliTPCseed * track1 = (AliTPCseed*)array->At(j);
+ //
+ //
+ Int_t row0 = GetRowNumber(rkink);
+ if (row0<10) continue;
+ if (row0>150) continue;
+ //
+ //
+ Float_t dens00=-1,dens01=-1;
+ Float_t dens10=-1,dens11=-1;
+ //
+ Int_t found,foundable,shared;
+ track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
+ if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
+ track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
+ if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
+ //
+ track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
+ if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
+ track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
+ if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
+ //
+ if (dens00<dens10 && dens01<dens11) continue;
+ if (dens00>dens10 && dens01>dens11) continue;
+ if (TMath::Max(dens00,dens10)<0.1) continue;
+ if (TMath::Max(dens01,dens11)<0.3) continue;
+ //
+ if (TMath::Min(dens00,dens10)>0.6) continue;
+ if (TMath::Min(dens01,dens11)>0.6) continue;
+
+ //
+ AliTPCseed * ktrack0, *ktrack1;
+ if (dens00>dens10){
+ ktrack0 = track0;
+ ktrack1 = track1;
+ }
+ else{
+ ktrack0 = track1;
+ ktrack1 = track0;
+ }
+ if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
+ AliExternalTrackParam paramm(*ktrack0);
+ AliExternalTrackParam paramd(*ktrack1);
+ if (row0>60&&ktrack1->GetReference().GetX()>90.) new (¶md) AliExternalTrackParam(ktrack1->GetReference());
+ //
+ //
+ kink->SetMother(paramm);
+ kink->SetDaughter(paramd);
+ kink->Update();
+
+ Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
+ Int_t index[4];
+ fParam->Transform0to1(x,index);
+ fParam->Transform1to2(x,index);
+ row0 = GetRowNumber(x[0]);
+
+ if (kink->GetR()<100) continue;
+ if (kink->GetR()>240) continue;
+ if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
+ if (kink->GetDistance()>cdist3) continue;
+ Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
+ if (dird<0) continue;
+
+ Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
+ if (dirm<0) continue;
+ Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
+ if (mpt<0.2) continue;
+
+ if (mpt<1){
+ //for high momenta momentum not defined well in first iteration
+ Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->P();
+ if (qt>0.35) continue;
+ }
+
+ kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
+ kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
+ if (dens00>dens10){
+ kink->SetTPCDensity(dens00,0,0);
+ kink->SetTPCDensity(dens01,0,1);
+ kink->SetTPCDensity(dens10,1,0);
+ kink->SetTPCDensity(dens11,1,1);
+ kink->SetIndex(i,0);
+ kink->SetIndex(j,1);
+ }
+ else{
+ kink->SetTPCDensity(dens10,0,0);
+ kink->SetTPCDensity(dens11,0,1);
+ kink->SetTPCDensity(dens00,1,0);
+ kink->SetTPCDensity(dens01,1,1);
+ kink->SetIndex(j,0);
+ kink->SetIndex(i,1);
+ }
+
+ if (mpt<1||kink->GetAngle(2)>0.1){
+ // angle and densities not defined yet
+ if (kink->GetTPCDensityFactor()<0.8) continue;
+ if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
+ if (kink->GetAngle(2)*ktrack0->P()<0.003) continue; //too small angle
+ if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
+ if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
+
+ Float_t criticalangle = track0->fC22+track0->fC33;
+ criticalangle+= track1->fC22+track1->fC33;
+ criticalangle= 3*TMath::Sqrt(criticalangle);
+ if (criticalangle>0.02) criticalangle=0.02;
+ if (kink->GetAngle(2)<criticalangle) continue;
+ }
+ //
+ Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
+ Float_t shapesum =0;
+ Float_t sum = 0;
+ for ( Int_t row = row0-drow; row<row0+drow;row++){
+ if (row<0) continue;
+ if (row>155) continue;
+ if (ktrack0->fClusterPointer[row]){
+ AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
+ shapesum+=point->GetSigmaY()+point->GetSigmaZ();
+ sum++;
+ }
+ if (ktrack1->fClusterPointer[row]){
+ AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
+ shapesum+=point->GetSigmaY()+point->GetSigmaZ();
+ sum++;
+ }
+ }
+ if (sum<4){
+ kink->SetShapeFactor(-1.);
+ }
+ else{
+ kink->SetShapeFactor(shapesum/sum);
+ }
+ // esd->AddKink(kink);
+ kinks->AddLast(kink);
+ kink = new AliESDkink;
+ ncandidates++;
+ }
+ }
+ //
+ // sort the kinks according quality - and refit them towards vertex
+ //
+ Int_t nkinks = kinks->GetEntriesFast();
+ Float_t *quality = new Float_t[nkinks];
+ Int_t *indexes = new Int_t[nkinks];
+ AliTPCseed *mothers = new AliTPCseed[nkinks];
+ AliTPCseed *daughters = new AliTPCseed[nkinks];
+ //
+ //
+ for (Int_t i=0;i<nkinks;i++){
+ quality[i] =100000;
+ AliESDkink *kink = (AliESDkink*)kinks->At(i);
+ //
+ // refit kinks towards vertex
+ //
+ Int_t index0 = kink->GetIndex(0);
+ Int_t index1 = kink->GetIndex(1);
+ AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
+ AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
+ //
+ Int_t sumn=ktrack0->fN+ktrack1->fN;
+ //
+ // Refit Kink under if too small angle
+ //
+ if (kink->GetAngle(2)<0.05){
+ kink->SetTPCRow0(GetRowNumber(kink->GetR()));
+ Int_t row0 = kink->GetTPCRow0();
+ Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
+ //
+ //
+ Int_t last = row0-drow;
+ if (last<40) last=40;
+ if (last<ktrack0->fFirstPoint+25) last = ktrack0->fFirstPoint+25;
+ AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
+ //
+ //
+ Int_t first = row0+drow;
+ if (first>130) first=130;
+ if (first>ktrack1->fLastPoint-25) first = TMath::Max(ktrack1->fLastPoint-25,30);
+ AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
+ //
+ if (seed0 && seed1){
+ kink->SetStatus(1,8);
+ if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
+ row0 = GetRowNumber(kink->GetR());
+ sumn = seed0->fN+seed1->fN;
+ new (&mothers[i]) AliTPCseed(*seed0);
+ new (&daughters[i]) AliTPCseed(*seed1);
+ }
+ else{
+ delete kinks->RemoveAt(i);
+ if (seed0) delete seed0;
+ if (seed1) delete seed1;
+ continue;
+ }
+ if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
+ delete kinks->RemoveAt(i);
+ if (seed0) delete seed0;
+ if (seed1) delete seed1;
+ continue;
+ }
+ //
+ delete seed0;
+ delete seed1;
+ }
+ //
+ if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
+ }
+ TMath::Sort(nkinks,quality,indexes,kFALSE);
+ //
+ //remove double find kinks
+ //
+ for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
+ AliESDkink * kink0 = (AliESDkink*) kinks->At(indexes[ikink0]);
+ if (!kink0) continue;
+ //
+ for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
+ if (!kink0) continue;
+ AliESDkink * kink1 = (AliESDkink*) kinks->At(indexes[ikink1]);
+ if (!kink1) continue;
+ // if not close kink continue
+ if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
+ if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
+ if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
+ //
+ AliTPCseed &mother0 = mothers[indexes[ikink0]];
+ AliTPCseed &daughter0 = daughters[indexes[ikink0]];
+ AliTPCseed &mother1 = mothers[indexes[ikink1]];
+ AliTPCseed &daughter1 = daughters[indexes[ikink1]];
+ Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
+ //
+ Int_t same = 0;
+ Int_t both = 0;
+ Int_t samem = 0;
+ Int_t bothm = 0;
+ Int_t samed = 0;
+ Int_t bothd = 0;
+ //
+ for (Int_t i=0;i<row0;i++){
+ if (mother0.fIndex[i]>0 && mother1.fIndex[i]>0){
+ both++;
+ bothm++;
+ if (mother0.fIndex[i]==mother1.fIndex[i]){
+ same++;
+ samem++;
+ }
+ }
+ }
+
+ for (Int_t i=row0;i<158;i++){
+ if (daughter0.fIndex[i]>0 && daughter0.fIndex[i]>0){
+ both++;
+ bothd++;
+ if (mother0.fIndex[i]==mother1.fIndex[i]){
+ same++;
+ samed++;
+ }
+ }
+ }
+ Float_t ratio = Float_t(same+1)/Float_t(both+1);
+ Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
+ Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
+ if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
+ Int_t sum0 = mother0.fN+daughter0.fN;
+ Int_t sum1 = mother1.fN+daughter1.fN;
+ if (sum1>sum0){
+ shared[kink0->GetIndex(0)]= kTRUE;
+ shared[kink0->GetIndex(1)]= kTRUE;
+ delete kinks->RemoveAt(indexes[ikink0]);
+ }
+ else{
+ shared[kink1->GetIndex(0)]= kTRUE;
+ shared[kink1->GetIndex(1)]= kTRUE;
+ delete kinks->RemoveAt(indexes[ikink1]);
+ }
+ }
+ }
+ }
+
+
+ for (Int_t i=0;i<nkinks;i++){
+ AliESDkink * kink = (AliESDkink*) kinks->At(indexes[i]);
+ if (!kink) continue;
+ kink->SetTPCRow0(GetRowNumber(kink->GetR()));
+ Int_t index0 = kink->GetIndex(0);
+ Int_t index1 = kink->GetIndex(1);
+ if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
+ kink->SetMultiple(usage[index0],0);
+ kink->SetMultiple(usage[index1],1);
+ if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
+ if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
+ if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
+ if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
+
+ AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
+ AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
+ if (!ktrack0 || !ktrack1) continue;
+ Int_t index = esd->AddKink(kink);
+ //
+ //
+ if ( ktrack0->fKinkIndexes[0]==0 && ktrack1->fKinkIndexes[0]==0) { //best kink
+ if (mothers[indexes[i]].fN>20 && daughters[indexes[i]].fN>20 && (mothers[indexes[i]].fN+daughters[indexes[i]].fN)>100){
+ new (ktrack0) AliTPCseed(mothers[indexes[i]]);
+ new (ktrack1) AliTPCseed(daughters[indexes[i]]);
+ }
+ }
+ //
+ ktrack0->fKinkIndexes[usage[index0]] = -(index+1);
+ ktrack1->fKinkIndexes[usage[index1]] = (index+1);
+ usage[index0]++;
+ usage[index1]++;
+ }
+ //
+ // Remove tracks corresponding to shared kink's
+ //
+ for (Int_t i=0;i<nentries;i++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+ if (!track0) continue;
+ if (track0->fKinkIndexes[0]!=0) continue;
+ if (shared[i]) delete array->RemoveAt(i);
+ }
+
+ //
+ //
+ RemoveUsed2(array,0.5,0.4,30);
+ UnsignClusters();
+ for (Int_t i=0;i<nentries;i++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+ if (!track0) continue;
+ track0->CookdEdx(0.02,0.6);
+ track0->CookPID();
+ }
+ //
+ for (Int_t i=0;i<nentries;i++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+ if (!track0) continue;
+ if (track0->Pt()<1.4) continue;
+ //remove double high momenta tracks - overlapped with kink candidates
+ Int_t shared=0;
+ Int_t all =0;
+ for (Int_t icl=track0->fFirstPoint;icl<track0->fLastPoint; icl++){
+ if (track0->fClusterPointer[icl]!=0){
+ all++;
+ if (track0->fClusterPointer[icl]->IsUsed(10)) shared++;
+ }
+ }
+ if (Float_t(shared+1)/Float_t(nall+1)>0.5) {
+ delete array->RemoveAt(i);
+ continue;
+ }
+ //
+ if (track0->fKinkIndexes[0]!=0) continue;
+ if (track0->GetNumberOfClusters()<80) continue;
+
+ AliTPCseed *pmother = new AliTPCseed();
+ AliTPCseed *pdaughter = new AliTPCseed();
+ AliESDkink *pkink = new AliESDkink;
+
+ AliTPCseed & mother = *pmother;
+ AliTPCseed & daughter = *pdaughter;
+ AliESDkink & kink = *pkink;
+ if (CheckKinkPoint(track0,mother,daughter, kink)){
+ if (mother.fN<30||daughter.fN<20) {
+ delete pmother;
+ delete pdaughter;
+ delete pkink;
+ continue; //too short tracks
+ }
+ if (mother.Pt()<1.4) {
+ delete pmother;
+ delete pdaughter;
+ delete pkink;
+ continue;
+ }
+ Int_t row0= kink.GetTPCRow0();
+ if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
+ delete pmother;
+ delete pdaughter;
+ delete pkink;
+ continue;
+ }
+ //
+ Int_t index = esd->AddKink(&kink);
+ mother.fKinkIndexes[0] = -(index+1);
+ daughter.fKinkIndexes[0] = index+1;
+ if (mother.fN>50) {
+ delete array->RemoveAt(i);
+ array->AddAt(new AliTPCseed(mother),i);
+ }
+ else{
+ array->AddLast(new AliTPCseed(mother));
+ }
+ array->AddLast(new AliTPCseed(daughter));
+ for (Int_t icl=0;icl<row0;icl++) {
+ if (mother.fClusterPointer[icl]) mother.fClusterPointer[icl]->Use(20);
+ }
+ //
+ for (Int_t icl=row0;icl<158;icl++) {
+ if (daughter.fClusterPointer[icl]) daughter.fClusterPointer[icl]->Use(20);
+ }
+ //
+ }
+ delete pmother;
+ delete pdaughter;
+ delete pkink;
+ }
+
+ delete [] daughters;
+ delete [] mothers;
+ //
+ //
+ delete [] dca;
+ delete []circular;
+ delete []shared;
+ delete []quality;
+ delete []indexes;
+ //
+ delete kink;
+ delete[]fim;
+ delete[] zm;
+ delete[] z0;
+ delete [] usage;
+ delete[] alpha;
+ delete[] nclusters;
+ delete[] sign;
+ delete[] helixes;
+ kinks->Delete();
+ delete kinks;
+
+ printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
+ timer.Print();
+}
+
+void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESD *esd)
+{
+ //
+ // find V0s
+ //
+ //
+ TObjArray *tpcv0s = new TObjArray(100000);
+ Int_t nentries = array->GetEntriesFast();
+ AliHelix *helixes = new AliHelix[nentries];
+ Int_t *sign = new Int_t[nentries];
+ Float_t *alpha = new Float_t[nentries];
+ Float_t *z0 = new Float_t[nentries];
+ Float_t *dca = new Float_t[nentries];
+ Float_t *sdcar = new Float_t[nentries];
+ Float_t *cdcar = new Float_t[nentries];
+ Float_t *pulldcar = new Float_t[nentries];
+ Float_t *pulldcaz = new Float_t[nentries];
+ Float_t *pulldca = new Float_t[nentries];
+ Bool_t *isPrim = new Bool_t[nentries];
+ const AliESDVertex * primvertex = esd->GetVertex();
+ Double_t zvertex = primvertex->GetZv();
+ //
+ // nentries = array->GetEntriesFast();
+ //
+ for (Int_t i=0;i<nentries;i++){
+ sign[i]=0;
+ isPrim[i]=0;
+ AliTPCseed* track = (AliTPCseed*)array->At(i);
+ if (!track) continue;
+ track->GetV0Indexes()[0] = 0; //rest v0 indexes
+ track->GetV0Indexes()[1] = 0; //rest v0 indexes
+ track->GetV0Indexes()[2] = 0; //rest v0 indexes
+ //
+ alpha[i] = track->GetAlpha();
+ new (&helixes[i]) AliHelix(*track);
+ Double_t xyz[3];
+ helixes[i].Evaluate(0,xyz);
+ sign[i] = (track->GetC()>0) ? -1:1;
+ Double_t x,y,z;
+ x=160;
+ z0[i]=1000;
+ if (track->GetProlongation(0,y,z)) z0[i] = z;
+ dca[i] = track->GetD(0,0);
+ //
+ // dca error parrameterezation + pulls
+ //
+ sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->fP4)*(100*track->fP4));
+ if (TMath::Abs(track->fP3)>1) sdcar[i]*=2.5;
+ cdcar[i] = TMath::Exp((TMath::Abs(track->fP4)-0.0106)*525.3);
+ pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
+ pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
+ pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
+ if (track->fTPCr[1]+track->fTPCr[2]+track->fTPCr[3]>0.5) {
+ if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
+ }
+ if (track->fTPCr[4]>0.5) {
+ if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
+ }
+ if (track->fTPCr[0]>0.4) {
+ isPrim[i]=kFALSE; //electron no sigma cut
+ }
+ }
+ //
+ //
+ TStopwatch timer;
+ timer.Start();
+ Int_t ncandidates =0;
+ Int_t nall =0;
+ Int_t ntracks=0;
+ Double_t phase[2][2],radius[2];
+ //
+ // Finf V0s loop
+ //
+ //
+ // //
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
+ AliESDV0MI vertex;
+ Double_t cradius0 = 10*10;
+ Double_t cradius1 = 200*200;
+ Double_t cdist1=3.;
+ Double_t cdist2=4.;
+ Double_t cpointAngle = 0.95;
+ //
+ Double_t delta[2]={10000,10000};
+ for (Int_t i =0;i<nentries;i++){
+ if (sign[i]==0) continue;
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+ if (!track0) continue;
+ if (AliTPCReconstructor::StreamLevel()>0){
+ cstream<<"Tracks"<<
+ "Tr0.="<<track0<<
+ "dca="<<dca[i]<<
+ "z0="<<z0[i]<<
+ "zvertex="<<zvertex<<
+ "sdcar0="<<sdcar[i]<<
+ "cdcar0="<<cdcar[i]<<
+ "pulldcar0="<<pulldcar[i]<<
+ "pulldcaz0="<<pulldcaz[i]<<
+ "pulldca0="<<pulldca[i]<<
+ "isPrim="<<isPrim[i]<<
+ "\n";
+ }
+ //
+ if (track0->fP4<0) continue;
+ if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
+ //
+ if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
+ ntracks++;
+ // debug output
+
+
+ for (Int_t j =0;j<nentries;j++){
+ AliTPCseed * track1 = (AliTPCseed*)array->At(j);
+ if (!track1) continue;
+ if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
+ if (sign[j]*sign[i]>0) continue;
+ if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
+ if (track0->fCircular+track1->fCircular>1) continue; //circling -returning track
+ nall++;
+ //
+ // DCA to prim vertex cut
+ //
+ //
+ delta[0]=10000;
+ delta[1]=10000;
+
+ Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
+ if (npoints<1) continue;
+ Int_t iclosest=0;
+ // cuts on radius
+ if (npoints==1){
+ if (radius[0]<cradius0||radius[0]>cradius1) continue;
+ helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
+ if (delta[0]>cdist1) continue;
+ }
+ else{
+ if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
+ helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
+ helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
+ if (delta[1]<delta[0]) iclosest=1;
+ if (delta[iclosest]>cdist1) continue;
+ }
+ helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
+ if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
+ //
+ Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
+ if (pointAngle<cpointAngle) continue;
+ //
+ Bool_t isGamma = kFALSE;
+ vertex.SetP(*track0); //track0 - plus
+ vertex.SetM(*track1); //track1 - minus
+ vertex.Update(fprimvertex);
+ if (track0->fTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
+ Double_t pointAngle2 = vertex.GetPointAngle();
+ //continue;
+ if (vertex.GetPointAngle()<cpointAngle && (!isGamma)) continue; // point angle cut
+ if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
+ Float_t sigmae = 0.15*0.15;
+ if (vertex.GetRr()<80)
+ sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
+ sigmae+= TMath::Sqrt(sigmae);
+ if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
+ Float_t densb0=0,densb1=0,densa0=0,densa1=0;
+ Int_t row0 = GetRowNumber(vertex.GetRr());
+ if (row0>15){
+ if (vertex.GetDist2()>0.2) continue;
+ densb0 = track0->Density2(0,row0-5);
+ densb1 = track1->Density2(0,row0-5);
+ if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
+ densa0 = track0->Density2(row0+5,row0+40);
+ densa1 = track1->Density2(row0+5,row0+40);
+ if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
+ }
+ else{
+ densa0 = track0->Density2(0,40); //cluster density
+ densa1 = track1->Density2(0,40); //cluster density
+ if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
+ }
+ vertex.SetLab(0,track0->GetLabel());
+ vertex.SetLab(1,track1->GetLabel());
+ vertex.SetChi2After((densa0+densa1)*0.5);
+ vertex.SetChi2Before((densb0+densb1)*0.5);
+ vertex.SetIndex(0,i);
+ vertex.SetIndex(1,j);
+ vertex.SetStatus(1); // TPC v0 candidate
+ vertex.SetRp(track0->fTPCr);
+ vertex.SetRm(track1->fTPCr);
+ tpcv0s->AddLast(new AliESDV0MI(vertex));
+ ncandidates++;
+ {
+ Int_t eventNr = esd->GetEventNumber();
+ Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
+ Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
+ if (AliTPCReconstructor::StreamLevel()>0)
+ cstream<<"V0"<<
+ "Event="<<eventNr<<
+ "vertex.="<<&vertex<<
+ "Tr0.="<<track0<<
+ "lab0="<<track0->fLab<<
+ "Helix0.="<<&helixes[i]<<
+ "Tr1.="<<track1<<
+ "lab1="<<track1->fLab<<
+ "Helix1.="<<&helixes[j]<<
+ "pointAngle="<<pointAngle<<
+ "pointAngle2="<<pointAngle2<<
+ "dca0="<<dca[i]<<
+ "dca1="<<dca[j]<<
+ "z0="<<z0[i]<<
+ "z1="<<z0[j]<<
+ "zvertex="<<zvertex<<
+ "circular0="<<track0->fCircular<<
+ "circular1="<<track1->fCircular<<
+ "npoints="<<npoints<<
+ "radius0="<<radius[0]<<
+ "delta0="<<delta[0]<<
+ "radius1="<<radius[1]<<
+ "delta1="<<delta[1]<<
+ "radiusm="<<radiusm<<
+ "deltam="<<deltam<<
+ "sdcar0="<<sdcar[i]<<
+ "sdcar1="<<sdcar[j]<<
+ "cdcar0="<<cdcar[i]<<
+ "cdcar1="<<cdcar[j]<<
+ "pulldcar0="<<pulldcar[i]<<
+ "pulldcar1="<<pulldcar[j]<<
+ "pulldcaz0="<<pulldcaz[i]<<
+ "pulldcaz1="<<pulldcaz[j]<<
+ "pulldca0="<<pulldca[i]<<
+ "pulldca1="<<pulldca[j]<<
+ "densb0="<<densb0<<
+ "densb1="<<densb1<<
+ "densa0="<<densa0<<
+ "densa1="<<densa1<<
+ "sigmae="<<sigmae<<
+ "\n";
+ }
+ }
+ }
+ Float_t *quality = new Float_t[ncandidates];
+ Int_t *indexes = new Int_t[ncandidates];
+ Int_t naccepted =0;
+ for (Int_t i=0;i<ncandidates;i++){
+ quality[i] = 0;
+ AliESDV0MI *v0 = (AliESDV0MI*)tpcv0s->At(i);
+ quality[i] = 1./(1.00001-v0->GetPointAngle()); //base point angle
+ // quality[i] /= (0.5+v0->GetDist2());
+ // quality[i] *= v0->GetChi2After(); //density factor
+ Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
+ Int_t index0 = v0->GetIndex(0);
+ Int_t index1 = v0->GetIndex(1);
+ AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
+ AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
+ if (track0->fTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
+ if (track0->fTPCr[4]>0.9||track1->fTPCr[4]>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
+ }
+
+ TMath::Sort(ncandidates,quality,indexes,kTRUE);
+ //
+ //
+ for (Int_t i=0;i<ncandidates;i++){
+ AliESDV0MI * v0 = (AliESDV0MI*)tpcv0s->At(indexes[i]);
+ if (!v0) continue;
+ Int_t index0 = v0->GetIndex(0);
+ Int_t index1 = v0->GetIndex(1);
+ AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
+ AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
+ if (!track0||!track1) {
+ printf("Bug\n");
+ continue;
+ }
+ Bool_t accept =kTRUE; //default accept
+ Int_t *v0indexes0 = track0->GetV0Indexes();
+ Int_t *v0indexes1 = track1->GetV0Indexes();
+ //
+ Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
+ Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
+ if (v0indexes0[1]!=0) order0 =2;
+ if (v0indexes1[1]!=0) order1 =2;
+ //
+ if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
+ if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
+ //
+ AliESDV0MI * v02 = v0;
+ if (accept){
+ v0->SetOrder(0,order0);
+ v0->SetOrder(1,order1);
+ v0->SetOrder(1,order0+order1);
+ Int_t index = esd->AddV0MI(v0);
+ v02 = esd->GetV0MI(index);
+ v0indexes0[order0]=index;
+ v0indexes1[order1]=index;
+ naccepted++;
+ }
+ {
+ Int_t eventNr = esd->GetEventNumber();
+ if (AliTPCReconstructor::StreamLevel()>0)
+ cstream<<"V02"<<
+ "Event="<<eventNr<<
+ "vertex.="<<v0<<
+ "vertex2.="<<v02<<
+ "Tr0.="<<track0<<
+ "lab0="<<track0->fLab<<
+ "Tr1.="<<track1<<
+ "lab1="<<track1->fLab<<
+ "dca0="<<dca[index0]<<
+ "dca1="<<dca[index1]<<
+ "order0="<<order0<<
+ "order1="<<order1<<
+ "accept="<<accept<<
+ "quality="<<quality[i]<<
+ "pulldca0="<<pulldca[index0]<<
+ "pulldca1="<<pulldca[index1]<<
+ "index="<<i<<
+ "\n";
+ }
+ }
+
+
+ //
+ //
+ delete []quality;
+ delete []indexes;
+//
+ delete [] isPrim;
+ delete [] pulldca;
+ delete [] pulldcaz;
+ delete [] pulldcar;
+ delete [] cdcar;
+ delete [] sdcar;
+ delete [] dca;
+ //
+ delete[] z0;
+ delete[] alpha;
+ delete[] sign;
+ delete[] helixes;
+ printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
+ timer.Print();
+}
+
+Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink)
+{
+ //
+ // refit kink towards to the vertex
+ //
+ //
+ Int_t row0 = GetRowNumber(kink.GetR());
+ FollowProlongation(mother,0);
+ mother.Reset(kFALSE);
+ //
+ FollowProlongation(daughter,row0);
+ daughter.Reset(kFALSE);
+ FollowBackProlongation(daughter,158);
+ daughter.Reset(kFALSE);
+ Int_t first = TMath::Max(row0-20,30);
+ Int_t last = TMath::Min(row0+20,140);
+ //
+ const Int_t kNdiv =5;
+ AliTPCseed param0[kNdiv]; // parameters along the track
+ AliTPCseed param1[kNdiv]; // parameters along the track
+ AliESDkink kinks[kNdiv]; // corresponding kink parameters
+ //
+ Int_t rows[kNdiv];
+ for (Int_t irow=0; irow<kNdiv;irow++){
+ rows[irow] = first +((last-first)*irow)/(kNdiv-1);
+ }
+ // store parameters along the track
+ //
+ for (Int_t irow=0;irow<kNdiv;irow++){
+ FollowBackProlongation(mother, rows[irow]);
+ FollowProlongation(daughter,rows[kNdiv-1-irow]);
+ new(¶m0[irow]) AliTPCseed(mother);
+ new(¶m1[kNdiv-1-irow]) AliTPCseed(daughter);
+ }
+ //
+ // define kinks
+ for (Int_t irow=0; irow<kNdiv-1;irow++){
+ if (param0[irow].fN<kNdiv||param1[irow].fN<kNdiv) continue;
+ kinks[irow].SetMother(param0[irow]);
+ kinks[irow].SetDaughter(param1[irow]);
+ kinks[irow].Update();
+ }
+ //
+ // choose kink with best "quality"
+ Int_t index =-1;
+ Double_t mindist = 10000;
+ for (Int_t irow=0;irow<kNdiv;irow++){
+ if (param0[irow].fN<20||param1[irow].fN<20) continue;
+ if (TMath::Abs(kinks[irow].GetR())>240.) continue;
+ if (TMath::Abs(kinks[irow].GetR())<100.) continue;
+ //
+ Float_t normdist = TMath::Abs(param0[irow].fX-kinks[irow].GetR())*(0.1+kink.GetDistance());
+ normdist/= (param0[irow].fN+param1[irow].fN+40.);
+ if (normdist < mindist){
+ mindist = normdist;
+ index = irow;
+ }
+ }
+ //
+ if (index==-1) return 0;
+ //
+ //
+ param0[index].Reset(kTRUE);
+ FollowProlongation(param0[index],0);
+ //
+ new (&mother) AliTPCseed(param0[index]);
+ new (&daughter) AliTPCseed(param1[index]); // daughter in vertex
+ //
+ kink.SetMother(mother);
+ kink.SetDaughter(daughter);
+ kink.Update();
+ kink.SetTPCRow0(GetRowNumber(kink.GetR()));
+ kink.SetTPCncls(param0[index].fN,0);
+ kink.SetTPCncls(param1[index].fN,1);
+ kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
+ kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
+ mother.SetLabel(kink.GetLabel(0));
+ daughter.SetLabel(kink.GetLabel(1));
+
+ return 1;
+}
+
+
+void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
+ //
+ // update Kink quality information for mother after back propagation
+ //
+ if (seed->GetKinkIndex(0)>=0) return;
+ for (Int_t ikink=0;ikink<3;ikink++){
+ Int_t index = seed->GetKinkIndex(ikink);
+ if (index>=0) break;
+ index = TMath::Abs(index)-1;
+ AliESDkink * kink = fEvent->GetKink(index);
+ //kink->fTPCdensity2[0][0]=-1;
+ //kink->fTPCdensity2[0][1]=-1;
+ kink->SetTPCDensity2(-1,0,0);
+ kink->SetTPCDensity2(1,0,1);
+ //
+ Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
+ if (row0<15) row0=15;
+ //
+ Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
+ if (row1>145) row1=145;
+ //
+ Int_t found,foundable,shared;
+ seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
+ if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),0,0);
+ seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
+ if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),0,1);
+ }
+
+}
+
+void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
+ //
+ // update Kink quality information for daughter after refit
+ //
+ if (seed->GetKinkIndex(0)<=0) return;
+ for (Int_t ikink=0;ikink<3;ikink++){
+ Int_t index = seed->GetKinkIndex(ikink);
+ if (index<=0) break;
+ index = TMath::Abs(index)-1;
+ AliESDkink * kink = fEvent->GetKink(index);
+ kink->SetTPCDensity2(-1,1,0);
+ kink->SetTPCDensity2(-1,1,1);
+ //
+ Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
+ if (row0<15) row0=15;
+ //
+ Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
+ if (row1>145) row1=145;
+ //
+ Int_t found,foundable,shared;
+ seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
+ if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),1,0);
+ seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
+ if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),1,1);
+ }
+
+}
+
+
+Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink)
+{
+ //
+ // check kink point for given track
+ // if return value=0 kink point not found
+ // otherwise seed0 correspond to mother particle
+ // seed1 correspond to daughter particle
+ // kink parameter of kink point
+
+ Int_t middlerow = (seed->fFirstPoint+seed->fLastPoint)/2;
+ Int_t first = seed->fFirstPoint;
+ Int_t last = seed->fLastPoint;
+ if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
+
+
+ AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
+ if (!seed1) return 0;
+ FollowProlongation(*seed1,seed->fLastPoint-20);
+ seed1->Reset(kTRUE);
+ FollowProlongation(*seed1,158);
+ seed1->Reset(kTRUE);
+ last = seed1->fLastPoint;
+ //
+ AliTPCseed *seed0 = new AliTPCseed(*seed);
+ seed0->Reset(kFALSE);
+ seed0->Reset();
+ //
+ AliTPCseed param0[20]; // parameters along the track
+ AliTPCseed param1[20]; // parameters along the track
+ AliESDkink kinks[20]; // corresponding kink parameters
+ Int_t rows[20];
+ for (Int_t irow=0; irow<20;irow++){
+ rows[irow] = first +((last-first)*irow)/19;
+ }
+ // store parameters along the track
+ //
+ for (Int_t irow=0;irow<20;irow++){
+ FollowBackProlongation(*seed0, rows[irow]);
+ FollowProlongation(*seed1,rows[19-irow]);
+ new(¶m0[irow]) AliTPCseed(*seed0);
+ new(¶m1[19-irow]) AliTPCseed(*seed1);
+ }
+ //
+ // define kinks
+ for (Int_t irow=0; irow<19;irow++){
+ kinks[irow].SetMother(param0[irow]);
+ kinks[irow].SetDaughter(param1[irow]);
+ kinks[irow].Update();
+ }
+ //
+ // choose kink with biggest change of angle
+ Int_t index =-1;
+ Double_t maxchange= 0;
+ for (Int_t irow=1;irow<19;irow++){
+ if (TMath::Abs(kinks[irow].GetR())>240.) continue;
+ if (TMath::Abs(kinks[irow].GetR())<110.) continue;
+ Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].fX));
+ if ( quality > maxchange){
+ maxchange = quality;
+ index = irow;
+ //
+ }
+ }
+ delete seed0;
+ delete seed1;
+ if (index<0) return 0;
+ //
+ Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
+ seed0 = new AliTPCseed(param0[index]);
+ seed1 = new AliTPCseed(param1[index]);
+ seed0->Reset(kFALSE);
+ seed1->Reset(kFALSE);
+ seed0->ResetCovariance();
+ seed1->ResetCovariance();
+ FollowProlongation(*seed0,0);
+ FollowBackProlongation(*seed1,158);
+ new (&mother) AliTPCseed(*seed0); // backup mother at position 0
+ seed0->Reset(kFALSE);
+ seed1->Reset(kFALSE);
+ seed0->ResetCovariance();
+ seed1->ResetCovariance();
+ //
+ first = TMath::Max(row0-20,0);
+ last = TMath::Min(row0+20,158);
+ //
+ for (Int_t irow=0; irow<20;irow++){
+ rows[irow] = first +((last-first)*irow)/19;
+ }
+ // store parameters along the track
+ //
+ for (Int_t irow=0;irow<20;irow++){
+ FollowBackProlongation(*seed0, rows[irow]);
+ FollowProlongation(*seed1,rows[19-irow]);
+ new(¶m0[irow]) AliTPCseed(*seed0);
+ new(¶m1[19-irow]) AliTPCseed(*seed1);
+ }
+ //
+ // define kinks
+ for (Int_t irow=0; irow<19;irow++){
+ kinks[irow].SetMother(param0[irow]);
+ kinks[irow].SetDaughter(param1[irow]);
+ // param0[irow].Dump();
+ //param1[irow].Dump();
+ kinks[irow].Update();
+ }
+ //
+ // choose kink with biggest change of angle
+ index =-1;
+ maxchange= 0;
+ for (Int_t irow=0;irow<20;irow++){
+ if (TMath::Abs(kinks[irow].GetR())>250.) continue;
+ if (TMath::Abs(kinks[irow].GetR())<90.) continue;
+ Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].fX));
+ if ( quality > maxchange){
+ maxchange = quality;
+ index = irow;
+ //
+ }
+ }
+ //
+ //
+ if (index==-1 || param0[index].fN+param1[index].fN<100){
+ delete seed0;
+ delete seed1;
+ return 0;
+ }
+ // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
+
+ kink.SetMother(param0[index]);
+ kink.SetDaughter(param1[index]);
+ kink.Update();
+ row0 = GetRowNumber(kink.GetR());
+ kink.SetTPCRow0(row0);
+ kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
+ kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
+ kink.SetIndex(-10,0);
+ kink.SetIndex(int(param0[index].fN+param1[index].fN),1);
+ kink.SetTPCncls(param0[index].fN,0);
+ kink.SetTPCncls(param1[index].fN,1);
+ //
+ //
+ // new (&mother) AliTPCseed(param0[index]);
+ new (&daughter) AliTPCseed(param1[index]);
+ daughter.SetLabel(kink.GetLabel(1));
+ param0[index].Reset(kTRUE);
+ FollowProlongation(param0[index],0);
+ new (&mother) AliTPCseed(param0[index]);
+ mother.SetLabel(kink.GetLabel(0));
+ delete seed0;
+ delete seed1;
+ //
+ return 1;
+}
+
+
+
+
+AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
+{
+ //
+ // reseed - refit - track
+ //
+ Int_t first = 0;
+ // Int_t last = fSectors->GetNRows()-1;
+ //
+ if (fSectors == fOuterSec){
+ first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows());
+ //last =
+ }
+ else
+ first = t->fFirstPoint;
+ //
+ AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
+ FollowBackProlongation(*t,fSectors->GetNRows()-1);
+ t->Reset(kFALSE);
+ FollowProlongation(*t,first);
+ return seed;
+}
+
+
+
+
+
+
+
+//_____________________________________________________________________________