//
// The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
// To enable storage of the TPC tracks in the ESD friend track
-// use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
+// use AliTPCReconstructor::SetStreamLevel(n);
//
// The debug level - different procedure produce tree for numerical debugging
// To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
// There are several places in the code which can be numerically debuged
// This code is keeped in order to enable code development and to check the calibration implementtion
//
-// 1. ErrParam stream (Log level 9) - dump information about
+// 1. ErrParam stream - dump information about
// 1.a) cluster
// 2.a) cluster error estimate
// 3.a) cluster shape estimate
//
//
+// Debug streamer levels:
+//
//-------------------------------------------------------
#include <TFile.h>
#include <TObjArray.h>
#include <TTree.h>
+#include <TGraphErrors.h>
#include "AliLog.h"
#include "AliComplexCluster.h"
#include "AliESDEvent.h"
#include "AliTPCtrackerMI.h"
#include "TStopwatch.h"
#include "AliTPCReconstructor.h"
-#include "AliPID.h"
-#include "TTreeStream.h"
#include "AliAlignObj.h"
#include "AliTrackPointArray.h"
#include "TRandom.h"
#include "AliTPCcalibDB.h"
+#include "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;
+ }
+
}
//_____________________________________________________________________
Double_t sdistancey2 = sy2+seed->GetSigmaY2();
Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
- Double_t dz=seed->GetCurrentCluster()->GetY()-zt;
+ 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)*
Double_t rdistance2 = rdistancey2+rdistancez2;
//Int_t accept =0;
- if (AliTPCReconstructor::StreamLevel()>1 && seed->GetNumberOfClusters()>20) {
+ if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
Float_t rmsy2 = seed->GetCurrentSigmaY2();
Float_t rmsz2 = seed->GetCurrentSigmaZ2();
Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
- if (AliTPCReconstructor::StreamLevel()>0) {
+ if (AliTPCReconstructor::StreamLevel()>2) {
(*fDebugStreamer)<<"ErrParam"<<
"Cl.="<<cluster<<
"T.="<<¶m<<
"rmsz2p30="<<rmsz2p30<<
"rmsy2p30R="<<rmsy2p30R<<
"rmsz2p30R="<<rmsz2p30R<<
+ // normalize distance -
+ "rdisty="<<rdistancey2<<
+ "rdistz="<<rdistancez2<<
+ "rdist="<<rdistance2<< //
"\n";
}
}
+
//_____________________________________________________________________________
AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
AliTracker(),
// 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)
- }
- printf("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()));
+
}
erry2+=0.5; // edge cluster
}
erry2*=erry2;
- erry2+=0.015*0.015;
seed->SetErrorY2(erry2);
//
return erry2;
if (ctype<0) {
errz2+=0.5; // edge cluster
}
- errz2*=errz2;
- errz2+=0.015*0.015;
+ errz2*=errz2;
seed->SetErrorZ2(errz2);
//
return errz2;
y2 -=y1;
//
Double_t det = x3*y2-x2*y3;
- if (det==0) {
+ if (TMath::Abs(det)<1e-10){
return 100;
}
//
y2 -=y1;
//
Double_t det = x3*y2-x2*y3;
- if (det==0) {
+ if (TMath::Abs(det)<1e-10) {
return 100;
}
//
return angle2;
}
-Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
+Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
{//-----------------------------------------------------------------
// This function find proloncation of a track to a reference plane x=x2.
//-----------------------------------------------------------------
return kTRUE;
}
-Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
+Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
{
//
//
{
//
// load clusters to the memory
- AliTPCClustersRow *clrow = 0x0;
+ AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
Int_t lower = arr->LowerBound();
Int_t entries = arr->GetEntriesFast();
for (Int_t j=0;j<tpcrow->GetN2();++j)
tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
}
+ clrow->GetArray()->Clear("C");
}
//
delete clrow;
{
//
// load clusters to the memory
- AliTPCClustersRow *clrow= new AliTPCClustersRow;
- clrow->SetClass("AliTPCclusterMI");
- clrow->SetArray(0);
- clrow->GetArray()->ExpandCreateFast(10000);
+ AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
//
// TTree * tree = fClustersArray.GetTree();
//
//
//
- 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()};
//
// in debug mode check the transformation
//
- if (AliTPCReconstructor::StreamLevel()>1) {
+ if (AliTPCReconstructor::StreamLevel()>2) {
Float_t gx[3];
cluster->GetGlobalXYZ(gx);
Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
- Double_t xt=GetXrow(cluster->GetRow());
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Transform"<<
"event="<<event<<
- "xt="<<xt<<
"x0="<<x[0]<<
"x1="<<x[1]<<
"x2="<<x[2]<<
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
- }
- if (AliTPCReconstructor::StreamLevel()>1) {
- Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
- Double_t xt=GetXrow(cluster->GetRow());
- TTreeSRedirector &cstream = *fDebugStreamer;
- cstream<<"Transform2"<<
- "event="<<event<<
- "xt="<<xt<<
- "x0="<<pos[0]<<
- "x1="<<pos[1]<<
- "x2="<<pos[2]<<
- "ax0="<<posC[0]<<
- "ax1="<<posC[1]<<
- "ax2="<<posC[2]<<
- "Cl.="<<cluster<<
- "\n";
+ 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]);
}
-
- cluster->SetX(posC[0]);
- cluster->SetY(posC[1]);
- cluster->SetZ(posC[2]);
}
//_____________________________________________________________________________
Int_t lastpoint = 160;
//
// if (firstpoint>=lastpoint-5) return;;
- Int_t firstShared=lastpoint, lastShared=firstpoint;
- Int_t firstRow=lastpoint, lastRow=firstpoint;
+
for (Int_t i=firstpoint;i<lastpoint;i++){
// if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
//
for (Int_t i=firstpoint;i<lastpoint;i++){
// if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
- 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) {
AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
- if (i<firstShared) firstShared=i;
- if (i>lastShared) lastShared=i;
if (s1->IsActive()&&s2->IsActive()){
p1->SetShared(kTRUE);
p2->SetShared(kTRUE);
}
//
if (sumshared>cutN0){
- Float_t dedx1S=s1->CookdEdx(0.05,0.7, firstShared, lastShared);
- Float_t dedx2S=s2->CookdEdx(0.05,0.7, firstShared, lastShared);
- Float_t dedx1A=s1->CookdEdx(0.05,0.7, 0, 160);
- Float_t dedx2A=s2->CookdEdx(0.05,0.7, 0, 160);
- Int_t rowM= (firstShared+lastShared)/2;
- Double_t xmid = GetXrow(rowM);
- AliExternalTrackParam p1(*s1);
- AliExternalTrackParam p2(*s2);
- p2.Rotate(p1.GetAlpha());
- p1.PropagateTo(xmid,GetBz());
- p2.PropagateTo(xmid,GetBz());
for (Int_t i=0;i<4;i++){
if (s1->GetOverlapLabel(3*i)==0){
s1->SetOverlapLabel(3*i, s2->GetLabel());
break;
}
}
- if (fDebugStreamer){
- TTreeSRedirector &cstream = *fDebugStreamer;
- //
- cstream<<"Shared"<<
- "iter="<<fIteration<<
- "shared="<<sumshared<<
- "firstS="<<firstShared<<
- "lastS="<<lastShared<<
- "firstR="<<firstRow<<
- "lastR="<<lastRow<<
- "dedx1S="<<dedx1S<<
- "dedx2S="<<dedx2S<<
- "dedx1A="<<dedx1A<<
- "dedx2A="<<dedx2A<<
- "p0.="<<&p1<<
- "p1.="<<&p2<<
- "Tr0.="<<s1<<
- "Tr1.="<<s2<<
- "\n";
- }
}
}
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->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())%18)>1) continue;
+ if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
// if (pt2){
if (pt2->GetRemoval()<=10) {
//if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
if (!pt) continue;
//
if (quality[trackindex]<0){
- if (pt) {
delete arr->RemoveAt(trackindex);
- }
- else{
- arr->RemoveAt(trackindex);
- }
- continue;
+ continue;
}
//
//
//
if (Float_t(shared+1)/Float_t(found+1)>factor){
if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
- if( AliTPCReconstructor::StreamLevel()>15){
+ if( AliTPCReconstructor::StreamLevel()>3){
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"RemoveUsed"<<
"iter="<<fIteration<<
}
if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
- if( AliTPCReconstructor::StreamLevel()>15){
+ if( AliTPCReconstructor::StreamLevel()>3){
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"RemoveShort"<<
"iter="<<fIteration<<
delete []indexes;
}
+void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
+{
+ //
+ // Dump clusters after reco
+ // signed and unsigned cluster can be visualized
+ // 1. Unsign all cluster
+ // 2. Sign all used clusters
+ // 3. Dump clusters
+ UnsignClusters();
+ Int_t nseed = trackArray->GetEntries();
+ for (Int_t i=0; 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); // kink
+ c->Use(10); // by default usage 10
+ }
+ }
+ //
+
+ for (Int_t sec=0;sec<fkNIS;sec++){
+ for (Int_t row=0;row<fInnerSec->GetNRows();row++){
+ AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
+ for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
+ Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
+ (*fDebugStreamer)<<"clDump"<<
+ "iter="<<iter<<
+ "cl.="<<&cl[icl]<<
+ "gx0="<<gx[0]<<
+ "gx1="<<gx[1]<<
+ "gx2="<<gx[2]<<
+ "\n";
+ }
+ cl = fInnerSec[sec][row].GetClusters2();
+ for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
+ Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
+ (*fDebugStreamer)<<"clDump"<<
+ "iter="<<iter<<
+ "cl.="<<&cl[icl]<<
+ "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++){
+ AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
+ for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
+ Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
+ (*fDebugStreamer)<<"clDump"<<
+ "iter="<<iter<<
+ "cl.="<<&cl[icl]<<
+ "gx0="<<gx[0]<<
+ "gx1="<<gx[1]<<
+ "gx2="<<gx[2]<<
+ "\n";
+ }
+ cl = fOuterSec[sec][row].GetClusters2();
+ for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
+ Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
+ (*fDebugStreamer)<<"clDump"<<
+ "iter="<<iter<<
+ "cl.="<<&cl[icl]<<
+ "gx0="<<gx[0]<<
+ "gx1="<<gx[1]<<
+ "gx2="<<gx[2]<<
+ "\n";
+ }
+ }
+ }
+
+}
void AliTPCtrackerMI::UnsignClusters()
{
//
// 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);
PropagateForward2(fSeeds);
RemoveUsed2(fSeeds,0.4,0.4,20);
- FindSplitted(fSeeds); // find multi found tracks
- //
+
TObjArray arraySeed(fSeeds->GetEntries());
for (Int_t i=0;i<fSeeds->GetEntries();i++) {
arraySeed.AddAt(fSeeds->At(i),i);
}
SignShared(&arraySeed);
// FindCurling(fSeeds, event,2); // find multi found tracks
- if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
+ FindSplitted(fSeeds, event,2); // find multi found tracks
+ if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
Int_t ntracks=0;
Int_t nseed = fSeeds->GetEntriesFast();
AliExternalTrackParam paramIn;
AliExternalTrackParam paramOut;
Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
- if (AliTPCReconstructor::StreamLevel()>0) {
+ if (AliTPCReconstructor::StreamLevel()>2) {
(*fDebugStreamer)<<"RecoverIn"<<
"seed.="<<seed<<
"esd.="<<esd<<
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
}
}
//FindKinks(fSeeds,event);
+ if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
Info("RefitInward","Number of refitted tracks %d",ntracks);
return 0;
}
//
// back propagation of ESD tracks
//
-
+ if (!event) return 0;
fEvent = event;
fIteration = 1;
ReadSeeds(event,1);
PropagateBack(fSeeds);
RemoveUsed2(fSeeds,0.4,0.4,20);
//FindCurling(fSeeds, fEvent,1);
- FindSplitted(fSeeds); // find multi found tracks
- if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
- SignShared(fSeeds);
+ FindSplitted(fSeeds, event,1); // find multi found tracks
+ if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
+
//
Int_t nseed = fSeeds->GetEntriesFast();
Int_t ntracks=0;
seed->UpdatePoints();
AddCovariance(seed);
AliESDtrack *esd=event->GetTrack(i);
+ if (!esd) continue; //never happen
if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
AliExternalTrackParam paramIn;
AliExternalTrackParam paramOut;
Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
- if (AliTPCReconstructor::StreamLevel()>0) {
+ if (AliTPCReconstructor::StreamLevel()>2) {
(*fDebugStreamer)<<"RecoverBack"<<
"seed.="<<seed<<
"esd.="<<esd<<
}
}
}
+ if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
//FindKinks(fSeeds,event);
Info("PropagateBack","Number of back propagated tracks %d",ntracks);
fEvent =0;
fSeeds =0;
}
-void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
+void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
{
//
//read seeds from the event
AliESDtrack *esd=event->GetTrack(i);
ULong_t status=esd->GetStatus();
if (!(status&AliESDtrack::kTPCin)) continue;
- AliTPCtrack t(*esd,fDebugStreamer);
+ AliTPCtrack t(*esd);
t.SetNumberOfClusters(0);
// AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
}
if (fDebug>3){
- Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
+ Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
}
delete seed;
}
}
-AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
+AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
{
//
//
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;
//
// General cuts - for splitted tracks and for curling tracks
//
- return;
const Float_t kMaxdPhi = 0.2; // maximal distance in phi
//
// Curling tracks cuts
}
}
//
- if (AliTPCReconstructor::StreamLevel()>0) {
+ if (AliTPCReconstructor::StreamLevel()>5) {
TTreeSRedirector &cstream = *fDebugStreamer;
cstream<<"Multi"<<
"iter="<<iter<<
}
delete [] helixes;
delete [] xm;
+ delete [] dz0;
+ delete [] dz1;
if (AliTPCReconstructor::StreamLevel()>1) {
AliInfo("Time for curling tracks removal DEBUGGING MC");
timer.Print();
}
-void AliTPCtrackerMI::FindSplitted(TObjArray * array){
+
+void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
//
// Find Splitted tracks and remove the one with worst quality
// Corresponding debug streamer to tune selections - "Splitted2"
//0. Sort tracks according quality
//1. Propagate the ext. param to reference radius
Int_t nseed = array->GetEntriesFast();
+ if (nseed<=0) return;
Float_t * quality = new Float_t[nseed];
Int_t * indexes = new Int_t[nseed];
for (Int_t i=0; i<nseed; i++) {
Float_t * points = pt->GetPoints();
if (points[3]<0.8) quality[i] =-1;
quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
+ //prefer high momenta tracks if overlaps
+ 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);
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()>2){
+ if( AliTPCReconstructor::StreamLevel()>1){
TTreeSRedirector &cstream = *fDebugStreamer;
Int_t n0=s1->GetNumberOfClusters();
Int_t n1=s2->GetNumberOfClusters();
// 4. Delete temporary array
//
delete [] params;
+ delete [] quality;
+ delete [] indexes;
+
}
//
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
//
track0->SetCircular(track0->GetCircular()+2);
}
}
- if (AliTPCReconstructor::StreamLevel()>1){
+ if (AliTPCReconstructor::StreamLevel()>2){
//
//debug stream to tune "fine" cuts
Int_t lab0=track0->GetLabel();
Int_t 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
kink->SetShapeFactor(shapesum/sum);
}
// esd->AddKink(kink);
+ //
+ // kink->SetMother(paramm);
+ //kink->SetDaughter(paramd);
+
+ Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
+ chi2P2*=chi2P2;
+ chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
+ Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
+ chi2P3*=chi2P3;
+ chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
+ //
+ if (AliTPCReconstructor::StreamLevel()>1) {
+ (*fDebugStreamer)<<"kinkLpt"<<
+ "chi2P2="<<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++;
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;
shared[kink0->GetIndex(0)]= kTRUE;
shared[kink0->GetIndex(1)]= kTRUE;
delete kinks->RemoveAt(indexes[ikink0]);
+ break;
}
else{
shared[kink1->GetIndex(0)]= kTRUE;
kinks->Delete();
delete kinks;
- printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
+ AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
timer.Print();
}
-void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *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, AliESDkink &knk)
+Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
{
//
// refit kink towards to the vertex
}
-Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
+Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
{
//
// check kink point for given track
delete seed1;
return 0;
}
+
// Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
kink.SetMother(param0[index]);
kink.SetDaughter(param1[index]);
kink.Update();
+
+ Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
+ chi2P2*=chi2P2;
+ chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
+ Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
+ chi2P3*=chi2P3;
+ chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
+ //
+ if (AliTPCReconstructor::StreamLevel()>1) {
+ (*fDebugStreamer)<<"kinkHpt"<<
+ "chi2P2="<<chi2P2<<
+ "chi2P3="<<chi2P3<<
+ "p0.="<<¶m0[index]<<
+ "p1.="<<¶m1[index]<<
+ "k.="<<&kink<<
+ "\n";
+ }
+ if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
+ delete seed0;
+ delete seed1;
+ return 0;
+ }
+
+
row0 = GetRowNumber(kink.GetR());
kink.SetTPCRow0(row0);
kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
daughter = param1[index];
daughter.SetLabel(kink.GetLabel(1));
param0[index].Reset(kTRUE);
- FollowProlongation(param0[index],0);
+ FollowProlongation(param0[index],0);
mother = param0[index];
mother.SetLabel(kink.GetLabel(0));
+ if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
+ mother=*seed;
+ }
delete seed0;
delete seed1;
//
return 0;
}
-Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
+Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
{
//
+
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);
+ if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
return 0;
//
}
RemoveUsed2(fSeeds,0.85,0.85,0);
if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
//FindCurling(fSeeds, fEvent,0);
- if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
+ if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
RemoveUsed2(fSeeds,0.5,0.4,20);
- FindSplitted(fSeeds); // find multi found tracks
- if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
+ FindSplitted(fSeeds, fEvent,0); // find multi found tracks
+ if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
// //
// // refit short tracks
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];
//
//find primaries
cuts[0]=0.0066;
- for (Int_t delta = 0; delta<18; delta+=6){
+ for (Int_t delta = 0; delta<18; delta+=gapPrim){
//
cuts[0]=0.0070;
cuts[1] = 1.5;
//find primaries
cuts[0]=0.0077;
- for (Int_t delta = 20; delta<120; delta+=10){
+ for (Int_t delta = 20; delta<120; delta+=gapPrim){
//
// seed high pt tracks
cuts[0]=0.0060;
SumTracks(seeds,arr);
SignClusters(seeds,fnumber,fdensity);
//
+ arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
+ SumTracks(seeds,arr);
+ SignClusters(seeds,fnumber,fdensity);
+ //
+ arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
+ SumTracks(seeds,arr);
+ SignClusters(seeds,fnumber,fdensity);
+ //
+ //
+ arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
+ SumTracks(seeds,arr);
+ SignClusters(seeds,fnumber,fdensity);
+ //
- for (Int_t delta = 3; delta<30; delta+=5){
+ for (Int_t delta = 9; delta<30; delta+=gapSec){
//
cuts[0] = 0.3;
cuts[1] = 1.5;
fdensity = 2.;
cuts[0]=0.0080;
- Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
// find secondaries
- for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
+ for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
//
cuts[0] = 0.3;
cuts[1] = 3.5;
-void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
+void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
{
//
// try to track in parralel
}
}
-void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
+void AliTPCtrackerMI::PrepareForBackProlongation(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)
}
-void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
+void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
{
//
//
}
-Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
+Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
{
//
// make back propagation
}
-Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
+Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
{
//
// make forward propagation
-Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
+Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
{
//
// make back propagation, in between row0 and row1
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);
if (TMath::Abs(c->GetLabel(1)) == lab ||
TMath::Abs(c->GetLabel(2)) == lab ) max++;
}
-
- if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
+ if (noc<=0) { lab=-1; return;}
+ if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
else {
Int_t tail=Int_t(0.10*noc);
max=0;
Int_t ind=0;
- for (i=1; i<=160&&ind<tail; i++) {
+ for (i=1; i<160&&ind<tail; i++) {
// AliTPCclusterMI *c=clusters[noc-i];
AliTPCclusterMI *c=clusters[i];
if (!c) continue;
//__________________________________________________________________________
-Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
+Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
//--------------------------------------------------------------------
//This function "cooks" a track label. If label<0, this track is fake.
//--------------------------------------------------------------------
}
}
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];
if (TMath::Abs(c->GetLabel(1)) == lab ||
TMath::Abs(c->GetLabel(2)) == lab ) max++;
}
-
- if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
+ if (noc<=0) { lab=-1; return -1;}
+ if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
else {
Int_t tail=Int_t(0.10*noc);
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];
for (Int_t i=0;i<15;i++) covar[i]=0;
// 0
covar[9] = param[3]*param[3];
Double_t facC = AliTracker::GetBz()*kB2C;
covar[14]= param[4]*param[4]*facC*facC;
+ //
+ covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
+ //
+ covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
+ covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
+ //
+ covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
+ covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
+ covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
+ //
+ covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
+ covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
+ covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
+ covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
+ //
seed->AddCovariance(covar);
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff