[u/mrichter/AliRoot.git] / TPC / AliTPCkalmanAlign.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
  TPC Kalman filter implementation for TPC sector alignment.
  Output of the AliTPCcalibAlign is used as a input for TPC global alignment.
  In AliTPCcalibAlign histograms - track parameter matching at sector boundaries are created.
  Each sector is aligned with 5 neighborhoud (sectors)
    1. Up-down    - 1
    2. Left-right - 4

  Sector alignment parameters are obtained finding the alignment parameters, minimizing
  misalignmnet for all piars fo sectors.

  Global minimization-  MakeGlobalAlign
  

  Example usage:
  gSystem->Load("libANALYSIS");
  gSystem->Load("libTPCcalib");
  //
  Int_t run=117092;
  .x ConfigCalibTrain.C(run)

  AliTPCkalmanAlign kalmanAlign("TPC align", "TPC align");  // create the object
  kalmanAlign.ReadAlign("CalibObjects.root");               // read the calibration form file         
  kalmanAlign.MakeGlobalAlign();                            // do kalman alignment
  kalmanAlign.DrawDeltaAlign();                             // make QA plot
  //
  

*/
#include "TMath.h"
#include "TTreeStream.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TVectorD.h"
#include "TClonesArray.h"


#include "AliCDBMetaData.h"
#include "AliCDBId.h"
#include "AliCDBManager.h"
#include "AliCDBStorage.h"
#include "AliCDBEntry.h"
#include "AliAlignObjParams.h"
#include "AliTPCROC.h"
#include "TFile.h"
#include "TLinearFitter.h"
#include "AliTPCcalibAlign.h"
#include "TH1.h"
#include "AliTPCCalPad.h"
#include "AliTPCkalmanAlign.h"
#include "TLegend.h"
#include "TCanvas.h"

AliTPCkalmanAlign::AliTPCkalmanAlign():
  TNamed(),
  fCalibAlign(0),     // kalman alignnmnt
  fOriginalAlign(0),   // original alignment 0 read for the OCDB
  fNewAlign(0)
{
  //
  // Default constructor
  //
  for (Int_t i=0; i<4; i++){
    fDelta1D[i]=0;
    fCovar1D[i]=0;
  }
}

AliTPCkalmanAlign::AliTPCkalmanAlign(const char* name, const char* title): 
  TNamed(name,title),
  fCalibAlign(0),     // kalman alignnmnt
  fOriginalAlign(0),   // original alignment 0 read for the OCDB
  fNewAlign(0)     // kalman alignnmnt
{
  //
  // Default constructor
  //
  for (Int_t i=0; i<4; i++){
    fDelta1D[i]=0;
    fCovar1D[i]=0;
  }
}

void AliTPCkalmanAlign::ReadAlign(const char *fname){
  //
  // Read old alignment used in the reco
  // and the residual histograms
  // WE ASSUME that the OCDB path is set in the same way as done in the calibration
  //
  TFile fcalib(fname);
  TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  fCalibAlign=0;
  if (array) fCalibAlign=( AliTPCcalibAlign *)array->FindObject("alignTPC");
  fCalibAlign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
  //
  // old alignment used
  AliCDBEntry * cdbEntry= AliCDBManager::Instance()->Get("TPC/Align/Data");
  fOriginalAlign =0;
  if (cdbEntry){
    fOriginalAlign = (TClonesArray*)cdbEntry->GetObject();
  }
  
}

void AliTPCkalmanAlign::BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t mean, Double_t sigma){
  //
  // Book Align 1D parameters and covariance
  //
  param.ResizeTo(72,1);
  covar.ResizeTo(72,72);
  for (Int_t i=0;i<72;i++){
    param(i,0)=mean;
    for (Int_t j=0;j<72;j++) covar(i,j)=0;      
    covar(i,i)=sigma*sigma;
  }
}


void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &vecXk, TMatrixD &covXk){
  //
  // Update 1D kalman filter
  //
  const Int_t knMeas=1;
  const Int_t knElem=72;
  TMatrixD mat1(72,72);            // update covariance matrix
  TMatrixD matHk(1,knElem);        // vector to mesurement
  TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
  TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
  TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
  TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
  TMatrixD covXk2(knElem,knElem);  // helper matrix
  TMatrixD covXk3(knElem,knElem);  // helper matrix
  TMatrixD vecZk(1,1);
  TMatrixD measR(1,1);
  vecZk(0,0)=delta;
  measR(0,0)=sigma*sigma;
  //
  // reset matHk
  for (Int_t iel=0;iel<knElem;iel++) 
    for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
  //mat1
  for (Int_t iel=0;iel<knElem;iel++) {
    for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
    mat1(iel,iel)=1;
  }
  //
  matHk(0, s1)=1;
  matHk(0, s2)=-1;
  vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
  matHkT=matHk.T(); matHk.T();
  matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
  matSk.Invert();
  matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
  vecXk += matKk*vecYk;                    //  updated vector 
  covXk2= (mat1-(matKk*matHk));
  covXk3 =  covXk2*covXk;          
  covXk = covXk3;  
}


void AliTPCkalmanAlign::MakeGlobalAlign(){
  //
  // Combine all pairs of fitters and make global alignemnt
  //
  AliTPCkalmanAlign &kalmanAlign=*this;
  TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
  DumpOldAlignment(pcstream);
  const Int_t kMinEntries=400;
  TMatrixD vec[5];
  TMatrixD cov[5];
  kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
  kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
  kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
  kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
  //
  TVectorD delta(5);
  TVectorD rms(5);
  TH1 * his=0;
  for (Int_t is0=0;is0<72;is0++)
    for (Int_t is1=0;is1<72;is1++){
      //
      //
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
      if (!his) continue;
      if (his->GetEntries()<kMinEntries) continue;
      delta[0]=his->GetMean();
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
      //     
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
      if (!his) continue;
      delta[1]=his->GetMean();
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
      //
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
      if (!his) continue;
      delta[2] = his->GetMean();
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
      //
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
      if (!his) continue;
      delta[3] = his->GetMean();       
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
    }  
  for (Int_t is0=0;is0<72;is0++)
    for (Int_t is1=0;is1<72;is1++){
      //
      //
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
      if (!his) continue;
      if (his->GetEntries()<kMinEntries) continue;
      delta[0]=his->GetMean();
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
      //     
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
      if (!his) continue;
      delta[1]=his->GetMean();
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
      //
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
      if (!his) continue;
      delta[2] = his->GetMean();
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
      //
      his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
      if (!his) continue;
      delta[3] = his->GetMean();       
      kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
      (*pcstream)<<"kalmanAlignDebug"<<
	"is0="<<is0<<
	"is1="<<is1<<
	"delta.="<<&delta<<
	"vec0.="<<&vec[0]<<
	"vec1.="<<&vec[1]<<
	"vec2.="<<&vec[2]<<
	"vec3.="<<&vec[3]<<
	"\n";
    }
  fDelta1D[0] = (TMatrixD*)vec[0].Clone();
  fDelta1D[1] = (TMatrixD*)vec[1].Clone();
  fDelta1D[2] = (TMatrixD*)vec[2].Clone();
  fDelta1D[3] = (TMatrixD*)vec[3].Clone();
  //
  fCovar1D[0] = (TMatrixD*)cov[0].Clone();
  fCovar1D[1] = (TMatrixD*)cov[1].Clone();
  fCovar1D[2] = (TMatrixD*)cov[2].Clone();
  fCovar1D[3] = (TMatrixD*)cov[3].Clone();
  delete pcstream;
}


void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad  *pad, Int_t ustartRun, Int_t uendRun,  const char* storagePath ){
  //
  // Update OCDB
  // .x ConfigCalibTrain.C(117116)
  // AliTPCcalibDB::Instance()->GetPulserTmean()
  // pad->Add(-pad->GetMean())
  AliCDBMetaData *metaData= new AliCDBMetaData();
  metaData->SetObjectClassName("TObjArray");
  metaData->SetResponsible("Marian Ivanov");
  metaData->SetBeamPeriod(1);
  metaData->SetAliRootVersion("05-25-01"); //root version
  metaData->SetComment("Calibration of the z - time misalignment");
  AliCDBId* id1=NULL;
  id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
  AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
  gStorage->Put(pad, (*id1), metaData);
}


