X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCtrackerMI.cxx;h=6edb2eb50e2060d93a34f5c196444b6586231d2d;hp=7c24512ea6b635d39e731dd2e9085750b2a2d0fc;hb=4ce766eb4226c070c454ba72807fbf481119a5cb;hpb=56c6de9ece6056f4774e7934abbd20d30465c456 diff --git a/TPC/AliTPCtrackerMI.cxx b/TPC/AliTPCtrackerMI.cxx index 7c24512ea6b..6edb2eb50e2 100644 --- a/TPC/AliTPCtrackerMI.cxx +++ b/TPC/AliTPCtrackerMI.cxx @@ -67,12 +67,42 @@ // // 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 // +// +// Adding systematic errors to the covariance: +// +// 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 default values are 0. +// +// The sytematic errors are added to the covariance matrix in following places: +// +// 1. During fisrt itteration - AliTPCtrackerMI::FillESD +// 2. Second iteration - +// 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds +// 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack +// 3. Third iteration - +// 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds +// 3.b TPC->ITS - AliTPCtrackerMI::RefitInward +// +// 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 - dump information about +// 1.a) cluster +// 2.a) cluster error estimate +// 3.a) cluster shape estimate +// +// +// Debug streamer levels: +// //------------------------------------------------------- @@ -97,23 +127,25 @@ #include "AliTPCReconstructor.h" #include "AliTPCpolyTrack.h" #include "AliTPCreco.h" -#include "AliTPCseed.h" +#include "AliTPCseed.h" + +#include "AliTPCtrackerSector.h" #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 "AliTPCTransform.h" +#include "AliTPCClusterParam.h" // ClassImp(AliTPCtrackerMI) + class AliTPCFastMath { public: AliTPCFastMath(); @@ -164,7 +196,7 @@ AliTPCtrackerMI::AliTPCtrackerMI() fNtracks(0), fSeeds(0), fIteration(0), - fParam(0), + fkParam(0), fDebugStreamer(0) { // @@ -190,7 +222,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; @@ -233,10 +265,10 @@ Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){ Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster()); // - track->SetErrorY2(track->GetErrorY2()*1.3); - track->SetErrorY2(track->GetErrorY2()+0.01); - track->SetErrorZ2(track->GetErrorZ2()*1.3); - track->SetErrorZ2(track->GetErrorZ2()+0.005); +// track->SetErrorY2(track->GetErrorY2()*1.3); +// track->SetErrorY2(track->GetErrorY2()+0.01); +// track->SetErrorZ2(track->GetErrorZ2()*1.3); +// track->SetErrorZ2(track->GetErrorZ2()+0.005); //} if (accept>0) return 0; if (track->GetNumberOfClusters()%20==0){ @@ -252,39 +284,80 @@ Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){ -Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor, - Float_t cory, Float_t corz) +Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster) { // // decide according desired precision to accept given // cluster for tracking - Double_t sy2=ErrY2(seed,cluster)*cory; - Double_t sz2=ErrZ2(seed,cluster)*corz; - //sy2=ErrY2(seed,cluster)*cory; - //sz2=ErrZ2(seed,cluster)*cory; + 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 (rdistance2>16) return 3; + 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); + 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.="<32) return 3; - if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0) + if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0) return 2; //suspisiouce - will be changed - if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0) + if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0) // strict cut on overlaped cluster return 2; //suspisiouce - will be changed - if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor ) + if ( (rdistancey2>1. || rdistancez2>6.25 ) && cluster->GetType()<0){ seed->SetNFoundable(seed->GetNFoundable()-1); return 2; @@ -295,6 +368,7 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste + //_____________________________________________________________________________ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par): AliTracker(), @@ -314,38 +388,40 @@ AliTracker(), fNtracks(0), fSeeds(0), fIteration(0), - fParam(0), + fkParam(0), fDebugStreamer(0) { //--------------------------------------------------------------------- // The main TPC tracker constructor //--------------------------------------------------------------------- - fInnerSec=new AliTPCSector[fkNIS]; - fOuterSec=new AliTPCSector[fkNOS]; + fInnerSec=new AliTPCtrackerSector[fkNIS]; + fOuterSec=new AliTPCtrackerSector[fkNOS]; Int_t i; 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"); + } } //________________________________________________________________________ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): @@ -366,7 +442,7 @@ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): fNtracks(0), fSeeds(0), fIteration(0), - fParam(0), + fkParam(0), fDebugStreamer(0) { //------------------------------------ @@ -374,7 +450,8 @@ AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t): //------------------------------------------------------------------ fOutput=t.fOutput; } -AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){ +AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/) +{ //------------------------------ // dummy //-------------------------------------------------------------- @@ -394,70 +471,8 @@ AliTPCtrackerMI::~AliTPCtrackerMI() { if (fDebugStreamer) delete fDebugStreamer; } -void AliTPCtrackerMI::SetIO() -{ - // - fNewIO = kTRUE; - fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName()); - - fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName()); - if (fOutput){ - AliTPCtrack *iotrack= new AliTPCtrack; - fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100); - delete iotrack; - } -} - -void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESDEvent * event) -{ - - // set input - fNewIO = kFALSE; - fInput = 0; - fOutput = 0; - fSeedTree = 0; - fTreeDebug =0; - fInput = input; - if (input==0){ - return; - } - //set output - fOutput = output; - if (output){ - AliTPCtrack *iotrack= new AliTPCtrack; - // iotrack->fHelixIn = new TClonesArray("AliHelix"); - //iotrack->fHelixOut = new TClonesArray("AliHelix"); - fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100); - delete iotrack; - } - if (output && (fDebug&2)){ - //write the full seed information if specified in debug mode - // - fSeedTree = new TTree("Seeds","Seeds"); - AliTPCseed * vseed = new AliTPCseed; - // - TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160); - arrtr->ExpandCreateFast(160); - TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160); - // - vseed->SetPoints(arrtr); - vseed->SetEPoints(arre); - // vseed->fClusterPoints = arrcl; - fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99); - delete arrtr; - delete arre; - fTreeDebug = new TTree("trackDebug","trackDebug"); - TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0); - fTreeDebug->Branch("debug",&arrd,32000,99); - } - - - //set ESD event - fEvent = event; -} - -void AliTPCtrackerMI::FillESD(TObjArray* arr) +void AliTPCtrackerMI::FillESD(const TObjArray* arr) { // // @@ -471,9 +486,15 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i); if (!pt) continue; pt->UpdatePoints(); - // pt->PropagateTo(fParam->GetInnerRadiusLow()); + AddCovariance(pt); + 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){ @@ -571,361 +592,338 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr) continue; } } + // >> 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) } - printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()); + 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){ // - //seed->SetErrorY2(0.1); - //return 0.1; - //calculate look-up table at the beginning - static Bool_t ginit = kFALSE; - static Float_t gnoise1,gnoise2,gnoise3; - static Float_t ggg1[10000]; - static Float_t ggg2[10000]; - static Float_t ggg3[10000]; - static Float_t glandau1[10000]; - static Float_t glandau2[10000]; - static Float_t glandau3[10000]; // - static Float_t gcor01[500]; - static Float_t gcor02[500]; - static Float_t gcorp[500]; + // Use calibrated cluster error from OCDB // - + AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam(); // - if (ginit==kFALSE){ - for (Int_t i=1;i<500;i++){ - Float_t rsigma = float(i)/100.; - gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6); - gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6); - gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2); - } - - // - for (Int_t i=3;i<10000;i++){ - // - // - // inner sector - Float_t amp = float(i); - Float_t padlength =0.75; - 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; - 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.; - // - // outer short - padlength =1.; - 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; - 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.; - // - // - // outer long - padlength =1.5; - 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; - 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.; - // - } - ginit = kTRUE; - } - // - // - // - Int_t amp = int(TMath::Abs(cl->GetQ())); - if (amp>9999) { - seed->SetErrorY2(1.); - 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(); - angle2 = angle2/(1-angle2); - // - //cluster "quality" - Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2())); - Float_t res; - // - if (fSectors==fInnerSec){ - snoise2 = gnoise1; - res = ggg1[amp]*z+glandau1[amp]*angle2; - if (ctype==0) res *= gcor01[rsigmay]; - if ((ctype>0)){ - res+=0.002; - res*= gcorp[rsigmay]; - } - } - else { - if (padlength<1.1){ - snoise2 = gnoise2; - res = ggg2[amp]*z+glandau2[amp]*angle2; - if (ctype==0) res *= gcor02[rsigmay]; - if ((ctype>0)){ - res+=0.002; - res*= gcorp[rsigmay]; - } - } - else{ - snoise2 = gnoise3; - res = ggg3[amp]*z+glandau3[amp]*angle2; - if (ctype==0) res *= gcor02[rsigmay]; - if ((ctype>0)){ - res+=0.002; - res*= gcorp[rsigmay]; - } - } - } + Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2; + Double_t angle = seed->GetSnp()*seed->GetSnp(); + angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle))); + Double_t erry2 = clparam->GetError0Par(0,type, z,angle); + if (ctype<0) { + erry2+=0.5; // edge cluster + } + erry2*=erry2; + seed->SetErrorY2(erry2); + // + return erry2; + +//calculate look-up table at the beginning +// static Bool_t ginit = kFALSE; +// static Float_t gnoise1,gnoise2,gnoise3; +// static Float_t ggg1[10000]; +// static Float_t ggg2[10000]; +// static Float_t ggg3[10000]; +// static Float_t glandau1[10000]; +// static Float_t glandau2[10000]; +// static Float_t glandau3[10000]; +// // +// static Float_t gcor01[500]; +// static Float_t gcor02[500]; +// static Float_t gcorp[500]; +// // - if (ctype<0){ - res+=0.005; - res*=2.4; // overestimate error 2 times - } - res+= snoise2; +// // +// if (ginit==kFALSE){ +// for (Int_t i=1;i<500;i++){ +// Float_t rsigma = float(i)/100.; +// gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6); +// gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6); +// gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2); +// } + +// // +// for (Int_t i=3;i<10000;i++){ +// // +// // +// // inner sector +// Float_t amp = float(i); +// Float_t padlength =0.75; +// gnoise1 = 0.0004/padlength; +// Float_t nel = 0.268*amp; +// Float_t nprim = 0.155*amp; +// 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.; +// // +// // outer short +// padlength =1.; +// gnoise2 = 0.0004/padlength; +// nel = 0.3*amp; +// nprim = 0.133*amp; +// 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.; +// // +// // +// // outer long +// padlength =1.5; +// gnoise3 = 0.0004/padlength; +// nel = 0.3*amp; +// nprim = 0.133*amp; +// 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.; +// // +// } +// ginit = kTRUE; +// } +// // +// // +// // +// Int_t amp = int(TMath::Abs(cl->GetQ())); +// if (amp>9999) { +// seed->SetErrorY2(1.); +// return 1.; +// } +// Float_t snoise2; +// 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(); +// angle2 = angle2/(1-angle2); +// // +// //cluster "quality" +// Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2())); +// Float_t res; +// // +// if (fSectors==fInnerSec){ +// snoise2 = gnoise1; +// res = ggg1[amp]*z+glandau1[amp]*angle2; +// if (ctype==0) res *= gcor01[rsigmay]; +// if ((ctype>0)){ +// res+=0.002; +// res*= gcorp[rsigmay]; +// } +// } +// else { +// if (padlength<1.1){ +// snoise2 = gnoise2; +// res = ggg2[amp]*z+glandau2[amp]*angle2; +// if (ctype==0) res *= gcor02[rsigmay]; +// if ((ctype>0)){ +// res+=0.002; +// res*= gcorp[rsigmay]; +// } +// } +// else{ +// snoise2 = gnoise3; +// res = ggg3[amp]*z+glandau3[amp]*angle2; +// if (ctype==0) res *= gcor02[rsigmay]; +// if ((ctype>0)){ +// res+=0.002; +// res*= gcorp[rsigmay]; +// } +// } +// } + +// if (ctype<0){ +// res+=0.005; +// res*=2.4; // overestimate error 2 times +// } +// res+= snoise2; - if (res<2*snoise2) - res = 2*snoise2; +// if (res<2*snoise2) +// res = 2*snoise2; - seed->SetErrorY2(res); - return res; +// seed->SetErrorY2(res); +// return res; } -Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ +Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){ // // - //seed->SetErrorY2(0.1); - //return 0.1; - //calculate look-up table at the beginning - static Bool_t ginit = kFALSE; - static Float_t gnoise1,gnoise2,gnoise3; - static Float_t ggg1[10000]; - static Float_t ggg2[10000]; - static Float_t ggg3[10000]; - static Float_t glandau1[10000]; - static Float_t glandau2[10000]; - static Float_t glandau3[10000]; + // Use calibrated cluster error from OCDB // - static Float_t gcor01[1000]; - static Float_t gcor02[1000]; - static Float_t gcorp[1000]; + AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam(); // - - // - if (ginit==kFALSE){ - for (Int_t i=1;i<1000;i++){ - Float_t rsigma = float(i)/100.; - gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6); - gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6); - gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2); - } - - // - for (Int_t i=3;i<10000;i++){ - // - // - // inner sector - Float_t amp = float(i); - Float_t padlength =0.75; - 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; - 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.; - // - // outer short - padlength =1.; - 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; - 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.; - // - // - // outer long - padlength =1.5; - 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; - 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.; - // - } - ginit = kTRUE; - } - // - // - // - Int_t amp = int(TMath::Abs(cl->GetQ())); - if (amp>9999) { - seed->SetErrorY2(1.); - 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()); + Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2; // Double_t angle2 = seed->GetSnp()*seed->GetSnp(); - // if (angle2<0.6) angle2 = 0.6; angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2)); - // - //cluster "quality" - Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2())); - Float_t res; - // - if (fSectors==fInnerSec){ - snoise2 = gnoise1; - res = ggg1[amp]*z+glandau1[amp]*angle2; - if (ctype==0) res *= gcor01[rsigmaz]; - if ((ctype>0)){ - res+=0.002; - res*= gcorp[rsigmaz]; - } - } - else { - if (padlength<1.1){ - snoise2 = gnoise2; - res = ggg2[amp]*z+glandau2[amp]*angle2; - if (ctype==0) res *= gcor02[rsigmaz]; - if ((ctype>0)){ - res+=0.002; - res*= gcorp[rsigmaz]; - } - } - else{ - snoise2 = gnoise3; - res = ggg3[amp]*z+glandau3[amp]*angle2; - if (ctype==0) res *= gcor02[rsigmaz]; - if ((ctype>0)){ - res+=0.002; - res*= gcorp[rsigmaz]; - } - } - } - - if (ctype<0){ - res+=0.002; - res*=1.3; - } - if ((ctype<0) &&<70){ - res+=0.002; - res*=1.3; + Double_t angle = TMath::Sqrt(TMath::Abs(angle2)); + Double_t errz2 = clparam->GetError0Par(1,type, z,angle); + if (ctype<0) { + errz2+=0.5; // edge cluster } - res += snoise2; - if (res<2*snoise2) - res = 2*snoise2; - if (res>3) res =3; - seed->SetErrorZ2(res); - return res; -} - - - -/* -Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){ - // + errz2*=errz2; + seed->SetErrorZ2(errz2); // - //seed->SetErrorZ2(0.1); - //return 0.1; + return errz2; - Float_t snoise2; - Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ())); - // - Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2); - Int_t ctype = cl->GetType(); - Float_t amp = TMath::Abs(cl->GetQ()); - - Float_t nel; - Float_t nprim; - // - Float_t landau=2 ; //landau fluctuation part - Float_t gg=2; // gg fluctuation part - Float_t padlength= GetPadPitchLength(seed->GetX()); - - if (fSectors==fInnerSec){ - snoise2 = 0.0004/padlength; - nel = 0.268*amp; - nprim = 0.155*amp; - gg = (2+0.001*nel/(padlength*padlength))/nel; - landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; - if (landau>1) landau=1; - } - else { - snoise2 = 0.0004/padlength; - nel = 0.3*amp; - nprim = 0.133*amp; - gg = (2+0.0008*nel/(padlength*padlength))/nel; - landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim; - if (landau>1) landau=1; - } - Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z; - - // - Float_t angle2 = seed->GetSnp()*seed->GetSnp(); - angle2 = TMath::Sqrt((1-angle2)); - if (angle2<0.6) angle2 = 0.6; - //angle2 = 1; - - Float_t angle = seed->GetTgl()/angle2; - Float_t angular = landau*angle*angle*padlength*padlength/12.; - Float_t res = sdiff + angular; - - if ((ctype==0) && (fSectors ==fOuterSec)) - res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2)); - if ((ctype==0) && (fSectors ==fInnerSec)) - res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2)); - - if ((ctype>0)){ - res+=0.005; - res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype; - } - if (ctype<0){ - res+=0.002; - res*=1.3; - } - if ((ctype<0) &&<70){ - res+=0.002; - res*=1.3; - } - res += snoise2; - if (res<2*snoise2) - res = 2*snoise2; +// //seed->SetErrorY2(0.1); +// //return 0.1; +// //calculate look-up table at the beginning +// static Bool_t ginit = kFALSE; +// static Float_t gnoise1,gnoise2,gnoise3; +// static Float_t ggg1[10000]; +// static Float_t ggg2[10000]; +// static Float_t ggg3[10000]; +// static Float_t glandau1[10000]; +// static Float_t glandau2[10000]; +// static Float_t glandau3[10000]; +// // +// static Float_t gcor01[1000]; +// static Float_t gcor02[1000]; +// static Float_t gcorp[1000]; +// // - seed->SetErrorZ2(res); - return res; +// // +// if (ginit==kFALSE){ +// for (Int_t i=1;i<1000;i++){ +// Float_t rsigma = float(i)/100.; +// gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6); +// gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6); +// gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2); +// } + +// // +// for (Int_t i=3;i<10000;i++){ +// // +// // +// // inner sector +// Float_t amp = float(i); +// Float_t padlength =0.75; +// gnoise1 = 0.0004/padlength; +// Float_t nel = 0.268*amp; +// Float_t nprim = 0.155*amp; +// 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.; +// // +// // outer short +// padlength =1.; +// gnoise2 = 0.0004/padlength; +// nel = 0.3*amp; +// nprim = 0.133*amp; +// 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.; +// // +// // +// // outer long +// padlength =1.5; +// gnoise3 = 0.0004/padlength; +// nel = 0.3*amp; +// nprim = 0.133*amp; +// 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.; +// // +// } +// ginit = kTRUE; +// } +// // +// // +// // +// Int_t amp = int(TMath::Abs(cl->GetQ())); +// if (amp>9999) { +// seed->SetErrorY2(1.); +// return 1.; +// } +// Float_t snoise2; +// 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(); +// // if (angle2<0.6) angle2 = 0.6; +// angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2)); +// // +// //cluster "quality" +// Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2())); +// Float_t res; +// // +// if (fSectors==fInnerSec){ +// snoise2 = gnoise1; +// res = ggg1[amp]*z+glandau1[amp]*angle2; +// if (ctype==0) res *= gcor01[rsigmaz]; +// if ((ctype>0)){ +// res+=0.002; +// res*= gcorp[rsigmaz]; +// } +// } +// else { +// if (padlength<1.1){ +// snoise2 = gnoise2; +// res = ggg2[amp]*z+glandau2[amp]*angle2; +// if (ctype==0) res *= gcor02[rsigmaz]; +// if ((ctype>0)){ +// res+=0.002; +// res*= gcorp[rsigmaz]; +// } +// } +// else{ +// snoise2 = gnoise3; +// res = ggg3[amp]*z+glandau3[amp]*angle2; +// if (ctype==0) res *= gcor02[rsigmaz]; +// if ((ctype>0)){ +// res+=0.002; +// res*= gcorp[rsigmaz]; +// } +// } +// } + +// if (ctype<0){ +// res+=0.002; +// res*=1.3; +// } +// if ((ctype<0) &&<70){ +// res+=0.002; +// res*=1.3; +// } +// res += snoise2; +// if (res<2*snoise2) +// res = 2*snoise2; +// if (res>3) res =3; +// seed->SetErrorZ2(res); +// return res; } -*/ + @@ -954,7 +952,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 @@ -975,7 +973,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 @@ -986,7 +984,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; } // @@ -1001,7 +999,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 @@ -1012,7 +1010,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; } // @@ -1031,7 +1029,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 @@ -1050,7 +1048,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 @@ -1061,7 +1059,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 @@ -1081,7 +1079,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. //----------------------------------------------------------------- @@ -1092,8 +1090,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]; @@ -1116,7 +1114,7 @@ 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) { // // @@ -1124,14 +1122,123 @@ Int_t AliTPCtrackerMI::LoadClusters (TTree *tree) return LoadClusters(); } + +Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr) +{ + // + // load clusters to the memory + 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; + + // + Int_t sec,row; + fkParam->AdjustSectorRow(clrow->GetID(),sec,row); + + for (Int_t icl=0; iclGetArray()->GetEntriesFast(); icl++){ + Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl))); + } + // + if (clrow->GetArray()->GetEntriesFast()<=0) continue; + AliTPCtrackerRow * tpcrow=0; + Int_t left=0; + if (secSetN1(clrow->GetArray()->GetEntriesFast()); + tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]); + 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 j=0;jGetN2();++j) + tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j))); + } + clrow->GetArray()->Clear("C"); + } + // + delete clrow; + LoadOuterSectors(); + LoadInnerSectors(); + 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); + AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI"); // // TTree * tree = fClustersArray.GetTree(); @@ -1144,12 +1251,12 @@ Int_t AliTPCtrackerMI::LoadClusters() br->GetEntry(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))); } // - AliTPCRow * tpcrow=0; + AliTPCtrackerRow * tpcrow=0; Int_t left=0; if (secSetN1(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))); } } // @@ -1188,7 +1295,7 @@ void AliTPCtrackerMI::UnloadClusters() Int_t nrows = fOuterSec->GetNRows(); for (Int_t sec = 0;secfClusters1) delete []tpcrow->fClusters1; // if (tpcrow->fClusters2) delete []tpcrow->fClusters2; @@ -1199,7 +1306,7 @@ void AliTPCtrackerMI::UnloadClusters() nrows = fInnerSec->GetNRows(); for (Int_t sec = 0;secfClusters1) delete []tpcrow->fClusters1; //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; @@ -1210,6 +1317,32 @@ void AliTPCtrackerMI::UnloadClusters() return ; } +void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{ + // + // Filling cluster to the array - For visualization purposes + // + Int_t nrows=0; + nrows = fOuterSec->GetNRows(); + for (Int_t sec = 0;secGetN();icl++){ + array->AddLast((TObject*)((*tpcrow)[icl])); + } + } + nrows = fInnerSec->GetNRows(); + for (Int_t sec = 0;secGetN();icl++){ + array->AddLast((TObject*)(*tpcrow)[icl]); + } + } +} + + void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){ // // @@ -1218,24 +1351,30 @@ void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){ if (!transform) { AliFatal("Tranformations not in calibDB"); } + 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(); TTreeSRedirector &cstream = *fDebugStreamer; cstream<<"Transform"<< + "event="<IsGeoRead()) fParam->ReadGeoMatrices(); - TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector()); + //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices(); + 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()}; @@ -1270,7 +1409,7 @@ Int_t AliTPCtrackerMI::LoadOuterSectors() { UInt_t index=0; for (Int_t sec = 0;secGetN1(); @@ -1318,7 +1457,7 @@ Int_t AliTPCtrackerMI::LoadInnerSectors() { UInt_t index=0; for (Int_t sec = 0;secGetN1(); @@ -1370,7 +1509,7 @@ AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const { Int_t row=(index&0x00ff0000)>>16; Int_t ncl=(index&0x00007fff)>>00; - const AliTPCRow * tpcrow=0; + const AliTPCtrackerRow * tpcrow=0; AliTPCclusterMI * clrow =0; if (sec<0 || sec>=fkNIS*4) { @@ -1379,7 +1518,9 @@ AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const { } if (secAliTPCReconstructor::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; @@ -1507,7 +1652,7 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { if (!isActive || !isActive2) return 0; - const AliTPCRow &krow=GetRow(t.GetRelativeSector(),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.; @@ -1519,7 +1664,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())IsUsed(10)) return 0; - Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.); + Int_t accept = AcceptCluster(&t,t.GetCurrentCluster()); if (fIteration==2&&cl->IsUsed(11)) { t.SetErrorY2(t.GetErrorY2()+0.03); t.SetErrorZ2(t.GetErrorZ2()+0.03); @@ -1562,83 +1708,6 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) { return 1; } -Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) { - //----------------------------------------------------------------- - // This function tries to find a track prolongation to next pad row - //----------------------------------------------------------------- - // - Double_t x= GetXrow(nr), ymax=GetMaxY(nr); - Double_t y,z; - if (!t.GetProlongation(x,y,z)) { - t.SetRemoval(10); - return 0; - } - // - // - if (TMath::Abs(y)>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 3 pad-row - if ( (nr%6==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 AliTPCRow &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,1.); - //if (accept<3){ - t.SetClusterIndex2(row,index); - t.SetClusterPointer(row, cl); - //} - } - return 1; -} - //_________________________________________________________________________ @@ -1651,19 +1720,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.; @@ -1673,13 +1742,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); @@ -1746,7 +1815,7 @@ Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) { } //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha())); - AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr); + AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr); if (TMath::Abs(TMath::Abs(y)-ymax)IsUsed(10)){ // @@ -1892,27 +1963,6 @@ 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; -} - @@ -2154,8 +2204,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(); @@ -2212,11 +2265,25 @@ 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 + 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 + if( AliTPCReconstructor::StreamLevel()>3){ + TTreeSRedirector &cstream = *fDebugStreamer; + cstream<<"RemoveShort"<< + "iter="<RemoveAt(trackindex); continue; } @@ -2245,6 +2312,88 @@ 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++){ + AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1(); + for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){ + Float_t gx[3]; cl[icl].GetGlobalXYZ(gx); + (*fDebugStreamer)<<"clDump"<< + "iter="<GetNRows();row++){ + AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1(); + for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){ + Float_t gx[3]; cl[icl].GetGlobalXYZ(gx); + (*fDebugStreamer)<<"clDump"<< + "iter="<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); @@ -2428,11 +2577,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); @@ -2447,7 +2596,7 @@ 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 +void AliTPCtrackerMI::StopNotActive(const 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 @@ -2522,6 +2671,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(); @@ -2529,11 +2679,31 @@ 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); seed->CookdEdx(0.02,0.6); CookLabel(seed,0.1); //For comparison only esd->SetTPCClusterMap(seed->GetClusterMap()); @@ -2551,8 +2721,8 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event) 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(); esd->SetTPCsignal(dedx, sdedx, ndedx); // @@ -2570,8 +2740,8 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event) } } //FindKinks(fSeeds,event); + if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds); Info("RefitInward","Number of refitted tracks %d",ntracks); - fEvent =0; return 0; } @@ -2589,6 +2759,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(); @@ -2598,28 +2769,49 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event) if (!seed) continue; if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks seed->UpdatePoints(); + AddCovariance(seed); 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)<<"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; @@ -2643,7 +2835,7 @@ void AliTPCtrackerMI::DeleteSeeds() fSeeds =0; } -void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction) +void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction) { // //read seeds from the event @@ -2668,6 +2860,7 @@ void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction) // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha()); AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/); seed->SetUniqueID(esd->GetID()); + AddCovariance(seed); //add systematic ucertainty for (Int_t ikink=0;ikink<3;ikink++) { Int_t index = esd->GetKinkIndex(ikink); seed->GetKinkIndexes()[ikink] = index; @@ -2686,11 +2879,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); + // delete seed; + // continue; + //} if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) { Double_t par0[5],par1[5],alpha,x; esd->GetInnerExternalParameters(alpha,x,par0); @@ -2783,11 +2976,11 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Int_t imiddle = (i2+i1)/2; //middle pad row index Double_t xm = GetXrow(imiddle); // radius of middle pad-row - const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row + const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row // Int_t ns =sec; - const AliTPCRow& kr1=GetRow(ns,i1); + const AliTPCtrackerRow& kr1=GetRow(ns,i1); Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5; Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5; @@ -2839,10 +3032,10 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){ Int_t sec2 = sec + dsec; // - // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2]; - //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle]; - AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2); - AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle); + // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2]; + //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle]; + AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2); + AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle); Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0); Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2); @@ -3095,25 +3288,25 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, // first 3 padrows Double_t x1 = GetXrow(i1-1); - const AliTPCRow& kr1=GetRow(sec,i1-1); + const AliTPCtrackerRow& kr1=GetRow(sec,i1-1); Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5; // Double_t x1p = GetXrow(i1); - const AliTPCRow& kr1p=GetRow(sec,i1); + const AliTPCtrackerRow& kr1p=GetRow(sec,i1); // Double_t x1m = GetXrow(i1-2); - const AliTPCRow& kr1m=GetRow(sec,i1-2); + const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2); // //last 3 padrow for seeding - AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7); + AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7); Double_t x3 = GetXrow(i1-7); // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5; // - AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6); + AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6); Double_t x3p = GetXrow(i1-6); // - AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8); + AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8); Double_t x3m = GetXrow(i1-8); // @@ -3121,7 +3314,7 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, // middle padrow Int_t im = i1-4; //middle pad row index Double_t xm = GetXrow(im); // radius of middle pad-row - const AliTPCRow& krm=GetRow(sec,im); //middle pad -row + const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5; // // @@ -3267,7 +3460,7 @@ 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); + index=kr1.GetIndex(is); seed->~AliTPCseed(); AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index); @@ -3316,7 +3509,7 @@ 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; } @@ -3345,8 +3538,8 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, // Double_t cs=cos(alpha), sn=sin(alpha); Int_t row0 = (i1+i2)/2; Int_t drow = (i1-i2)/2; - const AliTPCRow& kr0=fSectors[sec][row0]; - AliTPCRow * kr=0; + const AliTPCtrackerRow& kr0=fSectors[sec][row0]; + AliTPCtrackerRow * kr=0; AliTPCpolyTrack polytrack; Int_t nclusters=fSectors[sec][row0]; @@ -3359,8 +3552,8 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Int_t nfound =0; Int_t nfoundable =0; for (Int_t iter =1; iter<2; iter++){ //iterations - const AliTPCRow& krm=fSectors[sec][row0-iter]; - const AliTPCRow& krp=fSectors[sec][row0+iter]; + const AliTPCtrackerRow& krm=fSectors[sec][row0-iter]; + const AliTPCtrackerRow& krp=fSectors[sec][row0+iter]; const AliTPCclusterMI * cl= kr0[is]; if (cl->IsUsed(10)) { @@ -3429,9 +3622,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) { @@ -3580,7 +3773,7 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, } -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) { // // @@ -3712,7 +3905,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) { // // @@ -3979,7 +4172,7 @@ AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward) -void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter) +void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter) { // // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES @@ -4048,7 +4241,6 @@ void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_ } if (ncl>0) xm[i]/=Float_t(ncl); } - TTreeSRedirector &cstream = *fDebugStreamer; // for (Int_t i0=0;i0At(i0); @@ -4109,6 +4301,8 @@ void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_ } } // + if (AliTPCReconstructor::StreamLevel()>5) { + TTreeSRedirector &cstream = *fDebugStreamer; cstream<<"Multi"<< "iter="<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++; + // 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; + // + Int_t sumShared=0; + Int_t nall0=0; + Int_t nall1=0; + Int_t firstShared=lastpoint, lastShared=firstpoint; + Int_t firstRow=lastpoint, lastRow=firstpoint; + // + 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++; + } } - } - 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; - // - Float_t xc1 = helixes[i1].GetHelix(6); - Float_t yc1 = helixes[i1].GetHelix(7); - Float_t fi1 = TMath::ATan2(yc1,xc1); - // - Float_t dfi = fi0-fi1; - if (dfi>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; + 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); + delete 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 @@ -4416,7 +4569,6 @@ void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_ // // Find tracks // - TTreeSRedirector &cstream = *fDebugStreamer; // for (Int_t i0=0;i0At(i0); @@ -4509,11 +4661,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++; @@ -4962,12 +5139,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); // @@ -4975,10 +5152,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; @@ -4993,9 +5170,9 @@ 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; @@ -5006,7 +5183,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) if (seed1) delete 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; @@ -5017,7 +5194,7 @@ void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd) delete 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); // @@ -5092,23 +5269,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 @@ -5149,15 +5326,15 @@ 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) { + if (Float_t(ishared+1)/Float_t(all+1)>0.5) { delete array->RemoveAt(i); continue; } @@ -5171,8 +5348,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; @@ -5185,15 +5362,15 @@ 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) { @@ -5240,11 +5417,11 @@ 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::FindV0s(const TObjArray * array, AliESDEvent *const esd) { // // find V0s @@ -5319,7 +5496,6 @@ void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd) // // // // - TTreeSRedirector &cstream = *fDebugStreamer; Float_t fprimvertex[3]={GetX(),GetY(),GetZ()}; AliV0 vertex; Double_t cradius0 = 10*10; @@ -5334,6 +5510,7 @@ void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd) AliTPCseed * track0 = (AliTPCseed*)array->At(i); if (!track0) continue; if (AliTPCReconstructor::StreamLevel()>1){ + TTreeSRedirector &cstream = *fDebugStreamer; cstream<<"Tracks"<< "Tr0.="<GetKinkAngleCutChi2(0)){ + delete seed0; + delete seed1; + return 0; + } + + row0 = GetRowNumber(kink.GetR()); kink.SetTPCRow0(row0); kink.SetLabel(CookLabel(seed0,0.5,0,row0),0); @@ -5888,9 +6092,12 @@ 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)); + if ( chi2P2+chi2P3GetKinkAngleCutChi2(1)){ + mother=*seed; + } delete seed0; delete seed1; // @@ -5965,7 +6172,7 @@ Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) { return 0; } -Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd) +Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd) { // if (fSeeds) DeleteSeeds(); @@ -5973,6 +6180,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd) Clusters2Tracks(); if (!fSeeds) return 1; FillESD(fSeeds); + if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds); return 0; // } @@ -6015,9 +6223,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 @@ -6187,6 +6396,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]; @@ -6200,7 +6412,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; @@ -6225,7 +6437,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; @@ -6278,9 +6490,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; @@ -6303,10 +6528,9 @@ TObjArray * AliTPCtrackerMI::Tracking() fdensity = 2.; cuts[0]=0.0080; - Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec(); // find secondaries - for (Int_t delta = 30; deltaGetNRowLow()>rfirst+1) ) + if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) ) FollowProlongation(t, rfirst+1); } @@ -6452,7 +6676,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; } @@ -6463,7 +6687,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; } @@ -6472,7 +6696,7 @@ void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rla } } -void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const +void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const { // // @@ -6508,7 +6732,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 { // // @@ -6527,7 +6751,7 @@ void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const } -Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) +Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr) { // // make back propagation @@ -6558,7 +6782,7 @@ Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr) } -Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr) +Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr) { // // make forward propagation @@ -6606,7 +6830,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 @@ -6641,19 +6865,27 @@ Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1) void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row) { // - // - 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 padlength = GetPadPitchLength(row); - // - Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3; + // + AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam(); + 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); - seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy); - // - Float_t sresz = fParam->GetZSigma(); - Float_t angularz = seed->GetTgl(); - seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz); + 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((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() < fkParam->GetNSector()/2) ? 0.2 :0.3; +// seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy); +// // +// Float_t sresz = fkParam->GetZSigma(); +// seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz); /* Float_t wy = GetSigmaY(seed); Float_t wz = GetSigmaZ(seed); @@ -6666,32 +6898,6 @@ void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row) } -Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed) -{ - // - // - Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); - Float_t padlength = fParam->GetPadPitchLength(seed->GetSector()); - Float_t sres = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3; - Float_t angular = seed->GetSnp(); - angular = angular*angular/(1-angular*angular); - // angular*=angular; - //angular = TMath::Sqrt(angular/(1-angular)); - Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres); - return res; -} -Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed) -{ - // - // - Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL(); - Float_t padlength = fParam->GetPadPitchLength(seed->GetSector()); - Float_t sres = fParam->GetZSigma(); - Float_t angular = seed->GetTgl(); - Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres); - return res; -} - //__________________________________________________________________________ void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const { @@ -6755,8 +6961,8 @@ 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); @@ -6783,7 +6989,7 @@ void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const { //__________________________________________________________________________ -Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const { +Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const { //-------------------------------------------------------------------- //This function "cooks" a track label. If label<0, this track is fake. //-------------------------------------------------------------------- @@ -6839,8 +7045,8 @@ Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t 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 -1;} + if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab; else { Int_t tail=Int_t(0.10*noc); @@ -6866,29 +7072,6 @@ Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t } -Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const -{ - //return pad row number for this x - Double_t r; - if (fN < 64){ - r=fRow[fN-1].GetX(); - if (x > r) return fN; - r=fRow[0].GetX(); - if (x < r) return -1; - return Int_t((x-r)/fPadPitchLength + 0.5);} - else{ - r=fRow[fN-1].GetX(); - if (x > r) return fN; - r=fRow[0].GetX(); - if (x < r) return -1; - Double_t r1=fRow[64].GetX(); - if(xGetInnerAngle(); - fAlphaShift=par->GetInnerAngleShift(); - fPadPitchWidth=par->GetInnerPadPitchWidth(); - fPadPitchLength=par->GetInnerPadPitchLength(); - fN=par->GetNRowLow(); - if(fRow)delete [] fRow;fRow = 0; - fRow=new AliTPCRow[fN]; - for (Int_t i=0; iGetPadRowRadiiLow(i)); - fRow[i].SetDeadZone(1.5); //1.5 cm of dead zone - } - } else { - fAlpha=par->GetOuterAngle(); - fAlphaShift=par->GetOuterAngleShift(); - fPadPitchWidth = par->GetOuterPadPitchWidth(); - fPadPitchLength = par->GetOuter1PadPitchLength(); - f1PadPitchLength = par->GetOuter1PadPitchLength(); - f2PadPitchLength = par->GetOuter2PadPitchLength(); - fN=par->GetNRowUp(); - if(fRow)delete [] fRow;fRow = 0; - fRow=new AliTPCRow[fN]; - for (Int_t i=0; iGetPadRowRadiiUp(i)); - fRow[i].SetDeadZone(1.5); // 1.5 cm of dead zone - } - } -} - -AliTPCtrackerMI::AliTPCRow::AliTPCRow(): - fDeadZone(0.), - fClusters1(0), - fN1(0), - fClusters2(0), - fN2(0), - fN(0), - fX(0.) -{ - // - // default constructor - // -} - -AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){ - // - for (Int_t i = 0; i < fN1; i++) - fClusters1[i].~AliTPCclusterMI(); - delete [] fClusters1; - for (Int_t i = 0; i < fN2; i++) - fClusters2[i].~AliTPCclusterMI(); - delete [] fClusters2; -} - - - -//_________________________________________________________________________ -void -AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) { - //----------------------------------------------------------------------- - // Insert a cluster into this pad row in accordence with its y-coordinate - //----------------------------------------------------------------------- - if (fN==kMaxClusterPerRow) { - //AliInfo("AliTPCRow::InsertCluster(): Too many clusters"); - return; - } - if (fN>=fN1+fN2) { - //AliInfo("AliTPCRow::InsertCluster(): Too many clusters !"); - } - - if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;} - Int_t i=Find(c->GetZ()); - memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*)); - memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t)); - fIndex[i]=index; fClusters[i]=c; fN++; -} - -void AliTPCtrackerMI::AliTPCRow::ResetClusters() { - // - // reset clusters - // MvL: Need to call destructors for AliTPCclusterMI, to delete fInfo - for (Int_t i = 0; i < fN1; i++) - fClusters1[i].~AliTPCclusterMI(); - delete [] fClusters1; fClusters1=0; - for (Int_t i = 0; i < fN2; i++) - fClusters2[i].~AliTPCclusterMI(); - delete [] fClusters2; fClusters2=0; - - fN = 0; - fN1 = 0; - fN2 = 0; - //delete[] fClusterArray; - - //fClusterArray=0; -} - - -//___________________________________________________________________ -Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const { - //----------------------------------------------------------------------- - // Return the index of the nearest cluster - //----------------------------------------------------------------------- - if (fN==0) return 0; - if (z <= fClusters[0]->GetZ()) return 0; - if (z > fClusters[fN-1]->GetZ()) return fN; - Int_t b=0, e=fN-1, m=(b+e)/2; - for (; b fClusters[m]->GetZ()) b=m+1; - else e=m; - } - return m; -} - - - -//___________________________________________________________________ -AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const { - //----------------------------------------------------------------------- - // Return the index of the nearest cluster in z y - //----------------------------------------------------------------------- - Float_t maxdistance = roady*roady + roadz*roadz; - - AliTPCclusterMI *cl =0; - for (Int_t i=Find(z-roadz); iGetZ() > z+roadz) break; - if ( (c->GetY()-y) > roady ) continue; - Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y); - if (maxdistance>distance) { - maxdistance = distance; - cl=c; - } - } - return cl; -} - -AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const -{ - //----------------------------------------------------------------------- - // Return the index of the nearest cluster in z y - //----------------------------------------------------------------------- - Float_t maxdistance = roady*roady + roadz*roadz; - AliTPCclusterMI *cl =0; - - //PH Check boundaries. 510 is the size of fFastCluster - Int_t iz1 = Int_t(z-roadz+254.5); - if (iz1<0 || iz1>=510) return cl; - iz1 = TMath::Max(GetFastCluster(iz1)-1,0); - Int_t iz2 = Int_t(z+roadz+255.5); - if (iz2<0 || iz2>=510) return cl; - iz2 = TMath::Min(GetFastCluster(iz2)+1,fN); - Bool_t skipUsed = !(AliTPCReconstructor::GetRecoParam()->GetClusterSharing()); - //FindNearest3(y,z,roady,roadz,index); - // for (Int_t i=Find(z-roadz); iGetZ() > z+roadz) break; - if ( c->GetY()-y > roady ) continue; - if ( y-c->GetY() > roady ) continue; - if (skipUsed && c->IsUsed(11)) continue; - Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y); - if (maxdistance>distance) { - maxdistance = distance; - cl=c; - index =i; - //roady = TMath::Sqrt(maxdistance); - } - } - return cl; -} - - -void AliTPCtrackerMI::AliTPCRow::SetFastCluster(Int_t i, Short_t cl){ - // - // Set cluster info for fast navigation - // - if (i>510|| i<0){ - }else{ - fFastCluster[i]=cl; - } -} - - void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t) @@ -7099,11 +7095,14 @@ void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t) Int_t firstpoint = 0; Int_t lastpoint = 159; AliTPCTrackerPoint *point; + AliTPCclusterMI *cluster; for (int iter=firstpoint; iterGetTrackPoint(iter); - if (point) { + // Change to cluster pointers to see if we have a cluster at given padrow + cluster = t->GetClusterPointer(iter); + if (cluster) { t->SetClusterMapBit(iter, kTRUE); + point = t->GetTrackPoint(iter); if (point->IsShared()) t->SetSharedMapBit(iter,kTRUE); else @@ -7115,3 +7114,56 @@ void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t) } } } + +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(); + Double_t *covarIn= (Double_t*)seed->GetCovariance(); + Double_t covar[15]; + for (Int_t i=0;i<15;i++) covar[i]=0; + // 0 + // 1 2 + // 3 4 5 + // 6 7 8 9 + // 10 11 12 13 14 + covar[0] = param[0]*param[0]; + 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[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); +}