#include "AliTPCTracklet.h"
#include "TH1D.h"
#include "TH2F.h"
+#include "THnSparse.h"
#include "TVectorD.h"
#include "TTreeStream.h"
#include "TFile.h"
#include "AliTPCseed.h"
#include "AliTracker.h"
#include "TClonesArray.h"
-#include "AliExternalComparison.h"
#include "AliLog.h"
#include "TFile.h"
#include "TProfile.h"
#include "TCanvas.h"
-#include "TLegend.h"
-
+#include "TDatabasePDG.h"
#include "TTreeStream.h"
-#include <iostream>
+#include "Riostream.h"
#include <sstream>
using namespace std;
fMatrixArray9(72*72),
fMatrixArray6(72*72),
fCombinedMatrixArray6(72),
- fCompTracklet(0), // tracklet comparison
fNoField(kFALSE),
fXIO(0),
fXmiddle(0),
fXquadrant = roc->GetPadRowRadii(36,53);
fXmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
+ fClusterDelta[0]=0; // cluster residuals - Y
+ fClusterDelta[1]=0; // cluster residuals - Z
+
+
+ fTrackletDelta[0]=0; // tracklet residuals
+ fTrackletDelta[1]=0; // tracklet residuals
+ fTrackletDelta[2]=0; // tracklet residuals
+ fTrackletDelta[3]=0; // tracklet residuals
}
AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
fMatrixArray9(72*72),
fMatrixArray6(72*72),
fCombinedMatrixArray6(72),
- fCompTracklet(0), // tracklet comparison
fNoField(kFALSE),
fXIO(0),
fXmiddle(0),
fXquadrant = roc->GetPadRowRadii(36,53);
fXmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
fXIO = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
+ fClusterDelta[0]=0; // cluster residuals
+ fClusterDelta[1]=0; // cluster residuals
-
+ fTrackletDelta[0]=0; // tracklet residuals
+ fTrackletDelta[1]=0; // tracklet residuals
+ fTrackletDelta[2]=0; // tracklet residuals
+ fTrackletDelta[3]=0; // tracklet residuals
}
fMatrixArray9(align.fMatrixArray9),
fMatrixArray6(align.fMatrixArray6),
fCombinedMatrixArray6(align.fCombinedMatrixArray6),
- fCompTracklet(align.fCompTracklet), // tracklet comparison
fNoField(align.fNoField),
fXIO(align.fXIO),
fXmiddle(align.fXmiddle),
fSectorParamC = (TMatrixD*)align.fSectorParamA->Clone();
fSectorParamC = (TMatrixD*)align.fSectorCovarA->Clone();
}
+ fClusterDelta[0]=0; // cluster residuals
+ fClusterDelta[1]=0; // cluster residuals
+
+ fTrackletDelta[0]=0; // tracklet residuals
+ fTrackletDelta[1]=0; // tracklet residuals
+ fTrackletDelta[2]=0; // tracklet residuals
+ fTrackletDelta[3]=0; // tracklet residuals
}
fMatrixArray9.Delete(); // array of transnformtation matrix
fMatrixArray6.Delete(); // array of transnformtation matrix
- if (fCompTracklet) delete fCompTracklet;
fArraySectorIntParam.SetOwner(kTRUE); // array of sector alignment parameters
fArraySectorIntCovar.SetOwner(kTRUE); // array of sector alignment covariances
fArraySectorIntParam.Delete(); // array of sector alignment parameters
fArraySectorIntCovar.Delete(); // array of sector alignment covariances
+ for (Int_t i=0; i<2; i++){
+ delete fClusterDelta[i]; // cluster residuals
+ }
+
+ for (Int_t i=0; i<4; i++){
+ delete fTrackletDelta[i]; // tracklet residuals
+ }
+
}
//
// Process pairs of cosmic tracks
//
+ if (!fClusterDelta[0]) MakeResidualHistos();
+ if (!fTrackletDelta[0]) MakeResidualHistosTracklet();
+ //
+ fCurrentEvent=event;
ExportTrackPoints(event); // export track points for external calibration
- const Int_t kMaxTracks =50;
+ const Int_t kMaxTracks =6;
const Int_t kminCl = 40;
- AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
- if (!ESDfriend) return;
+ AliESDfriend *
eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+ if (!eESDfriend) return;
Int_t ntracks=event->GetNumberOfTracks();
Float_t dca0[2];
Float_t dca1[2];
//
//
+ // process seeds
//
+ for (Int_t i0=0;i0<ntracks;++i0) {
+ AliESDtrack *track0 = event->GetTrack(i0);
+ AliESDfriendTrack *friendTrack = 0;
+ TObject *calibObject=0;
+ AliTPCseed *seed0 = 0;
+ //
+ friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
+ if (!friendTrack) continue;
+ for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
+ if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
+ }
+ if (!seed0) continue;
+ fCurrentTrack=track0;
+ fCurrentFriendTrack=friendTrack;
+ fCurrentSeed=seed0;
+ fCurrentEvent=event;
+ ProcessSeed(seed0);
+ }
+ //
+ // process cosmic pairs
//
if (ntracks>kMaxTracks) return;
//
// fast cuts on dca and theta
// if (TMath::Abs(dca1[0]+dca0[0])>15) continue;
// if (TMath::Abs(dca1[1]-dca0[1])>15) continue;
- // if (TMath::Abs(track0->GetParameter()[3]+track1->GetParameter()[3])>0.1) continue;
+ if (TMath::Abs(track0->GetParameter()[3]+track1->GetParameter()[3])>0.1) continue;
//
AliESDfriendTrack *friendTrack = 0;
TObject *calibObject=0;
AliTPCseed *seed0 = 0,*seed1=0;
//
- friendTrack = (AliESDfriendTrack *)ESDfriend->GetTrack(i0);;
+ friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i0);;
+ if (!friendTrack) continue;
for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
}
- friendTrack = (AliESDfriendTrack *)ESDfriend->GetTrack(i1);;
+ friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(i1);;
+ if (!friendTrack) continue;
for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
}
-
if (!seed0) continue;
+ //
+ //
if (!seed1) continue;
Int_t nclsectors0[72], nclsectors1[72];
for (Int_t isec=0;isec<72;isec++){
// Export track points for alignment - calibration
// export space points for pairs of tracks if possible
//
- AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
- if (!ESDfriend) return;
+ AliESDfriend *
eESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+ if (!eESDfriend) return;
Int_t ntracks=event->GetNumberOfTracks();
+ Int_t kMaxTracks=4; // maximal number of tracks for cosmic pairs
+ Int_t kMinVertexTracks=5; // maximal number of tracks for vertex mesurement
+
//cuts
const Int_t kminCl = 60;
const Int_t kminClSum = 120;
- // const Double_t kDistY = 5;
+ //const Double_t kDistY = 5;
// const Double_t kDistZ = 40;
const Double_t kDistTh = 0.05;
+ const Double_t kDistThS = 0.002;
const Double_t kDist1Pt = 0.1;
+ const Double_t kMaxD0 =3; // max distance to the primary vertex
+ const Double_t kMaxD1 =5; // max distance to the primary vertex
+ const AliESDVertex *tpcVertex = 0;
+ // get the primary vertex TPC
+ if (ntracks>kMinVertexTracks) {
+ tpcVertex = event->GetPrimaryVertexSPD();
+ if (tpcVertex->GetNContributors()<kMinVertexTracks) tpcVertex=0;
+ }
//
Float_t dca0[2];
- Float_t dca1[2];
+ // Float_t dca1[2];
Int_t index0=0,index1=0;
//
for (Int_t i0=0;i0<ntracks;++i0) {
if (!track0) continue;
if ((track0->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
if (track0->GetOuterParam()==0) continue;
+ if (track0->GetInnerParam()==0) continue;
+ if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt()-track0->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
+ if (TMath::Abs(track0->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
+ if (TMath::Abs(track0->GetInnerParam()->GetTgl()-track0->GetOuterParam()->GetTgl())>kDistThS) continue;
AliESDtrack *track1P = 0;
if (track0->GetTPCNcls()<kminCl) continue;
track0->GetImpactParameters(dca0[0],dca0[1]);
index0=i0;
index1=-1;
//
- for (Int_t i1=0;i1<ntracks;++i1) {
+ if (ntracks<kMaxTracks) for (Int_t i1=i0+1;i1<ntracks;++i1) {
if (i0==i1) continue;
AliESDtrack *track1 = event->GetTrack(i1);
if (!track1) continue;
if ((track1->GetStatus() & AliESDtrack::kTPCrefit)==0) continue;
if (track1->GetOuterParam()==0) continue;
+ if (track1->GetInnerParam()==0) continue;
if (track1->GetTPCNcls()<kminCl) continue;
- track1->GetImpactParameters(dca1[0],dca1[1]);
+ if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt()-track1->GetOuterParam()->GetSigned1Pt())>kDist1Pt) continue;
+ if (TMath::Abs(track1->GetInnerParam()->GetTgl()-track1->GetOuterParam()->GetTgl())>kDistThS) continue;
+ if (TMath::Abs(track1->GetInnerParam()->GetSigned1Pt())>kDist1Pt) continue;
+ //track1->GetImpactParameters(dca1[0],dca1[1]);
//if (TMath::Abs(dca1[0]-dca0[0])>kDistY) continue;
//if (TMath::Abs(dca1[1]-dca0[1])>kDistZ) continue;
if (TMath::Abs(track0->GetTgl()+track1->GetTgl())>kDistTh) continue;
TObject *calibObject=0;
AliTPCseed *seed0 = 0,*seed1=0;
//
- friendTrack = (AliESDfriendTrack *)ESDfriend->GetTrack(index0);;
+ friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index0);;
+ if (!friendTrack) continue;
for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
}
if (index1>0){
- friendTrack = (AliESDfriendTrack *)ESDfriend->GetTrack(index1);;
+ friendTrack = (AliESDfriendTrack *)eESDfriend->GetTrack(index1);;
+ if (!friendTrack) continue;
for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
}
if (npoints<kminClSum) continue;
Int_t cpoint=0;
AliTrackPointArray array(npoints);
+ if (tpcVertex){
+ Double_t dxyz[3]={0,0,0};
+ Double_t dc[6]={0,0,0};
+ tpcVertex->GetXYZ(dxyz);
+ tpcVertex->GetCovarianceMatrix(dc);
+ Float_t xyz[3]={dxyz[0],dxyz[1],dxyz[2]};
+ Float_t cov[6]={dc[0]+1,dc[1],dc[2]+1,dc[3],dc[4],dc[5]+100.};
+ AliTrackPoint point(xyz,cov,73); // add point to not existing volume
+ Float_t dtpc[2],dcov[3];
+ track0->GetImpactParametersTPC(dtpc,dcov);
+ if (TMath::Abs(dtpc[0])>kMaxD0) continue;
+ if (TMath::Abs(dtpc[1])>kMaxD1) continue;
+ }
+
if (seed0) for (Int_t icl = 0; icl<160; icl++){
AliTPCclusterMI *cluster=seed0->GetClusterPointer(icl);
if (!cluster) continue;
cluster->GetGlobalXYZ(xyz);
cluster->GetGlobalCov(cov);
AliTrackPoint point(xyz,cov,cluster->GetDetector());
- array.AddPoint(npoints, &point);
+ array.AddPoint(npoints, &point);
if (cpoint>=npoints) continue; //shoul not happen
array.AddPoint(cpoint, &point);
cpoint++;
//
TTreeSRedirector *cstream = GetDebugStreamer();
if (cstream){
+ Bool_t isVertex=(tpcVertex)? kTRUE:kFALSE;
+ Double_t tof0=track0->GetTOFsignal();
+ Double_t tof1=(track1P) ? track1P->GetTOFsignal(): 0;
static AliExternalTrackParam dummy;
AliExternalTrackParam *p0In = &dummy;
AliExternalTrackParam *p1In = &dummy;
AliExternalTrackParam *p0Out = &dummy;
AliExternalTrackParam *p1Out = &dummy;
+ AliESDVertex vdummy;
+ AliESDVertex *pvertex= (tpcVertex)? (AliESDVertex *)tpcVertex: &vdummy;
if (track0) {
p0In= new AliExternalTrackParam(*track0);
p0Out=new AliExternalTrackParam(*(track0->GetOuterParam()));
"trigger="<<fTrigger<< // trigger
"triggerClass="<<&fTriggerClass<< // trigger
"mag="<<fMagF<< // magnetic field
+ "pvertex.="<<pvertex<< // vertex
//
+ "isVertex="<<isVertex<< // flag is with prim vertex
+ "tof0="<<tof0<< // tof signal 0
+ "tof1="<<tof1<< // tof signal 1
+ "seed0.="<<seed0<< // track info
"ntracks="<<ntracks<< // number of tracks
"ncont="<<ncont<< // number of contributors
"p0In.="<<p0In<< // track parameters0
-void AliTPCcalibAlign::Process(AliTPCseed *seed) {
+void AliTPCcalibAlign::ProcessSeed(AliTPCseed *seed) {
//
//
//
// make a kalman tracklets out of seed
//
+ UpdateClusterDeltaField(seed);
TObjArray tracklets=
AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
kFALSE,20,4);
AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[i2]);
AliExternalTrackParam *common1=0,*common2=0;
if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2)){
- ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
+ ProcessTracklets(*common1,*common2,seed, t1->GetSector(),t2->GetSector());
UpdateAlignSector(seed,t1->GetSector());
}
delete common1;
//
// Process function to fill fitters
//
+ if (!seed) return;
Double_t t1[10],t2[10];
Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
Double_t &x2=t2[0], &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
//
Int_t accept = AcceptTracklet(tp1,tp2);
Int_t acceptLinear = AcceptTracklet(parLine1,parLine2);
-
- if (fStreamLevel>1 && seed){
+
+
+ if (fStreamLevel>1){
TTreeSRedirector *cstream = GetDebugStreamer();
if (cstream){
static TVectorD vec1(5);
// use Linear fit
//
if (nl1>10 && nl2>10 &&(acceptLinear==0)){
- if (seed) ProcessDiff(tp1,tp2, seed,s1,s2);
+ ProcessDiff(tp1,tp2, seed,s1,s2);
if (TMath::Abs(parLine1[2])<0.8 &&TMath::Abs(parLine1[2])<0.8 ){ //angular cut
FillHisto(parLine1,parLine2,s1,s2);
ProcessAlign(parLine1,parLine2,s1,s2);
- UpdateKalman(s1,s2,par1, cov1, par2, cov2);
+ FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
+ FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
+ //UpdateKalman(s1,s2,par1, cov1, par2, cov2); - OBSOLETE to be removed - 50 % of time here
}
}
}
//
// use Kalman if mag field
//
- if (seed) ProcessDiff(tp1,tp2, seed,s1,s2);
+ ProcessDiff(tp1,tp2, seed,s1,s2);
+ FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
+ FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
FillHisto(t1,t2,s1,s2);
ProcessAlign(t1,t2,s1,s2);
}
//
// Do intersector alignment
//
- Process12(t1,t2,GetOrMakeFitter12(s1,s2));
- Process9(t1,t2,GetOrMakeFitter9(s1,s2));
+ //Process12(t1,t2,GetOrMakeFitter12(s1,s2));
+ //Process9(t1,t2,GetOrMakeFitter9(s1,s2));
Process6(t1,t2,GetOrMakeFitter6(s1,s2));
++fPoints[GetIndex(s1,s2)];
}
-void AliTPCcalibAlign::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){
- //
- // Process the debug streamer tree
- // Possible to modify selection criteria
- // Used with entry list
- //
- TTreeSRedirector * cstream = new TTreeSRedirector("aligndump.root");
-
- AliTPCcalibAlign *align = this;
- //
- TVectorD * vec1 = 0;
- TVectorD * vec2 = 0;
- AliExternalTrackParam * tp1 = 0;
- AliExternalTrackParam * tp2 = 0;
- Int_t s1 = 0;
- Int_t s2 = 0;
- Int_t npoints =0;
- {
- Int_t entries=chainTracklet->GetEntries();
- for (Int_t i=0; i< entries; i++){
- chainTracklet->GetBranch("tp1.")->SetAddress(&tp1);
- chainTracklet->GetBranch("tp2.")->SetAddress(&tp2);
- chainTracklet->GetBranch("v1.")->SetAddress(&vec1);
- chainTracklet->GetBranch("v2.")->SetAddress(&vec2);
- chainTracklet->GetBranch("s1")->SetAddress(&s1);
- chainTracklet->GetBranch("s2")->SetAddress(&s2);
- chainTracklet->GetEntry(i);
- if (!vec1) continue;
- if (!vec2) continue;
- if (!tp1) continue;
- if (!tp2) continue;
- if (!vec1->GetMatrixArray()) continue;
- if (!vec2->GetMatrixArray()) continue;
- // make a local copy
- AliExternalTrackParam par1(*tp1);
- AliExternalTrackParam par2(*tp2);
- TVectorD svec1(*vec1);
- TVectorD svec2(*vec2);
- //
- if (s1==s2) continue;
- if (i%100==0) printf("%d\t%d\t%d\t%d\t\n",i, npoints,s1,s2);
- AliExternalTrackParam cpar1(par1);
- AliExternalTrackParam cpar2(par2);
- Constrain1Pt(cpar1,par2,fNoField);
- Constrain1Pt(cpar2,par1,fNoField);
- Bool_t acceptComp = kFALSE;
- if (comp) acceptComp=comp->AcceptPair(&par1,&par2);
- if (comp) acceptComp&=comp->AcceptPair(&cpar1,&cpar2);
- //
- Int_t reject = align->AcceptTracklet(par1,par2);
- Int_t rejectC =align->AcceptTracklet(cpar1,cpar2);
-
- if (1||fStreamLevel>0){
- (*cstream)<<"Tracklet"<<
- "s1="<<s1<<
- "s2="<<s2<<
- "reject="<<reject<<
- "rejectC="<<rejectC<<
- "acceptComp="<<acceptComp<<
- "tp1.="<<&par1<<
- "tp2.="<<&par2<<
- "ctp1.="<<&cpar1<<
- "ctp2.="<<&cpar2<<
- "v1.="<<&svec1<<
- "v2.="<<&svec2<<
- "\n";
- }
- //
- if (fNoField){
- //
- //
- }
- if (acceptComp) comp->Process(&cpar1,&cpar2);
- //
- if (reject>0 || rejectC>0) continue;
- npoints++;
- align->ProcessTracklets(cpar1,cpar2,0,s1,s2);
- align->ProcessTracklets(cpar2,cpar1,0,s2,s1);
- }
- }
- delete cstream;
-}
Int_t AliTPCcalibAlign::AcceptTracklet(const AliExternalTrackParam &p1,
- const AliExternalTrackParam &p2){
+ const AliExternalTrackParam &p2) const
+{
//
// Accept pair of tracklets?
}
-Int_t AliTPCcalibAlign::AcceptTracklet(const Double_t *t1, const Double_t *t2){
+Int_t AliTPCcalibAlign::AcceptTracklet(const Double_t *t1, const Double_t *t2) const
+{
//
// accept tracklet -
// dist cut + 6 sigma cut
void AliTPCcalibAlign::Process12(const Double_t *t1,
const Double_t *t2,
- TLinearFitter *fitter) {
+ TLinearFitter *fitter) const
+{
// x2 = a00*x1 + a01*y1 + a02*z1 + a03
// y2 = a10*x1 + a11*y1 + a12*z1 + a13
// z2 = a20*x1 + a21*y1 + a22*z1 + a23
fitter->AddPoint(p,value,sdzdx);
}
-void AliTPCcalibAlign::Process9(Double_t *t1,
- Double_t *t2,
- TLinearFitter *fitter) {
+void AliTPCcalibAlign::Process9(const Double_t * const t1,
+ const Double_t * const t2,
+ TLinearFitter *fitter) const
+{
// x2 = a00*x1 + a01*y1 + a02*z1 + a03
// y2 = a10*x1 + a11*y1 + a12*z1 + a13
// z2 = a20*x1 + a21*y1 + a22*z1 + a23
// a20 a21 a21 a23 p[4] p[5] 1 p[8]
- Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
- Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
+ const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
+ const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
//
Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
Double_t sdydx = TMath::Sqrt(t1[7]*t1[7]+t2[7]*t2[7]);
fitter->AddPoint(p,value,sdzdx);
}
-void AliTPCcalibAlign::Process6(Double_t *t1,
- Double_t *t2,
- TLinearFitter *fitter) {
+void AliTPCcalibAlign::Process6(const Double_t *const t1,
+ const Double_t *const t2,
+ TLinearFitter *fitter) const
+{
// x2 = 1 *x1 +-a01*y1 + 0 +a03
// y2 = a01*x1 + 1 *y1 + 0 +a13
// z2 = a20*x1 + a21*y1 + 1 *z1 +a23
// a10 a11 a12 a13 ==> p[0] 1 0 p[4]
// a20 a21 a21 a23 p[1] p[2] 1 p[5]
- Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
- Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
+ const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[3], &dydx1=t1[2], &dzdx1=t1[4];
+ const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[3], &dydx2=t2[2], &dzdx2=t2[4];
//
Double_t sy = TMath::Sqrt(t1[6]*t1[6]+t2[6]*t2[6]);
fitter->StoreData(kFALSE);
fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
counter12++;
- if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter12<<endl;
+ // if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter12<<endl;
return fitter;
}
fitter->StoreData(kFALSE);
fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
counter9++;
- if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter9<<endl;
+ // if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter9<<endl;
return fitter;
}
fitter->StoreData(kFALSE);
fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
counter6++;
- if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<" : "<<counter6<<endl;
+ // if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<" : "<<counter6<<endl;
return fitter;
}
}
}
+void AliTPCcalibAlign::MakeResidualHistos(){
+ //
+ // Make cluster residual histograms
+ //
+ Double_t xminTrack[9], xmaxTrack[9];
+ Int_t binsTrack[9];
+ TString axisName[9],axisTitle[9];
+ //
+ // 0 - delta of interest
+ // 1 - global phi in sector number as float
+ // 2 - local x
+ // 3 - local ky
+ // 4 - local kz
+ //
+ axisName[0]="delta"; axisTitle[0]="#Delta (cm)";
+ if (TMath::Abs(AliTracker::GetBz())<0.01){
+ binsTrack[0]=60; xminTrack[0]=-1.2; xmaxTrack[0]=1.2;
+ }else{
+ binsTrack[0]=60; xminTrack[0]=-0.6; xmaxTrack[0]=0.6;
+ }
+ //
+ axisName[1]="sector"; axisTitle[1]="Sector Number";
+ binsTrack[1]=180; xminTrack[1]=0; xmaxTrack[1]=18;
+ //
+ axisName[2]="R"; axisTitle[2]="r (cm)";
+ binsTrack[2]=53; xminTrack[2]=85.; xmaxTrack[2]=245.;
+ //
+ //
+ axisName[3]="kZ"; axisTitle[3]="dz/dx";
+ binsTrack[3]=36; xminTrack[3]=-1.8; xmaxTrack[3]=1.8;
+ //
+ fClusterDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
+ fClusterDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
+ //
+ //
+ //
+ for (Int_t ivar=0;ivar<2;ivar++){
+ for (Int_t ivar2=0;ivar2<4;ivar2++){
+ fClusterDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
+ fClusterDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
+ }
+ }
+
+}
+
+
+void AliTPCcalibAlign::MakeResidualHistosTracklet(){
+ //
+ // Make tracklet residual histograms
+ //
+ Double_t xminTrack[9], xmaxTrack[9];
+ Int_t binsTrack[9];
+ TString axisName[9],axisTitle[9];
+ //
+ // 0 - delta of interest
+ // 1 - global phi in sector number as float
+ // 2 - local x
+ // 3 - local ky
+ // 4 - local kz
+ // 5 - sector 1
+ // 6 - sector 0
+ // 7 - z position 0
+
+ axisName[0]="delta"; axisTitle[0]="#Delta (cm)";
+ binsTrack[0]=60; xminTrack[0]=-0.6; xmaxTrack[0]=0.6;
+ //
+ axisName[1]="phi"; axisTitle[1]="#phi";
+ binsTrack[1]=180; xminTrack[1]=-TMath::Pi(); xmaxTrack[1]=TMath::Pi();
+ //
+ axisName[2]="localX"; axisTitle[2]="x (cm)";
+ binsTrack[2]=10; xminTrack[2]=120.; xmaxTrack[2]=200.;
+ //
+ axisName[3]="kY"; axisTitle[3]="dy/dx";
+ binsTrack[3]=10; xminTrack[3]=-0.5; xmaxTrack[3]=0.5;
+ //
+ axisName[4]="kZ"; axisTitle[4]="dz/dx";
+ binsTrack[4]=22; xminTrack[4]=-1.1; xmaxTrack[4]=1.1;
+ //
+ axisName[5]="is1"; axisTitle[5]="is1";
+ binsTrack[5]=72; xminTrack[5]=0; xmaxTrack[5]=72;
+ //
+ axisName[6]="is0"; axisTitle[6]="is0";
+ binsTrack[6]=72; xminTrack[6]=0; xmaxTrack[6]=72;
+ //
+ axisName[7]="z"; axisTitle[7]="z(cm)";
+ binsTrack[7]=12; xminTrack[7]=-240; xmaxTrack[7]=240;
+ //
+ axisName[8]="IsPrimary"; axisTitle[8]="Is Primary";
+ binsTrack[8]=2; xminTrack[8]=-0.1; xmaxTrack[8]=1.1;
+
+ //
+ xminTrack[0]=-0.25; xmaxTrack[0]=0.25;
+ fTrackletDelta[0] = new THnSparseF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0]=-0.5; xmaxTrack[0]=0.5;
+ fTrackletDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0]=-0.005; xmaxTrack[0]=0.005;
+ fTrackletDelta[2] = new THnSparseF("#Delta_{kY}","#Delta_{kY}", 9, binsTrack,xminTrack, xmaxTrack);
+ xminTrack[0]=-0.008; xmaxTrack[0]=0.008;
+ fTrackletDelta[3] = new THnSparseF("#Delta_{kZ}","#Delta_{kZ}", 9, binsTrack,xminTrack, xmaxTrack);
+ //
+ //
+ //
+ for (Int_t ivar=0;ivar<4;ivar++){
+ for (Int_t ivar2=0;ivar2<9;ivar2++){
+ fTrackletDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
+ fTrackletDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
+ }
+ }
+
+}
+
+
+
void AliTPCcalibAlign::FillHisto(const Double_t *t1,
const Double_t *t2,
Int_t s1,Int_t s2) {
}
+void AliTPCcalibAlign::FillHisto(AliExternalTrackParam *tp1,
+ AliExternalTrackParam *tp2,
+ Int_t s1,Int_t s2) {
+ //
+ // Fill residual histograms
+ // Track2-Track1
+ if (s2<s1) return;//
+ const Double_t kEpsilon=0.001;
+ Double_t x[9]={0,0,0,0,0,0,0,0,0};
+ AliExternalTrackParam p1(*tp1);
+ AliExternalTrackParam p2(*tp2);
+ if (s1%18==s2%18) {
+ // inner outer - match at the IROC-OROC boundary
+ if (!p1.PropagateTo(fXIO, AliTrackerBase::GetBz())) return;
+ }
+ if (!p2.Rotate(p1.GetAlpha())) return;
+ if (!p2.PropagateTo(p1.GetX(),AliTrackerBase::GetBz())) return;
+ if (TMath::Abs(p1.GetX()-p2.GetX())>kEpsilon) return;
+ Double_t xyz[3];
+ p1.GetXYZ(xyz);
+ x[1]=TMath::ATan2(xyz[1],xyz[0]);
+ x[2]=p1.GetX();
+ x[3]=0.5*(p1.GetSnp()+p2.GetSnp()); // mean snp
+ x[4]=0.5*(p1.GetTgl()+p2.GetTgl()); // mean tgl
+ x[5]=s2;
+ x[6]=s1;
+ x[7]=0.5*(p1.GetZ()+p2.GetZ());
+ // is primary ?
+ Int_t isPrimary = (TMath::Abs(p1.GetTgl()-p1.GetZ()/p1.GetX())<0.1) ? 1:0;
+ x[8]= isPrimary;
+ //
+ x[0]=p2.GetY()-p1.GetY();
+ fTrackletDelta[0]->Fill(x);
+ x[0]=p2.GetZ()-p1.GetZ();
+ fTrackletDelta[1]->Fill(x);
+ x[0]=p2.GetSnp()-p1.GetSnp();
+ fTrackletDelta[2]->Fill(x);
+ x[0]=p2.GetTgl()-p1.GetTgl();
+ fTrackletDelta[3]->Fill(x);
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ if (cstream){
+ (*cstream)<<"trackletMatch"<<
+ "tp1.="<<tp1<< // input tracklet
+ "tp2.="<<tp2<<
+ "p1.="<<&p1<< // tracklet in the ref frame
+ "p2.="<<&p2<<
+ "s1="<<s1<<
+ "s2="<<s2<<
+ "\n";
+ }
+
+}
+
+
TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
{
}
TGraphErrors *gr = new TGraphErrors(npoints,xsec,ysec,0,0);
Char_t name[1000];
- sprintf(name,"Mat[%d,%d] Type=%d",i0,i1,type);
+ snprintf(name,100,"Mat[%d,%d] Type=%d",i0,i1,type);
gr->SetName(name);
return gr;
}
//
//
- TMatrixD * kpar = fSectorParamA;
- TMatrixD * kcov = fSectorCovarA;
- if (s1%36>=18){
- kpar = fSectorParamC;
- kcov = fSectorCovarC;
- }
- for (Int_t ipar=0;ipar<6;ipar++){
- Int_t isec1 = s1%18;
- Int_t isec2 = s2%18;
- if (s1>35) isec1+=18;
- if (s2>35) isec2+=18;
- param6s1(ipar)=(*kpar)(6*isec1+ipar,0);
- param6s2(ipar)=(*kpar)(6*isec2+ipar,0);
+ if (fSectorParamA){
+ TMatrixD * kpar = fSectorParamA;
+ TMatrixD * kcov = fSectorCovarA;
+ if (s1%36>=18){
+ kpar = fSectorParamC;
+ kcov = fSectorCovarC;
+ }
+ for (Int_t ipar=0;ipar<6;ipar++){
+ Int_t isec1 = s1%18;
+ Int_t isec2 = s2%18;
+ if (s1>35) isec1+=18;
+ if (s2>35) isec2+=18;
+ param6s1(ipar)=(*kpar)(6*isec1+ipar,0);
+ param6s2(ipar)=(*kpar)(6*isec2+ipar,0);
+ }
}
-
-
+
Double_t dy=0, dz=0, dphi=0,dtheta=0;
Double_t sy=0, sz=0, sphi=0,stheta=0;
Double_t ny=0, nz=0, nphi=0,ntheta=0;
Double_t chi2v12=0, chi2v9=0, chi2v6=0;
- Int_t npoints=0;
- TLinearFitter * fitter = 0;
+ // Int_t npoints=0;
+ // TLinearFitter * fitter = 0;
if (fPoints[GetIndex(s1,s2)]>minPoints){
//
//
//
- fitter = GetFitter12(s1,s2);
- npoints = fitter->GetNpoints();
- chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
+// fitter = GetFitter12(s1,s2);
+// npoints = fitter->GetNpoints();
+// chi2v12 = TMath::Sqrt(fitter->GetChisquare()/npoints);
- //
- fitter = GetFitter9(s1,s2);
- npoints = fitter->GetNpoints();
- chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
- //
- fitter = GetFitter6(s1,s2);
- npoints = fitter->GetNpoints();
- chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
- fitter->GetParameters(param6Diff);
- //
- GetTransformation6(s1,s2,m6);
- GetTransformation9(s1,s2,m9);
- GetTransformation12(s1,s2,m12);
- //
- fitter = GetFitter6(s1,s2);
- //fitter->FixParameter(3,0);
- //fitter->Eval();
- GetTransformation6(s1,s2,m6FX);
+// //
+// fitter = GetFitter9(s1,s2);
+// npoints = fitter->GetNpoints();
+// chi2v9 = TMath::Sqrt(fitter->GetChisquare()/npoints);
+// //
+// fitter = GetFitter6(s1,s2);
+// npoints = fitter->GetNpoints();
+// chi2v6 = TMath::Sqrt(fitter->GetChisquare()/npoints);
+// fitter->GetParameters(param6Diff);
+// //
+// GetTransformation6(s1,s2,m6);
+// GetTransformation9(s1,s2,m9);
+// GetTransformation12(s1,s2,m12);
+// //
+// fitter = GetFitter6(s1,s2);
+// //fitter->FixParameter(3,0);
+// //fitter->Eval();
+// GetTransformation6(s1,s2,m6FX);
//
TH1 * his=0;
his = GetHisto(kY,s1,s2);
//
}
+ AliMagF* magF= (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+ if (!magF) AliError("Magneticd field - not initialized");
+ Double_t bz = magF->SolenoidField()/10.; //field in T
- // x2 = a00*x1 + a01*y1 + a02*z1 + a03
- // y2 = a10*x1 + a11*y1 + a12*z1 + a13
- // z2 = a20*x1 + a21*y1 + a22*z1 + a23
- // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
- // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/(a00 + a01*dydx1 + a02*dzdx1)
- //
- // a00 a01 a02 a03 p[0] p[1] p[2] p[9]
- // a10 a11 a12 a13 ==> p[3] p[4] p[5] p[10]
- // a20 a21 a22 a23 p[6] p[7] p[8] p[11]
-
- //
- //
- // dy:-(fXIO*m6.fElements[4]+m6.fElements[7])
- //
- // dphi:-(m6.fElements[4])
- //
- // dz:fXIO*m6.fElements[8]+m6.fElements[11]
- //
- // dtheta:m6.fElements[8]
- //
cstream<<"Align"<<
+ "run="<<fRun<< // run
+ "bz="<<bz<<
"s1="<<s1<< // reference sector
"s2="<<s2<< // sector to align
"m6FX.="<<&m6FX<< // tranformation matrix
//_____________________________________________________________________
-Long64_t AliTPCcalibAlign::Merge(TCollection* list) {
+Long64_t AliTPCcalibAlign::Merge(TCollection* const list) {
//
// merge function
//
UpdateKalman(*fSectorParamA,*fSectorCovarA,*align->fSectorParamA,*align->fSectorCovarA);
UpdateKalman(*fSectorParamC,*fSectorCovarC,*align->fSectorParamC,*align->fSectorCovarC);
}
+ if (!fClusterDelta[0]) MakeResidualHistos();
+
+ for (Int_t i=0; i<2; i++){
+ if (align->fClusterDelta[i]){
+ fClusterDelta[i]->Add(align->fClusterDelta[i]);
+ // align->fClusterDelta[i]->GetAxis(0)->SetRangeUser(-0.87,0.87);
+// align->fClusterDelta[i]->GetAxis(3)->SetRangeUser(-0.87,0.87);
+// fClusterDelta[i]->GetAxis(0)->SetRangeUser(-0.87,0.87);
+// fClusterDelta[i]->GetAxis(3)->SetRangeUser(-0.87,0.87);
+// Int_t idim[4]={0,1,2,3};
+// THnSparse *htemp=align->fClusterDelta[i]->Projection(4,idim);
+// THnSparse *htemp1=fClusterDelta[i]->Projection(4,idim);
+// htemp1->Add(htemp);
+// delete fClusterDelta[i];
+// fClusterDelta[i]=htemp1;
+// delete htemp;
+ }
+ }
+
+ for (Int_t i=0; i<4; i++){
+ if (!fTrackletDelta[i] && align->fTrackletDelta[i]) {
+ fTrackletDelta[i]= (THnSparse*)(align->fTrackletDelta[i]->Clone());
+ continue;
+ }
+ if (align->fTrackletDelta[i]) {
+ fTrackletDelta[i]->Add(align->fTrackletDelta[i]);
+ //
+ // align->fTrackletDelta[i]->GetAxis(3)->SetRangeUser(-0.36,0.36);
+// align->fTrackletDelta[i]->GetAxis(4)->SetRangeUser(-0.87,0.87);
+// fTrackletDelta[i]->GetAxis(3)->SetRangeUser(-0.36,0.36);
+// fTrackletDelta[i]->GetAxis(4)->SetRangeUser(-0.87,0.87);
+// //
+// Int_t idim[9]={0,1,2,3,4,5,6,7,8};
+// THnSparse *htemp=align->fTrackletDelta[i]->Projection(9,idim);
+// THnSparse *htemp1=fTrackletDelta[i]->Projection(9,idim);
+// htemp1->Add(htemp);
+// delete fTrackletDelta[i];
+// fTrackletDelta[i]=htemp1;
+// delete htemp;
+ }
+ }
+
}
Double_t AliTPCcalibAlign::Correct(Int_t type, Int_t value, Int_t s1, Int_t s2, Double_t x1, Double_t y1, Double_t z1, Double_t dydx1,Double_t dzdx1){
//
// full track fit parameters
//
- TLinearFitter fyf(2,"pol1");
- TLinearFitter fzf(2,"pol1");
+ static TLinearFitter fyf(2,"pol1"); // change to static - suggestion of calgrind - 30 % of time
+ static TLinearFitter fzf(2,"pol1"); // relative to time of given class
TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
TMatrixD covY(4,4),covZ(4,4);
Double_t chi2FacY =1;
return nf;
}
+void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){
+ //
+ // Update the cluster residula histograms for setup with field
+ // Kalman track fitting is used
+ // Only high momenta primary tracks used
+ //
+ // 1. Apply selection
+ // 2. Refit the track - in-out
+ // 3. Refit the track - out-in
+ // 4. Combine In and Out track - - fil cluster residuals
+ //
+ const Double_t kPtCut=1.0; // pt
+ const Double_t kSnpCut=0.2; // snp cut
+ const Double_t kNclCut=120; //
+ const Double_t kVertexCut=1;
+ const Double_t kMaxDist=0.5; // max distance between tracks and cluster
+ const Double_t kEdgeCut = 2.5;
+ const Double_t kDelta2=0.2*0.2; // initial increase in covar matrix
+ const Double_t kSigma=0.3; // error increase towards edges of TPC
+ const Double_t kSkipBoundary=7.5; // skip track updates in the boundary IFC,OFC, IO
+ //
+ if (!fCurrentTrack) return;
+ if (!fCurrentFriendTrack) return;
+ Float_t vertexXY=0,vertexZ=0;
+ fCurrentTrack->GetImpactParameters(vertexXY,vertexZ);
+ if (TMath::Abs(vertexXY)>kVertexCut) return;
+ if (TMath::Abs(vertexZ)>kVertexCut) return;
+ if (TMath::Abs(seed->Pt())<kPtCut) return;
+ if (seed->GetNumberOfClusters()<kNclCut) return;
+ if (TMath::Abs(seed->GetSnp())>kSnpCut) return;
+ if (!fClusterDelta[0]) MakeResidualHistos();
+ //
+ AliExternalTrackParam fitIn[160];
+ AliExternalTrackParam fitOut[160];
+ AliTPCROC * roc = AliTPCROC::Instance();
+ Double_t xmiddle = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+ Double_t xDiff = ( -roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+ Double_t xIFC = ( roc->GetPadRowRadii(0,0));
+ Double_t xOFC = ( roc->GetPadRowRadii(36,roc->GetNRows(36)-1));
+ //
+ Int_t detector=-1;
+ //
+ //
+ AliExternalTrackParam trackIn = *(fCurrentTrack->GetInnerParam());
+ AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut());
+ trackIn.ResetCovariance(10);
+ trackOut.ResetCovariance(10);
+ Double_t *covarIn = (Double_t*)trackIn.GetCovariance();
+ Double_t *covarOut = (Double_t*)trackOut.GetCovariance();
+ covarIn[0]+=kDelta2; covarIn[2]+=kDelta2;
+ covarIn[5]+=kDelta2/(100.*100.); covarIn[9]=kDelta2/(100.*100.);
+ covarIn[14]+=kDelta2/(5.*5.);
+ covarOut[0]+=kDelta2; covarOut[2]+=kDelta2;
+ covarOut[5]+=kDelta2/(100.*100.); covarOut[9]=kDelta2/(100.*100.);
+ covarOut[14]+=kDelta2/(5.*5.);
+ //
+ static Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+ //
+ Int_t ncl=0;
+ for (Int_t irow=0; irow<160; irow++){
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ if (detector!=cl->GetDetector()%36) return; // cluster from different sectors
+ // skip such tracks
+ ncl++;
+ }
+ if (ncl<kNclCut) return;
+ Int_t nclIn=0,nclOut=0;
+ Double_t xyz[3];
+ //
+ // Refit out - store residual maps
+ //
+ for (Int_t irow=0; irow<160; irow++){
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ Int_t sector = cl->GetDetector();
+ Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
+ if (cl->GetDetector()%36!=detector) continue;
+ if (TMath::Abs(dalpha)>0.01){
+ if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
+ }
+ Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
+ Double_t cov[3]={0.1,0.,0.1};
+ Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackOut.GetY());
+ Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
+ dmiddle*=dmiddle;
+ //
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma/dedge; // bigger error close to the boundary
+ cov[2]+=kSigma/dedge; // bigger error close to the boundary
+ cov[0]*=cov[0];
+ cov[2]*=cov[2];
+ if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
+ if (TMath::Abs(dedge)<kEdgeCut) continue;
+ //
+ Bool_t doUpdate=kTRUE;
+ if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
+ //
+ if (TMath::Abs(cl->GetY()-trackOut.GetY())<kMaxDist){
+ nclOut++;
+ if (doUpdate) trackOut.Update(&r[1],cov);
+ }
+ fitOut[irow]=trackOut;
+ }
+
+ //
+ // Refit In - store residual maps
+ //
+ for (Int_t irow=159; irow>=0; irow--){
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ Int_t sector = cl->GetDetector();
+ Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
+ if (cl->GetDetector()%36!=detector) continue;
+ if (TMath::Abs(dalpha)>0.01){
+ if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
+ }
+ Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
+ Double_t cov[3]={0.1,0.,0.1};
+ Double_t dedge = cl->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(trackIn.GetY());
+ Double_t dmiddle = TMath::Abs(cl->GetX()-xmiddle)/xDiff;
+ dmiddle*=dmiddle;
+ //
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[2]+=kSigma*dmiddle; // bigger error at boundary
+ cov[0]+=kSigma/dedge; // bigger error close to the boundary
+ cov[2]+=kSigma/dedge; // bigger error close to the boundary
+ cov[0]*=cov[0];
+ cov[2]*=cov[2];
+ if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
+ if (TMath::Abs(dedge)<kEdgeCut) continue;
+ Bool_t doUpdate=kTRUE;
+ if (TMath::Abs(cl->GetX()-xIFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-xOFC)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetX()-fXIO)<kSkipBoundary) doUpdate=kFALSE;
+ if (TMath::Abs(cl->GetY()-trackIn.GetY())<kMaxDist){
+ nclIn++;
+ if (doUpdate) trackIn.Update(&r[1],cov);
+ }
+ fitIn[irow]=trackIn;
+ }
+ //
+ //
+ for (Int_t irow=159; irow>=0; irow--){
+ //
+ // Update kalman - +- direction
+ // Store cluster residuals
+ AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+ if (!cl) continue;
+ if (cl->GetX()<80) continue;
+ if (detector<0) detector=cl->GetDetector()%36;
+ if (cl->GetDetector()%36!=detector) continue;
+ AliExternalTrackParam trackSmooth = fitIn[irow];
+ AliTrackerBase::UpdateTrack(trackSmooth, fitOut[irow]);
+ //
+ Double_t resVector[5];
+ trackSmooth.GetXYZ(xyz);
+ resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
+ if (resVector[1]<0) resVector[1]+=18;
+ resVector[2]= TMath::Sqrt(cl->GetX()*cl->GetX()+cl->GetY()*cl->GetY());
+ resVector[3]= cl->GetZ()/resVector[2];
+ //
+ resVector[0]= cl->GetY()-trackSmooth.GetY();
+ fClusterDelta[0]->Fill(resVector);
+ resVector[0]= cl->GetZ()-trackSmooth.GetZ();
+ fClusterDelta[1]->Fill(resVector);
+ }
+
+}
+
+
void AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
//
- // Update Kalman filter of Alignment
+ // Update Kalman filter of Alignment - only setup without filed
// IROC - OROC quadrants
//
-
- const Int_t kMinClusterF=40;
+ if (TMath::Abs(AliTracker::GetBz())>0.5) return;
+ if (!fClusterDelta[0]) MakeResidualHistos();
+ // const Int_t kMinClusterF=40;
+ const Int_t kMinClusterFit=10;
const Int_t kMinClusterQ=10;
//
- const Int_t kdrow1 =8; // rows to skip at the end
- const Int_t kdrow0 =2; // rows to skip at beginning
+ const Int_t kdrow1Fit =5; // rows to skip from fit at the end
+ const Int_t kdrow0Fit =10; // rows to skip from fit at beginning
const Float_t kedgey=3.0;
- const Float_t kMaxDist=0.5;
+ const Float_t kMaxDist=1;
const Float_t kMaxCorrY=0.05;
const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF
isec = isec%36; // use the hardware numbering
//
// full track fit parameters
//
- TLinearFitter fyf(2,"pol1");
- TLinearFitter fzf(2,"pol1");
+ static TLinearFitter fyf(2,"pol1"); // make it static - too much time for comiling of formula
+ static TLinearFitter fzf(2,"pol1"); // calgrind recomendation
TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
+ TVectorD pyfc(2),pzfc(2);
Int_t nf=0;
//
// make full fit as reference
//
for (Int_t iter=0; iter<2; iter++){
fyf.ClearPoints();
- for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ fzf.ClearPoints();
+ for (Int_t irow=kdrow0Fit;irow<159-kdrow1Fit;irow++) {
AliTPCclusterMI *c=track->GetClusterPointer(irow);
if (!c) continue;
if ((c->GetDetector()%36)!=isec) continue;
- if (c->GetRow()<kdrow0) continue;
+ if (c->GetRow()<kdrow0Fit) continue;
Double_t dx = c->GetX()-fXmiddle;
Double_t x[2]={dx, dx*dx};
if (iter==0 &&c->GetType()<0) continue;
if (iter==1){
Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t zfit = pzf[0]+pzf[1]*dx;
Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
if (dedge<kedgey) continue;
Double_t corrtrY =
corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
}
- fyf.AddPoint(x,c->GetY(),0.1);
- fzf.AddPoint(x,c->GetZ(),0.1);
+ if (TMath::Abs(x[0])<10){
+ fyf.AddPoint(x,c->GetY(),0.1); //use only middle rows+-10cm
+ fzf.AddPoint(x,c->GetZ(),0.1);
+ }
}
nf = fyf.GetNpoints();
- if (nf<kMinClusterF) return; // not enough points - skip
+ if (fyf.GetNpoints()<kMinClusterFit) return; // not enough points - skip
+ if (fzf.GetNpoints()<kMinClusterFit) return; // not enough points - skip
fyf.Eval();
fyf.GetParameters(pyf);
fyf.GetErrors(peyf);
fzf.GetErrors(pezf);
}
//
+ //
+ //
+ TVectorD vecX(160); // x vector
+ TVectorD vecY(160); // residuals vector
+ TVectorD vecZ(160); // residuals vector
+ TVectorD vPosG(3); //vertex position
+ TVectorD vPosL(3); // vertex position in the TPC local system
+ TVectorF vImpact(2); //track impact parameter
+ // Double_t tofSignal= fCurrentTrack->GetTOFsignal(); // tof signal
+ TVectorF tpcPosG(3); // global position of track at the middle of fXmiddle
+ Double_t lphi = TMath::ATan2(pyf[0],fXmiddle); // expected local phi angle - if vertex at 0
+ Double_t gphi = 2.*TMath::Pi()*(isec%18+0.5)/18.+lphi; // expected global phi if vertex at 0
+ Double_t ky = pyf[0]/fXmiddle;
+ Double_t kyE =0, kzE=0; // ky and kz expected
+ Double_t alpha =2.*TMath::Pi()*(isec%18+0.5)/18.;
+ Double_t scos=TMath::Cos(alpha);
+ Double_t ssin=TMath::Sin(alpha);
+ const AliESDVertex* vertex = fCurrentEvent->GetPrimaryVertexTracks();
+ vertex->GetXYZ(vPosG.GetMatrixArray());
+ fCurrentTrack->GetImpactParameters(vImpact[0],vImpact[1]); // track impact parameters
+ //
+ tpcPosG[0]= scos*fXmiddle-ssin*pyf[0];
+ tpcPosG[1]=+ssin*fXmiddle+scos*pyf[0];
+ tpcPosG[2]=pzf[0];
+ vPosL[0]= scos*vPosG[0]+ssin*vPosG[1];
+ vPosL[1]=-ssin*vPosG[0]+scos*vPosG[1];
+ vPosL[2]= vPosG[2];
+ kyE = (pyf[0]-vPosL[1])/(fXmiddle-vPosL[0]);
+ kzE = (pzf[0]-vPosL[2])/(fXmiddle-vPosL[0]);
+ //
+ // get constrained parameters
+ //
+ Double_t xvertex=vPosL[0]-fXmiddle;
+ fyf.AddPoint(&xvertex,vPosL[1], 0.00001);
+ fzf.AddPoint(&xvertex,vPosL[2], 2.);
+ fyf.Eval();
+ fyf.GetParameters(pyfc);
+ fzf.Eval();
+ fzf.GetParameters(pzfc);
+ //
+ //
// Make Fitters and params for 5 fitters
// 1-4 OROC quadrants
// 0 IROC
//
- TLinearFitter *fittersY[5];
- TLinearFitter *fittersZ[5];
+ static TLinearFitter *fittersY[5]={0,0,0,0,0}; // calgrind recomendation - fater to clear points
+ static TLinearFitter *fittersZ[5]={0,0,0,0,0}; // than create the fitter
+ if (fittersY[0]==0){
+ for (Int_t i=0;i<5;i++) {
+ fittersY[i] = new TLinearFitter(2,"pol1");
+ fittersZ[i] = new TLinearFitter(2,"pol1");
+ }
+ }
+ //
Int_t npoints[5];
TVectorD paramsY[5];
TVectorD errorsY[5];
Double_t chi2CZ[5];
for (Int_t i=0;i<5;i++) {
npoints[i]=0;
- fittersY[i] = new TLinearFitter(2,"pol1");
paramsY[i].ResizeTo(2);
errorsY[i].ResizeTo(2);
covY[i].ResizeTo(2,2);
- fittersZ[i] = new TLinearFitter(2,"pol1");
paramsZ[i].ResizeTo(2);
errorsZ[i].ResizeTo(2);
covZ[i].ResizeTo(2,2);
+ fittersY[i]->ClearPoints();
+ fittersZ[i]->ClearPoints();
}
//
// Update fitters
//
- for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ Int_t countRes=0;
+ for (Int_t irow=0;irow<159;irow++) {
AliTPCclusterMI *c=track->GetClusterPointer(irow);
if (!c) continue;
if ((c->GetDetector()%36)!=isec) continue;
- if (c->GetRow()<kdrow0) continue;
Double_t dx = c->GetX()-fXmiddle;
Double_t x[2]={dx, dx*dx};
Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t zfit = pzf[0]+pzf[1]*dx;
+ Double_t yfitC = pyfc[0]+pyfc[1]*dx;
+ Double_t zfitC = pzfc[0]+pzfc[1]*dx;
Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (TMath::Abs(c->GetZ()-zfit)>kMaxDist) continue;
if (dedge<kedgey) continue;
Double_t corrtrY =
corr->RPhiCOGCorrection(c->GetDetector(),c->GetRow(), c->GetPad(),
c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
-
+ //
+ vecX[countRes]=c->GetX();
+ vecY[countRes]=c->GetY()-yfit;
+ vecZ[countRes]=c->GetZ()-zfit;
+ countRes++;
+ //
+ // Fill THnSparse cluster residuals
+ // use only primary candidates with ITS signal
+ if (fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){
+ Double_t resVector[5];
+ resVector[1]= 9.*gphi/TMath::Pi();
+ resVector[2]= TMath::Sqrt(c->GetX()*c->GetX()+c->GetY()*c->GetY());
+ resVector[3]= c->GetZ()/resVector[2];
+ //
+ //
+ resVector[0]= c->GetY()-yfitC;
+ fClusterDelta[0]->Fill(resVector);
+ resVector[0]= c->GetZ()-zfitC;
+ fClusterDelta[1]->Fill(resVector);
+ }
+ if (c->GetRow()<kdrow0Fit) continue;
+ if (c->GetRow()>159-kdrow1Fit) continue;
+ //
+
if (c->GetDetector()>35){
if (c->GetX()<fXquadrant){
if (yfit<-kPRFWidth) fittersY[1]->AddPoint(x,c->GetY(),0.1);
covZ[i](1,1)*=chi2FacZ*chi2FacZ;
}
}
- for (Int_t i=0;i<5;i++){
- delete fittersY[i];
- delete fittersZ[i];
- }
-
//
// void UpdateSectorKalman
//
//
// Dump debug information
//
- if (fStreamLevel>0){
- TTreeSRedirector *cstream = GetDebugStreamer();
+ if (fStreamLevel>0){
+ // get vertex position
+ //
+ TTreeSRedirector *cstream = GetDebugStreamer();
+
+
if (cstream){
for (Int_t i0=0;i0<5;i0++){
for (Int_t i1=i0;i1<5;i1++){
"triggerClass="<<&fTriggerClass<< // trigger
"mag="<<fMagF<< // magnetic field
"isec="<<isec<< // current sector
+ //
+ "xref="<<fXmiddle<< // reference X
+ "vPosG.="<<&vPosG<< // vertex position in global system
+ "vPosL.="<<&vPosL<< // vertex position in local system
+ "vImpact.="<<&vImpact<< // track impact parameter
+ //"tofSignal="<<tofSignal<< // tof signal
+ "tpcPosG.="<<&tpcPosG<< // global position of track at the middle of fXmiddle
+ "lphi="<<lphi<< // expected local phi angle - if vertex at 0
+ "gphi="<<gphi<< // expected global phi if vertex at 0
+ "ky="<<ky<<
+ "kyE="<<kyE<< // expect ky - assiming pirmary track
+ "kzE="<<kzE<< // expected kz - assuming primary tracks
+ "salpha="<<alpha<< // sector alpha
+ "scos="<<scos<< // tracking cosinus
+ "ssin="<<ssin<< // tracking sinus
+ //
// full fit
+ //
"nf="<<nf<< // total number of points
"pyf.="<<&pyf<< // full OROC fit y
"pzf.="<<&pzf<< // full OROC fit z
+ "vX.="<<&vecX<< // x cluster
+ "vY.="<<&vecY<< // y residual cluster
+ "vZ.="<<&vecZ<< // z residual cluster
// quadrant and IROC fit
"i0="<<i0<< // quadrant number
"i1="<<i1<<
}
}
-void AliTPCcalibAlign::UpdateSectorKalman(Int_t sector, Int_t quadrant0, Int_t quadrant1, TMatrixD *p0, TMatrixD *c0, TMatrixD *p1, TMatrixD *c1 ){
+void AliTPCcalibAlign::UpdateSectorKalman(Int_t sector, Int_t quadrant0, Int_t quadrant1, TMatrixD *const p0, TMatrixD *const c0, TMatrixD *const p1, TMatrixD *const c1 ){
//
//
// tracks are refitted at sector middle
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff