X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCtrackerMI.cxx;h=4d5f0631fccb62922d7615bc6d5d7f61accafa9f;hp=91545affb09dbb7e1b46340678922f15164dabac;hb=93c541b7bdc7ed8c3f775c0de773420c52f3b108;hpb=92f513f5eea1f935fc8c1e9d7cd0e69c7fd9e308 diff --git a/TPC/AliTPCtrackerMI.cxx b/TPC/AliTPCtrackerMI.cxx index 91545affb09..4d5f0631fcc 100644 --- a/TPC/AliTPCtrackerMI.cxx +++ b/TPC/AliTPCtrackerMI.cxx @@ -67,7 +67,7 @@ // // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction // To enable storage of the TPC tracks in the ESD friend track -// use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0 +// use AliTPCReconstructor::SetStreamLevel(n); // // The debug level - different procedure produce tree for numerical debugging // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1 @@ -79,7 +79,7 @@ // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix // of the tracks (not to the clusters as they are dependent): // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError -// The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm) +// The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/GeV) // The default values are 0. // // The sytematic errors are added to the covariance matrix in following places: @@ -95,12 +95,14 @@ // There are several places in the code which can be numerically debuged // This code is keeped in order to enable code development and to check the calibration implementtion // -// 1. ErrParam stream (Log level 9) - dump information about +// 1. ErrParam stream - dump information about // 1.a) cluster // 2.a) cluster error estimate // 3.a) cluster shape estimate // // +// Debug streamer levels: +// //------------------------------------------------------- @@ -111,6 +113,7 @@ #include #include #include +#include #include "AliLog.h" #include "AliComplexCluster.h" #include "AliESDEvent.h" @@ -131,14 +134,18 @@ #include "AliTPCtrackerMI.h" #include "TStopwatch.h" #include "AliTPCReconstructor.h" -#include "AliPID.h" -#include "TTreeStream.h" #include "AliAlignObj.h" #include "AliTrackPointArray.h" #include "TRandom.h" #include "AliTPCcalibDB.h" +#include "AliTPCcalibDButil.h" #include "AliTPCTransform.h" #include "AliTPCClusterParam.h" +#include "AliTPCdEdxInfo.h" +#include "AliDCSSensorArray.h" +#include "AliDCSSensor.h" +#include "AliDAQ.h" +#include "AliCosmicTracker.h" // @@ -196,12 +203,23 @@ AliTPCtrackerMI::AliTPCtrackerMI() fNtracks(0), fSeeds(0), fIteration(0), - fParam(0), - fDebugStreamer(0) + fkParam(0), + fDebugStreamer(0), + fUseHLTClusters(4), + fSeedsPool(0), + fFreeSeedsID(500), + fNFreeSeeds(0), + fLastSeedID(-1) { // // default constructor // + for (Int_t irow=0; irow<200; irow++){ + fXRow[irow]=0; + fYMax[irow]=0; + fPadLength[irow]=0; + } + } //_____________________________________________________________________ @@ -213,8 +231,10 @@ Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){ AliTPCclusterMI* c =track->GetCurrentCluster(); - if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters - + if (accept > 0) //sign not accepted clusters + track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); + else // unsign accpeted clusters + track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff); UInt_t i = track->GetCurrentClusterIndex1(); Int_t sec=(i&0xff000000)>>24; @@ -222,7 +242,7 @@ Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){ track->SetRow((i&0x00ff0000)>>16); track->SetSector(sec); // Int_t index = i&0xFFFF; - if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow()); + if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow()); track->SetClusterIndex2(track->GetRow(), i); //track->fFirstPoint = row; //if ( track->fLastPointfLastPoint =row; @@ -289,37 +309,46 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste // // decide according desired precision to accept given // cluster for tracking + Double_t yt=0,zt=0; + seed->GetProlongation(cluster->GetX(),yt,zt); Double_t sy2=ErrY2(seed,cluster); Double_t sz2=ErrZ2(seed,cluster); Double_t sdistancey2 = sy2+seed->GetSigmaY2(); Double_t sdistancez2 = sz2+seed->GetSigmaZ2(); - - Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())* - (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2; - Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())* - (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2; + Double_t dy=seed->GetCurrentCluster()->GetY()-yt; + Double_t dz=seed->GetCurrentCluster()->GetZ()-zt; + Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)* + (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2; + Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)* + (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2; Double_t rdistance2 = rdistancey2+rdistancez2; //Int_t accept =0; - if (AliTPCReconstructor::StreamLevel()>5) { + if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) { Float_t rmsy2 = seed->GetCurrentSigmaY2(); Float_t rmsz2 = seed->GetCurrentSigmaZ2(); Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30(); Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30(); Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R(); Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R(); - - AliExternalTrackParam param(*seed); + AliExternalTrackParam param(*seed); static TVectorD gcl(3),gtr(3); Float_t gclf[3]; param.GetXYZ(gcl.GetMatrixArray()); cluster->GetGlobalXYZ(gclf); gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2]; + + + if (AliTPCReconstructor::StreamLevel()>2) { (*fDebugStreamer)<<"ErrParam"<< "Cl.="<16) return 3; + //return 0; // temporary + if (rdistance2>32) return 3; if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0) @@ -348,12 +382,23 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste seed->SetNFoundable(seed->GetNFoundable()-1); return 2; } + + if (fUseHLTClusters == 3 || fUseHLTClusters == 4) { + if (fIteration==2){ + if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) { + if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0) + return 2; + } + } + } + return 0; } + //_____________________________________________________________________________ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par): AliTracker(), @@ -373,8 +418,13 @@ AliTracker(), fNtracks(0), fSeeds(0), fIteration(0), - fParam(0), - fDebugStreamer(0) + fkParam(0), + fDebugStreamer(0), + fUseHLTClusters(4), + fSeedsPool(0), + fFreeSeedsID(500), + fNFreeSeeds(0), + fLastSeedID(-1) { //--------------------------------------------------------------------- // The main TPC tracker constructor @@ -386,25 +436,29 @@ AliTracker(), for (i=0; iGetNRowLow(); Int_t nrowup = par->GetNRowUp(); - for (Int_t i=0;iGetPadRowRadiiLow(i); fPadLength[i]= par->GetPadPitchLength(0,i); fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle()); } - for (Int_t i=0;iGetPadRowRadiiUp(i); fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i); fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle()); } - fDebugStreamer = new TTreeSRedirector("TPCdebug.root"); + if (AliTPCReconstructor::StreamLevel()>0) { + fDebugStreamer = new TTreeSRedirector("TPCdebug.root"); + } + // + fSeedsPool = new TClonesArray("AliTPCseed",1000); } //________________________________________________________________________ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): @@ -425,15 +479,27 @@ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): fNtracks(0), fSeeds(0), fIteration(0), - fParam(0), - fDebugStreamer(0) + fkParam(0), + fDebugStreamer(0), + fUseHLTClusters(4), + fSeedsPool(0), + fFreeSeedsID(500), + fNFreeSeeds(0), + fLastSeedID(-1) { //------------------------------------ // dummy copy constructor //------------------------------------------------------------------ fOutput=t.fOutput; + for (Int_t irow=0; irow<200; irow++){ + fXRow[irow]=0; + fYMax[irow]=0; + fPadLength[irow]=0; + } + } -AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){ +AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/) +{ //------------------------------ // dummy //-------------------------------------------------------------- @@ -447,19 +513,21 @@ AliTPCtrackerMI::~AliTPCtrackerMI() { delete[] fInnerSec; delete[] fOuterSec; if (fSeeds) { - fSeeds->Delete(); + fSeeds->Clear(); delete fSeeds; } if (fDebugStreamer) delete fDebugStreamer; + if (fSeedsPool) delete fSeedsPool; } -void AliTPCtrackerMI::FillESD(TObjArray* arr) +void AliTPCtrackerMI::FillESD(const TObjArray* arr) { // // //fill esds using updated tracks - if (fEvent){ + if (!fEvent) return; + // write tracks to the event // store index of the track Int_t nseed=arr->GetEntriesFast(); @@ -469,9 +537,14 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) if (!pt) continue; pt->UpdatePoints(); AddCovariance(pt); - // pt->PropagateTo(fParam->GetInnerRadiusLow()); + if (AliTPCReconstructor::StreamLevel()>1) { + (*fDebugStreamer)<<"Track0"<< + "Tr0.="<PropagateTo(fkParam->GetInnerRadiusLow()); if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks - pt->PropagateTo(fParam->GetInnerRadiusLow()); + pt->PropagateTo(fkParam->GetInnerRadiusLow()); } if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){ @@ -569,22 +642,39 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) continue; } } - } - printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()); + // >> account for suppressed tracks in the kink indices (RS) + int nESDtracks = fEvent->GetNumberOfTracks(); + for (int it=nESDtracks;it--;) { + AliESDtrack* esdTr = fEvent->GetTrack(it); + if (!esdTr || !esdTr->GetKinkIndex(0)) continue; + for (int ik=0;ik<3;ik++) { + int knkId=0; + if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track + AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1); + if (!kink) { + AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1)); + continue; + } + kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1 + } + } + // << account for suppressed tracks in the kink indices (RS) + AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks())); + } -Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ +Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){ // // // Use calibrated cluster error from OCDB // AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam(); // - Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ())); + Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())); Int_t ctype = cl->GetType(); Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2; Double_t angle = seed->GetSnp()*seed->GetSnp(); @@ -632,7 +722,7 @@ Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ // gnoise1 = 0.0004/padlength; // Float_t nel = 0.268*amp; // Float_t nprim = 0.155*amp; -// ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel; +// ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel; // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim; // if (glandau1[i]>1) glandau1[i]=1; // glandau1[i]*=padlength*padlength/12.; @@ -642,7 +732,7 @@ Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ // gnoise2 = 0.0004/padlength; // nel = 0.3*amp; // nprim = 0.133*amp; -// ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; +// ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; // if (glandau2[i]>1) glandau2[i]=1; // glandau2[i]*=padlength*padlength/12.; @@ -653,7 +743,7 @@ Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ // gnoise3 = 0.0004/padlength; // nel = 0.3*amp; // nprim = 0.133*amp; -// ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; +// ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; // if (glandau3[i]>1) glandau3[i]=1; // glandau3[i]*=padlength*padlength/12.; @@ -670,7 +760,7 @@ Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ // return 1.; // } // Float_t snoise2; -// Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ())); +// Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())); // Int_t ctype = cl->GetType(); // Float_t padlength= GetPadPitchLength(seed->GetRow()); // Double_t angle2 = seed->GetSnp()*seed->GetSnp(); @@ -727,14 +817,14 @@ Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){ -Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ +Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){ // // // Use calibrated cluster error from OCDB // AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam(); // - Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ())); + Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())); Int_t ctype = cl->GetType(); Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2; // @@ -788,7 +878,7 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ // gnoise1 = 0.0004/padlength; // Float_t nel = 0.268*amp; // Float_t nprim = 0.155*amp; -// ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel; +// ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel; // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim; // if (glandau1[i]>1) glandau1[i]=1; // glandau1[i]*=padlength*padlength/12.; @@ -798,7 +888,7 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ // gnoise2 = 0.0004/padlength; // nel = 0.3*amp; // nprim = 0.133*amp; -// ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; +// ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; // if (glandau2[i]>1) glandau2[i]=1; // glandau2[i]*=padlength*padlength/12.; @@ -809,7 +899,7 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ // gnoise3 = 0.0004/padlength; // nel = 0.3*amp; // nprim = 0.133*amp; -// ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; +// ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel; // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; // if (glandau3[i]>1) glandau3[i]=1; // glandau3[i]*=padlength*padlength/12.; @@ -826,7 +916,7 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ // return 1.; // } // Float_t snoise2; -// Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ())); +// Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())); // Int_t ctype = cl->GetType(); // Float_t padlength= GetPadPitchLength(seed->GetRow()); // // @@ -912,7 +1002,7 @@ void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed) //_____________________________________________________________________________ Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1, Double_t x2,Double_t y2, - Double_t x3,Double_t y3) + Double_t x3,Double_t y3) const { //----------------------------------------------------------------- // Initial approximation of the track curvature @@ -933,7 +1023,7 @@ Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1, //_____________________________________________________________________________ Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1, Double_t x2,Double_t y2, - Double_t x3,Double_t y3) + Double_t x3,Double_t y3) const { //----------------------------------------------------------------- // Initial approximation of the track curvature @@ -944,7 +1034,7 @@ Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1, y2 -=y1; // Double_t det = x3*y2-x2*y3; - if (det==0) { + if (TMath::Abs(det)<1e-10){ return 100; } // @@ -959,7 +1049,7 @@ Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1, Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1, Double_t x2,Double_t y2, - Double_t x3,Double_t y3) + Double_t x3,Double_t y3) const { //----------------------------------------------------------------- // Initial approximation of the track curvature @@ -970,7 +1060,7 @@ Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1, y2 -=y1; // Double_t det = x3*y2-x2*y3; - if (det==0) { + if (TMath::Abs(det)<1e-10) { return 100; } // @@ -989,7 +1079,7 @@ Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1, //_____________________________________________________________________________ Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1, Double_t x2,Double_t y2, - Double_t x3,Double_t y3) + Double_t x3,Double_t y3) const { //----------------------------------------------------------------- // Initial approximation of the track curvature times center of curvature @@ -1008,7 +1098,7 @@ Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1, //_____________________________________________________________________________ Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1, Double_t x2,Double_t y2, - Double_t z1,Double_t z2) + Double_t z1,Double_t z2) const { //----------------------------------------------------------------- // Initial approximation of the tangent of the track dip angle @@ -1019,7 +1109,7 @@ Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1, Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1, Double_t x2,Double_t y2, - Double_t z1,Double_t z2, Double_t c) + Double_t z1,Double_t z2, Double_t c) const { //----------------------------------------------------------------- // Initial approximation of the tangent of the track dip angle @@ -1039,7 +1129,7 @@ Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1, return angle2; } -Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) +Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const {//----------------------------------------------------------------- // This function find proloncation of a track to a reference plane x=x2. //----------------------------------------------------------------- @@ -1050,8 +1140,8 @@ Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5 return kFALSE; } - Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1); - Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2); + Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1)); + Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2)); y = x[0]; z = x[1]; @@ -1074,32 +1164,32 @@ Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5 return kTRUE; } -Int_t AliTPCtrackerMI::LoadClusters (TTree *tree) +Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree) { - // + // load clusters // fInput = tree; return LoadClusters(); } -Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr) +Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr) { // // load clusters to the memory - AliTPCClustersRow *clrow = 0x0; + AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI"); Int_t lower = arr->LowerBound(); Int_t entries = arr->GetEntriesFast(); - + for (Int_t i=lower; iAt(i); - if(!clrow) continue; - if(!clrow->GetArray()) continue; - + if(!clrow) continue; + if(!clrow->GetArray()) continue; + // Int_t sec,row; - fParam->AdjustSectorRow(clrow->GetID(),sec,row); - + fkParam->AdjustSectorRow(clrow->GetID(),sec,row); + for (Int_t icl=0; iclGetArray()->GetEntriesFast(); icl++){ Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl))); } @@ -1117,16 +1207,15 @@ Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr) } if (left ==0){ tpcrow->SetN1(clrow->GetArray()->GetEntriesFast()); - tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]); - for (Int_t i=0;iGetN1();i++) - tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i))); + for (Int_t j=0;jGetN1();++j) + tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j))); } if (left ==1){ tpcrow->SetN2(clrow->GetArray()->GetEntriesFast()); - tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]); - for (Int_t i=0;iGetN2();i++) - tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i))); + for (Int_t j=0;jGetN2();++j) + tpcrow->SetCluster2(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j))); } + clrow->GetArray()->Clear("C"); } // delete clrow; @@ -1135,28 +1224,85 @@ Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr) return 0; } +Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr) +{ + // + // load clusters to the memory from one + // TClonesArray + // + AliTPCclusterMI *clust=0; + Int_t count[72][96] = { {0} , {0} }; + + // loop over clusters + for (Int_t icl=0; iclGetEntriesFast(); icl++) { + clust = (AliTPCclusterMI*)arr->At(icl); + if(!clust) continue; + //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow()); + + // transform clusters + Transform(clust); + + // count clusters per pad row + count[clust->GetDetector()][clust->GetRow()]++; + } + + // insert clusters to sectors + for (Int_t icl=0; iclGetEntriesFast(); icl++) { + clust = (AliTPCclusterMI*)arr->At(icl); + if(!clust) continue; + + Int_t sec = clust->GetDetector(); + Int_t row = clust->GetRow(); + + // filter overlapping pad rows needed by HLT + if(secSetClass("AliTPCclusterMI"); - clrow->SetArray(0); - clrow->GetArray()->ExpandCreateFast(10000); + static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI"); // // TTree * tree = fClustersArray.GetTree(); TTree * tree = fInput; TBranch * br = tree->GetBranch("Segment"); br->SetAddress(&clrow); - // + + // Conversion of pad, row coordinates in local tracking coords. + // Could be skipped here; is already done in clusterfinder + Int_t j=Int_t(tree->GetEntries()); for (Int_t i=0; iGetEntry(i); // Int_t sec,row; - fParam->AdjustSectorRow(clrow->GetID(),sec,row); + fkParam->AdjustSectorRow(clrow->GetID(),sec,row); for (Int_t icl=0; iclGetArray()->GetEntriesFast(); icl++){ Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl))); } @@ -1173,19 +1319,17 @@ Int_t AliTPCtrackerMI::LoadClusters() } if (left ==0){ tpcrow->SetN1(clrow->GetArray()->GetEntriesFast()); - tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]); - for (Int_t i=0;iGetN1();i++) - tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i))); + for (Int_t k=0;kGetN1();++k) + tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k))); } if (left ==1){ tpcrow->SetN2(clrow->GetArray()->GetEntriesFast()); - tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]); - for (Int_t i=0;iGetN2();i++) - tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i))); + for (Int_t k=0;kGetN2();++k) + tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k))); } } // - delete clrow; + clrow->Clear("C"); LoadOuterSectors(); LoadInnerSectors(); return 0; @@ -1250,25 +1394,26 @@ void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{ void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){ // + // transformation // - // - AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance(); + AliTPCTransform *transform = calibDB->GetTransform() ; if (!transform) { AliFatal("Tranformations not in calibDB"); + return; } + transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam()); Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()}; Int_t i[1]={cluster->GetDetector()}; transform->Transform(x,i,0,1); - if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){ - if (cluster->GetDetector()%36>17){ - x[1]*=-1; - } - } + // if (cluster->GetDetector()%36>17){ + // x[1]*=-1; + //} // // in debug mode check the transformation // - if (AliTPCReconstructor::StreamLevel()>1) { + if (AliTPCReconstructor::StreamLevel()>2) { Float_t gx[3]; cluster->GetGlobalXYZ(gx); Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile(); @@ -1288,22 +1433,23 @@ void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){ cluster->SetY(x[1]); cluster->SetZ(x[2]); // The old stuff: - // // // - //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices(); - TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector()); - //TGeoHMatrix mat; - Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()}; - Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; - if (mat) mat->LocalToMaster(pos,posC); - else{ - // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX - } - cluster->SetX(posC[0]); - cluster->SetY(posC[1]); - cluster->SetZ(posC[2]); + //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices(); + if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){ + TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector()); + //TGeoHMatrix mat; + Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()}; + Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()}; + if (mat) mat->LocalToMaster(pos,posC); + else{ + // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX + } + cluster->SetX(posC[0]); + cluster->SetY(posC[1]); + cluster->SetZ(posC[2]); + } } //_____________________________________________________________________________ @@ -1416,7 +1562,7 @@ AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const { Int_t ncl=(index&0x00007fff)>>00; const AliTPCtrackerRow * tpcrow=0; - AliTPCclusterMI * clrow =0; + TClonesArray * clrow =0; if (sec<0 || sec>=fkNIS*4) { AliWarning(Form("Wrong sector %d",sec)); @@ -1424,7 +1570,9 @@ AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const { } if (secAt(ncl); } @@ -1462,6 +1612,8 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { //----------------------------------------------------------------- // Double_t x= GetXrow(nr), ymax=GetMaxY(nr); + // + // AliTPCclusterMI *cl=0; Int_t tpcindex= t.GetClusterIndex2(nr); // @@ -1473,9 +1625,10 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { if (fIteration>0 && tpcindex>=-1){ //if we have already clusters // if (tpcindex==-1) return 0; //track in dead zone - if (tpcindex>0){ // + if (tpcindex >= 0){ // cl = t.GetClusterPointer(nr); - if ( (cl==0) ) cl = GetClusterMI(tpcindex); + //if (cl==0) cl = GetClusterMI(tpcindex); + if (!cl) cl = GetClusterMI(tpcindex); t.SetCurrentClusterIndex1(tpcindex); } if (cl){ @@ -1488,9 +1641,15 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { if (TMath::Abs(angle-t.GetAlpha())>0.001){ Double_t rotation = angle-t.GetAlpha(); t.SetRelativeSector(relativesector); - if (!t.Rotate(rotation)) return 0; + if (!t.Rotate(rotation)) { + t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000); + return 0; + } + } + if (!t.PropagateTo(x)) { + t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000); + return 0; } - if (!t.PropagateTo(x)) return 0; // t.SetCurrentCluster(cl); t.SetRow(nr); @@ -1507,11 +1666,15 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { t.SetNFoundable(t.GetNFoundable()+1); UpdateTrack(&t,accept); return 1; - } + } + else { // Remove old cluster from track + t.SetClusterIndex(nr, -3); + t.SetClusterPointer(nr, 0); + } } } if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle - if (fIteration>1){ + if (fIteration>1 && IsFindable(t)){ // not look for new cluster during refitting t.SetNFoundable(t.GetNFoundable()+1); return 0; @@ -1519,7 +1682,11 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { // UInt_t index=0; // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06 - Double_t y=t.GetYat(x); + if (!t.PropagateTo(x)) { + if (fIteration==0) t.SetRemoval(10); + return 0; + } + Double_t y = t.GetY(); if (TMath::Abs(y)>ymax){ if (y > ymax) { t.SetRelativeSector((t.GetRelativeSector()+1) % fN); @@ -1530,14 +1697,13 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { if (!t.Rotate(-fSectors->GetAlpha())) return 0; } - //return 1; + if (!t.PropagateTo(x)) { + if (fIteration==0) t.SetRemoval(10); + return 0; + } + y = t.GetY(); } // - if (!t.PropagateTo(x)) { - if (fIteration==0) t.SetRemoval(10); - return 0; - } - y=t.GetY(); Double_t z=t.GetZ(); // @@ -1564,7 +1730,8 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { } else { - if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)GetZLength(0) && (TMath::Abs(t.GetSnp())GetZLength(0) && (TMath::Abs(t.GetSnp())ymax){ - - if (y > ymax) { - t.SetRelativeSector((t.GetRelativeSector()+1) % fN); - if (!t.Rotate(fSectors->GetAlpha())) - return 0; - } else if (y <-ymax) { - t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN); - if (!t.Rotate(-fSectors->GetAlpha())) - return 0; - } - if (!t.PropagateTo(x)) { - return 0; - } - t.GetProlongation(x,y,z); - } - // - // update current shape info every 2 pad-row - if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){ - // t.fCurrentSigmaY = GetSigmaY(&t); - //t.fCurrentSigmaZ = GetSigmaZ(&t); - GetShape(&t,nr); - } - // - AliTPCclusterMI *cl=0; - UInt_t index=0; - - - //Int_t nr2 = nr; - const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr); - if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0; - Double_t roady =1.; - Double_t roadz = 1.; - // - Int_t row = nr; - if (TMath::Abs(TMath::Abs(y)-ymax)(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1); - } - //calculate - - if ((cl==0)&&(krow)) { - // cl = krow.FindNearest2(y+10,z,roady,roadz,index); - cl = krow.FindNearest2(y,z,roady,roadz,index); - - if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index)); - } - - if (cl) { - t.SetCurrentCluster(cl); - // Int_t accept = AcceptCluster(&t,t.fCurrentCluster); - //if (accept<3){ - t.SetClusterIndex2(row,index); - t.SetClusterPointer(row, cl); - //} - } - return 1; -} - //_________________________________________________________________________ @@ -1696,19 +1786,19 @@ Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const Int_t sector = (index&0xff000000)>>24; // Int_t row = (index&0x00ff0000)>>16; Float_t xyz[3]; - // xyz[0] = fParam->GetPadRowRadii(sector,row); + // xyz[0] = fkParam->GetPadRowRadii(sector,row); xyz[0] = cl->GetX(); xyz[1] = cl->GetY(); xyz[2] = cl->GetZ(); Float_t sin,cos; - fParam->AdjustCosSin(sector,cos,sin); + fkParam->AdjustCosSin(sector,cos,sin); Float_t x = cos*xyz[0]-sin*xyz[1]; Float_t y = cos*xyz[1]+sin*xyz[0]; Float_t cov[6]; Float_t sigmaY2 = 0.027*cl->GetSigmaY2(); - if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07; + if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07; Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2(); - if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77; + if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77; cov[0] = sin*sin*sigmaY2; cov[1] = -sin*cos*sigmaY2; cov[2] = 0.; @@ -1718,13 +1808,13 @@ Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const p.SetXYZ(x,y,xyz[2],cov); AliGeomManager::ELayerID iLayer; Int_t idet; - if (sector < fParam->GetNInnerSector()) { + if (sector < fkParam->GetNInnerSector()) { iLayer = AliGeomManager::kTPC1; idet = sector; } else { iLayer = AliGeomManager::kTPC2; - idet = sector - fParam->GetNInnerSector(); + idet = sector - fkParam->GetNInnerSector(); } UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet); p.SetVolumeID(volid); @@ -1738,7 +1828,7 @@ Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) { // This function tries to find a track prolongation to next pad row //----------------------------------------------------------------- t.SetCurrentCluster(0); - t.SetCurrentClusterIndex1(0); + t.SetCurrentClusterIndex1(-3); Double_t xt=t.GetX(); Int_t row = GetRowNumber(xt)-1; @@ -1800,7 +1890,9 @@ Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) { } else { - if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())0) && (index&0x8000)==0){ + if ( (index >= 0) && (index&0x8000)==0){ cl = t.GetClusterPointer(nr); - if ( (cl==0) && (index>0)) cl = GetClusterMI(index); + if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index); t.SetCurrentClusterIndex1(index); if (cl) { t.SetCurrentCluster(cl); @@ -1874,8 +1966,8 @@ Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) { } else { if (fIteration==0){ - if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10); - if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10); + if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10); + if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10); if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10); } @@ -1886,19 +1978,23 @@ Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) { //_____________________________________________________________________________ -Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) { +Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) { //----------------------------------------------------------------- // This function tries to find a track prolongation. //----------------------------------------------------------------- Double_t xt=t.GetX(); // - Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift(); + Double_t alpha=t.GetAlpha(); if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); if (alpha < 0. ) alpha += 2.*TMath::Pi(); // t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN); - Int_t first = GetRowNumber(xt)-1; + Int_t first = GetRowNumber(xt); + if (!fromSeeds) + first -= step; + if (first < 0) + first = 0; for (Int_t nr= first; nr>=rf; nr-=step) { // update kink info if (t.GetKinkIndexes()[0]>0){ @@ -1937,44 +2033,27 @@ Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) { } -//_____________________________________________________________________________ -Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) { - //----------------------------------------------------------------- - // This function tries to find a track prolongation. - //----------------------------------------------------------------- - Double_t xt=t.GetX(); - // - Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift(); - if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); - if (alpha < 0. ) alpha += 2.*TMath::Pi(); - t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN); - - for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) { - - if (FollowToNextFast(t,nr)==0) - if (!t.IsActive()) return 0; - - } - return 1; -} - -Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) { +Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) { //----------------------------------------------------------------- // This function tries to find a track prolongation. //----------------------------------------------------------------- // Double_t xt=t.GetX(); - Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift(); + Double_t alpha=t.GetAlpha(); if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); if (alpha < 0. ) alpha += 2.*TMath::Pi(); t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN); Int_t first = t.GetFirstPoint(); - if (firstGetSigmaTglZ()+ s2->GetSigmaTglZ()); if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return; @@ -2081,11 +2160,7 @@ void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2) // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) { if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) { sumshared++; - s1->SetSharedMapBit(i, kTRUE); - s2->SetSharedMapBit(i, kTRUE); } - if (s1->GetClusterIndex2(i)>0) - s1->SetClusterMapBit(i, kTRUE); } if (sumshared>cutN0){ // sign clusters @@ -2143,7 +2218,7 @@ void AliTPCtrackerMI::SignShared(TObjArray * arr) for (Int_t i=0; iUncheckedAt(i); if (!pt) continue; - for (Int_t j=0;j<=12;j++){ + for (Int_t j=0;j<12;j++){ pt->SetOverlapLabel(j,0); } } @@ -2199,8 +2274,11 @@ void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t fac // - AliTPCtrackerMI::PropagateBack() // - AliTPCtrackerMI::RefitInward() // - - + // Arguments: + // factor1 - factor for constrained + // factor2 - for non constrained tracks + // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE + // UnsignClusters(); // Int_t nseed = arr->GetEntriesFast(); @@ -2231,12 +2309,7 @@ void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t fac if (!pt) continue; // if (quality[trackindex]<0){ - if (pt) { - delete arr->RemoveAt(trackindex); - } - else{ - arr->RemoveAt(trackindex); - } + MarkSeedFree( arr->RemoveAt(trackindex) ); continue; } // @@ -2257,12 +2330,26 @@ void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t fac // if (Float_t(shared+1)/Float_t(found+1)>factor){ if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks - delete arr->RemoveAt(trackindex); + if( AliTPCReconstructor::StreamLevel()>3){ + TTreeSRedirector &cstream = *fDebugStreamer; + cstream<<"RemoveUsed"<< + "iter="<RemoveAt(trackindex) ); continue; } if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks - delete arr->RemoveAt(trackindex); + if( AliTPCReconstructor::StreamLevel()>3){ + TTreeSRedirector &cstream = *fDebugStreamer; + cstream<<"RemoveShort"<< + "iter="<RemoveAt(trackindex) ); continue; } } @@ -2290,6 +2377,94 @@ void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t fac delete []indexes; } +void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray) +{ + // + // Dump clusters after reco + // signed and unsigned cluster can be visualized + // 1. Unsign all cluster + // 2. Sign all used clusters + // 3. Dump clusters + UnsignClusters(); + Int_t nseed = trackArray->GetEntries(); + for (Int_t i=0; iUncheckedAt(i); + if (!pt) { + continue; + } + Bool_t isKink=pt->GetKinkIndex(0)!=0; + for (Int_t j=0; j<160; ++j) { + Int_t index=pt->GetClusterIndex2(j); + if (index<0) continue; + AliTPCclusterMI *c= pt->GetClusterPointer(j); + if (!c) continue; + if (isKink) c->Use(100); // kink + c->Use(10); // by default usage 10 + } + } + // + + for (Int_t sec=0;secGetNRows();row++){ + TClonesArray *cla = fInnerSec[sec][row].GetClusters1(); + for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){ + AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl); + Float_t gx[3]; cl->GetGlobalXYZ(gx); + (*fDebugStreamer)<<"clDump"<< + "iter="<At(icl); + Float_t gx[3]; cl->GetGlobalXYZ(gx); + (*fDebugStreamer)<<"clDump"<< + "iter="<GetNRows();row++){ + TClonesArray *cla = fOuterSec[sec][row].GetClusters1(); + for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){ + Float_t gx[3]; + AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl); + cl->GetGlobalXYZ(gx); + (*fDebugStreamer)<<"clDump"<< + "iter="<At(icl); + cl->GetGlobalXYZ(gx); + (*fDebugStreamer)<<"clDump"<< + "iter="<GetNRows();row++){ - AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1(); + TClonesArray *cla = fInnerSec[sec][row].GetClusters1(); for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++) // if (cl[icl].IsUsed(10)) - cl[icl].Use(-1); - cl = fInnerSec[sec][row].GetClusters2(); + ((AliTPCclusterMI*) cla->At(icl))->Use(-1); + cla = fInnerSec[sec][row].GetClusters2(); for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++) //if (cl[icl].IsUsed(10)) - cl[icl].Use(-1); + ((AliTPCclusterMI*) cla->At(icl))->Use(-1); } } for (Int_t sec=0;secGetNRows();row++){ - AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1(); + TClonesArray *cla = fOuterSec[sec][row].GetClusters1(); for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++) //if (cl[icl].IsUsed(10)) - cl[icl].Use(-1); - cl = fOuterSec[sec][row].GetClusters2(); + ((AliTPCclusterMI*) cla->At(icl))->Use(-1); + cla = fOuterSec[sec][row].GetClusters2(); for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++) //if (cl[icl].IsUsed(10)) - cl[icl].Use(-1); + ((AliTPCclusterMI*) cla->At(icl))->Use(-1); } } @@ -2326,7 +2501,7 @@ void AliTPCtrackerMI::UnsignClusters() -void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity) +void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity) { // //sign clusters to be "used" @@ -2415,7 +2590,7 @@ void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fde pt->GetClusterStatistic(0,160,found, foundable,shared); if (shared/float(found)>0.3) { if (shared/float(found)>0.9 ){ - //delete arr->RemoveAt(i); + //MarkSeedFree( arr->RemoveAt(i) ); } continue; } @@ -2431,10 +2606,10 @@ void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fde isok =kTRUE; if (isok) - for (Int_t i=0; i<160; i++) { - Int_t index=pt->GetClusterIndex2(i); + for (Int_t j=0; j<160; ++j) { + Int_t index=pt->GetClusterIndex2(j); if (index<0) continue; - AliTPCclusterMI *c= pt->GetClusterPointer(i); + AliTPCclusterMI *c= pt->GetClusterPointer(j); if (!c) continue; //if (!(c->IsUsed(10))) c->Use(); c->Use(10); @@ -2473,11 +2648,11 @@ void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fde if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chiGetNumberOfClusters(); pt->SetBSigned(kTRUE); - for (Int_t i=0; i<160; i++) { + for (Int_t j=0; j<160; ++j) { - Int_t index=pt->GetClusterIndex2(i); + Int_t index=pt->GetClusterIndex2(j); if (index<0) continue; - AliTPCclusterMI *c= pt->GetClusterPointer(i); + AliTPCclusterMI *c= pt->GetClusterPointer(j); if (!c) continue; // if (!(c->IsUsed(10))) c->Use(); c->Use(10); @@ -2492,59 +2667,6 @@ void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fde } -void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const -{ - // stop not active tracks - // take th1 as threshold for number of founded to number of foundable on last 10 active rows - // take th2 as threshold for number of founded to number of foundable on last 20 active rows - Int_t nseed = arr->GetEntriesFast(); - // - for (Int_t i=0; iUncheckedAt(i); - if (!pt) { - continue; - } - if (!(pt->IsActive())) continue; - StopNotActive(pt,row0,th0, th1,th2); - } -} - - - -void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1, - Float_t th2) const -{ - // stop not active tracks - // take th1 as threshold for number of founded to number of foundable on last 10 active rows - // take th2 as threshold for number of founded to number of foundable on last 20 active rows - Int_t sumgood1 = 0; - Int_t sumgood2 = 0; - Int_t foundable = 0; - Int_t maxindex = seed->GetLastPoint(); //last foundable row - if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) { - seed->Desactivate(10) ; - return; - } - - for (Int_t i=row0; iGetClusterIndex2(i); - if (index!=-1) foundable++; - //if (!c) continue; - if (foundable<=30) sumgood1++; - if (foundable<=50) { - sumgood2++; - } - else{ - break; - } - } - if (foundable>=30.){ - if (sumgood1<(th1*30.)) seed->Desactivate(10); - } - if (foundable>=50) - if (sumgood2<(th2*50.)) seed->Desactivate(10); -} - Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event) { @@ -2552,8 +2674,29 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event) // back propagation of ESD tracks // //return 0; + if (!event) return 0; const Int_t kMaxFriendTracks=2000; fEvent = event; + // extract correction object for multiplicity dependence of dEdx + TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber()); + + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + if (!transform) { + AliFatal("Tranformations not in RefitInward"); + return 0; + } + transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam()); + const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam(); + Int_t nContribut = event->GetNumberOfTracks(); + TGraphErrors * graphMultDependenceDeDx = 0x0; + if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) { + if (recoParam->GetUseTotCharge()) { + graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL"); + } else { + graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL"); + } + } + // ReadSeeds(event,2); fIteration=2; //PrepareForProlongation(fSeeds,1); @@ -2567,6 +2710,7 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event) SignShared(&arraySeed); // FindCurling(fSeeds, event,2); // find multi found tracks FindSplitted(fSeeds, event,2); // find multi found tracks + if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks Int_t ntracks=0; Int_t nseed = fSeeds->GetEntriesFast(); @@ -2574,18 +2718,35 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event) AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i); if (!seed) continue; if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks + AliESDtrack *esd=event->GetTrack(i); + + if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){ + AliExternalTrackParam paramIn; + AliExternalTrackParam paramOut; + Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut); + if (AliTPCReconstructor::StreamLevel()>2) { + (*fDebugStreamer)<<"RecoverIn"<< + "seed.="<15) { + seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance()); + seed->SetNumberOfClusters(ncl); + } + } - seed->PropagateTo(fParam->GetInnerRadiusLow()); + seed->PropagateTo(fkParam->GetInnerRadiusLow()); seed->UpdatePoints(); AddCovariance(seed); - MakeBitmaps(seed); - AliESDtrack *esd=event->GetTrack(i); + MakeESDBitmaps(seed, esd); seed->CookdEdx(0.02,0.6); CookLabel(seed,0.1); //For comparison only - esd->SetTPCClusterMap(seed->GetClusterMap()); - esd->SetTPCSharedMap(seed->GetSharedMap()); // - if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) { + if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) { TTreeSRedirector &cstream = *fDebugStreamer; cstream<<"Crefit"<< "Esd.="<UpdateTrackParams(seed,AliESDtrack::kTPCrefit); esd->SetTPCPoints(seed->GetPoints()); esd->SetTPCPointsF(seed->GetNFoundable()); - Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3); - Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25; + Int_t ndedx = seed->GetNCDEDX(0); + Float_t sdedx = seed->GetSDEDX(0); Float_t dedx = seed->GetdEdx(); + // apply mutliplicity dependent dEdx correction if available + if (graphMultDependenceDeDx) { + Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut); + dedx += (1 - corrGain)*50.; // MIP is normalized to 50 + } esd->SetTPCsignal(dedx, sdedx, ndedx); // + // fill new dEdx information + // + Double32_t signal[4]; + Char_t ncl[3]; + Char_t nrows[3]; + // + for(Int_t iarr=0;iarr<3;iarr++) { + signal[iarr] = seed->GetDEDXregion(iarr+1); + ncl[iarr] = seed->GetNCDEDX(iarr+1); + nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1); + } + signal[3] = seed->GetDEDXregion(4); + // + AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo(); + infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows); + esd->SetTPCdEdxInfo(infoTpcPid); + // // add seed to the esd track in Calib level // Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed); if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){ + // RS: this is the only place where the seed is created not in the pool, + // since it should belong to ESDevent AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE); esd->AddCalibObject(seedCopy); } @@ -2616,7 +2801,11 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event) } } //FindKinks(fSeeds,event); + if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds); Info("RefitInward","Number of refitted tracks %d",ntracks); + + AliCosmicTracker::FindCosmic(event, kTRUE); + return 0; } @@ -2626,7 +2815,7 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event) // // back propagation of ESD tracks // - + if (!event) return 0; fEvent = event; fIteration = 1; ReadSeeds(event,1); @@ -2634,6 +2823,7 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event) RemoveUsed2(fSeeds,0.4,0.4,20); //FindCurling(fSeeds, fEvent,1); FindSplitted(fSeeds, event,1); // find multi found tracks + if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks // Int_t nseed = fSeeds->GetEntriesFast(); @@ -2645,27 +2835,48 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event) seed->UpdatePoints(); AddCovariance(seed); AliESDtrack *esd=event->GetTrack(i); + if (!esd) continue; //never happen + if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){ + AliExternalTrackParam paramIn; + AliExternalTrackParam paramOut; + Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut); + if (AliTPCReconstructor::StreamLevel()>2) { + (*fDebugStreamer)<<"RecoverBack"<< + "seed.="<15) { + seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance()); + seed->SetNumberOfClusters(ncl); + } + } seed->CookdEdx(0.02,0.6); CookLabel(seed,0.1); //For comparison only if (seed->GetNumberOfClusters()>15){ esd->UpdateTrackParams(seed,AliESDtrack::kTPCout); esd->SetTPCPoints(seed->GetPoints()); esd->SetTPCPointsF(seed->GetNFoundable()); - Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3); - Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25; + Int_t ndedx = seed->GetNCDEDX(0); + Float_t sdedx = seed->GetSDEDX(0); Float_t dedx = seed->GetdEdx(); esd->SetTPCsignal(dedx, sdedx, ndedx); ntracks++; Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number - if (AliTPCReconstructor::StreamLevel()>1) { + if (AliTPCReconstructor::StreamLevel()>1 && esd) { (*fDebugStreamer)<<"Cback"<< "Tr0.="<3) DumpClusters(1,fSeeds); //FindKinks(fSeeds,event); Info("PropagateBack","Number of back propagated tracks %d",ntracks); fEvent =0; @@ -2674,22 +2885,39 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event) } -void AliTPCtrackerMI::DeleteSeeds() +Int_t AliTPCtrackerMI::PostProcess(AliESDEvent *event) { // - //delete Seeds + // Post process events + // + if (!event) return 0; - Int_t nseed = fSeeds->GetEntriesFast(); - for (Int_t i=0;iAt(i); - if (seed) delete fSeeds->RemoveAt(i); + // + // Set TPC event status + // + + // event affected by HV dip + // reset TPC status + if(IsTPCHVDipEvent(event)) { + event->ResetDetectorStatus(AliDAQ::kTPC); } - delete fSeeds; + + //printf("Status %d \n", event->IsDetectorOn(AliDAQ::kTPC)); + return 0; +} + + + void AliTPCtrackerMI::DeleteSeeds() +{ + // + fSeeds->Clear(); + ResetSeedsPool(); + delete fSeeds; fSeeds =0; } -void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction) +void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction) { // //read seeds from the event @@ -2712,7 +2940,8 @@ void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction) AliTPCtrack t(*esd); t.SetNumberOfClusters(0); // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha()); - AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/); + AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t/*,t.GetAlpha()*/); + seed->SetPoolID(fLastSeedID); seed->SetUniqueID(esd->GetID()); AddCovariance(seed); //add systematic ucertainty for (Int_t ikink=0;ikink<3;ikink++) { @@ -2733,11 +2962,11 @@ void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction) } if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.); if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.); - if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) { - fSeeds->AddAt(0,i); - delete seed; - continue; - } + //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) { + // fSeeds->AddAt(0,i); + // MarkSeedFree( seed ); + // continue; + //} if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) { Double_t par0[5],par1[5],alpha,x; esd->GetInnerExternalParameters(alpha,x,par0); @@ -2756,24 +2985,27 @@ void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction) // rotate to the local coordinate system // fSectors=fInnerSec; fN=fkNIS; - Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift(); + Double_t alpha=seed->GetAlpha(); if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi(); if (alpha < 0. ) alpha += 2.*TMath::Pi(); - Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN; + Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN; alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift(); + alpha-=seed->GetAlpha(); if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi(); if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi(); - alpha-=seed->GetAlpha(); - if (!seed->Rotate(alpha)) { - delete seed; - continue; + if (TMath::Abs(alpha) > 0.001) { // This should not happen normally + AliWarning(Form("Rotating track over %f",alpha)); + if (!seed->Rotate(alpha)) { + MarkSeedFree( seed ); + continue; + } } seed->SetESD(esd); // sign clusters if (esd->GetKinkIndex(0)<=0){ for (Int_t irow=0;irow<160;irow++){ Int_t index = seed->GetClusterIndex2(irow); - if (index>0){ + if (index >= 0){ // AliTPCclusterMI * cl = GetClusterMI(index); seed->SetClusterPointer(irow,cl); @@ -2816,7 +3048,8 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Double_t x[5], c[15]; // Int_t di = i1-i2; // - AliTPCseed * seed = new AliTPCseed(); + AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed(); + seed->SetPoolID(fLastSeedID); Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift(); Double_t cs=cos(alpha), sn=sin(alpha); // @@ -2985,12 +3218,12 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, } - Double_t dym = 0; - Double_t dzm = 0; - if (cm){ - dym = ym - cm->GetY(); - dzm = zm - cm->GetZ(); - } + // Double_t dym = 0; + // Double_t dzm = 0; + // if (cm){ + // dym = ym - cm->GetY(); + // dzm = zm - cm->GetZ(); + // } nin2++; @@ -3026,9 +3259,9 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue; UInt_t index=kr1.GetIndex(is); - seed->~AliTPCseed(); // this does not set the pointer to 0... - AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index); - + if (seed) {MarkSeedFree(seed); seed = 0;} + AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, ns*alpha+shift, x, c, index); + seed->SetPoolID(fLastSeedID); track->SetIsSeeding(kTRUE); track->SetSeed1(i1); track->SetSeed2(i2); @@ -3040,8 +3273,7 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Int_t foundable,found,shared; track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE); if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){ - seed->Reset(); - seed->~AliTPCseed(); + MarkSeedFree(seed); seed = 0; continue; } //} @@ -3059,8 +3291,7 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, if (track->GetNumberOfClusters()<(i1-i2)*0.5 || track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || track->GetNShared()>0.4*track->GetNumberOfClusters() ) { - seed->Reset(); - seed->~AliTPCseed(); + MarkSeedFree(seed); seed = 0; continue; } nout1++; @@ -3080,13 +3311,11 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, track2->SetBConstrain(kFALSE); track2->SetSeedType(1); arr->AddLast(track2); - seed->Reset(); - seed->~AliTPCseed(); + MarkSeedFree( seed ); seed = 0; continue; } else{ - seed->Reset(); - seed->~AliTPCseed(); + MarkSeedFree( seed ); seed = 0; continue; } @@ -3094,8 +3323,9 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, } track->SetSeedType(0); - arr->AddLast(track); - seed = new AliTPCseed; + arr->AddLast(track); // note, track is seed, don't free the seed + seed = new( NextFreeSeed() ) AliTPCseed; + seed->SetPoolID(fLastSeedID); nout2++; // don't consider other combinations if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8) @@ -3106,7 +3336,7 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, if (fDebug>3){ Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2); } - delete seed; + if (seed) MarkSeedFree( seed ); } @@ -3134,7 +3364,8 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Double_t x[5], c[15]; // // make temporary seed - AliTPCseed * seed = new AliTPCseed; + AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed; + seed->SetPoolID(fLastSeedID); Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift(); // Double_t cs=cos(alpha), sn=sin(alpha); // @@ -3197,7 +3428,7 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, y3 = kcl->GetY(); // apply angular cuts if (TMath::Abs(y1-y3)>dymax) continue; - x3 = x3; + //x3 = x3; z3 = kcl->GetZ(); if (TMath::Abs(z1-z3)>dzmax) continue; // @@ -3314,9 +3545,10 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue; - UInt_t index=kr1.GetIndex(is); - seed->~AliTPCseed(); - AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index); + index=kr1.GetIndex(is); + if (seed) {MarkSeedFree( seed ); seed = 0;} + AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, sec*alpha+shift, x, c, index); + seed->SetPoolID(fLastSeedID); track->SetIsSeeding(kTRUE); @@ -3324,8 +3556,7 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, FollowProlongation(*track, i1-7,1); if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 || track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){ - seed->Reset(); - seed->~AliTPCseed(); + MarkSeedFree( seed ); seed = 0; continue; } nout1++; @@ -3339,8 +3570,7 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, if (track->GetNumberOfClusters()<(i1-i2)*0.5 || track->GetNumberOfClusters()GetNFoundable()*0.7 || track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) { - seed->Reset(); - seed->~AliTPCseed(); + MarkSeedFree( seed ); seed = 0; continue; } @@ -3350,9 +3580,8 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, FollowProlongation(*track2, i2,1); track2->SetBConstrain(kFALSE); track2->SetSeedType(4); - arr->AddLast(track2); - seed->Reset(); - seed->~AliTPCseed(); + arr->AddLast(track2); + MarkSeedFree( seed ); seed = 0; } @@ -3363,9 +3592,9 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, } if (fDebug>3){ - Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3); + Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3); } - delete seed; + if (seed) MarkSeedFree(seed); } @@ -3397,7 +3626,8 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, AliTPCpolyTrack polytrack; Int_t nclusters=fSectors[sec][row0]; - AliTPCseed * seed = new AliTPCseed; + AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed; + seed->SetPoolID(fLastSeedID); Int_t sumused=0; Int_t cused=0; @@ -3476,9 +3706,9 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Int_t row = row0+ddrow*delta; kr = &(fSectors[sec][row]); Double_t xn = kr->GetX(); - Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5; + Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5; polytrack.GetFitPoint(xn,yn,zn); - if (TMath::Abs(yn)>ymax) continue; + if (TMath::Abs(yn)>ymax1) continue; nfoundable++; AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz); if (cln) { @@ -3595,8 +3825,9 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, UInt_t index=0; //kr0.GetIndex(is); - seed->~AliTPCseed(); - AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index); + if (seed) {MarkSeedFree( seed ); seed = 0;} + AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1,sec*alpha+shift,x,c,index); + seed->SetPoolID(fLastSeedID); track->SetIsSeeding(kTRUE); Int_t rc=FollowProlongation(*track, i2); if (constrain) track->SetBConstrain(1); @@ -3608,13 +3839,12 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 || track->GetNumberOfClusters() < track->GetNFoundable()*0.6 || track->GetNShared()>0.4*track->GetNumberOfClusters()) { - //delete track; - seed->Reset(); - seed->~AliTPCseed(); + MarkSeedFree( seed ); seed = 0; } else { - arr->AddLast(track); - seed = new AliTPCseed; + arr->AddLast(track); // track IS seed, don't free seed + seed = new( NextFreeSeed() ) AliTPCseed; + seed->SetPoolID(fLastSeedID); } nin3++; } @@ -3623,11 +3853,11 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, if (fDebug>3){ Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3); } - delete seed; + if (seed) MarkSeedFree( seed ); } -AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2) +AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2) { // // @@ -3654,7 +3884,7 @@ AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, if ( (indexGetX()>1){ clindex = track->GetClusterIndex2(i); - if (clindex>0){ + if (clindex >= 0){ x0[0] = trpoint->GetX(); x0[1] = trpoint->GetY(); x0[2] = trpoint->GetZ(); @@ -3665,7 +3895,7 @@ AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, if ( (indexGetX()>1)){ clindex = track->GetClusterIndex2(i); - if (clindex>0){ + if (clindex >= 0){ x1[0] = trpoint->GetX(); x1[1] = trpoint->GetY(); x1[2] = trpoint->GetZ(); @@ -3674,7 +3904,7 @@ AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, } if ( (indexGetX()>1)){ clindex = track->GetClusterIndex2(i); - if (clindex>0){ + if (clindex >= 0){ x2[0] = trpoint->GetX(); x2[1] = trpoint->GetY(); x2[2] = trpoint->GetZ(); @@ -3745,7 +3975,8 @@ AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43; // Int_t row1 = fSectors->GetRowNumber(x2[0]); - AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + seed->SetPoolID(fLastSeedID); // Double_t y0,z0,y1,z1, y2,z2; //seed->GetProlongation(x0[0],y0,z0); // seed->GetProlongation(x1[0],y1,z1); @@ -3759,7 +3990,7 @@ AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, } -AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2) +AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2) { // // @@ -3777,8 +4008,8 @@ AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, F ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point // // - Double_t xyz[3][3]; - Int_t row[3],sec[3]={0,0,0}; + Double_t xyz[3][3]={{0}}; + Int_t row[3]={0},sec[3]={0,0,0}; // // find track row position at given ratio of the length Int_t index=-1; @@ -3865,7 +4096,8 @@ AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, F c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43; // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]); - AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + seed->SetPoolID(fLastSeedID); seed->SetLastPoint(row[2]); seed->SetFirstPoint(row[2]); return seed; @@ -4014,7 +4246,8 @@ AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward) c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43; // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]); - AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0); + seed->SetPoolID(fLastSeedID); seed->SetLastPoint(row[2]); seed->SetFirstPoint(row[2]); for (Int_t i=row[0];i0) xm[i]/=Float_t(ncl); } - TTreeSRedirector &cstream = *fDebugStreamer; // for (Int_t i0=0;i0At(i0); @@ -4156,6 +4388,8 @@ void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_ } } // + if (AliTPCReconstructor::StreamLevel()>5) { + TTreeSRedirector &cstream = *fDebugStreamer; cstream<<"Multi"<< "iter="<1) { AliInfo("Time for curling tracks removal DEBUGGING MC"); timer.Print(); @@ -4202,44 +4439,38 @@ void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_ } -void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter) -{ - // - // - // Two reasons to have multiple find tracks - // 1. Curling tracks can be find more than once - // 2. Splitted tracks - // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached) - // b.) Edge effect on the sector boundaries - // - // This function tries to find tracks closed in the parametric space - // - // cut logic if distance is bigger than cut continue - Do Nothing - const Float_t kMaxdTheta = 0.05; // maximal distance in theta - const Float_t kMaxdPhi = 0.05; // maximal deistance in phi - const Float_t kdelta = 40.; //delta r to calculate track distance - // - // const Float_t kMaxDist0 = 1.; // maximal distance 0 - //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows - // - /* - TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2"); - TCut cdtheta("cdtheta","abs(dtheta)<0.05"); - */ - // + +void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){ // + // Find Splitted tracks and remove the one with worst quality + // Corresponding debug streamer to tune selections - "Splitted2" + // Algorithm: + // 0. Sort tracks according quility + // 1. Propagate the tracks to the reference radius + // 2. Double_t loop to select close tracks (only to speed up process) + // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold + // 4. Delete temporary parameters + // + const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary + // rough cuts + const Double_t kCutP1=10; // delta Z cut 10 cm + const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15 + const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15 + const Double_t kCutAlpha=0.15; // delta alpha cut + Int_t firstpoint = 0; + Int_t lastpoint = 160; // Int_t nentries = array->GetEntriesFast(); - AliHelix *helixes = new AliHelix[nentries]; - Float_t *xm = new Float_t[nentries]; + AliExternalTrackParam *params = new AliExternalTrackParam[nentries]; // // TStopwatch timer; timer.Start(); // - //Sort tracks according quality - // + //0. Sort tracks according quality + //1. Propagate the ext. param to reference radius Int_t nseed = array->GetEntriesFast(); + if (nseed<=0) return; Float_t * quality = new Float_t[nseed]; Int_t * indexes = new Int_t[nseed]; for (Int_t i=0; iGetNumberOfClusters(); //prefer high momenta tracks if overlaps - quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5); + quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5); + params[i]=(*pt); + AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE); + AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE); } TMath::Sort(nseed,quality,indexes); - - - // - // Find track COG in x direction - point with best defined parameters - // - for (Int_t i=0;iAt(i); - if (!track) continue; - track->SetCircular(0); - new (&helixes[i]) AliHelix(*track); - Int_t ncl=0; - xm[i]=0; - for (Int_t icl=0; icl<160; icl++){ - AliTPCclusterMI * cl = track->GetClusterPointer(icl); - if (cl) { - xm[i]+=cl->GetX(); - ncl++; - } - } - if (ncl>0) xm[i]/=Float_t(ncl); - } - TTreeSRedirector &cstream = *fDebugStreamer; // - for (Int_t is0=0;is0At(i0); - if (!track0) continue; - if (track0->GetKinkIndexes()[0]!=0) continue; - Float_t xc0 = helixes[i0].GetHelix(6); - Float_t yc0 = helixes[i0].GetHelix(7); - Float_t fi0 = TMath::ATan2(yc0,xc0); - - for (Int_t is1=is0+1;is1At(i1); - if (!track1) continue; - // - if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue; - if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue; - if (track1->GetKinkIndexes()[0]!=0) continue; - - Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl(); - if (TMath::Abs(dtheta)>kMaxdTheta) continue; + // 3. Loop over pair of tracks + // + for (Int_t i0=0; i0UncheckedAt(index0))) continue; + AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0); + if (!s1->IsActive()) continue; + AliExternalTrackParam &par0=params[index0]; + for (Int_t i1=i0+1; i1UncheckedAt(index1))) continue; + AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1); + if (!s2->IsActive()) continue; + if (s2->GetKinkIndexes()[0]!=0) + if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue; + AliExternalTrackParam &par1=params[index1]; + if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue; + if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue; + if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue; + Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha()); + if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi(); + if (TMath::Abs(dAlpha)>kCutAlpha) continue; // - Float_t xc1 = helixes[i1].GetHelix(6); - Float_t yc1 = helixes[i1].GetHelix(7); - Float_t fi1 = TMath::ATan2(yc1,xc1); + Int_t sumShared=0; + Int_t nall0=0; + Int_t nall1=0; + Int_t firstShared=lastpoint, lastShared=firstpoint; + Int_t firstRow=lastpoint, lastRow=firstpoint; // - Float_t dfi = fi0-fi1; - if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect - if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); // - if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){ - // - // if short tracks with undefined sign - fi1 = -TMath::ATan2(yc1,-xc1); - dfi = fi0-fi1; + for (Int_t i=firstpoint;iGetClusterIndex2(i)>0) nall0++; + if (s2->GetClusterIndex2(i)>0) nall1++; + if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) { + if (ilastRow) lastRow=i; + } + if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) { + if (ilastShared) lastShared=i; + sumShared++; + } + } + Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1)); + Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1)); + + if( AliTPCReconstructor::StreamLevel()>1){ + TTreeSRedirector &cstream = *fDebugStreamer; + Int_t n0=s1->GetNumberOfClusters(); + Int_t n1=s2->GetNumberOfClusters(); + Int_t n0F=s1->GetNFoundable(); + Int_t n1F=s2->GetNFoundable(); + Int_t lab0=s1->GetLabel(); + Int_t lab1=s2->GetLabel(); + + cstream<<"Splitted2"<< + "iter="<kMaxdPhi) continue; - // - // - Float_t sum =0; - Float_t sums=0; - Float_t sum0=0; - Float_t sum1=0; - for (Int_t icl=0; icl<160; icl++){ - Int_t index0=track0->GetClusterIndex2(icl); - Int_t index1=track1->GetClusterIndex2(icl); - Bool_t used0 = (index0>0 && !(index0&0x8000)); - Bool_t used1 = (index1>0 && !(index1&0x8000)); - // - if (used0) sum0++; // used cluster0 - if (used1) sum1++; // used clusters1 - if (used0&&used1) sum++; - if (index0==index1 && used0 && used1) sums++; - } - - // - if (sums<10) continue; - if (sum<40) continue; - if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue; - // - Double_t dist[5][4]; // distance at X - Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta - - // - // - track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz()); - for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]); - track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz()); - for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]); - // - track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz()); - for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]); - track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz()); - for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]); - // - track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz()); - for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]); - for (Int_t i=0;i<3;i++) mdist[i]*=0.2; // + // remove track with lower quality // - Int_t lab0=track0->GetLabel(); - Int_t lab1=track1->GetLabel(); - cstream<<"Splitted"<< - "iter="<RemoveAt(i1); + MarkSeedFree( array->RemoveAt(index1) ); + } } - } - delete [] helixes; - delete [] xm; - AliInfo("Time for splitted tracks removal"); - timer.Print(); + } + // + // 4. Delete temporary array + // + delete [] params; + delete [] quality; + delete [] indexes; + } -void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter) +void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter) { // // find Curling tracks @@ -4459,11 +4654,11 @@ void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_ // TStopwatch timer; timer.Start(); - Double_t phase[2][2],radius[2]; + Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0}; + // // Find tracks // - TTreeSRedirector &cstream = *fDebugStreamer; // for (Int_t i0=0;i0At(i0); @@ -4556,11 +4751,12 @@ void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_ track0->SetCircular(track0->GetCircular()+2); } } - if (AliTPCReconstructor::StreamLevel()>1){ + if (AliTPCReconstructor::StreamLevel()>2){ // //debug stream to tune "fine" cuts Int_t lab0=track0->GetLabel(); Int_t lab1=track1->GetLabel(); + TTreeSRedirector &cstream = *fDebugStreamer; cstream<<"Curling2"<< "iter="<GetKinkAngleCutChi2(0)){ + continue; + } + // kinks->AddLast(kink); kink = new AliKink; ncandidates++; @@ -5009,12 +5228,12 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) // for (Int_t i=0;iAt(i); + AliKink *kinkl = (AliKink*)kinks->At(i); // // refit kinks towards vertex // - Int_t index0 = kink->GetIndex(0); - Int_t index1 = kink->GetIndex(1); + Int_t index0 = kinkl->GetIndex(0); + Int_t index1 = kinkl->GetIndex(1); AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0); AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1); // @@ -5022,10 +5241,10 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) // // Refit Kink under if too small angle // - if (kink->GetAngle(2)<0.05){ - kink->SetTPCRow0(GetRowNumber(kink->GetR())); - Int_t row0 = kink->GetTPCRow0(); - Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); + if (kinkl->GetAngle(2)<0.05){ + kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR())); + Int_t row0 = kinkl->GetTPCRow0(); + Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2))); // // Int_t last = row0-drow; @@ -5040,31 +5259,31 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE); // if (seed0 && seed1){ - kink->SetStatus(1,8); - if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9); - row0 = GetRowNumber(kink->GetR()); + kinkl->SetStatus(1,8); + if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9); + row0 = GetRowNumber(kinkl->GetR()); sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters(); mothers[i] = *seed0; daughters[i] = *seed1; } else{ delete kinks->RemoveAt(i); - if (seed0) delete seed0; - if (seed1) delete seed1; + if (seed0) MarkSeedFree( seed0 ); + if (seed1) MarkSeedFree( seed1 ); continue; } - if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) { + if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) { delete kinks->RemoveAt(i); - if (seed0) delete seed0; - if (seed1) delete seed1; + if (seed0) MarkSeedFree( seed0 ); + if (seed1) MarkSeedFree( seed1 ); continue; } // - delete seed0; - delete seed1; + MarkSeedFree( seed0 ); + MarkSeedFree( seed1 ); } // - if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win + if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win } TMath::Sort(nkinks,quality,indexes,kFALSE); // @@ -5074,7 +5293,8 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]); if (!kink0) continue; // - for (Int_t ikink1=0;ikink1At(indexes[ikink0]); if (!kink0) continue; AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]); if (!kink1) continue; @@ -5108,7 +5328,8 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) } for (Int_t i=row0;i<158;i++){ - if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ + //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug? + if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){ both++; bothd++; if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){ @@ -5127,6 +5348,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) shared[kink0->GetIndex(0)]= kTRUE; shared[kink0->GetIndex(1)]= kTRUE; delete kinks->RemoveAt(indexes[ikink0]); + break; } else{ shared[kink1->GetIndex(0)]= kTRUE; @@ -5139,23 +5361,23 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) for (Int_t i=0;iAt(indexes[i]); - if (!kink) continue; - kink->SetTPCRow0(GetRowNumber(kink->GetR())); - Int_t index0 = kink->GetIndex(0); - Int_t index1 = kink->GetIndex(1); - if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue; - kink->SetMultiple(usage[index0],0); - kink->SetMultiple(usage[index1],1); - if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue; - if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue; - if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue; - if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue; + AliKink * kinkl = (AliKink*) kinks->At(indexes[i]); + if (!kinkl) continue; + kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR())); + Int_t index0 = kinkl->GetIndex(0); + Int_t index1 = kinkl->GetIndex(1); + if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue; + kinkl->SetMultiple(usage[index0],0); + kinkl->SetMultiple(usage[index1],1); + if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue; + if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue; + if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue; + if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue; AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0); AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1); if (!ktrack0 || !ktrack1) continue; - Int_t index = esd->AddKink(kink); + Int_t index = esd->AddKink(kinkl); // // if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink @@ -5177,7 +5399,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) AliTPCseed * track0 = (AliTPCseed*)array->At(i); if (!track0) continue; if (track0->GetKinkIndex(0)!=0) continue; - if (shared[i]) delete array->RemoveAt(i); + if (shared[i]) MarkSeedFree( array->RemoveAt(i) ); } // @@ -5196,16 +5418,16 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) if (!track0) continue; if (track0->Pt()<1.4) continue; //remove double high momenta tracks - overlapped with kink candidates - Int_t shared=0; + Int_t ishared=0; Int_t all =0; for (Int_t icl=track0->GetFirstPoint();iclGetLastPoint(); icl++){ if (track0->GetClusterPointer(icl)!=0){ all++; - if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++; + if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++; } } - if (Float_t(shared+1)/Float_t(all+1)>0.5) { - delete array->RemoveAt(i); + if (Float_t(ishared+1)/Float_t(all+1)>0.5) { + MarkSeedFree( array->RemoveAt(i) ); continue; } // @@ -5218,8 +5440,8 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) AliTPCseed & mother = *pmother; AliTPCseed & daughter = *pdaughter; - AliKink & kink = *pkink; - if (CheckKinkPoint(track0,mother,daughter, kink)){ + AliKink & kinkl = *pkink; + if (CheckKinkPoint(track0,mother,daughter, kinkl)){ if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) { delete pmother; delete pdaughter; @@ -5232,25 +5454,31 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) delete pkink; continue; } - Int_t row0= kink.GetTPCRow0(); - if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) { + Int_t row0= kinkl.GetTPCRow0(); + if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) { delete pmother; delete pdaughter; delete pkink; continue; } // - Int_t index = esd->AddKink(&kink); + Int_t index = esd->AddKink(&kinkl); mother.SetKinkIndex(0,-(index+1)); daughter.SetKinkIndex(0,index+1); if (mother.GetNumberOfClusters()>50) { - delete array->RemoveAt(i); - array->AddAt(new AliTPCseed(mother),i); + MarkSeedFree( array->RemoveAt(i) ); + AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother); + mtc->SetPoolID(fLastSeedID); + array->AddAt(mtc,i); } else{ - array->AddLast(new AliTPCseed(mother)); + AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother); + mtc->SetPoolID(fLastSeedID); + array->AddLast(mtc); } - array->AddLast(new AliTPCseed(daughter)); + AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter); + dtc->SetPoolID(fLastSeedID); + array->AddLast(dtc); for (Int_t icl=0;iclUse(20); } @@ -5287,43 +5515,52 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) kinks->Delete(); delete kinks; - printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall); + AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall)); timer.Print(); } -void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd) + +/* +void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) { // - // find V0s + // find kinks // // - TObjArray *tpcv0s = new TObjArray(100000); - Int_t nentries = array->GetEntriesFast(); + + TObjArray *kinks= new TObjArray(10000); + // TObjArray *v0s= new TObjArray(10000); + Int_t nentries = array->GetEntriesFast(); AliHelix *helixes = new AliHelix[nentries]; Int_t *sign = new Int_t[nentries]; + Int_t *nclusters = new Int_t[nentries]; Float_t *alpha = new Float_t[nentries]; + AliKink *kink = new AliKink(); + Int_t * usage = new Int_t[nentries]; + Float_t *zm = new Float_t[nentries]; Float_t *z0 = new Float_t[nentries]; - Float_t *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(); + Float_t *fim = new Float_t[nentries]; + Float_t *shared = new Float_t[nentries]; + Bool_t *circular = new Bool_t[nentries]; + Float_t *dca = new Float_t[nentries]; + //const AliESDVertex * primvertex = esd->GetVertex(); // // nentries = array->GetEntriesFast(); // + + // + // for (Int_t i=0;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 - // + track->SetCircular(0); + shared[i] = kFALSE; + track->UpdatePoints(); + if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){ + } + nclusters[i]=track->GetNumberOfClusters(); alpha[i] = track->GetAlpha(); new (&helixes[i]) AliHelix(*track); Double_t xyz[3]; @@ -5331,27 +5568,18 @@ void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd) sign[i] = (track->GetC()>0) ? -1:1; Double_t x,y,z; x=160; + if (track->GetProlongation(x,y,z)){ + zm[i] = z; + fim[i] = alpha[i]+TMath::ATan2(y,x); + } + else{ + zm[i] = track->GetZ(); + fim[i] = alpha[i]; + } z0[i]=1000; + circular[i]= kFALSE; if (track->GetProlongation(0,y,z)) z0[i] = z; - dca[i] = track->GetD(0,0); - // - // dca error parrameterezation + pulls - // - sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC())); - if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5; - cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3); - pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i]; - pulldcaz[i] = (z0[i]-zvertex)/sdcar[i]; - pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]); - if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) { - if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut - } - if (track->TPCrPID(4)>0.5) { - if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut - } - if (track->TPCrPID(0)>0.4) { - isPrim[i]=kFALSE; //electron no sigma cut - } + dca[i] = track->GetD(0,0); } // // @@ -5360,294 +5588,666 @@ void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd) Int_t ncandidates =0; Int_t nall =0; Int_t ntracks=0; - Double_t phase[2][2],radius[2]; + Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0}; + // - // Finf V0s loop - // + // Find circling track // - // // - TTreeSRedirector &cstream = *fDebugStreamer; - Float_t fprimvertex[3]={GetX(),GetY(),GetZ()}; - AliV0 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; - if (AliTPCReconstructor::StreamLevel()>1){ - cstream<<"Tracks"<< - "Tr0.="<GetSigned1Pt()<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); + for (Int_t i0=0;i0At(i0); + if (!track0) continue; + if (track0->GetNumberOfClusters()<40) continue; + if (TMath::Abs(1./track0->GetC())>200) continue; + for (Int_t i1=i0+1;i1At(i1); 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->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track - nall++; + if (track1->GetNumberOfClusters()<40) continue; + if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue; + if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; + if (TMath::Abs(1./track1->GetC())>200) continue; + if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; + if (track1->GetTgl()*track0->GetTgl()>0) continue; + if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue; + if (track0->GetBConstrain()&&track1->OneOverPt()OneOverPt()) continue; //returning - lower momenta + if (track1->GetBConstrain()&&track0->OneOverPt()OneOverPt()) continue; //returning - lower momenta + // + Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1])); + if (mindcar<5) continue; + Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ())); + if (mindcaz<5) continue; + if (mindcar+mindcaz<20) continue; + // // - // DCA to prim vertex cut + Float_t xc0 = helixes[i0].GetHelix(6); + Float_t yc0 = helixes[i0].GetHelix(7); + Float_t r0 = helixes[i0].GetHelix(8); + Float_t xc1 = helixes[i1].GetHelix(6); + Float_t yc1 = helixes[i1].GetHelix(7); + Float_t r1 = helixes[i1].GetHelix(8); + + Float_t rmean = (r0+r1)*0.5; + Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0)); + //if (delta>30) continue; + if (delta>rmean*0.25) continue; + if (TMath::Abs(r0-r1)/rmean>0.3) continue; // + Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10); + if (npoints==0) continue; + helixes[i0].GetClosestPhases(helixes[i1], phase); // - delta[0]=10000; - delta[1]=10000; + Double_t xyz0[3]; + Double_t xyz1[3]; + Double_t hangles[3]; + helixes[i0].Evaluate(phase[0][0],xyz0); + helixes[i1].Evaluate(phase[0][1],xyz1); - Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2); + helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles); + Double_t deltah[2],deltabest; + if (hangles[2]<2.8) continue; + if (npoints>0){ + Int_t ibest=0; + helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2); + if (npoints==2){ + helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2); + if (deltah[1]6) continue; + if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue; + Bool_t lsign =kFALSE; + if (hangles[2]>3.06) lsign =kTRUE; + // + if (lsign){ + circular[i0] = kTRUE; + circular[i1] = kTRUE; + if (track0->OneOverPt()OneOverPt()){ + track0->SetCircular(track0->GetCircular()+1); + track1->SetCircular(track1->GetCircular()+2); + } + else{ + track1->SetCircular(track1->GetCircular()+1); + track0->SetCircular(track0->GetCircular()+2); + } + } + if (lsign&&AliTPCReconstructor::StreamLevel()>1){ + //debug stream + Int_t lab0=track0->GetLabel(); + Int_t lab1=track1->GetLabel(); + TTreeSRedirector &cstream = *fDebugStreamer; + cstream<<"Curling"<< + "lab0="<At(i); + if (track0==0) { + AliInfo("seed==0"); + continue; + } + ntracks++; + // + Double_t cradius0 = 40*40; + Double_t cradius1 = 270*270; + Double_t cdist1=8.; + Double_t cdist2=8.; + Double_t cdist3=0.55; + for (Int_t j =i+1;j200) continue; + if ( (nclusters[i]+nclusters[j])<80) continue; + if ( TMath::Abs(zm[i]-zm[j])>60.) continue; + if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue; + //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2]; + Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20); if (npoints<1) continue; - 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 (pointAngleTPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate - Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle(); - //continue; - if (vertex.GetV0CosineOfPointingAngle()2&&(!isGamma)) continue; // point angle cut - if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // 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); - //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue; - if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue; - Float_t densb0=0,densb1=0,densa0=0,densa1=0; - Int_t row0 = GetRowNumber(vertex.GetRr()); - if (row0>15){ - //Bo: if (vertex.GetDist2()>0.2) continue; - if (vertex.GetDcaV0Daughters()>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{ + if ( (radius[0]cradius1) && (radius[1]cradius1) ) continue; + } + // + Double_t delta1=10000,delta2=10000; + // cuts on the intersection radius + helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1); + if (radius[0]<20&&delta1<1) continue; //intersection at vertex + if (radius[0]<10&&delta1<3) continue; //intersection at vertex + if (npoints==2){ + helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2); + if (radius[1]<20&&delta2<1) continue; //intersection at vertex + if (radius[1]<10&&delta2<3) continue; //intersection at vertex + } + // + Double_t distance1 = TMath::Min(delta1,delta2); + if (distance1>cdist1) continue; // cut on DCA linear approximation + // + npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20); + helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1); + if (radius[0]<20&&delta1<1) continue; //intersection at vertex + if (radius[0]<10&&delta1<3) continue; //intersection at vertex + // + if (npoints==2){ + helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2); + if (radius[1]<20&&delta2<1) continue; //intersection at vertex + if (radius[1]<10&&delta2<3) continue; //intersection at vertex + } + distance1 = TMath::Min(delta1,delta2); + Float_t rkink =0; + if (delta1cdist2) continue; + // + // + AliTPCseed * track1 = (AliTPCseed*)array->At(j); + // + // + Int_t row0 = GetRowNumber(rkink); + if (row0<10) continue; + if (row0>150) continue; + // + // + Float_t dens00=-1,dens01=-1; + Float_t dens10=-1,dens11=-1; + // + Int_t found,foundable,ishared; + track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE); + if (foundable>5) dens00 = Float_t(found)/Float_t(foundable); + track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE); + if (foundable>5) dens01 = Float_t(found)/Float_t(foundable); + // + track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE); + if (foundable>10) dens10 = Float_t(found)/Float_t(foundable); + track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE); + if (foundable>10) dens11 = Float_t(found)/Float_t(foundable); + // + if (dens00dens10 && dens01>dens11) continue; + if (TMath::Max(dens00,dens10)<0.1) continue; + if (TMath::Max(dens01,dens11)<0.3) continue; + // + if (TMath::Min(dens00,dens10)>0.6) continue; + if (TMath::Min(dens01,dens11)>0.6) continue; + + // + AliTPCseed * ktrack0, *ktrack1; + if (dens00>dens10){ + ktrack0 = track0; + ktrack1 = track1; + } + else{ + ktrack0 = track1; + ktrack1 = track0; + } + if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle + AliExternalTrackParam paramm(*ktrack0); + AliExternalTrackParam paramd(*ktrack1); + if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference()); + // + // + kink->SetMother(paramm); + kink->SetDaughter(paramd); + kink->Update(); + + Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]}; + Int_t index[4]; + fkParam->Transform0to1(x,index); + fkParam->Transform1to2(x,index); + row0 = GetRowNumber(x[0]); + + if (kink->GetR()<100) continue; + if (kink->GetR()>240) continue; + if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume + if (kink->GetDistance()>cdist3) continue; + Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate + if (dird<0) continue; + + Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate + if (dirm<0) continue; + Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]); + if (mpt<0.2) continue; + + if (mpt<1){ + //for high momenta momentum not defined well in first iteration + Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP(); + if (qt>0.35) continue; + } + + kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0); + kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1); + if (dens00>dens10){ + kink->SetTPCDensity(dens00,0,0); + kink->SetTPCDensity(dens01,0,1); + kink->SetTPCDensity(dens10,1,0); + kink->SetTPCDensity(dens11,1,1); + kink->SetIndex(i,0); + kink->SetIndex(j,1); + } + else{ + kink->SetTPCDensity(dens10,0,0); + kink->SetTPCDensity(dens11,0,1); + kink->SetTPCDensity(dens00,1,0); + kink->SetTPCDensity(dens01,1,1); + kink->SetIndex(j,0); + kink->SetIndex(i,1); + } + + if (mpt<1||kink->GetAngle(2)>0.1){ + // angle and densities not defined yet + if (kink->GetTPCDensityFactor()<0.8) continue; + if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue; + if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle + if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue; + if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue; + + Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2(); + criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2(); + criticalangle= 3*TMath::Sqrt(criticalangle); + if (criticalangle>0.02) criticalangle=0.02; + if (kink->GetAngle(2)GetAngle(2))); // overlap region defined + Float_t shapesum =0; + Float_t sum = 0; + for ( Int_t row = row0-drow; row155) continue; + if (ktrack0->GetClusterPointer(row)){ + AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row); + shapesum+=point->GetSigmaY()+point->GetSigmaZ(); + sum++; + } + if (ktrack1->GetClusterPointer(row)){ + AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row); + shapesum+=point->GetSigmaY()+point->GetSigmaZ(); + sum++; + } + } + if (sum<4){ + kink->SetShapeFactor(-1.); + } + else{ + kink->SetShapeFactor(shapesum/sum); + } + // esd->AddKink(kink); + // + // kink->SetMother(paramm); + //kink->SetDaughter(paramd); + + Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2]; + chi2P2*=chi2P2; + chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5]; + Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3]; + chi2P3*=chi2P3; + chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9]; + // + if (AliTPCReconstructor::StreamLevel()>1) { + (*fDebugStreamer)<<"kinkLpt"<< + "chi2P2="<GetKinkAngleCutChi2(0)){ + continue; + } + // + kinks->AddLast(kink); + kink = new AliKink; + ncandidates++; + } + } + // + // sort the kinks according quality - and refit them towards vertex + // + Int_t nkinks = kinks->GetEntriesFast(); + Float_t *quality = new Float_t[nkinks]; + Int_t *indexes = new Int_t[nkinks]; + AliTPCseed **mothers = new AliTPCseed*[nkinks]; memset(mothers, 0, nkinks*sizeof(AliTPCseed*)); + AliTPCseed **daughters = new AliTPCseed*[nkinks]; memset(daughters, 0, nkinks*sizeof(AliTPCseed*)); + // + // + for (Int_t i=0;iAt(i); + // + // refit kinks towards vertex + // + Int_t index0 = kinkl->GetIndex(0); + Int_t index1 = kinkl->GetIndex(1); + AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0); + AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1); + // + Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters(); + // + // Refit Kink under if too small angle + // + if (kinkl->GetAngle(2)<0.05){ + kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR())); + Int_t row0 = kinkl->GetTPCRow0(); + Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2))); + // + // + Int_t last = row0-drow; + if (last<40) last=40; + if (lastGetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25; + AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE); + // + // + Int_t first = row0+drow; + if (first>130) first=130; + if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30); + AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE); + // + if (seed0 && seed1){ + kinkl->SetStatus(1,8); + if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9); + row0 = GetRowNumber(kinkl->GetR()); + sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters(); + mothers[i] = new ( NextFreeSeed() ) AliTPCseed(*seed0); + mothers[i]->SetPoolID(fLastSeedID); + daughters[i] = new (NextFreeSeed() ) AliTPCseed(*seed1); + daughters[i]->SetPoolID(fLastSeedID); + } + else{ + delete kinks->RemoveAt(i); + if (seed0) MarkSeedFree( seed0 ); + if (seed1) MarkSeedFree( seed1 ); + continue; + } + if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) { + delete kinks->RemoveAt(i); + if (seed0) MarkSeedFree( seed0 ); + if (seed1) MarkSeedFree( seed1 ); + continue; + } + // + MarkSeedFree( seed0 ); + MarkSeedFree( seed1 ); + } + // + if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win + } + TMath::Sort(nkinks,quality,indexes,kFALSE); + // + //remove double find kinks + // + for (Int_t ikink0=1;ikink0At(indexes[ikink0]); + if (!kink0) continue; + // + for (Int_t ikink1=0;ikink1At(indexes[ikink0]); + if (!kink0) continue; + AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]); + if (!kink1) continue; + // if not close kink continue + if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue; + if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue; + if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue; + // + AliTPCseed &mother0 = *mothers[indexes[ikink0]]; + AliTPCseed &daughter0 = *daughters[indexes[ikink0]]; + AliTPCseed &mother1 = *mothers[indexes[ikink1]]; + AliTPCseed &daughter1 = *daughters[indexes[ikink1]]; + Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2; + // + Int_t same = 0; + Int_t both = 0; + Int_t samem = 0; + Int_t bothm = 0; + Int_t samed = 0; + Int_t bothd = 0; + // + for (Int_t i=0;i0 && mother1.GetClusterIndex(i)>0){ + both++; + bothm++; + if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){ + same++; + samem++; + } + } + } + + for (Int_t i=row0;i<158;i++){ + //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug? + if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){ + both++; + bothd++; + if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){ + same++; + samed++; + } + } } - else{ - densa0 = track0->Density2(0,40); //cluster density - densa1 = track1->Density2(0,40); //cluster density - if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue; - } -//Bo: vertex.SetLab(0,track0->GetLabel()); -//Bo: 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); -//Bo: vertex.SetStatus(1); // TPC v0 candidate - vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate -//Bo: vertex.SetRp(track0->TPCrPIDs()); -//Bo: vertex.SetRm(track1->TPCrPIDs()); - tpcv0s->AddLast(new AliESDv0(vertex)); - ncandidates++; - { - Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number - Double_t radiusm= (delta[0]1) { - Int_t lab0=track0->GetLabel(); - Int_t lab1=track1->GetLabel(); - Char_t c0=track0->GetCircular(); - Char_t c1=track1->GetCircular(); - cstream<<"V0"<< - "Event="<0.3 && ratiom>0.5 &&ratiod>0.5) { + Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters(); + Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters(); + if (sum1>sum0){ + shared[kink0->GetIndex(0)]= kTRUE; + shared[kink0->GetIndex(1)]= kTRUE; + delete kinks->RemoveAt(indexes[ikink0]); + break; + } + else{ + shared[kink1->GetIndex(0)]= kTRUE; + shared[kink1->GetIndex(1)]= kTRUE; + delete kinks->RemoveAt(indexes[ikink1]); + } } } - } - Float_t *quality = new Float_t[ncandidates]; - Int_t *indexes = new Int_t[ncandidates]; - Int_t naccepted =0; - for (Int_t i=0;iAt(i); - quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle - // quality[i] /= (0.5+v0->GetDist2()); - // quality[i] *= v0->GetChi2After(); //density factor - - Int_t index0 = v0->GetIndex(0); - Int_t index1 = v0->GetIndex(1); - //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull - Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo: + } + for (Int_t i=0;iAt(indexes[i]); + if (!kinkl) continue; + kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR())); + Int_t index0 = kinkl->GetIndex(0); + Int_t index1 = kinkl->GetIndex(1); + if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue; + kinkl->SetMultiple(usage[index0],0); + kinkl->SetMultiple(usage[index1],1); + if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue; + if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue; + if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue; + if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue; - AliTPCseed * track0 = (AliTPCseed*)array->At(index0); - AliTPCseed * track1 = (AliTPCseed*)array->At(index1); - if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate - if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate + AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0); + AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1); + if (!ktrack0 || !ktrack1) continue; + Int_t index = esd->AddKink(kinkl); + // + // + if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink + if (mothers[indexes[i]]->GetNumberOfClusters()>20 && daughters[indexes[i]]->GetNumberOfClusters()>20 && + (mothers[indexes[i]]->GetNumberOfClusters()+daughters[indexes[i]]->GetNumberOfClusters())>100){ + *ktrack0 = *mothers[indexes[i]]; + *ktrack1 = *daughters[indexes[i]]; + } + } + // + ktrack0->SetKinkIndex(usage[index0],-(index+1)); + ktrack1->SetKinkIndex(usage[index1], (index+1)); + usage[index0]++; + usage[index1]++; + } + // + // Remove tracks corresponding to shared kink's + // + for (Int_t i=0;iAt(i); + if (!track0) continue; + if (track0->GetKinkIndex(0)!=0) continue; + if (shared[i]) MarkSeedFree( array->RemoveAt(i) ); } - 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"); + RemoveUsed2(array,0.5,0.4,30); + UnsignClusters(); + for (Int_t i=0;iAt(i); + if (!track0) continue; + track0->CookdEdx(0.02,0.6); + track0->CookPID(); + } + // + for (Int_t i=0;iAt(i); + if (!track0) continue; + if (track0->Pt()<1.4) continue; + //remove double high momenta tracks - overlapped with kink candidates + Int_t ishared=0; + Int_t all =0; + for (Int_t icl=track0->GetFirstPoint();iclGetLastPoint(); icl++){ + if (track0->GetClusterPointer(icl)!=0){ + all++; + if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++; + } + } + if (Float_t(ishared+1)/Float_t(all+1)>0.5) { + MarkSeedFree( array->RemoveAt(i) ); 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;} - // - AliESDv0 * v02 = v0; - if (accept){ - //Bo: v0->SetOrder(0,order0); - //Bo: v0->SetOrder(1,order1); - //Bo: v0->SetOrder(1,order0+order1); - v0->SetOnFlyStatus(kTRUE); - Int_t index = esd->AddV0(v0); - v02 = esd->GetV0(index); - v0indexes0[order0]=index; - v0indexes1[order1]=index; - naccepted++; - } - { - Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number - if (AliTPCReconstructor::StreamLevel()>1) { - Int_t lab0=track0->GetLabel(); - Int_t lab1=track1->GetLabel(); - cstream<<"V02"<< - "Event="<GetKinkIndex(0)!=0) continue; + if (track0->GetNumberOfClusters()<80) continue; + + AliTPCseed *pmother = new( NextFreeSeed() ) AliTPCseed(); + pmother->SetPoolID(fLastSeedID); + AliTPCseed *pdaughter = new( NextFreeSeed() ) AliTPCseed(); + pdaughter->SetPoolID(fLastSeedID); + AliKink *pkink = new AliKink; + AliTPCseed & mother = *pmother; + AliTPCseed & daughter = *pdaughter; + AliKink & kinkl = *pkink; + if (CheckKinkPoint(track0,mother,daughter, kinkl)){ + if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) { + MarkSeedFree( pmother ); + MarkSeedFree( pdaughter ); + delete pkink; + continue; //too short tracks + } + if (mother.Pt()<1.4) { + MarkSeedFree( pmother ); + MarkSeedFree( pdaughter ); + delete pkink; + continue; + } + Int_t row0= kinkl.GetTPCRow0(); + if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) { + MarkSeedFree( pmother ); + MarkSeedFree( pdaughter ); + delete pkink; + continue; + } + // + Int_t index = esd->AddKink(&kinkl); + mother.SetKinkIndex(0,-(index+1)); + daughter.SetKinkIndex(0,index+1); + if (mother.GetNumberOfClusters()>50) { + MarkSeedFree( array->RemoveAt(i) ); + AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother); + mtc->SetPoolID(fLastSeedID); + array->AddAt(mtc,i); + } + else{ + AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother); + mtc->SetPoolID(fLastSeedID); + array->AddLast(mtc); + } + AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter); + dtc->SetPoolID(fLastSeedID); + array->AddLast(dtc); + for (Int_t icl=0;iclUse(20); + } + // + for (Int_t icl=row0;icl<158;icl++) { + if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20); + } + // + } + MarkSeedFree( pmother ); + MarkSeedFree( pdaughter ); + delete pkink; + } + delete [] daughters; + delete [] mothers; // // + delete [] dca; + delete []circular; + delete []shared; delete []quality; delete []indexes; -// - delete [] isPrim; - delete [] pulldca; - delete [] pulldcaz; - delete [] pulldcar; - delete [] cdcar; - delete [] sdcar; - delete [] dca; // + delete kink; + delete[]fim; + delete[] zm; delete[] z0; + delete [] usage; delete[] alpha; + delete[] nclusters; delete[] sign; delete[] helixes; - printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall); + kinks->Delete(); + delete kinks; + + AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall)); timer.Print(); } +*/ -Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk) +Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk) { // // refit kink towards to the vertex @@ -5789,7 +6389,7 @@ void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){ } -Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk) +Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk) { // // check kink point for given track @@ -5813,7 +6413,8 @@ Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTP seed1->Reset(kTRUE); last = seed1->GetLastPoint(); // - AliTPCseed *seed0 = new AliTPCseed(*seed); + AliTPCseed *seed0 = new( NextFreeSeed() ) AliTPCseed(*seed); + seed0->SetPoolID(fLastSeedID); seed0->Reset(kFALSE); seed0->Reset(); // @@ -5853,13 +6454,15 @@ Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTP // } } - delete seed0; - delete seed1; + MarkSeedFree( seed0 ); + MarkSeedFree( seed1 ); if (index<0) return 0; // Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate - seed0 = new AliTPCseed(param0[index]); - seed1 = new AliTPCseed(param1[index]); + seed0 = new( NextFreeSeed() ) AliTPCseed(param0[index]); + seed0->SetPoolID(fLastSeedID); + seed1 = new( NextFreeSeed() ) AliTPCseed(param1[index]); + seed1->SetPoolID(fLastSeedID); seed0->Reset(kFALSE); seed1->Reset(kFALSE); seed0->ResetCovariance(10.); @@ -5912,15 +6515,40 @@ Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTP // // if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){ - delete seed0; - delete seed1; + MarkSeedFree( seed0 ); + MarkSeedFree( seed1 ); return 0; } + // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33); kink.SetMother(param0[index]); kink.SetDaughter(param1[index]); kink.Update(); + + Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2]; + chi2P2*=chi2P2; + chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5]; + Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3]; + chi2P3*=chi2P3; + chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9]; + // + if (AliTPCReconstructor::StreamLevel()>1) { + (*fDebugStreamer)<<"kinkHpt"<< + "chi2P2="<GetKinkAngleCutChi2(0)){ + MarkSeedFree( seed0 ); + MarkSeedFree( seed1 ); + return 0; + } + + row0 = GetRowNumber(kink.GetR()); kink.SetTPCRow0(row0); kink.SetLabel(CookLabel(seed0,0.5,0,row0),0); @@ -5935,11 +6563,14 @@ Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTP daughter = param1[index]; daughter.SetLabel(kink.GetLabel(1)); param0[index].Reset(kTRUE); - FollowProlongation(param0[index],0); + FollowProlongation(param0[index],0); mother = param0[index]; mother.SetLabel(kink.GetLabel(0)); - delete seed0; - delete seed1; + if ( chi2P2+chi2P3GetKinkAngleCutChi2(1)){ + mother=*seed; + } + MarkSeedFree( seed0 ); + MarkSeedFree( seed1 ); // return 1; } @@ -6003,23 +6634,33 @@ Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) { Int_t n=(Int_t)seedTree->GetEntries(); for (Int_t i=0; iGetEvent(i); - fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/)); + AliTPCseed* sdc = new( NextFreeSeed() ) AliTPCseed(*seed/*,seed->GetAlpha()*/); + sdc->SetPoolID(fLastSeedID); + fSeeds->AddLast(sdc); } - delete seed; + delete seed; // RS: this seed is not from the pool, delete it !!! delete seedTree; savedir->cd(); return 0; } -Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd) +Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd) { // + // clusters to tracks if (fSeeds) DeleteSeeds(); - fEvent = esd; + else ResetSeedsPool(); + fEvent = esd; + + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + transform->SetCurrentTimeStamp( esd->GetTimeStamp()); + transform->SetCurrentRun(esd->GetRunNumber()); + Clusters2Tracks(); if (!fSeeds) return 1; FillESD(fSeeds); + if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds); return 0; // } @@ -6045,16 +6686,16 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { if (!pt) continue; Int_t nc=t.GetNumberOfClusters(); if (nc<20) { - delete fSeeds->RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); continue; } CookLabel(pt,0.1); if (pt->GetRemoval()==10) { if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7) - pt->Desactivate(10); // make track again active + pt->Desactivate(10); // make track again active // MvL: should be 0 ? else{ pt->Desactivate(20); - delete fSeeds->RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); } } } @@ -6062,9 +6703,10 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { RemoveUsed2(fSeeds,0.85,0.85,0); if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent); //FindCurling(fSeeds, fEvent,0); - if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks + if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks RemoveUsed2(fSeeds,0.5,0.4,20); FindSplitted(fSeeds, fEvent,0); // find multi found tracks + if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks // // // // refit short tracks @@ -6077,7 +6719,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { if (!pt) continue; Int_t nc=t.GetNumberOfClusters(); if (nc<15) { - delete fSeeds->RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); continue; } CookLabel(pt,0.1); //For comparison only @@ -6088,7 +6730,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { pt->SetLab2(i); } else - delete fSeeds->RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); } @@ -6103,7 +6745,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { if (!pt) continue; Int_t nc=t.GetNumberOfClusters(); if (nc<15) { - delete fSeeds->RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); continue; } t.SetUniqueID(i); @@ -6118,7 +6760,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { pt->SetLab2(i); } else - delete fSeeds->RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1); //if (seed1){ // FollowProlongation(*seed1,0); @@ -6142,7 +6784,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { fIteration = 2; PrepareForProlongation(fSeeds,5.); - PropagateForward2(fSeeds); + PropagateForard2(fSeeds); printf("Time for FORWARD propagation: \t");timer.Print();timer.Start(); // RemoveUsed(fSeeds,0.7,0.7,6); @@ -6155,7 +6797,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { if (!pt) continue; Int_t nc=t.GetNumberOfClusters(); if (nc<15) { - delete fSeeds->RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); continue; } t.CookdEdx(0.02,0.6); @@ -6165,7 +6807,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() { cerr<RemoveAt(i); + MarkSeedFree( fSeeds->RemoveAt(i) ); pt->fLab2 = i; } */ @@ -6200,7 +6842,8 @@ TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_ // // //tracking routine - TObjArray * arr = new TObjArray; + static TObjArray arrTracks; + TObjArray * arr = &arrTracks; // fSectors = fOuterSec; TStopwatch timer; @@ -6225,7 +6868,7 @@ TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_ TObjArray * AliTPCtrackerMI::Tracking() { - // + // tracking // if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial(); TStopwatch timer; @@ -6234,6 +6877,9 @@ TObjArray * AliTPCtrackerMI::Tracking() TObjArray * seeds = new TObjArray; TObjArray * arr=0; + Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec(); + Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim(); + Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec(); Int_t gap =20; Float_t cuts[4]; @@ -6247,7 +6893,7 @@ TObjArray * AliTPCtrackerMI::Tracking() // //find primaries cuts[0]=0.0066; - for (Int_t delta = 0; delta<18; delta+=6){ + for (Int_t delta = 0; delta<18; delta+=gapPrim){ // cuts[0]=0.0070; cuts[1] = 1.5; @@ -6272,7 +6918,7 @@ TObjArray * AliTPCtrackerMI::Tracking() //find primaries cuts[0]=0.0077; - for (Int_t delta = 20; delta<120; delta+=10){ + for (Int_t delta = 20; delta<120; delta+=gapPrim){ // // seed high pt tracks cuts[0]=0.0060; @@ -6325,9 +6971,22 @@ TObjArray * AliTPCtrackerMI::Tracking() SumTracks(seeds,arr); SignClusters(seeds,fnumber,fdensity); // + arr = Tracking(4,nup-5,nup-5-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + arr = Tracking(4,nup-7,nup-7-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // + // + arr = Tracking(4,nup-9,nup-9-gap,cuts,-1); + SumTracks(seeds,arr); + SignClusters(seeds,fnumber,fdensity); + // - for (Int_t delta = 3; delta<30; delta+=5){ + for (Int_t delta = 9; delta<30; delta+=gapSec){ // cuts[0] = 0.3; cuts[1] = 1.5; @@ -6350,10 +7009,9 @@ TObjArray * AliTPCtrackerMI::Tracking() fdensity = 2.; cuts[0]=0.0080; - Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec(); // find secondaries - for (Int_t delta = 30; deltaGetEntriesFast(); for (Int_t i=0;iUncheckedAt(i); @@ -6436,12 +7095,12 @@ void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const // remove tracks with too big curvature // if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){ - delete arr2->RemoveAt(i); + MarkSeedFree( arr2->RemoveAt(i) ); continue; } // REMOVE VERY SHORT TRACKS if (pt->GetNumberOfClusters()<20){ - delete arr2->RemoveAt(i); + MarkSeedFree( arr2->RemoveAt(i) ); continue; }// patch 28 fev06 // NORMAL ACTIVE TRACK @@ -6451,7 +7110,7 @@ void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const } //remove not usable tracks if (pt->GetRemoval()!=10){ - delete arr2->RemoveAt(i); + MarkSeedFree( arr2->RemoveAt(i) ); continue; } @@ -6459,16 +7118,16 @@ void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7) arr1->AddLast(arr2->RemoveAt(i)); else{ - delete arr2->RemoveAt(i); + MarkSeedFree( arr2->RemoveAt(i) ); } } } - delete arr2; arr2 = 0; + // delete arr2; arr2 = 0; // RS: this is static array, don't delete it } -void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast) +void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast) { // // try to track in parralel @@ -6480,7 +7139,7 @@ void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rla if (!pt) continue; if (!t.IsActive()) continue; // follow prolongation to the first layer - if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) ) + if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) ) FollowProlongation(t, rfirst+1); } @@ -6499,7 +7158,7 @@ void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rla if (!pt) continue; if (nr==80) pt->UpdateReference(); if (!pt->IsActive()) continue; - // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())fFirstPoint-fkParam->GetNRowLow())GetRelativeSector()>17) { continue; } @@ -6510,7 +7169,7 @@ void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rla AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i); if (!pt) continue; if (!pt->IsActive()) continue; - // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())fFirstPoint-fkParam->GetNRowLow())GetRelativeSector()>17) { continue; } @@ -6519,7 +7178,7 @@ void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rla } } -void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const +void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const { // // @@ -6541,7 +7200,7 @@ void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) co Float_t angle2 = pt->GetAlpha(); if (TMath::Abs(angle1-angle2)>0.001){ - pt->Rotate(angle1-angle2); + if (!pt->Rotate(angle1-angle2)) return; //angle2 = pt->GetAlpha(); //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha(); //if (pt->GetAlpha()<0) @@ -6555,7 +7214,7 @@ void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) co } -void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const +void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const { // // @@ -6574,7 +7233,7 @@ void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const } -Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) +Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr) { // // make back propagation @@ -6587,7 +7246,7 @@ Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) fSectors = fInnerSec; //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1); //fSectors = fOuterSec; - FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); + FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1); //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){ // Error("PropagateBack","Not prolonged track %d",pt->GetLabel()); // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); @@ -6597,7 +7256,7 @@ Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1); pt->SetFirstPoint(kink->GetTPCRow0()); fSectors = fInnerSec; - FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1); + FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1); } CookLabel(pt,0.3); } @@ -6605,7 +7264,7 @@ Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) } -Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr) +Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr) { // // make forward propagation @@ -6615,12 +7274,12 @@ Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr) for (Int_t i=0;iUncheckedAt(i); if (pt) { - FollowProlongation(*pt,0); + FollowProlongation(*pt,0,1,1); CookLabel(pt,0.3); } } - return 0; + return 0; } @@ -6653,7 +7312,7 @@ Int_t AliTPCtrackerMI::PropagateForward() -Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1) +Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1) { // // make back propagation, in between row0 and row1 @@ -6687,27 +7346,32 @@ Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1) void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row) { - // + // gets cluster shape // AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam(); - Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ()))); - Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2; + Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()))); + Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2; Double_t angulary = seed->GetSnp(); - angulary = angulary*angulary/(1.-angulary*angulary); + + if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) { + angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary); + } + + angulary = angulary*angulary/((1.-angulary)*(1.+angulary)); Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary); Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary))); Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz))); seed->SetCurrentSigmaY2(sigmay*sigmay); seed->SetCurrentSigmaZ2(sigmaz*sigmaz); - // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); -// // Float_t padlength = fParam->GetPadPitchLength(seed->fSector); + // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL(); +// // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector); // Float_t padlength = GetPadPitchLength(row); // // -// Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3; +// Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3; // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy); // // -// Float_t sresz = fParam->GetZSigma(); +// Float_t sresz = fkParam->GetZSigma(); // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz); /* Float_t wy = GetSigmaY(seed); @@ -6784,14 +7448,14 @@ void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const { if (TMath::Abs(c->GetLabel(1)) == lab || TMath::Abs(c->GetLabel(2)) == lab ) max++; } - - if ((1.- Float_t(max)/noc) > wrong) lab=-lab; + if (noc<=0) { lab=-1; return;} + if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab; else { Int_t tail=Int_t(0.10*noc); max=0; Int_t ind=0; - for (i=1; i<=160&&indGetLabel(1)) == lab || TMath::Abs(c->GetLabel(2)) == lab ) max++; } - - if ((1.- Float_t(max)/noc) > wrong) lab=-lab; + if (noc<=0) { lab=-1; return -1;} + if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab; else { Int_t tail=Int_t(0.10*noc); max=0; Int_t ind=0; - for (i=1; i<=160&&indGetTrackPoint(iter); - if (point) { - t->SetClusterMapBit(iter, kTRUE); + // Change to cluster pointers to see if we have a cluster at given padrow + + cluster = t->GetClusterPointer(iter); + if (cluster) { + clusterMap.SetBitNumber(iter, kTRUE); + point = t->GetTrackPoint(iter); if (point->IsShared()) - t->SetSharedMapBit(iter,kTRUE); - else - t->SetSharedMapBit(iter, kFALSE); + sharedMap.SetBitNumber(iter,kTRUE); } - else { - t->SetClusterMapBit(iter, kFALSE); - t->SetSharedMapBit(iter, kFALSE); + if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) { + fitMap.SetBitNumber(iter, kTRUE); + nclsf++; } } + esd->SetTPCClusterMap(clusterMap); + esd->SetTPCSharedMap(sharedMap); + esd->SetTPCFitMap(fitMap); + if (nclsf != t->GetNumberOfClusters()) + AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits())); } +Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){ + // + // return flag if there is findable cluster at given position + // + Float_t kDeltaZ=10; + Float_t z = track.GetZ(); + + if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) && + TMath::Abs(z)GetZLength(0) && + (TMath::Abs(track.GetSnp())GetSystematicError(); + // + // use only the diagonal part if not specified otherwise + if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed); + // + Double_t *covarS= (Double_t*)seed->GetCovariance(); + Double_t factor[5]={1,1,1,1,1}; + factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0])); + factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2])); + factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5])); + factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9])); + factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] +param[4]*param[4])/covarS[14])); + // + factor[0]=factor[2]; + factor[4]=factor[2]; + // 0 + // 1 2 + // 3 4 5 + // 6 7 8 9 + // 10 11 12 13 14 + for (Int_t i=0; i<5; i++){ + for (Int_t j=i; j<5; j++){ + Int_t index=seed->GetIndex(i,j); + covarS[index]*=factor[i]*factor[j]; + } + } +} + + +void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){ + // + // Adding systematic error - as additive factor without correlation + // + // !!!! the systematic error for element 4 is in 1/GeV + // 03.03.2012 MI changed in respect to the previous versions const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError(); + Double_t *covarIn= (Double_t*)seed->GetCovariance(); Double_t covar[15]; for (Int_t i=0;i<15;i++) covar[i]=0; // 0 @@ -6954,7 +7679,108 @@ void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){ covar[2] = param[1]*param[1]; covar[5] = param[2]*param[2]; covar[9] = param[3]*param[3]; - Double_t facC = AliTracker::GetBz()*kB2C; - covar[14]= param[4]*param[4]*facC*facC; + covar[14]= param[4]*param[4]; + // + covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2])); + // + covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5])); + covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5])); + // + covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9])); + covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9])); + covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9])); + // + covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14])); + covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14])); + covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14])); + covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14])); + // seed->AddCovariance(covar); } + +//_____________________________________________________________________________ +Bool_t AliTPCtrackerMI::IsTPCHVDipEvent(AliESDEvent const *esdEvent) { +// +// check events affected by TPC HV dip +// +if(!esdEvent) return kFALSE; + +// Init TPC OCDB +if(!AliTPCcalibDB::Instance()) return kFALSE; +AliTPCcalibDB::Instance()->SetRun(esdEvent->GetRunNumber()); + +// Get HV TPC chamber sensors and calculate the median +AliDCSSensorArray *voltageArray= AliTPCcalibDB::Instance()->GetVoltageSensors(esdEvent->GetRunNumber()); +if(!voltageArray) return kFALSE; + +TString sensorName=""; +Double_t kTPCHVdip = 2.0; // allow for 2V dip as compared to median from given sensor + + + for(Int_t sector=0; sector<72; sector++) + { + Char_t sideName='A'; + if ((sector/18)%2==1) sideName='C'; + if (sector<36){ + //IROC + sensorName=Form("TPC_ANODE_I_%c%02d_VMEAS",sideName,sector%18); + } else { + //OROC + sensorName=Form("TPC_ANODE_O_%c%02d_0_VMEAS",sideName,sector%18); + } + + AliDCSSensor* sensor = voltageArray->GetSensor(sensorName.Data()); + if(!sensor) continue; + TGraph *graph = sensor->GetGraph(); + if(!graph) continue; + Double_t median = TMath::Median(graph->GetN(), graph->GetY()); + if(median == 0) continue; + + //printf("chamber %d, sensor %s, HV %f, median %f\n", sector, sensorName.Data(), sensor->GetValue(esdEvent->GetTimeStamp()), median); + + if(TMath::Abs(sensor->GetValue(esdEvent->GetTimeStamp())-median)>kTPCHVdip) { + return kTRUE; + } + } + + return kFALSE; +} + +//________________________________________ +void AliTPCtrackerMI::MarkSeedFree(TObject *sd) +{ + // account that this seed is "deleted" + AliTPCseed* seed = dynamic_cast(sd); + if (!seed) { + AliError(Form("Freeing of non-AliTPCseed %p from the pool is requested",sd)); + return; + } + int id = seed->GetPoolID(); + if (id<0) { + AliError(Form("Freeing of seed %p NOT from the pool is requested",sd)); + return; + } + // AliInfo(Form("%d %p",id, seed)); + fSeedsPool->RemoveAt(id); + if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 ); + fFreeSeedsID.GetArray()[fNFreeSeeds++] = id; +} + +//________________________________________ +TObject *&AliTPCtrackerMI::NextFreeSeed() +{ + // return next free slot where the seed can be created + fLastSeedID = fNFreeSeeds ? fFreeSeedsID.GetArray()[--fNFreeSeeds] : fSeedsPool->GetEntriesFast(); + // AliInfo(Form("%d",fLastSeedID)); + return (*fSeedsPool)[ fLastSeedID ]; + // +} + +//________________________________________ +void AliTPCtrackerMI::ResetSeedsPool() +{ + // mark all seeds in the pool as unused + AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool->GetSize(),fSeedsPool->GetEntriesFast(),fNFreeSeeds)); + fNFreeSeeds = 0; + fSeedsPool->Clear("C"); // RS: nominally the seeds may allocate memory... +}