void AliTPCkalmanAlign::DrawDeltaAlign(){
  //
  // Draw the reuslts of the alignment
  // Residual misalignment in respect with previous alignment shown
  //
  //
  TFile f("AliTPCkalmanAlign.root","update");
  TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
  TH1::AddDirectory(0);
  // 
  treeDelta->SetAlias("sector","is0");
  treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
  treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
  treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
  treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
  //
  treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
  treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
  treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
  treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
  //
  treeDelta->SetMarkerStyle(25);
  treeDelta->SetMarkerColor(4);
  treeDelta->SetLineColor(4);
  const char *type[3]={"up-down","left-right","right-left"};
  const char *gropt[3]={"alp","lp","lp"};
  const char *varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10*vec0.fElements[is0]:sector"};
  const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
  const char *varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10*vec1.fElements[is0]:sector"};
  const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
  const char *varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000*vec3.fElements[is0]:sector"};
  const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
  const char *varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000*vec2.fElements[is0]:sector"};
  const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
  TLegend *legend = 0;
  Int_t grcol[3]={2,1,4};
  Int_t entries=0;
  TGraph *grDelta[3]={0,0,0};
  TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
  TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
  TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
  TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
  canvasDy->Divide(2,2);
  canvasDz->Divide(2,2);
  canvasDtheta->Divide(2,2);
  canvasDphi->Divide(2,2);

  //
  // Dy
  //
  canvasDy->cd(1);
  treeDelta->Draw("dYmeas:dYfit");
  for (Int_t itype=0; itype<3; itype++){
    canvasDy->cd(itype+2);
    entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
    grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
    grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
    grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
    for (Int_t i=0; i<3; i++) {
      grDelta[i]->SetMaximum(1.5);
      grDelta[i]->SetMinimum(-1);
      grDelta[i]->SetTitle(type[i]);
      grDelta[i]->SetMarkerColor(grcol[i]); 
      grDelta[i]->SetLineColor(grcol[i]); 
      grDelta[i]->SetMarkerStyle(25+i); 
      grDelta[i]->GetXaxis()->SetTitle("sector"); 
      grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)"); 
      if (itype==2 && i>0) continue;
      grDelta[i]->Draw(gropt[i]); 
      legend->AddEntry(grDelta[i]);
    }
    legend->Draw();
  }
  //
  // Dz
  //
  canvasDz->cd();
  canvasDz->cd(1);
  treeDelta->Draw("dZmeas:dZfit");
  for (Int_t itype=0; itype<3; itype++){
    canvasDz->cd(itype+2);
    entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
    grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
    grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
    grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
    for (Int_t i=0; i<3; i++) {
      grDelta[i]->SetMaximum(1.5);
      grDelta[i]->SetMinimum(-1);
      grDelta[i]->SetTitle(type[i]);
      grDelta[i]->SetMarkerColor(grcol[i]); 
      grDelta[i]->SetLineColor(grcol[i]); 
      grDelta[i]->SetMarkerStyle(25+i); 
      grDelta[i]->GetXaxis()->SetTitle("sector"); 
      grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)"); 
      if (itype==2 && i>0) continue;
      grDelta[i]->Draw(gropt[i]); 
      legend->AddEntry(grDelta[i]);
    }
    legend->Draw();
  }

  //
  // Dtheta
  //
  canvasDtheta->cd(1);
  treeDelta->Draw("dThetameas:dThetafit");
  for (Int_t itype=0; itype<3; itype++){
    canvasDtheta->cd(itype+2);
    entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
    grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
    grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
    grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
    for (Int_t i=0; i<3; i++) {
      grDelta[i]->SetMaximum(4.);
      grDelta[i]->SetMinimum(-3.);
      grDelta[i]->SetTitle(type[i]);
      grDelta[i]->SetMarkerColor(grcol[i]); 
      grDelta[i]->SetLineColor(grcol[i]); 
      grDelta[i]->SetMarkerStyle(25+i); 
      grDelta[i]->GetXaxis()->SetTitle("sector"); 
      grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)"); 
      if (itype==2 && i>0) continue;
      grDelta[i]->Draw(gropt[i]); 
      legend->AddEntry(grDelta[i]);
    }
    legend->Draw();
  }

  //
  // Dphi
  //
  canvasDphi->cd();
  canvasDphi->cd(1);
  treeDelta->Draw("dPhimeas:dPhifit");
  for (Int_t itype=0; itype<3; itype++){
    canvasDphi->cd(itype+2);
    entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
    grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
    grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
    grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
    legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
    for (Int_t i=0; i<3; i++) {
      grDelta[i]->SetMaximum(4.);
      grDelta[i]->SetMinimum(-3.);
      grDelta[i]->SetTitle(type[i]);
      grDelta[i]->SetMarkerColor(grcol[i]); 
      grDelta[i]->SetLineColor(grcol[i]); 
      grDelta[i]->SetMarkerStyle(25+i); 
      grDelta[i]->GetXaxis()->SetTitle("sector"); 
      grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)"); 
      if (itype==2 && i>0) continue;
      grDelta[i]->Draw(gropt[i]); 
      legend->AddEntry(grDelta[i]);
    }
    legend->Draw();
  }
  //
  //
  f.cd();
  canvasDy->Write();
  canvasDz->Write();
  canvasDtheta->Write();
  canvasDphi->Write();
  //
  canvasDy->SaveAs("alignDy.pdf");
  canvasDz->SaveAs("alignDz.pdf");
  canvasDtheta->SaveAs("alignDtheta.pdf");
  canvasDphi->SaveAs("alignDphi.pdf");
}


void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
  // Dump the content of old alignemnt
  // Expected that the old alignmnet is loaded
  //
  if (!fOriginalAlign) return;
  //
  TVectorD localTrans(3);
  TVectorD globalTrans(3);
  TVectorD localRot(3);
  TVectorD globalRot(3);
  AliGeomManager::ELayerID idLayer;
  Int_t idModule=0;
  //
  for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
    AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
    params->GetVolUID(idLayer,idModule);
    params->GetLocalTranslation(localTrans.GetMatrixArray());
    params->GetLocalAngles(localRot.GetMatrixArray());
    params->GetTranslation(globalTrans.GetMatrixArray());
    params->GetAngles(globalRot.GetMatrixArray());
    Int_t sector=idModule;
    if (idLayer>7) sector+=36;
    (*pcstream)<<"oldAlign"<<
      //"idLayer="<<idLayer<<
      "idModule="<<idModule<<
      "sector="<<sector<<
      "lT.="<<&localTrans<<
      "gT.="<<&localTrans<<
      "lR.="<<&localRot<<
      "gR.="<<&globalRot<<
      "\n";
  }
}


void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
  //
  // make a new Alignment entry
  //
  if (!fOriginalAlign) return;
  //
  TVectorD localTrans(3);
  TVectorD globalTrans(3);
  TVectorD localRot(3);
  TVectorD globalRot(3);
  //
  TVectorD localTransNew(3);   // new entries
  TVectorD globalTransNew(3);
  TVectorD localRotNew(3);
  TVectorD globalRotNew(3);
  //
  AliGeomManager::ELayerID idLayer;
  Int_t idModule=0;
  //
  fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
  for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
    AliAlignObjParams *params = (AliAlignObjParams*)fOriginalAlign->At(i);
    AliAlignObjParams *paramsNew = (AliAlignObjParams*)fNewAlign->At(i);
    params->GetVolUID(idLayer,idModule);
    Int_t sector=(Int_t)idModule;
    if (idLayer>7) sector+=36;
    params->GetLocalTranslation(localTrans.GetMatrixArray());
    params->GetLocalAngles(localRot.GetMatrixArray());
    params->GetTranslation(globalTrans.GetMatrixArray());
    params->GetAngles(globalRot.GetMatrixArray());
    //
    //
    //
    if (badd){ // addition if 
      localTransNew=localTrans;
      localRotNew=localRot;
    }
//     localTrans[1]=localTrans[1]-(*fDelta1D[0])[sector];
//     localRot[0]  =localRot[0]-(*fDelta1D[0])[sector];
    //
    if (pcstream) (*pcstream)<<"newAlign"<<
      //"idLayer="<<idLayer<<
      "idModule="<<idModule<<
      "sector="<<sector<<
      "olT.="<<&localTrans<<
      "ogT.="<<&localTrans<<
      "olR.="<<&localRot<<
      "ogR.="<<&globalRot<<
      "\n";
  }
  

}
Commit	Line	Data
125d3a38	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17	TPC Kalman filter implementation for TPC sector alignment.
	18	Output of the AliTPCcalibAlign is used as a input for TPC global alignment.
	19	In AliTPCcalibAlign histograms - track parameter matching at sector boundaries are created.
	20	Each sector is aligned with 5 neighborhoud (sectors)
	21	1. Up-down - 1
	22	2. Left-right - 4
	23
	24	Sector alignment parameters are obtained finding the alignment parameters, minimizing
	25	misalignmnet for all piars fo sectors.
	26
	27	Global minimization- MakeGlobalAlign
	28
	29
	30	Example usage:
	31	gSystem->Load("libANALYSIS");
	32	gSystem->Load("libTPCcalib");
	33	//
	34	Int_t run=117092;
	35	.x ConfigCalibTrain.C(run)
	36
	37	AliTPCkalmanAlign kalmanAlign("TPC align", "TPC align"); // create the object
	38	kalmanAlign.ReadAlign("CalibObjects.root"); // read the calibration form file
	39	kalmanAlign.MakeGlobalAlign(); // do kalman alignment
	40	kalmanAlign.DrawDeltaAlign(); // make QA plot
	41	//
	42
	43
	44	*/
	45	#include "TMath.h"
	46	#include "TTreeStream.h"
	47	#include "TGraph.h"
	48	#include "TGraphErrors.h"
	49	#include "TVectorD.h"
	50	#include "TClonesArray.h"
	51
	52
	53	#include "AliCDBMetaData.h"
	54	#include "AliCDBId.h"
	55	#include "AliCDBManager.h"
	56	#include "AliCDBStorage.h"
	57	#include "AliCDBEntry.h"
	58	#include "AliAlignObjParams.h"
	59	#include "AliTPCROC.h"
	60	#include "TFile.h"
	61	#include "TLinearFitter.h"
	62	#include "AliTPCcalibAlign.h"
	63	#include "TH1.h"
	64	#include "AliTPCCalPad.h"
65	#include "AliTPCkalmanAlign.h"
66	#include "TLegend.h"
67	#include "TCanvas.h"
68
69	AliTPCkalmanAlign::AliTPCkalmanAlign():
70	TNamed(),
71	fCalibAlign(0), // kalman alignnmnt
72	fOriginalAlign(0), // original alignment 0 read for the OCDB
73	fNewAlign(0)
74	{
75	//
76	// Default constructor
77	//
78	for (Int_t i=0; i<4; i++){
79	fDelta1D[i]=0;
80	fCovar1D[i]=0;
81	}
82	}
83
84	AliTPCkalmanAlign::AliTPCkalmanAlign(const char* name, const char* title):
85	TNamed(name,title),
86	fCalibAlign(0), // kalman alignnmnt
87	fOriginalAlign(0), // original alignment 0 read for the OCDB
88	fNewAlign(0) // kalman alignnmnt
89	{
90	//
91	// Default constructor
92	//
93	for (Int_t i=0; i<4; i++){
94	fDelta1D[i]=0;
95	fCovar1D[i]=0;
96	}
97	}
98
99	void AliTPCkalmanAlign::ReadAlign(const char *fname){
100	//
101	// Read old alignment used in the reco
102	// and the residual histograms
103	// WE ASSUME that the OCDB path is set in the same way as done in the calibration
104	//
105	TFile fcalib(fname);
106	TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
107	fCalibAlign=0;
108	if (array) fCalibAlign=( AliTPCcalibAlign *)array->FindObject("alignTPC");
109	fCalibAlign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
110	//
111	// old alignment used
112	AliCDBEntry * cdbEntry= AliCDBManager::Instance()->Get("TPC/Align/Data");
113	fOriginalAlign =0;
114	if (cdbEntry){
115	fOriginalAlign = (TClonesArray*)cdbEntry->GetObject();
116	}
117
118	}
119
120	void AliTPCkalmanAlign::BookAlign1D(TMatrixD &param, TMatrixD &covar, Double_t mean, Double_t sigma){
121	//
122	// Book Align 1D parameters and covariance
123	//
124	param.ResizeTo(72,1);
125	covar.ResizeTo(72,72);
126	for (Int_t i=0;i<72;i++){
127	param(i,0)=mean;
128	for (Int_t j=0;j<72;j++) covar(i,j)=0;
129	covar(i,i)=sigma*sigma;
130	}
131	}
132
133
134	void AliTPCkalmanAlign::UpdateAlign1D(Double_t delta, Double_t sigma, Int_t s1, Int_t s2, TMatrixD &vecXk, TMatrixD &covXk){
135	//
136	// Update 1D kalman filter
137	//
138	const Int_t knMeas=1;
139	const Int_t knElem=72;
140	TMatrixD mat1(72,72); // update covariance matrix
141	TMatrixD matHk(1,knElem); // vector to mesurement
142	TMatrixD vecYk(knMeas,1); // Innovation or measurement residual
143	TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose
144	TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance
145	TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain
146	TMatrixD covXk2(knElem,knElem); // helper matrix
147	TMatrixD covXk3(knElem,knElem); // helper matrix
148	TMatrixD vecZk(1,1);
149	TMatrixD measR(1,1);
150	vecZk(0,0)=delta;
151	measR(0,0)=sigma*sigma;
152	//
153	// reset matHk
154	for (Int_t iel=0;iel<knElem;iel++)
155	for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
156	//mat1
157	for (Int_t iel=0;iel<knElem;iel++) {
158	for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
159	mat1(iel,iel)=1;
160	}
161	//
162	matHk(0, s1)=1;
163	matHk(0, s2)=-1;
164	vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
165	matHkT=matHk.T(); matHk.T();
166	matSk = (matHk(covXkmatHkT))+measR; // Innovation (or residual) covariance
167	matSk.Invert();
168	matKk = (covXkmatHkT)matSk; // Optimal Kalman gain
169	vecXk += matKk*vecYk; // updated vector
170	covXk2= (mat1-(matKk*matHk));
171	covXk3 = covXk2*covXk;
172	covXk = covXk3;
173	}
174
175
176
177
178	void AliTPCkalmanAlign::MakeGlobalAlign(){
179	//
180	// Combine all pairs of fitters and make global alignemnt
181	//
182	AliTPCkalmanAlign &kalmanAlign=*this;
183	TTreeSRedirector *pcstream = new TTreeSRedirector("AliTPCkalmanAlign.root");
184	DumpOldAlignment(pcstream);
185	const Int_t kMinEntries=400;
186	TMatrixD vec[5];
187	TMatrixD cov[5];
188	kalmanAlign.BookAlign1D(vec[0],cov[0], 0,0.5);
189	kalmanAlign.BookAlign1D(vec[1],cov[1], 0,0.5);
190	kalmanAlign.BookAlign1D(vec[2],cov[2], 0,0.01);
191	kalmanAlign.BookAlign1D(vec[3],cov[3], 0,0.01);
192	//
193	TVectorD delta(5);
194	TVectorD rms(5);
195	TH1 * his=0;
196	for (Int_t is0=0;is0<72;is0++)
197	for (Int_t is1=0;is1<72;is1++){
198	//
199	//
200	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
201	if (!his) continue;
202	if (his->GetEntries()<kMinEntries) continue;
203	delta[0]=his->GetMean();
204	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
205	//
206	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
207	if (!his) continue;
208	delta[1]=his->GetMean();
209	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
210	//
211	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
212	if (!his) continue;
213	delta[2] = his->GetMean();
214	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
215	//
216	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
217	if (!his) continue;
218	delta[3] = his->GetMean();
219	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
220	}
221	for (Int_t is0=0;is0<72;is0++)
222	for (Int_t is1=0;is1<72;is1++){
223	//
224	//
225	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kY,is0,is1);
226	if (!his) continue;
227	if (his->GetEntries()<kMinEntries) continue;
228	delta[0]=his->GetMean();
229	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[0],cov[0]);
230	//
231	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kZ,is0,is1);
232	if (!his) continue;
233	delta[1]=his->GetMean();
234	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[1],cov[1]);
235	//
236	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kPhi,is0,is1);
237	if (!his) continue;
238	delta[2] = his->GetMean();
239	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[2],cov[2]);
240	//
241	his = kalmanAlign.fCalibAlign->GetHisto(AliTPCcalibAlign::kTheta,is0,is1);
242	if (!his) continue;
243	delta[3] = his->GetMean();
244	kalmanAlign.UpdateAlign1D(his->GetMean(),his->GetRMS(),is0,is1, vec[3],cov[3]);
245	(*pcstream)<<"kalmanAlignDebug"<<
246	"is0="<<is0<<
247	"is1="<<is1<<
248	"delta.="<<&delta<<
249	"vec0.="<<&vec[0]<<
250	"vec1.="<<&vec[1]<<
251	"vec2.="<<&vec[2]<<
252	"vec3.="<<&vec[3]<<
253	"\n";
254	}
255	fDelta1D[0] = (TMatrixD*)vec[0].Clone();
256	fDelta1D[1] = (TMatrixD*)vec[1].Clone();
257	fDelta1D[2] = (TMatrixD*)vec[2].Clone();
258	fDelta1D[3] = (TMatrixD*)vec[3].Clone();
259	//
260	fCovar1D[0] = (TMatrixD*)cov[0].Clone();
261	fCovar1D[1] = (TMatrixD*)cov[1].Clone();
262	fCovar1D[2] = (TMatrixD*)cov[2].Clone();
263	fCovar1D[3] = (TMatrixD*)cov[3].Clone();
264	delete pcstream;
265	}
266
267
268
269
270
271
272	void AliTPCkalmanAlign::UpdateOCDBTime0( AliTPCCalPad pad, Int_t ustartRun, Int_t uendRun, const char storagePath ){
273	//
274	// Update OCDB
275	// .x ConfigCalibTrain.C(117116)
276	// AliTPCcalibDB::Instance()->GetPulserTmean()
277	// pad->Add(-pad->GetMean())
278	AliCDBMetaData *metaData= new AliCDBMetaData();
279	metaData->SetObjectClassName("TObjArray");
280	metaData->SetResponsible("Marian Ivanov");
281	metaData->SetBeamPeriod(1);
282	metaData->SetAliRootVersion("05-25-01"); //root version
283	metaData->SetComment("Calibration of the z - time misalignment");
284	AliCDBId* id1=NULL;
285	id1=new AliCDBId("TPC/Calib/PadTime0", ustartRun, uendRun);
286	AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage(storagePath);
287	gStorage->Put(pad, (*id1), metaData);
288	}
289
290
291
292	void AliTPCkalmanAlign::DrawDeltaAlign(){
293	//
294	// Draw the reuslts of the alignment
295	// Residual misalignment in respect with previous alignment shown
296	//
297	//
298	TFile f("AliTPCkalmanAlign.root","update");
299	TTree * treeDelta = (TTree*)f.Get("kalmanAlignDebug");
300	TH1::AddDirectory(0);
301	//
302	treeDelta->SetAlias("sector","is0");
303	treeDelta->SetAlias("dYmeas","10*(delta.fElements[0])");
304	treeDelta->SetAlias("dZmeas","10*(delta.fElements[1])");
305	treeDelta->SetAlias("dPhimeas","1000*(delta.fElements[2])");
306	treeDelta->SetAlias("dThetameas","1000*(delta.fElements[3])");
307	//
308	treeDelta->SetAlias("dYfit","10*(vec0.fElements[is0]-vec0.fElements[is1])");
309	treeDelta->SetAlias("dZfit","10*(vec1.fElements[is0]-vec1.fElements[is1])");
310	treeDelta->SetAlias("dPhifit","1000*(vec2.fElements[is0]-vec2.fElements[is1])");
311	treeDelta->SetAlias("dThetafit","1000*(vec3.fElements[is0]-vec3.fElements[is1])");
312	//
313	treeDelta->SetMarkerStyle(25);
314	treeDelta->SetMarkerColor(4);
315	treeDelta->SetLineColor(4);
316	const char *type[3]={"up-down","left-right","right-left"};
317	const char *gropt[3]={"alp","lp","lp"};
318	const char varTypeY[3]={"dYmeas-dYfit:sector","dYfit:sector","10vec0.fElements[is0]:sector"};
319	const char *varLegendY[3]={"#Delta_{y} fit residual","#Delta_{y} fit value difference","#Delta_{y} sector"};
320	const char varTypeZ[3]={"dZmeas-dZfit:sector","dZfit:sector","10vec1.fElements[is0]:sector"};
321	const char *varLegendZ[3]={"#Delta_{Z} fit residual","#Delta_{Z} fit value difference","#Delta_{Z} sector"};
322	const char varTypeT[3]={"dThetameas-dThetafit:sector","dThetafit:sector","1000vec3.fElements[is0]:sector"};
323	const char *varLegendT[3]={"#Delta_{#theta} fit residual","#Delta_{#theta} fit value difference","#Delta_{#theta} sector"};
324	const char varTypeP[3]={"dPhimeas-dPhifit:sector","dPhifit:sector","1000vec2.fElements[is0]:sector"};
325	const char *varLegendP[3]={"#Delta_{#phi} fit residual","#Delta_{#phi} fit value difference","#Delta_{#phi} sector"};
326	TLegend *legend = 0;
327	Int_t grcol[3]={2,1,4};
328	Int_t entries=0;
329	TGraph *grDelta[3]={0,0,0};
330	TCanvas * canvasDy=new TCanvas("canvasDy","canvasDy",1200,800);
331	TCanvas * canvasDz=new TCanvas("canvasDz","canvasDz",1200,800);
332	TCanvas * canvasDphi=new TCanvas("canvasDphi","canvasDphi",1200,800);
333	TCanvas * canvasDtheta=new TCanvas("canvasDtheta","canvasDtheta",1200,800);
334	canvasDy->Divide(2,2);
335	canvasDz->Divide(2,2);
336	canvasDtheta->Divide(2,2);
337	canvasDphi->Divide(2,2);
338
339	//
340	// Dy
341	//
342	canvasDy->cd(1);
343	treeDelta->Draw("dYmeas:dYfit");
344	for (Int_t itype=0; itype<3; itype++){
345	canvasDy->cd(itype+2);
346	entries=treeDelta->Draw(varTypeY[itype],"is1==is0+36","goff");
347	grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
348	entries=treeDelta->Draw(varTypeY[itype],"is1==is0+35","goff");
349	grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
350	entries=treeDelta->Draw(varTypeY[itype],"is1==is0+37","goff");
351	grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
352	legend = new TLegend(0.5,0.7,0.9,0.9, varLegendY[itype]);
353	for (Int_t i=0; i<3; i++) {
354	grDelta[i]->SetMaximum(1.5);
355	grDelta[i]->SetMinimum(-1);
356	grDelta[i]->SetTitle(type[i]);
357	grDelta[i]->SetMarkerColor(grcol[i]);
358	grDelta[i]->SetLineColor(grcol[i]);
359	grDelta[i]->SetMarkerStyle(25+i);
360	grDelta[i]->GetXaxis()->SetTitle("sector");
361	grDelta[i]->GetYaxis()->SetTitle("#Delta_{y} (mm)");
362	if (itype==2 && i>0) continue;
363	grDelta[i]->Draw(gropt[i]);
364	legend->AddEntry(grDelta[i]);
365	}
366	legend->Draw();
367	}
368	//
369	// Dz
370	//
371	canvasDz->cd();
372	canvasDz->cd(1);
373	treeDelta->Draw("dZmeas:dZfit");
374	for (Int_t itype=0; itype<3; itype++){
375	canvasDz->cd(itype+2);
376	entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+36","goff");
377	grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
378	entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+35","goff");
379	grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
380	entries=treeDelta->Draw(varTypeZ[itype],"is1==is0+37","goff");
381	grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
382	legend = new TLegend(0.5,0.7,0.9,0.9, varLegendZ[itype]);
383	for (Int_t i=0; i<3; i++) {
384	grDelta[i]->SetMaximum(1.5);
385	grDelta[i]->SetMinimum(-1);
386	grDelta[i]->SetTitle(type[i]);
387	grDelta[i]->SetMarkerColor(grcol[i]);
388	grDelta[i]->SetLineColor(grcol[i]);
389	grDelta[i]->SetMarkerStyle(25+i);
390	grDelta[i]->GetXaxis()->SetTitle("sector");
391	grDelta[i]->GetYaxis()->SetTitle("#Delta_{z} (mm)");
392	if (itype==2 && i>0) continue;
393	grDelta[i]->Draw(gropt[i]);
394	legend->AddEntry(grDelta[i]);
395	}
396	legend->Draw();
397	}
398
399	//
400	// Dtheta
401	//
402	canvasDtheta->cd(1);
403	treeDelta->Draw("dThetameas:dThetafit");
404	for (Int_t itype=0; itype<3; itype++){
405	canvasDtheta->cd(itype+2);
406	entries=treeDelta->Draw(varTypeT[itype],"is1==is0+36","goff");
407	grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
408	entries=treeDelta->Draw(varTypeT[itype],"is1==is0+35","goff");
409	grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
410	entries=treeDelta->Draw(varTypeT[itype],"is1==is0+37","goff");
411	grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
412	legend = new TLegend(0.5,0.7,0.9,0.9, varLegendT[itype]);
413	for (Int_t i=0; i<3; i++) {
414	grDelta[i]->SetMaximum(4.);
415	grDelta[i]->SetMinimum(-3.);
416	grDelta[i]->SetTitle(type[i]);
417	grDelta[i]->SetMarkerColor(grcol[i]);
418	grDelta[i]->SetLineColor(grcol[i]);
419	grDelta[i]->SetMarkerStyle(25+i);
420	grDelta[i]->GetXaxis()->SetTitle("sector");
421	grDelta[i]->GetYaxis()->SetTitle("#Delta_{#theta} (mrad)");
422	if (itype==2 && i>0) continue;
423	grDelta[i]->Draw(gropt[i]);
424	legend->AddEntry(grDelta[i]);
425	}
426	legend->Draw();
427	}
428
429	//
430	// Dphi
431	//
432	canvasDphi->cd();
433	canvasDphi->cd(1);
434	treeDelta->Draw("dPhimeas:dPhifit");
435	for (Int_t itype=0; itype<3; itype++){
436	canvasDphi->cd(itype+2);
437	entries=treeDelta->Draw(varTypeP[itype],"is1==is0+36","goff");
438	grDelta[0]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
439	entries=treeDelta->Draw(varTypeP[itype],"is1==is0+35","goff");
440	grDelta[1]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
441	entries=treeDelta->Draw(varTypeP[itype],"is1==is0+37","goff");
442	grDelta[2]=new TGraph(entries,treeDelta->GetV2(),treeDelta->GetV1());
443	legend = new TLegend(0.5,0.7,0.9,0.9, varLegendP[itype]);
444	for (Int_t i=0; i<3; i++) {
445	grDelta[i]->SetMaximum(4.);
446	grDelta[i]->SetMinimum(-3.);
447	grDelta[i]->SetTitle(type[i]);
448	grDelta[i]->SetMarkerColor(grcol[i]);
449	grDelta[i]->SetLineColor(grcol[i]);
450	grDelta[i]->SetMarkerStyle(25+i);
451	grDelta[i]->GetXaxis()->SetTitle("sector");
452	grDelta[i]->GetYaxis()->SetTitle("#Delta_{#phi} (mrad)");
453	if (itype==2 && i>0) continue;
454	grDelta[i]->Draw(gropt[i]);
455	legend->AddEntry(grDelta[i]);
456	}
457	legend->Draw();
458	}
459	//
460	//
461	f.cd();
462	canvasDy->Write();
463	canvasDz->Write();
464	canvasDtheta->Write();
465	canvasDphi->Write();
466	//
467	canvasDy->SaveAs("alignDy.pdf");
468	canvasDz->SaveAs("alignDz.pdf");
469	canvasDtheta->SaveAs("alignDtheta.pdf");
470	canvasDphi->SaveAs("alignDphi.pdf");
471	}
472
473
474
475	void AliTPCkalmanAlign::DumpOldAlignment(TTreeSRedirector *pcstream){
476	// Dump the content of old alignemnt
477	// Expected that the old alignmnet is loaded
478	//
479	if (!fOriginalAlign) return;
480	//
481	TVectorD localTrans(3);
482	TVectorD globalTrans(3);
483	TVectorD localRot(3);
484	TVectorD globalRot(3);
485	AliGeomManager::ELayerID idLayer;
486	Int_t idModule=0;
487	//
488	for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
489	AliAlignObjParams params = (AliAlignObjParams)fOriginalAlign->At(i);
490	params->GetVolUID(idLayer,idModule);
491	params->GetLocalTranslation(localTrans.GetMatrixArray());
492	params->GetLocalAngles(localRot.GetMatrixArray());
493	params->GetTranslation(globalTrans.GetMatrixArray());
494	params->GetAngles(globalRot.GetMatrixArray());
495	Int_t sector=idModule;
496	if (idLayer>7) sector+=36;
497	(*pcstream)<<"oldAlign"<<
498	//"idLayer="<<idLayer<<
499	"idModule="<<idModule<<
500	"sector="<<sector<<
501	"lT.="<<&localTrans<<
502	"gT.="<<&localTrans<<
503	"lR.="<<&localRot<<
504	"gR.="<<&globalRot<<
505	"\n";
506	}
507	}
508
509
510	void AliTPCkalmanAlign::MakeNewAlignment(Bool_t badd, TTreeSRedirector * pcstream){
511	//
512	// make a new Alignment entry
513	//
514	if (!fOriginalAlign) return;
515	//
516	TVectorD localTrans(3);
517	TVectorD globalTrans(3);
518	TVectorD localRot(3);
519	TVectorD globalRot(3);
520	//
521	TVectorD localTransNew(3); // new entries
522	TVectorD globalTransNew(3);
523	TVectorD localRotNew(3);
524	TVectorD globalRotNew(3);
525	//
526	AliGeomManager::ELayerID idLayer;
527	Int_t idModule=0;
528	//
529	fNewAlign = (TClonesArray*)fOriginalAlign->Clone();
530	for (Int_t i=0; i<fOriginalAlign->GetEntries();i++){
531	AliAlignObjParams params = (AliAlignObjParams)fOriginalAlign->At(i);
532	AliAlignObjParams paramsNew = (AliAlignObjParams)fNewAlign->At(i);
533	params->GetVolUID(idLayer,idModule);
534	Int_t sector=(Int_t)idModule;
535	if (idLayer>7) sector+=36;
536	params->GetLocalTranslation(localTrans.GetMatrixArray());
537	params->GetLocalAngles(localRot.GetMatrixArray());
538	params->GetTranslation(globalTrans.GetMatrixArray());
539	params->GetAngles(globalRot.GetMatrixArray());
540	//
541	//
542	//
543	if (badd){ // addition if
544	localTransNew=localTrans;
545	localRotNew=localRot;
546	}
547	// localTrans[1]=localTrans[1]-(*fDelta1D[0])[sector];
548	// localRot[0] =localRot[0]-(*fDelta1D[0])[sector];
549	//
550	if (pcstream) (*pcstream)<<"newAlign"<<
551	//"idLayer="<<idLayer<<
552	"idModule="<<idModule<<
553	"sector="<<sector<<
554	"olT.="<<&localTrans<<
555	"ogT.="<<&localTrans<<
556	"olR.="<<&localRot<<
557	"ogR.="<<&globalRot<<
558	"\n";
559	}
560
561
562	}