+ Double_t alpha, cs,sn, xx2,yy2;
+ //
+ alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
+ cs = TMath::Cos(alpha);
+ sn = TMath::Sin(alpha);
+ xx2= xyz[1][0]*cs-xyz[1][1]*sn;
+ yy2= xyz[1][0]*sn+xyz[1][1]*cs;
+ xyz[1][0] = xx2;
+ xyz[1][1] = yy2;
+ //
+ alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
+ cs = TMath::Cos(alpha);
+ sn = TMath::Sin(alpha);
+ xx2= xyz[0][0]*cs-xyz[0][1]*sn;
+ yy2= xyz[0][0]*sn+xyz[0][1]*cs;
+ xyz[0][0] = xx2;
+ xyz[0][1] = yy2;
+ //
+ //
+ //
+ Double_t x[5],c[15];
+ //
+ x[0]=xyz[2][1];
+ x[1]=xyz[2][2];
+ x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
+ x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
+ x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
+ //
+ Double_t sy =0.1, sz =0.1;
+ //
+ Double_t sy1=0.2, sz1=0.2;
+ Double_t sy2=0.2, sz2=0.2;
+ Double_t sy3=0.2;
+ //
+ Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
+ Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
+ Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
+ Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
+ Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
+ Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
+ //
+ Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
+ Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
+ Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
+ Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
+
+
+ c[0]=sy1;
+ c[1]=0.; c[2]=sz1;
+ c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
+ c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
+ c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
+ c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
+ c[13]=f30*sy1*f40+f32*sy2*f42;
+ c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
+
+ // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
+ AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
+ seed->SetLastPoint(row[2]);
+ seed->SetFirstPoint(row[2]);
+ for (Int_t i=row[0];i<row[2];i++){
+ seed->SetClusterIndex(i, track->GetClusterIndex(i));
+ }
+
+ return seed;
+}
+
+
+
+void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
+{
+ //
+ // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
+ // USES MC LABELS
+ // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
+ //
+ // Two reasons to have multiple find tracks
+ // 1. Curling tracks can be find more than once
+ // 2. Splitted tracks
+ // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
+ // b.) Edge effect on the sector boundaries
+ //
+ //
+ // Algorithm done in 2 phases - because of CPU consumption
+ // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
+ //
+ // Algorihm for curling tracks sign:
+ // 1 phase -makes a very rough fast cuts to minimize combinatorics
+ // a.) opposite sign
+ // b.) one of the tracks - not pointing to the primary vertex -
+ // c.) delta tan(theta)
+ // d.) delta phi
+ // 2 phase - calculates DCA between tracks - time consument
+
+ //
+ // fast cuts
+ //
+ // General cuts - for splitted tracks and for curling tracks
+ //
+ const Float_t kMaxdPhi = 0.2; // maximal distance in phi
+ //
+ // Curling tracks cuts
+ //
+ //
+ //
+ //
+ Int_t nentries = array->GetEntriesFast();
+ AliHelix *helixes = new AliHelix[nentries];
+ Float_t *xm = new Float_t[nentries];
+ Float_t *dz0 = new Float_t[nentries];
+ Float_t *dz1 = new Float_t[nentries];
+ //
+ //
+ TStopwatch timer;
+ timer.Start();
+ //
+ // Find track COG in x direction - point with best defined parameters
+ //
+ for (Int_t i=0;i<nentries;i++){
+ AliTPCseed* track = (AliTPCseed*)array->At(i);
+ if (!track) continue;
+ track->SetCircular(0);
+ new (&helixes[i]) AliHelix(*track);
+ Int_t ncl=0;
+ xm[i]=0;
+ Float_t dz[2];
+ track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
+ dz0[i]=dz[0];
+ dz1[i]=dz[1];
+ for (Int_t icl=0; icl<160; icl++){
+ AliTPCclusterMI * cl = track->GetClusterPointer(icl);
+ if (cl) {
+ xm[i]+=cl->GetX();
+ ncl++;
+ }
+ }
+ if (ncl>0) xm[i]/=Float_t(ncl);
+ }
+ //
+ for (Int_t i0=0;i0<nentries;i0++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
+ if (!track0) continue;
+ Float_t xc0 = helixes[i0].GetHelix(6);
+ Float_t yc0 = helixes[i0].GetHelix(7);
+ Float_t r0 = helixes[i0].GetHelix(8);
+ Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
+ Float_t fi0 = TMath::ATan2(yc0,xc0);
+
+ for (Int_t i1=i0+1;i1<nentries;i1++){
+ AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
+ if (!track1) continue;
+ Int_t lab0=track0->GetLabel();
+ Int_t lab1=track1->GetLabel();
+ if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
+ //
+ Float_t xc1 = helixes[i1].GetHelix(6);
+ Float_t yc1 = helixes[i1].GetHelix(7);
+ Float_t r1 = helixes[i1].GetHelix(8);
+ Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
+ Float_t fi1 = TMath::ATan2(yc1,xc1);
+ //
+ Float_t dfi = fi0-fi1;
+ //
+ //
+ if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
+ if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
+ if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
+ //
+ // if short tracks with undefined sign
+ fi1 = -TMath::ATan2(yc1,-xc1);
+ dfi = fi0-fi1;
+ }
+ Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
+
+ //
+ // debug stream to tune "fast cuts"
+ //
+ Double_t dist[3]; // distance at X
+ Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
+ track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
+ track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
+ track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
+ for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
+
+ Float_t sum =0;
+ Float_t sums=0;
+ for (Int_t icl=0; icl<160; icl++){
+ AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
+ AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
+ if (cl0&&cl1) {
+ sum++;
+ if (cl0==cl1) sums++;
+ }
+ }
+ //
+ if (AliTPCReconstructor::StreamLevel()>0) {
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ cstream<<"Multi"<<
+ "iter="<<iter<<
+ "lab0="<<lab0<<
+ "lab1="<<lab1<<
+ "Tr0.="<<track0<< // seed0
+ "Tr1.="<<track1<< // seed1
+ "h0.="<<&helixes[i0]<<
+ "h1.="<<&helixes[i1]<<
+ //
+ "sum="<<sum<< //the sum of rows with cl in both
+ "sums="<<sums<< //the sum of shared clusters
+ "xm0="<<xm[i0]<< // the center of track
+ "xm1="<<xm[i1]<< // the x center of track
+ // General cut variables
+ "dfi="<<dfi<< // distance in fi angle
+ "dtheta="<<dtheta<< // distance int theta angle
+ //
+ "dz00="<<dz0[i0]<<
+ "dz01="<<dz0[i1]<<
+ "dz10="<<dz1[i1]<<
+ "dz11="<<dz1[i1]<<
+ "dist0="<<dist[0]<< //distance x
+ "dist1="<<dist[1]<< //distance y
+ "dist2="<<dist[2]<< //distance z
+ "mdist0="<<mdist[0]<< //distance x
+ "mdist1="<<mdist[1]<< //distance y
+ "mdist2="<<mdist[2]<< //distance z
+ //
+ "r0="<<r0<<
+ "rc0="<<rc0<<
+ "fi0="<<fi0<<
+ "fi1="<<fi1<<
+ "r1="<<r1<<
+ "rc1="<<rc1<<
+ "\n";
+ }
+ }
+ }
+ delete [] helixes;
+ delete [] xm;
+ if (AliTPCReconstructor::StreamLevel()>1) {
+ AliInfo("Time for curling tracks removal DEBUGGING MC");
+ timer.Print();
+ }
+}
+
+
+void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
+{
+ //
+ //
+ // Two reasons to have multiple find tracks
+ // 1. Curling tracks can be find more than once
+ // 2. Splitted tracks
+ // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
+ // b.) Edge effect on the sector boundaries
+ //
+ // This function tries to find tracks closed in the parametric space
+ //
+ // cut logic if distance is bigger than cut continue - Do Nothing
+ const Float_t kMaxdTheta = 0.05; // maximal distance in theta
+ const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
+ const Float_t kdelta = 40.; //delta r to calculate track distance
+ //
+ // const Float_t kMaxDist0 = 1.; // maximal distance 0
+ //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
+ //
+ /*
+ TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
+ TCut cdtheta("cdtheta","abs(dtheta)<0.05");
+ */
+ //
+ //
+ //
+ Int_t nentries = array->GetEntriesFast();
+ AliHelix *helixes = new AliHelix[nentries];
+ Float_t *xm = new Float_t[nentries];
+ //
+ //
+ TStopwatch timer;
+ timer.Start();
+ //
+ //Sort tracks according quality
+ //
+ Int_t nseed = array->GetEntriesFast();
+ Float_t * quality = new Float_t[nseed];
+ Int_t * indexes = new Int_t[nseed];
+ for (Int_t i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
+ if (!pt){
+ quality[i]=-1;
+ continue;
+ }
+ pt->UpdatePoints(); //select first last max dens points
+ Float_t * points = pt->GetPoints();
+ if (points[3]<0.8) quality[i] =-1;
+ quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
+ //prefer high momenta tracks if overlaps
+ quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
+ }
+ TMath::Sort(nseed,quality,indexes);
+
+
+ //
+ // Find track COG in x direction - point with best defined parameters
+ //
+ for (Int_t i=0;i<nentries;i++){
+ AliTPCseed* track = (AliTPCseed*)array->At(i);
+ if (!track) continue;
+ track->SetCircular(0);
+ new (&helixes[i]) AliHelix(*track);
+ Int_t ncl=0;
+ xm[i]=0;
+ for (Int_t icl=0; icl<160; icl++){
+ AliTPCclusterMI * cl = track->GetClusterPointer(icl);
+ if (cl) {
+ xm[i]+=cl->GetX();
+ ncl++;
+ }
+ }
+ if (ncl>0) xm[i]/=Float_t(ncl);
+ }
+ //
+ for (Int_t is0=0;is0<nentries;is0++){
+ Int_t i0 = indexes[is0];
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
+ if (!track0) continue;
+ Float_t xc0 = helixes[i0].GetHelix(6);
+ Float_t yc0 = helixes[i0].GetHelix(7);
+ Float_t fi0 = TMath::ATan2(yc0,xc0);
+
+ for (Int_t is1=is0+1;is1<nentries;is1++){
+ Int_t i1 = indexes[is1];
+ AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
+ if (!track1) continue;
+ //
+ Int_t dsec = TMath::Abs((track0->GetRelativeSector()%18)-(track1->GetRelativeSector()%18)); // sector distance
+ if (dsec>1 && dsec<17) continue;
+
+ if (track1->GetKinkIndexes()[0] == -track0->GetKinkIndexes()[0]) continue;
+
+ Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
+ if (TMath::Abs(dtheta)>kMaxdTheta) continue;
+ //
+ Float_t xc1 = helixes[i1].GetHelix(6);
+ Float_t yc1 = helixes[i1].GetHelix(7);
+ Float_t fi1 = TMath::ATan2(yc1,xc1);
+ //
+ Float_t dfi = fi0-fi1;
+ if (dfi>TMath::Pi()) dfi-=TMath::TwoPi(); // take care about edge effect
+ if (dfi<-TMath::Pi()) dfi+=TMath::TwoPi(); //
+ if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
+ //
+ // if short tracks with undefined sign
+ fi1 = -TMath::ATan2(yc1,-xc1);
+ dfi = fi0-fi1;
+ if (dfi>TMath::Pi()) dfi-=TMath::TwoPi(); // take care about edge effect
+ if (dfi<-TMath::Pi()) dfi+=TMath::TwoPi(); //
+ }
+ if (TMath::Abs(dfi)>kMaxdPhi) continue;
+ //
+ //
+ Float_t sum =0;
+ Float_t sums=0;
+ Float_t sum0=0;
+ Float_t sum1=0;
+ for (Int_t icl=0; icl<160; icl++){
+ Int_t index0=track0->GetClusterIndex2(icl);
+ Int_t index1=track1->GetClusterIndex2(icl);
+ Bool_t used0 = (index0>0 && !(index0&0x8000));
+ Bool_t used1 = (index1>0 && !(index1&0x8000));
+ //
+ if (used0) sum0++; // used cluster0
+ if (used1) sum1++; // used clusters1
+ if (used0&&used1) sum++;
+ if (index0==index1 && used0 && used1) sums++;
+ }
+
+ //
+ if (sums<10) continue;
+ if (sum<40) continue;
+ if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
+ //
+ Double_t dist[5][4]; // distance at X
+ Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
+
+ //
+ //
+ track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
+ track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
+ //
+ track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
+ track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
+ //
+ track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
+ for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
+ for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
+ //
+ //
+ Int_t lab0=track0->GetLabel();
+ Int_t lab1=track1->GetLabel();
+ if( AliTPCReconstructor::StreamLevel()>5){
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ cstream<<"Splitted"<<
+ "iter="<<iter<<
+ "lab0="<<lab0<<
+ "lab1="<<lab1<<
+ "Tr0.="<<track0<< // seed0
+ "Tr1.="<<track1<< // seed1
+ "h0.="<<&helixes[i0]<<
+ "h1.="<<&helixes[i1]<<
+ //
+ "sum="<<sum<< //the sum of rows with cl in both
+ "sum0="<<sum0<< //the sum of rows with cl in 0 track
+ "sum1="<<sum1<< //the sum of rows with cl in 1 track
+ "sums="<<sums<< //the sum of shared clusters
+ "xm0="<<xm[i0]<< // the center of track
+ "xm1="<<xm[i1]<< // the x center of track
+ // General cut variables
+ "dfi="<<dfi<< // distance in fi angle
+ "dtheta="<<dtheta<< // distance int theta angle
+ //
+ //
+ "dist0="<<dist[4][0]<< //distance x
+ "dist1="<<dist[4][1]<< //distance y
+ "dist2="<<dist[4][2]<< //distance z
+ "mdist0="<<mdist[0]<< //distance x
+ "mdist1="<<mdist[1]<< //distance y
+ "mdist2="<<mdist[2]<< //distance z
+ "\n";
+ }
+ delete array->RemoveAt(i1);
+ }
+ }
+ delete [] helixes;
+ delete [] xm;
+ delete [] quality;
+ delete [] indexes;
+ AliInfo("Time for splitted tracks removal");
+ timer.Print();
+}
+
+
+
+void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
+{
+ //
+ // find Curling tracks
+ // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
+ //
+ //
+ // Algorithm done in 2 phases - because of CPU consumption
+ // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
+ // see detal in MC part what can be used to cut
+ //
+ //
+ //
+ const Float_t kMaxC = 400; // maximal curvature to of the track
+ const Float_t kMaxdTheta = 0.15; // maximal distance in theta
+ const Float_t kMaxdPhi = 0.15; // maximal distance in phi
+ const Float_t kPtRatio = 0.3; // ratio between pt
+ const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
+
+ //
+ // Curling tracks cuts
+ //
+ //
+ const Float_t kMaxDeltaRMax = 40; // distance in outer radius
+ const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
+ const Float_t kMinAngle = 2.9; // angle between tracks
+ const Float_t kMaxDist = 5; // biggest distance
+ //
+ // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
+ /*
+ Fast cuts:
+ TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
+ TCut cmax("cmax","abs(Tr0.GetC())>1/400");
+ TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
+ TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
+ TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
+ //
+ TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
+ TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
+ //
+ Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
+ Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
+ //
+ Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
+
+ */
+ //
+ //
+ //
+ Int_t nentries = array->GetEntriesFast();
+ AliHelix *helixes = new AliHelix[nentries];
+ for (Int_t i=0;i<nentries;i++){
+ AliTPCseed* track = (AliTPCseed*)array->At(i);
+ if (!track) continue;
+ track->SetCircular(0);
+ new (&helixes[i]) AliHelix(*track);
+ }
+ //
+ //
+ TStopwatch timer;
+ timer.Start();
+ Double_t phase[2][2],radius[2];
+ //
+ // Find tracks
+ //
+ //
+ for (Int_t i0=0;i0<nentries;i0++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
+ if (!track0) continue;
+ if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
+ Float_t xc0 = helixes[i0].GetHelix(6);
+ Float_t yc0 = helixes[i0].GetHelix(7);
+ Float_t r0 = helixes[i0].GetHelix(8);
+ Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
+ Float_t fi0 = TMath::ATan2(yc0,xc0);
+
+ for (Int_t i1=i0+1;i1<nentries;i1++){
+ AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
+ if (!track1) continue;
+ if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
+ Float_t xc1 = helixes[i1].GetHelix(6);
+ Float_t yc1 = helixes[i1].GetHelix(7);
+ Float_t r1 = helixes[i1].GetHelix(8);
+ Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
+ Float_t fi1 = TMath::ATan2(yc1,xc1);
+ //
+ Float_t dfi = fi0-fi1;
+ //
+ //
+ if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
+ if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
+ Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
+ //
+ //
+ // FIRST fast cuts
+ if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
+ if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
+ if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
+ if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
+ if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
+ //
+ Float_t pt0 = track0->GetSignedPt();
+ Float_t pt1 = track1->GetSignedPt();
+ if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
+ if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
+ if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
+ if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
+ //
+ //
+ // Now find closest approach
+ //
+ //
+ //
+ 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 (TMath::Abs(hangles[2])<kMinAngle) continue;
+
+ 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>kMaxDist) continue;
+ // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
+ Bool_t sign =kFALSE;
+ if (hangles[2]>kMinAngle) sign =kTRUE;
+ //
+ if (sign){
+ // circular[i0] = kTRUE;
+ // circular[i1] = kTRUE;
+ if (track0->OneOverPt()<track1->OneOverPt()){
+ track0->SetCircular(track0->GetCircular()+1);
+ track1->SetCircular(track1->GetCircular()+2);
+ }
+ else{
+ track1->SetCircular(track1->GetCircular()+1);
+ track0->SetCircular(track0->GetCircular()+2);
+ }
+ }
+ if (AliTPCReconstructor::StreamLevel()>1){
+ //
+ //debug stream to tune "fine" cuts
+ Int_t lab0=track0->GetLabel();
+ Int_t lab1=track1->GetLabel();
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ cstream<<"Curling2"<<
+ "iter="<<iter<<
+ "lab0="<<lab0<<
+ "lab1="<<lab1<<
+ "Tr0.="<<track0<<
+ "Tr1.="<<track1<<
+ //
+ "r0="<<r0<<
+ "rc0="<<rc0<<
+ "fi0="<<fi0<<
+ "r1="<<r1<<
+ "rc1="<<rc1<<
+ "fi1="<<fi1<<
+ "dfi="<<dfi<<
+ "dtheta="<<dtheta<<
+ //
+ "npoints="<<npoints<<
+ "hangles0="<<hangles[0]<<
+ "hangles1="<<hangles[1]<<
+ "hangles2="<<hangles[2]<<
+ "xyz0="<<xyz0[2]<<
+ "xyzz1="<<xyz1[2]<<
+ "radius="<<radiusbest<<
+ "deltabest="<<deltabest<<
+ "phase0="<<phase[ibest][0]<<
+ "phase1="<<phase[ibest][1]<<
+ "\n";
+
+ }
+ }
+ }
+ }
+ delete [] helixes;
+ if (AliTPCReconstructor::StreamLevel()>1) {
+ AliInfo("Time for curling tracks removal");
+ timer.Print();
+ }
+}
+
+
+
+
+
+void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
+{
+ //
+ // find kinks
+ //
+ //
+
+ TObjArray *kinks= new TObjArray(10000);
+ // TObjArray *v0s= new TObjArray(10000);
+ Int_t nentries = array->GetEntriesFast();
+ AliHelix *helixes = new AliHelix[nentries];
+ Int_t *sign = new Int_t[nentries];
+ Int_t *nclusters = new Int_t[nentries];
+ Float_t *alpha = new Float_t[nentries];
+ AliKink *kink = new AliKink();
+ Int_t * usage = new Int_t[nentries];
+ Float_t *zm = new Float_t[nentries];
+ Float_t *z0 = new Float_t[nentries];
+ Float_t *fim = new Float_t[nentries];
+ Float_t *shared = new Float_t[nentries];
+ Bool_t *circular = new Bool_t[nentries];
+ Float_t *dca = new Float_t[nentries];
+ //const AliESDVertex * primvertex = esd->GetVertex();
+ //
+ // nentries = array->GetEntriesFast();
+ //
+
+ //
+ //
+ for (Int_t i=0;i<nentries;i++){
+ sign[i]=0;
+ usage[i]=0;
+ AliTPCseed* track = (AliTPCseed*)array->At(i);
+ if (!track) continue;
+ track->SetCircular(0);
+ shared[i] = kFALSE;
+ track->UpdatePoints();
+ if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
+ }
+ nclusters[i]=track->GetNumberOfClusters();
+ 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;
+ if (track->GetProlongation(x,y,z)){
+ zm[i] = z;
+ fim[i] = alpha[i]+TMath::ATan2(y,x);
+ }
+ else{
+ zm[i] = track->GetZ();
+ fim[i] = alpha[i];
+ }
+ z0[i]=1000;
+ circular[i]= kFALSE;
+ if (track->GetProlongation(0,y,z)) z0[i] = z;
+ dca[i] = track->GetD(0,0);
+ }
+ //
+ //
+ 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
+ //
+ for (Int_t i0=0;i0<nentries;i0++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
+ if (!track0) continue;
+ if (track0->GetNumberOfClusters()<40) continue;
+ if (TMath::Abs(1./track0->GetC())>200) continue;
+ for (Int_t i1=i0+1;i1<nentries;i1++){
+ AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
+ if (!track1) continue;
+ if (track1->GetNumberOfClusters()<40) continue;
+ if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
+ if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
+ if (TMath::Abs(1./track1->GetC())>200) continue;
+ if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
+ if (track1->GetTgl()*track0->GetTgl()>0) continue;
+ if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
+ if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
+ if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) 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;
+ 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 lsign =kFALSE;
+ if (hangles[2]>3.06) lsign =kTRUE;
+ //
+ if (lsign){
+ circular[i0] = kTRUE;
+ circular[i1] = kTRUE;
+ if (track0->OneOverPt()<track1->OneOverPt()){
+ track0->SetCircular(track0->GetCircular()+1);
+ track1->SetCircular(track1->GetCircular()+2);
+ }
+ else{
+ track1->SetCircular(track1->GetCircular()+1);
+ track0->SetCircular(track0->GetCircular()+2);
+ }
+ }
+ if (lsign&&AliTPCReconstructor::StreamLevel()>1){
+ //debug stream
+ Int_t lab0=track0->GetLabel();
+ Int_t lab1=track1->GetLabel();
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ cstream<<"Curling"<<
+ "lab0="<<lab0<<
+ "lab1="<<lab1<<
+ "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);
+ if (track0==0) {
+ AliInfo("seed==0");
+ continue;
+ }
+ 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,ishared;
+ track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
+ if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
+ track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
+ if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
+ //
+ track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
+ if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
+ track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,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];
+ fkParam->Transform0to1(x,index);
+ fkParam->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->GetP();
+ 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->GetP()<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->GetSigmaSnp2()+track0->GetSigmaTgl2();
+ criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
+ 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->GetClusterPointer(row)){
+ AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
+ shapesum+=point->GetSigmaY()+point->GetSigmaZ();
+ sum++;
+ }
+ if (ktrack1->GetClusterPointer(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 AliKink;
+ 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;
+ AliKink *kinkl = (AliKink*)kinks->At(i);
+ //
+ // refit kinks towards vertex
+ //
+ Int_t index0 = kinkl->GetIndex(0);
+ Int_t index1 = kinkl->GetIndex(1);
+ AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
+ AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
+ //
+ Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
+ //
+ // Refit Kink under if too small angle
+ //
+ if (kinkl->GetAngle(2)<0.05){
+ kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
+ Int_t row0 = kinkl->GetTPCRow0();
+ Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
+ //
+ //
+ Int_t last = row0-drow;
+ if (last<40) last=40;
+ if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
+ AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
+ //
+ //
+ Int_t first = row0+drow;
+ if (first>130) first=130;
+ if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
+ AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
+ //
+ if (seed0 && seed1){
+ kinkl->SetStatus(1,8);
+ if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
+ row0 = GetRowNumber(kinkl->GetR());
+ sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
+ mothers[i] = *seed0;
+ daughters[i] = *seed1;
+ }
+ else{
+ delete kinks->RemoveAt(i);
+ if (seed0) delete seed0;
+ if (seed1) delete seed1;
+ continue;
+ }
+ if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
+ delete kinks->RemoveAt(i);
+ if (seed0) delete seed0;
+ if (seed1) delete seed1;
+ continue;
+ }
+ //
+ delete seed0;
+ delete seed1;
+ }
+ //
+ if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->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++){
+ AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
+ if (!kink0) continue;
+ //
+ for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
+ if (!kink0) continue;
+ AliKink * kink1 = (AliKink*) 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.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
+ both++;
+ bothm++;
+ if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
+ same++;
+ samem++;
+ }
+ }
+ }
+
+ for (Int_t i=row0;i<158;i++){
+ if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
+ both++;
+ bothd++;
+ if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(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.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
+ Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
+ 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++){
+ AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
+ if (!kinkl) continue;
+ kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
+ Int_t index0 = kinkl->GetIndex(0);
+ Int_t index1 = kinkl->GetIndex(1);
+ if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
+ kinkl->SetMultiple(usage[index0],0);
+ kinkl->SetMultiple(usage[index1],1);
+ if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
+ if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
+ if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
+ if (circular[index0]||(circular[index1]&&kinkl->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(kinkl);
+ //
+ //
+ if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
+ if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
+ *ktrack0 = mothers[indexes[i]];
+ *ktrack1 = daughters[indexes[i]];
+ }
+ }
+ //
+ ktrack0->SetKinkIndex(usage[index0],-(index+1));
+ ktrack1->SetKinkIndex(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->GetKinkIndex(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 ishared=0;
+ Int_t all =0;
+ for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
+ if (track0->GetClusterPointer(icl)!=0){
+ all++;
+ if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
+ }
+ }
+ if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
+ delete array->RemoveAt(i);
+ continue;
+ }
+ //
+ if (track0->GetKinkIndex(0)!=0) continue;
+ if (track0->GetNumberOfClusters()<80) continue;
+
+ AliTPCseed *pmother = new AliTPCseed();
+ AliTPCseed *pdaughter = new AliTPCseed();
+ AliKink *pkink = new AliKink;
+
+ AliTPCseed & mother = *pmother;
+ AliTPCseed & daughter = *pdaughter;
+ AliKink & kinkl = *pkink;
+ if (CheckKinkPoint(track0,mother,daughter, kinkl)){
+ if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<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= kinkl.GetTPCRow0();
+ if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
+ delete pmother;
+ delete pdaughter;
+ delete pkink;
+ continue;
+ }
+ //
+ Int_t index = esd->AddKink(&kinkl);
+ mother.SetKinkIndex(0,-(index+1));
+ daughter.SetKinkIndex(0,index+1);
+ if (mother.GetNumberOfClusters()>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.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
+ }
+ //
+ for (Int_t icl=row0;icl<158;icl++) {
+ if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(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(const TObjArray * array, AliESDEvent *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->GetC())*(100*track->GetC()));
+ if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
+ cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-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->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
+ if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut