#include <TFile.h>
#include <TObjArray.h>
#include <TTree.h>
+#include <TGraphErrors.h>
#include "AliLog.h"
#include "AliComplexCluster.h"
#include "AliESDEvent.h"
#include "AliTrackPointArray.h"
#include "TRandom.h"
#include "AliTPCcalibDB.h"
+#include "AliTPCcalibDButil.h"
#include "AliTPCTransform.h"
#include "AliTPCClusterParam.h"
//
// default constructor
//
+ for (Int_t irow=0; irow<200; irow++){
+ fXRow[irow]=0;
+ fYMax[irow]=0;
+ fPadLength[irow]=0;
+ }
+
}
//_____________________________________________________________________
AliTPCclusterMI* c =track->GetCurrentCluster();
- if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
-
+ if (accept > 0) //sign not accepted clusters
+ track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
+ else // unsign accpeted clusters
+ track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
UInt_t i = track->GetCurrentClusterIndex1();
Int_t sec=(i&0xff000000)>>24;
seed->SetNFoundable(seed->GetNFoundable()-1);
return 2;
}
+
+ if (AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 3 ||
+ AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 4 ) {
+ if(!AliTPCReconstructor::GetRecoParam()->GetUseOnePadCluster())
+ if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
+ return 2;
+ }
return 0;
}
// dummy copy constructor
//------------------------------------------------------------------
fOutput=t.fOutput;
+ for (Int_t irow=0; irow<200; irow++){
+ fXRow[irow]=0;
+ fYMax[irow]=0;
+ fPadLength[irow]=0;
+ }
+
}
AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
{
//
//
//fill esds using updated tracks
- if (fEvent){
+ if (!fEvent) return;
+
// write tracks to the event
// store index of the track
Int_t nseed=arr->GetEntriesFast();
kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
}
}
- // << account for suppressed tracks in the kink indices (RS)
- }
- AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
+ // << account for suppressed tracks in the kink indices (RS)
+ AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
+
}
Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
{
- //
+ // load clusters
//
fInput = tree;
return LoadClusters();
void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
//
+ // transformation
//
- //
- AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+ AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
+ AliTPCTransform *transform = calibDB->GetTransform() ;
if (!transform) {
AliFatal("Tranformations not in calibDB");
+ return;
}
transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
cluster->SetY(x[1]);
cluster->SetZ(x[2]);
// The old stuff:
-
//
//
//
//if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
- TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
- //TGeoHMatrix mat;
- Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
- Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
- if (mat) mat->LocalToMaster(pos,posC);
- else{
- // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
- }
- cluster->SetX(posC[0]);
- cluster->SetY(posC[1]);
- cluster->SetZ(posC[2]);
+ if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
+ TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
+ //TGeoHMatrix mat;
+ Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
+ Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
+ if (mat) mat->LocalToMaster(pos,posC);
+ else{
+ // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
+ }
+ cluster->SetX(posC[0]);
+ cluster->SetY(posC[1]);
+ cluster->SetZ(posC[2]);
+ }
}
//_____________________________________________________________________________
if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
//
if (tpcindex==-1) return 0; //track in dead zone
- if (tpcindex>0){ //
+ if (tpcindex >= 0){ //
cl = t.GetClusterPointer(nr);
if ( (cl==0) ) cl = GetClusterMI(tpcindex);
t.SetCurrentClusterIndex1(tpcindex);
t.SetNFoundable(t.GetNFoundable()+1);
UpdateTrack(&t,accept);
return 1;
- }
+ }
+ else { // Remove old cluster from track
+ t.SetClusterIndex(nr, -3);
+ t.SetClusterPointer(nr, 0);
+ }
}
}
if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
//
UInt_t index=0;
// if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
- Double_t y=t.GetYat(x);
+ if (!t.PropagateTo(x)) {
+ if (fIteration==0) t.SetRemoval(10);
+ return 0;
+ }
+ Double_t y = t.GetY();
if (TMath::Abs(y)>ymax){
if (y > ymax) {
t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
if (!t.Rotate(-fSectors->GetAlpha()))
return 0;
}
- //return 1;
+ if (!t.PropagateTo(x)) {
+ if (fIteration==0) t.SetRemoval(10);
+ return 0;
+ }
+ y = t.GetY();
}
//
- if (!t.PropagateTo(x)) {
- if (fIteration==0) t.SetRemoval(10);
- return 0;
- }
- y=t.GetY();
Double_t z=t.GetZ();
//
// This function tries to find a track prolongation to next pad row
//-----------------------------------------------------------------
t.SetCurrentCluster(0);
- t.SetCurrentClusterIndex1(0);
+ t.SetCurrentClusterIndex1(-3);
Double_t xt=t.GetX();
Int_t row = GetRowNumber(xt)-1;
} else {
if (fIteration==0){
- if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
- if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
+ if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
+ if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
}
//_____________________________________________________________________________
-Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
+Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
//-----------------------------------------------------------------
// This function tries to find a track prolongation.
//-----------------------------------------------------------------
Double_t xt=t.GetX();
//
- Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
+ Double_t alpha=t.GetAlpha();
if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
if (alpha < 0. ) alpha += 2.*TMath::Pi();
//
- t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
+ t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha())%fN);
- Int_t first = GetRowNumber(xt)-1;
+ Int_t first = GetRowNumber(xt);
+ if (!fromSeeds)
+ first -= step;
for (Int_t nr= first; nr>=rf; nr-=step) {
// update kink info
if (t.GetKinkIndexes()[0]>0){
-Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
+Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
//-----------------------------------------------------------------
// This function tries to find a track prolongation.
//-----------------------------------------------------------------
//
Double_t xt=t.GetX();
- Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
+ Double_t alpha=t.GetAlpha();
if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
if (alpha < 0. ) alpha += 2.*TMath::Pi();
- t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
+ t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha())%fN);
Int_t first = t.GetFirstPoint();
- if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
+ 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++) {
Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
{
- //
+ // overlapping factor
//
sum1=0;
sum2=0;
void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
{
- //
+ // shared clusters
//
Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
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++){
+ for (Int_t j=0;j<12;j++){
pt->SetOverlapLabel(j,0);
}
}
if (!pt) continue;
//
if (quality[trackindex]<0){
- if (pt) {
delete arr->RemoveAt(trackindex);
- }
- else{
- arr->RemoveAt(trackindex);
- }
- continue;
+ continue;
}
//
//
// back propagation of ESD tracks
//
//return 0;
+ if (!event) return 0;
const Int_t kMaxFriendTracks=2000;
fEvent = event;
+ // extract correction object for multiplicity dependence of dEdx
+ TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
+
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+ if (!transform) {
+ AliFatal("Tranformations not in RefitInward");
+ return 0;
+ }
+ transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
+ const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
+ Int_t nContribut = event->GetNumberOfTracks();
+ TGraphErrors * graphMultDependenceDeDx = 0x0;
+ if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
+ if (recoParam->GetUseTotCharge()) {
+ graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
+ } else {
+ graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
+ }
+ }
+ //
ReadSeeds(event,2);
fIteration=2;
//PrepareForProlongation(fSeeds,1);
AliESDtrack *esd=event->GetTrack(i);
if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
+ AliInfo("Refitting track");
AliExternalTrackParam paramIn;
AliExternalTrackParam paramOut;
Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
esd->SetTPCClusterMap(seed->GetClusterMap());
esd->SetTPCSharedMap(seed->GetSharedMap());
//
- if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
+ if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Crefit"<<
"Esd.="<<esd<<
Int_t ndedx = seed->GetNCDEDX(0);
Float_t sdedx = seed->GetSDEDX(0);
Float_t dedx = seed->GetdEdx();
+ // apply mutliplicity dependent dEdx correction if available
+ if (graphMultDependenceDeDx) {
+ Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
+ dedx += (1 - corrGain)*50.; // MIP is normalized to 50
+ }
esd->SetTPCsignal(dedx, sdedx, ndedx);
//
// add seed to the esd track in Calib level
//
// back propagation of ESD tracks
//
-
+ if (!event) return 0;
fEvent = event;
fIteration = 1;
ReadSeeds(event,1);
seed->UpdatePoints();
AddCovariance(seed);
AliESDtrack *esd=event->GetTrack(i);
+ if (!esd) continue; //never happen
if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
AliExternalTrackParam paramIn;
AliExternalTrackParam paramOut;
fSeeds =0;
}
-void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
+void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
{
//
//read seeds from the event
// rotate to the local coordinate system
//
fSectors=fInnerSec; fN=fkNIS;
- Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
+ Double_t alpha=seed->GetAlpha();
if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
if (alpha < 0. ) alpha += 2.*TMath::Pi();
Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
+ alpha-=seed->GetAlpha();
if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
- alpha-=seed->GetAlpha();
- if (!seed->Rotate(alpha)) {
- delete seed;
- continue;
+ if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
+ AliWarning(Form("Rotating track over %f",alpha));
+ if (!seed->Rotate(alpha)) {
+ delete seed;
+ continue;
+ }
}
seed->SetESD(esd);
// sign clusters
ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
//
//
- Double_t xyz[3][3];
- Int_t row[3],sec[3]={0,0,0};
+ Double_t xyz[3][3]={{0}};
+ Int_t row[3]={0},sec[3]={0,0,0};
//
// find track row position at given ratio of the length
Int_t index=-1;
}
delete [] helixes;
delete [] xm;
+ delete [] dz0;
+ delete [] dz1;
if (AliTPCReconstructor::StreamLevel()>1) {
AliInfo("Time for curling tracks removal DEBUGGING MC");
timer.Print();
//
TStopwatch timer;
timer.Start();
- Double_t phase[2][2],radius[2];
+ Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
+
//
// Find tracks
//
Int_t ncandidates =0;
Int_t nall =0;
Int_t ntracks=0;
- Double_t phase[2][2],radius[2];
+ Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
//
// Find circling track
AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
if (!kink0) continue;
//
- for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
+ 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;
timer.Print();
}
-void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *const esd)
-{
- //
- // find V0s
- //
- //
- TObjArray *tpcv0s = new TObjArray(100000);
- Int_t nentries = array->GetEntriesFast();
- AliHelix *helixes = new AliHelix[nentries];
- Int_t *sign = new Int_t[nentries];
- Float_t *alpha = new Float_t[nentries];
- Float_t *z0 = new Float_t[nentries];
- Float_t *dca = new Float_t[nentries];
- Float_t *sdcar = new Float_t[nentries];
- Float_t *cdcar = new Float_t[nentries];
- Float_t *pulldcar = new Float_t[nentries];
- Float_t *pulldcaz = new Float_t[nentries];
- Float_t *pulldca = new Float_t[nentries];
- Bool_t *isPrim = new Bool_t[nentries];
- const AliESDVertex * primvertex = esd->GetVertex();
- Double_t zvertex = primvertex->GetZv();
- //
- // nentries = array->GetEntriesFast();
- //
- for (Int_t i=0;i<nentries;i++){
- sign[i]=0;
- isPrim[i]=0;
- AliTPCseed* track = (AliTPCseed*)array->At(i);
- if (!track) continue;
- track->GetV0Indexes()[0] = 0; //rest v0 indexes
- track->GetV0Indexes()[1] = 0; //rest v0 indexes
- track->GetV0Indexes()[2] = 0; //rest v0 indexes
- //
- alpha[i] = track->GetAlpha();
- new (&helixes[i]) AliHelix(*track);
- Double_t xyz[3];
- helixes[i].Evaluate(0,xyz);
- sign[i] = (track->GetC()>0) ? -1:1;
- Double_t x,y,z;
- x=160;
- z0[i]=1000;
- if (track->GetProlongation(0,y,z)) z0[i] = z;
- dca[i] = track->GetD(0,0);
- //
- // dca error parrameterezation + pulls
- //
- sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
- if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
- cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
- pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
- pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
- pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
- if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
- if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
- }
- if (track->TPCrPID(4)>0.5) {
- if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
- }
- if (track->TPCrPID(0)>0.4) {
- isPrim[i]=kFALSE; //electron no sigma cut
- }
- }
- //
- //
- TStopwatch timer;
- timer.Start();
- Int_t ncandidates =0;
- Int_t nall =0;
- Int_t ntracks=0;
- Double_t phase[2][2],radius[2];
- //
- // Finf V0s loop
- //
- //
- // //
- Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
- AliV0 vertex;
- Double_t cradius0 = 10*10;
- Double_t cradius1 = 200*200;
- Double_t cdist1=3.;
- Double_t cdist2=4.;
- Double_t cpointAngle = 0.95;
- //
- Double_t delta[2]={10000,10000};
- for (Int_t i =0;i<nentries;i++){
- if (sign[i]==0) continue;
- AliTPCseed * track0 = (AliTPCseed*)array->At(i);
- if (!track0) continue;
- if (AliTPCReconstructor::StreamLevel()>1){
- TTreeSRedirector &cstream = *fDebugStreamer;
- cstream<<"Tracks"<<
- "Tr0.="<<track0<<
- "dca="<<dca[i]<<
- "z0="<<z0[i]<<
- "zvertex="<<zvertex<<
- "sdcar0="<<sdcar[i]<<
- "cdcar0="<<cdcar[i]<<
- "pulldcar0="<<pulldcar[i]<<
- "pulldcaz0="<<pulldcaz[i]<<
- "pulldca0="<<pulldca[i]<<
- "isPrim="<<isPrim[i]<<
- "\n";
- }
- //
- if (track0->GetSigned1Pt()<0) continue;
- if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
- //
- if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
- ntracks++;
- // debug output
-
-
- for (Int_t j =0;j<nentries;j++){
- AliTPCseed * track1 = (AliTPCseed*)array->At(j);
- if (!track1) continue;
- if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
- if (sign[j]*sign[i]>0) continue;
- if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
- if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
- nall++;
- //
- // DCA to prim vertex cut
- //
- //
- delta[0]=10000;
- delta[1]=10000;
-
- Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
- if (npoints<1) continue;
- Int_t iclosest=0;
- // cuts on radius
- if (npoints==1){
- if (radius[0]<cradius0||radius[0]>cradius1) continue;
- helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
- if (delta[0]>cdist1) continue;
- }
- else{
- if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
- helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
- helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
- if (delta[1]<delta[0]) iclosest=1;
- if (delta[iclosest]>cdist1) continue;
- }
- helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
- if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
- //
- Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
- if (pointAngle<cpointAngle) continue;
- //
- Bool_t isGamma = kFALSE;
- vertex.SetParamP(*track0); //track0 - plus
- vertex.SetParamN(*track1); //track1 - minus
- vertex.Update(fprimvertex);
- if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
- Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
- //continue;
- if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
- //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
- if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
- Float_t sigmae = 0.15*0.15;
- if (vertex.GetRr()<80)
- sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
- sigmae+= TMath::Sqrt(sigmae);
- //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
- if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
- Float_t densb0=0,densb1=0,densa0=0,densa1=0;
- Int_t row0 = GetRowNumber(vertex.GetRr());
- if (row0>15){
- //Bo: if (vertex.GetDist2()>0.2) continue;
- if (vertex.GetDcaV0Daughters()>0.2) continue;
- densb0 = track0->Density2(0,row0-5);
- densb1 = track1->Density2(0,row0-5);
- if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
- densa0 = track0->Density2(row0+5,row0+40);
- densa1 = track1->Density2(row0+5,row0+40);
- if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
- }
- else{
- densa0 = track0->Density2(0,40); //cluster density
- densa1 = track1->Density2(0,40); //cluster density
- if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
- }
-//Bo: vertex.SetLab(0,track0->GetLabel());
-//Bo: vertex.SetLab(1,track1->GetLabel());
- vertex.SetChi2After((densa0+densa1)*0.5);
- vertex.SetChi2Before((densb0+densb1)*0.5);
- vertex.SetIndex(0,i);
- vertex.SetIndex(1,j);
-//Bo: vertex.SetStatus(1); // TPC v0 candidate
- vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
-//Bo: vertex.SetRp(track0->TPCrPIDs());
-//Bo: vertex.SetRm(track1->TPCrPIDs());
- tpcv0s->AddLast(new AliESDv0(vertex));
- ncandidates++;
- {
- Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
- Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
- Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
- if (AliTPCReconstructor::StreamLevel()>1) {
- Int_t lab0=track0->GetLabel();
- Int_t lab1=track1->GetLabel();
- Char_t c0=track0->GetCircular();
- Char_t c1=track1->GetCircular();
- TTreeSRedirector &cstream = *fDebugStreamer;
- cstream<<"V0"<<
- "Event="<<eventNr<<
- "vertex.="<<&vertex<<
- "Tr0.="<<track0<<
- "lab0="<<lab0<<
- "Helix0.="<<&helixes[i]<<
- "Tr1.="<<track1<<
- "lab1="<<lab1<<
- "Helix1.="<<&helixes[j]<<
- "pointAngle="<<pointAngle<<
- "pointAngle2="<<pointAngle2<<
- "dca0="<<dca[i]<<
- "dca1="<<dca[j]<<
- "z0="<<z0[i]<<
- "z1="<<z0[j]<<
- "zvertex="<<zvertex<<
- "circular0="<<c0<<
- "circular1="<<c1<<
- "npoints="<<npoints<<
- "radius0="<<radius[0]<<
- "delta0="<<delta[0]<<
- "radius1="<<radius[1]<<
- "delta1="<<delta[1]<<
- "radiusm="<<radiusm<<
- "deltam="<<deltam<<
- "sdcar0="<<sdcar[i]<<
- "sdcar1="<<sdcar[j]<<
- "cdcar0="<<cdcar[i]<<
- "cdcar1="<<cdcar[j]<<
- "pulldcar0="<<pulldcar[i]<<
- "pulldcar1="<<pulldcar[j]<<
- "pulldcaz0="<<pulldcaz[i]<<
- "pulldcaz1="<<pulldcaz[j]<<
- "pulldca0="<<pulldca[i]<<
- "pulldca1="<<pulldca[j]<<
- "densb0="<<densb0<<
- "densb1="<<densb1<<
- "densa0="<<densa0<<
- "densa1="<<densa1<<
- "sigmae="<<sigmae<<
- "\n";}
- }
- }
- }
- Float_t *quality = new Float_t[ncandidates];
- Int_t *indexes = new Int_t[ncandidates];
- Int_t naccepted =0;
- for (Int_t i=0;i<ncandidates;i++){
- quality[i] = 0;
- AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
- quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
- // quality[i] /= (0.5+v0->GetDist2());
- // quality[i] *= v0->GetChi2After(); //density factor
-
- Int_t index0 = v0->GetIndex(0);
- Int_t index1 = v0->GetIndex(1);
- //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
- Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
-
-
-
- AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
- AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
- if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
- if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
- }
-
- TMath::Sort(ncandidates,quality,indexes,kTRUE);
- //
- //
- for (Int_t i=0;i<ncandidates;i++){
- AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
- if (!v0) continue;
- Int_t index0 = v0->GetIndex(0);
- Int_t index1 = v0->GetIndex(1);
- AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
- AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
- if (!track0||!track1) {
- printf("Bug\n");
- continue;
- }
- Bool_t accept =kTRUE; //default accept
- Int_t *v0indexes0 = track0->GetV0Indexes();
- Int_t *v0indexes1 = track1->GetV0Indexes();
- //
- Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
- Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
- if (v0indexes0[1]!=0) order0 =2;
- if (v0indexes1[1]!=0) order1 =2;
- //
- if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
- if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
- //
- AliESDv0 * v02 = v0;
- if (accept){
- //Bo: v0->SetOrder(0,order0);
- //Bo: v0->SetOrder(1,order1);
- //Bo: v0->SetOrder(1,order0+order1);
- v0->SetOnFlyStatus(kTRUE);
- Int_t index = esd->AddV0(v0);
- v02 = esd->GetV0(index);
- v0indexes0[order0]=index;
- v0indexes1[order1]=index;
- naccepted++;
- }
- {
- Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
- if (AliTPCReconstructor::StreamLevel()>1) {
- Int_t lab0=track0->GetLabel();
- Int_t lab1=track1->GetLabel();
- TTreeSRedirector &cstream = *fDebugStreamer;
- cstream<<"V02"<<
- "Event="<<eventNr<<
- "vertex.="<<v0<<
- "vertex2.="<<v02<<
- "Tr0.="<<track0<<
- "lab0="<<lab0<<
- "Tr1.="<<track1<<
- "lab1="<<lab1<<
- "dca0="<<dca[index0]<<
- "dca1="<<dca[index1]<<
- "order0="<<order0<<
- "order1="<<order1<<
- "accept="<<accept<<
- "quality="<<quality[i]<<
- "pulldca0="<<pulldca[index0]<<
- "pulldca1="<<pulldca[index1]<<
- "index="<<i<<
- "\n";}
- }
- }
-
-
- //
- //
- delete []quality;
- delete []indexes;
-//
- delete [] isPrim;
- delete [] pulldca;
- delete [] pulldcaz;
- delete [] pulldcar;
- delete [] cdcar;
- delete [] sdcar;
- delete [] dca;
- //
- delete[] z0;
- delete[] alpha;
- delete[] sign;
- delete[] helixes;
- printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
- timer.Print();
-}
Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
{
Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
{
//
+ // clusters to tracks
+
if (fSeeds) DeleteSeeds();
- fEvent = esd;
+ fEvent = esd;
+
+ AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
+ transform->SetCurrentTimeStamp( esd->GetTimeStamp());
+ transform->SetCurrentRun(esd->GetRunNumber());
+
Clusters2Tracks();
if (!fSeeds) return 1;
FillESD(fSeeds);
fIteration = 2;
PrepareForProlongation(fSeeds,5.);
- PropagateForward2(fSeeds);
+ PropagateForard2(fSeeds);
printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
// RemoveUsed(fSeeds,0.7,0.7,6);
TObjArray * AliTPCtrackerMI::Tracking()
{
- //
+ // tracking
//
if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
TStopwatch timer;
}
}
-void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
+void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
{
//
//
Float_t angle2 = pt->GetAlpha();
if (TMath::Abs(angle1-angle2)>0.001){
- pt->Rotate(angle1-angle2);
+ if (!pt->Rotate(angle1-angle2)) return;
//angle2 = pt->GetAlpha();
//pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
//if (pt->GetAlpha()<0)
}
-Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
+Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
{
//
// make back propagation
fSectors = fInnerSec;
//FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
//fSectors = fOuterSec;
- FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
+ FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
//if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
// Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
// FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
pt->SetFirstPoint(kink->GetTPCRow0());
fSectors = fInnerSec;
- FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
+ FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
}
CookLabel(pt,0.3);
}
}
-Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
+Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
{
//
// make forward propagation
for (Int_t i=0;i<nseed;i++){
AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
if (pt) {
- FollowProlongation(*pt,0);
+ FollowProlongation(*pt,0,1,1);
CookLabel(pt,0.3);
}
}
- return 0;
+ return 0;
}
void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
{
- //
+ // gets cluster shape
//
AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
Int_t tail=Int_t(0.10*noc);
max=0;
Int_t ind=0;
- for (i=1; i<=160&&ind<tail; i++) {
+ for (i=1; i<160&&ind<tail; i++) {
// AliTPCclusterMI *c=clusters[noc-i];
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
}
}
noc = current;
- if (noc<5) return -1;
+ //if (noc<5) return -1;
Int_t lab=123456789;
for (i=0; i<noc; i++) {
AliTPCclusterMI *c=clusters[i];
Int_t tail=Int_t(0.10*noc);
max=0;
Int_t ind=0;
- for (i=1; i<=160&&ind<tail; i++) {
+ for (i=1; i<160&&ind<tail; i++) {
// AliTPCclusterMI *c=clusters[noc-i];
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
//
- // Adding systematic error
+ // Adding systematic error estimate to the covariance matrix
// !!!! the systematic error for element 4 is in 1/cm not in pt
+ //
+ const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
+ //
+ // use only the diagonal part if not specified otherwise
+ if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
+ //
+ Double_t *covarS= (Double_t*)seed->GetCovariance();
+ Double_t factor[5]={1,1,1,1,1};
+ Double_t facC = AliTracker::GetBz()*kB2C;
+ 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] + facC*facC*param[4]*param[4])/covarS[14]));
+ // 0
+ // 1 2
+ // 3 4 5
+ // 6 7 8 9
+ // 10 11 12 13 14
+ for (Int_t i=0; i<5; i++){
+ for (Int_t j=i; j<5; j++){
+ Int_t index=seed->GetIndex(i,j);
+ covarS[index]*=factor[i]*factor[j];
+ }
+ }
+}
+
+void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
+ //
+ // Adding systematic error - as additive factor without correlation
+ //
const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
Double_t *covarIn= (Double_t*)seed->GetCovariance();
Double_t covar[15];
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff