// The debug level - different procedure produce tree for numerical debugging
// To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
//
+
+//
+// Adding systematic errors to the covariance:
+//
+// The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
+// of the tracks (not to the clusters as they are dependent):
+// The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
+// The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
+// The default values are 0.
+//
+// The sytematic errors are added to the covariance matrix in following places:
+//
+// 1. During fisrt itteration - AliTPCtrackerMI::FillESD
+// 2. Second iteration -
+// 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
+// 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
+// 3. Third iteration -
+// 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
+// 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
+//
// There are several places in the code which can be numerically debuged
// This code is keeped in order to enable code development and to check the calibration implementtion
//
#include "AliTPCReconstructor.h"
#include "AliTPCpolyTrack.h"
#include "AliTPCreco.h"
-#include "AliTPCseed.h"
+#include "AliTPCseed.h"
+
+#include "AliTPCtrackerSector.h"
#include "AliTPCtrackerMI.h"
#include "TStopwatch.h"
#include "AliTPCReconstructor.h"
//
ClassImp(AliTPCtrackerMI)
-ClassImp(AliTPCtrackerRow)
+
class AliTPCFastMath {
Double_t rdistance2 = rdistancey2+rdistancez2;
//Int_t accept =0;
- if (AliTPCReconstructor::StreamLevel()>5) {
+ if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
Float_t rmsy2 = seed->GetCurrentSigmaY2();
Float_t rmsz2 = seed->GetCurrentSigmaZ2();
Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
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];
(*fDebugStreamer)<<"ErrParam"<<
"Cl.="<<cluster<<
"T.="<<¶m<<
+ "gcl.="<<&gcl<<
+ "gtr.="<<>r<<
"erry2="<<sy2<<
"errz2="<<sz2<<
"rmsy2="<<rmsy2<<
//---------------------------------------------------------------------
// The main TPC tracker constructor
//---------------------------------------------------------------------
- fInnerSec=new AliTPCSector[fkNIS];
- fOuterSec=new AliTPCSector[fkNOS];
+ fInnerSec=new AliTPCtrackerSector[fkNIS];
+ fOuterSec=new AliTPCtrackerSector[fkNOS];
Int_t i;
for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
if (!pt) continue;
pt->UpdatePoints();
+ AddCovariance(pt);
+ if (AliTPCReconstructor::StreamLevel()>1) {
+ (*fDebugStreamer)<<"Track0"<<
+ "Tr0.="<<pt<<
+ "\n";
+ }
// pt->PropagateTo(fParam->GetInnerRadiusLow());
if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
pt->PropagateTo(fParam->GetInnerRadiusLow());
return ;
}
+void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
+ //
+ // Filling cluster to the array - For visualization purposes
+ //
+ Int_t nrows=0;
+ nrows = fOuterSec->GetNRows();
+ for (Int_t sec = 0;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 AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
//
//
if (AliTPCReconstructor::StreamLevel()>1) {
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]<<
seed->PropagateTo(fParam->GetInnerRadiusLow());
seed->UpdatePoints();
+ AddCovariance(seed);
MakeBitmaps(seed);
AliESDtrack *esd=event->GetTrack(i);
seed->CookdEdx(0.02,0.6);
}
//FindKinks(fSeeds,event);
Info("RefitInward","Number of refitted tracks %d",ntracks);
- fEvent =0;
return 0;
}
if (!seed) continue;
if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
seed->UpdatePoints();
+ AddCovariance(seed);
AliESDtrack *esd=event->GetTrack(i);
seed->CookdEdx(0.02,0.6);
CookLabel(seed,0.1); //For comparison only
ntracks++;
Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
// This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
- if (AliTPCReconstructor::StreamLevel()>1) {
+ if (AliTPCReconstructor::StreamLevel()>1 && esd) {
(*fDebugStreamer)<<"Cback"<<
"Tr0.="<<seed<<
+ "esd.="<<esd<<
"EventNrInFile="<<eventnumber<<
- "\n"; // patch 28 fev 06
+ "\n";
}
}
}
// AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
seed->SetUniqueID(esd->GetID());
+ AddCovariance(seed); //add systematic ucertainty
for (Int_t ikink=0;ikink<3;ikink++) {
Int_t index = esd->GetKinkIndex(ikink);
seed->GetKinkIndexes()[ikink] = index;
}
-Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
-{
- //return pad row number for this x
- Double_t r;
- if (fN < 64){
- r=fRow[fN-1].GetX();
- if (x > r) return fN;
- r=fRow[0].GetX();
- if (x < r) return -1;
- return Int_t((x-r)/fPadPitchLength + 0.5);}
- else{
- r=fRow[fN-1].GetX();
- if (x > r) return fN;
- r=fRow[0].GetX();
- if (x < r) return -1;
- Double_t r1=fRow[64].GetX();
- if(x<r1){
- return Int_t((x-r)/f1PadPitchLength + 0.5);}
- else{
- return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
- }
-}
-
Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
{
//return pad row number for given x vector
return GetRowNumber(localx);
}
-//_________________________________________________________________________
-void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
- //-----------------------------------------------------------------------
- // Setup inner sector
- //-----------------------------------------------------------------------
- if (f==0) {
- fAlpha=par->GetInnerAngle();
- fAlphaShift=par->GetInnerAngleShift();
- fPadPitchWidth=par->GetInnerPadPitchWidth();
- fPadPitchLength=par->GetInnerPadPitchLength();
- fN=par->GetNRowLow();
- if(fRow)delete [] fRow;fRow = 0;
- fRow=new AliTPCtrackerRow[fN];
- for (Int_t i=0; i<fN; i++) {
- fRow[i].SetX(par->GetPadRowRadiiLow(i));
- fRow[i].SetDeadZone(1.5); //1.5 cm of dead zone
- }
- } else {
- fAlpha=par->GetOuterAngle();
- fAlphaShift=par->GetOuterAngleShift();
- fPadPitchWidth = par->GetOuterPadPitchWidth();
- fPadPitchLength = par->GetOuter1PadPitchLength();
- f1PadPitchLength = par->GetOuter1PadPitchLength();
- f2PadPitchLength = par->GetOuter2PadPitchLength();
- fN=par->GetNRowUp();
- if(fRow)delete [] fRow;fRow = 0;
- fRow=new AliTPCtrackerRow[fN];
- for (Int_t i=0; i<fN; i++) {
- fRow[i].SetX(par->GetPadRowRadiiUp(i));
- fRow[i].SetDeadZone(1.5); // 1.5 cm of dead zone
- }
- }
-}
-
-AliTPCtrackerRow::AliTPCtrackerRow():
- fDeadZone(0.),
- fClusters1(0),
- fN1(0),
- fClusters2(0),
- fN2(0),
- fN(0),
- fX(0.)
-{
- //
- // default constructor
- //
-}
-
-AliTPCtrackerRow::~AliTPCtrackerRow(){
- //
- for (Int_t i = 0; i < fN1; i++)
- fClusters1[i].~AliTPCclusterMI();
- delete [] fClusters1;
- for (Int_t i = 0; i < fN2; i++)
- fClusters2[i].~AliTPCclusterMI();
- delete [] fClusters2;
-}
-
-
-
-//_________________________________________________________________________
-void
-AliTPCtrackerRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
- //-----------------------------------------------------------------------
- // Insert a cluster into this pad row in accordence with its y-coordinate
- //-----------------------------------------------------------------------
- if (fN==kMaxClusterPerRow) {
- //AliInfo("AliTPCtrackerRow::InsertCluster(): Too many clusters");
- return;
- }
- if (fN>=fN1+fN2) {
- //AliInfo("AliTPCtrackerRow::InsertCluster(): Too many clusters !");
- }
-
- if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
- Int_t i=Find(c->GetZ());
- memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
- memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
- fIndex[i]=index; fClusters[i]=c; fN++;
-}
-
-void AliTPCtrackerRow::ResetClusters() {
- //
- // reset clusters
- // MvL: Need to call destructors for AliTPCclusterMI, to delete fInfo
- for (Int_t i = 0; i < fN1; i++)
- fClusters1[i].~AliTPCclusterMI();
- delete [] fClusters1; fClusters1=0;
- for (Int_t i = 0; i < fN2; i++)
- fClusters2[i].~AliTPCclusterMI();
- delete [] fClusters2; fClusters2=0;
-
- fN = 0;
- fN1 = 0;
- fN2 = 0;
- //delete[] fClusterArray;
-
- //fClusterArray=0;
-}
-
-
-//___________________________________________________________________
-Int_t AliTPCtrackerRow::Find(Double_t z) const {
- //-----------------------------------------------------------------------
- // Return the index of the nearest cluster
- //-----------------------------------------------------------------------
- if (fN==0) return 0;
- if (z <= fClusters[0]->GetZ()) return 0;
- if (z > fClusters[fN-1]->GetZ()) return fN;
- Int_t b=0, e=fN-1, m=(b+e)/2;
- for (; b<e; m=(b+e)/2) {
- if (z > fClusters[m]->GetZ()) b=m+1;
- else e=m;
- }
- return m;
-}
-
-
-
-//___________________________________________________________________
-AliTPCclusterMI * AliTPCtrackerRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
- //-----------------------------------------------------------------------
- // Return the index of the nearest cluster in z y
- //-----------------------------------------------------------------------
- Float_t maxdistance = roady*roady + roadz*roadz;
-
- AliTPCclusterMI *cl =0;
- for (Int_t i=Find(z-roadz); i<fN; i++) {
- AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
- if (c->GetZ() > z+roadz) break;
- if ( (c->GetY()-y) > roady ) continue;
- Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
- if (maxdistance>distance) {
- maxdistance = distance;
- cl=c;
- }
- }
- return cl;
-}
-
-AliTPCclusterMI * AliTPCtrackerRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
-{
- //-----------------------------------------------------------------------
- // Return the index of the nearest cluster in z y
- //-----------------------------------------------------------------------
- Float_t maxdistance = roady*roady + roadz*roadz;
- AliTPCclusterMI *cl =0;
-
- //PH Check boundaries. 510 is the size of fFastCluster
- Int_t iz1 = Int_t(z-roadz+254.5);
- if (iz1<0 || iz1>=510) return cl;
- iz1 = TMath::Max(GetFastCluster(iz1)-1,0);
- Int_t iz2 = Int_t(z+roadz+255.5);
- if (iz2<0 || iz2>=510) return cl;
- iz2 = TMath::Min(GetFastCluster(iz2)+1,fN);
- Bool_t skipUsed = !(AliTPCReconstructor::GetRecoParam()->GetClusterSharing());
- //FindNearest3(y,z,roady,roadz,index);
- // for (Int_t i=Find(z-roadz); i<fN; i++) {
- for (Int_t i=iz1; i<iz2; i++) {
- AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
- if (c->GetZ() > z+roadz) break;
- if ( c->GetY()-y > roady ) continue;
- if ( y-c->GetY() > roady ) continue;
- if (skipUsed && c->IsUsed(11)) continue;
- Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
- if (maxdistance>distance) {
- maxdistance = distance;
- cl=c;
- index =i;
- //roady = TMath::Sqrt(maxdistance);
- }
- }
- return cl;
-}
-
-
-void AliTPCtrackerRow::SetFastCluster(Int_t i, Short_t cl){
- //
- // Set cluster info for fast navigation
- //
- if (i>510|| i<0){
- }else{
- fFastCluster[i]=cl;
- }
-}
-
-
void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
}
}
}
+
+
+
+void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
+ //
+ // Adding systematic error
+ // !!!! the systematic error for element 4 is in 1/cm not in pt
+
+ const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
+ Double_t covar[15];
+ for (Int_t i=0;i<15;i++) covar[i]=0;
+ // 0
+ // 1 2
+ // 3 4 5
+ // 6 7 8 9
+ // 10 11 12 13 14
+ covar[0] = param[0]*param[0];
+ covar[2] = param[1]*param[1];
+ covar[5] = param[2]*param[2];
+ covar[9] = param[3]*param[3];
+ Double_t facC = AliTracker::GetBz()*kB2C;
+ covar[14]= param[4]*param[4]*facC*facC;
+ seed->AddCovariance(covar);
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff