From 81e97e0da7f50c2a0bc70f39a14d845607ba361a Mon Sep 17 00:00:00 2001 From: hristov Date: Fri, 29 Apr 2005 03:23:08 +0000 Subject: [PATCH] New version of the V0 finder (M.Ivanov) --- ITS/AliITStrackV2.cxx | 24 +- ITS/AliITStrackerMI.cxx | 1336 ++++++++++++++++++++--------------- ITS/AliITStrackerMI.h | 15 +- ITS/AliITStrackerV2.cxx | 4 + STEER/AliESDComparisonMI.C | 49 +- STEER/AliESDV0MI.cxx | 51 +- STEER/AliESDV0MI.h | 16 +- STEER/AliESDtrack.cxx | 34 +- STEER/AliESDtrack.h | 9 +- STEER/AliESDv0.h | 2 +- STEER/AliGenInfo.C | 24 +- STEER/AliGenInfo.h | 1 + STEER/AliHelix.cxx | 25 +- STEER/AliHelix.h | 3 +- STEER/AliMagFMaps.cxx | 2 +- STEER/TTreeStream.cxx | 4 +- TPC/AliTPCseed.cxx | 922 ++++++++++++++++++++++++ TPC/AliTPCseed.h | 120 ++++ TPC/AliTPCtrack.cxx | 5 +- TPC/AliTPCtrack.h | 5 +- TPC/AliTPCtrackerMI.cxx | 1367 +++++++++++++----------------------- TPC/AliTPCtrackerMI.h | 90 +-- TPC/libTPCrec.pkg | 2 +- 23 files changed, 2498 insertions(+), 1612 deletions(-) create mode 100644 TPC/AliTPCseed.cxx create mode 100644 TPC/AliTPCseed.h diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx index 43a394e43ca..af3a54b691a 100644 --- a/ITS/AliITStrackV2.cxx +++ b/ITS/AliITStrackV2.cxx @@ -107,9 +107,9 @@ AliKalmanTrack() { } fESDtrack=&t; + // if (!Invariant()) throw "AliITStrackV2: conversion failed !\n"; for(Int_t i=0; i<4; i++) fdEdxSample[i]=0; - if (!Invariant()) throw "AliITStrackV2: conversion failed !\n"; } void AliITStrackV2::UpdateESDtrack(ULong_t flags) const { @@ -278,10 +278,12 @@ Int_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) { //------------------------------------------------------------------ Double_t x1=fX, x2=xk, dx=x2-x1; Double_t f1=fP2, f2=f1 + fP4*dx; - if (TMath::Abs(f2) >= 0.9999) { - Int_t n=GetNumberOfClusters(); - if (n>kWARN) - Warning("PropagateTo","Propagation failed !\n",n); + if (TMath::Abs(f2) >= 0.98) { + // MI change - don't propagate highly inclined tracks + // covariance matrix distorted + // Int_t n=GetNumberOfClusters(); + // if (n>kWARN) + // Warning("PropagateTo","Propagation failed !\n",n); return 0; } @@ -504,12 +506,16 @@ Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) { { Double_t dx=xk-fX; Double_t f1=fP2, f2=f1 + fP4*dx; - if (TMath::Abs(f2) >= 0.9999) { - Int_t n=GetNumberOfClusters(); - if (n>kWARN) - Warning("Propagate","Propagation failed (%d) !\n",n); + if (TMath::Abs(f2) >= 0.98) { + // don't propagate highly inclined tracks MI return 0; } + // if (TMath::Abs(f2) >= 0.9999) { +// Int_t n=GetNumberOfClusters(); +// if (n>kWARN) +// Warning("Propagate","Propagation failed (%d) !\n",n); +// return 0; +// } Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2); fX=xk; diff --git a/ITS/AliITStrackerMI.cxx b/ITS/AliITStrackerMI.cxx index 998b8c6f25a..647ef97b812 100644 --- a/ITS/AliITStrackerMI.cxx +++ b/ITS/AliITStrackerMI.cxx @@ -100,6 +100,9 @@ AliITStrackerMI::AliITStrackerMI(const AliITSgeom *geom) : AliTracker() { for (Int_t i=0;i<100000;i++){ fBestTrackIndex[i]=0; } + // + fDebugStreamer = new TTreeSRedirector("ITSdebug.root"); + } AliITStrackerMI::~AliITStrackerMI() @@ -108,6 +111,10 @@ AliITStrackerMI::~AliITStrackerMI() //destructor // if (fCoeficients) delete []fCoeficients; + if (fDebugStreamer) { + //fDebugStreamer->Close(); + delete fDebugStreamer; + } } void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) { @@ -231,6 +238,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) { //-------------------------------------------------------------------- TObjArray itsTracks(15000); fOriginal.Clear(); + fEsd = event; // store pointer to the esd {/* Read ESD tracks */ Int_t nentr=event->GetNumberOfTracks(); Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr); @@ -256,25 +264,30 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) { // write expected q t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.); - if (TMath::Abs(t->fD[0])>10) { - delete t; - continue; - } - - if (TMath::Abs(vdist)>20) { - delete t; - continue; - } - if (TMath::Abs(1/t->Get1Pt())<0.120) { - delete t; - continue; + if (esd->GetV0Index(0)>0 && t->fD[0]<30){ + //track - can be V0 according to TPC } - - if (CorrectForDeadZoneMaterial(t)!=0) { - //Warning("Clusters2Tracks", - // "failed to correct for the material in the dead zone !\n"); - delete t; - continue; + else{ + if (TMath::Abs(t->fD[0])>10) { + delete t; + continue; + } + + if (TMath::Abs(vdist)>20) { + delete t; + continue; + } + if (TMath::Abs(1/t->Get1Pt())<0.120) { + delete t; + continue; + } + + if (CorrectForDeadZoneMaterial(t)!=0) { + //Warning("Clusters2Tracks", + // "failed to correct for the material in the dead zone !\n"); + delete t; + continue; + } } t->fReconstructed = kFALSE; itsTracks.AddLast(t); @@ -286,6 +299,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) { fOriginal.Sort(); Int_t nentr=itsTracks.GetEntriesFast(); fTrackHypothesys.Expand(nentr); + fBestHypothesys.Expand(nentr); MakeCoeficients(nentr); Int_t ntrk=0; for (fPass=0; fPass<2; fPass++) { @@ -303,9 +317,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) { fI = 6; ResetTrackToFollow(*t); ResetBestTrack(); - - FollowProlongationTree(t,i); - + FollowProlongationTree(t,i,fConstraint[fPass]); SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change // @@ -336,7 +348,8 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) { } //GetBestHypothesysMIP(itsTracks); - FindV0(event); + UpdateTPCV0(event); + FindV02(event); fAfterV0 = kTRUE; //GetBestHypothesysMIP(itsTracks); // @@ -350,6 +363,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) { } fTrackHypothesys.Delete(); + fBestHypothesys.Delete(); fOriginal.Clear(); delete []fCoeficients; fCoeficients=0; @@ -483,9 +497,6 @@ Int_t AliITStrackerMI::RefitInward(AliESD *event) { if (fTrackToFollow.Propagate(fv+a,xv)) { fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit); - Float_t d=fTrackToFollow.GetD(GetX(),GetY()); - Float_t z=fTrackToFollow.GetZ()-GetZ(); - fTrackToFollow.GetESDtrack()->SetImpactParameters(d,z); //UseClusters(&fTrackToFollow); { AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ()); @@ -520,12 +531,40 @@ AliCluster *AliITStrackerMI::GetCluster(Int_t index) const { } -void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex) +void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain) { //-------------------------------------------------------------------- // Follow prolongation tree //-------------------------------------------------------------------- + // + AliESDtrack * esd = otrack->fESDtrack; + if (esd->GetV0Index(0)>0){ + // + // TEMPORARY SOLLUTION: map V0 indexes to point to proper track + // mapping of esd track is different as its track in Containers + // Need something more stable + // Indexes are set back againg to the ESD track indexes in UpdateTPCV0 + for (Int_t i=0;i<3;i++){ + Int_t index = esd->GetV0Index(i); + if (index==0) break; + AliESDV0MI * vertex = fEsd->GetV0MI(index); + if (vertex->GetStatus()<0) continue; // rejected V0 + // + if (esd->GetSign()>0) { + vertex->SetIndex(0,esdindex); + } + else{ + vertex->SetIndex(1,esdindex); + } + } + } + TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex); + if (!bestarray){ + bestarray = new TObjArray(5); + fBestHypothesys.AddAt(bestarray,esdindex); + } + // //setup tree of the prolongations // static AliITStrackMI tracks[7][100]; @@ -540,7 +579,7 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin otrack->fNSkipped=0; new (&(tracks[6][0])) AliITStrackMI(*otrack); ntracks[6]=1; - nindexes[6][0]=0; + for (Int_t i=0;i<7;i++) nindexes[i][0]=0; // // // follow prolongations @@ -626,7 +665,7 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin // Double_t msz=1./((currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer])); Double_t msy=1./((currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer])); - if (fConstraint[fPass]){ + if (constrain){ msy/=60; msz/=60.; } else{ @@ -687,8 +726,8 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin Double_t x0; Double_t d=layer.GetThickness(updatetrack->GetY(),updatetrack->GetZ(),x0); updatetrack->CorrectForMaterial(d,x0); - if (fConstraint[fPass]) { - updatetrack->fConstrain = fConstraint[fPass]; + if (constrain) { + updatetrack->fConstrain = constrain; fI = ilayer; Double_t d=GetEffectiveThickness(0,0); //Think of this !!!! Double_t xyz[]={GetX(),GetY(),GetZ()}; @@ -705,8 +744,8 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin ntracks[ilayer]++; } // create new hypothesy } // loop over possible cluster prolongation - // if (fConstraint[fPass]&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0){ - if (fConstraint[fPass]&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){ + // if (constrain&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0){ + if (constrain&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){ AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1); vtrack->fClIndex[ilayer]=0; fI = ilayer; @@ -718,7 +757,7 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin ntracks[ilayer]++; } - if (fConstraint[fPass]&&itrack<1&&TMath::Abs(currenttrack1.fP3)>1.1){ //big theta -- for low mult. runs + if (constrain&&itrack<1&&TMath::Abs(currenttrack1.fP3)>1.1){ //big theta -- for low mult. runs AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1); vtrack->fClIndex[ilayer]=0; fI = ilayer; @@ -742,8 +781,8 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin if ( normalizedchi2[itrack]<3+0.5*ilayer) golds++; if (ilayer>4) accepted++; else{ - if ( fConstraint[fPass] && normalizedchi2[itrack]90) ntracks[ilayer]=90; } //loop over layers //printf("%d\t%d\t%d\t%d\t%d\t%d\n",ntracks[0],ntracks[1],ntracks[2],ntracks[3],ntracks[4],ntracks[5]); - Int_t max = fConstraint[fPass]? 20: 5; + Int_t max = constrain? 20: 5; for (Int_t i=0;i7.)continue; + if (!constrain&&track.fNormChi2[0]>7.)continue; AddTrackHypothesys(new AliITStrackMI(track), esdindex); } for (Int_t i=0;i7.)continue; - if (fConstraint[fPass]) track.fNSkipped+=1; - if (!fConstraint[fPass]) { + if (!constrain&&track.fNormChi2[1]>7.)continue; + if (constrain) track.fNSkipped+=1; + if (!constrain) { track.fD[0] = track.GetD(GetX(),GetY()); track.fNSkipped+=4./(4.+8.*TMath::Abs(track.fD[0])); if (track.fN+track.fNDeadZone+track.fNSkipped>6) { @@ -776,13 +815,13 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin } //} - if (!fConstraint[fPass]){ + if (!constrain){ for (Int_t i=0;i7.)continue; - if (fConstraint[fPass]) track.fNSkipped+=2; - if (!fConstraint[fPass]){ + if (!constrain&&track.fNormChi2[2]>7.)continue; + if (constrain) track.fNSkipped+=2; + if (!constrain){ track.fD[0] = track.GetD(GetX(),GetY()); track.fNSkipped+= 7./(7.+8.*TMath::Abs(track.fD[0])); if (track.fN+track.fNDeadZone+track.fNSkipped>6) { @@ -792,6 +831,76 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin AddTrackHypothesys(new AliITStrackMI(track), esdindex); } } + + if (!constrain){ + // + // register best tracks - important for V0 finder + // + for (Int_t ilayer=0;ilayer<5;ilayer++){ + if (ntracks[ilayer]==0) continue; + AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]]; + if (track.GetNumberOfClusters()<1) continue; + CookLabel(&track,0); + bestarray->AddAt(new AliITStrackMI(track),ilayer); + } + } + // + // update TPC V0 information + // + if (otrack->fESDtrack->GetV0Index(0)>0){ + Float_t fprimvertex[3]={GetX(),GetY(),GetZ()}; + for (Int_t i=0;i<3;i++){ + Int_t index = otrack->fESDtrack->GetV0Index(i); + if (index==0) break; + AliESDV0MI * vertex = fEsd->GetV0MI(index); + if (vertex->GetStatus()<0) continue; // rejected V0 + // + if (otrack->fP4>0) { + vertex->SetIndex(0,esdindex); + } + else{ + vertex->SetIndex(1,esdindex); + } + //find nearest layer with track info + Int_t nearestold = GetNearestLayer(vertex->GetXrp()); + Int_t nearest = nearestold; + for (Int_t ilayer =nearest;ilayer<8;ilayer++){ + if (ntracks[nearest]==0){ + nearest = ilayer; + } + } + // + AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]]; + if (nearestold<5&&nearest<5){ + Bool_t accept = track.fNormChi2[nearest]<10; + if (accept){ + if (track.fP4>0) { + vertex->SetP(track); + vertex->Update(fprimvertex); + // vertex->SetIndex(0,track.fESDtrack->GetID()); + if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex); + }else{ + vertex->SetM(track); + vertex->Update(fprimvertex); + //vertex->SetIndex(1,track.fESDtrack->GetID()); + if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex); + } + vertex->SetStatus(vertex->GetStatus()+1); + }else{ + // vertex->SetStatus(-2); // reject V0 - not enough clusters + } + } + // if (nearestold>3){ +// Int_t indexlayer = (ntracks[0]>0)? 0:1; +// if (ntracks[indexlayer]>0){ +// AliITStrackMI & track= tracks[indexlayer][nindexes[indexlayer][0]]; +// if (track.GetNumberOfClusters()>4&&track.fNormChi2[indexlayer]<4){ +// vertex->SetStatus(-1); // reject V0 - clusters before +// } +// } +// } + } + } } @@ -2389,7 +2498,7 @@ AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI if (!backtrack->Improve(0,xyzv,ersv)) continue; if (!backtrack->PropagateToVertex()) continue; backtrack->ResetCovariance(); - if (!backtrack->Improve(0,xyzv,ersv)) continue; + if (!backtrack->Improve(0,xyzv,ersv)) continue; }else{ backtrack->ResetCovariance(); } @@ -2975,7 +3084,6 @@ void AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Fl // sigmarfi = 0.004+1.4 *TMath::Abs(track->fP4)+332.*track->fP4*track->fP4; sigmaz = 0.011+4.37*TMath::Abs(track->fP4); - } @@ -3064,607 +3172,685 @@ void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const +Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{ + // + //Get nearest upper layer close to the point xr. + // rough approximation + // + const Float_t radiuses[6]={4,6.5,15.03,24.,38.5,43.7}; + Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]); + Int_t res =6; + for (Int_t i=0;i<6;i++){ + if (radiusGetNumberOfV0MIs(); + Int_t nitstracks = fTrackHypothesys.GetEntriesFast(); + for (Int_t i=0;iGetV0MI(i); + Int_t ip = vertex->GetIndex(0); + Int_t im = vertex->GetIndex(1); + // + TObjArray * arrayp = (ipAt(0):0; + AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0; + // + // + if (trackp){ + if (trackp->fN+trackp->fNDeadZone>5.5){ + if (trackp->fConstrain&&trackp->fChi2MIP[0]<3) vertex->SetStatus(-100); + if (!trackp->fConstrain&&trackp->fChi2MIP[0]<2) vertex->SetStatus(-100); + } + } + if (trackm){ + if (trackm->fN+trackm->fNDeadZone>5.5){ + if (trackm->fConstrain&&trackm->fChi2MIP[0]<3) vertex->SetStatus(-100); + if (!trackm->fConstrain&&trackm->fChi2MIP[0]<2) vertex->SetStatus(-100); + } + } + if (vertex->GetStatus()==-100) continue; + // + Int_t clayer = GetNearestLayer(vertex->GetXrp()); + vertex->SetNBefore(clayer); // + vertex->SetChi2Before(9*clayer); // + vertex->SetNAfter(6-clayer); // + vertex->SetChi2After(0); // + // + if (clayer >1 ){ // calculate chi2 before vertex + Float_t chi2p = 0, chi2m=0; + // + if (trackp){ + for (Int_t ilayer=0;ilayerfClIndex[ilayer]>0){ + chi2p+=trackp->fDy[ilayer]*trackp->fDy[ilayer]/(trackp->fSigmaY[ilayer]*trackp->fSigmaY[ilayer])+ + trackp->fDz[ilayer]*trackp->fDz[ilayer]/(trackp->fSigmaZ[ilayer]*trackp->fSigmaZ[ilayer]); + } + else{ + chi2p+=9; + } + } + }else{ + chi2p = 9*clayer; + } + // + if (trackm){ + for (Int_t ilayer=0;ilayerfClIndex[ilayer]>0){ + chi2m+=trackm->fDy[ilayer]*trackm->fDy[ilayer]/(trackm->fSigmaY[ilayer]*trackm->fSigmaY[ilayer])+ + trackm->fDz[ilayer]*trackm->fDz[ilayer]/(trackm->fSigmaZ[ilayer]*trackm->fSigmaZ[ilayer]); + } + else{ + chi2m+=9; + } + } + }else{ + chi2m = 9*clayer; + } + vertex->SetChi2Before(TMath::Min(chi2p,chi2m)); + if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10); // track exist before vertex + } + + if (clayer < 5 ){ // calculate chi2 after vertex + Float_t chi2p = 0, chi2m=0; + // + if (trackp&&TMath::Abs(trackp->fP3)<1.){ + for (Int_t ilayer=clayer;ilayer<6;ilayer++){ + if (trackp->fClIndex[ilayer]>0){ + chi2p+=trackp->fDy[ilayer]*trackp->fDy[ilayer]/(trackp->fSigmaY[ilayer]*trackp->fSigmaY[ilayer])+ + trackp->fDz[ilayer]*trackp->fDz[ilayer]/(trackp->fSigmaZ[ilayer]*trackp->fSigmaZ[ilayer]); + } + else{ + chi2p+=9; + } + } + }else{ + chi2p = 0; + } + // + if (trackm&&TMath::Abs(trackm->fP3)<1.){ + for (Int_t ilayer=clayer;ilayer<6;ilayer++){ + if (trackm->fClIndex[ilayer]>0){ + chi2m+=trackm->fDy[ilayer]*trackm->fDy[ilayer]/(trackm->fSigmaY[ilayer]*trackm->fSigmaY[ilayer])+ + trackm->fDz[ilayer]*trackm->fDz[ilayer]/(trackm->fSigmaZ[ilayer]*trackm->fSigmaZ[ilayer]); + } + else{ + chi2m+=9; + } + } + }else{ + chi2m = 0; + } + vertex->SetChi2After(TMath::Max(chi2p,chi2m)); + if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20); // track not found in ITS + } + } + // +} -void AliITStrackerMI::FindV0(AliESD *event) +void AliITStrackerMI::FindV02(AliESD *event) { // - // fast V0 finder + // V0 finder + // + // Cuts on DCA - R dependent + // max distance DCA between 2 tracks cut + // maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist); + // + const Float_t kMaxDist0 = 0.1; + const Float_t kMaxDist1 = 0.1; + const Float_t kMaxDist = 1; + const Float_t kMinPointAngle = 0.85; + const Float_t kMinPointAngle2 = 0.99; + const Float_t kMinR = 0.5; + const Float_t kMaxR = 220; + //const Float_t kCausality0Cut = 0.19; + //const Float_t kLikelihood01Cut = 0.25; + //const Float_t kPointAngleCut = 0.9996; + const Float_t kCausality0Cut = 0.19; + const Float_t kLikelihood01Cut = 0.45; + const Float_t kLikelihood1Cut = 0.5; + const Float_t kCombinedCut = 0.55; + + // + // + TTreeSRedirector &cstream = *fDebugStreamer; + Int_t ntracks = event->GetNumberOfTracks(); + Int_t nitstracks = fTrackHypothesys.GetEntriesFast(); + fOriginal.Expand(ntracks); + fTrackHypothesys.Expand(ntracks); + fBestHypothesys.Expand(ntracks); + // + AliHelix * helixes = new AliHelix[ntracks+2]; + TObjArray trackarray(ntracks+2); //array with tracks - with vertex constrain + TObjArray trackarrayc(ntracks+2); //array of "best tracks" - without vertex constrain + TObjArray trackarrayl(ntracks+2); //array of "longest tracks" - without vertex constrain + Bool_t * forbidden = new Bool_t [ntracks+2]; + Int_t *itsmap = new Int_t [ntracks+2]; + Float_t *dist = new Float_t[ntracks+2]; + Float_t *normdist0 = new Float_t[ntracks+2]; + Float_t *normdist1 = new Float_t[ntracks+2]; + Float_t *normdist = new Float_t[ntracks+2]; + Float_t *norm = new Float_t[ntracks+2]; + Float_t *maxr = new Float_t[ntracks+2]; + Float_t *minr = new Float_t[ntracks+2]; + Float_t *minPointAngle= new Float_t[ntracks+2]; + // + AliESDV0MI *pvertex = new AliESDV0MI; + AliITStrackMI * dummy= new AliITStrackMI; + dummy->SetLabel(0); + AliITStrackMI trackat0; //temporary track for DCA calculation // - //TTreeSRedirector cstream("itsv0.root"); - Int_t centries=0; - AliHelix * helixes = new AliHelix[30000]; - TObjArray trackarray(30000); - TObjArray trackarrayc(30000); - Float_t * dist = new Float_t[30000]; - Float_t * normdist0 = new Float_t[30000]; - Float_t * normdist1 = new Float_t[30000]; - Float_t * normdist = new Float_t[30000]; - Float_t * norm = new Float_t[30000]; - AliESDV0MI *vertexarray = new AliESDV0MI[100000]; - AliESDV0MI *pvertex = &vertexarray[0]; - AliITStrackMI * dummy=0; + Float_t primvertex[3]={GetX(),GetY(),GetZ()}; // + // make its - esd map // - Int_t entries = fTrackHypothesys.GetEntriesFast(); - for (Int_t i=0;iGetEntriesFast(); + for (Int_t itrack=0;itrackfESDtrack->GetID(); + itsmap[esdindex] = itrack; + } + // + // create its tracks from esd tracks if not done before + // + for (Int_t itrack=0;itrack=0) continue; + AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack))); + tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY()); + tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ(); + if (tpctrack->fD[0]<20 && tpctrack->fD[1]<20){ + // tracks which can reach inner part of ITS + // propagate track to outer its volume - with correction for material + CorrectForDeadZoneMaterial(tpctrack); + } + itsmap[itrack] = nitstracks; + fOriginal.AddAt(tpctrack,nitstracks); + nitstracks++; + } + // + // fill temporary arrays + // + for (Int_t itrack=0;itrackGetTrack(itrack); + Int_t itsindex = itsmap[itrack]; + AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]); + if (!original) continue; + AliITStrackMI *bestConst = 0; + AliITStrackMI *bestLong = 0; + AliITStrackMI *best = 0; + // // - // best with vertex constrain - AliITStrackMI * trackc = (AliITStrackMI*)array->At(0); - if (trackc&&trackc->fConstrain&&trackc->fN==6&&trackc->fNormChi2[0]<2.) continue; - trackc=0; + TObjArray * array = (TObjArray*) fTrackHypothesys.At(itsindex); + Int_t hentries = (array==0) ? 0 : array->GetEntriesFast(); + // Get best track with vertex constrain for (Int_t ih=0;ihAt(ih); if (!trackh->fConstrain) continue; - if (trackh->fN<6) continue; - trackc = trackh; - if (!dummy) dummy = trackc; - dummy->SetLabel(0); + if (!bestConst) bestConst = trackh; + if (trackh->fN>5.0){ + bestConst = trackh; // full track - with minimal chi2 + break; + } + if (trackh->fN+trackh->fNDeadZone<=bestConst->fN+bestConst->fNDeadZone) continue; + bestConst = trackh; break; - } - // - // best without vertex - AliITStrackMI * track = 0; + } + // Get best long track without vertex constrain and best track without vertex constrain for (Int_t ih=0;ihAt(ih); if (trackh->fConstrain) continue; - track = trackh; - break; + if (!best) best = trackh; + if (!bestLong) bestLong = trackh; + if (trackh->fN>5.0){ + bestLong = trackh; // full track - with minimal chi2 + break; + } + if (trackh->fN+trackh->fNDeadZone<=bestLong->fN+bestLong->fNDeadZone) continue; + bestLong = trackh; } - if (trackc&&track){ - if (trackc->fChi2MIP[1]<2.) continue; - if (trackc->fChi2MIP[0]<2. && trackc->fChi2MIP[1]<2.) continue; - trackarrayc.AddAt(trackc,i); - if (trackc->fN==6&&track->fN&&trackc->fNormChi2[0] < track->fNormChi2[0]-2) continue; + if (!best) { + best = original; + bestLong = original; } + trackat0 = *bestLong; + Double_t xx,yy,zz,alpha; + bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz); + alpha = TMath::ATan2(yy,xx); + trackat0.Propagate(alpha,0); + // calculate normalized distances to the vertex // + if ( bestLong->fN>3 ){ + dist[itsindex] = trackat0.fP0; + norm[itsindex] = TMath::Sqrt(trackat0.fC00); + normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]); + normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/TMath::Sqrt(trackat0.fC11)); + normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]); + } + else{ + if (bestConst&&bestConst->fN+bestConst->fNDeadZone>4.5){ + dist[itsindex] = bestConst->fD[0]; + norm[itsindex] = bestConst->fDnorm[0]; + normdist0[itsindex] = TMath::Abs(bestConst->fD[0]/norm[itsindex]); + normdist1[itsindex] = TMath::Abs(bestConst->fD[0]/norm[itsindex]); + normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]); + }else{ + dist[itsindex] = trackat0.fP0; + norm[itsindex] = TMath::Sqrt(trackat0.fC00); + normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]); + normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/TMath::Sqrt(trackat0.fC11)); + normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]); + if (TMath::Abs(trackat0.fP3)>1.1){ + if (normdist[itsindex]<10) forbidden[itsindex]=kTRUE; + } + } + } // + //----------------------------------------------------------- + //Forbid primary track candidates - // - if (track){ - dist[i] = TMath::Sqrt(track->fD[0]*track->fD[0]+track->fD[1]*track->fD[1]); - norm[i] = track->fDnorm[0]; - normdist0[i] = TMath::Abs(track->fD[0]/track->fDnorm[0]); - normdist1[i] = TMath::Abs(track->fD[1]/track->fDnorm[1]); - normdist[i] = TMath::Sqrt(normdist0[i]*normdist0[i]+normdist1[i]*normdist1[i]); - if (track->IsGoldPrimary()) continue; //primary track - if (track->fD[0]<0.02 && (track->fN+track->fNDeadZone>5.8)){ - if (normdist[i]<3.) continue; // primary track - cutoff 3 sigma - if (normdist0[i]<2.) continue; //DCA normalized cut 2 sigma + //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5"); + //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5"); + //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0"); + //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0"); + //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2"); + //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1"); + //----------------------------------------------------------- + if (bestConst){ + if (bestLong->fN<4 && bestConst->fN+bestConst->fNDeadZone>4.5) forbidden[itsindex]=kTRUE; + if (normdist[itsindex]<3 && bestConst->fN+bestConst->fNDeadZone>5.5) forbidden[itsindex]=kTRUE; + if (normdist[itsindex]<2 && bestConst->fClIndex[0]>0 && bestConst->fClIndex[1]>0 ) forbidden[itsindex]=kTRUE; + if (normdist[itsindex]<1 && bestConst->fClIndex[0]>0) forbidden[itsindex]=kTRUE; + if (normdist[itsindex]<4 && bestConst->fNormChi2[0]<2) forbidden[itsindex]=kTRUE; + if (normdist[itsindex]<5 && bestConst->fNormChi2[0]<1) forbidden[itsindex]=kTRUE; + if (bestConst->fNormChi2[0]<2.5) { + minPointAngle[itsindex]= 0.9999; + maxr[itsindex] = 10; } - trackarray.AddAt(track,i); - new (&helixes[i]) AliHelix(*track); } + // + //forbid daughter kink candidates + // + if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE; + Bool_t isElectron = kTRUE; + Double_t pid[5]; + esdtrack->GetESDpid(pid); + for (Int_t i=1;i<5;i++){ + if (pid[0]GetTPCdensity(0,10) >0.6) maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]); + if (esdtrack->GetTPCdensity(10,30)>0.6) maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]); + if (esdtrack->GetTPCdensity(20,40)>0.6) maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]); + if (esdtrack->GetTPCdensity(30,50)>0.6) maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]); + // + if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->fN<3) minr[itsindex]=100; + // + // + if (0){ + cstream<<"Track"<< + "Tr0.="<fP4>0) continue; + // + // first iterration of V0 finder + // + for (Int_t iesd0=0;iesd0fP4>0) continue; AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0); // - TObjArray * array0 = (TObjArray*)fTrackHypothesys.At(itrack0); - // - Int_t vertexes =0; - for (Int_t itrack1=0;itrack1fP4<0) continue; - AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1); - if (trackc0&&trackc1){ - if (TMath::Min(trackc0->fChi2MIP[1],trackc1->fChi2MIP[1])<2.) continue; + for (Int_t iesd1=0;iesd1fP4<0) continue; + Bool_t isGold = kFALSE; + if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){ + isGold = kTRUE; } - if (track1->fNDeadZone+track0->fNDeadZone>1.1) continue; - TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(itrack1); + AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1); + AliHelix &h1 = helixes[itrack0]; + AliHelix &h2 = helixes[itrack1]; // - //if (normdist0[itrack0]+normdist0[itrack1]<3) continue; - //if (normdist[itrack0]+normdist[itrack1]<4) continue; + // find linear distance + Double_t rmin =0; // // - AliHelix *h1 = &helixes[itrack0]; - AliHelix *h2 = &helixes[itrack1]; - Double_t rmin =0; - Double_t distance = TestV0(h1,h2,pvertex,rmin); // - if (distance>0.4) continue; - if (pvertex->GetRr()<0.3) continue; - if (pvertex->GetRr()>20.) continue; - pvertex->SetM(*track0); - pvertex->SetP(*track1); - pvertex->Update(primvertex); - if (pvertex->GetRr()<0.3) continue; - if (pvertex->GetRr()>20.) continue; - if (track1->fNDeadZone+track0->fNDeadZone>0.5 &&distance>0.12) continue; - + Double_t phase[2][2],radius[2]; + Int_t points = h1.GetRPHIintersections(h2, phase, radius); + if (points==0) continue; + Double_t delta[2]={1000,1000}; + rmin = radius[0]; + h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]); + if (points==2){ + if (radius[1]TMath::Min(maxr[itrack0],maxr[itrack1])) continue; + Float_t maxDist = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1)); + if (distance>maxDist) continue; + Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex); + if (pointAngleGetLabel())-TMath::Abs(track1->GetLabel())))<2 - ||(centries<5000&&(pvertex->GetPointAngle()>0.95))){ - //cstream<<"Iter0"<GetPointAngle()<0.85) continue; - // if (normdist[itrack0]+normdist[itrack1]<6&&pvertex->GetPointAngle()<0.99) continue; + // if (distance>maxDist) continue; + // if (pvertex->GetRr()GetRr()>kMaxR) continue; + AliITStrackMI * track0=btrack0; + AliITStrackMI * track1=btrack1; + // if (pvertex->GetRr()<3.5){ + if (radiusC<3.5){ + //use longest tracks inside the pipe + track0 = (AliITStrackMI*)trackarrayl.At(itrack0); + track1 = (AliITStrackMI*)trackarrayl.At(itrack1); + } // // + pvertex->SetM(*track0); + pvertex->SetP(*track1); + pvertex->Update(primvertex); + if (pvertex->GetRr()GetRr()>kMaxR) continue; + if (pvertex->GetPointAngle()GetDist2()>maxDist) continue; pvertex->SetLab(0,track0->GetLabel()); pvertex->SetLab(1,track1->GetLabel()); - pvertex->SetIndex(0,track0->GetESDtrack()->GetID()); - pvertex->SetIndex(1,track1->GetESDtrack()->GetID()); - // calculate chi2s - // - pvertex->SetChi2After(0); - pvertex->SetChi2Before(0); - pvertex->SetNBefore(0); - pvertex->SetNAfter(0); - for (Int_t i=0;i<6;i++){ - Float_t radius = fgLayers[i].GetR(); - if (pvertex->GetRr()>radius+0.5){ - pvertex->SetNBefore(pvertex->GetNBefore()+2.); - // - if (track0->fClIndex[i]<=0) { - pvertex->SetChi2Before(pvertex->GetChi2Before()+9); - }else{ - Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+ - track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]); - pvertex->SetChi2Before(pvertex->GetChi2Before()+chi2); - } + // + AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy; + AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy; - if (track1->fClIndex[i]<=0) { - pvertex->SetChi2Before(pvertex->GetChi2Before()+9); - - }else{ - Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+ - track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]); - // pvertex->fChi2Before+=chi2; - pvertex->SetChi2Before(pvertex->GetChi2Before()+chi2); + // + // + TObjArray * array0b = (TObjArray*)fBestHypothesys.At(itrack0); + if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->fP3)<1.1) + FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE); + TObjArray * array1b = (TObjArray*)fBestHypothesys.At(itrack1); + if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->fP3)<1.1) + FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE); + // + AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0); + AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1); + AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0); + AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1); + + Float_t minchi2before0=16; + Float_t minchi2before1=16; + Float_t minchi2after0 =16; + Float_t minchi2after1 =16; + Int_t maxLayer = GetNearestLayer(pvertex->GetXrp()); + + if (array0b) for (Int_t i=0;i<5;i++){ + // best track after vertex + AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i); + if (!btrack) continue; + if (btrack->fN>track0l->fN) track0l = btrack; + if (btrack->fXGetRr()-2) { + if (maxLayer>i+2 && btrack->fN==(6-i)&&i<2){ + Float_t sumchi2= 0; + Float_t sumn = 0; + for (Int_t ilayer=i;ilayerfClIndex[ilayer]){ + sumchi2+=25; + continue; + }else{ + sumchi2+=btrack->fDy[ilayer]*btrack->fDy[ilayer]/(btrack->fSigmaY[ilayer]*btrack->fSigmaY[ilayer]); + sumchi2+=btrack->fDz[ilayer]*btrack->fDz[ilayer]/(btrack->fSigmaZ[ilayer]*btrack->fSigmaZ[ilayer]); + } + } + sumchi2/=sumn; + if (sumchi2GetRr()SetNAfter(pvertex->GetNAfter()+2.); - // - if (track0->fClIndex[i]<=0) { - pvertex->SetChi2After(pvertex->GetChi2After()+9); - }else{ - Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+ - track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]); - pvertex->SetChi2After(pvertex->GetChi2After()+chi2); - } - - if (track1->fClIndex[i]<=0) { - pvertex->SetChi2After(pvertex->GetChi2After()+9.); - }else{ - Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+ - track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]); - pvertex->SetChi2After(pvertex->GetChi2After()+chi2); + track0b = btrack; + minchi2after0 = btrack->fNormChi2[i]; + break; + } + if (array1b) for (Int_t i=0;i<5;i++){ + // best track after vertex + AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i); + if (!btrack) continue; + if (btrack->fN>track1l->fN) track1l = btrack; + if (btrack->fXGetRr()-2){ + if (maxLayer>i+2&&btrack->fN==(6-i)&&(i<2)){ + Float_t sumchi2= 0; + Float_t sumn = 0; + for (Int_t ilayer=i;ilayerfClIndex[ilayer]){ + sumchi2+=30; + continue; + }else{ + sumchi2+=btrack->fDy[ilayer]*btrack->fDy[ilayer]/(btrack->fSigmaY[ilayer]*btrack->fSigmaY[ilayer]); + sumchi2+=btrack->fDz[ilayer]*btrack->fDz[ilayer]/(btrack->fSigmaZ[ilayer]*btrack->fSigmaZ[ilayer]); + } + } + sumchi2/=sumn; + if (sumchi2fNormChi2[i]; + break; } - if (pvertex->GetNBefore()>2){ - if (pvertex->GetChi2Before()/pvertex->GetNBefore()<4.) continue; //have clusters before vetex - } - Int_t ibest0=0,ibest1=0; - AliITStrackMI * ntrack0 = track0; - AliITStrackMI * ntrack1 = track1; - // - // - //PH Float_t oldistance = pvertex->GetDist2(); - Bool_t improve = FindBestPair(itrack0,itrack1,pvertex,ibest0,ibest1); // try to improve vertex - if (pvertex->GetPointAngle()<0.5) continue; - distance = pvertex->GetDist2(); - if (improve){ - ntrack0 = (AliITStrackMI*)array0->At(ibest0); - ntrack1 = (AliITStrackMI*)array1->At(ibest1); - } - Bool_t accept0 = kFALSE; // accept ==> because of pointing angle - if (pvertex->GetPointAngle()>0.999){ - if (pvertex->GetRr()<3.5 && (ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone)<11.5) continue; - if (pvertex->GetRr()>3.5&& pvertex->GetDistNorm()<12) accept0 = kTRUE; - if (pvertex->GetRr()>1 && pvertex->GetDist2()<0.1 && pvertex->GetDistNorm()<12) accept0 = kTRUE; - if (pvertex->GetPointAngle()>0.9995&&pvertex->GetRr()>5.) accept0 = kTRUE; - } - Bool_t reject1= kFALSE; // reject ==> bad kinematic - // - reject1 |= TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>1.02 || - TMath::Abs(ntrack0->fN-ntrack1->fN)>1.02; // cut1 - reject1 |= ntrack0->fNUsed+ntrack1->fNUsed>1.01; // cut2 - reject1 |= pvertex->GetDistNorm()>12; // cut3 - reject1 |= pvertex->GetDist2()>0.1 && improve; // cut4 - reject1 |= (TMath::Abs(ntrack0->fD[0])+TMath::Abs(ntrack1->fD[0]))/pvertex->GetDist2()<5; //cut5 - reject1 |= TMath::Abs(ntrack0->fD[0]/pvertex->GetDist2())<2 || TMath::Abs(ntrack1->fD[0]/pvertex->GetDist2())<2; //cut 6 // - // small radii cuts - Bool_t reject2 = kFALSE; - if (pvertex->GetRr()<3.6){ - reject2 |= (TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>0.01); // cut7 - reject2 |= ntrack0->fNUsed+ntrack1->fNUsed>0.01; // cut8 - reject2 |= ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone<11.5; // cut9 - reject2 |= (ntrack0->fN+ntrack1->fN<11.5)&&pvertex->GetRr()<2; // cut10 - reject2 |= pvertex->GetDist2()>0.1; // cut11 - } - //PH AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy; - //PH AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy; + // position resolution - used for DCA cut + Float_t sigmad = track0b->fC00+track0b->fC11+track1b->fC00+track1b->fC11+ + (track0b->fX-pvertex->GetRr())*(track0b->fX-pvertex->GetRr())*(track0b->fC22+track0b->fC33)+ + (track1b->fX-pvertex->GetRr())*(track1b->fX-pvertex->GetRr())*(track1b->fC22+track1b->fC33); + sigmad =TMath::Sqrt(sigmad)+0.04; + if (pvertex->GetRr()>50){ + Double_t cov0[15],cov1[15]; + track0b->fESDtrack->GetInnerExternalCovariance(cov0); + track1b->fESDtrack->GetInnerExternalCovariance(cov1); + sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+ + (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+ + (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]); + sigmad =TMath::Sqrt(sigmad)+0.05; + } + // + AliESDV0MI vertex2; + vertex2.SetM(*track0b); + vertex2.SetP(*track1b); + vertex2.Update(primvertex); + if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetPointAngle()>=pvertex->GetPointAngle())){ + pvertex->SetM(*track0b); + pvertex->SetP(*track1b); + pvertex->Update(primvertex); + } + pvertex->SetDistSigma(sigmad); + pvertex->SetDistNorm(pvertex->GetDist2()/sigmad); // + // define likelihhod and causalities // - // - if ( TMath::Abs((TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel())))<2 - ||(centries<500000)){ - /* - cstream<<"It1"<<"Tr0.="<GetAnglep()[2]>0.2){ + pb0 = TMath::Exp(-TMath::Min(normdist[itrack0],Float_t(16.))/12.); + pb1 = TMath::Exp(-TMath::Min(normdist[itrack1],Float_t(16.))/12.); + } + pvertex->SetChi2Before(normdist[itrack0]); + pvertex->SetChi2After(normdist[itrack1]); + pvertex->SetNAfter(0); + pvertex->SetNBefore(0); + }else{ + pvertex->SetChi2Before(minchi2before0); + pvertex->SetChi2After(minchi2before1); + if (pvertex->GetAnglep()[2]>0.2){ + pb0 = TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.); + pb1 = TMath::Exp(-TMath::Min(minchi2before1,float_t(16))/12.); + } + pvertex->SetNAfter(maxLayer); + pvertex->SetNBefore(maxLayer); } - - if (!accept0 && (reject1 || reject2)) continue; - -// if (distance>0.5) continue; -// distance = pvertex->GetDist2(); -// if (pvertex->GetRr()>25 || pvertex->GetRr()<0.2) continue; -// if (pvertex->GetRr()/pvertex->fDistSigma<1) continue; -// if (pvertex->GetDistNorm()>10) continue; -// if (pvertex->GetPointAngle()<0.85) continue; -// if ((normdist[itrack0]<3||normdist[itrack1]<3)){ -// if (pvertex->GetPointAngle()<0.99||pvertex->GetDist2()>0.15) continue; -// } -// if (distance>0.05*(0.8+0.2*(0.5+pvertex->GetRr()))) continue; -// if (pvertex->GetRr()<0.3) continue; -// if (pvertex->GetRr()>27.) continue; - - - - // - if (distance<0.3 &&pvertex->GetPointAngle()>0.998){ - track0->fGoldV0 = kTRUE; - track1->fGoldV0 = kTRUE; + if (pvertex->GetRr()<90){ + pa0 *= TMath::Min(track0->fESDtrack->GetTPCdensity(0,60),Float_t(1.)); + pa1 *= TMath::Min(track1->fESDtrack->GetTPCdensity(0,60),Float_t(1.)); } - vertexes++; - vertexall++; - if (vertexall>=100000) break; - pvertex = &vertexarray[vertexall]; - } - } - // printf("\n\n\nMultifound\t%d\n\n\n",multifound); - // - // sort vertices according quality - Float_t quality[40000]; - Int_t indexes[40000]; - Int_t trackvertices[40000]; - for (Int_t i=0;iAddV0MI(&vertexarray[indexes[i]]); - // - if (trackvertices[index1]+trackvertices[index0]>5) continue; - if (trackvertices[index0]>2) continue; - if (trackvertices[index1]>2) continue; - - if (trackvertices[index1]+trackvertices[index0]>0) { - // if (pvertex->GetPointAngle()<0.995) continue; - } - trackvertices[index0]++; - trackvertices[index1]++; - - AliESDtrack * ptrack0 = event->GetTrack(vertexarray[indexes[i]].GetIndex(0)); - AliESDtrack * ptrack1 = event->GetTrack(vertexarray[indexes[i]].GetIndex(1)); - if (!ptrack0 || !ptrack1){ - printf("BBBBBBBUUUUUUUUUUGGGGGGGGGG\n"); - } - Int_t v0index0[3]={ptrack0->GetV0Index(0),ptrack0->GetV0Index(1),ptrack0->GetV0Index(2)}; - Int_t v0index1[3]={ptrack1->GetV0Index(0),ptrack1->GetV0Index(1),ptrack1->GetV0Index(2)}; - for (Int_t i=0;i<3;i++){ - if (v0index0[i]<0) { - v0index0[i]=v0index; - ptrack0->SetV0Indexes(v0index0); - break; - } - } - for (Int_t i=0;i<3;i++){ - if (v0index1[i]<0) { - v0index1[i]=v0index; - ptrack1->SetV0Indexes(v0index1); - break; + if (pvertex->GetRr()<20){ + pa0 *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2; + pa1 *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2; } - } - } - - delete [] helixes; - delete [] dist; - delete [] normdist0; - delete [] normdist1; - delete [] normdist; - delete [] norm; - - delete[] vertexarray; -} - - - -Bool_t AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1, AliESDV0MI *vertex, Int_t &i0, Int_t &i1) -{ - // - // try to find best pair from the tree of track hyp. - // - TObjArray * array0 = (TObjArray*)fTrackHypothesys.At(esdtrack0); - Int_t entries0 = array0->GetEntriesFast(); - TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(esdtrack1); - Int_t entries1 = array1->GetEntriesFast(); - // AliITStrackMI *orig0 = (AliITStrackMI*)fOriginal.At(esdtrack0); - //AliITStrackMI *orig1 = (AliITStrackMI*)fOriginal.At(esdtrack1); - Double_t criticalradius = vertex->GetRr(); - AliITStrackMI * track0= 0; - AliITStrackMI * track1= 0; - i0 = -1; - i1 = -1; - // - // - Float_t rfirst0[2000]; //radius position of the first cluster - track0 - Float_t rfirst1[2000]; // - track1 - Float_t maxlocalx0=0; //local x for first track - Float_t maxlocalx1=0; //local x for second track - Float_t cs0=1, sn0=0; //rotations - Float_t cs1=1, sn1=0; //rotations - - // - for (Int_t itrack0=0;itrack0At(itrack0); - if (!htrack0) continue; - if (htrack0->fConstrain) continue; - if (i0<0){ - i0 = itrack0; - track0 = htrack0; - } - Double_t cs = TMath::Cos(htrack0->fAlpha); - Double_t sn = TMath::Sin(htrack0->fAlpha); - Double_t x = htrack0->fX*cs - htrack0->fP0*sn; - Double_t y = htrack0->fX*sn + htrack0->fP0*cs; - Double_t radius = TMath::Sqrt(x*x+y*y); - if (criticalradius<3 && radius>6&&htrack0->fNDeadZone<0.2) continue; // all cluster required - if (criticalradius>10 && radius<6) continue; // causality - Double_t localx = TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn); - if (localx>maxlocalx0) { - maxlocalx0=localx; - cs0 = cs; sn0=sn; - } - rfirst0[itrack0] = radius; - } - for (Int_t itrack1=0;itrack1At(itrack1); - if (!htrack1) continue; - if (htrack1->fConstrain) continue; - if (i1<0){ - i1 = itrack1; - track1 = htrack1; - } - Double_t cs = TMath::Cos(htrack1->fAlpha); - Double_t sn = TMath::Sin(htrack1->fAlpha); - // - // - Double_t x = htrack1->fX*cs - htrack1->fP0*sn; - Double_t y = htrack1->fX*sn + htrack1->fP0*cs; - Double_t radius = TMath::Sqrt(x*x+y*y); - if (criticalradius<3 && radius>6&&htrack1->fNDeadZone<0.2) continue; //all clusters required - if (criticalradius>10 && radius<6) continue; //causality - Double_t localx = TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn); - if (localx>maxlocalx1) { - maxlocalx1=localx; - cs1 = cs; sn1=sn; - } - rfirst1[itrack1] = radius; - } - // - // - // - const Float_t radiuses[4]={4,6.5,15.03,24.}; - // - // - // find the best tracks after decay point - Float_t bestquality =100000; - Float_t bestradius=0; - AliESDV0MI v0; - Double_t vpos[3]; - Float_t v[3]={GetX(),GetY(),GetZ()}; - // - // - for (Int_t itrack0=0;itrack0At(itrack0); - if (!htrack0) continue; - // - for (Int_t itrack1=0;itrack1At(itrack1); - if (!htrack1) continue; - if (TMath::Abs(rfirst0[itrack0]-rfirst1[itrack1])>1.) continue; - if (htrack0->fClIndex[6-htrack0->fN]==htrack1->fClIndex[6-htrack1->fN]) continue; //sharing of last cluster not allowe // - if (htrack0->fNUsed+htrack1->fNUsed>0) continue; //sharing of clusters not allowed + pvertex->SetCausality(pb0,pb1,pa0,pa1); // - // - v0.SetM(*htrack0); - v0.SetP(*htrack1); - // if (v0.Update(v)==0) continue; - v0.Update(v); - if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4]; + p12 += pvertex->GetParamM()->GetParameter()[4]*pvertex->GetParamM()->GetParameter()[4]; + p12 = TMath::Sqrt(p12); // "mean" momenta + Float_t sigmap0 = 0.0001+0.001/(0.1+pvertex->GetRr()); + Float_t sigmap = 0.5*sigmap0*(0.6+0.4*p12); // "resolution: of point angle - as a function of radius and momenta + + Float_t CausalityA = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]); + Float_t CausalityB = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Float_t(0.7))* + TMath::Min(pvertex->GetCausalityP()[3],Float_t(0.7))); // + Float_t Likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5); + + Float_t Likelihood1 = TMath::Exp(-(1.0001-pvertex->GetPointAngle())/sigmap)+ + 0.4*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/(4.*sigmap))+ + 0.4*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/(8.*sigmap))+ + 0.1*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/0.01); // - if (v0.GetDist2()>0.3) continue; - if (v0.GetRr()fN+htrack0->fNDeadZone-htrack1->fN-htrack1->fNDeadZone)>0.5)) continue; + if (CausalityAfN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<11.7) continue; - //if (v0.GetRr()<6. && (htrack0->fN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<9.7) continue; - if (v0.GetRr()<3. && (htrack0->fN+htrack1->fN)<11.7) continue; - if (v0.GetRr()<6. && (htrack0->fN+htrack1->fN)<9.7) continue; - Double_t localx0=v0.GetXr(0)*cs0+v0.GetXr(1)*sn0; - Double_t localx1=v0.GetXr(0)*cs1+v0.GetXr(1)*sn1; - Double_t maxlocalx = TMath::Max(localx0,localx1); - if (maxlocalx<3.4 && (htrack0->fN+htrack1->fN)<11.7) continue; - if (maxlocalx<6.1 && (htrack0->fN+htrack1->fN)<9.7) continue; // - Float_t fnormdist = v0.GetDist2()/0.05; - fnormdist +=(htrack0->fNormChi2[6-htrack0->fN]+htrack1->fNormChi2[6-htrack1->fN]); - fnormdist +=3*(htrack0->fNUsed+htrack1->fNUsed); - if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) GetLabel())-TMath::Abs(track1->GetLabel()))==1; + cstream<<"It0"<< + "Tr0.="<GetDist2()){ - i0=itrack0; - i1=itrack1; - track0 =htrack0; - track1 =htrack1; - bestquality = quality; - bestradius = v0.GetRr(); - vpos[0] = v0.GetXr(0); - vpos[1] = v0.GetXr(1); - vpos[2] = v0.GetXr(2); + pvertex->SetStatus(0); + // if (rejectBase) { + // pvertex->SetStatus(-100); + //} + if (pvertex->GetPointAngle()>kMinPointAngle2) { + if (v0OK){ + AliV0vertex vertexjuri(*track0,*track1); + vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID()); + pvertex->SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID()); + event->AddV0(&vertexjuri); + pvertex->SetStatus(100); + } + event->AddV0MI(pvertex); } } } // // - if (!track0||!track1) return kFALSE; - if (track0->fNUsed+track1->fNUsed>0) return kFALSE; // sharing of clusters not allowed - // - //propagate to vertex - - Double_t alpha = TMath::ATan2(vpos[1],vpos[0]); - // - AliITStrackMI track0p = *track0; - AliITStrackMI track1p = *track1; - // - // RefitAt(bestradius+0.5,&track0p,track0); - //RefitAt(bestradius+0.5,&track1p,track1); - - if ( track0p.Propagate(alpha,bestradius+0.2)) { - track0= &track0p; - } - if (track1p.Propagate(alpha,bestradius+0.2)){ - track1 = &track1p; - } - // - v0.SetM(*track0); - v0.SetP(*track1); - v0.Update(v); - if (v0.GetDist2()GetDist2() && v0.GetRr()<20){ - vertex->SetM(*track0); - vertex->SetP(*track1); - vertex->Update(v); - return kTRUE; - } - return kFALSE; + // delete temporary arrays + // + delete[] minPointAngle; + delete[] maxr; + delete[] minr; + delete[] norm; + delete[] normdist; + delete[] normdist1; + delete[] normdist0; + delete[] dist; + delete[] itsmap; + delete[] helixes; + delete pvertex; } -Double_t AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliESDV0MI *vertex, Double_t &rmin) -{ - // - // test the helixes for the distnce calculate vertex characteristic - // - rmin =0; - Float_t distance1,distance2; - AliHelix & dhelix1 = *helix1; - Double_t pp[3],xx[3]; - dhelix1.GetMomentum(0,pp,0); - dhelix1.Evaluate(0,xx); - AliHelix &mhelix = *helix2; - // - //find intersection linear - // - Double_t phase[2][2],radius[2]; - Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius); - Double_t delta1=10000,delta2=10000; - - if (points>0){ - rmin = radius[0]; - dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); - dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); - dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); - } - if (points==2){ - if (radius[1]SetDist1(TMath::Sqrt(distance1)); - - // - //find intersection parabolic - // - points = dhelix1.GetRPHIintersections(mhelix, phase, radius); - delta1=10000,delta2=10000; - - if (points>0){ - dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); - } - if (points==2){ - dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); - } - - distance2 = TMath::Min(delta1,delta2); - vertex->SetDist2(TMath::Sqrt(distance2)); - if (distance2<100){ - if (delta1GetXrp()); - dhelix1.GetMomentum(phase[0][0],vertex->GetPPp()); - mhelix.GetMomentum(phase[0][1],vertex->GetPMp()); - dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],vertex->GetAnglep()); - vertex->SetRr(TMath::Sqrt(radius[0])); - } - else{ - dhelix1.Evaluate(phase[1][0],vertex->GetXrp()); - dhelix1.GetMomentum(phase[1][0],vertex->GetPPp()); - mhelix.GetMomentum(phase[1][1],vertex->GetPMp()); - dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],vertex->GetAnglep()); - vertex->SetRr(TMath::Sqrt(radius[1])); - } - } - // - // - return vertex->GetDist2(); -} + diff --git a/ITS/AliITStrackerMI.h b/ITS/AliITStrackerMI.h index 79369a33faf..b6f3d9ad790 100644 --- a/ITS/AliITStrackerMI.h +++ b/ITS/AliITStrackerMI.h @@ -23,6 +23,8 @@ class AliHelix; class AliV0vertex; class AliESDV0MI; +class TTreeSRedirector; + @@ -110,7 +112,7 @@ public: Int_t GetSkip() const {return fSkip;} void SetSkip(Int_t skip){fSkip=skip;} void IncAccepted(){fAccepted++;} - Int_t GetAccepted() const {return fAccepted;} + Int_t GetAccepted() const {return fAccepted;} protected: AliITSlayer(const AliITSlayer& /*layer*/){;} Double_t fR; // mean radius of this layer @@ -172,13 +174,13 @@ public: AliITStrackerMI::AliITSdetector & GetDetector(Int_t layer, Int_t n) const {return GetLayer(layer).GetDetector(n); } protected: - void FindV0(AliESD *event); //try to find V0 - Double_t TestV0(AliHelix *h1, AliHelix *h2, AliESDV0MI *vertex, Double_t &rmin); //try to find V0 - return DCA - Bool_t FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliESDV0MI *vertex, Int_t &ibest0, Int_t &ibest1); // try to find best pair from the tree of track hyp. + Int_t GetNearestLayer(const Double_t *xr) const; //get nearest upper layer close to the point xr + void FindV02(AliESD *event); //try to find V0 + void UpdateTPCV0(AliESD *event); //try to update, or reject TPC V0s void CookLabel(AliKalmanTrack *t,Float_t wrong) const; void CookLabel(AliITStrackMI *t,Float_t wrong) const; Double_t GetEffectiveThickness(Double_t y, Double_t z) const; - void FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex); + void FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain); void ResetBestTrack() { fBestTrack.~AliITStrackMI(); new(&fBestTrack) AliITStrackMI(fTrackToFollow); @@ -218,6 +220,7 @@ protected: AliITStrackMI fBestTrack; // "best" track AliITStrackMI fTrackToFollow; // followed track TObjArray fTrackHypothesys; // ! array with track hypothesys- ARRAY is the owner of tracks- MI + TObjArray fBestHypothesys; // ! array with track hypothesys- ARRAY is the owner of tracks- MI TObjArray fOriginal; // ! array with seeds from the TPC Int_t fBestTrackIndex[100000]; // ! index of the best track Int_t fCurrentEsdTrack; // ! current esd track - MI @@ -227,6 +230,8 @@ protected: Int_t fLayersNotToSkip[kMaxLayer]; // layer masks Int_t fLastLayerToTrackTo; // the innermost layer to track to Float_t * fCoeficients; //! working array with errors and mean cluser shape + AliESD * fEsd; //! pointer to the ESD event + TTreeSRedirector *fDebugStreamer; //!debug streamer private: AliITStrackerMI(const AliITStrackerMI * /*tracker*/){;} ClassDef(AliITStrackerMI,2) //ITS tracker MI diff --git a/ITS/AliITStrackerV2.cxx b/ITS/AliITStrackerV2.cxx index efced05e745..be13f576463 100644 --- a/ITS/AliITStrackerV2.cxx +++ b/ITS/AliITStrackerV2.cxx @@ -670,6 +670,10 @@ Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) { //-------------------------------------------------------------------- Float_t circ=TMath::TwoPi()*fR; Int_t sec=Int_t(kNsector*c->GetPhiR()/circ); + if (sec>=kNsector) { + ::Error("InsertCluster","Wrong sector !\n"); + return 1; + } Int_t &n=fN[sec]; if (n>=kMaxClusterPerSector) { ::Error("InsertCluster","Too many clusters !\n"); diff --git a/STEER/AliESDComparisonMI.C b/STEER/AliESDComparisonMI.C index 74dc37d633d..ffaa2833878 100644 --- a/STEER/AliESDComparisonMI.C +++ b/STEER/AliESDComparisonMI.C @@ -189,6 +189,49 @@ TProfile prof("prof","prof",10,0.5,5); +TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01"); +TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5"); +TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5"); +TCut crec("crec","fReconstructed==1"); +TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25"); +TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1"); + +TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1"); +TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50"); +TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)"); +TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6."); + + + +void MakeAliases(AliESDComparisonDraw&comp) +{ + // + // aliases definition + // + comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)"); + comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy"); + comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy"); + comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()"); + comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)"); + comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); + comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN"); + + comp.fTree->SetAlias("trddedx","(RC.fESDTrack.fTRDsignals[0]+RC.fESDTrack.fTRDsignals[1]+RC.fESDTrack.fTRDsignals[2]+RC.fESDTrack.fTRDsignals[3]+RC.fESDTrack.fTRDsignals[4]+RC.fESDTrack.fTRDsignals[5])/6."); + + comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)"); + comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)"); + + comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]"); + comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum"); + comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum"); + comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum"); + comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum"); + comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum"); + comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)"); +} + + + void AliESDRecInfo::UpdatePoints(AliESDtrack*track) { @@ -410,7 +453,7 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconst } - if (fStatus[1]>0 &&info->fNTPCRef>0){ + if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){ //TPC fESDTrack.GetInnerXYZ(fTPCinR1); fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]); @@ -930,6 +973,7 @@ Int_t ESDCmpTr::Exec() fNextTreeGenEntryToRead = 0; fNextKinkToRead = 0; + fNextV0ToRead =0; cerr<<"fFirstEventNr, fNEvents: "<fESDTrack =*track; if (track->GetITStrack()) fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack()); + else{ + fRecInfo->fITStrack = *track; + } if (track->GetTRDtrack()){ fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack()); } diff --git a/STEER/AliESDV0MI.cxx b/STEER/AliESDV0MI.cxx index 0fad89192bd..13d839543cb 100644 --- a/STEER/AliESDV0MI.cxx +++ b/STEER/AliESDV0MI.cxx @@ -51,23 +51,50 @@ AliESDV0MI::AliESDV0MI() : // //Dafault constructor // + for (Int_t i=0;i<4;i++){fCausality[i]=0;} +} + + +void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1) +{ + // + // set probabilities + // + fCausality[0] = pb0; // probability - track 0 exist before vertex + fCausality[1] = pb1; // probability - track 1 exist before vertex + fCausality[2] = pa0; // probability - track 0 exist close after vertex + fCausality[3] = pa1; // probability - track 1 exist close after vertex } void AliESDV0MI::SetP(const AliExternalTrackParam & paramp) { // - // set mother + // set track + // fParamP = paramp; } void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){ // - //set daughter + //set track - // fParamM = paramm; - } +void AliESDV0MI::SetRp(const Double_t *rp){ + // + // set pid + + // + for (Int_t i=0;i<5;i++) fRP[i]=rp[i]; +} + +void AliESDV0MI::SetRm(const Double_t *rm){ + // + // set pid - + // + for (Int_t i=0;i<5;i++) fRM[i]=rm[i]; +} + + void AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5]) { // @@ -136,7 +163,8 @@ void AliESDV0MI::Update(Float_t vertex[3]) // // updates Kink Info // - Float_t distance1,distance2; + // Float_t distance1,distance2; + Float_t distance2; // AliHelix phelix(fParamP); AliHelix mhelix(fParamM); @@ -146,7 +174,7 @@ void AliESDV0MI::Update(Float_t vertex[3]) Double_t phase[2][2],radius[2]; Int_t points = phelix.GetRPHIintersections(mhelix, phase, radius,200); Double_t delta1=10000,delta2=10000; - + /* if (points<=0) return; if (points>0){ phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1); @@ -159,6 +187,7 @@ void AliESDV0MI::Update(Float_t vertex[3]) phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2); } distance1 = TMath::Min(delta1,delta2); + */ // //find intersection parabolic // @@ -215,13 +244,13 @@ void AliESDV0MI::Update(Float_t vertex[3]) fDist2 = TMath::Sqrt(distance2); // // - Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]}; - Float_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]}; - Float_t vnorm2 = v[0]*v[0]+v[1]*v[1]; - Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2); + Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]}; + Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]}; + Double_t vnorm2 = v[0]*v[0]+v[1]*v[1]; + Double_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2); vnorm2 = TMath::Sqrt(vnorm2); - Float_t pnorm2 = p[0]*p[0]+p[1]*p[1]; - Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2); + Double_t pnorm2 = p[0]*p[0]+p[1]*p[1]; + Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2); pnorm2 = TMath::Sqrt(pnorm2); fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2); fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3); diff --git a/STEER/AliESDV0MI.h b/STEER/AliESDV0MI.h index a0c9878c458..4828d08af18 100644 --- a/STEER/AliESDV0MI.h +++ b/STEER/AliESDV0MI.h @@ -22,9 +22,15 @@ public: // friend class AliITStrackerMI; AliESDV0MI(); //constructor // + const AliExternalTrackParam *GetParamP() const {return &fParamP;} + const AliExternalTrackParam *GetParamM() const {return &fParamM;} void SetP(const AliExternalTrackParam & paramp); void SetM(const AliExternalTrackParam & paramd); + void SetRp(const Double_t *rp); + void SetRm(const Double_t *rm); void UpdatePID(Double_t pidp[5], Double_t pidm[5]); + void SetStatus(Int_t status){fStatus=status;} + Int_t GetStatus() const {return fStatus;} Float_t GetEffMass(UInt_t p1, UInt_t p2); Float_t GetProb(UInt_t p1, UInt_t p2); void Update(Float_t vertex[3]); //update @@ -58,7 +64,8 @@ public: Float_t GetNBefore() const {return fNBefore;} void SetNBefore(Float_t nb) {fNBefore=nb;} void SetLab(Int_t i, Int_t lab) {fLab[i]=lab;} - + void SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1); + const Float_t * GetCausalityP() const {return fCausality;} private: AliExternalTrackParam fParamP; AliExternalTrackParam fParamM; @@ -78,12 +85,13 @@ private: Double_t fXr[3]; //rec. position according helix Double_t fAngle[3]; //three angles Double_t fRr; //rec position of the vertex - Int_t fStatus; //status + Int_t fStatus; //status - 1 - TPC V0 2- ITS V0 4- accepted - 0 -rejected Int_t fRow0; // critical layer Int_t fOrder[3]; //order of the vertex // quality information Double_t fDistNorm; //normalized DCA Double_t fDistSigma; //sigma of distance + Float_t fCausality[4]; // causality information - see comments in SetCausality Float_t fChi2Before; //chi2 of the tracks before V0 Float_t fNBefore; // number of possible points before V0 Float_t fChi2After; // chi2 of the tracks after V0 @@ -92,9 +100,7 @@ private: Float_t fPointAngleTh; //point angle theta Float_t fPointAngle; //point angle full - - - ClassDef(AliESDV0MI,1) // ESD V0 vertex + ClassDef(AliESDV0MI,2) // ESD V0 vertex }; diff --git a/STEER/AliESDtrack.cxx b/STEER/AliESDtrack.cxx index 86262cbf460..de33b519e3a 100644 --- a/STEER/AliESDtrack.cxx +++ b/STEER/AliESDtrack.cxx @@ -833,10 +833,42 @@ Int_t AliESDtrack::GetTPCclusters(Int_t *idx) const { return fTPCncls; } +Float_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{ + // + // GetDensity of the clusters on given region between row0 and row1 + // Dead zone effect takin into acoount + // + Int_t good = 0; + Int_t found = 0; + // + for (Int_t i=row0;i<=row1;i++){ + Int_t index = fTPCindex[i]; + if (index!=-1) good++; // track outside of dead zone + if (index>0) found++; + } + Float_t density=0.5; + if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good); + return density; +} //_______________________________________________________________________ void AliESDtrack::SetTPCpid(const Double_t *p) { // Sets values for the probability of each particle type (in TPC) - for (Int_t i=0; i0){ +// fTPCr[i]=p[i]/sump; +// } +// else{ +// fTPCr[i]=p[i]; +// } +// } + for (Int_t i=0; ifMCm = (*info); v0info->fMCd = (*dinfo); + v0info->fMotherP = (*motherparticle); v0info->Update(fVPrim); br->SetAddress(&v0info); fTreeV0->Fill(); diff --git a/STEER/AliGenInfo.h b/STEER/AliGenInfo.h index 1341b6057f8..bc4acb1e784 100644 --- a/STEER/AliGenInfo.h +++ b/STEER/AliGenInfo.h @@ -80,6 +80,7 @@ class AliGenV0Info: public TObject { public: AliMCInfo fMCd; //info about daughter particle - second particle for V0 AliMCInfo fMCm; //info about mother particle - first particle for V0 + TParticle fMotherP; //particle info about mother particle void Update(Float_t vertex[3]); // put some derived info to special field Double_t fMCDist1; //info about closest distance according closest MC - linear DCA Double_t fMCDist2; //info about closest distance parabolic DCA diff --git a/STEER/AliHelix.cxx b/STEER/AliHelix.cxx index 670d0c74666..1810a0f68f3 100644 --- a/STEER/AliHelix.cxx +++ b/STEER/AliHelix.cxx @@ -156,7 +156,7 @@ AliHelix::AliHelix(Double_t x[3], Double_t p[3], Double_t charge, Double_t conve } -void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion) +void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, Double_t *xr) { // return momentum at given phase Double_t x[3],g[3],gg[3]; @@ -166,6 +166,9 @@ void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion) p[0] = fHelix[8]*g[0]/(mt*conversion); p[1] = fHelix[8]*g[1]/(mt*conversion); p[2] = fHelix[8]*g[2]/(mt*conversion); + if (xr){ + xr[0] = x[0]; xr[1] = x[1]; xr[2] = x[2]; + } } void AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3]) @@ -267,6 +270,26 @@ Int_t AliHelix::GetClosestPhases(AliHelix &h, Double_t phase[2][2]) return 1; } +Double_t AliHelix::GetPointAngle(AliHelix &h, Double_t phase[2], const Float_t * vertex) +{ + // + // get point angle bettwen two helixes + // + Double_t r0[3],p0[4]; + Double_t r1[3],p1[4]; + GetMomentum(phase[0],p0,1,r0); + h.GetMomentum(phase[1],p1,1,r1); + // + Double_t r[3] = {(r0[0]+r1[0])*0.5-vertex[0],(r0[1]+r1[1])*0.5-vertex[1],(r0[2]+r1[2])*0.5-vertex[2]}; + //intersection point - relative to the prim vertex + Double_t p[3] = { p0[0]+p1[0], p0[1]+p1[1],p0[2]+p1[2]}; + // derivation vector + Double_t normr = TMath::Sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]); + Double_t normp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]); + Double_t pointAngle = (r[0]*p[0]+r[1]*p[1]+r[2]*p[2])/(normr*normp); + return pointAngle; +} + Double_t AliHelix::GetPhase(Double_t x, Double_t y ) { diff --git a/STEER/AliHelix.h b/STEER/AliHelix.h index 804caaeb490..2a21b04c13a 100644 --- a/STEER/AliHelix.h +++ b/STEER/AliHelix.h @@ -31,7 +31,7 @@ public: void Evaluate(Double_t t, Double_t r[3], //radius vector Double_t g[3], //first defivatives Double_t gg[3]); //second derivatives - void GetMomentum(Double_t phase, Double_t p[4], Double_t conversion=0.); // return momentum + void GetMomentum(Double_t phase, Double_t p[4], Double_t conversion=0., Double_t *xr=0); // return momentum void GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3]); inline Double_t GetHelixR(Double_t phase=0); inline Double_t GetHelixZ(Double_t phase=0); @@ -41,6 +41,7 @@ public: Int_t GetPhase(Double_t r0, Double_t t[2]); //return phase for the nearest point Int_t GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Double_t ri[2], Double_t cut=3.); Int_t GetClosestPhases(AliHelix &h, Double_t phase[2][2]); + Double_t GetPointAngle(AliHelix &h, Double_t phase[2],const Float_t *vertex); Int_t LinearDCA(AliHelix &h, Double_t &t1, Double_t &t2, Double_t &R, Double_t &dist); // diff --git a/STEER/AliMagFMaps.cxx b/STEER/AliMagFMaps.cxx index 1c78a65184d..c4b2d6dd375 100644 --- a/STEER/AliMagFMaps.cxx +++ b/STEER/AliMagFMaps.cxx @@ -207,7 +207,7 @@ void AliMagFMaps::Field(Float_t *x, Float_t *b) const // // Constant L3 field , if this option was selected // - b[2] = - fSolenoid; + b[2] = (- fSolenoid)*fFactor; return; } } else if (fFieldMap[1]->Inside(xm[0], xm[1], xm[2])) { diff --git a/STEER/TTreeStream.cxx b/STEER/TTreeStream.cxx index a6cdeba07ec..f2511810265 100644 --- a/STEER/TTreeStream.cxx +++ b/STEER/TTreeStream.cxx @@ -232,7 +232,7 @@ TTreeDataElement:: TTreeDataElement(Char_t type) : TTreeDataElement:: TTreeDataElement(TDataType* type) : fName(), - fType(' '), + fType(0), fDType(type), fClass(0), fPointer(0) @@ -244,7 +244,7 @@ TTreeDataElement:: TTreeDataElement(TDataType* type) : TTreeDataElement:: TTreeDataElement(TClass* cl) : fName(), - fType(' '), + fType(0), fDType(0), fClass(cl), fPointer(0) diff --git a/TPC/AliTPCseed.cxx b/TPC/AliTPCseed.cxx new file mode 100644 index 00000000000..3d802591f63 --- /dev/null +++ b/TPC/AliTPCseed.cxx @@ -0,0 +1,922 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + + + + +//----------------------------------------------------------------- +// Implementation of the TPC seed class +// This class is used by the AliTPCtrackerMI class +// Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch +//----------------------------------------------------------------- +#include "TClonesArray.h" +#include "AliTPCseed.h" + +ClassImp(AliTPCseed) + + + +AliTPCseed::AliTPCseed():AliTPCtrack(){ + // + fRow=0; + fRemoval =0; + for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); + for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; + for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0; + for (Int_t i=0;i<5;i++) fTPCr[i]=0.2; + fPoints = 0; + fEPoints = 0; + fNFoundable =0; + fNShared =0; + fRemoval = 0; + fSort =0; + fFirstPoint =0; + fNoCluster =0; + fBSigned = kFALSE; + fSeed1 =-1; + fSeed2 =-1; + fCurrentCluster =0; + fCurrentSigmaY2=0; + fCurrentSigmaZ2=0; + fCircular = 0; // not curling track +} +AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){ + //--------------------- + // dummy copy constructor + //------------------------- + for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i]; + for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i]; + + fPoints = 0; + fEPoints = 0; + fCircular =0; +} +AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){ + // + //copy constructor + fPoints = 0; + fEPoints = 0; + fNShared =0; + // fTrackPoints =0; + fRemoval =0; + fSort =0; + for (Int_t i=0;i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i); + + for (Int_t i=0;i<160;i++) { + fClusterPointer[i] = 0; + Int_t index = t.GetClusterIndex(i); + if (index>=-1){ + SetClusterIndex2(i,index); + } + else{ + SetClusterIndex2(i,-3); + } + } + fFirstPoint =0; + fNoCluster =0; + fBSigned = kFALSE; + fSeed1 =-1; + fSeed2 =-1; + fCurrentCluster =0; + fCurrentSigmaY2=0; + fCurrentSigmaZ2=0; + fCircular =0; +} + +AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15], + Double_t xr, Double_t alpha): + AliTPCtrack(index, xx, cc, xr, alpha) { + // + // + //constructor + fRow =0; + for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); + for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; + for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0; + for (Int_t i=0;i<5;i++) fTPCr[i]=0.2; + + fPoints = 0; + fEPoints = 0; + fNFoundable =0; + fNShared = 0; + // fTrackPoints =0; + fRemoval =0; + fSort =0; + fFirstPoint =0; + // fHelixIn = new TClonesArray("AliHelix",0); + //fHelixOut = new TClonesArray("AliHelix",0); + fNoCluster =0; + fBSigned = kFALSE; + fSeed1 =-1; + fSeed2 =-1; + fCurrentCluster =0; + fCurrentSigmaY2=0; + fCurrentSigmaZ2=0; +} + +AliTPCseed::~AliTPCseed(){ + // + // destructor + if (fPoints) delete fPoints; + fPoints =0; + if (fEPoints) delete fEPoints; + fEPoints = 0; + fNoCluster =0; +} + +AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i) +{ + // + // + return &fTrackPoints[i]; +} + +void AliTPCseed::RebuildSeed() +{ + // + // rebuild seed to be ready for storing + AliTPCclusterMI cldummy; + cldummy.SetQ(0); + AliTPCTrackPoint pdummy; + pdummy.GetTPoint().fIsShared = 10; + for (Int_t i=0;i<160;i++){ + AliTPCclusterMI * cl0 = fClusterPointer[i]; + AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i); + if (cl0){ + trpoint->GetTPoint() = *(GetTrackPoint(i)); + trpoint->GetCPoint() = *cl0; + trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ())); + } + else{ + *trpoint = pdummy; + trpoint->GetCPoint()= cldummy; + } + + } + +} + + +Double_t AliTPCseed::GetDensityFirst(Int_t n) +{ + // + // + // return cluster for n rows bellow first point + Int_t nfoundable = 1; + Int_t nfound = 1; + for (Int_t i=fLastPoint-1;i>0&&nfoundable0) nfound++; + } + if (nfoundableIsUsed(10)) { + shared++; + continue; + } + if (!plus2) continue; //take also neighborhoud + // + if ( (i>0) && fClusterPointer[i-1]){ + if (fClusterPointer[i-1]->IsUsed(10)) { + shared++; + continue; + } + } + if ( fClusterPointer[i+1]){ + if (fClusterPointer[i+1]->IsUsed(10)) { + shared++; + continue; + } + } + + } + //if (shared>found){ + //Error("AliTPCseed::GetClusterStatistic","problem\n"); + //} +} + + + + + +void AliTPCseed::Reset(Bool_t all) +{ + // + // + SetNumberOfClusters(0); + fNFoundable = 0; + SetChi2(0); + ResetCovariance(); + /* + if (fTrackPoints){ + for (Int_t i=0;i<8;i++){ + delete [] fTrackPoints[i]; + } + delete fTrackPoints; + fTrackPoints =0; + } + */ + + if (all){ + for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); + for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; + } + +} + + +void AliTPCseed::Modify(Double_t factor) +{ + + //------------------------------------------------------------------ + //This function makes a track forget its history :) + //------------------------------------------------------------------ + if (factor<=0) { + ResetCovariance(); + return; + } + fC00*=factor; + fC10*=0; fC11*=factor; + fC20*=0; fC21*=0; fC22*=factor; + fC30*=0; fC31*=0; fC32*=0; fC33*=factor; + fC40*=0; fC41*=0; fC42*=0; fC43*=0; fC44*=factor; + SetNumberOfClusters(0); + fNFoundable =0; + SetChi2(0); + fRemoval = 0; + fCurrentSigmaY2 = 0.000005; + fCurrentSigmaZ2 = 0.000005; + fNoCluster = 0; + //fFirstPoint = 160; + //fLastPoint = 0; +} + + + + +Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const +{ + //----------------------------------------------------------------- + // This function find proloncation of a track to a reference plane x=xk. + // doesn't change internal state of the track + //----------------------------------------------------------------- + + Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1; + + if (TMath::Abs(fP4*xk - fP2) >= 0.999) { + return 0; + } + + // Double_t y1=fP0, z1=fP1; + Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1); + Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2); + + y = fP0; + z = fP1; + //y += dx*(c1+c2)/(r1+r2); + //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3; + + Double_t dy = dx*(c1+c2)/(r1+r2); + Double_t dz = 0; + // + Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1); + /* + if (TMath::Abs(delta)>0.0001){ + dz = fP3*TMath::ASin(delta)/fP4; + }else{ + dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1); + } + */ + // dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4; + dz = fP3*TMath::ASin(delta)/fP4; + // + y+=dy; + z+=dz; + + + return 1; +} + + +//_____________________________________________________________________________ +Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const +{ + //----------------------------------------------------------------- + // This function calculates a predicted chi2 increment. + //----------------------------------------------------------------- + //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2(); + Double_t r00=fErrorY2, r01=0., r11=fErrorZ2; + r00+=fC00; r01+=fC10; r11+=fC11; + + Double_t det=r00*r11 - r01*r01; + if (TMath::Abs(det) < 1.e-10) { + //Int_t n=GetNumberOfClusters(); + //if (n>4) cerr<GetY() - fP0, dz=c->GetZ() - fP1; + + return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; +} + + +//_________________________________________________________________________________________ + + +Int_t AliTPCseed::Compare(const TObject *o) const { + //----------------------------------------------------------------- + // This function compares tracks according to the sector - for given sector according z + //----------------------------------------------------------------- + AliTPCseed *t=(AliTPCseed*)o; + + if (fSort == 0){ + if (t->fRelativeSector>fRelativeSector) return -1; + if (t->fRelativeSectorGetZ(); + Double_t z1 = GetZ(); + if (z2>z1) return 1; + if (z2fC44)/(TMath::Abs(t->GetC())+0.0066); + if (t->fBConstrain) f2=1.2; + + Float_t f1 =1; + f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066); + + if (fBConstrain) f1=1.2; + + if (t->GetNumberOfClusters()*f2 GetY() - fP0, dz=c->GetZ() - fP1; + Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz; + if (TMath::Abs(cur*fX-eta) >= 0.9) { + return 0; + } + + fP0 += k00*dy + k01*dz; + fP1 += k10*dy + k11*dz; + fP2 = eta; + fP3 += k30*dy + k31*dz; + fP4 = cur; + + Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40; + Double_t c12=fC21, c13=fC31, c14=fC41; + + fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11; + fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13; + fC40-=k00*c04+k01*c14; + + fC11-=k10*c01+k11*fC11; + fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13; + fC41-=k10*c04+k11*c14; + + fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13; + fC42-=k20*c04+k21*c14; + + fC33-=k30*c03+k31*c13; + fC43-=k40*c03+k41*c13; + + fC44-=k40*c04+k41*c14; + + Int_t n=GetNumberOfClusters(); + // fIndex[n]=index; + SetNumberOfClusters(n+1); + SetChi2(GetChi2()+chisq); + + return 1; +} + + + +//_____________________________________________________________________________ +void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) { + //----------------------------------------------------------------- + // This funtion calculates dE/dX within the "low" and "up" cuts. + //----------------------------------------------------------------- + + Float_t amp[200]; + Float_t angular[200]; + Float_t weight[200]; + Int_t index[200]; + //Int_t nc = 0; + // TClonesArray & arr = *fPoints; + Float_t meanlog = 100.; + + Float_t mean[4] = {0,0,0,0}; + Float_t sigma[4] = {1000,1000,1000,1000}; + Int_t nc[4] = {0,0,0,0}; + Float_t norm[4] = {1000,1000,1000,1000}; + // + // + fNShared =0; + + for (Int_t of =0; of<4; of++){ + for (Int_t i=of+i1;iIsUsed(10))) continue; + if (cl->IsUsed(11)) { + fNShared++; + continue; + } + Int_t type = cl->GetType(); + //if (point->fIsShared){ + // fNShared++; + // continue; + //} + //if (pointm) + // if (pointm->fIsShared) continue; + //if (pointp) + // if (pointp->fIsShared) continue; + + if (type<0) continue; + //if (type>10) continue; + //if (point->GetErrY()==0) continue; + //if (point->GetErrZ()==0) continue; + + //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY(); + //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ(); + //if ((ddy*ddy+ddz*ddz)>10) continue; + + + // if (point->GetCPoint().GetMax()<5) continue; + if (cl->GetMax()<5) continue; + Float_t angley = point->GetAngleY(); + Float_t anglez = point->GetAngleZ(); + + Float_t rsigmay2 = point->GetSigmaY(); + Float_t rsigmaz2 = point->GetSigmaZ(); + /* + Float_t ns = 1.; + if (pointm){ + rsigmay += pointm->GetTPoint().GetSigmaY(); + rsigmaz += pointm->GetTPoint().GetSigmaZ(); + ns+=1.; + } + if (pointp){ + rsigmay += pointp->GetTPoint().GetSigmaY(); + rsigmaz += pointp->GetTPoint().GetSigmaZ(); + ns+=1.; + } + rsigmay/=ns; + rsigmaz/=ns; + */ + + Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2); + + Float_t ampc = 0; // normalization to the number of electrons + if (i>64){ + // ampc = 1.*point->GetCPoint().GetMax(); + ampc = 1.*cl->GetMax(); + //ampc = 1.*point->GetCPoint().GetQ(); + // AliTPCClusterPoint & p = point->GetCPoint(); + // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5); + // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; + //Float_t dz = + // TMath::Abs( Int_t(iz) - iz + 0.5); + //ampc *= 1.15*(1-0.3*dy); + //ampc *= 1.15*(1-0.3*dz); + // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ())); + //ampc *=zfactor; + } + else{ + //ampc = 1.0*point->GetCPoint().GetMax(); + ampc = 1.0*cl->GetMax(); + //ampc = 1.0*point->GetCPoint().GetQ(); + //AliTPCClusterPoint & p = point->GetCPoint(); + // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5); + //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; + //Float_t dz = + // TMath::Abs( Int_t(iz) - iz + 0.5); + + //ampc *= 1.15*(1-0.3*dy); + //ampc *= 1.15*(1-0.3*dz); + // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ())); + //ampc *=zfactor; + + } + ampc *= 2.0; // put mean value to channel 50 + //ampc *= 0.58; // put mean value to channel 50 + Float_t w = 1.; + // if (type>0) w = 1./(type/2.-0.5); + // Float_t z = TMath::Abs(cl->GetZ()); + if (i<64) { + ampc /= 0.6; + //ampc /= (1+0.0008*z); + } else + if (i>128){ + ampc /=1.5; + //ampc /= (1+0.0008*z); + }else{ + //ampc /= (1+0.0008*z); + } + + if (type<0) { //amp at the border - lower weight + // w*= 2.; + + continue; + } + if (rsigma>1.5) ampc/=1.3; // if big backround + amp[nc[of]] = ampc; + angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez); + weight[nc[of]] = w; + nc[of]++; + } + + TMath::Sort(nc[of],amp,index,kFALSE); + Float_t sumamp=0; + Float_t sumamp2=0; + Float_t sumw=0; + //meanlog = amp[index[Int_t(nc[of]*0.33)]]; + meanlog = 50; + for (Int_t i=int(nc[of]*low+0.5);i0.1) + sigma[of] = TMath::Sqrt(sigma[of]); + else + sigma[of] = 1000; + + mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; + //mean *=(1-0.02*(sigma/(mean*0.17)-1.)); + //mean *=(1-0.1*(norm-1.)); + } + } + + Float_t dedx =0; + fSdEdx =0; + fMAngular =0; + // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1)); + // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1)); + + + // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ + // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1]))); + + Int_t norm2 = 0; + Int_t norm3 = 0; + for (Int_t i =0;i<4;i++){ + if (nc[i]>2&&nc[i]<1000){ + dedx += mean[i] *nc[i]; + fSdEdx += sigma[i]*(nc[i]-2); + fMAngular += norm[i] *nc[i]; + norm2 += nc[i]; + norm3 += nc[i]-2; + } + fDEDX[i] = mean[i]; + fSDEDX[i] = sigma[i]; + fNCDEDX[i]= nc[i]; + } + + if (norm3>0){ + dedx /=norm2; + fSdEdx /=norm3; + fMAngular/=norm2; + } + else{ + SetdEdx(0); + return; + } + // Float_t dedx1 =dedx; + /* + dedx =0; + for (Int_t i =0;i<4;i++){ + if (nc[i]>2&&nc[i]<1000){ + mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.)); + dedx += mean[i] *nc[i]; + } + fDEDX[i] = mean[i]; + } + dedx /= norm2; + */ + + + SetdEdx(dedx); + + //mi deDX + + + + //Very rough PID + Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt())); + + if (p<0.6) { + if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;} + if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;} + SetMass(0.93827); return; + } + + if (p<1.2) { + if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;} + SetMass(0.93827); return; + } + + SetMass(0.13957); return; + +} +Double_t AliTPCseed::Bethe(Double_t bg){ + // + // This is the Bethe-Bloch function normalised to 1 at the minimum + // + Double_t bg2=bg*bg; + Double_t bethe; + if (bg<3.5e1) + bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2)); + else // Density effect ( approximately :) + bethe=1.15*(1.+ bg2)/bg2*(log(3.5*5940*bg) - bg2/(1.+ bg2)); + return bethe/11.091; +} + +void AliTPCseed::CookPID() +{ + // + // cook PID information according dEdx + // + Double_t fRange = 10.; + Double_t fRes = 0.1; + Double_t fMIP = 47.; + // + Int_t ns=AliPID::kSPECIES; + Double_t sumr =0; + for (Int_t j=0; j0.001){ + if (TMath::Abs(dedx-bethe) > fRange*sigma) { + fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; + sumr+=fTPCr[j]; + continue; + } + fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; + sumr+=fTPCr[j]; + } + else{ + fTPCr[j]=1.; + sumr+=fTPCr[j]; + } + } + for (Int_t j=0; jfIsShared){ + fNShared++; + continue; + } + Int_t type = point->GetCPoint().GetType(); + if (type<0) continue; + if (point->GetCPoint().GetMax()<5) continue; + Float_t angley = point->GetTPoint().GetAngleY(); + Float_t anglez = point->GetTPoint().GetAngleZ(); + Float_t rsigmay = point->GetCPoint().GetSigmaY(); + Float_t rsigmaz = point->GetCPoint().GetSigmaZ(); + Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz); + + Float_t ampc = 0; // normalization to the number of electrons + if (i>64){ + ampc = point->GetCPoint().GetMax(); + } + else{ + ampc = point->GetCPoint().GetMax(); + } + ampc *= 2.0; // put mean value to channel 50 + // ampc *= 0.565; // put mean value to channel 50 + + Float_t w = 1.; + Float_t z = TMath::Abs(point->GetCPoint().GetZ()); + if (i<64) { + ampc /= 0.63; + } else + if (i>128){ + ampc /=1.51; + } + if (type<0) { //amp at the border - lower weight + continue; + } + if (rsigma>1.5) ampc/=1.3; // if big backround + angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez); + amp[i] = ampc/angular[i]; + weight[i] = w; + anc++; + } + + TMath::Sort(159,amp,index,kFALSE); + for (Int_t i=int(anc*low+0.5);i0.1) + sigma[of] = TMath::Sqrt(sigma[of]); + else + sigma[of] = 1000; + mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; + } + } + + Float_t dedx =0; + fSdEdx =0; + fMAngular =0; + // + Int_t norm2 = 0; + Int_t norm3 = 0; + Float_t www[3] = {12.,14.,17.}; + //Float_t www[3] = {1.,1.,1.}; + + for (Int_t i =0;i<3;i++){ + if (nc[i]>2&&nc[i]<1000){ + dedx += mean[i] *nc[i]*www[i]/sigma[i]; + fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i]; + fMAngular += norm[i] *nc[i]; + norm2 += nc[i]*www[i]/sigma[i]; + norm3 += (nc[i]-2)*www[i]/sigma[i]; + } + fDEDX[i] = mean[i]; + fSDEDX[i] = sigma[i]; + fNCDEDX[i]= nc[i]; + } + + if (norm3>0){ + dedx /=norm2; + fSdEdx /=norm3; + fMAngular/=norm2; + } + else{ + SetdEdx(0); + return; + } + // Float_t dedx1 =dedx; + + dedx =0; + Float_t norm4 = 0; + for (Int_t i =0;i<3;i++){ + if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){ + //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.)); + dedx += mean[i] *(nc[i])/(sigma[i]); + norm4 += (nc[i])/(sigma[i]); + } + fDEDX[i] = mean[i]; + } + if (norm4>0) dedx /= norm4; + + + + SetdEdx(dedx); + + //mi deDX + +} +*/ diff --git a/TPC/AliTPCseed.h b/TPC/AliTPCseed.h new file mode 100644 index 00000000000..1b01fbb15d2 --- /dev/null +++ b/TPC/AliTPCseed.h @@ -0,0 +1,120 @@ +#ifndef ALITPCSEED_H +#define ALITPCSEED_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + +/* $Id$ */ + +//------------------------------------------------------- +// TPC seed +// Class needed for TPC parallel tracking +// +// Origin: +//------------------------------------------------------- + +#include + +#include "AliTPCtrack.h" +#include "AliComplexCluster.h" +#include "AliPID.h" + +class TFile; +class AliTPCParam; +class AliTPCseed; +class AliTPCclusterMI; +class AliTPCTrackerPoint; +class AliESD; +class TClonesArray; + +class AliTPCseed : public AliTPCtrack { + friend class AliTPCtrackerMI; + public: + AliTPCseed(); + virtual ~AliTPCseed(); + AliTPCseed(const AliTPCtrack &t); + AliTPCseed(const AliTPCseed &s); + //AliTPCseed(const AliTPCseed &t, Double_t a); + AliTPCseed(UInt_t index, const Double_t xx[5], + const Double_t cc[15], Double_t xr, Double_t alpha); + Int_t Compare(const TObject *o) const; + void Reset(Bool_t all = kTRUE); + Int_t GetProlongation(Double_t xr, Double_t &y, Double_t & z) const; + virtual Double_t GetPredictedChi2(const AliTPCclusterMI *cluster) const; + virtual Int_t Update(const AliTPCclusterMI* c, Double_t chi2, UInt_t i); + AliTPCTrackerPoint * GetTrackPoint(Int_t i); + void RebuildSeed(); // rebuild seed to be ready for storing + Double_t GetDensityFirst(Int_t n); + Double_t GetSigma2C() const {return fC44;}; + void GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2); + + void Modify(Double_t factor); + void SetClusterIndex2(Int_t row, Int_t index) { + fIndex[row] = index; + } + Int_t GetClusterIndex2(Int_t row) const { + return fIndex[row]; + } + Int_t GetClusterSector(Int_t row) const { + Int_t pica = -1; + if (fIndex[row]>=0) pica = ((fIndex[row]&0xff000000)>>24); + return pica; + } + + void SetErrorY2(Float_t sy2){fErrorY2=sy2;} + void SetErrorZ2(Float_t sz2){fErrorZ2=sz2;} + void CookdEdx(Double_t low=0.05, Double_t up=0.70, Int_t i1=0, Int_t i2=159, Bool_t onlyused = kFALSE); + void CookPID(); + Double_t Bethe(Double_t bg); // return bethe-bloch + // void CookdEdx2(Double_t low=0.05, Double_t up=0.70); + Bool_t IsActive() const { return !(fRemoval);} + void Desactivate(Int_t reason){ fRemoval = reason;} + // + // + private: + // AliTPCseed & operator = (const AliTPCseed &) + // {::Fatal("= operator","Not Implemented\n");return *this;} + AliESDtrack * fEsd; //! + AliTPCclusterMI* fClusterPointer[160]; //! array of cluster pointers - + TClonesArray * fPoints; // array with points along the track + TClonesArray * fEPoints; // array with exact points - calculated in special macro not used in tracking + //---CURRENT VALUES + Int_t fRow; //!current row number + Int_t fSector; //!current sector number + Int_t fRelativeSector; //! index of current relative sector + Float_t fCurrentSigmaY2; //!expected current cluster sigma Y + Float_t fCurrentSigmaZ2; //!expected current cluster sigma Z + Float_t fErrorY2; //!sigma of current cluster + Float_t fErrorZ2; //!sigma of current cluster + AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation + Int_t fCurrentClusterIndex1; //! index of the current cluster + Bool_t fInDead; //! indicate if the track is in dead zone + Bool_t fIsSeeding; //!indicates if it is proces of seeading + Int_t fNoCluster; //!indicates number of rows without clusters + Int_t fSort; //!indicate criteria for sorting + Bool_t fBSigned; //indicates that clusters of this trackes are signed to be used + // + // + Float_t fDEDX[4]; // dedx according padrows + Float_t fSDEDX[4]; // sdedx according padrows + Int_t fNCDEDX[4]; // number of clusters for dedx measurment + Double_t fTPCr[AliPID::kSPECIES]; // rough PID according TPC + // + Int_t fSeedType; //seeding type + Int_t fSeed1; //first row for seeding + Int_t fSeed2; //last row for seeding + Int_t fOverlapLabels[12]; //track labels and the length of the overlap + Float_t fMAngular; // mean angular factor + Char_t fCircular; // indicates curlin track + AliTPCTrackerPoint fTrackPoints[160]; //!track points - array track points + + ClassDef(AliTPCseed,1) +}; + + + + + +#endif + + diff --git a/TPC/AliTPCtrack.cxx b/TPC/AliTPCtrack.cxx index f514d185c27..32c780dff81 100644 --- a/TPC/AliTPCtrack.cxx +++ b/TPC/AliTPCtrack.cxx @@ -38,7 +38,7 @@ AliTPCtrack::AliTPCtrack(): AliKalmanTrack() fX = fP0 = fP1 = fP2 = fP3 = fP3 = fP4 = 0.0; fAlpha = fdEdx = 0.0; fNumber = 0; // [SR, 01.04.2003] - for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0; + for (Int_t i=0; i<3;i++) {fKinkIndexes[i]=0; fV0Indexes[i]=0;} } //_________________________________________________________________________ @@ -77,6 +77,7 @@ const Double_t cc[15], Double_t xref, Double_t alpha) : AliKalmanTrack() { fTrackType = 0; fLab2 = 0; for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0; + for (Int_t i=0; i<3;i++) fV0Indexes[i]=0; } //_____________________________________________________________________________ @@ -88,6 +89,7 @@ AliTPCtrack::AliTPCtrack(const AliESDtrack& t) : AliKalmanTrack() { SetLabel(t.GetLabel()); SetMass(t.GetMass()); for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i); + for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.GetV0Index(i); fdEdx = t.GetTPCsignal(); fAlpha = t.GetAlpha(); @@ -165,6 +167,7 @@ AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : AliKalmanTrack(t) { fTrackType = t.fTrackType; fLab2 = t.fLab2; for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.fKinkIndexes[i]; + for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.fV0Indexes[i]; } //_____________________________________________________________________________ diff --git a/TPC/AliTPCtrack.h b/TPC/AliTPCtrack.h index c339ce07d00..a99cdaddc97 100644 --- a/TPC/AliTPCtrack.h +++ b/TPC/AliTPCtrack.h @@ -108,8 +108,10 @@ public: Double_t GetD(Double_t x=0, Double_t y=0) const; AliExternalTrackParam & GetReference(){ return fReference;} void UpdateReference(){ new (&fReference) AliExternalTrackParam(*this);} - Int_t GetKinkIndex(Int_t i) const{ return fKinkIndexes[i];} + Int_t GetKinkIndex(Int_t i) const{ return fKinkIndexes[i];} Int_t* GetKinkIndexes() { return fKinkIndexes;} + Int_t GetV0Index(Int_t i) const{ return fV0Indexes[i];} + Int_t* GetV0Indexes() { return fV0Indexes;} protected: Double_t fX; // X-coordinate of this track (reference plane) Double_t fAlpha; // Rotation angle the local (TPC sector) @@ -147,6 +149,7 @@ protected: AliExternalTrackParam fReference; // track parameters at the middle of the chamber Float_t fKinkPoint[12]; //radius, of kink, dfi and dtheta Int_t fKinkIndexes[3]; // kink indexes - minus = mother + daughter + Int_t fV0Indexes[3]; // kink indexes - minus = mother + daughter ClassDef(AliTPCtrack,2) // Time Projection Chamber reconstructed tracks }; diff --git a/TPC/AliTPCtrackerMI.cxx b/TPC/AliTPCtrackerMI.cxx index f35ec1290eb..5ef0349b538 100644 --- a/TPC/AliTPCtrackerMI.cxx +++ b/TPC/AliTPCtrackerMI.cxx @@ -46,13 +46,16 @@ #include "AliTPCReconstructor.h" #include "AliTPCclusterMI.h" #include "AliTPCpolyTrack.h" -#include "AliTPCreco.h" +#include "AliTPCreco.h" +#include "AliTPCseed.h" #include "AliTPCtrackerMI.h" #include "TStopwatch.h" +#include "AliTPCReconstructor.h" +#include "AliESDkink.h" +#include "AliPID.h" #include "TTreeStream.h" // -ClassImp(AliTPCseed) ClassImp(AliTPCtrackerMI) @@ -255,6 +258,7 @@ AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2) fNewIO =0; fDebug =0; fEvent =0; + fDebugStreamer = new TTreeSRedirector("TPCdebug.root"); } //________________________________________________________________________ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): @@ -283,6 +287,7 @@ AliTPCtrackerMI::~AliTPCtrackerMI() { fSeeds->Delete(); delete fSeeds; } + if (fDebugStreamer) delete fDebugStreamer; } void AliTPCtrackerMI::SetIO() @@ -372,6 +377,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); iotrack.SetTPCPoints(pt->GetPoints()); iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + // iotrack.SetTPCpid(pt->fTPCr); //iotrack.SetTPCindex(i); fEvent->AddTrack(&iotrack); continue; @@ -383,6 +390,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) iotrack.SetTPCPoints(pt->GetPoints()); //iotrack.SetTPCindex(i); iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + // iotrack.SetTPCpid(pt->fTPCr); fEvent->AddTrack(&iotrack); continue; } @@ -398,6 +407,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) //iotrack.SetTPCindex(i); iotrack.SetTPCPoints(pt->GetPoints()); iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + //iotrack.SetTPCpid(pt->fTPCr); fEvent->AddTrack(&iotrack); continue; } @@ -413,6 +424,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); iotrack.SetTPCPoints(pt->GetPoints()); iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + //iotrack.SetTPCpid(pt->fTPCr); //iotrack.SetTPCindex(i); fEvent->AddTrack(&iotrack); continue; @@ -427,6 +440,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); iotrack.SetTPCPoints(pt->GetPoints()); iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + //iotrack.SetTPCpid(pt->fTPCr); //iotrack.SetTPCindex(i); fEvent->AddTrack(&iotrack); continue; @@ -444,6 +459,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin); iotrack.SetTPCPoints(pt->GetPoints()); iotrack.SetKinkIndexes(pt->GetKinkIndexes()); + iotrack.SetV0Indexes(pt->GetV0Indexes()); + // iotrack.SetTPCpid(pt->fTPCr); //iotrack.SetTPCindex(i); fEvent->AddTrack(&iotrack); continue; @@ -895,160 +912,6 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ -void AliTPCseed::Reset(Bool_t all) -{ - // - // - SetNumberOfClusters(0); - fNFoundable = 0; - SetChi2(0); - ResetCovariance(); - /* - if (fTrackPoints){ - for (Int_t i=0;i<8;i++){ - delete [] fTrackPoints[i]; - } - delete fTrackPoints; - fTrackPoints =0; - } - */ - - if (all){ - for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); - for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; - } - -} - - -void AliTPCseed::Modify(Double_t factor) -{ - - //------------------------------------------------------------------ - //This function makes a track forget its history :) - //------------------------------------------------------------------ - if (factor<=0) { - ResetCovariance(); - return; - } - fC00*=factor; - fC10*=0; fC11*=factor; - fC20*=0; fC21*=0; fC22*=factor; - fC30*=0; fC31*=0; fC32*=0; fC33*=factor; - fC40*=0; fC41*=0; fC42*=0; fC43*=0; fC44*=factor; - SetNumberOfClusters(0); - fNFoundable =0; - SetChi2(0); - fRemoval = 0; - fCurrentSigmaY2 = 0.000005; - fCurrentSigmaZ2 = 0.000005; - fNoCluster = 0; - //fFirstPoint = 160; - //fLastPoint = 0; -} - - - - -Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const -{ - //----------------------------------------------------------------- - // This function find proloncation of a track to a reference plane x=xk. - // doesn't change internal state of the track - //----------------------------------------------------------------- - - Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1; - - if (TMath::Abs(fP4*xk - fP2) >= 0.999) { - return 0; - } - - // Double_t y1=fP0, z1=fP1; - Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1); - Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2); - - y = fP0; - z = fP1; - //y += dx*(c1+c2)/(r1+r2); - //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3; - - Double_t dy = dx*(c1+c2)/(r1+r2); - Double_t dz = 0; - // - Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1); - /* - if (TMath::Abs(delta)>0.0001){ - dz = fP3*TMath::ASin(delta)/fP4; - }else{ - dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1); - } - */ - dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4; - // - y+=dy; - z+=dz; - - - return 1; -} - - -//_____________________________________________________________________________ -Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const -{ - //----------------------------------------------------------------- - // This function calculates a predicted chi2 increment. - //----------------------------------------------------------------- - //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2(); - Double_t r00=fErrorY2, r01=0., r11=fErrorZ2; - r00+=fC00; r01+=fC10; r11+=fC11; - - Double_t det=r00*r11 - r01*r01; - if (TMath::Abs(det) < 1.e-10) { - Int_t n=GetNumberOfClusters(); - if (n>4) cerr<GetY() - fP0, dz=c->GetZ() - fP1; - - return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; -} - - -//_________________________________________________________________________________________ - - -Int_t AliTPCseed::Compare(const TObject *o) const { - //----------------------------------------------------------------- - // This function compares tracks according to the sector - for given sector according z - //----------------------------------------------------------------- - AliTPCseed *t=(AliTPCseed*)o; - - if (fSort == 0){ - if (t->fRelativeSector>fRelativeSector) return -1; - if (t->fRelativeSectorGetZ(); - Double_t z1 = GetZ(); - if (z2>z1) return 1; - if (z2fC44)/(TMath::Abs(t->GetC())+0.0066); - if (t->fBConstrain) f2=1.2; - - Float_t f1 =1; - f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066); - - if (fBConstrain) f1=1.2; - - if (t->GetNumberOfClusters()*f2 GetY() - fP0, dz=c->GetZ() - fP1; - Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz; - if (TMath::Abs(cur*fX-eta) >= 0.9) { - return 0; - } - - fP0 += k00*dy + k01*dz; - fP1 += k10*dy + k11*dz; - fP2 = eta; - fP3 += k30*dy + k31*dz; - fP4 = cur; - - Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40; - Double_t c12=fC21, c13=fC31, c14=fC41; - - fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11; - fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13; - fC40-=k00*c04+k01*c14; - - fC11-=k10*c01+k11*fC11; - fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13; - fC41-=k10*c04+k11*c14; - - fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13; - fC42-=k20*c04+k21*c14; - - fC33-=k30*c03+k31*c13; - fC43-=k40*c03+k41*c13; - - fC44-=k40*c04+k41*c14; - - Int_t n=GetNumberOfClusters(); - // fIndex[n]=index; - SetNumberOfClusters(n+1); - SetChi2(GetChi2()+chisq); - - return 1; -} - - - //_____________________________________________________________________________ Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1, Double_t x2,Double_t y2, @@ -4083,6 +3887,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) // 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]; @@ -4095,6 +3900,8 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) 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(); // @@ -4106,6 +3913,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) usage[i]=0; AliTPCseed* track = (AliTPCseed*)array->At(i); if (!track) continue; + track->fCircular =0; shared[i] = kFALSE; track->UpdatePoints(); if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){ @@ -4128,7 +3936,8 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) } z0[i]=1000; circular[i]= kFALSE; - if (track->GetProlongation(0,y,z)) z0[i] = TMath::Abs(z); + if (track->GetProlongation(0,y,z)) z0[i] = z; + dca[i] = track->GetD(0,0); } // // @@ -4141,7 +3950,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) // // Find circling track - // TTreeSRedirector cstream("circling.root"); + TTreeSRedirector &cstream = *fDebugStreamer; // for (Int_t i0=0;i0At(i0); @@ -4151,17 +3960,23 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) for (Int_t i1=i0+1;i1At(i1); if (!track1) continue; - if (TMath::Abs(1./track1->fP4)>200) continue; - if (track1->fN<40) continue; - if (z0[i0]<20&&z0[i1]<20) continue; - if (track1->fP4*track0->fP4>0) continue; - if (track1->fP3*track0->fP3>0) 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)fP4)) continue; //returning - lower momenta if (track1->fBConstrain&&TMath::Abs(track0->fP4)fP4)) continue; //returning - lower momenta // - if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue; + 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); @@ -4171,7 +3986,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) 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>30) continue; if (delta>rmean*0.25) continue; if (TMath::Abs(r0-r1)/rmean>0.3) continue; // @@ -4198,44 +4013,66 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) */ if (npoints>0){ Int_t ibest=0; - helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],5); + 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],5); + helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2); if (deltah[1]fLab<fLab<< //0-1 - track0->fP3<fP3<< //2-3 - track0->fP4<fP4<< //4-5 - delta<fFirstPoint<fFirstPoint<< //15-16 - track0->fLastPoint<fLastPoint<< //17-18 - track0<6) continue; + if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue; Bool_t sign =kFALSE; - if (deltabest<3&&hangles[2]>3.06) sign =kTRUE; - if (TMath::Min(TMath::Abs(track0->GetD()),TMath::Abs(track1->GetD()))>10&&deltabest<6&&hangles[2]>3.) sign= kTRUE; + if (hangles[2]>3.06) sign =kTRUE; + // if (sign){ circular[i0] = kTRUE; - circular[i1] = kTRUE; + circular[i1] = kTRUE; + if (TMath::Abs(track0->fP4)fP4)){ + track0->fCircular += 1; + track1->fCircular += 2; + } + else{ + track1->fCircular += 1; + track0->fCircular += 2; + } + } + if (sign){ + //debug stream + cstream<<"Curling"<< + "lab0="<fLab<< + "lab1="<fLab<< + "Tr0.="<At(i); + if (!track0) continue; + track0->CookdEdx(0.02,0.6); + track0->CookPID(); + } // for (Int_t i=0;iAt(i); @@ -4681,6 +4524,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) delete [] mothers; // // + delete [] dca; delete []circular; delete []shared; delete []quality; @@ -4702,6 +4546,342 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd) 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;iAt(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;iAt(i); + if (!track0) continue; + cstream<<"Tracks"<< + "Tr0.="<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;jAt(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]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])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]cdist1) continue; + } + helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]); + if (radius[iclosest]cradius1 || delta[iclosest]>cdist1) continue; + // + Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex); + if (pointAnglefTPCr[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()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]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;iAt(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(); + cstream<<"V02"<< + "Event="<GetEntriesFast(); +// for (Int_t i=0; iUncheckedAt(i), &t=*pt; +// if (!pt) continue; +// Int_t nc=t.GetNumberOfClusters(); +// if (nc<15) { +// delete fSeeds->RemoveAt(i); +// continue; +// } +// if (pt->GetKinkIndexes()[0]!=0) continue; // ignore kink candidates +// if (nc>100) continue; // hopefully, enough for ITS +// AliTPCseed *seed2 = new AliTPCseed(*pt); +// //seed2->Reset(kFALSE); +// //pt->ResetCovariance(); +// seed2->Modify(1); +// FollowBackProlongation(*seed2,158); +// //seed2->Reset(kFALSE); +// seed2->Modify(10); +// FollowProlongation(*seed2,0); +// TTreeSRedirector &cstream = *fDebugStreamer; +// cstream<<"Crefit"<< +// "Tr0.="<fN>pt->fN){ +// delete fSeeds->RemoveAt(i); +// fSeeds->AddAt(seed2,i); +// }else{ +// delete seed2; +// } +// } +// RemoveUsed2(fSeeds,0.6,0.6,50); + +// FindV0s(fSeeds,fEvent); //RemoveDouble(fSeeds,0.2,0.6,11); - // RemoveUsed(fSeeds,0.5,0.5,6); // - Int_t nseed=fSeeds->GetEntriesFast(); Int_t found = 0; for (Int_t i=0; iUncheckedAt(i), &t=*pt; @@ -6116,654 +6331,12 @@ AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t -AliTPCseed::AliTPCseed():AliTPCtrack(){ - // - fRow=0; - fRemoval =0; - for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); - for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; - for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0; - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared =0; - fRemoval = 0; - fSort =0; - fFirstPoint =0; - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; - fCurrentSigmaY2=0; - fCurrentSigmaZ2=0; -} -AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){ - //--------------------- - // dummy copy constructor - //------------------------- - for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i]; - for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i]; - fPoints = 0; - fEPoints = 0; -} -AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){ - // - //copy constructor - fPoints = 0; - fEPoints = 0; - fNShared =0; - // fTrackPoints =0; - fRemoval =0; - fSort =0; - for (Int_t i=0;i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i); - - for (Int_t i=0;i<160;i++) { - fClusterPointer[i] = 0; - Int_t index = t.GetClusterIndex(i); - if (index>=-1){ - SetClusterIndex2(i,index); - } - else{ - SetClusterIndex2(i,-3); - } - } - fFirstPoint =0; - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; - fCurrentSigmaY2=0; - fCurrentSigmaZ2=0; -} -/* -AliTPCseed::AliTPCseed(const AliTPCtrack &t, Double_t a):AliTPCtrack(t,a){ - // - //copy constructor - fRow=0; - for (Int_t i=0;i<160;i++) { - fClusterPointer[i] = 0; - Int_t index = t.GetClusterIndex(i); - SetClusterIndex2(i,index); - } - for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0; - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared =0; - // fTrackPoints =0; - fRemoval =0; - fSort = 0; - fFirstPoint =0; - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; - fCurrentSigmaY2=0; - fCurrentSigmaZ2=0; - -} -*/ - -AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15], - Double_t xr, Double_t alpha): - AliTPCtrack(index, xx, cc, xr, alpha) { - // - // - //constructor - fRow =0; - for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3); - for (Int_t i=0;i<160;i++) fClusterPointer[i]=0; - for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0; - fPoints = 0; - fEPoints = 0; - fNFoundable =0; - fNShared = 0; - // fTrackPoints =0; - fRemoval =0; - fSort =0; - fFirstPoint =0; - // fHelixIn = new TClonesArray("AliHelix",0); - //fHelixOut = new TClonesArray("AliHelix",0); - fNoCluster =0; - fBSigned = kFALSE; - fSeed1 =-1; - fSeed2 =-1; - fCurrentCluster =0; - fCurrentSigmaY2=0; - fCurrentSigmaZ2=0; -} - -AliTPCseed::~AliTPCseed(){ - // - // destructor - if (fPoints) delete fPoints; - fPoints =0; - if (fEPoints) delete fEPoints; - fEPoints = 0; - fNoCluster =0; -} - -AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i) -{ - // - // - return &fTrackPoints[i]; -} - -void AliTPCseed::RebuildSeed() -{ - // - // rebuild seed to be ready for storing - AliTPCclusterMI cldummy; - cldummy.SetQ(0); - AliTPCTrackPoint pdummy; - pdummy.GetTPoint().fIsShared = 10; - for (Int_t i=0;i<160;i++){ - AliTPCclusterMI * cl0 = fClusterPointer[i]; - AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i); - if (cl0){ - trpoint->GetTPoint() = *(GetTrackPoint(i)); - trpoint->GetCPoint() = *cl0; - trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ())); - } - else{ - *trpoint = pdummy; - trpoint->GetCPoint()= cldummy; - } - - } - -} - - -Double_t AliTPCseed::GetDensityFirst(Int_t n) -{ - // - // - // return cluster for n rows bellow first point - Int_t nfoundable = 1; - Int_t nfound = 1; - for (Int_t i=fLastPoint-1;i>0&&nfoundable0) nfound++; - } - if (nfoundableIsUsed(10)) { - shared++; - continue; - } - if (!plus2) continue; //take also neighborhoud - // - if ( (i>0) && fClusterPointer[i-1]){ - if (fClusterPointer[i-1]->IsUsed(10)) { - shared++; - continue; - } - } - if ( fClusterPointer[i+1]){ - if (fClusterPointer[i+1]->IsUsed(10)) { - shared++; - continue; - } - } - - } - //if (shared>found){ - //Error("AliTPCseed::GetClusterStatistic","problem\n"); - //} -} - -//_____________________________________________________________________________ -void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) { - //----------------------------------------------------------------- - // This funtion calculates dE/dX within the "low" and "up" cuts. - //----------------------------------------------------------------- - - Float_t amp[200]; - Float_t angular[200]; - Float_t weight[200]; - Int_t index[200]; - //Int_t nc = 0; - // TClonesArray & arr = *fPoints; - Float_t meanlog = 100.; - - Float_t mean[4] = {0,0,0,0}; - Float_t sigma[4] = {1000,1000,1000,1000}; - Int_t nc[4] = {0,0,0,0}; - Float_t norm[4] = {1000,1000,1000,1000}; - // - // - fNShared =0; - - for (Int_t of =0; of<4; of++){ - for (Int_t i=of+i1;iIsUsed(10))) continue; - if (cl->IsUsed(11)) { - fNShared++; - continue; - } - Int_t type = cl->GetType(); - //if (point->fIsShared){ - // fNShared++; - // continue; - //} - //if (pointm) - // if (pointm->fIsShared) continue; - //if (pointp) - // if (pointp->fIsShared) continue; - - if (type<0) continue; - //if (type>10) continue; - //if (point->GetErrY()==0) continue; - //if (point->GetErrZ()==0) continue; - - //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY(); - //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ(); - //if ((ddy*ddy+ddz*ddz)>10) continue; - - - // if (point->GetCPoint().GetMax()<5) continue; - if (cl->GetMax()<5) continue; - Float_t angley = point->GetAngleY(); - Float_t anglez = point->GetAngleZ(); - - Float_t rsigmay2 = point->GetSigmaY(); - Float_t rsigmaz2 = point->GetSigmaZ(); - /* - Float_t ns = 1.; - if (pointm){ - rsigmay += pointm->GetTPoint().GetSigmaY(); - rsigmaz += pointm->GetTPoint().GetSigmaZ(); - ns+=1.; - } - if (pointp){ - rsigmay += pointp->GetTPoint().GetSigmaY(); - rsigmaz += pointp->GetTPoint().GetSigmaZ(); - ns+=1.; - } - rsigmay/=ns; - rsigmaz/=ns; - */ - - Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2); - - Float_t ampc = 0; // normalization to the number of electrons - if (i>64){ - // ampc = 1.*point->GetCPoint().GetMax(); - ampc = 1.*cl->GetMax(); - //ampc = 1.*point->GetCPoint().GetQ(); - // AliTPCClusterPoint & p = point->GetCPoint(); - // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5); - // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; - //Float_t dz = - // TMath::Abs( Int_t(iz) - iz + 0.5); - //ampc *= 1.15*(1-0.3*dy); - //ampc *= 1.15*(1-0.3*dz); - // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ())); - //ampc *=zfactor; - } - else{ - //ampc = 1.0*point->GetCPoint().GetMax(); - ampc = 1.0*cl->GetMax(); - //ampc = 1.0*point->GetCPoint().GetQ(); - //AliTPCClusterPoint & p = point->GetCPoint(); - // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5); - //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566; - //Float_t dz = - // TMath::Abs( Int_t(iz) - iz + 0.5); - - //ampc *= 1.15*(1-0.3*dy); - //ampc *= 1.15*(1-0.3*dz); - // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ())); - //ampc *=zfactor; - - } - ampc *= 2.0; // put mean value to channel 50 - //ampc *= 0.58; // put mean value to channel 50 - Float_t w = 1.; - // if (type>0) w = 1./(type/2.-0.5); - // Float_t z = TMath::Abs(cl->GetZ()); - if (i<64) { - ampc /= 0.6; - //ampc /= (1+0.0008*z); - } else - if (i>128){ - ampc /=1.5; - //ampc /= (1+0.0008*z); - }else{ - //ampc /= (1+0.0008*z); - } - - if (type<0) { //amp at the border - lower weight - // w*= 2.; - - continue; - } - if (rsigma>1.5) ampc/=1.3; // if big backround - amp[nc[of]] = ampc; - angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez); - weight[nc[of]] = w; - nc[of]++; - } - - TMath::Sort(nc[of],amp,index,kFALSE); - Float_t sumamp=0; - Float_t sumamp2=0; - Float_t sumw=0; - //meanlog = amp[index[Int_t(nc[of]*0.33)]]; - meanlog = 50; - for (Int_t i=int(nc[of]*low+0.5);i0.1) - sigma[of] = TMath::Sqrt(sigma[of]); - else - sigma[of] = 1000; - - mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; - //mean *=(1-0.02*(sigma/(mean*0.17)-1.)); - //mean *=(1-0.1*(norm-1.)); - } - } - - Float_t dedx =0; - fSdEdx =0; - fMAngular =0; - // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1)); - // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1)); - - - // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ - // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1]))); - - Int_t norm2 = 0; - Int_t norm3 = 0; - for (Int_t i =0;i<4;i++){ - if (nc[i]>2&&nc[i]<1000){ - dedx += mean[i] *nc[i]; - fSdEdx += sigma[i]*(nc[i]-2); - fMAngular += norm[i] *nc[i]; - norm2 += nc[i]; - norm3 += nc[i]-2; - } - fDEDX[i] = mean[i]; - fSDEDX[i] = sigma[i]; - fNCDEDX[i]= nc[i]; - } - - if (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; - } - else{ - SetdEdx(0); - return; - } - // Float_t dedx1 =dedx; - /* - dedx =0; - for (Int_t i =0;i<4;i++){ - if (nc[i]>2&&nc[i]<1000){ - mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.)); - dedx += mean[i] *nc[i]; - } - fDEDX[i] = mean[i]; - } - dedx /= norm2; - */ - - - SetdEdx(dedx); - - //mi deDX - - - - //Very rough PID - Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt())); - - if (p<0.6) { - if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;} - if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;} - SetMass(0.93827); return; - } - - if (p<1.2) { - if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;} - SetMass(0.93827); return; - } - - SetMass(0.13957); return; - -} - - - -/* - - -void AliTPCseed::CookdEdx2(Double_t low, Double_t up) { - //----------------------------------------------------------------- - // This funtion calculates dE/dX within the "low" and "up" cuts. - //----------------------------------------------------------------- - - Float_t amp[200]; - Float_t angular[200]; - Float_t weight[200]; - Int_t index[200]; - Bool_t inlimit[200]; - for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE; - for (Int_t i=0;i<200;i++) amp[i]=10000; - for (Int_t i=0;i<200;i++) angular[i]= 1;; - - - // - Float_t meanlog = 100.; - Int_t indexde[4]={0,64,128,160}; +// AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i) +// { +// // +// // +// return &fTrackPoints[i]; +// } - Float_t amean =0; - Float_t asigma =0; - Float_t anc =0; - Float_t anorm =0; - - Float_t mean[4] = {0,0,0,0}; - Float_t sigma[4] = {1000,1000,1000,1000}; - Int_t nc[4] = {0,0,0,0}; - Float_t norm[4] = {1000,1000,1000,1000}; - // - // - fNShared =0; - - // for (Int_t of =0; of<3; of++){ - // for (Int_t i=indexde[of];ifIsShared){ - fNShared++; - continue; - } - Int_t type = point->GetCPoint().GetType(); - if (type<0) continue; - if (point->GetCPoint().GetMax()<5) continue; - Float_t angley = point->GetTPoint().GetAngleY(); - Float_t anglez = point->GetTPoint().GetAngleZ(); - Float_t rsigmay = point->GetCPoint().GetSigmaY(); - Float_t rsigmaz = point->GetCPoint().GetSigmaZ(); - Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz); - - Float_t ampc = 0; // normalization to the number of electrons - if (i>64){ - ampc = point->GetCPoint().GetMax(); - } - else{ - ampc = point->GetCPoint().GetMax(); - } - ampc *= 2.0; // put mean value to channel 50 - // ampc *= 0.565; // put mean value to channel 50 - - Float_t w = 1.; - Float_t z = TMath::Abs(point->GetCPoint().GetZ()); - if (i<64) { - ampc /= 0.63; - } else - if (i>128){ - ampc /=1.51; - } - if (type<0) { //amp at the border - lower weight - continue; - } - if (rsigma>1.5) ampc/=1.3; // if big backround - angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez); - amp[i] = ampc/angular[i]; - weight[i] = w; - anc++; - } - TMath::Sort(159,amp,index,kFALSE); - for (Int_t i=int(anc*low+0.5);i0.1) - sigma[of] = TMath::Sqrt(sigma[of]); - else - sigma[of] = 1000; - mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog; - } - } - - Float_t dedx =0; - fSdEdx =0; - fMAngular =0; - // - Int_t norm2 = 0; - Int_t norm3 = 0; - Float_t www[3] = {12.,14.,17.}; - //Float_t www[3] = {1.,1.,1.}; - - for (Int_t i =0;i<3;i++){ - if (nc[i]>2&&nc[i]<1000){ - dedx += mean[i] *nc[i]*www[i]/sigma[i]; - fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i]; - fMAngular += norm[i] *nc[i]; - norm2 += nc[i]*www[i]/sigma[i]; - norm3 += (nc[i]-2)*www[i]/sigma[i]; - } - fDEDX[i] = mean[i]; - fSDEDX[i] = sigma[i]; - fNCDEDX[i]= nc[i]; - } - - if (norm3>0){ - dedx /=norm2; - fSdEdx /=norm3; - fMAngular/=norm2; - } - else{ - SetdEdx(0); - return; - } - // Float_t dedx1 =dedx; - - dedx =0; - Float_t norm4 = 0; - for (Int_t i =0;i<3;i++){ - if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){ - //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.)); - dedx += mean[i] *(nc[i])/(sigma[i]); - norm4 += (nc[i])/(sigma[i]); - } - fDEDX[i] = mean[i]; - } - if (norm4>0) dedx /= norm4; - - - - SetdEdx(dedx); - - //mi deDX - -} - -*/ diff --git a/TPC/AliTPCtrackerMI.h b/TPC/AliTPCtrackerMI.h index 211e941497e..53836a25209 100644 --- a/TPC/AliTPCtrackerMI.h +++ b/TPC/AliTPCtrackerMI.h @@ -13,9 +13,11 @@ // Origin: //------------------------------------------------------- +#include #include "AliTracker.h" -#include "AliTPCtrack.h" -#include "AliComplexCluster.h" +#include "AliTPCreco.h" +#include "AliPID.h" + class TFile; class AliTPCParam; @@ -25,87 +27,7 @@ class AliTPCTrackerPoint; class AliESD; class TTree; class AliESDkink; - -class AliTPCseed : public AliTPCtrack { - friend class AliTPCtrackerMI; - public: - AliTPCseed(); - virtual ~AliTPCseed(); - AliTPCseed(const AliTPCtrack &t); - AliTPCseed(const AliTPCseed &s); - //AliTPCseed(const AliTPCseed &t, Double_t a); - AliTPCseed(UInt_t index, const Double_t xx[5], - const Double_t cc[15], Double_t xr, Double_t alpha); - Int_t Compare(const TObject *o) const; - void Reset(Bool_t all = kTRUE); - Int_t GetProlongation(Double_t xr, Double_t &y, Double_t & z) const; - virtual Double_t GetPredictedChi2(const AliTPCclusterMI *cluster) const; - virtual Int_t Update(const AliTPCclusterMI* c, Double_t chi2, UInt_t i); - AliTPCTrackerPoint * GetTrackPoint(Int_t i); - void RebuildSeed(); // rebuild seed to be ready for storing - Double_t GetDensityFirst(Int_t n); - Double_t GetSigma2C() const {return fC44;}; - void GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2); - - void Modify(Double_t factor); - void SetClusterIndex2(Int_t row, Int_t index) { - fIndex[row] = index; - } - Int_t GetClusterIndex2(Int_t row) const { - return fIndex[row]; - } - Int_t GetClusterSector(Int_t row) const { - Int_t pica = -1; - if (fIndex[row]>=0) pica = ((fIndex[row]&0xff000000)>>24); - return pica; - } - - void SetErrorY2(Float_t sy2){fErrorY2=sy2;} - void SetErrorZ2(Float_t sz2){fErrorZ2=sz2;} - void CookdEdx(Double_t low=0.05, Double_t up=0.70, Int_t i1=0, Int_t i2=159, Bool_t onlyused = kFALSE); - // void CookdEdx2(Double_t low=0.05, Double_t up=0.70); - Bool_t IsActive() const { return !(fRemoval);} - void Desactivate(Int_t reason){ fRemoval = reason;} - // - // - private: - AliTPCseed & operator = (const AliTPCseed &); - AliESDtrack * fEsd; //! - AliTPCclusterMI* fClusterPointer[160]; //! array of cluster pointers - - TClonesArray * fPoints; // array with points along the track - TClonesArray * fEPoints; // array with exact points - calculated in special macro not used in tracking - //---CURRENT VALUES - Int_t fRow; //!current row number - Int_t fSector; //!current sector number - Int_t fRelativeSector; //! index of current relative sector - Float_t fCurrentSigmaY2; //!expected current cluster sigma Y - Float_t fCurrentSigmaZ2; //!expected current cluster sigma Z - Float_t fErrorY2; //!sigma of current cluster - Float_t fErrorZ2; //!sigma of current cluster - AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation - Int_t fCurrentClusterIndex1; //! index of the current cluster - Bool_t fInDead; //! indicate if the track is in dead zone - Bool_t fIsSeeding; //!indicates if it is proces of seeading - Int_t fNoCluster; //!indicates number of rows without clusters - Int_t fSort; //!indicate criteria for sorting - Bool_t fBSigned; //indicates that clusters of this trackes are signed to be used - // - // - Float_t fDEDX[4]; // dedx according padrows - Float_t fSDEDX[4]; // sdedx according padrows - Int_t fNCDEDX[4]; // number of clusters for dedx measurment - // - Int_t fSeedType; //seeding type - Int_t fSeed1; //first row for seeding - Int_t fSeed2; //last row for seeding - Int_t fOverlapLabels[12]; //track labels and the length of the overlap - Float_t fMAngular; // mean angular factor - AliTPCTrackerPoint fTrackPoints[160]; //!track points - array track points - - ClassDef(AliTPCseed,1) -}; - - +class TTreeSRedirector; class AliTPCtrackerMI : public AliTracker { @@ -134,6 +56,7 @@ public: void DeleteSeeds(); void SetDebug(Int_t debug){ fDebug = debug;} void FindKinks(TObjArray * array, AliESD * esd); + void FindV0s(TObjArray * array, AliESD * esd); void UpdateKinkQualityM(AliTPCseed * seed); void UpdateKinkQualityD(AliTPCseed * seed); Int_t CheckKinkPoint(AliTPCseed*seed, AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink); @@ -338,6 +261,7 @@ private: Double_t fYMax[200]; // max y for given pad row Double_t fPadLength[200]; // max y for given pad row const AliTPCParam *fParam; //pointer to the parameters + TTreeSRedirector *fDebugStreamer; //!debug streamer ClassDef(AliTPCtrackerMI,1) }; diff --git a/TPC/libTPCrec.pkg b/TPC/libTPCrec.pkg index a61fe3a3606..2bf4acb6217 100644 --- a/TPC/libTPCrec.pkg +++ b/TPC/libTPCrec.pkg @@ -4,7 +4,7 @@ SRCS:= AliTPCcluster.cxx \ AliClustersArray.cxx AliTPCClustersArray.cxx \ AliTPCclusterer.cxx AliTPCclustererMI.cxx \ AliTPCtrack.cxx AliTPCtracker.cxx \ - AliTPCpolyTrack.cxx AliTPCtrackerMI.cxx \ + AliTPCpolyTrack.cxx AliTPCseed.cxx AliTPCtrackerMI.cxx \ AliTPCPid.cxx AliTPCtrackPid.cxx AliTPCpidESD.cxx \ AliTPCReconstructor.cxx -- 2.43.0