#include "Riostream.h"
#include <TClonesArray.h>
#include <TFile.h>
#include <TObjArray.h>
#include <TTree.h>
#include <TMatrixD.h>
#include <TGraphErrors.h>
#include <TTimeStamp.h>
#include "AliLog.h"
#include "AliComplexCluster.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliKink.h"
#include "AliV0.h"
#include "AliHelix.h"
#include "AliRunLoader.h"
#include "AliTPCClustersRow.h"
#include "AliTPCParam.h"
#include "AliTPCReconstructor.h"
#include "AliTPCpolyTrack.h"
#include "AliTPCreco.h"
#include "AliTPCseed.h"
#include "AliTPCtrackerSector.h"
#include "AliTPCtracker.h"
#include "TStopwatch.h"
#include "AliTPCReconstructor.h"
#include "AliAlignObj.h"
#include "AliTrackPointArray.h"
#include "TRandom.h"
#include "AliTPCcalibDB.h"
#include "AliTPCcalibDButil.h"
#include "AliTPCTransform.h"
#include "AliTPCClusterParam.h"
#include "AliTPCdEdxInfo.h"
#include "AliDCSSensorArray.h"
#include "AliDCSSensor.h"
#include "AliDAQ.h"
#include "AliCosmicTracker.h"
#include "AliTPCROC.h"
#include "AliMathBase.h"
using std::cerr;
using std::endl;
ClassImp(AliTPCtracker)
class AliTPCFastMath {
public:
AliTPCFastMath();
static Double_t FastAsin(Double_t x);
private:
static Double_t fgFastAsin[20000];
};
Double_t AliTPCFastMath::fgFastAsin[20000];
AliTPCFastMath gAliTPCFastMath;
AliTPCFastMath::AliTPCFastMath(){
for (Int_t i=0;i<10000;i++){
fgFastAsin[2*i] = TMath::ASin(i/10000.);
fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
}
}
Double_t AliTPCFastMath::FastAsin(Double_t x){
if (x>0){
Int_t index = int(x*10000);
return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
}
x*=-1;
Int_t index = int(x*10000);
return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
}
AliTPCtracker::AliTPCtracker()
:AliTracker(),
fkNIS(0),
fInnerSec(0),
fkNOS(0),
fOuterSec(0),
fN(0),
fSectors(0),
fInput(0),
fOutput(0),
fSeedTree(0),
fTreeDebug(0),
fEvent(0),
fEventHLT(0),
fDebug(0),
fNewIO(kFALSE),
fNtracks(0),
fSeeds(0),
fIteration(0),
fkParam(0),
fDebugStreamer(0),
fUseHLTClusters(4),
fCrossTalkSignalArray(0),
fSeedsPool(0),
fFreeSeedsID(500),
fNFreeSeeds(0),
fLastSeedID(-1)
{
for (Int_t irow=0; irow<200; irow++){
fXRow[irow]=0;
fYMax[irow]=0;
fPadLength[irow]=0;
}
}
Int_t AliTPCtracker::UpdateTrack(AliTPCseed * track, Int_t accept){
AliTPCclusterMI* c =track->GetCurrentCluster();
if (accept > 0)
track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
else
track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
UInt_t i = track->GetCurrentClusterIndex1();
Int_t sec=(i&0xff000000)>>24;
track->SetRow((i&0x00ff0000)>>16);
track->SetSector(sec);
if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
track->SetClusterIndex2(track->GetRow(), i);
if (track->GetFirstPoint()>track->GetRow())
track->SetFirstPoint(track->GetRow());
if (track->GetLastPoint()<track->GetRow())
track->SetLastPoint(track->GetRow());
track->SetClusterPointer(track->GetRow(),c);
Double_t angle2 = track->GetSnp()*track->GetSnp();
if (angle2<1)
{
angle2 = TMath::Sqrt(angle2/(1-angle2));
AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
point.SetErrY(sqrt(track->GetErrorY2()));
point.SetErrZ(sqrt(track->GetErrorZ2()));
point.SetX(track->GetX());
point.SetY(track->GetY());
point.SetZ(track->GetZ());
point.SetAngleY(angle2);
point.SetAngleZ(track->GetTgl());
if (point.IsShared()){
track->SetErrorY2(track->GetErrorY2()*4);
track->SetErrorZ2(track->GetErrorZ2()*4);
}
}
Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
if (accept>0) return 0;
if (track->GetNumberOfClusters()%20==0){
}
if (AliTPCReconstructor::StreamLevel()&kStreamUpdateTrack) {
Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
AliExternalTrackParam param(*track);
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"UpdateTrack"<<
"cl.="<<c<<
"event="<<event<<
"track.="<<¶m<<
"\n";
}
track->SetNoCluster(0);
return track->Update(c,chi2,i);
}
Int_t AliTPCtracker::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
{
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 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;
if (AliTPCReconstructor::StreamLevel()>2 && ( (fIteration>0)|| (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];
Int_t nclSeed=seed->GetNumberOfClusters();
if (AliTPCReconstructor::StreamLevel()&kStreamErrParam) {
Int_t eventNr = fEvent->GetEventNumberInFile();
(*fDebugStreamer)<<"ErrParam"<<
"iter="<<fIteration<<
"eventNr="<<eventNr<<
"Cl.="<<cluster<<
"nclSeed="<<nclSeed<<
"T.="<<¶m<<
"dy="<<dy<<
"dz="<<dz<<
"yt="<<yt<<
"zt="<<zt<<
"gcl.="<<&gcl<<
"gtr.="<<>r<<
"erry2="<<sy2<<
"errz2="<<sz2<<
"rmsy2="<<rmsy2<<
"rmsz2="<<rmsz2<<
"rmsy2p30="<<rmsy2p30<<
"rmsz2p30="<<rmsz2p30<<
"rmsy2p30R="<<rmsy2p30R<<
"rmsz2p30R="<<rmsz2p30R<<
"rdisty="<<rdistancey2<<
"rdistz="<<rdistancez2<<
"rdist="<<rdistance2<<
"\n";
}
}
if (rdistance2>32) return 3;
if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
return 2;
if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
return 2;
if ( (rdistancey2>1. || rdistancez2>6.25 )
&& cluster->GetType()<0){
seed->SetNFoundable(seed->GetNFoundable()-1);
return 2;
}
if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
if (fIteration==2){
if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
return 2;
}
}
}
return 0;
}
AliTPCtracker::AliTPCtracker(const AliTPCParam *par):
AliTracker(),
fkNIS(par->GetNInnerSector()/2),
fInnerSec(0),
fkNOS(par->GetNOuterSector()/2),
fOuterSec(0),
fN(0),
fSectors(0),
fInput(0),
fOutput(0),
fSeedTree(0),
fTreeDebug(0),
fEvent(0),
fEventHLT(0),
fDebug(0),
fNewIO(0),
fNtracks(0),
fSeeds(0),
fIteration(0),
fkParam(0),
fDebugStreamer(0),
fUseHLTClusters(4),
fCrossTalkSignalArray(0),
fSeedsPool(0),
fFreeSeedsID(500),
fNFreeSeeds(0),
fLastSeedID(-1)
{
fInnerSec=new AliTPCtrackerSector[fkNIS];
fOuterSec=new AliTPCtrackerSector[fkNOS];
Int_t i;
for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
fkParam = par;
Int_t nrowlow = par->GetNRowLow();
Int_t nrowup = par->GetNRowUp();
for (i=0;i<nrowlow;i++){
fXRow[i] = par->GetPadRowRadiiLow(i);
fPadLength[i]= par->GetPadPitchLength(0,i);
fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
}
for (i=0;i<nrowup;i++){
fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
}
if (AliTPCReconstructor::StreamLevel()>0) {
fDebugStreamer = new TTreeSRedirector("TPCdebug.root","recreate");
}
fSeedsPool = new TClonesArray("AliTPCseed",1000);
Int_t nROCs = 72;
Int_t nTimeBinsAll = par->GetMaxTBin();
Int_t nWireSegments = 11;
fCrossTalkSignalArray = new TObjArray(nROCs*4);
fCrossTalkSignalArray->SetOwner(kTRUE);
for (Int_t isector=0; isector<4*nROCs; isector++){
TMatrixD * crossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll);
for (Int_t imatrix = 0; imatrix<11; imatrix++)
for (Int_t jmatrix = 0; jmatrix<nTimeBinsAll; jmatrix++){
(*crossTalkSignal)[imatrix][jmatrix]=0.;
}
fCrossTalkSignalArray->AddAt(crossTalkSignal,isector);
}
}
AliTPCtracker::AliTPCtracker(const AliTPCtracker &t):
AliTracker(t),
fkNIS(t.fkNIS),
fInnerSec(0),
fkNOS(t.fkNOS),
fOuterSec(0),
fN(0),
fSectors(0),
fInput(0),
fOutput(0),
fSeedTree(0),
fTreeDebug(0),
fEvent(0),
fEventHLT(0),
fDebug(0),
fNewIO(kFALSE),
fNtracks(0),
fSeeds(0),
fIteration(0),
fkParam(0),
fDebugStreamer(0),
fUseHLTClusters(4),
fCrossTalkSignalArray(0),
fSeedsPool(0),
fFreeSeedsID(500),
fNFreeSeeds(0),
fLastSeedID(-1)
{
fOutput=t.fOutput;
for (Int_t irow=0; irow<200; irow++){
fXRow[irow]=0;
fYMax[irow]=0;
fPadLength[irow]=0;
}
}
AliTPCtracker & AliTPCtracker::operator=(const AliTPCtracker& )
{
return *this;
}
AliTPCtracker::~AliTPCtracker() {
delete[] fInnerSec;
delete[] fOuterSec;
if (fSeeds) {
fSeeds->Clear();
delete fSeeds;
}
if (fCrossTalkSignalArray) delete fCrossTalkSignalArray;
if (fDebugStreamer) delete fDebugStreamer;
if (fSeedsPool) delete fSeedsPool;
}
void AliTPCtracker::FillESD(const TObjArray* arr)
{
if (!fEvent) return;
AliESDtrack iotrack;
Int_t nseed=arr->GetEntriesFast();
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) continue;
pt->UpdatePoints();
AddCovariance(pt);
if (AliTPCReconstructor::StreamLevel()&kStreamFillESD) {
(*fDebugStreamer)<<"FillESD"<<
"Tr0.="<<pt<<
"\n";
}
if (pt->GetKinkIndex(0)<=0){
pt->PropagateTo(fkParam->GetInnerRadiusLow());
}
if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
iotrack.~AliESDtrack();
new(&iotrack) AliESDtrack;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
iotrack.SetV0Indexes(pt->GetV0Indexes());
MakeESDBitmaps(pt, &iotrack);
fEvent->AddTrack(&iotrack);
continue;
}
if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
iotrack.~AliESDtrack();
new(&iotrack) AliESDtrack;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
iotrack.SetV0Indexes(pt->GetV0Indexes());
MakeESDBitmaps(pt, &iotrack);
fEvent->AddTrack(&iotrack);
continue;
}
if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
Int_t found,foundable,shared;
pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
iotrack.~AliESDtrack();
new(&iotrack) AliESDtrack;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
iotrack.SetV0Indexes(pt->GetV0Indexes());
MakeESDBitmaps(pt, &iotrack);
fEvent->AddTrack(&iotrack);
continue;
}
}
if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
Int_t found,foundable,shared;
pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
if (found<20) continue;
if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
iotrack.~AliESDtrack();
new(&iotrack) AliESDtrack;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
iotrack.SetV0Indexes(pt->GetV0Indexes());
MakeESDBitmaps(pt, &iotrack);
fEvent->AddTrack(&iotrack);
continue;
}
if ( (pt->GetNumberOfClusters()>30) ) {
Int_t found,foundable,shared;
pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
iotrack.~AliESDtrack();
new(&iotrack) AliESDtrack;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
iotrack.SetV0Indexes(pt->GetV0Indexes());
MakeESDBitmaps(pt, &iotrack);
fEvent->AddTrack(&iotrack);
continue;
}
}
if ( (pt->GetNumberOfClusters()>15)) {
Int_t found,foundable,shared;
pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
if (found<15) continue;
if (foundable<=0) continue;
if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
if (float(found)/float(foundable)<0.8) continue;
iotrack.~AliESDtrack();
new(&iotrack) AliESDtrack;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCsignal(pt->GetdEdx(), pt->GetSDEDX(0), pt->GetNCDEDX(0));
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
iotrack.SetV0Indexes(pt->GetV0Indexes());
MakeESDBitmaps(pt, &iotrack);
fEvent->AddTrack(&iotrack);
continue;
}
}
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;
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);
}
}
AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
}
Double_t AliTPCtracker::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
Int_t ctype = cl->GetType();
Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
Double_t angle = seed->GetSnp()*seed->GetSnp();
angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
if (ctype<0) {
erry2+=0.5;
}
erry2*=erry2;
Double_t addErr=0;
const Double_t *errInner = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorClusterInner();
addErr=errInner[0]*TMath::Exp(-TMath::Abs((cl->GetX()-85.)/errInner[1]));
erry2+=addErr*addErr;
const Double_t *errCluster = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorCluster();
erry2+=errCluster[0]*errCluster[0];
seed->SetErrorY2(erry2);
return erry2;
}
Double_t AliTPCtracker::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
Int_t ctype = cl->GetType();
Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
Double_t angle2 = seed->GetSnp()*seed->GetSnp();
angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
if (ctype<0) {
errz2+=0.5;
}
errz2*=errz2;
Double_t addErr=0;
const Double_t *errInner = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorClusterInner();
addErr=errInner[0]*TMath::Exp(-TMath::Abs((cl->GetX()-85.)/errInner[1]));
errz2+=addErr*addErr;
const Double_t *errCluster = AliTPCReconstructor::GetRecoParam()->GetSystematicErrorCluster();
errz2+=errCluster[1]*errCluster[1];
seed->SetErrorZ2(errz2);
return errz2;
}
void AliTPCtracker::RotateToLocal(AliTPCseed *seed)
{
Float_t x = seed->GetX();
Float_t y = seed->GetY();
Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
if (y > ymax) {
seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
if (!seed->Rotate(fSectors->GetAlpha()))
return;
} else if (y <-ymax) {
seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
if (!seed->Rotate(-fSectors->GetAlpha()))
return;
}
}
Double_t AliTPCtracker::F1old(Double_t x1,Double_t y1,
Double_t x2,Double_t y2,
Double_t x3,Double_t y3) const
{
Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
(y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
(x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
if ( xr*xr+yr*yr<=0.00000000000001) return 100;
return -xr*yr/sqrt(xr*xr+yr*yr);
}
Double_t AliTPCtracker::F1(Double_t x1,Double_t y1,
Double_t x2,Double_t y2,
Double_t x3,Double_t y3) const
{
x3 -=x1;
x2 -=x1;
y3 -=y1;
y2 -=y1;
Double_t det = x3*y2-x2*y3;
if (TMath::Abs(det)<1e-10){
return 100;
}
Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
Double_t x0 = x3*0.5-y3*u;
Double_t y0 = y3*0.5+x3*u;
Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
if (det<0) c2*=-1;
return c2;
}
Double_t AliTPCtracker::F2(Double_t x1,Double_t y1,
Double_t x2,Double_t y2,
Double_t x3,Double_t y3) const
{
x3 -=x1;
x2 -=x1;
y3 -=y1;
y2 -=y1;
Double_t det = x3*y2-x2*y3;
if (TMath::Abs(det)<1e-10) {
return 100;
}
Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
Double_t x0 = x3*0.5-y3*u;
Double_t y0 = y3*0.5+x3*u;
Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
if (det<0) c2*=-1;
x0+=x1;
x0*=c2;
return x0;
}
Double_t AliTPCtracker::F2old(Double_t x1,Double_t y1,
Double_t x2,Double_t y2,
Double_t x3,Double_t y3) const
{
Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
(y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
(x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
}
Double_t AliTPCtracker::F3(Double_t x1,Double_t y1,
Double_t x2,Double_t y2,
Double_t z1,Double_t z2) const
{
return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}
Double_t AliTPCtracker::F3n(Double_t x1,Double_t y1,
Double_t x2,Double_t y2,
Double_t z1,Double_t z2, Double_t c) const
{
Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
if (TMath::Abs(d*c*0.5)>1) return 0;
Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
angle2 = (z1-z2)*c/(angle2*2.);
return angle2;
}
Bool_t AliTPCtracker::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
{
Double_t dx=x2-x1;
if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
return kFALSE;
}
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];
Double_t dy = dx*(c1+c2)/(r1+r2);
Double_t dz = 0;
Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
if (TMath::Abs(delta)>0.01){
dz = x[3]*TMath::ASin(delta)/x[4];
}else{
dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
}
y+=dy;
z+=dz;
return kTRUE;
}
Int_t AliTPCtracker::LoadClusters (TTree *const tree)
{
fInput = tree;
return LoadClusters();
}
Int_t AliTPCtracker::LoadClusters(const TObjArray *arr)
{
AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
Int_t lower = arr->LowerBound();
Int_t entries = arr->GetEntriesFast();
for (Int_t i=lower; i<entries; i++) {
clrow = (AliTPCClustersRow*) arr->At(i);
if(!clrow) continue;
if(!clrow->GetArray()) continue;
Int_t sec,row;
fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
}
if (clrow->GetArray()->GetEntriesFast()<=0) continue;
AliTPCtrackerRow * tpcrow=0;
Int_t left=0;
if (sec<fkNIS*2){
tpcrow = &(fInnerSec[sec%fkNIS][row]);
left = sec/fkNIS;
}
else{
tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
left = (sec-fkNIS*2)/fkNOS;
}
if (left ==0){
tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
for (Int_t j=0;j<tpcrow->GetN1();++j)
tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
}
if (left ==1){
tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
for (Int_t j=0;j<tpcrow->GetN2();++j)
tpcrow->SetCluster2(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
}
clrow->GetArray()->Clear("C");
}
delete clrow;
LoadOuterSectors();
LoadInnerSectors();
return 0;
}
Int_t AliTPCtracker::LoadClusters(const TClonesArray *arr)
{
AliTPCclusterMI *clust=0;
Int_t count[72][96] = { {0} , {0} };
for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
clust = (AliTPCclusterMI*)arr->At(icl);
if(!clust) continue;
Transform(clust);
count[clust->GetDetector()][clust->GetRow()]++;
}
for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
clust = (AliTPCclusterMI*)arr->At(icl);
if(!clust) continue;
Int_t sec = clust->GetDetector();
Int_t row = clust->GetRow();
if(sec<fkNIS*2) {
if(row == 30) continue;
}
else {
if(row == 27 || row == 76) continue;
}
if (sec<fkNIS*2){
fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
}
else{
fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
}
}
return 0;
}
Int_t AliTPCtracker::LoadClusters()
{
static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
AliInfo("LoadClusters()\n");
TTree * tree = fInput;
TBranch * br = tree->GetBranch("Segment");
br->SetAddress(&clrow);
Int_t j=Int_t(tree->GetEntries());
for (Int_t i=0; i<j; i++) {
br->GetEntry(i);
Int_t sec,row;
fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
}
AliTPCtrackerRow * tpcrow=0;
Int_t left=0;
if (sec<fkNIS*2){
tpcrow = &(fInnerSec[sec%fkNIS][row]);
left = sec/fkNIS;
}
else{
tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
left = (sec-fkNIS*2)/fkNOS;
}
if (left ==0){
tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
for (Int_t k=0;k<tpcrow->GetN1();++k)
tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
}
if (left ==1){
tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
for (Int_t k=0;k<tpcrow->GetN2();++k)
tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
}
}
clrow->Clear("C");
LoadOuterSectors();
LoadInnerSectors();
cout << " =================================================================================================================================== " << endl;
cout << " AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection() << endl;
cout << " AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() = " << AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection() << endl;
cout << " =================================================================================================================================== " << endl;
if (AliTPCReconstructor::GetRecoParam()->GetUseIonTailCorrection()) ApplyTailCancellation();
if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) CalculateXtalkCorrection();
if (AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrection()!=0.) ApplyXtalkCorrection();
return 0;
}
void AliTPCtracker::CalculateXtalkCorrection(){
const Int_t nROCs = 72;
const Int_t nIterations=3;
for (Int_t isector=0; isector<nROCs*4; isector++){
TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
if (crossTalkMatrix)(*crossTalkMatrix)*=0;
}
Double_t missingChargeFactor= AliTPCReconstructor::GetRecoParam()->GetCrosstalkCorrectionMissingCharge();
for (Int_t iter=0; iter<nIterations;iter++){
for (Int_t isector=0; isector<36; isector++){
for (Int_t iside=0; iside<2; iside++){
AliTPCtrackerSector §or= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
Int_t nrows = sector.GetNRows();
Int_t sec=0;
if (isector<18) sec=isector+18*iside;
if (isector>=18) sec=18+isector+18*iside;
for (Int_t row = 0;row<nrows;row++){
Int_t wireSegmentID = fkParam->GetWireSegment(sec,row);
Float_t nPadsPerSegment = (Float_t)(fkParam->GetNPadsPerSegment(wireSegmentID));
TMatrixD &crossTalkSignal = *((TMatrixD*)fCrossTalkSignalArray->At(sec));
TMatrixD &crossTalkSignalCache = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs*2));
TMatrixD &crossTalkSignalBelow = *((TMatrixD*)fCrossTalkSignalArray->At(sec+nROCs));
Int_t nCols=crossTalkSignal.GetNcols();
AliTPCtrackerRow& tpcrow = sector[row];
Int_t ncl = tpcrow.GetN1();
if (iside>0) ncl=tpcrow.GetN2();
for (Int_t i=0;i<ncl;i++) {
AliTPCclusterMI *clXtalk= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
Int_t timeBinXtalk = clXtalk->GetTimeBin();
Double_t rmsPadMin2=0.5*0.5+(fkParam->GetDiffT()*fkParam->GetDiffT())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec));
Double_t rmsTimeMin2=1+(fkParam->GetDiffL()*fkParam->GetDiffL())*(TMath::Abs((clXtalk->GetZ()-fkParam->GetZLength())))/(fkParam->GetZWidth()*fkParam->GetZWidth());
Double_t rmsTime2 = clXtalk->GetSigmaZ2()/(fkParam->GetZWidth()*fkParam->GetZWidth());
Double_t rmsPad2 = clXtalk->GetSigmaY2()/(fkParam->GetPadPitchWidth(sec)*fkParam->GetPadPitchWidth(sec));
if (rmsPadMin2>rmsPad2){
rmsPad2=rmsPadMin2;
}
if (rmsTimeMin2>rmsTime2){
rmsTime2=rmsTimeMin2;
}
Double_t norm= 2.*TMath::Exp(-1.0/(2.*rmsTime2))+2.*TMath::Exp(-4.0/(2.*rmsTime2))+1.;
Double_t qTotXtalk = 0.;
Double_t qTotXtalkMissing = 0.;
for (Int_t itb=timeBinXtalk-2, idelta=-2; itb<=timeBinXtalk+2; itb++,idelta++) {
if (itb<0 || itb>=nCols) continue;
Double_t missingCharge=0;
Double_t trf= TMath::Exp(-idelta*idelta/(2.*rmsTime2));
if (missingChargeFactor>0) {
for (Int_t dpad=-2; dpad<=2; dpad++){
Double_t qPad = clXtalk->GetMax()*TMath::Exp(-dpad*dpad/(2.*rmsPad2))*trf;
if (TMath::Nint(qPad-crossTalkSignalCache[wireSegmentID][itb])<=fkParam->GetZeroSup()){
missingCharge+=qPad+crossTalkSignalCache[wireSegmentID][itb];
}else{
missingCharge+=crossTalkSignalCache[wireSegmentID][itb];
}
}
}
qTotXtalk = clXtalk->GetQ()*trf/norm+missingCharge*missingChargeFactor;
qTotXtalkMissing = missingCharge;
crossTalkSignal[wireSegmentID][itb]+= qTotXtalk/nPadsPerSegment;
crossTalkSignalBelow[wireSegmentID][itb]+= qTotXtalkMissing/nPadsPerSegment;
}
}
}
}
}
if (AliTPCReconstructor::StreamLevel()&kStreamCrosstalkMatrix) {
for (Int_t isector=0; isector<nROCs; isector++){
TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
TMatrixD * crossTalkMatrixCache = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs*2);
TVectorD vecAll(crossTalkMatrix->GetNrows());
TVectorD vecBelow(crossTalkMatrix->GetNrows());
TVectorD vecCache(crossTalkMatrixCache->GetNrows());
for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
vecCache[iwire]=(*crossTalkMatrixCache)(iwire,itime);
}
(*fDebugStreamer)<<"crosstalkMatrix"<<
"iter="<<iter<<
"sector="<<isector<<
"itime="<<itime<<
"vecAll.="<<&vecAll<<
"vecCache.="<<&vecCache<<
"vecBelow.="<<&vecBelow<<
"\n";
}
}
}
if (iter<nIterations-1) for (Int_t isector=0; isector<nROCs*2; isector++){
TMatrixD * crossTalkMatrix = (TMatrixD*)fCrossTalkSignalArray->At(isector);
TMatrixD * crossTalkMatrixCache = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs*2);
if (crossTalkMatrix){
(*crossTalkMatrixCache)*=0;
(*crossTalkMatrixCache)+=(*crossTalkMatrix);
(*crossTalkMatrix)*=0;
}
}
}
}
void AliTPCtracker::FilterOutlierClusters(){
Int_t nSectors=AliTPCROC::Instance()->GetNSectors();
Int_t nTimeBins= 1100;
Int_t nPadRows=(AliTPCROC::Instance()->GetNRows(0) + AliTPCROC::Instance()->GetNRows(36));
const Double_t nSigmaCut=9.;
const Double_t offsetTime=100;
const Double_t offsetPadRow=300;
const Double_t offsetTimeAccept=8;
TH2F hisTime("hisSectorTime","hisSectorTime", nSectors,0,nSectors, nTimeBins,0,nTimeBins);
TH2F hisPadRow("hisSectorRow","hisSectorRow", nSectors,0,nSectors, nPadRows,0,nPadRows);
for (Int_t isector=0; isector<36; isector++){
for (Int_t iside=0; iside<2; iside++){
AliTPCtrackerSector §or= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
Int_t nrows = sector.GetNRows();
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow& tpcrow = sector[row];
Int_t ncl = tpcrow.GetN1();
if (iside>0) ncl=tpcrow.GetN2();
for (Int_t i=0;i<ncl;i++) {
AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
hisTime.Fill(cluster->GetDetector(),cluster->GetTimeBin());
hisPadRow.Fill(cluster->GetDetector(),cluster->GetRow());
}
}
}
}
TVectorD vecTime(nTimeBins);
TVectorD vecPadRow(nPadRows);
TVectorD vecMedianSectorTime(nSectors);
TVectorD vecRMSSectorTime(nSectors);
TVectorD vecMedianSectorTimeOut6(nSectors);
TVectorD vecMedianSectorTimeOut9(nSectors);
TVectorD vecMedianSectorTimeOut(nSectors);
TVectorD vecMedianSectorPadRow(nSectors);
TVectorD vecRMSSectorPadRow(nSectors);
TVectorD vecMedianSectorPadRowOut6(nSectors);
TVectorD vecMedianSectorPadRowOut9(nSectors);
TVectorD vecMedianSectorPadRowOut(nSectors);
TVectorD vecSectorOut6(nSectors);
TVectorD vecSectorOut9(nSectors);
TMatrixD matSectorCluster(nSectors,2);
for (Int_t isec=0; isec<nSectors; isec++){
vecMedianSectorTimeOut6[isec]=0;
vecMedianSectorTimeOut9[isec]=0;
for (Int_t itime=0; itime<nTimeBins; itime++){
vecTime[itime]=hisTime.GetBinContent(isec+1, itime+1);
}
Double_t median= TMath::Mean(nTimeBins,vecTime.GetMatrixArray());
Double_t rms= TMath::RMS(nTimeBins,vecTime.GetMatrixArray());
vecMedianSectorTime[isec]=median;
vecRMSSectorTime[isec]=rms;
if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector TimeStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
for (Int_t itime=0; itime<nTimeBins; itime++){
Double_t entries= hisTime.GetBinContent(isec+1, itime+1);
if (entries>median+6.*rms+offsetTime) {
vecMedianSectorTimeOut6[isec]+=1;
}
if (entries>median+9.*rms+offsetTime) {
vecMedianSectorTimeOut9[isec]+=1;
}
}
}
for (Int_t isec=0; isec<nSectors; isec++){
vecMedianSectorPadRowOut6[isec]=0;
vecMedianSectorPadRowOut9[isec]=0;
for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
vecPadRow[ipadrow]=hisPadRow.GetBinContent(isec+1, ipadrow+1);
}
Int_t nPadRowsSector= AliTPCROC::Instance()->GetNRows(isec);
Double_t median= TMath::Mean(nPadRowsSector,vecPadRow.GetMatrixArray());
Double_t rms= TMath::RMS(nPadRowsSector,vecPadRow.GetMatrixArray());
vecMedianSectorPadRow[isec]=median;
vecRMSSectorPadRow[isec]=rms;
if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector PadRowStat: %d\t%8.0f\t%8.0f",isec,median,rms).Data());
for (Int_t ipadrow=0; ipadrow<nPadRows; ipadrow++){
Double_t entries= hisPadRow.GetBinContent(isec+1, ipadrow+1);
if (entries>median+6.*rms+offsetPadRow) {
vecMedianSectorPadRowOut6[isec]+=1;
}
if (entries>median+9.*rms+offsetPadRow) {
vecMedianSectorPadRowOut9[isec]+=1;
}
}
}
Double_t medianSectorTime = TMath::Median(nSectors, vecTime.GetMatrixArray());
Double_t mean69SectorTime, rms69SectorTime=0;
AliMathBase::EvaluateUni(nSectors, vecTime.GetMatrixArray(), mean69SectorTime,rms69SectorTime,69);
for (Int_t isec=0; isec<nSectors; isec++){
vecSectorOut6[isec]=0;
vecSectorOut9[isec]=0;
matSectorCluster(isec,0)=0;
matSectorCluster(isec,1)=0;
if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+6.*(rms69SectorTime+ offsetTimeAccept))) {
vecSectorOut6[isec]=1;
}
if (TMath::Abs(vecMedianSectorTime[isec])>(mean69SectorTime+9.*(rms69SectorTime+ offsetTimeAccept))){
vecSectorOut9[isec]=1;
}
}
Int_t filteredSector= vecSectorOut9.Sum();
Int_t filteredSectorTime= vecMedianSectorTimeOut9.Sum();
Int_t filteredSectorPadRow= vecMedianSectorPadRowOut9.Sum();
if (fEvent) if (fEvent->GetHeader()){
fEvent->GetHeader()->SetTPCNoiseFilterCounter(0,TMath::Min(filteredSector,255));
fEvent->GetHeader()->SetTPCNoiseFilterCounter(1,TMath::Min(filteredSectorTime,255));
fEvent->GetHeader()->SetTPCNoiseFilterCounter(2,TMath::Min(filteredSectorPadRow,255));
}
Int_t counterAll=0;
Int_t counterOut=0;
for (Int_t isector=0; isector<36; isector++){
for (Int_t iside=0; iside<2; iside++){
AliTPCtrackerSector §or= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
Int_t nrows = sector.GetNRows();
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow& tpcrow = sector[row];
Int_t ncl = tpcrow.GetN1();
if (iside>0) ncl=tpcrow.GetN2();
for (Int_t i=0;i<ncl;i++) {
AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
Double_t medianTime=vecMedianSectorTime[cluster->GetDetector()];
Double_t medianPadRow=vecMedianSectorPadRow[cluster->GetDetector()];
Double_t rmsTime=vecRMSSectorTime[cluster->GetDetector()];
Double_t rmsPadRow=vecRMSSectorPadRow[cluster->GetDetector()];
Int_t entriesPadRow=hisPadRow.GetBinContent(cluster->GetDetector()+1, cluster->GetRow()+1);
Int_t entriesTime=hisTime.GetBinContent(cluster->GetDetector()+1, cluster->GetTimeBin()+1);
Bool_t isOut=kFALSE;
if (vecSectorOut9[cluster->GetDetector()]>0.5) {
isOut=kTRUE;
}
if (entriesTime>medianTime+nSigmaCut*rmsTime+offsetTime) {
isOut=kTRUE;
vecMedianSectorTimeOut[cluster->GetDetector()]++;
}
if (entriesPadRow>medianPadRow+nSigmaCut*rmsPadRow+offsetPadRow) {
isOut=kTRUE;
vecMedianSectorPadRowOut[cluster->GetDetector()]++;
}
counterAll++;
matSectorCluster(cluster->GetDetector(),0)+=1;
if (isOut){
cluster->Disable();
counterOut++;
matSectorCluster(cluster->GetDetector(),1)+=1;
}
}
}
}
}
for (Int_t isec=0; isec<nSectors; isec++){
if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) AliInfo(TString::Format("Sector Stat: %d\t%8.0f\t%8.0f",isec,matSectorCluster(isec,1),matSectorCluster(isec,0)).Data());
}
if ((AliTPCReconstructor::StreamLevel()&kStreamFilterClusterInfo)>0) {
AliLog::Flush();
AliInfo(TString::Format("Cluster counter: (%d/%d) (Filtered/All)",counterOut,counterAll).Data());
for (Int_t iSec=0; iSec<nSectors; iSec++){
if (vecSectorOut9[iSec]>0 || matSectorCluster(iSec,1)>0) {
AliInfo(TString::Format("Filtered sector\t%d",iSec).Data());
Double_t vecMedTime =TMath::Median(72,vecMedianSectorTime.GetMatrixArray());
Double_t vecMedPadRow =TMath::Median(72,vecMedianSectorPadRow.GetMatrixArray());
Double_t vecMedCluster=(counterAll-counterOut)/72;
AliInfo(TString::Format("VecMedianSectorTime\t(%4.4f/%4.4f/%4.4f)", vecMedianSectorTimeOut[iSec],vecMedianSectorTime[iSec],vecMedTime).Data());
AliInfo(TString::Format("VecMedianSectorPadRow\t(%4.4f/%4.4f/%4.4f)", vecMedianSectorPadRowOut[iSec],vecMedianSectorPadRow[iSec],vecMedPadRow).Data());
AliInfo(TString::Format("MatSectorCluster\t(%4.4f/%4.4f/%4.4f)\n", matSectorCluster(iSec,1), matSectorCluster(iSec,0), vecMedCluster).Data());
AliLog::Flush();
}
}
AliLog::Flush();
Int_t eventNr = fEvent->GetEventNumberInFile();
(*fDebugStreamer)<<"filterClusterInfo"<<
"eventNr="<<eventNr<<
"counterAll="<<counterAll<<
"counterOut="<<counterOut<<
"matSectotCluster.="<<&matSectorCluster<<
"filteredSector="<<filteredSector<<
"filteredSectorTime="<<filteredSectorTime<<
"filteredSectorPadRow="<<filteredSectorPadRow<<
"medianSectorTime="<<medianSectorTime<<
"mean69SectorTime="<<mean69SectorTime<<
"rms69SectorTime="<<rms69SectorTime<<
"vecSectorOut6.="<<&vecSectorOut6<<
"vecSectorOut9.="<<&vecSectorOut9<<
"vecMedianSectorTime.="<<&vecMedianSectorTime<<
"vecRMSSectorTime.="<<&vecRMSSectorTime<<
"vecMedianSectorTimeOut6.="<<&vecMedianSectorTimeOut6<<
"vecMedianSectorTimeOut9.="<<&vecMedianSectorTimeOut9<<
"vecMedianSectorTimeOut0.="<<&vecMedianSectorTimeOut<<
"vecMedianSectorPadRow.="<<&vecMedianSectorPadRow<<
"vecRMSSectorPadRow.="<<&vecRMSSectorPadRow<<
"vecMedianSectorPadRowOut6.="<<&vecMedianSectorPadRowOut6<<
"vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut9<<
"vecMedianSectorPadRowOut9.="<<&vecMedianSectorPadRowOut<<
"\n";
((*fDebugStreamer)<<"filterClusterInfo").GetTree()->Write();
fDebugStreamer->GetFile()->Flush();
}
}
void AliTPCtracker::UnloadClusters()
{
Int_t nrows = fOuterSec->GetNRows();
for (Int_t sec = 0;sec<fkNOS;sec++)
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
tpcrow->ResetClusters();
}
nrows = fInnerSec->GetNRows();
for (Int_t sec = 0;sec<fkNIS;sec++)
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
tpcrow->ResetClusters();
}
return ;
}
void AliTPCtracker::FillClusterArray(TObjArray* array) const{
Int_t nrows=0;
nrows = fOuterSec->GetNRows();
for (Int_t sec = 0;sec<fkNOS;sec++)
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
if (!tpcrow) continue;
for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
array->AddLast((TObject*)((*tpcrow)[icl]));
}
}
nrows = fInnerSec->GetNRows();
for (Int_t sec = 0;sec<fkNIS;sec++)
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
if (!tpcrow) continue;
for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
array->AddLast((TObject*)(*tpcrow)[icl]);
}
}
}
void AliTPCtracker::Transform(AliTPCclusterMI * cluster){
AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
AliTPCTransform *transform = calibDB->GetTransform() ;
if (!transform) {
AliFatal("Tranformations not in calibDB");
return;
}
transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),static_cast<Double_t>(cluster->GetPad()),static_cast<Double_t>(cluster->GetTimeBin())};
Int_t i[1]={cluster->GetDetector()};
transform->Transform(x,i,0,1);
if ((AliTPCReconstructor::StreamLevel()&kStreamTransform)>0) {
Float_t gx[3];
cluster->GetGlobalXYZ(gx);
Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Transform"<<
"event="<<event<<
"x0="<<x[0]<<
"x1="<<x[1]<<
"x2="<<x[2]<<
"gx0="<<gx[0]<<
"gx1="<<gx[1]<<
"gx2="<<gx[2]<<
"Cl.="<<cluster<<
"\n";
}
cluster->SetX(x[0]);
cluster->SetY(x[1]);
cluster->SetZ(x[2]);
if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
if (mat) mat->LocalToMaster(pos,posC);
else{
}
cluster->SetX(posC[0]);
cluster->SetY(posC[1]);
cluster->SetZ(posC[2]);
}
}
void AliTPCtracker::ApplyXtalkCorrection(){
for (Int_t isector=0; isector<36; isector++){
for (Int_t iside=0; iside<2; iside++){
AliTPCtrackerSector §or= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
Int_t nrows = sector.GetNRows();
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow& tpcrow = sector[row];
Int_t ncl = tpcrow.GetN1();
if (iside>0) ncl=tpcrow.GetN2();
Int_t xSector=0;
if (isector<18){
xSector=isector+(iside>0)*18;
}else{
xSector=isector+18;
if (iside>0) xSector+=18;
}
TMatrixD &crossTalkMatrix= *((TMatrixD*)fCrossTalkSignalArray->At(xSector));
Int_t wireSegmentID = fkParam->GetWireSegment(xSector,row);
for (Int_t i=0;i<ncl;i++) {
AliTPCclusterMI *cluster= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
Int_t iTimeBin=TMath::Nint(cluster->GetTimeBin());
Double_t xTalk= crossTalkMatrix[wireSegmentID][iTimeBin];
cluster->SetMax(cluster->GetMax()+xTalk);
const Double_t kDummy=4;
Double_t sumxTalk=xTalk*kDummy;
cluster->SetQ(cluster->GetQ()+sumxTalk);
if ((AliTPCReconstructor::StreamLevel()&kStreamXtalk)>0) {
TTreeSRedirector &cstream = *fDebugStreamer;
if (gRandom->Rndm() > 0.){
cstream<<"Xtalk"<<
"isector=" << isector <<
"iside=" << iside <<
"row=" << row <<
"i=" << i <<
"xSector=" << xSector <<
"wireSegmentID=" << wireSegmentID <<
"iTimeBin=" << iTimeBin <<
"xTalk=" << xTalk <<
"sumxTalk=" << sumxTalk <<
"cluster.=" << cluster <<
"\n";
}
}
}
}
}
}
}
void AliTPCtracker::ApplyTailCancellation(){
TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
TObject *rocFactorIROC = ionTailArr->FindObject("factorIROC");
TObject *rocFactorOROC = ionTailArr->FindObject("factorOROC");
Float_t factorIROC = (atof(rocFactorIROC->GetTitle()));
Float_t factorOROC = (atof(rocFactorOROC->GetTitle()));
Int_t nclALL=0;
for (Int_t isector=0; isector<36; isector++){
AliTPCtrackerSector §or= (isector<18)?fInnerSec[isector%18]:fOuterSec[isector%18];
nclALL += sector.GetNClInSector(0);
nclALL += sector.GetNClInSector(1);
}
for (Int_t iside=0; iside<2; iside++){
for (Int_t secType=0; secType<2; secType++){
const Float_t ampfactor = (secType==0)?factorIROC:factorOROC;
std::cout << " ampfactor = " << ampfactor << std::endl;
for (Int_t sec = 0;sec<fkNOS;sec++){
TGraphErrors ** graphRes = new TGraphErrors *[20];
Float_t * indexAmpGraphs = new Float_t[20];
for (Int_t icache=0; icache<20; icache++)
{
graphRes[icache] = NULL;
indexAmpGraphs[icache] = 0;
}
if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(sec+36*secType+18*iside,graphRes,indexAmpGraphs))
{
continue;
}
AliTPCtrackerSector §or= (secType==0)?fInnerSec[sec]:fOuterSec[sec];
Int_t nrows = sector.GetNRows();
Int_t nclSector = sector.GetNClInSector(iside);
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow& tpcrow = sector[row];
Int_t ncl = tpcrow.GetN1();
if (iside>0) ncl=tpcrow.GetN2();
Float_t qTotArray[ncl];
Float_t qMaxArray[ncl];
Int_t sortedClusterIndex[ncl];
Float_t sortedClusterTimeBin[ncl];
TObjArray *rowClusterArray = new TObjArray(ncl);
for (Int_t i=0;i<ncl;i++)
{
qTotArray[i]=0;
qMaxArray[i]=0;
sortedClusterIndex[i]=i;
AliTPCclusterMI *rowcl= (iside>0)?(tpcrow.GetCluster2(i)):(tpcrow.GetCluster1(i));
if (rowcl) {
rowClusterArray->AddAt(rowcl,i);
} else {
rowClusterArray->RemoveAt(i);
}
if (!rowcl) {
sortedClusterTimeBin[i]=0.0;
} else {
sortedClusterTimeBin[i] = rowcl->GetTimeBin();
}
}
TMath::Sort(ncl,sortedClusterTimeBin,sortedClusterIndex,kFALSE);
for (Int_t icl0=0; icl0<ncl;icl0++){
AliTPCclusterMI *cl0= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl0]));
if (!cl0) continue;
Int_t nclPad=0;
for (Int_t icl1=0; icl1<ncl;icl1++){
AliTPCclusterMI *cl1= static_cast<AliTPCclusterMI*>(rowClusterArray->At(sortedClusterIndex[icl1]));
if (!cl1) continue;
if (TMath::Abs(cl0->GetPad()-cl1->GetPad())>4) continue;
if (cl0->GetTimeBin()<= cl1->GetTimeBin()) continue;
if (TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin())>600) continue;
if (TMath::Abs(cl0->GetPad()-cl1->GetPad())<4) nclPad++;
Double_t ionTailMax=0.;
Double_t ionTailTotal=0.;
GetTailValue(ampfactor,ionTailMax,ionTailTotal,graphRes,indexAmpGraphs,cl0,cl1);
ionTailMax=TMath::Abs(ionTailMax);
ionTailTotal=TMath::Abs(ionTailTotal);
qTotArray[icl0]+=ionTailTotal;
qMaxArray[icl0]+=ionTailMax;
if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {
TTreeSRedirector &cstream = *fDebugStreamer;
if (gRandom->Rndm() > 0.999){
cstream<<"IonTail"<<
"cl0.=" <<cl0 <<
"cl1.=" <<cl1 <<
"ionTailTotal=" <<ionTailTotal <<
"ionTailMax=" <<ionTailMax <<
"\n";
}
}
}
cl0->SetQ(TMath::Nint(Float_t(cl0->GetQ())+Float_t(qTotArray[icl0])));
cl0->SetMax(TMath::Nint(Float_t(cl0->GetMax())+qMaxArray[icl0]));
if ((AliTPCReconstructor::StreamLevel()&kStreamIonTail)>0) {
TTreeSRedirector &cstream = *fDebugStreamer;
if (gRandom->Rndm() > 0.999){
cstream<<"IonTailCorrected"<<
"cl0.=" << cl0 <<
"ionTailTotalPerCluster=" << qTotArray[icl0] <<
"ionTailMaxPerCluster=" << qMaxArray[icl0] <<
"nclALL=" << nclALL <<
"nclSector=" << nclSector <<
"nclRow=" << ncl <<
"nclPad=" << nclPad <<
"row=" << row <<
"sector=" << sec <<
"icl0=" << icl0 <<
"\n";
}
}
}
delete rowClusterArray;
}
for (int i=0; i<20; i++) delete graphRes[i];
delete [] graphRes;
delete [] indexAmpGraphs;
}
}
}
}
void AliTPCtracker::GetTailValue(Float_t ampfactor,Double_t &ionTailMax, Double_t &ionTailTotal,TGraphErrors **graphRes,Float_t *indexAmpGraphs,AliTPCclusterMI *cl0,AliTPCclusterMI *cl1){
const Double_t kMinPRF = 0.5;
ionTailTotal = 0.;
ionTailMax = 0.;
Float_t qTot0 = cl0->GetQ();
Float_t qTot1 = cl1->GetQ();
Int_t sectorPad = cl1->GetDetector();
Int_t padcl0 = TMath::Nint(cl0->GetPad());
Int_t padcl1 = TMath::Nint(cl1->GetPad());
Float_t padWidth = (sectorPad < 36)?0.4:0.6;
const Int_t deltaTimebin = TMath::Nint(TMath::Abs(cl1->GetTimeBin()-cl0->GetTimeBin()))+12;
Double_t rmsPad1 = (cl1->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl1->GetSigmaY2())/padWidth);
Double_t rmsPad0 = (cl0->GetSigmaY2()==0)?kMinPRF:(TMath::Sqrt(cl0->GetSigmaY2())/padWidth);
Double_t sumAmp1=0.;
for (Int_t idelta =-2; idelta<=2;idelta++){
sumAmp1+=TMath::Exp(-idelta*idelta/(2*rmsPad1));
}
Double_t sumAmp0=0.;
for (Int_t idelta =-2; idelta<=2;idelta++){
sumAmp0+=TMath::Exp(-idelta*idelta/(2*rmsPad0));
}
Int_t padScan=2;
for (Int_t ipad1=padcl1-padScan; ipad1<=padcl1+padScan; ipad1++) {
Float_t deltaPad1 = TMath::Abs(cl1->GetPad()-(Float_t)ipad1);
Double_t amp1 = (TMath::Exp(-(deltaPad1*deltaPad1)/(2*rmsPad1)))/sumAmp1;
Float_t qTotPad1 = amp1*qTot1;
Int_t ampIndex = 0;
Float_t diffAmp = TMath::Abs(deltaPad1-indexAmpGraphs[0]);
for (Int_t j=0;j<20;j++) {
if (diffAmp > TMath::Abs(deltaPad1-indexAmpGraphs[j]) && indexAmpGraphs[j]!=0)
{
diffAmp = TMath::Abs(deltaPad1-indexAmpGraphs[j]);
ampIndex = j;
}
}
if (!graphRes[ampIndex]) continue;
if (deltaTimebin+2 >= graphRes[ampIndex]->GetN()) continue;
if (graphRes[ampIndex]->GetY()[deltaTimebin+2]>=0) continue;
for (Int_t ipad0=padcl0-padScan; ipad0<=padcl0+padScan; ipad0++) {
if (ipad1!=ipad0) continue;
Float_t deltaPad0 = TMath::Abs(cl0->GetPad()-(Float_t)ipad0);
Double_t amp0 = (TMath::Exp(-(deltaPad0*deltaPad0)/(2*rmsPad0)))/sumAmp0;
Float_t qMaxPad0 = amp0*qTot0;
for (Int_t itb=deltaTimebin-2; itb<=deltaTimebin+2; itb++) {
if (itb<0) continue;
if (itb>=graphRes[ampIndex]->GetN()) continue;
Float_t tailCorr = TMath::Abs((qTotPad1*ampfactor)*(graphRes[ampIndex])->GetY()[itb]);
if (ipad1!=padcl0) {
ionTailTotal += TMath::Min(qMaxPad0,tailCorr);
} else {
ionTailTotal += tailCorr;
}
if (itb == deltaTimebin && ipad1 == padcl0) ionTailMax += tailCorr;
}
}
}
}
Int_t AliTPCtracker::LoadOuterSectors() {
Int_t nrows = fOuterSec->GetNRows();
UInt_t index=0;
for (Int_t sec = 0;sec<fkNOS;sec++)
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
Int_t sec2 = sec+2*fkNIS;
Int_t ncl = tpcrow->GetN1();
while (ncl--) {
AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
index=(((sec2<<8)+row)<<16)+ncl;
tpcrow->InsertCluster(c,index);
}
ncl = tpcrow->GetN2();
while (ncl--) {
AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
tpcrow->InsertCluster(c,index);
}
for (Int_t i=0;i<510;i++)
tpcrow->SetFastCluster(i,-1);
for (Int_t i=0;i<tpcrow->GetN();i++){
Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
tpcrow->SetFastCluster(zi,i);
}
Int_t last = 0;
for (Int_t i=0;i<510;i++){
if (tpcrow->GetFastCluster(i)<0)
tpcrow->SetFastCluster(i,last);
else
last = tpcrow->GetFastCluster(i);
}
}
fN=fkNOS;
fSectors=fOuterSec;
return 0;
}
Int_t AliTPCtracker::LoadInnerSectors() {
Int_t nrows = fInnerSec->GetNRows();
UInt_t index=0;
for (Int_t sec = 0;sec<fkNIS;sec++)
for (Int_t row = 0;row<nrows;row++){
AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
Int_t ncl = tpcrow->GetN1();
while (ncl--) {
AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
index=(((sec<<8)+row)<<16)+ncl;
tpcrow->InsertCluster(c,index);
}
ncl = tpcrow->GetN2();
while (ncl--) {
AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
tpcrow->InsertCluster(c,index);
}
for (Int_t i=0;i<510;i++)
tpcrow->SetFastCluster(i,-1);
for (Int_t i=0;i<tpcrow->GetN();i++){
Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
tpcrow->SetFastCluster(zi,i);
}
Int_t last = 0;
for (Int_t i=0;i<510;i++){
if (tpcrow->GetFastCluster(i)<0)
tpcrow->SetFastCluster(i,last);
else
last = tpcrow->GetFastCluster(i);
}
}
fN=fkNIS;
fSectors=fInnerSec;
return 0;
}
AliTPCclusterMI *AliTPCtracker::GetClusterMI(Int_t index) const {
if (index<0) return 0;
Int_t sec=(index&0xff000000)>>24;
Int_t row=(index&0x00ff0000)>>16;
Int_t ncl=(index&0x00007fff)>>00;
const AliTPCtrackerRow * tpcrow=0;
TClonesArray * clrow =0;
if (sec<0 || sec>=fkNIS*4) {
AliWarning(Form("Wrong sector %d",sec));
return 0x0;
}
if (sec<fkNIS*2){
AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
if (tracksec.GetNRows()<=row) return 0;
tpcrow = &(tracksec[row]);
if (tpcrow==0) return 0;
if (sec<fkNIS) {
if (tpcrow->GetN1()<=ncl) return 0;
clrow = tpcrow->GetClusters1();
}
else {
if (tpcrow->GetN2()<=ncl) return 0;
clrow = tpcrow->GetClusters2();
}
}
else {
AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
if (tracksec.GetNRows()<=row) return 0;
tpcrow = &(tracksec[row]);
if (tpcrow==0) return 0;
if (sec-2*fkNIS<fkNOS) {
if (tpcrow->GetN1()<=ncl) return 0;
clrow = tpcrow->GetClusters1();
}
else {
if (tpcrow->GetN2()<=ncl) return 0;
clrow = tpcrow->GetClusters2();
}
}
return (AliTPCclusterMI*)clrow->At(ncl);
}
Int_t AliTPCtracker::FollowToNext(AliTPCseed& t, Int_t nr) {
Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
AliTPCclusterMI *cl=0;
Int_t tpcindex= t.GetClusterIndex2(nr);
GetShape(&t,nr);
if (fIteration>0 && tpcindex>=-1){
if (tpcindex==-1) return 0;
if (tpcindex >= 0){
cl = t.GetClusterPointer(nr);
if (!cl) cl = GetClusterMI(tpcindex);
t.SetCurrentClusterIndex1(tpcindex);
}
if (cl){
Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;
Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
if (TMath::Abs(angle-t.GetAlpha())>0.001){
Double_t rotation = angle-t.GetAlpha();
t.SetRelativeSector(relativesector);
if (!t.Rotate(rotation)) {
t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
return 0;
}
}
if (!t.PropagateTo(x)) {
t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
return 0;
}
t.SetCurrentCluster(cl);
t.SetRow(nr);
Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
if ((tpcindex&0x8000)==0) accept =0;
if (accept<3) {
if (cl->IsUsed(11)) {
t.SetErrorY2(t.GetErrorY2()+0.03);
t.SetErrorZ2(t.GetErrorZ2()+0.03);
t.SetErrorY2(t.GetErrorY2()*3);
t.SetErrorZ2(t.GetErrorZ2()*3);
}
t.SetNFoundable(t.GetNFoundable()+1);
UpdateTrack(&t,accept);
return 1;
}
else {
t.SetClusterIndex(nr, -3);
t.SetClusterPointer(nr, 0);
}
}
}
if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
if (fIteration>1 && IsFindable(t)){
t.SetNFoundable(t.GetNFoundable()+1);
return 0;
}
UInt_t index=0;
if (!t.PropagateTo(x)) {
if (fIteration==0) t.SetRemoval(10);
return 0;
}
Double_t y = t.GetY();
if (TMath::Abs(y)>ymax){
if (y > ymax) {
t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
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)) {
if (fIteration==0) t.SetRemoval(10);
return 0;
}
y = t.GetY();
}
Double_t z=t.GetZ();
if (!IsActive(t.GetRelativeSector(),nr)) {
t.SetInDead(kTRUE);
t.SetClusterIndex2(nr,-1);
return 0;
}
Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
if (!isActive || !isActive2) return 0;
const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
Double_t roady =1.;
Double_t roadz = 1.;
if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
t.SetInDead(kTRUE);
t.SetClusterIndex2(nr,-1);
return 0;
}
else
{
if (IsFindable(t))
t.SetNFoundable(t.GetNFoundable()+1);
else
return 0;
}
if (krow) {
cl = krow.FindNearest2(y,z,roady,roadz,index);
if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
}
if (cl) {
t.SetCurrentCluster(cl);
t.SetRow(nr);
if (fIteration==2&&cl->IsUsed(10)) return 0;
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);
t.SetErrorY2(t.GetErrorY2()*3);
t.SetErrorZ2(t.GetErrorZ2()*3);
}
if (accept<3) UpdateTrack(&t,accept);
} else {
if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
}
return 1;
}
Bool_t AliTPCtracker::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
{
AliTPCclusterMI *cl = GetClusterMI(index);
if (!cl) return kFALSE;
Int_t sector = (index&0xff000000)>>24;
Float_t xyz[3];
xyz[0] = cl->GetX();
xyz[1] = cl->GetY();
xyz[2] = cl->GetZ();
Float_t sin,cos;
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 < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
cov[0] = sin*sin*sigmaY2;
cov[1] = -sin*cos*sigmaY2;
cov[2] = 0.;
cov[3] = cos*cos*sigmaY2;
cov[4] = 0.;
cov[5] = sigmaZ2;
p.SetXYZ(x,y,xyz[2],cov);
AliGeomManager::ELayerID iLayer;
Int_t idet;
if (sector < fkParam->GetNInnerSector()) {
iLayer = AliGeomManager::kTPC1;
idet = sector;
}
else {
iLayer = AliGeomManager::kTPC2;
idet = sector - fkParam->GetNInnerSector();
}
UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
p.SetVolumeID(volid);
return kTRUE;
}
Int_t AliTPCtracker::UpdateClusters(AliTPCseed& t, Int_t nr) {
t.SetCurrentCluster(0);
t.SetCurrentClusterIndex1(-3);
Double_t xt=t.GetX();
Int_t row = GetRowNumber(xt)-1;
Double_t ymax= GetMaxY(nr);
if (row < nr) return 1;
Double_t x= GetXrow(nr);
Double_t y,z;
if (!t.PropagateTo(x)){
return 0;
}
y=t.GetY();
z=t.GetZ();
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;
}
return 1;
}
if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
if (!IsActive(t.GetRelativeSector(),nr)) {
t.SetInDead(kTRUE);
t.SetClusterIndex2(nr,-1);
return 0;
}
AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
t.SetInDead(kTRUE);
t.SetClusterIndex2(nr,-1);
return 0;
}
else
{
if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
else
return 0;
}
if ( (nr%2==0) || t.GetNumberOfClusters()<2){
GetShape(&t,nr);
}
AliTPCclusterMI *cl=0;
Int_t index=0;
Double_t roady = 1.;
Double_t roadz = 1.;
if (!cl){
index = t.GetClusterIndex2(nr);
if ( (index >= 0) && (index&0x8000)==0){
cl = t.GetClusterPointer(nr);
if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
t.SetCurrentClusterIndex1(index);
if (cl) {
t.SetCurrentCluster(cl);
return 1;
}
}
}
UInt_t uindex = TMath::Abs(index);
if (krow) {
cl = krow.FindNearest2(y,z,roady,roadz,uindex);
}
if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
t.SetCurrentCluster(cl);
return 1;
}
Int_t AliTPCtracker::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
if (t.GetCurrentCluster()) {
t.SetRow(nr);
Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
if (t.GetCurrentCluster()->IsUsed(10)){
t.SetNShared(t.GetNShared()+1);
if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
t.SetRemoval(10);
return 0;
}
}
if (fIteration>0) accept = 0;
if (accept<3) UpdateTrack(&t,accept);
} else {
if (fIteration==0){
if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
}
}
return 1;
}
Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
Double_t xt=t.GetX();
Double_t alpha=t.GetAlpha();
if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
if (alpha < 0. ) alpha += 2.*TMath::Pi();
t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
Int_t first = GetRowNumber(xt);
if (!fromSeeds)
first -= step;
if (first < 0)
first = 0;
for (Int_t nr= first; nr>=rf; nr-=step) {
if (t.GetKinkIndexes()[0]>0){
for (Int_t i=0;i<3;i++){
Int_t index = t.GetKinkIndexes()[i];
if (index==0) break;
if (index<0) continue;
AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
if (!kink){
printf("PROBLEM\n");
}
else{
Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
if (kinkrow==nr){
AliExternalTrackParam paramd(t);
kink->SetDaughter(paramd);
kink->SetStatus(2,5);
kink->Update();
}
}
}
}
if (nr==80) t.UpdateReference();
if (nr<fInnerSec->GetNRows())
fSectors = fInnerSec;
else
fSectors = fOuterSec;
if (FollowToNext(t,nr)==0)
if (!t.IsActive())
return 0;
}
return 1;
}
Int_t AliTPCtracker::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
Double_t xt=t.GetX();
Double_t alpha=t.GetAlpha();
if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
if (alpha < 0. ) alpha += 2.*TMath::Pi();
t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
Int_t first = t.GetFirstPoint();
Int_t ri = GetRowNumber(xt);
if (!fromSeeds)
ri += 1;
if (first<ri) first = ri;
if (first<0) first=0;
for (Int_t nr=first; nr<=rf; nr++) {
if (t.GetKinkIndexes()[0]<0){
for (Int_t i=0;i<3;i++){
Int_t index = t.GetKinkIndexes()[i];
if (index==0) break;
if (index>0) continue;
index = TMath::Abs(index);
AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
if (!kink){
printf("PROBLEM\n");
}
else{
Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
if (kinkrow==nr){
AliExternalTrackParam paramm(t);
kink->SetMother(paramm);
kink->SetStatus(2,1);
kink->Update();
}
}
}
}
if (nr<fInnerSec->GetNRows())
fSectors = fInnerSec;
else
fSectors = fOuterSec;
FollowToNext(t,nr);
}
return 1;
}
Float_t AliTPCtracker::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
{
sum1=0;
sum2=0;
Int_t sum=0;
Float_t dz2 =(s1->GetZ() - s2->GetZ());
dz2*=dz2;
Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
dy2*=dy2;
Float_t distance = TMath::Sqrt(dz2+dy2);
if (distance>4.) return 0;
Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
if (lastpoint>160)
lastpoint =160;
if (firstpoint<0)
firstpoint = 0;
if (firstpoint>lastpoint) {
firstpoint =lastpoint;
}
for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
if (s1->GetClusterIndex2(i)>0) sum1++;
if (s2->GetClusterIndex2(i)>0) sum2++;
if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
sum++;
}
}
if (sum<5) return 0;
Float_t summin = TMath::Min(sum1+1,sum2+1);
Float_t ratio = (sum+1)/Float_t(summin);
return ratio;
}
void AliTPCtracker::SignShared(AliTPCseed * s1, AliTPCseed * s2)
{
Float_t thetaCut = 0.2;
if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
Int_t sumshared=0;
Int_t firstpoint = 0;
Int_t lastpoint = 160;
for (Int_t i=firstpoint;i<lastpoint;i++){
if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
sumshared++;
}
}
if (sumshared>cutN0){
for (Int_t i=firstpoint;i<lastpoint;i++){
if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
if (s1->IsActive()&&s2->IsActive()){
p1->SetShared(kTRUE);
p2->SetShared(kTRUE);
}
}
}
}
if (sumshared>cutN0){
for (Int_t i=0;i<4;i++){
if (s1->GetOverlapLabel(3*i)==0){
s1->SetOverlapLabel(3*i, s2->GetLabel());
s1->SetOverlapLabel(3*i+1,sumshared);
s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
break;
}
}
for (Int_t i=0;i<4;i++){
if (s2->GetOverlapLabel(3*i)==0){
s2->SetOverlapLabel(3*i, s1->GetLabel());
s2->SetOverlapLabel(3*i+1,sumshared);
s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
break;
}
}
}
}
void AliTPCtracker::SignShared(TObjArray * arr)
{
for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) continue;
pt->SetSort(0);
}
arr->UnSort();
arr->Sort();
arr->Expand(arr->GetEntries());
Int_t nseed=arr->GetEntriesFast();
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) continue;
for (Int_t j=0;j<12;j++){
pt->SetOverlapLabel(j,0);
}
}
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) continue;
if (pt->GetRemoval()>10) continue;
for (Int_t j=i+1; j<nseed; j++){
AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
if (pt2->GetRemoval()<=10) {
SignShared(pt,pt2);
}
}
}
}
void AliTPCtracker::SortTracks(TObjArray * arr, Int_t mode) const
{
Int_t nseed = arr->GetEntriesFast();
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) {
continue;
}
pt->SetSort(mode);
}
arr->UnSort();
arr->Sort();
}
void AliTPCtracker::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
{
UnsignClusters();
Int_t nseed = arr->GetEntriesFast();
Float_t * quality = new Float_t[nseed];
Int_t * indexes = new Int_t[nseed];
Int_t good =0;
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt){
quality[i]=-1;
continue;
}
pt->UpdatePoints();
Float_t * points = pt->GetPoints();
if (points[3]<0.8) quality[i] =-1;
quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
}
TMath::Sort(nseed,quality,indexes);
for (Int_t itrack=0; itrack<nseed; itrack++) {
Int_t trackindex = indexes[itrack];
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
if (!pt) continue;
if (quality[trackindex]<0){
MarkSeedFree( arr->RemoveAt(trackindex) );
continue;
}
Int_t first = Int_t(pt->GetPoints()[0]);
Int_t last = Int_t(pt->GetPoints()[2]);
Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
Int_t found,foundable,shared;
pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
Bool_t itsgold =kFALSE;
if (pt->GetESD()){
Int_t dummy[12];
if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
}
if (!itsgold){
if (Float_t(shared+1)/Float_t(found+1)>factor){
if (pt->GetKinkIndexes()[0]!=0) continue;
if( (AliTPCReconstructor::StreamLevel()&kStreamRemoveUsed)>0){
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"RemoveUsed"<<
"iter="<<fIteration<<
"pt.="<<pt<<
"\n";
}
MarkSeedFree( arr->RemoveAt(trackindex) );
continue;
}
if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){
if (pt->GetKinkIndexes()[0]!=0) continue;
if( (AliTPCReconstructor::StreamLevel()&kStreamRemoveShort)>0){
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"RemoveShort"<<
"iter="<<fIteration<<
"pt.="<<pt<<
"\n";
}
MarkSeedFree( arr->RemoveAt(trackindex) );
continue;
}
}
good++;
if (pt->GetKinkIndexes()[0]>0) continue;
if (pt->GetSigmaY2()<kAlmost0) continue;
for (Int_t i=first; i<last; i++) {
Int_t index=pt->GetClusterIndex2(i);
if (index<0 || index&0x8000 ) continue;
AliTPCclusterMI *c= pt->GetClusterPointer(i);
if (!c) continue;
c->Use(10);
}
}
fNtracks = good;
if (fDebug>0){
Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
}
delete []quality;
delete []indexes;
}
void AliTPCtracker::DumpClusters(Int_t iter, TObjArray *trackArray)
{
UnsignClusters();
Int_t nseed = trackArray->GetEntries();
for (Int_t i=0; i<nseed; i++){
AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(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);
c->Use(10);
}
}
Int_t eventNr = fEvent->GetEventNumberInFile();
for (Int_t sec=0;sec<fkNIS;sec++){
for (Int_t row=0;row<fInnerSec->GetNRows();row++){
TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
Float_t gx[3]; cl->GetGlobalXYZ(gx);
(*fDebugStreamer)<<"clDump"<<
"eventNr="<<eventNr<<
"iter="<<iter<<
"cl.="<<cl<<
"gx0="<<gx[0]<<
"gx1="<<gx[1]<<
"gx2="<<gx[2]<<
"\n";
}
cla = fInnerSec[sec][row].GetClusters2();
for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
Float_t gx[3]; cl->GetGlobalXYZ(gx);
(*fDebugStreamer)<<"clDump"<<
"eventNr="<<eventNr<<
"iter="<<iter<<
"cl.="<<cl<<
"gx0="<<gx[0]<<
"gx1="<<gx[1]<<
"gx2="<<gx[2]<<
"\n";
}
}
}
for (Int_t sec=0;sec<fkNOS;sec++){
for (Int_t row=0;row<fOuterSec->GetNRows();row++){
TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
Float_t gx[3];
AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
cl->GetGlobalXYZ(gx);
(*fDebugStreamer)<<"clDump"<<
"eventNr="<<eventNr<<
"iter="<<iter<<
"cl.="<<cl<<
"gx0="<<gx[0]<<
"gx1="<<gx[1]<<
"gx2="<<gx[2]<<
"\n";
}
cla = fOuterSec[sec][row].GetClusters2();
for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
Float_t gx[3];
AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
cl->GetGlobalXYZ(gx);
(*fDebugStreamer)<<"clDump"<<
"eventNr="<<eventNr<<
"iter="<<iter<<
"cl.="<<cl<<
"gx0="<<gx[0]<<
"gx1="<<gx[1]<<
"gx2="<<gx[2]<<
"\n";
}
}
}
}
void AliTPCtracker::UnsignClusters()
{
for (Int_t sec=0;sec<fkNIS;sec++){
for (Int_t row=0;row<fInnerSec->GetNRows();row++){
TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
((AliTPCclusterMI*) cla->At(icl))->Use(-1);
cla = fInnerSec[sec][row].GetClusters2();
for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
((AliTPCclusterMI*) cla->At(icl))->Use(-1);
}
}
for (Int_t sec=0;sec<fkNOS;sec++){
for (Int_t row=0;row<fOuterSec->GetNRows();row++){
TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
((AliTPCclusterMI*) cla->At(icl))->Use(-1);
cla = fOuterSec[sec][row].GetClusters2();
for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
((AliTPCclusterMI*) cla->At(icl))->Use(-1);
}
}
}
void AliTPCtracker::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
{
Float_t sumdens=0;
Float_t sumdens2=0;
Float_t sumn =0;
Float_t sumn2 =0;
Float_t sumchi =0;
Float_t sumchi2 =0;
Float_t sum =0;
TStopwatch timer;
timer.Start();
Int_t nseed = arr->GetEntriesFast();
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) {
continue;
}
if (!(pt->IsActive())) continue;
Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
sumdens += dens;
sumdens2+= dens*dens;
sumn += pt->GetNumberOfClusters();
sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
if (chi2>5) chi2=5;
sumchi +=chi2;
sumchi2 +=chi2*chi2;
sum++;
}
}
Float_t mdensity = 0.9;
Float_t meann = 130;
Float_t meanchi = 1;
Float_t sdensity = 0.1;
Float_t smeann = 10;
Float_t smeanchi =0.4;
if (sum>20){
mdensity = sumdens/sum;
meann = sumn/sum;
meanchi = sumchi/sum;
sdensity = sumdens2/sum-mdensity*mdensity;
if (sdensity >= 0)
sdensity = TMath::Sqrt(sdensity);
else
sdensity = 0.1;
smeann = sumn2/sum-meann*meann;
if (smeann >= 0)
smeann = TMath::Sqrt(smeann);
else
smeann = 10;
smeanchi = sumchi2/sum - meanchi*meanchi;
if (smeanchi >= 0)
smeanchi = TMath::Sqrt(smeanchi);
else
smeanchi = 0.4;
}
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) {
continue;
}
if (pt->GetBSigned()) continue;
if (pt->GetBConstrain()) continue;
Bool_t isok =kFALSE;
if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
isok = kTRUE;
if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
isok =kTRUE;
if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
isok =kTRUE;
if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
isok =kTRUE;
if (isok)
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;
c->Use(10);
}
}
Double_t maxchi = meanchi+2.*smeanchi;
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) {
continue;
}
if (pt->GetBSigned()) continue;
Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
if (chi>maxchi) continue;
Float_t bfactor=1;
Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
minn=0;
if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
pt->SetBSigned(kTRUE);
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;
c->Use(10);
}
}
}
if (fDebug>0){
timer.Print();
}
}
Int_t AliTPCtracker::RefitInward(AliESDEvent *event)
{
if (!event) return 0;
const Int_t kMaxFriendTracks=2000;
fEvent = event;
fEventHLT = 0;
TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
if (!transform) {
AliFatal("Tranformations not in RefitInward");
return 0;
}
transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
Int_t nContribut = event->GetNumberOfTracks();
TGraphErrors * graphMultDependenceDeDx = 0x0;
if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
if (recoParam->GetUseTotCharge()) {
graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
} else {
graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
}
}
ReadSeeds(event,2);
fIteration=2;
PropagateForward2(fSeeds);
RemoveUsed2(fSeeds,0.4,0.4,20);
Int_t entriesSeed=fSeeds->GetEntries();
TObjArray arraySeed(entriesSeed);
for (Int_t i=0;i<entriesSeed;i++) {
arraySeed.AddAt(fSeeds->At(i),i);
}
SignShared(&arraySeed);
FindSplitted(fSeeds, event,2);
if ((AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC)>0) FindMultiMC(fSeeds, fEvent,2);
Int_t ntracks=0;
Int_t nseed = fSeeds->GetEntriesFast();
for (Int_t i=0;i<nseed;i++){
AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
if (!seed) continue;
if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(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() & kStreamRecoverIn)>0) {
(*fDebugStreamer)<<"RecoverIn"<<
"seed.="<<seed<<
"esd.="<<esd<<
"pin.="<<¶mIn<<
"pout.="<<¶mOut<<
"ncl="<<ncl<<
"\n";
}
if (ncl>15) {
seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
seed->SetNumberOfClusters(ncl);
}
}
seed->PropagateTo(fkParam->GetInnerRadiusLow());
seed->UpdatePoints();
AddCovariance(seed);
MakeESDBitmaps(seed, esd);
seed->CookdEdx(0.02,0.6);
CookLabel(seed,0.1);
if (((AliTPCReconstructor::StreamLevel()&kStreamRefitInward)>0) && seed!=0) {
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"RefitInward"<<
"Esd.="<<esd<<
"Track.="<<seed<<
"\n";
}
if (seed->GetNumberOfClusters()>15){
esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
esd->SetTPCPoints(seed->GetPoints());
esd->SetTPCPointsF(seed->GetNFoundable());
Int_t ndedx = seed->GetNCDEDX(0);
Float_t sdedx = seed->GetSDEDX(0);
Float_t dedx = seed->GetdEdx();
if (graphMultDependenceDeDx) {
Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
dedx += (1 - corrGain)*50.;
}
esd->SetTPCsignal(dedx, sdedx, ndedx);
Double32_t signal[4];
Double32_t signalMax[4];
Char_t ncl[3];
Char_t nrows[3];
for(Int_t iarr=0;iarr<3;iarr++) {
signal[iarr] = seed->GetDEDXregion(iarr+1);
signalMax[iarr] = seed->GetDEDXregion(iarr+5);
ncl[iarr] = seed->GetNCDEDX(iarr+1);
nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
}
signal[3] = seed->GetDEDXregion(4);
signalMax[3] = seed->GetDEDXregion(8);
AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
infoTpcPid->SetTPCSignalsQmax(signalMax);
esd->SetTPCdEdxInfo(infoTpcPid);
Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
esd->AddCalibObject(seedCopy);
}
ntracks++;
}
else{
}
}
if ((AliTPCReconstructor::StreamLevel()&kStreamClDump)>0) DumpClusters(2,fSeeds);
Info("RefitInward","Number of refitted tracks %d",ntracks);
AliCosmicTracker::FindCosmic(event, kTRUE);
FillClusterOccupancyInfo();
return 0;
}
Int_t AliTPCtracker::PropagateBack(AliESDEvent *event)
{
if (!event) return 0;
fEvent = event;
fEventHLT = 0;
fIteration = 1;
ReadSeeds(event,1);
PropagateBack(fSeeds);
RemoveUsed2(fSeeds,0.4,0.4,20);
FindSplitted(fSeeds, event,1);
if ((AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC)>0) FindMultiMC(fSeeds, fEvent,1);
Int_t nseed = fSeeds->GetEntriesFast();
Int_t ntracks=0;
for (Int_t i=0;i<nseed;i++){
AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
if (!seed) continue;
if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed);
seed->UpdatePoints();
AddCovariance(seed);
AliESDtrack *esd=event->GetTrack(i);
if (!esd) continue;
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()&kStreamRecoverBack)>0) {
(*fDebugStreamer)<<"RecoverBack"<<
"seed.="<<seed<<
"esd.="<<esd<<
"pin.="<<¶mIn<<
"pout.="<<¶mOut<<
"ncl="<<ncl<<
"\n";
}
if (ncl>15) {
seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
seed->SetNumberOfClusters(ncl);
}
}
seed->CookdEdx(0.02,0.6);
CookLabel(seed,0.1);
if (seed->GetNumberOfClusters()>15){
esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
esd->SetTPCPoints(seed->GetPoints());
esd->SetTPCPointsF(seed->GetNFoundable());
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();
if (((AliTPCReconstructor::StreamLevel()&kStreamCPropagateBack)>0) && esd) {
(*fDebugStreamer)<<"PropagateBack"<<
"Tr0.="<<seed<<
"esd.="<<esd<<
"EventNrInFile="<<eventnumber<<
"\n";
}
}
}
if (AliTPCReconstructor::StreamLevel()&kStreamClDump) DumpClusters(1,fSeeds);
Info("PropagateBack","Number of back propagated tracks %d",ntracks);
fEvent =0;
fEventHLT = 0;
return 0;
}
Int_t AliTPCtracker::PostProcess(AliESDEvent *event)
{
if (!event) return 0;
if(IsTPCHVDipEvent(event)) {
event->ResetDetectorStatus(AliDAQ::kTPC);
}
return 0;
}
void AliTPCtracker::DeleteSeeds()
{
fSeeds->Clear();
ResetSeedsPool();
delete fSeeds;
fSeeds =0;
}
void AliTPCtracker::ReadSeeds(const AliESDEvent *const event, Int_t direction)
{
Int_t nentr=event->GetNumberOfTracks();
if (fDebug>0){
Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
}
if (fSeeds)
DeleteSeeds();
if (!fSeeds){
fSeeds = new TObjArray(nentr);
}
UnsignClusters();
for (Int_t i=0; i<nentr; i++) {
AliESDtrack *esd=event->GetTrack(i);
ULong_t status=esd->GetStatus();
if (!(status&AliESDtrack::kTPCin)) continue;
AliTPCtrack t(*esd);
t.SetNumberOfClusters(0);
AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t);
seed->SetPoolID(fLastSeedID);
seed->SetUniqueID(esd->GetID());
AddCovariance(seed);
for (Int_t ikink=0;ikink<3;ikink++) {
Int_t index = esd->GetKinkIndex(ikink);
seed->GetKinkIndexes()[ikink] = index;
if (index==0) continue;
index = TMath::Abs(index);
AliESDkink * kink = fEvent->GetKink(index-1);
if (kink&&esd->GetKinkIndex(ikink)<0){
if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
}
if (kink&&esd->GetKinkIndex(ikink)>0){
if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
}
}
if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
fSectors=fInnerSec; fN=fkNIS;
Double_t alpha=seed->GetAlpha();
if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
if (alpha < 0. ) alpha += 2.*TMath::Pi();
Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
alpha-=seed->GetAlpha();
if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
if (TMath::Abs(alpha) > 0.001) {
AliWarning(Form("Rotating track over %f",alpha));
if (!seed->Rotate(alpha)) {
MarkSeedFree( seed );
continue;
}
}
seed->SetESD(esd);
if (esd->GetKinkIndex(0)<=0){
for (Int_t irow=0;irow<160;irow++){
Int_t index = seed->GetClusterIndex2(irow);
if (index >= 0){
AliTPCclusterMI * cl = GetClusterMI(index);
seed->SetClusterPointer(irow,cl);
if (cl){
if ((index & 0x8000)==0){
cl->Use(10);
}else{
cl->Use(6);
}
}else{
Info("ReadSeeds","Not found cluster");
}
}
}
}
fSeeds->AddAt(seed,i);
}
}
void AliTPCtracker::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
Float_t deltay, Int_t ddsec) {
Int_t nin0 = 0;
Int_t nin1 = 0;
Int_t nin2 = 0;
Int_t nin = 0;
Int_t nout1 = 0;
Int_t nout2 = 0;
Double_t x[5], c[15];
AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed();
seed->SetPoolID(fLastSeedID);
Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
Double_t cs=cos(alpha), sn=sin(alpha);
Double_t x1 =GetXrow(i1);
Double_t xx2=GetXrow(i2);
Double_t x3=GetX(), y3=GetY(), z3=GetZ();
Int_t imiddle = (i2+i1)/2;
Double_t xm = GetXrow(imiddle);
const AliTPCtrackerRow& krm=GetRow(sec,imiddle);
Int_t ns =sec;
const AliTPCtrackerRow& kr1=GetRow(ns,i1);
Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
if (dvertexmax*0.5*cuts[0]>0.85){
cuts[0] = 0.85/(dvertexmax*0.5+1.);
}
Double_t r2min = 1/(cuts[0]*cuts[0]);
if (deltay>0) ddsec = 0;
for (Int_t is=0; is < kr1; is++) {
if (kr1[is]->IsUsed(10)) continue;
if (kr1[is]->IsDisabled()) {
continue;
}
Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;
Float_t anglez = (z1-z3)/(x1-x3);
Float_t extraz = z1 - anglez*(x1-xx2);
Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
Double_t dvertex = TMath::Sqrt(dvertex2);
Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
Int_t dsec1=-ddsec;
Int_t dsec2= ddsec;
if (y1<0) dsec2= 0;
if (y1>0) dsec1= 0;
Double_t dddz1=0;
Double_t dddz2=0;
if ( (z1-z3)>0)
dddz1 =1;
else
dddz2 =1;
for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
Int_t sec2 = sec + dsec;
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);
Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
Double_t x2, y2, z2;
Double_t dxx0 = (xx2-x3)*cs13r;
Double_t dyy0 = (xx2-x3)*sn13r;
for (Int_t js=index1; js < index2; js++) {
const AliTPCclusterMI *kcl = kr2[js];
if (kcl->IsUsed(10)) continue;
if (kcl->IsDisabled()) {
continue;
}
Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
if (TMath::Abs(yy0)<0.000001) continue;
Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
Double_t r02 = (0.25+y0*y0)*dvertex2;
if (r02<r2min) continue;
nin0++;
Double_t c0 = 1/TMath::Sqrt(r02);
if (yy0>0) c0*=-1.;
Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
Double_t z0 = kcl->GetZ();
Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
if (TMath::Abs(zzzz2-z0)>0.5) continue;
nin1++;
Double_t dip = (z1-z0)*c0/dfi1;
Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
y2 = kcl->GetY();
if (dsec==0){
x2 = xx2;
z2 = kcl->GetZ();
}
else
{
z2 = kcl->GetZ();
x2= xx2*cs-y2*sn*dsec;
y2=+xx2*sn*dsec+y2*cs;
}
x[0] = y1;
x[1] = z1;
x[2] = x0;
x[3] = dip;
x[4] = c0;
Double_t ym,zm;
GetProlongation(x1,xm,x,ym,zm);
UInt_t dummy;
AliTPCclusterMI * cm=0;
if (TMath::Abs(ym)-ymaxm<0){
cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
if ((!cm) || (cm->IsUsed(10))) {
continue;
}
}
else{
Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
Double_t xr2 = x0*cs+yr1*sn*dsec;
Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
GetProlongation(xx2,xm,xr,ym,zm);
if (TMath::Abs(ym)-ymaxm<0){
cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
if ((!cm) || (cm->IsUsed(10))) {
continue;
}
}
}
nin2++;
Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
c[0]=sy1;
c[1]=0.; c[2]=sz1;
c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
c[13]=f30*sy1*f40+f32*sy2*f42;
c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
UInt_t index=kr1.GetIndex(is);
if (seed) {MarkSeedFree(seed); seed = 0;}
AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, ns*alpha+shift, x, c, index);
seed->SetPoolID(fLastSeedID);
track->SetIsSeeding(kTRUE);
track->SetSeed1(i1);
track->SetSeed2(i2);
track->SetSeedType(3);
FollowProlongation(*track, (i1+i2)/2,1);
Int_t foundable,found,shared;
track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
MarkSeedFree(seed); seed = 0;
continue;
}
nin++;
FollowProlongation(*track, i2,1);
track->SetBConstrain(1);
track->SetLastPoint(i1);
track->SetFirstPoint(track->GetLastPoint());
if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
MarkSeedFree(seed); seed = 0;
continue;
}
nout1++;
Double_t zv, bz=GetBz();
if ( !track->GetZAt(0.,bz,zv) ) continue;
if (TMath::Abs(zv-z3)>cuts[2]) {
FollowProlongation(*track, TMath::Max(i2-20,0));
if ( !track->GetZAt(0.,bz,zv) ) continue;
if (TMath::Abs(zv-z3)>cuts[2]){
FollowProlongation(*track, TMath::Max(i2-40,0));
if ( !track->GetZAt(0.,bz,zv) ) continue;
if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
FollowProlongation(*track2, i2,1);
track2->SetBConstrain(kFALSE);
track2->SetSeedType(1);
arr->AddLast(track2);
MarkSeedFree( seed ); seed = 0;
continue;
}
else{
MarkSeedFree( seed ); seed = 0;
continue;
}
}
}
track->SetSeedType(0);
arr->AddLast(track);
seed = new( NextFreeSeed() ) AliTPCseed;
seed->SetPoolID(fLastSeedID);
nout2++;
if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
break;
}
}
}
if (fDebug>3){
Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
}
if (seed) MarkSeedFree( seed );
}
void AliTPCtracker::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
Float_t deltay) {
Int_t nin0 = 0;
Int_t nin1 = 0;
Int_t nin2 = 0;
Int_t nin = 0;
Int_t nout1 = 0;
Int_t nout2 = 0;
Int_t nout3 =0;
Double_t x[5], c[15];
AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
seed->SetPoolID(fLastSeedID);
Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
Double_t x1 = GetXrow(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 AliTPCtrackerRow& kr1p=GetRow(sec,i1);
Double_t x1m = GetXrow(i1-2);
const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
Double_t x3 = GetXrow(i1-7);
AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
Double_t x3p = GetXrow(i1-6);
AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
Double_t x3m = GetXrow(i1-8);
Int_t im = i1-4;
Double_t xm = GetXrow(im);
const AliTPCtrackerRow& krm=GetRow(sec,im);
Double_t deltax = x1-x3;
Double_t dymax = deltax*cuts[1];
Double_t dzmax = deltax*cuts[3];
for (Int_t is=0; is < kr1; is++) {
if (kr1[is]->IsUsed(10)) continue;
if (kr1[is]->IsDisabled()) {
continue;
}
Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;
Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
Double_t y3, z3;
UInt_t index;
for (Int_t js=index1; js < index2; js++) {
const AliTPCclusterMI *kcl = kr3[js];
if (kcl->IsDisabled()) {
continue;
}
if (kcl->IsUsed(10)) continue;
y3 = kcl->GetY();
if (TMath::Abs(y1-y3)>dymax) continue;
z3 = kcl->GetZ();
if (TMath::Abs(z1-z3)>dzmax) continue;
Double_t angley = (y1-y3)/(x1-x3);
Double_t anglez = (z1-z3)/(x1-x3);
Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
Double_t yyym = angley*(xm-x1)+y1;
Double_t zzzm = anglez*(xm-x1)+z1;
const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
if (!kcm) continue;
if (kcm->IsUsed(10)) continue;
if (kcm->IsDisabled()) {
continue;
}
erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
Int_t used =0;
Int_t found =0;
const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
anglez*(x1m-x1)+z1,
erry,errz,index);
if (kc1m){
found++;
if (kc1m->IsUsed(10)) used++;
}
const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
anglez*(x1p-x1)+z1,
erry,errz,index);
if (kc1p){
found++;
if (kc1p->IsUsed(10)) used++;
}
if (used>1) continue;
if (found<1) continue;
const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
anglez*(x3m-x3)+z3,
erry,errz,index);
if (kc3m){
found++;
if (kc3m->IsUsed(10)) used++;
}
else
continue;
const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
anglez*(x3p-x3)+z3,
erry,errz,index);
if (kc3p){
found++;
if (kc3p->IsUsed(10)) used++;
}
else
continue;
if (used>1) continue;
if (found<3) continue;
Double_t x2,y2,z2;
x2 = xm;
y2 = kcm->GetY();
z2 = kcm->GetZ();
x[0]=y1;
x[1]=z1;
x[4]=F1(x1,y1,x2,y2,x3,y3);
nin0++;
x[2]=F2(x1,y1,x2,y2,x3,y3);
nin1++;
x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
nin2++;
Double_t sy1=0.1, sz1=0.1;
Double_t sy2=0.1, sz2=0.1;
Double_t sy3=0.1, sy=0.1, sz=0.1;
Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
c[0]=sy1;
c[1]=0.; c[2]=sz1;
c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
c[13]=f30*sy1*f40+f32*sy2*f42;
c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
index=kr1.GetIndex(is);
if (seed) {MarkSeedFree( seed ); seed = 0;}
AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, sec*alpha+shift, x, c, index);
seed->SetPoolID(fLastSeedID);
track->SetIsSeeding(kTRUE);
nin++;
FollowProlongation(*track, i1-7,1);
if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
MarkSeedFree( seed ); seed = 0;
continue;
}
nout1++;
nout2++;
FollowProlongation(*track, i2,1);
track->SetBConstrain(0);
track->SetLastPoint(i1+fInnerSec->GetNRows());
track->SetFirstPoint(track->GetLastPoint());
if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
MarkSeedFree( seed ); seed = 0;
continue;
}
{
FollowProlongation(*track, TMath::Max(i2-10,0),1);
AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
FollowProlongation(*track2, i2,1);
track2->SetBConstrain(kFALSE);
track2->SetSeedType(4);
arr->AddLast(track2);
MarkSeedFree( seed ); seed = 0;
}
nout3++;
}
}
if (fDebug>3){
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);
}
if (seed) MarkSeedFree(seed);
}
void AliTPCtracker::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t *,
Float_t deltay, Bool_t ) {
Int_t nin0=0;
Int_t nin1=0;
Int_t nin2=0;
Int_t nin3=0;
Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
Int_t row0 = (i1+i2)/2;
Int_t drow = (i1-i2)/2;
const AliTPCtrackerRow& kr0=fSectors[sec][row0];
AliTPCtrackerRow * kr=0;
AliTPCpolyTrack polytrack;
Int_t nclusters=fSectors[sec][row0];
AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
seed->SetPoolID(fLastSeedID);
Int_t sumused=0;
Int_t cused=0;
Int_t cnused=0;
for (Int_t is=0; is < nclusters; is++) {
Int_t nfound =0;
Int_t nfoundable =0;
for (Int_t iter =1; iter<2; iter++){
const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
const AliTPCclusterMI * cl= kr0[is];
if (cl->IsDisabled()) {
continue;
}
if (cl->IsUsed(10)) {
cused++;
}
else{
cnused++;
}
Double_t x = kr0.GetX();
nfound =0;
nfoundable =0;
polytrack.Reset();
Double_t y0= cl->GetY();
Double_t z0= cl->GetZ();
Float_t erry = 0;
Float_t errz = 0;
Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;
erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
polytrack.AddPoint(x,y0,z0,erry, errz);
sumused=0;
if (cl->IsUsed(10)) sumused++;
Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
x = krm.GetX();
AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
if (cl1->IsUsed(10)) sumused++;
polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
}
x = krp.GetX();
AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
if (cl2) {
erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
if (cl2->IsUsed(10)) sumused++;
polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
}
if (sumused>0) continue;
nin0++;
polytrack.UpdateParameters();
roadz = 1.2;
roady = 1.2;
Double_t yn,zn;
nfoundable = polytrack.GetN();
nfound = nfoundable;
for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
Float_t maxdist = 0.8*(1.+3./(ddrow));
for (Int_t delta = -1;delta<=1;delta+=2){
Int_t row = row0+ddrow*delta;
kr = &(fSectors[sec][row]);
Double_t xn = kr->GetX();
Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
polytrack.GetFitPoint(xn,yn,zn);
if (TMath::Abs(yn)>ymax1) continue;
nfoundable++;
AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
if (cln) {
Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
if (dist<maxdist){
erry=0.1;
errz=0.1;
polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
nfound++;
}
}
}
if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
polytrack.UpdateParameters();
}
}
if ( (sumused>3) || (sumused>0.5*nfound)) {
continue;
}
nin1++;
Double_t dy,dz;
polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
AliTPCpolyTrack track2;
polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
if (track2.GetN()<0.5*nfoundable) continue;
nin2++;
if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
for (Int_t constrain=0; constrain<=0;constrain++){
Double_t x[5], c[15];
Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
track2.GetBoundaries(x3,x1);
x2 = (x1+x3)/2.;
track2.GetFitPoint(x1,y1,z1);
track2.GetFitPoint(x2,y2,z2);
track2.GetFitPoint(x3,y3,z3);
Double_t x0,y0,z0;
x0=0;
polytrack.GetFitPoint(x0,y0,z0);
if (constrain) {
x2 = x3;
y2 = y3;
z2 = z3;
x3 = 0;
y3 = 0;
z3 = 0;
}
x[0]=y1;
x[1]=z1;
x[4]=F1(x1,y1,x2,y2,x3,y3);
x[2]=F2(x1,y1,x2,y2,x3,y3);
x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
Double_t sy =0.1, sz =0.1;
Double_t sy1=0.02, sz1=0.02;
Double_t sy2=0.02, sz2=0.02;
Double_t sy3=0.02;
if (constrain){
sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
}
Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
c[0]=sy1;
c[1]=0.; c[2]=sz1;
c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
c[13]=f30*sy1*f40+f32*sy2*f42;
c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
Int_t row1 = GetRowNumber(x1);
UInt_t index=0;
if (seed) {MarkSeedFree( seed ); seed = 0;}
AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1,sec*alpha+shift,x,c,index);
seed->SetPoolID(fLastSeedID);
track->SetIsSeeding(kTRUE);
Int_t rc=FollowProlongation(*track, i2);
if (constrain) track->SetBConstrain(1);
else
track->SetBConstrain(0);
track->SetLastPoint(row1+fInnerSec->GetNRows());
track->SetFirstPoint(track->GetLastPoint());
if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
track->GetNShared()>0.4*track->GetNumberOfClusters()) {
MarkSeedFree( seed ); seed = 0;
}
else {
arr->AddLast(track);
seed = new( NextFreeSeed() ) AliTPCseed;
seed->SetPoolID(fLastSeedID);
}
nin3++;
}
}
}
if (fDebug>3){
Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
}
if (seed) MarkSeedFree( seed );
}
AliTPCseed *AliTPCtracker::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
{
Int_t p0 = int(r0*track->GetNumberOfClusters());
Int_t p1 = int(r1*track->GetNumberOfClusters());
Int_t p2 = int(r2*track->GetNumberOfClusters());
Int_t pp2=0;
Double_t x0[3],x1[3],x2[3];
for (Int_t i=0;i<3;i++){
x0[i]=-1;
x1[i]=-1;
x2[i]=-1;
}
Int_t sec0=0, sec1=0, sec2=0;
Int_t index=-1;
Int_t clindex;
for (Int_t i=0;i<160;i++){
if (track->GetClusterPointer(i)){
index++;
AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
if ( (index<p0) || x0[0]<0 ){
if (trpoint->GetX()>1){
clindex = track->GetClusterIndex2(i);
if (clindex >= 0){
x0[0] = trpoint->GetX();
x0[1] = trpoint->GetY();
x0[2] = trpoint->GetZ();
sec0 = ((clindex&0xff000000)>>24)%18;
}
}
}
if ( (index<p1) &&(trpoint->GetX()>1)){
clindex = track->GetClusterIndex2(i);
if (clindex >= 0){
x1[0] = trpoint->GetX();
x1[1] = trpoint->GetY();
x1[2] = trpoint->GetZ();
sec1 = ((clindex&0xff000000)>>24)%18;
}
}
if ( (index<p2) &&(trpoint->GetX()>1)){
clindex = track->GetClusterIndex2(i);
if (clindex >= 0){
x2[0] = trpoint->GetX();
x2[1] = trpoint->GetY();
x2[2] = trpoint->GetZ();
sec2 = ((clindex&0xff000000)>>24)%18;
pp2 = i;
}
}
}
}
Double_t alpha, cs,sn, xx2,yy2;
alpha = (sec1-sec2)*fSectors->GetAlpha();
cs = TMath::Cos(alpha);
sn = TMath::Sin(alpha);
xx2= x1[0]*cs-x1[1]*sn;
yy2= x1[0]*sn+x1[1]*cs;
x1[0] = xx2;
x1[1] = yy2;
alpha = (sec0-sec2)*fSectors->GetAlpha();
cs = TMath::Cos(alpha);
sn = TMath::Sin(alpha);
xx2= x0[0]*cs-x0[1]*sn;
yy2= x0[0]*sn+x0[1]*cs;
x0[0] = xx2;
x0[1] = yy2;
Double_t x[5],c[15];
x[0]=x2[1];
x[1]=x2[2];
x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
Double_t sy =0.1, sz =0.1;
Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
Double_t sy3=0.01+track->GetSigmaY2();
Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
c[0]=sy1;
c[1]=0.; c[2]=sz1;
c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
c[13]=f30*sy1*f40+f32*sy2*f42;
c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
seed->SetPoolID(fLastSeedID);
seed->SetLastPoint(pp2);
seed->SetFirstPoint(pp2);
return seed;
}
AliTPCseed *AliTPCtracker::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
{
Int_t nclusters = 0;
for (Int_t irow=0;irow<160;irow++){
if (track->GetClusterIndex(irow)>0) nclusters++;
}
Int_t ipos[3];
ipos[0] = TMath::Max(int(r0*nclusters),0);
ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1);
ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1);
Double_t xyz[3][3]={{0}};
Int_t row[3]={0},sec[3]={0,0,0};
Int_t index=-1;
for (Int_t irow=0;irow<160;irow++){
if (track->GetClusterIndex2(irow)<0) continue;
index++;
for (Int_t ipoint=0;ipoint<3;ipoint++){
if (index<=ipos[ipoint]) row[ipoint] = irow;
}
}
for (Int_t ipoint=0;ipoint<3;ipoint++){
Int_t clindex = track->GetClusterIndex2(row[ipoint]);
AliTPCclusterMI * cl = GetClusterMI(clindex);
if (cl==0) {
return 0;
}
sec[ipoint] = ((clindex&0xff000000)>>24)%18;
xyz[ipoint][0] = GetXrow(row[ipoint]);
xyz[ipoint][1] = cl->GetY();
xyz[ipoint][2] = cl->GetZ();
}
Double_t alpha, cs,sn, xx2,yy2;
alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
cs = TMath::Cos(alpha);
sn = TMath::Sin(alpha);
xx2= xyz[1][0]*cs-xyz[1][1]*sn;
yy2= xyz[1][0]*sn+xyz[1][1]*cs;
xyz[1][0] = xx2;
xyz[1][1] = yy2;
alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
cs = TMath::Cos(alpha);
sn = TMath::Sin(alpha);
xx2= xyz[0][0]*cs-xyz[0][1]*sn;
yy2= xyz[0][0]*sn+xyz[0][1]*cs;
xyz[0][0] = xx2;
xyz[0][1] = yy2;
Double_t x[5],c[15];
x[0]=xyz[2][1];
x[1]=xyz[2][2];
x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
Double_t sy =0.1, sz =0.1;
Double_t sy1=0.2, sz1=0.2;
Double_t sy2=0.2, sz2=0.2;
Double_t sy3=0.2;
Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
c[0]=sy1;
c[1]=0.; c[2]=sz1;
c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
c[13]=f30*sy1*f40+f32*sy2*f42;
c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
seed->SetPoolID(fLastSeedID);
seed->SetLastPoint(row[2]);
seed->SetFirstPoint(row[2]);
return seed;
}
AliTPCseed *AliTPCtracker::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
{
Double_t xyz[3][3];
Int_t row[3]={0,0,0};
Int_t sec[3]={0,0,0};
if (forward){
for (Int_t irow=r0;irow<160;irow++){
if (track->GetClusterIndex(irow)>0){
row[0] = irow;
break;
}
}
for (Int_t irow=160;irow>r0;irow--){
if (track->GetClusterIndex(irow)>0){
row[2] = irow;
break;
}
}
for (Int_t irow=row[2]-15;irow>row[0];irow--){
if (track->GetClusterIndex(irow)>0){
row[1] = irow;
break;
}
}
}
if (!forward){
for (Int_t irow=0;irow<r0;irow++){
if (track->GetClusterIndex(irow)>0){
row[0] = irow;
break;
}
}
for (Int_t irow=r0;irow>0;irow--){
if (track->GetClusterIndex(irow)>0){
row[2] = irow;
break;
}
}
for (Int_t irow=row[2]-15;irow>row[0];irow--){
if (track->GetClusterIndex(irow)>0){
row[1] = irow;
break;
}
}
}
if ((row[2]-row[0])<20) return 0;
if (row[1]==0) return 0;
for (Int_t ipoint=0;ipoint<3;ipoint++){
Int_t clindex = track->GetClusterIndex2(row[ipoint]);
AliTPCclusterMI * cl = GetClusterMI(clindex);
if (cl==0) {
return 0;
}
sec[ipoint] = ((clindex&0xff000000)>>24)%18;
xyz[ipoint][0] = GetXrow(row[ipoint]);
AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
if (point&&ipoint<2){
xyz[ipoint][1] = point->GetY();
xyz[ipoint][2] = point->GetZ();
}
else{
xyz[ipoint][1] = cl->GetY();
xyz[ipoint][2] = cl->GetZ();
}
}
Double_t alpha, cs,sn, xx2,yy2;
alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
cs = TMath::Cos(alpha);
sn = TMath::Sin(alpha);
xx2= xyz[1][0]*cs-xyz[1][1]*sn;
yy2= xyz[1][0]*sn+xyz[1][1]*cs;
xyz[1][0] = xx2;
xyz[1][1] = yy2;
alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
cs = TMath::Cos(alpha);
sn = TMath::Sin(alpha);
xx2= xyz[0][0]*cs-xyz[0][1]*sn;
yy2= xyz[0][0]*sn+xyz[0][1]*cs;
xyz[0][0] = xx2;
xyz[0][1] = yy2;
Double_t x[5],c[15];
x[0]=xyz[2][1];
x[1]=xyz[2][2];
x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
Double_t sy =0.1, sz =0.1;
Double_t sy1=0.2, sz1=0.2;
Double_t sy2=0.2, sz2=0.2;
Double_t sy3=0.2;
Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
c[0]=sy1;
c[1]=0.; c[2]=sz1;
c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
c[13]=f30*sy1*f40+f32*sy2*f42;
c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
seed->SetPoolID(fLastSeedID);
seed->SetLastPoint(row[2]);
seed->SetFirstPoint(row[2]);
for (Int_t i=row[0];i<row[2];i++){
seed->SetClusterIndex(i, track->GetClusterIndex(i));
}
return seed;
}
void AliTPCtracker::FindMultiMC(const TObjArray * array, AliESDEvent *, Int_t iter)
{
const Float_t kMaxdPhi = 0.2;
Int_t nentries = array->GetEntriesFast();
AliHelix *helixes = new AliHelix[nentries];
Float_t *xm = new Float_t[nentries];
Float_t *dz0 = new Float_t[nentries];
Float_t *dz1 = new Float_t[nentries];
TStopwatch timer;
timer.Start();
for (Int_t i=0;i<nentries;i++){
AliTPCseed* track = (AliTPCseed*)array->At(i);
if (!track) continue;
track->SetCircular(0);
new (&helixes[i]) AliHelix(*track);
Int_t ncl=0;
xm[i]=0;
Float_t dz[2];
track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
dz0[i]=dz[0];
dz1[i]=dz[1];
for (Int_t icl=0; icl<160; icl++){
AliTPCclusterMI * cl = track->GetClusterPointer(icl);
if (cl) {
xm[i]+=cl->GetX();
ncl++;
}
}
if (ncl>0) xm[i]/=Float_t(ncl);
}
for (Int_t i0=0;i0<nentries;i0++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
if (!track0) continue;
Float_t xc0 = helixes[i0].GetHelix(6);
Float_t yc0 = helixes[i0].GetHelix(7);
Float_t r0 = helixes[i0].GetHelix(8);
Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
Float_t fi0 = TMath::ATan2(yc0,xc0);
for (Int_t i1=i0+1;i1<nentries;i1++){
AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
if (!track1) continue;
Int_t lab0=track0->GetLabel();
Int_t lab1=track1->GetLabel();
if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
Float_t xc1 = helixes[i1].GetHelix(6);
Float_t yc1 = helixes[i1].GetHelix(7);
Float_t r1 = helixes[i1].GetHelix(8);
Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
Float_t fi1 = TMath::ATan2(yc1,xc1);
Float_t dfi = fi0-fi1;
if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi();
if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi();
if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
fi1 = -TMath::ATan2(yc1,-xc1);
dfi = fi0-fi1;
}
Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
Double_t dist[3];
Double_t mdist[3]={0,0,0};
track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
Float_t sum =0;
Float_t sums=0;
for (Int_t icl=0; icl<160; icl++){
AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
if (cl0&&cl1) {
sum++;
if (cl0==cl1) sums++;
}
}
if ((AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC)>0) {
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Multi"<<
"iter="<<iter<<
"lab0="<<lab0<<
"lab1="<<lab1<<
"Tr0.="<<track0<<
"Tr1.="<<track1<<
"h0.="<<&helixes[i0]<<
"h1.="<<&helixes[i1]<<
"sum="<<sum<<
"sums="<<sums<<
"xm0="<<xm[i0]<<
"xm1="<<xm[i1]<<
"dfi="<<dfi<<
"dtheta="<<dtheta<<
"dz00="<<dz0[i0]<<
"dz01="<<dz0[i1]<<
"dz10="<<dz1[i1]<<
"dz11="<<dz1[i1]<<
"dist0="<<dist[0]<<
"dist1="<<dist[1]<<
"dist2="<<dist[2]<<
"mdist0="<<mdist[0]<<
"mdist1="<<mdist[1]<<
"mdist2="<<mdist[2]<<
"r0="<<r0<<
"rc0="<<rc0<<
"fi0="<<fi0<<
"fi1="<<fi1<<
"r1="<<r1<<
"rc1="<<rc1<<
"\n";
}
}
}
delete [] helixes;
delete [] xm;
delete [] dz0;
delete [] dz1;
if (AliTPCReconstructor::StreamLevel()>0) {
AliInfo("Time for curling tracks removal DEBUGGING MC");
timer.Print();
}
}
void AliTPCtracker::FindSplitted(TObjArray * array, AliESDEvent *, Int_t ){
const Double_t xref=GetXrow(63);
const Double_t kCutP1=10;
const Double_t kCutP2=0.15;
const Double_t kCutP3=0.15;
const Double_t kCutAlpha=0.15;
Int_t firstpoint = 0;
Int_t lastpoint = 160;
Int_t nentries = array->GetEntriesFast();
AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
TStopwatch timer;
timer.Start();
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; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
if (!pt){
quality[i]=-1;
continue;
}
pt->UpdatePoints();
Float_t * points = pt->GetPoints();
if (points[3]<0.8) quality[i] =-1;
quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
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);
for (Int_t i0=0; i0<nseed; i0++) {
Int_t index0=indexes[i0];
if (!(array->UncheckedAt(index0))) continue;
AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
if (!s1->IsActive()) continue;
AliExternalTrackParam &par0=params[index0];
for (Int_t i1=i0+1; i1<nseed; i1++) {
Int_t index1=indexes[i1];
if (!(array->UncheckedAt(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;i<lastpoint;i++){
if (s1->GetClusterIndex2(i)>0) nall0++;
if (s2->GetClusterIndex2(i)>0) nall1++;
if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
if (i<firstRow) firstRow=i;
if (i>lastRow) lastRow=i;
}
if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
if (i<firstShared) firstShared=i;
if (i>lastShared) lastShared=i;
sumShared++;
}
}
Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
if ((AliTPCReconstructor::StreamLevel()&kStreamSplitted2)>0){
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="<<fIteration<<
"lab0="<<lab0<<
"lab1="<<lab1<<
"index0="<<index0<<
"index1="<<index1<<
"ratio0="<<ratio0<<
"ratio1="<<ratio1<<
"p0.="<<&par0<<
"p1.="<<&par1<<
"s0.="<<s1<<
"s1.="<<s2<<
"n0="<<n0<<
"n1="<<n1<<
"nall0="<<nall0<<
"nall1="<<nall1<<
"n0F="<<n0F<<
"n1F="<<n1F<<
"shared="<<sumShared<<
"firstS="<<firstShared<<
"lastS="<<lastShared<<
"firstRow="<<firstRow<<
"lastRow="<<lastRow<<
"\n";
}
if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
MarkSeedFree( array->RemoveAt(index1) );
}
}
}
delete [] params;
delete [] quality;
delete [] indexes;
}
void AliTPCtracker::FindCurling(const TObjArray * array, AliESDEvent *, Int_t iter)
{
const Float_t kMaxC = 400;
const Float_t kMaxdTheta = 0.15;
const Float_t kMaxdPhi = 0.15;
const Float_t kPtRatio = 0.3;
const Float_t kMinDCAR = 2.;
const Float_t kMaxDeltaRMax = 40;
const Float_t kMaxDeltaRMin = 5.;
const Float_t kMinAngle = 2.9;
const Float_t kMaxDist = 5;
Int_t nentries = array->GetEntriesFast();
AliHelix *helixes = new AliHelix[nentries];
for (Int_t i=0;i<nentries;i++){
AliTPCseed* track = (AliTPCseed*)array->At(i);
if (!track) continue;
track->SetCircular(0);
new (&helixes[i]) AliHelix(*track);
}
TStopwatch timer;
timer.Start();
Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
for (Int_t i0=0;i0<nentries;i0++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
if (!track0) continue;
if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
Float_t xc0 = helixes[i0].GetHelix(6);
Float_t yc0 = helixes[i0].GetHelix(7);
Float_t r0 = helixes[i0].GetHelix(8);
Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
Float_t fi0 = TMath::ATan2(yc0,xc0);
for (Int_t i1=i0+1;i1<nentries;i1++){
AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
if (!track1) continue;
if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
Float_t xc1 = helixes[i1].GetHelix(6);
Float_t yc1 = helixes[i1].GetHelix(7);
Float_t r1 = helixes[i1].GetHelix(8);
Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
Float_t fi1 = TMath::ATan2(yc1,xc1);
Float_t dfi = fi0-fi1;
if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi();
if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi();
Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue;
if ( TMath::Abs(dfi)>kMaxdPhi) continue;
if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue;
Float_t pt0 = track0->GetSignedPt();
Float_t pt1 = track1->GetSignedPt();
if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
if (npoints==0) continue;
helixes[i0].GetClosestPhases(helixes[i1], phase);
Double_t xyz0[3];
Double_t xyz1[3];
Double_t hangles[3];
helixes[i0].Evaluate(phase[0][0],xyz0);
helixes[i1].Evaluate(phase[0][1],xyz1);
helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
Double_t deltah[2],deltabest;
if (TMath::Abs(hangles[2])<kMinAngle) continue;
if (npoints>0){
Int_t ibest=0;
helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
if (npoints==2){
helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
if (deltah[1]<deltah[0]) ibest=1;
}
deltabest = TMath::Sqrt(deltah[ibest]);
helixes[i0].Evaluate(phase[ibest][0],xyz0);
helixes[i1].Evaluate(phase[ibest][1],xyz1);
helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
Double_t radiusbest = TMath::Sqrt(radius[ibest]);
if (deltabest>kMaxDist) continue;
Bool_t sign =kFALSE;
if (hangles[2]>kMinAngle) sign =kTRUE;
if (sign){
if (track0->OneOverPt()<track1->OneOverPt()){
track0->SetCircular(track0->GetCircular()+1);
track1->SetCircular(track1->GetCircular()+2);
}
else{
track1->SetCircular(track1->GetCircular()+1);
track0->SetCircular(track0->GetCircular()+2);
}
}
if ((AliTPCReconstructor::StreamLevel()&kStreamFindCurling)>0){
Int_t lab0=track0->GetLabel();
Int_t lab1=track1->GetLabel();
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Curling2"<<
"iter="<<iter<<
"lab0="<<lab0<<
"lab1="<<lab1<<
"Tr0.="<<track0<<
"Tr1.="<<track1<<
"r0="<<r0<<
"rc0="<<rc0<<
"fi0="<<fi0<<
"r1="<<r1<<
"rc1="<<rc1<<
"fi1="<<fi1<<
"dfi="<<dfi<<
"dtheta="<<dtheta<<
"npoints="<<npoints<<
"hangles0="<<hangles[0]<<
"hangles1="<<hangles[1]<<
"hangles2="<<hangles[2]<<
"xyz0="<<xyz0[2]<<
"xyzz1="<<xyz1[2]<<
"radius="<<radiusbest<<
"deltabest="<<deltabest<<
"phase0="<<phase[ibest][0]<<
"phase1="<<phase[ibest][1]<<
"\n";
}
}
}
}
delete [] helixes;
if (AliTPCReconstructor::StreamLevel()>1) {
AliInfo("Time for curling tracks removal");
timer.Print();
}
}
void AliTPCtracker::FindKinks(TObjArray * array, AliESDEvent *esd)
{
TObjArray *kinks= new TObjArray(10000);
Int_t nentries = array->GetEntriesFast();
AliHelix *helixes = new AliHelix[nentries];
Int_t *sign = new Int_t[nentries];
Int_t *nclusters = new Int_t[nentries];
Float_t *alpha = new Float_t[nentries];
AliKink *kink = new AliKink();
Int_t * usage = new Int_t[nentries];
Float_t *zm = new Float_t[nentries];
Float_t *z0 = new Float_t[nentries];
Float_t *fim = new Float_t[nentries];
Float_t *shared = new Float_t[nentries];
Bool_t *circular = new Bool_t[nentries];
Float_t *dca = new Float_t[nentries];
for (Int_t i=0;i<nentries;i++){
sign[i]=0;
usage[i]=0;
AliTPCseed* track = (AliTPCseed*)array->At(i);
if (!track) continue;
track->SetCircular(0);
shared[i] = kFALSE;
track->UpdatePoints();
if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
}
nclusters[i]=track->GetNumberOfClusters();
alpha[i] = track->GetAlpha();
new (&helixes[i]) AliHelix(*track);
Double_t xyz[3];
helixes[i].Evaluate(0,xyz);
sign[i] = (track->GetC()>0) ? -1:1;
Double_t x,y,z;
x=160;
if (track->GetProlongation(x,y,z)){
zm[i] = z;
fim[i] = alpha[i]+TMath::ATan2(y,x);
}
else{
zm[i] = track->GetZ();
fim[i] = alpha[i];
}
z0[i]=1000;
circular[i]= kFALSE;
if (track->GetProlongation(0,y,z)) z0[i] = z;
dca[i] = track->GetD(0,0);
}
TStopwatch timer;
timer.Start();
Int_t ncandidates =0;
Int_t nall =0;
Int_t ntracks=0;
Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
for (Int_t i0=0;i0<nentries;i0++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
if (!track0) continue;
if (track0->GetNumberOfClusters()<40) continue;
if (TMath::Abs(1./track0->GetC())>200) continue;
for (Int_t i1=i0+1;i1<nentries;i1++){
AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
if (!track1) continue;
if (track1->GetNumberOfClusters()<40) continue;
if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
if (TMath::Abs(1./track1->GetC())>200) continue;
if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
if (track1->GetTgl()*track0->GetTgl()>0) continue;
if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue;
if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue;
Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
if (mindcar<5) continue;
Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
if (mindcaz<5) continue;
if (mindcar+mindcaz<20) continue;
Float_t xc0 = helixes[i0].GetHelix(6);
Float_t yc0 = helixes[i0].GetHelix(7);
Float_t r0 = helixes[i0].GetHelix(8);
Float_t xc1 = helixes[i1].GetHelix(6);
Float_t yc1 = helixes[i1].GetHelix(7);
Float_t r1 = helixes[i1].GetHelix(8);
Float_t rmean = (r0+r1)*0.5;
Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
if (delta>rmean*0.25) continue;
if (TMath::Abs(r0-r1)/rmean>0.3) continue;
Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
if (npoints==0) continue;
helixes[i0].GetClosestPhases(helixes[i1], phase);
Double_t xyz0[3];
Double_t xyz1[3];
Double_t hangles[3];
helixes[i0].Evaluate(phase[0][0],xyz0);
helixes[i1].Evaluate(phase[0][1],xyz1);
helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
Double_t deltah[2],deltabest;
if (hangles[2]<2.8) continue;
if (npoints>0){
Int_t ibest=0;
helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
if (npoints==2){
helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
if (deltah[1]<deltah[0]) ibest=1;
}
deltabest = TMath::Sqrt(deltah[ibest]);
helixes[i0].Evaluate(phase[ibest][0],xyz0);
helixes[i1].Evaluate(phase[ibest][1],xyz1);
helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
Double_t radiusbest = TMath::Sqrt(radius[ibest]);
if (deltabest>6) continue;
if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
Bool_t lsign =kFALSE;
if (hangles[2]>3.06) lsign =kTRUE;
if (lsign){
circular[i0] = kTRUE;
circular[i1] = kTRUE;
if (track0->OneOverPt()<track1->OneOverPt()){
track0->SetCircular(track0->GetCircular()+1);
track1->SetCircular(track1->GetCircular()+2);
}
else{
track1->SetCircular(track1->GetCircular()+1);
track0->SetCircular(track0->GetCircular()+2);
}
}
if (lsign&&((AliTPCReconstructor::StreamLevel()&kStreamFindCurling)>0)){
Int_t lab0=track0->GetLabel();
Int_t lab1=track1->GetLabel();
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Curling"<<
"lab0="<<lab0<<
"lab1="<<lab1<<
"Tr0.="<<track0<<
"Tr1.="<<track1<<
"dca0="<<dca[i0]<<
"dca1="<<dca[i1]<<
"mindcar="<<mindcar<<
"mindcaz="<<mindcaz<<
"delta="<<delta<<
"rmean="<<rmean<<
"npoints="<<npoints<<
"hangles0="<<hangles[0]<<
"hangles2="<<hangles[2]<<
"xyz0="<<xyz0[2]<<
"xyzz1="<<xyz1[2]<<
"z0="<<z0[i0]<<
"z1="<<z0[i1]<<
"radius="<<radiusbest<<
"deltabest="<<deltabest<<
"phase0="<<phase[ibest][0]<<
"phase1="<<phase[ibest][1]<<
"\n";
}
}
}
}
for (Int_t i =0;i<nentries;i++){
if (sign[i]==0) continue;
AliTPCseed * track0 = (AliTPCseed*)array->At(i);
if (track0==0) {
AliInfo("seed==0");
continue;
}
ntracks++;
Double_t cradius0 = 40*40;
Double_t cradius1 = 270*270;
Double_t cdist1=8.;
Double_t cdist2=8.;
Double_t cdist3=0.55;
for (Int_t j =i+1;j<nentries;j++){
nall++;
if (sign[j]*sign[i]<1) continue;
if ( (nclusters[i]+nclusters[j])>200) continue;
if ( (nclusters[i]+nclusters[j])<80) continue;
if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
if (npoints<1) continue;
if (npoints==1){
if (radius[0]<cradius0||radius[0]>cradius1) continue;
}
else{
if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
}
Double_t delta1=10000,delta2=10000;
helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
if (radius[0]<20&&delta1<1) continue;
if (radius[0]<10&&delta1<3) continue;
if (npoints==2){
helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
if (radius[1]<20&&delta2<1) continue;
if (radius[1]<10&&delta2<3) continue;
}
Double_t distance1 = TMath::Min(delta1,delta2);
if (distance1>cdist1) continue;
npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
if (radius[0]<20&&delta1<1) continue;
if (radius[0]<10&&delta1<3) continue;
if (npoints==2){
helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
if (radius[1]<20&&delta2<1) continue;
if (radius[1]<10&&delta2<3) continue;
}
distance1 = TMath::Min(delta1,delta2);
Float_t rkink =0;
if (delta1<delta2){
rkink = TMath::Sqrt(radius[0]);
}
else{
rkink = TMath::Sqrt(radius[1]);
}
if (distance1>cdist2) continue;
AliTPCseed * track1 = (AliTPCseed*)array->At(j);
Int_t row0 = GetRowNumber(rkink);
if (row0<10) continue;
if (row0>150) continue;
Float_t dens00=-1,dens01=-1;
Float_t dens10=-1,dens11=-1;
Int_t found,foundable,ishared;
track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
if (dens00<dens10 && dens01<dens11) continue;
if (dens00>dens10 && dens01>dens11) continue;
if (TMath::Max(dens00,dens10)<0.1) continue;
if (TMath::Max(dens01,dens11)<0.3) continue;
if (TMath::Min(dens00,dens10)>0.6) continue;
if (TMath::Min(dens01,dens11)>0.6) continue;
AliTPCseed * ktrack0, *ktrack1;
if (dens00>dens10){
ktrack0 = track0;
ktrack1 = track1;
}
else{
ktrack0 = track1;
ktrack1 = track0;
}
if (TMath::Abs(ktrack0->GetC())>5) continue;
AliExternalTrackParam paramm(*ktrack0);
AliExternalTrackParam paramd(*ktrack1);
if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
kink->SetMother(paramm);
kink->SetDaughter(paramd);
kink->Update();
Float_t x[3] = { static_cast<Float_t>(kink->GetPosition()[0]),static_cast<Float_t>(kink->GetPosition()[1]),static_cast<Float_t>(kink->GetPosition()[2])};
Int_t index[4];
fkParam->Transform0to1(x,index);
fkParam->Transform1to2(x,index);
row0 = GetRowNumber(x[0]);
if (kink->GetR()<100) continue;
if (kink->GetR()>240) continue;
if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue;
if (kink->GetDistance()>cdist3) continue;
Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1];
if (dird<0) continue;
Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1];
if (dirm<0) continue;
Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
if (mpt<0.2) continue;
if (mpt<1){
Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
if (qt>0.35) continue;
}
kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
if (dens00>dens10){
kink->SetTPCDensity(dens00,0,0);
kink->SetTPCDensity(dens01,0,1);
kink->SetTPCDensity(dens10,1,0);
kink->SetTPCDensity(dens11,1,1);
kink->SetIndex(i,0);
kink->SetIndex(j,1);
}
else{
kink->SetTPCDensity(dens10,0,0);
kink->SetTPCDensity(dens11,0,1);
kink->SetTPCDensity(dens00,1,0);
kink->SetTPCDensity(dens01,1,1);
kink->SetIndex(j,0);
kink->SetIndex(i,1);
}
if (mpt<1||kink->GetAngle(2)>0.1){
if (kink->GetTPCDensityFactor()<0.8) continue;
if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue;
if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
criticalangle= 3*TMath::Sqrt(criticalangle);
if (criticalangle>0.02) criticalangle=0.02;
if (kink->GetAngle(2)<criticalangle) continue;
}
Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
Float_t shapesum =0;
Float_t sum = 0;
for ( Int_t row = row0-drow; row<row0+drow;row++){
if (row<0) continue;
if (row>155) continue;
if (ktrack0->GetClusterPointer(row)){
AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
shapesum+=point->GetSigmaY()+point->GetSigmaZ();
sum++;
}
if (ktrack1->GetClusterPointer(row)){
AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
shapesum+=point->GetSigmaY()+point->GetSigmaZ();
sum++;
}
}
if (sum<4){
kink->SetShapeFactor(-1.);
}
else{
kink->SetShapeFactor(shapesum/sum);
}
Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
chi2P2*=chi2P2;
chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
chi2P3*=chi2P3;
chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
if ((AliTPCReconstructor::StreamLevel()&kStreamFindKinks)>0) {
(*fDebugStreamer)<<"kinkLpt"<<
"chi2P2="<<chi2P2<<
"chi2P3="<<chi2P3<<
"p0.="<<¶mm<<
"p1.="<<¶md<<
"k.="<<kink<<
"\n";
}
if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
continue;
}
kinks->AddLast(kink);
kink = new AliKink;
ncandidates++;
}
}
Int_t nkinks = kinks->GetEntriesFast();
Float_t *quality = new Float_t[nkinks];
Int_t *indexes = new Int_t[nkinks];
AliTPCseed *mothers = new AliTPCseed[nkinks];
AliTPCseed *daughters = new AliTPCseed[nkinks];
for (Int_t i=0;i<nkinks;i++){
quality[i] =100000;
AliKink *kinkl = (AliKink*)kinks->At(i);
Int_t index0 = kinkl->GetIndex(0);
Int_t index1 = kinkl->GetIndex(1);
AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
if (kinkl->GetAngle(2)<0.05){
kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
Int_t row0 = kinkl->GetTPCRow0();
Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
Int_t last = row0-drow;
if (last<40) last=40;
if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
Int_t first = row0+drow;
if (first>130) first=130;
if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
if (seed0 && seed1){
kinkl->SetStatus(1,8);
if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
row0 = GetRowNumber(kinkl->GetR());
sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
mothers[i] = *seed0;
daughters[i] = *seed1;
}
else{
delete kinks->RemoveAt(i);
if (seed0) MarkSeedFree( seed0 );
if (seed1) MarkSeedFree( seed1 );
continue;
}
if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
delete kinks->RemoveAt(i);
if (seed0) MarkSeedFree( seed0 );
if (seed1) MarkSeedFree( seed1 );
continue;
}
MarkSeedFree( seed0 );
MarkSeedFree( seed1 );
}
if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.);
}
TMath::Sort(nkinks,quality,indexes,kFALSE);
for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
if (!kink0) continue;
for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
kink0 = (AliKink*) kinks->At(indexes[ikink0]);
if (!kink0) continue;
AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
if (!kink1) continue;
if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
AliTPCseed &mother0 = mothers[indexes[ikink0]];
AliTPCseed &daughter0 = daughters[indexes[ikink0]];
AliTPCseed &mother1 = mothers[indexes[ikink1]];
AliTPCseed &daughter1 = daughters[indexes[ikink1]];
Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
Int_t same = 0;
Int_t both = 0;
Int_t samem = 0;
Int_t bothm = 0;
Int_t samed = 0;
Int_t bothd = 0;
for (Int_t i=0;i<row0;i++){
if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
both++;
bothm++;
if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
same++;
samem++;
}
}
}
for (Int_t i=row0;i<158;i++){
if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
both++;
bothd++;
if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
same++;
samed++;
}
}
}
Float_t ratio = Float_t(same+1)/Float_t(both+1);
Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
if (sum1>sum0){
shared[kink0->GetIndex(0)]= kTRUE;
shared[kink0->GetIndex(1)]= kTRUE;
delete kinks->RemoveAt(indexes[ikink0]);
break;
}
else{
shared[kink1->GetIndex(0)]= kTRUE;
shared[kink1->GetIndex(1)]= kTRUE;
delete kinks->RemoveAt(indexes[ikink1]);
}
}
}
}
for (Int_t i=0;i<nkinks;i++){
AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
if (!kinkl) continue;
kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
Int_t index0 = kinkl->GetIndex(0);
Int_t index1 = kinkl->GetIndex(1);
if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
kinkl->SetMultiple(usage[index0],0);
kinkl->SetMultiple(usage[index1],1);
if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
if (!ktrack0 || !ktrack1) continue;
Int_t index = esd->AddKink(kinkl);
if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) {
if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
*ktrack0 = mothers[indexes[i]];
*ktrack1 = daughters[indexes[i]];
}
}
ktrack0->SetKinkIndex(usage[index0],-(index+1));
ktrack1->SetKinkIndex(usage[index1], (index+1));
usage[index0]++;
usage[index1]++;
}
for (Int_t i=0;i<nentries;i++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i);
if (!track0) continue;
if (track0->GetKinkIndex(0)!=0) continue;
if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
}
RemoveUsed2(array,0.5,0.4,30);
UnsignClusters();
for (Int_t i=0;i<nentries;i++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i);
if (!track0) continue;
track0->CookdEdx(0.02,0.6);
track0->CookPID();
}
for (Int_t i=0;i<nentries;i++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i);
if (!track0) continue;
if (track0->Pt()<1.4) continue;
Int_t ishared=0;
Int_t all =0;
for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
if (track0->GetClusterPointer(icl)!=0){
all++;
if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
}
}
if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
MarkSeedFree( array->RemoveAt(i) );
continue;
}
if (track0->GetKinkIndex(0)!=0) continue;
if (track0->GetNumberOfClusters()<80) continue;
AliTPCseed *pmother = new AliTPCseed();
AliTPCseed *pdaughter = new AliTPCseed();
AliKink *pkink = new AliKink;
AliTPCseed & mother = *pmother;
AliTPCseed & daughter = *pdaughter;
AliKink & kinkl = *pkink;
if (CheckKinkPoint(track0,mother,daughter, kinkl)){
if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
delete pmother;
delete pdaughter;
delete pkink;
continue;
}
if (mother.Pt()<1.4) {
delete pmother;
delete pdaughter;
delete pkink;
continue;
}
Int_t row0= kinkl.GetTPCRow0();
if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
delete pmother;
delete pdaughter;
delete pkink;
continue;
}
Int_t index = esd->AddKink(&kinkl);
mother.SetKinkIndex(0,-(index+1));
daughter.SetKinkIndex(0,index+1);
if (mother.GetNumberOfClusters()>50) {
MarkSeedFree( array->RemoveAt(i) );
AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
mtc->SetPoolID(fLastSeedID);
array->AddAt(mtc,i);
}
else{
AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
mtc->SetPoolID(fLastSeedID);
array->AddLast(mtc);
}
AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
dtc->SetPoolID(fLastSeedID);
array->AddLast(dtc);
for (Int_t icl=0;icl<row0;icl++) {
if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
}
for (Int_t icl=row0;icl<158;icl++) {
if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
}
}
delete pmother;
delete pdaughter;
delete pkink;
}
delete [] daughters;
delete [] mothers;
delete [] dca;
delete []circular;
delete []shared;
delete []quality;
delete []indexes;
delete kink;
delete[]fim;
delete[] zm;
delete[] z0;
delete [] usage;
delete[] alpha;
delete[] nclusters;
delete[] sign;
delete[] helixes;
kinks->Delete();
delete kinks;
AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
timer.Print();
}
Int_t AliTPCtracker::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
{
AliKink &kink=(AliKink &)knk;
Int_t row0 = GetRowNumber(kink.GetR());
FollowProlongation(mother,0);
mother.Reset(kFALSE);
FollowProlongation(daughter,row0);
daughter.Reset(kFALSE);
FollowBackProlongation(daughter,158);
daughter.Reset(kFALSE);
Int_t first = TMath::Max(row0-20,30);
Int_t last = TMath::Min(row0+20,140);
const Int_t kNdiv =5;
AliTPCseed param0[kNdiv];
AliTPCseed param1[kNdiv];
AliKink kinks[kNdiv];
Int_t rows[kNdiv];
for (Int_t irow=0; irow<kNdiv;irow++){
rows[irow] = first +((last-first)*irow)/(kNdiv-1);
}
for (Int_t irow=0;irow<kNdiv;irow++){
FollowBackProlongation(mother, rows[irow]);
FollowProlongation(daughter,rows[kNdiv-1-irow]);
param0[irow] = mother;
param1[kNdiv-1-irow] = daughter;
}
for (Int_t irow=0; irow<kNdiv-1;irow++){
if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
kinks[irow].SetMother(param0[irow]);
kinks[irow].SetDaughter(param1[irow]);
kinks[irow].Update();
}
Int_t index =-1;
Double_t mindist = 10000;
for (Int_t irow=0;irow<kNdiv;irow++){
if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
if (TMath::Abs(kinks[irow].GetR())>240.) continue;
if (TMath::Abs(kinks[irow].GetR())<100.) continue;
Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
if (normdist < mindist){
mindist = normdist;
index = irow;
}
}
if (index==-1) return 0;
param0[index].Reset(kTRUE);
FollowProlongation(param0[index],0);
mother = param0[index];
daughter = param1[index];
kink.SetMother(mother);
kink.SetDaughter(daughter);
kink.Update();
kink.SetTPCRow0(GetRowNumber(kink.GetR()));
kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
mother.SetLabel(kink.GetLabel(0));
daughter.SetLabel(kink.GetLabel(1));
return 1;
}
void AliTPCtracker::UpdateKinkQualityM(AliTPCseed * seed){
if (seed->GetKinkIndex(0)>=0) return;
for (Int_t ikink=0;ikink<3;ikink++){
Int_t index = seed->GetKinkIndex(ikink);
if (index>=0) break;
index = TMath::Abs(index)-1;
AliESDkink * kink = fEvent->GetKink(index);
kink->SetTPCDensity(-1,0,0);
kink->SetTPCDensity(1,0,1);
Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
if (row0<15) row0=15;
Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
if (row1>145) row1=145;
Int_t found,foundable,shared;
seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
}
}
void AliTPCtracker::UpdateKinkQualityD(AliTPCseed * seed){
if (seed->GetKinkIndex(0)<=0) return;
for (Int_t ikink=0;ikink<3;ikink++){
Int_t index = seed->GetKinkIndex(ikink);
if (index<=0) break;
index = TMath::Abs(index)-1;
AliESDkink * kink = fEvent->GetKink(index);
kink->SetTPCDensity(-1,1,0);
kink->SetTPCDensity(-1,1,1);
Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
if (row0<15) row0=15;
Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
if (row1>145) row1=145;
Int_t found,foundable,shared;
seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
}
}
Int_t AliTPCtracker::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
{
AliKink &kink=(AliKink &)knk;
Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
Int_t first = seed->GetFirstPoint();
Int_t last = seed->GetLastPoint();
if (last-first<20) return 0;
AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE);
if (!seed1) return 0;
FollowProlongation(*seed1,seed->GetLastPoint()-20);
seed1->Reset(kTRUE);
FollowProlongation(*seed1,158);
seed1->Reset(kTRUE);
last = seed1->GetLastPoint();
AliTPCseed *seed0 = new( NextFreeSeed() ) AliTPCseed(*seed);
seed0->SetPoolID(fLastSeedID);
seed0->Reset(kFALSE);
seed0->Reset();
AliTPCseed param0[20];
AliTPCseed param1[20];
AliKink kinks[20];
Int_t rows[20];
for (Int_t irow=0; irow<20;irow++){
rows[irow] = first +((last-first)*irow)/19;
}
for (Int_t irow=0;irow<20;irow++){
FollowBackProlongation(*seed0, rows[irow]);
FollowProlongation(*seed1,rows[19-irow]);
param0[irow] = *seed0;
param1[19-irow] = *seed1;
}
for (Int_t irow=0; irow<19;irow++){
kinks[irow].SetMother(param0[irow]);
kinks[irow].SetDaughter(param1[irow]);
kinks[irow].Update();
}
Int_t index =-1;
Double_t maxchange= 0;
for (Int_t irow=1;irow<19;irow++){
if (TMath::Abs(kinks[irow].GetR())>240.) continue;
if (TMath::Abs(kinks[irow].GetR())<110.) continue;
Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
if ( quality > maxchange){
maxchange = quality;
index = irow;
}
}
MarkSeedFree( seed0 );
MarkSeedFree( seed1 );
if (index<0) return 0;
Int_t row0 = GetRowNumber(kinks[index].GetR());
seed0 = new( NextFreeSeed() ) AliTPCseed(param0[index]);
seed0->SetPoolID(fLastSeedID);
seed1 = new( NextFreeSeed() ) AliTPCseed(param1[index]);
seed1->SetPoolID(fLastSeedID);
seed0->Reset(kFALSE);
seed1->Reset(kFALSE);
seed0->ResetCovariance(10.);
seed1->ResetCovariance(10.);
FollowProlongation(*seed0,0);
FollowBackProlongation(*seed1,158);
mother = *seed0;
seed0->Reset(kFALSE);
seed1->Reset(kFALSE);
seed0->ResetCovariance(10.);
seed1->ResetCovariance(10.);
first = TMath::Max(row0-20,0);
last = TMath::Min(row0+20,158);
for (Int_t irow=0; irow<20;irow++){
rows[irow] = first +((last-first)*irow)/19;
}
for (Int_t irow=0;irow<20;irow++){
FollowBackProlongation(*seed0, rows[irow]);
FollowProlongation(*seed1,rows[19-irow]);
param0[irow] = *seed0;
param1[19-irow] = *seed1;
}
for (Int_t irow=0; irow<19;irow++){
kinks[irow].SetMother(param0[irow]);
kinks[irow].SetDaughter(param1[irow]);
kinks[irow].Update();
}
index =-1;
maxchange= 0;
for (Int_t irow=0;irow<20;irow++){
if (TMath::Abs(kinks[irow].GetR())>250.) continue;
if (TMath::Abs(kinks[irow].GetR())<90.) continue;
Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
if ( quality > maxchange){
maxchange = quality;
index = irow;
}
}
if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
MarkSeedFree( seed0 );
MarkSeedFree( seed1 );
return 0;
}
kink.SetMother(param0[index]);
kink.SetDaughter(param1[index]);
kink.Update();
Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
chi2P2*=chi2P2;
chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
chi2P3*=chi2P3;
chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
if (AliTPCReconstructor::StreamLevel()&kStreamFindKinks) {
(*fDebugStreamer)<<"kinkHpt"<<
"chi2P2="<<chi2P2<<
"chi2P3="<<chi2P3<<
"p0.="<<¶m0[index]<<
"p1.="<<¶m1[index]<<
"k.="<<&kink<<
"\n";
}
if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
MarkSeedFree( seed0 );
MarkSeedFree( seed1 );
return 0;
}
row0 = GetRowNumber(kink.GetR());
kink.SetTPCRow0(row0);
kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
kink.SetIndex(-10,0);
kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
daughter = param1[index];
daughter.SetLabel(kink.GetLabel(1));
param0[index].Reset(kTRUE);
FollowProlongation(param0[index],0);
mother = param0[index];
mother.SetLabel(kink.GetLabel(0));
if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
mother=*seed;
}
MarkSeedFree( seed0 );
MarkSeedFree( seed1 );
return 1;
}
AliTPCseed* AliTPCtracker::ReSeed(AliTPCseed *t)
{
Int_t first = 0;
if (fSectors == fOuterSec){
first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
}
else
first = t->GetFirstPoint();
AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
FollowBackProlongation(*t,fSectors->GetNRows()-1);
t->Reset(kFALSE);
FollowProlongation(*t,first);
return seed;
}
Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
TDirectory *savedir=gDirectory;
TFile *in=(TFile*)inp;
if (!in->IsOpen()) {
cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
return 1;
}
in->cd();
TTree *seedTree=(TTree*)in->Get("Seeds");
if (!seedTree) {
cerr<<"AliTPCtracker::ReadSeeds(): ";
cerr<<"can't get a tree with track seeds !\n";
return 2;
}
AliTPCtrack *seed=new AliTPCtrack;
seedTree->SetBranchAddress("tracks",&seed);
if (fSeeds==0) fSeeds=new TObjArray(15000);
Int_t n=(Int_t)seedTree->GetEntries();
for (Int_t i=0; i<n; i++) {
seedTree->GetEvent(i);
AliTPCseed* sdc = new( NextFreeSeed() ) AliTPCseed(*seed);
sdc->SetPoolID(fLastSeedID);
fSeeds->AddLast(sdc);
}
delete seed;
delete seedTree;
savedir->cd();
return 0;
}
Int_t AliTPCtracker::Clusters2TracksHLT (AliESDEvent *const esd, const AliESDEvent *hltEvent)
{
if (fSeeds) DeleteSeeds();
else ResetSeedsPool();
fEvent = esd;
fEventHLT = hltEvent;
if (AliTPCReconstructor::GetRecoParam()->GetUseOulierClusterFilter()) FilterOutlierClusters();
AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
transform->SetCurrentTimeStamp( esd->GetTimeStamp());
transform->SetCurrentRun(esd->GetRunNumber());
Clusters2Tracks();
fEventHLT = 0;
if (!fSeeds) return 1;
FillESD(fSeeds);
if ((AliTPCReconstructor::StreamLevel()&kStreamClDump)>0) DumpClusters(0,fSeeds);
return 0;
}
Int_t AliTPCtracker::Clusters2Tracks(AliESDEvent *const esd)
{
return Clusters2TracksHLT( esd, 0);
}
Int_t AliTPCtracker::Clusters2Tracks() {
TDirectory *savedir=gDirectory;
TStopwatch timer;
fIteration = 0;
fSeeds = Tracking();
if (fDebug>0){
Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
}
for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
if (!pt) continue;
Int_t nc=t.GetNumberOfClusters();
if (nc<20) {
MarkSeedFree( fSeeds->RemoveAt(i) );
continue;
}
CookLabel(pt,0.1);
if (pt->GetRemoval()==10) {
if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
pt->Desactivate(10);
else{
pt->Desactivate(20);
MarkSeedFree( fSeeds->RemoveAt(i) );
}
}
}
RemoveUsed2(fSeeds,0.85,0.85,0);
if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
if (AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC) FindMultiMC(fSeeds, fEvent,-1);
RemoveUsed2(fSeeds,0.5,0.4,20);
FindSplitted(fSeeds, fEvent,0);
if (AliTPCReconstructor::StreamLevel()&kStreamFindMultiMC) FindMultiMC(fSeeds, fEvent,0);
Int_t nseed=fSeeds->GetEntriesFast();
Int_t found = 0;
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
if (!pt) continue;
Int_t nc=t.GetNumberOfClusters();
if (nc<15) {
MarkSeedFree( fSeeds->RemoveAt(i) );
continue;
}
CookLabel(pt,0.1);
if ((pt->IsActive() || (pt->GetRemoval()==10) )){
found++;
if (fDebug>0) cerr<<found<<'\r';
pt->SetLab2(i);
}
else
MarkSeedFree( fSeeds->RemoveAt(i) );
}
SignShared(fSeeds);
nseed=fSeeds->GetEntriesFast();
found = 0;
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
if (!pt) continue;
Int_t nc=t.GetNumberOfClusters();
if (nc<15) {
MarkSeedFree( fSeeds->RemoveAt(i) );
continue;
}
t.SetUniqueID(i);
t.CookdEdx(0.02,0.6);
if ((pt->IsActive() || (pt->GetRemoval()==10) )){
found++;
if (fDebug>0){
cerr<<found<<'\r';
}
pt->SetLab2(i);
}
else
MarkSeedFree( fSeeds->RemoveAt(i) );
}
SortTracks(fSeeds, 1);
if (fDebug>0){
Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
}
Info("Clusters2Tracks","Number of found tracks %d",found);
savedir->cd();
return 0;
}
void AliTPCtracker::Tracking(TObjArray * arr)
{
fSectors = fOuterSec;
ParallelTracking(arr,150,63);
fSectors = fOuterSec;
ParallelTracking(arr,63,0);
}
TObjArray * AliTPCtracker::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
{
static TObjArray arrTracks;
TObjArray * arr = &arrTracks;
fSectors = fOuterSec;
TStopwatch timer;
timer.Start();
for (Int_t sec=0;sec<fkNOS;sec++){
if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
}
if (fDebug>0){
Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
timer.Print();
timer.Start();
}
Tracking(arr);
if (fDebug>0){
timer.Print();
}
return arr;
}
TObjArray * AliTPCtracker::Tracking()
{
if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
TStopwatch timer;
timer.Start();
Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
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];
cuts[0] = 0.002;
cuts[1] = 1.5;
cuts[2] = 3.;
cuts[3] = 3.;
Float_t fnumber = 3.0;
Float_t fdensity = 3.0;
if (AliTPCReconstructor::GetRecoParam()->GetUseHLTPreSeeding()) {
arr = MakeSeedsHLT( fEventHLT );
if( arr ){
SumTracks(seeds,arr);
delete arr;
arr=0;
}
}
cuts[0]=0.0066;
for (Int_t delta = 0; delta<18; delta+=gapPrim){
cuts[0]=0.0070;
cuts[1] = 1.5;
arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
for (Int_t i=2;i<6;i+=2){
cuts[0]=0.0022;
cuts[1]=0.3;
arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
}
}
fnumber = 4;
fdensity = 4.;
cuts[0]=0.0077;
for (Int_t delta = 20; delta<120; delta+=gapPrim){
cuts[0]=0.0060;
cuts[1]=0.3;
cuts[2]=6.;
arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
cuts[0]=0.003;
cuts[1]=0.3;
cuts[2]=6.;
arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
}
cuts[0] = 0.01;
cuts[1] = 2.0;
cuts[2] = 3.;
cuts[3] = 2.0;
fnumber = 2.;
fdensity = 2.;
if (fDebug>0){
Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
timer.Print();
timer.Start();
}
cuts[0] = 0.3;
cuts[1] = 1.5;
cuts[2] = 3.;
cuts[3] = 1.5;
arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
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 = 9; delta<30; delta+=gapSec){
cuts[0] = 0.3;
cuts[1] = 1.5;
cuts[2] = 3.;
cuts[3] = 1.5;
arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
}
fnumber = 1;
fdensity = 1;
fnumber = 2.;
fdensity = 2.;
cuts[0]=0.0080;
for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
cuts[0] = 0.3;
cuts[1] = 3.5;
cuts[2] = 3.;
cuts[3] = 3.5;
arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
}
if (fDebug>0){
Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
timer.Print();
timer.Start();
}
return seeds;
}
TObjArray * AliTPCtracker::TrackingSpecial()
{
TStopwatch timer;
timer.Start();
Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
TObjArray * seeds = new TObjArray;
TObjArray * arr=0;
Int_t gap = 15;
Float_t cuts[4];
Float_t fnumber = 3.0;
Float_t fdensity = 3.0;
cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC();
cuts[1] = 3.5;
cuts[2] = 3.;
cuts[3] = 3.5;
for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
}
if (fDebug>0){
Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
timer.Print();
timer.Start();
}
return seeds;
}
void AliTPCtracker::SumTracks(TObjArray *arr1,TObjArray *&arr2)
{
Int_t nseed = arr2->GetEntriesFast();
for (Int_t i=0;i<nseed;i++){
AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
if (pt){
if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
MarkSeedFree( arr2->RemoveAt(i) );
continue;
}
if (pt->GetNumberOfClusters()<20){
MarkSeedFree( arr2->RemoveAt(i) );
continue;
}
if (pt->IsActive()){
arr1->AddLast(arr2->RemoveAt(i));
continue;
}
if (pt->GetRemoval()!=10){
MarkSeedFree( arr2->RemoveAt(i) );
continue;
}
if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
arr1->AddLast(arr2->RemoveAt(i));
else{
MarkSeedFree( arr2->RemoveAt(i) );
}
}
}
}
void AliTPCtracker::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
{
Int_t nseed=arr->GetEntriesFast();
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
if (!pt) continue;
if (!t.IsActive()) continue;
if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
FollowProlongation(t, rfirst+1);
}
for (Int_t nr=rfirst; nr>=rlast; nr--){
if (nr<fInnerSec->GetNRows())
fSectors = fInnerSec;
else
fSectors = fOuterSec;
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
if (!pt) continue;
if (nr==80) pt->UpdateReference();
if (!pt->IsActive()) continue;
if (pt->GetRelativeSector()>17) {
continue;
}
UpdateClusters(t,nr);
}
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) continue;
if (!pt->IsActive()) continue;
if (pt->GetRelativeSector()>17) {
continue;
}
FollowToNextCluster(*pt,nr);
}
}
}
void AliTPCtracker::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
{
Int_t nseed= arr->GetEntriesFast();
for (Int_t i=0;i<nseed;i++){
AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
if (pt) {
pt->Modify(fac);
Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
Int_t sec = (index&0xff000000)>>24;
sec = sec%18;
Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
if (angle1>TMath::Pi())
angle1-=2.*TMath::Pi();
Float_t angle2 = pt->GetAlpha();
if (TMath::Abs(angle1-angle2)>0.001){
if (!pt->Rotate(angle1-angle2)) return;
}
}
}
}
void AliTPCtracker::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
{
Int_t nseed= arr->GetEntriesFast();
for (Int_t i=0;i<nseed;i++){
AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
if (pt) {
pt->Modify(fac);
pt->SetFirstPoint(pt->GetLastPoint());
}
}
}
Int_t AliTPCtracker::PropagateBack(const TObjArray *const arr)
{
Int_t nseed= arr->GetEntriesFast();
for (Int_t i=0;i<nseed;i++){
AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
if (pt&& pt->GetKinkIndex(0)<=0) {
fSectors = fInnerSec;
FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
}
if (pt&& pt->GetKinkIndex(0)>0) {
AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
pt->SetFirstPoint(kink->GetTPCRow0());
fSectors = fInnerSec;
FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
}
CookLabel(pt,0.3);
}
return 0;
}
Int_t AliTPCtracker::PropagateForward2(const TObjArray *const arr)
{
Int_t nseed= arr->GetEntriesFast();
for (Int_t i=0;i<nseed;i++){
AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
if (pt) {
FollowProlongation(*pt,0,1,1);
CookLabel(pt,0.3);
}
}
return 0;
}
Int_t AliTPCtracker::PropagateForward()
{
Int_t nseed = fSeeds->GetEntriesFast();
for (Int_t i=0;i<nseed;i++){
AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
if (pt){
AliTPCseed &t = *pt;
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);
}
}
fSectors = fOuterSec;
ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
fSectors = fInnerSec;
ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
return 1;
}
Int_t AliTPCtracker::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
{
if (pt) {
fSectors = fInnerSec;
Int_t r1;
if (row1<fSectors->GetNRows())
r1 = row1;
else
r1 = fSectors->GetNRows()-1;
if (row0<fSectors->GetNRows()&& r1>0 )
FollowBackProlongation(*pt,r1);
if (row1<=fSectors->GetNRows())
return 0;
r1 = row1 - fSectors->GetNRows();
if (r1<=0) return 0;
if (r1>=fOuterSec->GetNRows()) return 0;
fSectors = fOuterSec;
return FollowBackProlongation(*pt,r1);
}
return 0;
}
void AliTPCtracker::GetShape(AliTPCseed * seed, Int_t row)
{
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();
if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
}
angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
seed->SetCurrentSigmaY2(sigmay*sigmay);
seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
}
void AliTPCtracker::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
if(!t){
printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
return;
}
Int_t noc=t->GetNumberOfClusters();
if (noc<10){
return ;
}
Int_t lb[160];
Int_t mx[160];
AliTPCclusterMI *clusters[160];
for (Int_t i=0;i<160;i++) {
clusters[i]=0;
lb[i]=mx[i]=0;
}
Int_t i;
Int_t current=0;
for (i=0; i<160 && current<noc; i++) {
Int_t index=t->GetClusterIndex2(i);
if (index<=0) continue;
if (index&0x8000) continue;
if (t->GetClusterPointer(i)){
clusters[current]=t->GetClusterPointer(i);
current++;
}
}
noc = current;
Int_t lab=123456789;
for (i=0; i<noc; i++) {
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
lab=TMath::Abs(c->GetLabel(0));
Int_t j;
for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
lb[j]=lab;
(mx[j])++;
}
Int_t max=0;
for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
for (i=0; i<noc; i++) {
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
if (TMath::Abs(c->GetLabel(1)) == lab ||
TMath::Abs(c->GetLabel(2)) == lab ) max++;
}
if (noc<=0) { lab=-1; return;}
if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
else {
Int_t tail=Int_t(0.10*noc);
max=0;
Int_t ind=0;
for (i=1; i<160&&ind<tail; i++) {
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
if (lab == TMath::Abs(c->GetLabel(0)) ||
lab == TMath::Abs(c->GetLabel(1)) ||
lab == TMath::Abs(c->GetLabel(2))) max++;
ind++;
}
if (max < Int_t(0.5*tail)) lab=-lab;
}
t->SetLabel(lab);
}
Int_t AliTPCtracker::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
Int_t noc=t->GetNumberOfClusters();
if (noc<10){
return -1;
}
Int_t lb[160];
Int_t mx[160];
AliTPCclusterMI *clusters[160];
for (Int_t i=0;i<160;i++) {
clusters[i]=0;
lb[i]=mx[i]=0;
}
Int_t i;
Int_t current=0;
for (i=0; i<160 && current<noc; i++) {
if (i<first) continue;
if (i>last) continue;
Int_t index=t->GetClusterIndex2(i);
if (index<=0) continue;
if (index&0x8000) continue;
if (t->GetClusterPointer(i)){
clusters[current]=t->GetClusterPointer(i);
current++;
}
}
noc = current;
Int_t lab=123456789;
for (i=0; i<noc; i++) {
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
lab=TMath::Abs(c->GetLabel(0));
Int_t j;
for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
lb[j]=lab;
(mx[j])++;
}
Int_t max=0;
for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
for (i=0; i<noc; i++) {
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
if (TMath::Abs(c->GetLabel(1)) == lab ||
TMath::Abs(c->GetLabel(2)) == lab ) max++;
}
if (noc<=0) { lab=-1; return -1;}
if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
else {
Int_t tail=Int_t(0.10*noc);
max=0;
Int_t ind=0;
for (i=1; i<160&&ind<tail; i++) {
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
if (lab == TMath::Abs(c->GetLabel(0)) ||
lab == TMath::Abs(c->GetLabel(1)) ||
lab == TMath::Abs(c->GetLabel(2))) max++;
ind++;
}
if (max < Int_t(0.5*tail)) lab=-lab;
}
return lab;
}
Int_t AliTPCtracker::GetRowNumber(Double_t x[3]) const
{
Float_t phi = TMath::ATan2(x[1],x[0]);
if(phi<0) phi=2.*TMath::Pi()+phi;
const Float_t kRaddeg = 180/3.14159265358979312;
Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
return GetRowNumber(localx);
}
void AliTPCtracker::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
{
Int_t firstpoint = 0;
Int_t lastpoint = 159;
AliTPCTrackerPoint *point;
AliTPCclusterMI *cluster;
Int_t nclsf = 0;
TBits clusterMap(159);
TBits sharedMap(159);
TBits fitMap(159);
for (int iter=firstpoint; iter<lastpoint; iter++) {
cluster = t->GetClusterPointer(iter);
if (cluster) {
clusterMap.SetBitNumber(iter, kTRUE);
point = t->GetTrackPoint(iter);
if (point->IsShared())
sharedMap.SetBitNumber(iter,kTRUE);
}
if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
fitMap.SetBitNumber(iter, kTRUE);
nclsf++;
}
}
esd->SetTPCClusterMap(clusterMap);
esd->SetTPCSharedMap(sharedMap);
esd->SetTPCFitMap(fitMap);
if (nclsf != t->GetNumberOfClusters())
AliDebug(3,Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
}
Bool_t AliTPCtracker::IsFindable(AliTPCseed & track){
Float_t kDeltaZ=10;
Float_t z = track.GetZ();
if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
TMath::Abs(z)<fkParam->GetZLength(0) &&
(TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
return kTRUE;
return kFALSE;
}
void AliTPCtracker::AddCovariance(AliTPCseed * seed){
const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
Double_t *covarS= (Double_t*)seed->GetCovariance();
Double_t factor[5]={1,1,1,1,1};
factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] +param[4]*param[4])/covarS[14]));
factor[0]=factor[2];
factor[4]=factor[2];
for (Int_t i=0; i<5; i++){
for (Int_t j=i; j<5; j++){
Int_t index=seed->GetIndex(i,j);
covarS[index]*=factor[i]*factor[j];
}
}
}
void AliTPCtracker::AddCovarianceAdd(AliTPCseed * seed){
const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
Double_t *covarIn= (Double_t*)seed->GetCovariance();
Double_t covar[15];
for (Int_t i=0;i<15;i++) covar[i]=0;
covar[0] = param[0]*param[0];
covar[2] = param[1]*param[1];
covar[5] = param[2]*param[2];
covar[9] = param[3]*param[3];
covar[14]= param[4]*param[4];
covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
seed->AddCovariance(covar);
}
Bool_t AliTPCtracker::IsTPCHVDipEvent(AliESDEvent const *esdEvent)
{
if(!esdEvent) return kFALSE;
AliTPCcalibDB *db=AliTPCcalibDB::Instance();
if(!db) return kFALSE;
db->SetRun(esdEvent->GetRunNumber());
const Double_t kTPCHVdip = db->GetParameters()->GetMaxDipVoltage();
const Double_t dipEventScanPeriod = db->GetParameters()->GetVoltageDipScanPeriod();
const Double_t tevSec = esdEvent->GetTimeStamp();
for(Int_t sector=0; sector<72; sector++)
{
if (!db->GetChamberHVStatus(sector)) continue;
AliDCSSensor *sensor = db->GetChamberHVSensor(sector);
if (!sensor) continue;
TGraph *grSensor=sensor->GetGraph();
if (!grSensor) continue;
if (grSensor->GetN()<1) continue;
const Double_t median = db->GetChamberHighVoltageMedian(sector);
if(median < 1.) continue;
for (Int_t ipoint=0; ipoint<grSensor->GetN()-1; ++ipoint){
Double_t nextTime=grSensor->GetX()[ipoint+1]*3600+sensor->GetStartTime();
if (tevSec-dipEventScanPeriod>nextTime) continue;
const Float_t deltaV=TMath::Abs(grSensor->GetY()[ipoint]-median);
if (deltaV>kTPCHVdip) {
AliDebug(3,Form("HV dip detected in ROC '%02d' with dV '%.2f' at time stamp '%.0f'",sector,deltaV,tevSec));
return kTRUE;
}
if (nextTime>tevSec+dipEventScanPeriod) break;
}
}
return kFALSE;
}
void AliTPCtracker::MarkSeedFree(TObject *sd)
{
AliTPCseed* seed = dynamic_cast<AliTPCseed*>(sd);
if (!seed) {
AliError(Form("Freeing of non-AliTPCseed %p from the pool is requested",sd));
return;
}
int id = seed->GetPoolID();
if (id<0) {
AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
return;
}
fSeedsPool->RemoveAt(id);
if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
}
TObject *&AliTPCtracker::NextFreeSeed()
{
fLastSeedID = fNFreeSeeds ? fFreeSeedsID.GetArray()[--fNFreeSeeds] : fSeedsPool->GetEntriesFast();
return (*fSeedsPool)[ fLastSeedID ];
}
void AliTPCtracker::ResetSeedsPool()
{
AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool->GetSize(),fSeedsPool->GetEntriesFast(),fNFreeSeeds));
fNFreeSeeds = 0;
fSeedsPool->Clear("C");
}
Int_t AliTPCtracker::PropagateToRowHLT(AliTPCseed *pt, int nrow)
{
AliTPCseed &t=*pt;
Double_t x= GetXrow(nrow);
Double_t ymax= GetMaxY(nrow);
Int_t rotate = 0;
Int_t nRotations=0;
int ret = 1;
do{
rotate = 0;
if (!t.PropagateTo(x) ){
ret = 0;
break;
}
t.SetRow(nrow);
Double_t y = t.GetY();
if( y > ymax ) {
if( rotate!=-1 ) rotate=1;
} else if (y <-ymax) {
if( rotate!=1 ) rotate = -1;
}
if( rotate==0 ) break;
if (!t.Rotate( rotate==1 ?fSectors->GetAlpha() :-fSectors->GetAlpha())) {
ret=0;
break;
}
nRotations+= rotate;
}while(rotate!=0);
if( nRotations!=0 ){
int newSec= t.GetRelativeSector()+nRotations;
if( newSec>=fN ) newSec-=fN;
else if( newSec<0 ) newSec +=fN;
t.SetRelativeSector(newSec);
}
return ret;
}
void AliTPCtracker::TrackFollowingHLT(TObjArray *const arr )
{
Int_t nRows=fOuterSec->GetNRows()+fInnerSec->GetNRows();
fSectors=fInnerSec;
Int_t nseed=arr->GetEntriesFast();
double shapeY2[160], shapeZ2[160];
Int_t clusterIndex[160];
for (Int_t iSeed=0; iSeed<nseed; iSeed++) {
AliTPCseed *pt=(AliTPCseed*) (arr->UncheckedAt(iSeed));
if (!pt) continue;
AliTPCseed &t=*pt;
for( int iter=0; iter<3; iter++ ){
t.Reset();
t.SetLastPoint(0);
t.SetFirstPoint(nRows-1);
t.ResetCovariance(.1);
t.SetNumberOfClusters(0);
for( int i=0; i<nRows; i++ ){
shapeY2[i]=1.;
shapeZ2[i]=1.;
clusterIndex[i]=-1;
t.SetClusterIndex2(i,-1);
t.SetClusterIndex(i,-1);
}
Double_t roady = 20.;
Double_t roadz = 20.;
double roadr = 5;
AliTPCseed t0(t);
t0.Reset();
int nClusters = 0;
{
t0.SetRelativeSector(t.GetRelativeSector());
t0.SetLastPoint(0);
t0.SetFirstPoint(159);
for (Int_t nr=0; nr<nRows; nr++){
if( nr<fInnerSec->GetNRows() ) fSectors=fInnerSec;
else fSectors=fOuterSec;
if( !PropagateToRowHLT(&t0, nr ) ){ break; }
if (TMath::Abs(t0.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()){
continue;
}
if (!IsActive(t0.GetRelativeSector(),nr)) {
continue;
}
if( iter==0 ){
GetShape(&t0,nr);
shapeY2[nr]=t0.GetCurrentSigmaY2();
shapeZ2[nr]=t0.GetCurrentSigmaZ2();
}
AliTPCtrackerRow &krow=GetRow(t0.GetRelativeSector(),nr);
if( !krow ) continue;
t.SetClusterIndex2(nr,-3);
t.SetClusterIndex(nr,-3);
AliTPCclusterMI *cl=0;
UInt_t uindex = 0;
cl = krow.FindNearest2(t0.GetY(),t0.GetZ(),roady,roadz,uindex);
if (!cl ) continue;
double dy = cl->GetY()-t0.GetY();
double dz = cl->GetZ()-t0.GetZ();
double dr = sqrt(dy*dy+dz*dz);
if( dr>roadr ){
continue;
}
t0.SetClusterPointer(nr, cl);
clusterIndex[nr] = krow.GetIndex(uindex);
if( t0.GetFirstPoint()>nr ) t0.SetFirstPoint(nr);
t0.SetLastPoint(nr);
nClusters++;
}
}
if( nClusters <3 ){
break;
}
Int_t basePoints[3] = {t0.GetFirstPoint(),t0.GetFirstPoint(),t0.GetLastPoint()};
{
Int_t midRow = (t0.GetLastPoint()-t0.GetFirstPoint())/2;
int dist=200;
for( int nr=t0.GetFirstPoint()+1; nr< t0.GetLastPoint(); nr++){
if( !t0.GetClusterPointer(nr) ) continue;
int d = TMath::Abs(nr-midRow);
if( d < dist ){
dist = d;
basePoints[1] = nr;
}
}
}
if( 1||iter<2 ){
for( int icl=0; icl<3; icl++){
int nr = basePoints[icl];
int lr=nr;
if( nr>=fInnerSec->GetNRows()){
lr = nr - fInnerSec->GetNRows();
fSectors=fOuterSec;
} else fSectors=fInnerSec;
AliTPCclusterMI *cl=t0.GetClusterPointer(nr);
if(!cl){
continue;
}
int iSec = cl->GetDetector() %fkNIS;
int rotate = iSec - t.GetRelativeSector();
if( rotate!=0 ){
if (!t.Rotate( rotate*fSectors->GetAlpha()) ) {
break;
}
t.SetRelativeSector(iSec);
}
Double_t x= cl->GetX();
if (!t.PropagateTo(x)){
break;
}
if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()){
break;
}
t.SetCurrentClusterIndex1(clusterIndex[nr]);
t.SetCurrentCluster(cl);
t.SetRow(lr);
t.SetErrorY2(shapeY2[nr]);
t.SetErrorZ2(shapeZ2[nr]);
if( icl==0 ){
double cov[15];
for( int j=0; j<15; j++ ) cov[j]=0;
cov[0]=10;
cov[2]=10;
cov[5]=.5;
cov[9]=.5;
cov[14]=1.;
t.AliExternalTrackParam::AddCovariance(cov);
}
if( !UpdateTrack(&t,0) ){
t.SetClusterIndex2(nr,-1);
t.SetClusterIndex(nr,-1);
t.SetClusterPointer(nr,0);
break;
}
}
for (Int_t nr=t0.GetLastPoint(); nr>=t0.GetFirstPoint(); nr-- ){
int lr=nr;
if( nr>=fInnerSec->GetNRows()){
lr = nr - fInnerSec->GetNRows();
fSectors=fOuterSec;
} else fSectors=fInnerSec;
if(1|| iter<2 ){
if( nr == basePoints[0] ) continue;
if( nr == basePoints[1] ) continue;
if( nr == basePoints[2] ) continue;
}
AliTPCclusterMI *cl=t0.GetClusterPointer(nr);
if(!cl) continue;
int iSec = cl->GetDetector() %fkNIS;
int rotate = iSec - t.GetRelativeSector();
if( rotate!=0 ){
if (!t.Rotate( rotate*fSectors->GetAlpha()) ) {
break;
}
t.SetRelativeSector(iSec);
}
Double_t x= cl->GetX();
if (!t.PropagateTo(x)){
break;
}
if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()){
break;
}
t.SetCurrentClusterIndex1(clusterIndex[nr]);
t.SetCurrentCluster(cl);
t.SetRow(lr);
t.SetErrorY2(shapeY2[nr]);
t.SetErrorZ2(shapeZ2[nr]);
if( !UpdateTrack(&t,0) ){
t.SetClusterIndex2(nr,-1);
t.SetClusterIndex(nr,-1);
break;
}
}
}
}
Int_t foundable,found,shared;
t.GetClusterStatistic(0,nRows, found, foundable, shared, kTRUE);
t.SetNFoundable(foundable);
}
}
TObjArray * AliTPCtracker::MakeSeedsHLT(const AliESDEvent *hltEvent)
{
if( !hltEvent ) return 0;
Int_t nentr=hltEvent->GetNumberOfTracks();
AliInfo(Form("Using %d HLT tracks for seeding",nentr));
TObjArray * seeds = new TObjArray(nentr);
Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
Int_t index = 0;
Int_t nTr=hltEvent->GetNumberOfTracks();
for( int itr=0; itr<nTr; itr++ ){
const AliExternalTrackParam *param = hltEvent->GetTrack(itr)->GetTPCInnerParam();
if( !param ) continue;
AliTPCtrack tr;
tr.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
tr.SetNumberOfClusters(0);
AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed(tr);
Double_t alpha=seed->GetAlpha();
if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
if (alpha < 0. ) alpha += 2.*TMath::Pi();
seed->SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
Double_t alphaSec = fSectors->GetAlpha() * seed->GetRelativeSector() + fSectors->GetAlphaShift();
if (alphaSec >= TMath::Pi()) alphaSec -= 2.*TMath::Pi();
if (alphaSec < -TMath::Pi()) alphaSec += 2.*TMath::Pi();
seed->Rotate(alphaSec - alpha);
seed->SetPoolID(fLastSeedID);
seed->SetIsSeeding(kTRUE);
seed->SetSeed1(nup-1);
seed->SetSeed2(nup-2);
seed->SetSeedType(0);
seed->SetFirstPoint(-1);
seed->SetLastPoint(-1);
seeds->AddLast(seed);
index++;
}
fSectors = fOuterSec;
TrackFollowingHLT(seeds );
nTr = seeds->GetEntriesFast();
for( int itr=0; itr<nTr; itr++ ){
AliTPCseed * seed = (AliTPCseed*) seeds->UncheckedAt(itr);
if( !seed ) continue;
Int_t foundable,found,shared;
seed->GetClusterStatistic(0,nup, found, foundable, shared, kTRUE);
seed->SetNFoundable(foundable);
if (seed->GetNumberOfClusters()<30 ||
seed->GetNumberOfClusters() < seed->GetNFoundable()*0.6 ||
seed->GetNShared()>0.4*seed->GetNumberOfClusters() ) {
MarkSeedFree(seeds->RemoveAt(itr));
continue;
}
for( int ir=0; ir<nup; ir++){
AliTPCclusterMI *c = seed->GetClusterPointer(ir);
if( c ) c->Use(10);
}
}
std::cout<<"\n\nHLT tracks left: "<<seeds->GetEntries()<<" out of "<<hltEvent->GetNumberOfTracks()<<endl<<endl;
return seeds;
}
void AliTPCtracker::FillClusterOccupancyInfo()
{
AliESDfriend* esdFriend = static_cast<AliESDfriend*>(fEvent->FindListObject("AliESDfriend"));
if (!esdFriend) return;
for (Int_t isector=0; isector<18; isector++){
AliTPCtrackerSector &iroc = fInnerSec[isector];
AliTPCtrackerSector &oroc = fOuterSec[isector];
esdFriend->SetNclustersTPC(isector, iroc.GetNClInSector(0));
esdFriend->SetNclustersTPC(isector+18,iroc.GetNClInSector(1));
esdFriend->SetNclustersTPC(isector+36,oroc.GetNClInSector(0));
esdFriend->SetNclustersTPC(isector+54,oroc.GetNClInSector(1));
esdFriend->SetNclustersTPCused(isector, iroc.GetNClUsedInSector(0));
esdFriend->SetNclustersTPCused(isector+18, iroc.GetNClUsedInSector(1));
esdFriend->SetNclustersTPCused(isector+36, oroc.GetNClUsedInSector(0));
esdFriend->SetNclustersTPCused(isector+54, oroc.GetNClUsedInSector(1));
}
}