e.g:
kalmanFit0->SetInstance(kalmanFit0); //
kalmanFit0->InitTransformation(); //
- TF2 fxRZdz("fxRZdz","sign(y)*1000*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y)-AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y-1))",85,245,-250,250);
+ //
+ TF2 fxRZdz("fxRZdz","sign(y)*1000*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,0,x,0,y)-AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,0,x,0,y-1))",85,245,-250,250);
fxRZdz->Draw("");
+ TF2 fxRZ("fxRZ","sign(y)*10*(AliTPCkalmanFit::SGetTPCDeltaXYZ(0,-1,x,0,y))",85,245,-250,250);
+ fxRZ->Draw("");
+
*/
fCalibCovar = new TMatrixD(ncalibs,ncalibs);
for (Int_t icalib=0;icalib<ncalibs; icalib++){
AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
- (*fCalibParam)(icalib,0) = transform->fParam;
+ (*fCalibParam)(icalib,0) = transform->GetParam();
for (Int_t jcalib=0;jcalib<ncalibs; jcalib++){
if (icalib!=jcalib) (*fCalibCovar)(icalib,jcalib)= 0;
- if (icalib==jcalib) (*fCalibCovar)(icalib,jcalib) = transform->fSigma*transform->fSigma;
+ if (icalib==jcalib) (*fCalibCovar)(icalib,jcalib) = transform->GetSigma()*transform->GetSigma();
}
}
//
-void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug){
+void AliTPCkalmanFit::FitTrackLinear(AliTrackPointArray &points, TTreeSRedirector *debug, Float_t scalingRMSY, Float_t scalingRMSZ){
//
//
Double_t errz2 = chi2Z;
//set error - it is hack
Float_t *cov = (Float_t*) point.GetCov();
- cov[1] = erry2+kOff0;
- cov[2] = errz2+kOff0;
+ cov[1] = (erry2+kOff0)*scalingRMSY;
+ cov[2] = (errz2+kOff0)*scalingRMSZ;
UpdateLinear(point,debug);
if (!points.GetPoint(point,npoints-1-ipoint)) continue;
//set error - it is hack
cov = (Float_t*) point.GetCov();
- cov[1] = erry2+kOff0;
- cov[2] = errz2+kOff0;
+ cov[1] = (erry2+kOff0)*scalingRMSY;
+ cov[2] = (errz2+kOff0)*scalingRMSZ;
UpdateLinear(point,debug);
}
//
// save current param and covariance
for (Int_t i=0; i<ncalibs;i++){
AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(i);
- transform->fParam= (*fLinearParam)(i,0);
+ transform->SetParam( (*fLinearParam)(i,0));
(*fCalibParam)(i,0) = (*fLinearParam)(i,0);
for (Int_t j=0; j<ncalibs;j++){
(*fCalibCovar)(i,j) = (*fLinearCovar)(i,j);
//
}
+void AliTPCkalmanFit::AddCovariance(const char * varName, Double_t sigma){
+ //
+ //
+ //
+ if (!fCalibCovar) return;
+ if (!fCalibration) return;
+ if (!fCalibration->FindObject(varName)) return;
+ Int_t ncalibs = fCalibration->GetEntries();
+ TString strVar(varName);
+ for (Int_t icalib=0;icalib<ncalibs; icalib++){
+ AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
+ if (strVar.CompareTo(transform->GetName())==0){
+ (*fCalibCovar)(icalib,icalib)+=sigma*sigma;
+ }
+ }
+}
+
void AliTPCkalmanFit::PropagateTime(Int_t time){
//
// Propagate the calibration in time
// - Increase covariance matrix
//
- if (fCalibCovar) return;
+ if (!fCalibCovar) return;
Int_t ncalibs = fCalibration->GetEntries();
Double_t deltaT = 0;
if (fLastTimeStamp>0) deltaT = (fLastTimeStamp-time)/(60*60.); // delta T in hours
fLastTimeStamp = time;
for (Int_t icalib=0;icalib<ncalibs; icalib++){
AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
- (*fCalibCovar)(icalib,icalib)+= transform->fSigma2Time*deltaT;
+ if ((*fCalibCovar)(icalib,icalib)<transform->GetSigmaMax()*transform->GetSigmaMax())
+ (*fCalibCovar)(icalib,icalib)+= transform->GetSigma2Time()*TMath::Abs(deltaT);
}
}
}
+AliTrackPointArray * AliTPCkalmanFit::SortPoints(AliTrackPointArray &points){
+ //
+ //Creates the array - points sorted according radius - neccessay for kalman fit
+ //
+ //
+ // 0. choose the frame - rotation angle
+ //
+ Int_t npoints = points.GetNPoints();
+ if (npoints<1) return 0;
+ Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
+ Double_t ca = TMath::Cos(currentAlpha);
+ Double_t sa = TMath::Sin(currentAlpha);
+ //
+ // 1. sort the points
+ //
+ Double_t *rxvector = new Double_t[npoints];
+ Int_t *indexes = new Int_t[npoints];
+ for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
+ rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
+ }
+ TMath::Sort(npoints, rxvector,indexes,kFALSE);
+ AliTrackPoint point;
+ AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
+ for (Int_t ipoint=0; ipoint<npoints; ipoint++){
+ if (!points.GetPoint(point,indexes[ipoint])) continue;
+ pointsSorted->AddPoint(ipoint,&point);
+ }
+ delete [] rxvector;
+ delete [] indexes;
+ return pointsSorted;
+}
+
+
+
AliTrackPointArray * AliTPCkalmanFit::MakePointArrayLinear(Double_t alpha, Double_t y0, Double_t z0, Double_t ky, Double_t kz, Double_t err){
//
//
Int_t npoints=0;
for (Int_t i=0;i<6;i++) cov[i]=0.001;
for (Int_t i=0;i<500;i++){
- AliTrackPoint point(0, 0, 0, cov, -1,0,0);
+ AliTrackPoint point(0, 0, 0, cov, 0,0,0);
array.AddPoint(npoints, &point);
npoints++;
}
fCalibration->AddLast(calib);
}
+Int_t AliTPCkalmanFit::GetTransformationIndex(const char * trName){
+ //
+ //
+ //
+ if (!fCalibration) return -1;
+ if (!fCalibration->FindObject(trName)) return -1;
+ //
+ Int_t ncalibs = fCalibration->GetEntries();
+ TString strVar(trName);
+ for (Int_t icalib=0;icalib<ncalibs; icalib++){
+ AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
+ if (strVar.CompareTo(transform->GetName())==0){
+ return icalib;
+ }
+ }
+ return -1;
+}
+
+
+
void AliTPCkalmanFit::ApplyCalibration(AliTrackPointArray *array, Double_t csign){
//
//
Double_t dxdydz[3]={0,0,0};
for (Int_t icalib=0; icalib<ncalibs; icalib++){
AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
- dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
- dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
- dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->fParam, xyz[0], xyz[1], xyz[2]);
+ dxdydz[0] += transform->GetDeltaXYZ(0,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
+ dxdydz[1] += transform->GetDeltaXYZ(1,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
+ dxdydz[2] += transform->GetDeltaXYZ(2,volId, transform->GetParam(), xyz[0], xyz[1], xyz[2]);
}
((Float_t*)array->GetX())[ipoint]+=csign*dxdydz[0];
((Float_t*)array->GetY())[ipoint]+=csign*dxdydz[1];
}
}
}
- return (nrow>0);
+ return (nrow>0);
}
-Bool_t AliTPCkalmanFit::DumpCalib(const char *mask){
+Bool_t AliTPCkalmanFit::DumpCalib(const char *mask, Float_t correlationCut){
//
// Print calibration entries - name, value, error
//
TMatrixD &mat = *fCalibCovar;
Int_t nrow= mat.GetNrows();
TString strMask(mask);
+ TVectorD vecCorrSum(nrow);
+ for (Int_t irow=0; irow<nrow; irow++){
+ vecCorrSum[irow]=0;
+ for (Int_t icol=0; icol<nrow; icol++){
+ if (icol==irow) continue;
+ Double_t diag = TMath::Sqrt(TMath::Abs(mat(irow,irow)*mat(icol,icol)));
+ Double_t corr0 = mat(irow,icol)/diag;
+ vecCorrSum[irow]+=TMath::Abs(corr0);
+ }
+ vecCorrSum[irow]*=0.5;
+ }
+
for (Int_t irow=0; irow<nrow; irow++){
AliTPCTransformation * trans0 = GetTransformation(irow);
TString strName(trans0->GetName());
if (mask){
if (!strName.Contains(mask)) continue;
}
- printf("%d\t%s\t%f\t%f\n",
+ if (vecCorrSum[irow]<correlationCut) continue;
+ printf("%d\t%s\t%f\t%f\t%f\n",
irow,
trans0->GetName(),
(*fCalibParam)(irow,0),
- TMath::Sqrt(mat(irow,irow)));
+ TMath::Sqrt(mat(irow,irow)), vecCorrSum[irow]);
}
return (nrow>0);
}
if (ipar0+ipar1==0) continue;
Double_t param = (gRandom->Rndm()-0.5)*0.5; // generate random parameters
char tname[100];
- sprintf(tname,"tscalingR%d%dSide%d",ipar0,ipar1,iside);
+ snprintf(tname,100,"tscalingR%d%dSide%d",ipar0,ipar1,iside);
transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,1);
transformation->SetParams(0,5*0.25,0,&fpar);
kalmanFit0->AddCalibration(transformation);
}
-Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
+// Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
+// //
+// // function for visualization purposes
+// //
+// if (!fCalibration) return 0;
+// Int_t ncalibs = fCalibration->GetEntries();
+// if (ncalibs==0) return 0;
+// Double_t dxdydz[3]={0,0,0};
+// //
+// if (volID<0){
+// Double_t alpha = TMath::ATan2(y,x);
+// Double_t r = TMath::Sqrt(y*y+x*x);
+// volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
+// if (volID<0) volID+=18;
+// if (z<0) volID+=18;
+// if (r>120) volID+=36;
+// }
+// for (Int_t icalib=0; icalib<ncalibs; icalib++){
+// AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
+// Double_t param = (*fCalibParam)(icalib,0);
+// dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z);
+// }
+// return dxdydz[coord];
+// }
+
+// Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
+// //
+// //
+// //
+// if (AliTPCkalmanFit::fgInstance==0) return 0;
+// return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
+// }
+
+
+
+
+
+Double_t AliTPCkalmanFit::GetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
//
// function for visualization purposes
//
+ // coord - coordinate for output
+ // - 0 -X
+ // 1 -Y
+ // 2 -Z
+ // 3 -R
+ // 4 -RPhi
+ // 5 -Z
+ //
+ //icoordsys - type of coordinate system for input
+ // - 0 - x,y,z
+ // - 1 - r,phi,z
+ //
if (!fCalibration) return 0;
Int_t ncalibs = fCalibration->GetEntries();
if (ncalibs==0) return 0;
- Double_t dxdydz[3]={0,0,0};
+ Double_t xyz[3]={0,0,0};
+ Double_t dxdydz[6]={0,0,0,0,0,0};
+ Double_t alpha=0;
+ Double_t r=0;
+ if(icoordsys==0){alpha=TMath::ATan2(y,x); r =TMath::Sqrt(y*y+x*x);}
+ if(icoordsys==1){alpha=y; r=x;}
+ Double_t ca = TMath::Cos(alpha);
+ Double_t sa = TMath::Sin(alpha);
+ if(icoordsys==0){xyz[0]=x; xyz[1]=y; xyz[2]=z;}
+ if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;}
//
if (volID<0){
- Double_t alpha = TMath::ATan2(y,x);
- Double_t r = TMath::Sqrt(y*y+x*x);
+ // Double_t alpha = TMath::ATan2(y,x);
+ //Double_t r = TMath::Sqrt(y*y+x*x);
volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
if (volID<0) volID+=18;
if (z<0) volID+=18;
for (Int_t icalib=0; icalib<ncalibs; icalib++){
AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
Double_t param = (*fCalibParam)(icalib,0);
- dxdydz[coord] += transform->GetDeltaXYZ(coord,volID, param, x, y,z);
+ for (Int_t icoord=0;icoord<6;icoord++){
+ dxdydz[icoord] += transform->GetDeltaXYZ(icoord,volID, param, xyz[0],xyz[1],xyz[2]);
+ }
}
+
return dxdydz[coord];
}
-Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Double_t x, Double_t y, Double_t z){
+
+Double_t AliTPCkalmanFit::SGetTPCDeltaXYZ(Int_t coord, Int_t volID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
//
//
//
if (AliTPCkalmanFit::fgInstance==0) return 0;
- return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID,x,y,z);
+ return AliTPCkalmanFit::fgInstance->GetTPCDeltaXYZ(coord, volID, icoordsys,x,y,z);
+}
+
+
+
+
+
+Double_t AliTPCkalmanFit::GetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
+
+ Int_t ncalibs = fCalibration->GetEntries();
+ if (calibID>=ncalibs) return 0;
+ //Int_t volID=-1;
+ //Double_t xyz[3]={x,y,z};
+ Double_t r=0;
+ Double_t alpha=0;
+ if(icoordsys==0){r=TMath::Sqrt(x*x+y*y); alpha = TMath::ATan2(y,x);}
+ if(icoordsys==1){r=x; alpha = y;}
+ Double_t ca = TMath::Cos(alpha);
+ Double_t sa = TMath::Sin(alpha);
+ Double_t xyz[3]={0,0,0};
+ if(icoordsys==0){xyz[0]=x;xyz[1]=y;xyz[2]=z;}
+ if(icoordsys==1){xyz[0]=x*ca; xyz[1]=x*sa; xyz[2]=z;}
+ //xyz[3]=param; xyz[4]=volID;
+
+ if (volID<0){
+ //Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
+ //Double_t r = TMath::Sqrt(xyz[1]*xyz[1]+xyz[0]*xyz[0]);
+ volID = TMath::Nint(9*alpha/TMath::Pi()-0.5);
+ if (volID<0) volID+=18;
+ if (xyz[2]<0) volID+=18;
+ if (r>120) volID+=36;
+ }
+ AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(calibID);
+ //transform->SetInstance(transform);
+ Double_t param = (*fCalibParam)(calibID,0);
+ Double_t delta = (Double_t)transform->GetDeltaXYZ(coord,volID, param, xyz[0],xyz[1],xyz[2]);
+
+ return delta;
}
+Double_t AliTPCkalmanFit::SGetTPCtransXYZ(Int_t coord, Int_t volID, Int_t calibID, Int_t icoordsys, Double_t x, Double_t y, Double_t z){
+ //
+ //
+ //
+ if (AliTPCkalmanFit::fgInstance==0) return 0;
+ return AliTPCkalmanFit::fgInstance->GetTPCtransXYZ(coord, volID, calibID,icoordsys,x,y,z);
+}
+
+
+void AliTPCkalmanFit::MakeTreeTrans(TTreeSRedirector *debug, const char *treeName){
+ //
+ // Make the Tree before and after current calibration
+ //
+ if (!fCalibParam) {
+ AliError("Kalman Fit not initialized");
+ return;
+ }
+ //
+ //
+ //
+ const Int_t ncalibs = fCalibration->GetEntries();
+ TMatrixD dxdydz(ncalibs,5);
+ Double_t * adx = new Double_t[ncalibs];
+ Double_t * ady = new Double_t[ncalibs];
+ Double_t * adz = new Double_t[ncalibs];
+ Double_t * adr = new Double_t[ncalibs];
+ Double_t * adrphi = new Double_t[ncalibs];
+
+ Double_t x[3];
+ for (x[0]=-250.;x[0]<=250.;x[0]+=10.){
+ for (x[1]=-250.;x[1]<=250.;x[1]+=10.){
+ for (x[2]=-250.;x[2]<=250.;x[2]+=20.) {
+ Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
+ if (r<20) continue;
+ if (r>260) continue;
+ //Double_t z = x[2];
+ Double_t phi=TMath::ATan2(x[1],x[0]);
+ Double_t ca=TMath::Cos(phi);
+ Double_t sa=TMath::Sin(phi);
+ Double_t dx=0;
+ Double_t dy=0;
+ Double_t dz=0;
+ Double_t dr=0;
+ Double_t rdphi=0;
+
+ Int_t volID= TMath::Nint(9*phi/TMath::Pi()-0.5);
+ if (volID<0) volID+=18;
+ if (x[2]<0) volID+=18; //C-side
+ if (r>120) volID+=36; //outer
+ Double_t volalpha=(volID+0.5)*TMath::Pi()/9;
+ Double_t cva=TMath::Cos(volalpha);
+ Double_t sva=TMath::Sin(volalpha);
+
+ Double_t lx=x[0]*cva+x[1]*sva;
+ Double_t ly=-x[0]*sva+x[1]*cva;
+
+
+ for(Int_t icalib=0;icalib<ncalibs;icalib++){
+ for(Int_t icoord=0;icoord<5;icoord++){
+ dxdydz(icalib,icoord)= GetTPCtransXYZ(icoord, -1, icalib, 0, x[0], x[1], x[2]);
+ }
+ }
+ dx=GetTPCDeltaXYZ(0, -1, 0, x[0], x[1], x[2]);
+ dy=GetTPCDeltaXYZ(1, -1, 0, x[0], x[1], x[2]);
+ dz=GetTPCDeltaXYZ(2, -1, 0, x[0], x[1], x[2]);
+ dr=GetTPCDeltaXYZ(3, -1, 0, x[0], x[1], x[2]);
+ rdphi=GetTPCDeltaXYZ(4, -1, 0, x[0], x[1], x[2]);
+
+
+ if(debug){
+ TTreeStream &cstream=
+ (*debug)<<treeName<<
+ "x="<<x[0]<<
+ "y="<<x[1]<<
+ "z="<<x[2]<<
+ "r="<<r<<
+ "ca="<<ca<<
+ "sa="<<sa<<
+ "lx="<<lx<<
+ "ly="<<ly<<
+ "sector="<<volID<<
+ "phi="<<phi<<
+ "dx="<<dx<<
+ "dy="<<dy<<
+ "dz="<<dz<<
+ "dr="<<dr<<
+ "rdphi="<<rdphi<<
+ "dxdydz.="<<&dxdydz;
+ for (Int_t icalib=0;icalib<ncalibs;icalib++){
+ AliTPCTransformation * transform = (AliTPCTransformation *)fCalibration->At(icalib);
+ char tname[1000];
+ //
+ snprintf(tname,1000,"dx%s=",transform->GetName());
+ adx[icalib] =dxdydz(icalib,0);
+ cstream<<tname<<adx[icalib];
+ snprintf(tname,1000,"dy%s=",transform->GetName());
+ ady[icalib] =dxdydz(icalib,1);
+ cstream<<tname<<ady[icalib];
+ snprintf(tname,1000,"dz%s=",transform->GetName());
+ adz[icalib] =dxdydz(icalib,2);
+ cstream<<tname<<adz[icalib];
+ //
+ snprintf(tname,1000,"dr%s=",transform->GetName());
+ adr[icalib] =dxdydz(icalib,3);
+ cstream<<tname<<adr[icalib];
+ snprintf(tname,1000,"rdphi%s=",transform->GetName());
+ adrphi[icalib] =dxdydz(icalib,4);
+ cstream<<tname<<adrphi[icalib];
+ }
+ cstream<<"\n";
+ }
+ }
+ }
+ Printf("x0=%f finished",x[0]);
+ }
+ delete [] adx;// = new Double_t[ncalibs];
+ delete [] ady;// = new Double_t[ncalibs];
+ delete [] adz;// = new Double_t[ncalibs];
+ delete [] adr;// = new Double_t[ncalibs];
+ delete [] adrphi;// = new Double_t[ncalibs];
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff