X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TPC%2FAliTPCcalibAlign.cxx;h=6e1cfa6aa211f3ccabfaf1257c920f75c32a9b06;hp=40d8eda8aefe8cf7411eb3bd35f90bc0f7932bdb;hb=7eaa723ef27b1c82a19ffce6ec66303df9438c31;hpb=e4042305eece1ea7ad68e8dda2b76a957e06df1e

diff --git a/TPC/AliTPCcalibAlign.cxx b/TPC/AliTPCcalibAlign.cxx
index 40d8eda8aef..6e1cfa6aa21 100644
--- a/TPC/AliTPCcalibAlign.cxx
+++ b/TPC/AliTPCcalibAlign.cxx
@@ -33,6 +33,9 @@
 #include "AliExternalTrackParam.h"
 #include "AliTPCTracklet.h"
 #include "TH1D.h"
+#include "TVectorD.h"
+#include "TTreeStream.h"
+#include "TFile.h"
 
 #include <iostream>
 #include <sstream>
@@ -57,6 +60,26 @@ AliTPCcalibAlign::AliTPCcalibAlign()
   }
 }
 
+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
@@ -67,16 +90,19 @@ void AliTPCcalibAlign::Process(AliTPCseed *seed) {
   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]);
@@ -91,60 +117,82 @@ void AliTPCcalibAlign::Process(AliTPCseed *seed) {
     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;
+  
 
 
 
@@ -153,7 +201,7 @@ void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
   //
   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();
@@ -161,23 +209,45 @@ void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
   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
@@ -188,14 +258,14 @@ void AliTPCcalibAlign::Process12(Double_t *t1,
   //       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;
@@ -346,10 +416,10 @@ void AliTPCcalibAlign::Process6(Double_t *t1,
   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;
@@ -407,16 +477,46 @@ void AliTPCcalibAlign::Process6(Double_t *t1,
   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();