Bug fix - Teminate - as virtula function
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibAlign.cxx
index 8b1e8f29b83b6b4149cb6b037d7e4bcdac831a63..6e1cfa6aa211f3ccabfaf1257c920f75c32a9b06 100644 (file)
 #include "TLinearFitter.h"
 #include "AliTPCcalibAlign.h"
 #include "AliExternalTrackParam.h"
+#include "AliTPCTracklet.h"
+#include "TH1D.h"
+#include "TVectorD.h"
+#include "TTreeStream.h"
+#include "TFile.h"
 
 #include <iostream>
+#include <sstream>
 using namespace std;
 
 ClassImp(AliTPCcalibAlign)
 
 AliTPCcalibAlign::AliTPCcalibAlign()
-  :fFitterArray12(72*72),fFitterArray9(72*72),fFitterArray6(72*72)
+  :  fDphiHistArray(72*72),
+     fDthetaHistArray(72*72),
+     fDyHistArray(72*72),
+     fDzHistArray(72*72),
+     fFitterArray12(72*72),
+     fFitterArray9(72*72),
+     fFitterArray6(72*72)
 {
   //
   // Constructor
@@ -48,15 +60,142 @@ 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
   //
 }
 
-void AliTPCcalibAlign::Process(const AliExternalTrackParam &tp1,
-                              const AliExternalTrackParam &tp2,
-                              Int_t s1,Int_t s2) {
+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);
+  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]);
+    if (t1->GetSector()>t2->GetSector()) {
+      AliTPCTracklet* tmp=t1;
+      t1=t2;
+      t2=tmp;
+    }
+    AliExternalTrackParam *common1=0,*common2=0;
+    if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
+      ProcessTracklets(*common1,*common2,t1->GetSector(),t2->GetSector());
+    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[GetIndex(s1,s2)]) {
+      stringstream name;
+      stringstream title;
+      name<<"hist_phi_"<<s1<<"_"<<s2;
+      title<<"Phi Missalignment for sectors "<<s1<<" and "<<s2;
+      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[GetIndex(s1,s2)]) {
+      stringstream name;
+      stringstream title;
+      name<<"hist_theta_"<<s1<<"_"<<s2;
+      title<<"Theta Missalignment for sectors "<<s1<<" and "<<s2;
+      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[GetIndex(s1,s2)]) {
+      stringstream name;
+      stringstream title;
+      name<<"hist_y_"<<s1<<"_"<<s2;
+      title<<"Y Missalignment for sectors "<<s1<<" and "<<s2;
+      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[GetIndex(s1,s2)]) {
+      stringstream name;
+      stringstream title;
+      name<<"hist_z_"<<s1<<"_"<<s2;
+      title<<"Z Missalignment for sectors "<<s1<<" and "<<s2;
+      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[GetIndex(s1,s2)])
+      ->Fill(TMath::ASin(tp1.GetSnp())
+            -TMath::ASin(tp2.GetSnp()));
+    static_cast<TH1D*>(fDthetaHistArray[GetIndex(s1,s2)])
+      ->Fill(TMath::ATan(tp1.GetTgl())
+            -TMath::ATan(tp2.GetTgl()));
+    static_cast<TH1D*>(fDyHistArray[GetIndex(s1,s2)])
+      ->Fill(tp1.GetY()
+            -tp2.GetY());
+    static_cast<TH1D*>(fDzHistArray[GetIndex(s1,s2)])
+      ->Fill(tp1.GetZ()
+            -tp2.GetZ());
+  }
+  
+
+
+
   //
   // Process function to fill fitters
   //
@@ -70,23 +209,45 @@ void AliTPCcalibAlign::Process(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
@@ -97,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;
@@ -175,7 +336,7 @@ void AliTPCcalibAlign::Process9(Double_t *t1,
   //       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];
+  Double_t /*&x2=t2[0],*/ &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
 
   // TODO:
   Double_t sy    = 1.;
@@ -252,13 +413,13 @@ void AliTPCcalibAlign::Process6(Double_t *t1,
   //       a10  1   0    a13 ==> p[0]   x     x     p[4]
   //       a20  a21 1    a23     p[1]   p[2]  x     p[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];
 
   // 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;
@@ -316,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();