// Class to make a internal alignemnt of TPC chambers //
//
// Different linear tranformation investigated
-// 12 parameters - arbitrary transformation
+// 12 parameters - arbitrary linear transformation
// 9 parameters - scaling fixed to 1
-// 6 parameters -
+// 6 parameters - x-y rotation x-z, y-z tiliting
////
//// Only the 12 parameter version works so far...
////
#include "AliExternalTrackParam.h"
#include "AliTPCTracklet.h"
#include "TH1D.h"
+#include "TVectorD.h"
+#include "TTreeStream.h"
+#include "TFile.h"
+#include "TF1.h"
#include <iostream>
#include <sstream>
}
}
+AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
+ :AliTPCcalibBase(),
+ fDphiHistArray(72*72),
+ fDthetaHistArray(72*72),
+ fDyHistArray(72*72),
+ fDzHistArray(72*72),
+ fFitterArray12(72*72),
+ fFitterArray9(72*72),
+ fFitterArray6(72*72)
+{
+ //
+ // Constructor
+ //
+ SetName(name);
+ SetTitle(title);
+ for (Int_t i=0;i<72*72;++i) {
+ fPoints[i]=0;
+ }
+}
+
AliTPCcalibAlign::~AliTPCcalibAlign() {
//
// destructor
}
-void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
- const AliExternalTrackParam &tp2,
- Int_t s1,Int_t s2) {
+void AliTPCcalibAlign::Analyze(){
+ //
+ // Analyze function
+ //
+ EvalFitters();
+}
- if (s2-s1==36) {//only inner-outer
- if (!fDphiHistArray[s1*72+s2]) {
- stringstream name;
- stringstream title;
- name<<"hist_phi_"<<s1<<"_"<<s2;
- title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
- fDphiHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.01,0.01); // +/- 10 mrad
- ((TH1D*)fDphiHistArray[s1*72+s2])->SetDirectory(0);
- }
- if (!fDthetaHistArray[s1*72+s2]) {
- stringstream name;
- stringstream title;
- name<<"hist_theta_"<<s1<<"_"<<s2;
- title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
- fDthetaHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.01,0.01); // +/- 10 mrad
- ((TH1D*)fDthetaHistArray[s1*72+s2])->SetDirectory(0);
- }
- if (!fDyHistArray[s1*72+s2]) {
- stringstream name;
- stringstream title;
- name<<"hist_y_"<<s1<<"_"<<s2;
- title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
- fDyHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.3,0.3); // +/- 3 mm
- ((TH1D*)fDyHistArray[s1*72+s2])->SetDirectory(0);
- }
- if (!fDzHistArray[s1*72+s2]) {
- stringstream name;
- stringstream title;
- name<<"hist_z_"<<s1<<"_"<<s2;
- title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
- fDzHistArray[s1*72+s2]=new TH1D(name.str().c_str(),title.str().c_str(),1024,-0.3,0.3); // +/- 3 mm
- ((TH1D*)fDzHistArray[s1*72+s2])->SetDirectory(0);
- }
- static_cast<TH1D*>(fDphiHistArray[s1*72+s2])
- ->Fill(TMath::ASin(tp1.GetSnp())
- -TMath::ASin(tp2.GetSnp()));
- static_cast<TH1D*>(fDthetaHistArray[s1*72+s2])
- ->Fill(TMath::ATan(tp1.GetTgl())
- -TMath::ATan(tp2.GetTgl()));
- static_cast<TH1D*>(fDyHistArray[s1*72+s2])
- ->Fill(tp1.GetY()
- -tp2.GetY());
- static_cast<TH1D*>(fDzHistArray[s1*72+s2])
- ->Fill(tp1.GetZ()
- -tp2.GetZ());
- }
- return;
+
+void AliTPCcalibAlign::Terminate(){
+ //
+ // Terminate function
+ // call base terminate + Eval of fitters
+ //
+ if (GetDebugLevel()>0) Info("AliTPCcalibAlign","Teminate");
+ EvalFitters();
+ AliTPCcalibBase::Terminate();
+}
+
+void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
+ const AliExternalTrackParam &tp2,
+ Int_t s1,Int_t s2) {
+
+ //
+ //
+ //
+ FillHisto(tp1,tp2,s1,s2);
//
// Process function to fill fitters
//
Double_t t1[5],t2[5];
Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
- Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3]/*, &dzdx2=t2[4]*/;
+ Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
x1 =tp1.GetX();
y1 =tp1.GetY();
z1 =tp1.GetZ();
dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
Double_t tgl1=tp1.GetTgl();
// dz/dx = 1/(cos(theta)*cos(phi))
- dzdx1=1./TMath::Sqrt((1.+tgl1*tgl1)*(1.-snp1*snp1));
- // x2 =tp1->GetX();
+ dzdx1=tgl1/TMath::Sqrt(1.-snp1*snp1);
+ x2 =tp2.GetX();
y2 =tp2.GetY();
z2 =tp2.GetZ();
Double_t snp2=tp2.GetSnp();
dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
Double_t tgl2=tp2.GetTgl();
- dzdx1=1./TMath::Sqrt((1.+tgl2*tgl2)*(1.-snp2*snp2));
-
+ dzdx2=tgl2/TMath::Sqrt(1.-snp2*snp2);
+ //
+ //
+ //
+ if (fStreamLevel>1){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ if (cstream){
+ static TVectorD vec1(5);
+ static TVectorD vec2(5);
+ vec1.SetElements(t1);
+ vec2.SetElements(t2);
+ AliExternalTrackParam *p1 = &((AliExternalTrackParam&)tp1);
+ AliExternalTrackParam *p2 = &((AliExternalTrackParam&)tp2);
+ (*cstream)<<"Tracklet"<<
+ "tp1.="<<p1<<
+ "tp2.="<<p2<<
+ "v1.="<<&vec1<<
+ "v2.="<<&vec2<<
+ "s1="<<s1<<
+ "s2="<<s2<<
+ "\n";
+ }
+ }
Process12(t1,t2,GetOrMakeFitter12(s1,s2));
- Process9(t1,t2,GetOrMakeFitter12(s1,s2));
- Process6(t1,t2,GetOrMakeFitter12(s1,s2));
- ++fPoints[72*s1+s2];
+ Process9(t1,t2,GetOrMakeFitter9(s1,s2));
+ Process6(t1,t2,GetOrMakeFitter6(s1,s2));
+ ++fPoints[GetIndex(s1,s2)];
}
-void AliTPCcalibAlign::Process12(Double_t *t1,
- Double_t *t2,
+void AliTPCcalibAlign::Process12(const Double_t *t1,
+ const Double_t *t2,
TLinearFitter *fitter) {
// x2 = a00*x1 + a01*y1 + a02*z1 + a03
// y2 = a10*x1 + a11*y1 + a12*z1 + a13
// 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]
- Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
- Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
+ const Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
+ const Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
// TODO:
- Double_t sy = 1.;
- Double_t sz = 1.;
- Double_t sdydx = 1.;
- Double_t sdzdx = 1.;
+ Double_t sy = 0.1;
+ Double_t sz = 0.1;
+ Double_t sdydx = 0.001;
+ Double_t sdzdx = 0.001;
Double_t p[12];
Double_t value;
Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
// TODO:
- Double_t sy = 1.;
- Double_t sz = 1.;
- Double_t sdydx = 1.;
- Double_t sdzdx = 1.;
+ Double_t sy = 0.1;
+ Double_t sz = 0.1;
+ Double_t sdydx = 0.001;
+ Double_t sdzdx = 0.001;
Double_t p[12];
Double_t value;
fitter->AddPoint(p,value,sdzdx);
}
-void AliTPCcalibAlign::Eval() {
+
+
+
+void AliTPCcalibAlign::EvalFitters() {
+ //
+ // Analyze function
+ //
+ // Perform the fitting using linear fitters
+ //
+ Int_t kMinPoints =50;
TLinearFitter *f;
+ TFile fff("alignDebug.root","recreate");
for (Int_t s1=0;s1<72;++s1)
- for (Int_t s2=0;s2<72;++s2)
- if ((f=GetFitter12(s1,s2))&&fPoints[72*s1+s2]>12) {
- // cerr<<s1<<","<<s2<<": "<<fPoints[72*s1+s2]<<endl;
- if (!f->Eval()) {
+ for (Int_t s2=0;s2<72;++s2){
+ if ((f=GetFitter12(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
+ // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
+ if (f->Eval()!=0) {
+ cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
+ f->Write(Form("f12_%d_%d",s1,s2));
+ }else{
+ f->Write(Form("f12_%d_%d",s1,s2));
+ }
+ }
+ if ((f=GetFitter9(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
+ // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
+ if (f->Eval()!=0) {
+ cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
+ }else{
+ f->Write(Form("f9_%d_%d",s1,s2));
+ }
+ }
+ if ((f=GetFitter6(s1,s2))&&fPoints[GetIndex(s1,s2)]>kMinPoints) {
+ // cerr<<s1<<","<<s2<<": "<<fPoints[GetIndex(s1,s2)]<<endl;
+ if (f->Eval()!=0) {
cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
+ }else{
+ f->Write(Form("f6_%d_%d",s1,s2));
}
}
+ }
+ this->Write("align");
/*
fitter->Eval();
*/
}
+TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter12(Int_t s1,Int_t s2) {
+ //
+ // get or make fitter - general linear transformation
+ //
+ static Int_t counter12=0;
+ static TF1 f12("f12","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]++x[9]++x[10]++x[11]");
+ TLinearFitter * fitter = GetFitter12(s1,s2);
+ if (fitter) return fitter;
+ // fitter =new TLinearFitter(12,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]++x[9]++x[10]++x[11]");
+ fitter =new TLinearFitter(&f12,"");
+ fitter->StoreData(kFALSE);
+ fFitterArray12.AddAt(fitter,GetIndex(s1,s2));
+ counter12++;
+ if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter12<<endl;
+ return fitter;
+}
+
+TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter9(Int_t s1,Int_t s2) {
+ //
+ //get or make fitter - general linear transformation - no scaling
+ //
+ static Int_t counter9=0;
+ static TF1 f9("f9","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
+ TLinearFitter * fitter = GetFitter9(s1,s2);
+ if (fitter) return fitter;
+ // fitter =new TLinearFitter(9,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]++x[6]++x[7]++x[8]");
+ fitter =new TLinearFitter(&f9,"");
+ fitter->StoreData(kFALSE);
+ fFitterArray9.AddAt(fitter,GetIndex(s1,s2));
+ counter9++;
+ if (GetDebugLevel()>0) cerr<<"Creating fitter12 "<<s1<<","<<s2<<" : "<<counter9<<endl;
+ return fitter;
+}
+
+TLinearFitter* AliTPCcalibAlign::GetOrMakeFitter6(Int_t s1,Int_t s2) {
+ //
+ // get or make fitter - 6 paramater linear tranformation
+ // - no scaling
+ // - rotation x-y
+ // - tilting x-z, y-z
+ static Int_t counter6=0;
+ static TF1 f6("f6","x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
+ TLinearFitter * fitter = GetFitter6(s1,s2);
+ if (fitter) return fitter;
+ // fitter=new TLinearFitter(6,"x[0]++x[1]++x[2]++x[3]++x[4]++x[5]");
+ fitter=new TLinearFitter(&f6,"");
+ fitter->StoreData(kFALSE);
+ fFitterArray6.AddAt(fitter,GetIndex(s1,s2));
+ counter6++;
+ if (GetDebugLevel()>0) cerr<<"Creating fitter6 "<<s1<<","<<s2<<" : "<<counter6<<endl;
+ return fitter;
+}
+
+
+
+
+
Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
+ //
+ // GetTransformation matrix - 12 paramaters - generael linear transformation
+ //
if (!GetFitter12(s1,s2))
return false;
else {
TVectorD p(12);
- cerr<<"foo"<<endl;
GetFitter12(s1,s2)->GetParameters(p);
- cerr<<"bar"<<endl;
a.ResizeTo(4,4);
- a[0][0]=p[0];
- a[0][1]=p[1];
- a[0][2]=p[2];
- a[0][3]=p[9];
- a[1][0]=p[3];
- a[1][1]=p[4];
- a[1][2]=p[5];
- a[1][3]=p[10];
- a[2][0]=p[6];
- a[2][1]=p[7];
- a[2][2]=p[8];
- a[2][3]=p[11];
- a[3][0]=0.;
- a[3][1]=0.;
- a[3][2]=0.;
- a[3][3]=1.;
+ a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
+ a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
+ a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
+ a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
return true;
}
}
Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
+ //
+ // GetTransformation matrix - 9 paramaters - general linear transformation
+ // No scaling
+ //
if (!GetFitter9(s1,s2))
return false;
else {
TVectorD p(9);
GetFitter9(s1,s2)->GetParameters(p);
a.ResizeTo(4,4);
- a[0][0]=p[0];
- a[0][1]=p[1];
- a[0][2]=p[2];
- a[0][3]=p[9];
- a[1][0]=p[3];
- a[1][1]=p[4];
- a[1][2]=p[5];
- a[1][3]=p[10];
- a[2][0]=p[6];
- a[2][1]=p[7];
- a[2][2]=p[8];
- a[2][3]=p[11];
- a[3][0]=0.;
- a[3][1]=0.;
- a[3][2]=0.;
- a[3][3]=1.;
+ a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
+ a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
+ a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
+ a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
return true;
}
}
Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
+ //
+ // GetTransformation matrix - 6 paramaters
+ // 3 translation
+ // 1 rotation -x-y
+ // 2 tilting x-z y-z
if (!GetFitter6(s1,s2))
return false;
else {
TVectorD p(6);
- cerr<<"foo"<<endl;
GetFitter6(s1,s2)->GetParameters(p);
- cerr<<"bar"<<endl;
a.ResizeTo(4,4);
- a[0][0]=p[0];
- a[0][1]=p[1];
- a[0][2]=p[2];
- a[0][3]=p[9];
- a[1][0]=p[3];
- a[1][1]=p[4];
- a[1][2]=p[5];
- a[1][3]=p[10];
- a[2][0]=p[6];
- a[2][1]=p[7];
- a[2][2]=p[8];
- a[2][3]=p[11];
- a[3][0]=0.;
- a[3][1]=0.;
- a[3][2]=0.;
- a[3][3]=1.;
+ a[0][0]=p[0]; a[0][1]=p[1]; a[0][2]=p[2]; a[0][3]=p[9];
+ a[1][0]=p[3]; a[1][1]=p[4]; a[1][2]=p[5]; a[1][3]=p[10];
+ a[2][0]=p[6]; a[2][1]=p[7]; a[2][2]=p[8]; a[2][3]=p[11];
+ a[3][0]=0.; a[3][1]=0.; a[3][2]=0.; a[3][3]=1.;
return true;
}
}
+
+void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam &tp1,
+ const AliExternalTrackParam &tp2,
+ Int_t s1,Int_t s2) {
+ //
+ // Fill residual histograms
+ //
+ if (s2-s1==36) {//only inner-outer
+ GetHisto(kPhi,s1,s2,kTRUE)->Fill(TMath::ASin(tp1.GetSnp())-TMath::ASin(tp2.GetSnp()));
+ GetHisto(kTheta,s1,s2,kTRUE)->Fill(TMath::ATan(tp1.GetTgl())-TMath::ATan(tp2.GetTgl()));
+ GetHisto(kY,s1,s2,kTRUE)->Fill(tp1.GetY()-tp2.GetY());
+ GetHisto(kZ,s1,s2,kTRUE)->Fill(tp1.GetZ()-tp2.GetZ());
+ }
+}
+
+
+
+TH1 * AliTPCcalibAlign::GetHisto(HistoType type, Int_t s1, Int_t s2, Bool_t force)
+{
+ //
+ // return specified residual histogram - it is only QA
+ // if force specified the histogram and given histogram is not existing
+ // new histogram is created
+ //
+ if (GetIndex(s1,s2)>=72*72) return 0;
+ TObjArray *histoArray=0;
+ switch (type) {
+ case kY:
+ histoArray = &fDyHistArray; break;
+ case kZ:
+ histoArray = &fDzHistArray; break;
+ case kPhi:
+ histoArray = &fDphiHistArray; break;
+ case kTheta:
+ histoArray = &fDthetaHistArray; break;
+ }
+ TH1 * histo= (TH1*)histoArray->At(GetIndex(s1,s2));
+ if (histo) return histo;
+ if (force==kFALSE) return 0;
+ //
+ stringstream name;
+ stringstream title;
+ switch (type) {
+ case kY:
+ name<<"hist_y_"<<s1<<"_"<<s2;
+ title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
+ histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
+ break;
+ case kZ:
+ name<<"hist_z_"<<s1<<"_"<<s2;
+ title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
+ histo = new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
+ break;
+ case kPhi:
+ name<<"hist_phi_"<<s1<<"_"<<s2;
+ title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
+ histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
+ break;
+ case kTheta:
+ name<<"hist_theta_"<<s1<<"_"<<s2;
+ title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
+ histo =new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
+ break;
+ }
+ histo->SetDirectory(0);
+ histoArray->AddAt(histo,GetIndex(s1,s2));
+ return histo;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff