#include "AliExternalTrackParam.h"
#include "AliTPCTracklet.h"
#include "TH1D.h"
+#include "TVectorD.h"
+#include "TTreeStream.h"
+#include "TFile.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
TObjArray tracklets=
AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
kFALSE,20,2);
- TObjArray trackletsL=
- AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
- kFALSE,20,2);
- TObjArray trackletsQ=
- AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
- kFALSE,20,2);
- TObjArray trackletsR=
- AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
- kFALSE,20,2);
+ // TObjArray trackletsL=
+// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kLinear,
+// kFALSE,20,2);
+// TObjArray trackletsQ=
+// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kQuadratic,
+// kFALSE,20,2);
+// TObjArray trackletsR=
+// AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kRiemann,
+// kFALSE,20,2);
tracklets.SetOwner();
+ // trackletsL.SetOwner();
+// trackletsQ.SetOwner();
+// trackletsR.SetOwner();
if (tracklets.GetEntries()==2) {
AliTPCTracklet *t1=static_cast<AliTPCTracklet*>(tracklets[0]);
AliTPCTracklet *t2=static_cast<AliTPCTracklet*>(tracklets[1]);
delete common1;
delete common2;
}
+
+}
+
+void AliTPCcalibAlign::Analyze(){
+ //
+ // Analyze function
+ //
+ EvalFitters();
+}
+
+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) {
+
if (s2-s1==36) {//only inner-outer
- if (!fDphiHistArray[s1*72+s2]) {
+ if (!fDphiHistArray[GetIndex(s1,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);
+ fDphiHistArray[GetIndex(s1,s2)]=new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
+ ((TH1D*)fDphiHistArray[GetIndex(s1,s2)])->SetDirectory(0);
}
- if (!fDthetaHistArray[s1*72+s2]) {
+ if (!fDthetaHistArray[GetIndex(s1,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);
+ fDthetaHistArray[GetIndex(s1,s2)]=new TH1D(name.str().c_str(),title.str().c_str(),512,-0.01,0.01); // +/- 10 mrad
+ ((TH1D*)fDthetaHistArray[GetIndex(s1,s2)])->SetDirectory(0);
}
- if (!fDyHistArray[s1*72+s2]) {
+ if (!fDyHistArray[GetIndex(s1,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);
+ fDyHistArray[GetIndex(s1,s2)]=new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
+ ((TH1D*)fDyHistArray[GetIndex(s1,s2)])->SetDirectory(0);
}
- if (!fDzHistArray[s1*72+s2]) {
+ if (!fDzHistArray[GetIndex(s1,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);
+ fDzHistArray[GetIndex(s1,s2)]=new TH1D(name.str().c_str(),title.str().c_str(),512,-0.3,0.3); // +/- 3 mm
+ ((TH1D*)fDzHistArray[GetIndex(s1,s2)])->SetDirectory(0);
}
- static_cast<TH1D*>(fDphiHistArray[s1*72+s2])
+ static_cast<TH1D*>(fDphiHistArray[GetIndex(s1,s2)])
->Fill(TMath::ASin(tp1.GetSnp())
-TMath::ASin(tp2.GetSnp()));
- static_cast<TH1D*>(fDthetaHistArray[s1*72+s2])
+ static_cast<TH1D*>(fDthetaHistArray[GetIndex(s1,s2)])
->Fill(TMath::ATan(tp1.GetTgl())
-TMath::ATan(tp2.GetTgl()));
- static_cast<TH1D*>(fDyHistArray[s1*72+s2])
+ static_cast<TH1D*>(fDyHistArray[GetIndex(s1,s2)])
->Fill(tp1.GetY()
-tp2.GetY());
- static_cast<TH1D*>(fDzHistArray[s1*72+s2])
+ static_cast<TH1D*>(fDzHistArray[GetIndex(s1,s2)])
->Fill(tp1.GetZ()
-tp2.GetZ());
}
- return;
+
//
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));
+ Process9(t1,t2,GetOrMakeFitter9(s1,s2));
+ Process6(t1,t2,GetOrMakeFitter6(s1,s2));
++fPoints[72*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) {
+ for (Int_t s2=0;s2<72;++s2){
+ if ((f=GetFitter12(s1,s2))&&fPoints[72*s1+s2]>kMinPoints) {
+ // cerr<<s1<<","<<s2<<": "<<fPoints[72*s1+s2]<<endl;
+ if (!f->Eval()) {
+ 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[72*s1+s2]>kMinPoints) {
// cerr<<s1<<","<<s2<<": "<<fPoints[72*s1+s2]<<endl;
if (!f->Eval()) {
cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
+ }else{
+ f->Write(Form("f9_%d_%d",s1,s2));
}
}
+ if ((f=GetFitter6(s1,s2))&&fPoints[72*s1+s2]>kMinPoints) {
+ // cerr<<s1<<","<<s2<<": "<<fPoints[72*s1+s2]<<endl;
+ if (!f->Eval()) {
+ cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
+ }else{
+ f->Write(Form("f6_%d_%d",s1,s2));
+ }
+ }
+ }
/*
fitter->Eval();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff