Update for extraction of 3D cluster residual maps
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Jul 2010 10:50:58 +0000 (10:50 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Jul 2010 10:50:58 +0000 (10:50 +0000)
                  - A lot of code in the AliTPCcaliAlign started to be obsolete
                  - To be cleaned later
+
MakeAlignCalPad.C - for E field distortion in 3D
                  - previously working in 2D
                  - name is obsolete - will be replaced in future

TPC/AliTPCcalibAlign.cxx
TPC/AliTPCcalibAlign.h
TPC/CalibMacros/MakeAlignCalPad.C [new file with mode: 0644]

index 5a44932..0486c8e 100644 (file)
 #include "TFile.h"
 #include "TProfile.h"
 #include "TCanvas.h"
-
+#include "TDatabasePDG.h"
 
 #include "TTreeStream.h"
 #include "Riostream.h"
@@ -206,12 +206,8 @@ AliTPCcalibAlign::AliTPCcalibAlign()
   fXquadrant = roc->GetPadRowRadii(36,53);
   fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
-  fClusterDelta[0]=0;   // cluster residuals
-  fClusterDelta[1]=0;   // cluster residuals
-  fClusterDelta[2]=0;   // cluster residuals - vertex constrained
-  fClusterDelta[3]=0;   // cluster residuals
-  fClusterDelta[4]=0;   // cluster residuals - ITS constrained
-  fClusterDelta[5]=0;   // cluster residuals
+  fClusterDelta[0]=0;   // cluster residuals -  Y
+  fClusterDelta[1]=0;   // cluster residuals -  Z
   
   
   fTrackletDelta[0]=0;   // tracklet residuals
@@ -269,10 +265,6 @@ AliTPCcalibAlign::AliTPCcalibAlign(const Text_t *name, const Text_t *title)
   fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
   fClusterDelta[0]=0;   // cluster residuals
   fClusterDelta[1]=0;   // cluster residuals
-  fClusterDelta[2]=0;   // cluster residuals - vertex constrained
-  fClusterDelta[3]=0;   // cluster residuals
-  fClusterDelta[4]=0;   // cluster residuals - ITS constrained
-  fClusterDelta[5]=0;   // cluster residuals
 
   fTrackletDelta[0]=0;   // tracklet residuals
   fTrackletDelta[1]=0;   // tracklet residuals
@@ -383,10 +375,6 @@ AliTPCcalibAlign::AliTPCcalibAlign(const AliTPCcalibAlign &align)
   }
   fClusterDelta[0]=0;   // cluster residuals
   fClusterDelta[1]=0;   // cluster residuals
-  fClusterDelta[2]=0;   // cluster residuals - vertex constrained
-  fClusterDelta[3]=0;   // cluster residuals
-  fClusterDelta[4]=0;   // cluster residuals - ITS constrained
-  fClusterDelta[5]=0;   // cluster residuals
 
   fTrackletDelta[0]=0;   // tracklet residuals
   fTrackletDelta[1]=0;   // tracklet residuals
@@ -446,7 +434,7 @@ AliTPCcalibAlign::~AliTPCcalibAlign() {
   fArraySectorIntCovar.SetOwner(kTRUE); // array of sector alignment covariances 
   fArraySectorIntParam.Delete(); // array of sector alignment parameters
   fArraySectorIntCovar.Delete(); // array of sector alignment covariances 
-  for (Int_t i=0; i<6; i++){
+  for (Int_t i=0; i<2; i++){
     delete fClusterDelta[i];   // cluster residuals
   }
 
@@ -490,6 +478,7 @@ void AliTPCcalibAlign::Process(AliESDEvent *event) {
     }
     if (!seed0) continue;
     fCurrentTrack=track0;
+    fCurrentFriendTrack=friendTrack;
     fCurrentSeed=seed0;
     fCurrentEvent=event;
     ProcessSeed(seed0);
@@ -798,6 +787,7 @@ void AliTPCcalibAlign::ProcessSeed(AliTPCseed *seed) {
   //
   // make a kalman tracklets out of seed
   //
+  UpdateClusterDeltaField(seed);
   TObjArray tracklets=
     AliTPCTracklet::CreateTracklets(seed,AliTPCTracklet::kKalman,
                                    kFALSE,20,4);
@@ -921,10 +911,6 @@ void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
   //
   Int_t accept       =   AcceptTracklet(tp1,tp2);  
   Int_t acceptLinear =   AcceptTracklet(parLine1,parLine2);
-  if (accept==0){
-    FillHisto(&tp1,&tp2, s1,s2);
-    FillHisto(&tp2,&tp1, s2,s1);
-  }
 
 
   if (fStreamLevel>1 && seed){
@@ -983,7 +969,11 @@ void AliTPCcalibAlign::ProcessTracklets(const AliExternalTrackParam &tp1,
     //
     // use Kalman if mag field
     //
-    if (seed) ProcessDiff(tp1,tp2, seed,s1,s2);
+    if (seed) {
+      ProcessDiff(tp1,tp2, seed,s1,s2);
+      FillHisto((AliExternalTrackParam*)&tp1,(AliExternalTrackParam*)&tp2,s1,s2);
+      FillHisto((AliExternalTrackParam*)&tp2,(AliExternalTrackParam*)&tp1,s2,s1);
+    }
     FillHisto(t1,t2,s1,s2);  
     ProcessAlign(t1,t2,s1,s2);
   }
@@ -1689,23 +1679,17 @@ void AliTPCcalibAlign::MakeResidualHistos(){
   axisName[2]="localX";   axisTitle[2]="x (cm)"; 
   binsTrack[2]=53;       xminTrack[2]=85.;        xmaxTrack[2]=245.; 
   //
-  axisName[3]="kY";      axisTitle[3]="dy/dx"; 
-  binsTrack[3]=1;       xminTrack[3]=-0.16;        xmaxTrack[3]=0.16; 
   //
-  axisName[4]="kZ";      axisTitle[4]="dz/dx"; 
-  binsTrack[4]=22;       xminTrack[4]=-1.1;        xmaxTrack[4]=1.1; 
+  axisName[3]="kZ";      axisTitle[3]="dz/dx"; 
+  binsTrack[3]=36;       xminTrack[3]=-1.8;        xmaxTrack[3]=1.8; 
   //
-  fClusterDelta[0] = new THnSparseF("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 5, binsTrack,xminTrack, xmaxTrack);
-  fClusterDelta[1] = new THnSparseF("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 5, binsTrack,xminTrack, xmaxTrack);
-  fClusterDelta[2] = new THnSparseF("#Delta_{Y} (cm) const","#Delta_{Y} (cm) const ", 5, binsTrack,xminTrack, xmaxTrack);
-  fClusterDelta[3] = new THnSparseF("#Delta_{Z} (cm) const","#Delta_{Z} (cm) const", 5, binsTrack,xminTrack, xmaxTrack);
-  fClusterDelta[4] = new THnSparseF("#Delta_{Y} (cm) ITS","#Delta_{Y} (cm) ITS", 5, binsTrack,xminTrack, xmaxTrack);
-  fClusterDelta[5] = new THnSparseF("#Delta_{Z} (cm) ITS","#Delta_{Z} (cm) ITS", 5, binsTrack,xminTrack, xmaxTrack);
+  fClusterDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
+  fClusterDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 4, binsTrack,xminTrack, xmaxTrack);
   //
   //
   //
-  for (Int_t ivar=0;ivar<6;ivar++){
-    for (Int_t ivar2=0;ivar2<5;ivar2++){
+  for (Int_t ivar=0;ivar<2;ivar++){
+    for (Int_t ivar2=0;ivar2<4;ivar2++){
       fClusterDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
       fClusterDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
     }
@@ -1806,21 +1790,23 @@ void AliTPCcalibAlign::FillHisto(const Double_t *t1,
 }
 
 
-void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam *tp1,
-                                const AliExternalTrackParam *tp2,
+void AliTPCcalibAlign::FillHisto(AliExternalTrackParam *tp1,
+                                AliExternalTrackParam *tp2,
                                 Int_t s1,Int_t s2) {
   //
   // Fill residual histograms
   // Track2-Track1
+  const Double_t kEpsilon=0.001;
   Double_t x[8]={0,0,0,0,0,0,0,0};
   AliExternalTrackParam p1(*tp1);
   AliExternalTrackParam p2(*tp2);
   if (s1%18==s2%18) {
     // inner outer - match at the IROC-OROC boundary
-    p1.PropagateTo(fXIO, AliTrackerBase::GetBz());
+    if (!p1.PropagateTo(fXIO, AliTrackerBase::GetBz())) return;
   }
-  p2.Rotate(p1.GetAlpha());
-  p2.PropagateTo(p1.GetX(),AliTrackerBase::GetBz());
+  if (!p2.Rotate(p1.GetAlpha())) return;
+  if (!p2.PropagateTo(p1.GetX(),AliTrackerBase::GetBz())) return;
+  if (TMath::Abs(p1.GetX()-p2.GetX())>kEpsilon) return;
   Double_t xyz[3];
   p1.GetXYZ(xyz);
   x[1]=TMath::ATan2(xyz[1],xyz[0]);
@@ -1838,7 +1824,18 @@ void AliTPCcalibAlign::FillHisto(const AliExternalTrackParam *tp1,
   fTrackletDelta[2]->Fill(x);
   x[0]=p2.GetTgl()-p1.GetTgl();
   fTrackletDelta[3]->Fill(x);
-
+  TTreeSRedirector *cstream = GetDebugStreamer();    
+  if (cstream){
+    (*cstream)<<"trackletMatch"<<
+      "tp1.="<<tp1<<   // input tracklet
+      "tp2.="<<tp2<<
+      "p1.="<<&p1<<    // tracklet in the ref frame
+      "p2.="<<&p2<<
+      "s1="<<s1<<
+      "s2="<<s2<<
+      "\n";
+  }      
+  
 }
 
 
@@ -2216,17 +2213,9 @@ void AliTPCcalibAlign::Add(AliTPCcalibAlign * align){
     UpdateKalman(*fSectorParamA,*fSectorCovarA,*align->fSectorParamA,*align->fSectorCovarA);
     UpdateKalman(*fSectorParamC,*fSectorCovarC,*align->fSectorParamC,*align->fSectorCovarC);
   }
-  if (!fClusterDelta[1]) MakeResidualHistos();
+  if (!fClusterDelta[0]) MakeResidualHistos();
 
-  for (Int_t i=0; i<6; i++){
-    if (i==0 || i==3){
-      delete fClusterDelta[i];   // memory problem do not fit into memory
-      fClusterDelta[i]=0;        // 
-      delete align->fClusterDelta[i];   // memory problem do not fit into memory
-      align->fClusterDelta[i]=0;        // 
-    }
-    if (i==3) continue;  // skip constrained histo z
-    if (i==0) continue;  // skip non constrained histo y
+  for (Int_t i=0; i<2; i++){
     if (align->fClusterDelta[i]) fClusterDelta[i]->Add(align->fClusterDelta[i]);
   }
 
@@ -2640,13 +2629,145 @@ void AliTPCcalibAlign::GlobalAlign6(Int_t minPoints, Float_t sysError, Int_t nit
   return nf;
 }
 
+void AliTPCcalibAlign::UpdateClusterDeltaField(const AliTPCseed * seed){
+  //
+  // Update the cluster residula histograms for setup with field
+  // Kalman track fitting is used
+  //  Only high momenta primary tracks used
+  //
+  // 1. Apply selection
+  // 2. Refit the track - in-out
+  //                    - update the cluster delta in upper part
+  // 3. Refit the track - out-in
+  //                    - update the cluster delta histo lower part
+  //
+  const Double_t kPtCut=1;    // pt
+  const Double_t kSnpCut=0.2; // snp cut
+  const Double_t kNclCut=120; //
+  const Double_t kVertexCut=1;
+  const Double_t kMaxDist=0.5; // max distance between tracks and cluster
+  if (!fCurrentTrack) return;
+  if (!fCurrentFriendTrack) return;
+  Float_t vertexXY=0,vertexZ=0;
+  fCurrentTrack->GetImpactParameters(vertexXY,vertexZ);
+  if (TMath::Abs(vertexXY)>kVertexCut) return;
+  if (TMath::Abs(vertexZ)>kVertexCut) return;
+  if (TMath::Abs(seed->Pt())<kPtCut) return;
+  if (seed->GetNumberOfClusters()<kNclCut) return;
+  if (TMath::Abs(seed->GetSnp())>kSnpCut) return; 
+  if (!fClusterDelta[0])  MakeResidualHistos();
+
+  Int_t detector=-1;
+  //
+  //
+  AliExternalTrackParam trackIn  = *(fCurrentTrack->GetInnerParam());
+  AliExternalTrackParam trackOut = *(fCurrentFriendTrack->GetTPCOut());
+  Double_t mass =    TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  //  
+  Int_t nclIn=0,nclOut=0;
+  Double_t xyz[3];
+  //
+  // Refit out - store residual maps
+  //
+  for (Int_t irow=0; irow<160; irow++){
+    AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+    if (!cl) continue;
+    if (cl->GetX()<80) continue;
+    if (detector<0) detector=cl->GetDetector()%36;
+    Int_t sector = cl->GetDetector();
+    Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();    
+    if (TMath::Abs(dalpha)>0.01){
+      if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
+    }
+    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};    
+    Double_t cov[3]={0.01,0.,0.01}; 
+    AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
+    cov[0]+=1./(irow+1.); // bigger error at boundary
+    cov[0]+=1./(160.-irow); // bigger error at boundary
+    cov[2]+=1./(irow+1.); // bigger error at boundary
+    cov[2]+=1./(160.-irow); // bigger error at boundary
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
+    if (TMath::Abs(cl->GetY()-trackOut.GetY())<kMaxDist){
+      nclOut++;
+      trackOut.Update(&r[1],cov);
+    }
+    if (nclOut<kNclCut/2) continue;
+    if (cl->GetDetector()%36!=detector) continue;
+    //
+    // fill residual histogram
+    //
+    Double_t resVector[5]; 
+    trackOut.GetXYZ(xyz);
+    resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
+    if (resVector[1]<0) resVector[1]+=18;
+    resVector[2]= cl->GetX();
+    resVector[3]= cl->GetZ()/cl->GetX();
+    //
+    resVector[0]= cl->GetY()-trackOut.GetY();
+    fClusterDelta[0]->Fill(resVector);
+    resVector[0]= cl->GetZ()-trackOut.GetZ();
+    fClusterDelta[1]->Fill(resVector);
+  }
+  //
+  // Refit in - store residual maps
+  //
+  for (Int_t irow=159; irow>=0; irow--){
+    AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
+    if (!cl) continue;
+    if (cl->GetX()<80) continue;
+    if (detector<0) detector=cl->GetDetector()%36;
+    Int_t sector = cl->GetDetector();
+    Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();    
+    if (TMath::Abs(dalpha)>0.01){
+      if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
+    }
+    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};    
+    Double_t cov[3]={0.01,0.,0.01}; 
+    AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
+    cov[0]+=1./(irow+1.); // bigger error at boundary
+    cov[0]+=1./(160.-irow); // bigger error at boundary
+    cov[2]+=1./(irow+1.); // bigger error at boundary
+    cov[2]+=1./(160.-irow); // bigger error at boundary
+    cov[0]*=cov[0];
+    cov[2]*=cov[2];
+    if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
+    if (TMath::Abs(cl->GetY()-trackIn.GetY())<kMaxDist){
+      nclIn++;
+      trackIn.Update(&r[1],cov);
+    }
+    if (nclIn<kNclCut/2) continue;
+    if (cl->GetDetector()%36!=detector) continue;
+    //
+    // fill residual histogram
+    //
+    Double_t resVector[5]; 
+    trackIn.GetXYZ(xyz);
+    resVector[1]= 9.*TMath::ATan2(xyz[1],xyz[0])/TMath::Pi();
+    if (resVector[1]<0) resVector[1]+=18;
+    resVector[2]= cl->GetX();
+    resVector[3]= cl->GetZ()/cl->GetX();
+    //
+    resVector[0]= cl->GetY()-trackIn.GetY();
+    fClusterDelta[0]->Fill(resVector);
+    resVector[0]= cl->GetZ()-trackIn.GetZ();
+    fClusterDelta[1]->Fill(resVector);
+  }
+
+
+}
+
+
 void  AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
   //
-  // Update Kalman filter of Alignment 
+  // Update Kalman filter of Alignment - only setup without filed
   //       IROC - OROC quadrants
   //
+  if (TMath::Abs(AliTracker::GetBz())>0.5) return;
   if (!fClusterDelta[0])  MakeResidualHistos();
   const Int_t kMinClusterF=40;
+  const Int_t kMinClusterFit=10;
   const Int_t kMinClusterQ=10;
   //
   const  Int_t     kdrow1Fit =5;         // rows to skip from fit at the end      
@@ -2693,11 +2814,14 @@ void  AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
                                  c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
        if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
       }
-      fyf.AddPoint(x,c->GetY(),0.1);
-      fzf.AddPoint(x,c->GetZ(),0.1);
+      if (TMath::Abs(x[0])<10){
+       fyf.AddPoint(x,c->GetY(),0.1); //use only middle rows+-10cm
+      }
+      fzf.AddPoint(x,c->GetZ(),0.1);      
     }
     nf = fyf.GetNpoints();
-    if (nf<kMinClusterF) return;   // not enough points - skip 
+    if (fyf.GetNpoints()<kMinClusterFit) return;   // not enough points - skip 
+    if (fzf.GetNpoints()<kMinClusterF) return;   // not enough points - skip 
     fyf.Eval(); 
     fyf.GetParameters(pyf); 
     fyf.GetErrors(peyf);
@@ -2708,9 +2832,9 @@ void  AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
   //
   //
   //
-  TVectorD vecX(2*nf+kdrow0Fit+kdrow1Fit+5);          // x         vector
-  TVectorD vecY(2*nf+kdrow0Fit+kdrow1Fit+5);          // residuals vector
-  TVectorD vecZ(2*nf+kdrow0Fit+kdrow1Fit+5);                              // residuals vector
+  TVectorD vecX(160);          // x         vector
+  TVectorD vecY(160);          // residuals vector
+  TVectorD vecZ(160);                              // residuals vector
   TVectorD vPosG(3);                  //vertex position
   TVectorD vPosL(3);                 // vertex position in the TPC local system
   TVectorF vImpact(2);               //track impact parameter
@@ -2739,8 +2863,8 @@ void  AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
   // get constrained parameters
   //
   Double_t xvertex=vPosL[0]-fXmiddle;
-  fyf.AddPoint(&xvertex,vPosL[1], 0.1+TMath::Abs(vImpact[0]));
-  fzf.AddPoint(&xvertex,vPosL[2], 0.1+TMath::Abs(vImpact[1]));
+  fyf.AddPoint(&xvertex,vPosL[1], 0.00001);
+  fzf.AddPoint(&xvertex,vPosL[2], 2.);
   fyf.Eval(); 
   fyf.GetParameters(pyfc); 
   fzf.Eval(); 
@@ -2810,23 +2934,17 @@ void  AliTPCcalibAlign::UpdateAlignSector(const AliTPCseed * track,Int_t isec){
     //
     // Fill THnSparse cluster residuals
     // use only primary candidates with ITS signal
-    if (nf>100&&fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){    
+    if (fCurrentTrack->IsOn(0x4)&&TMath::Abs(vImpact[0])<1&&TMath::Abs(vImpact[1])<1){    
       Double_t resVector[5];
       resVector[1]= 9.*gphi/TMath::Pi();
       resVector[2]= c->GetX();
-      resVector[3]= c->GetY()/c->GetX();
-      resVector[4]= c->GetZ()/c->GetX();
+      resVector[3]= c->GetZ()/c->GetX();
       //
-      resVector[0]= c->GetY()-yfit;
-      //fClusterDelta[0]->Fill(resVector);
-      resVector[0]= c->GetZ()-zfit;
-      fClusterDelta[1]->Fill(resVector);
       //
       resVector[0]= c->GetY()-yfitC;
-      fClusterDelta[2]->Fill(resVector);
+      fClusterDelta[0]->Fill(resVector);
       resVector[0]= c->GetZ()-zfitC;
-      //fClusterDelta[3]->Fill(resVector);
-
+      fClusterDelta[1]->Fill(resVector);
     }
     if (c->GetRow()<kdrow0Fit) continue;      
     if (c->GetRow()>159-kdrow1Fit) continue;      
index 56545de..372ee15 100644 (file)
@@ -60,6 +60,7 @@ public:
                        const AliTPCseed * seed,
                        Int_t s1,Int_t s2);
   
+  void UpdateClusterDeltaField(const AliTPCseed * seed);
   void UpdateAlignSector(const AliTPCseed * seed,Int_t isec); 
   Int_t GetIndex(Int_t s1,Int_t s2) const {return 72*s1+s2;}
   //
@@ -140,12 +141,12 @@ public:
   void FillHisto(const Double_t *t1,
                 const Double_t *t2,
                 Int_t s1,Int_t s2);
-  void FillHisto(const AliExternalTrackParam *tp1,
-                const AliExternalTrackParam *tp2,
+  void FillHisto(AliExternalTrackParam *tp1,
+                AliExternalTrackParam *tp2,
                 Int_t s1,Int_t s2);
 
 protected:
-  THnSparse *fClusterDelta[6];  //clusters residuals
+  THnSparse *fClusterDelta[2];  //clusters residuals
   THnSparse *fTrackletDelta[4]; //track residuals
 
   TObjArray fDphiHistArray;    // array of residual histograms  phi      -kPhi
@@ -203,7 +204,7 @@ protected:
 private:
   AliTPCcalibAlign&  operator=(const AliTPCcalibAlign&);// not implemented
 
-  ClassDef(AliTPCcalibAlign,4)
+  ClassDef(AliTPCcalibAlign,6)
 };
 
 
diff --git a/TPC/CalibMacros/MakeAlignCalPad.C b/TPC/CalibMacros/MakeAlignCalPad.C
new file mode 100644 (file)
index 0000000..1ad219b
--- /dev/null
@@ -0,0 +1,625 @@
+/*
+  marian.ivanov@cern.ch
+  Macro to create  alignment/distortion maps
+  As a input output of AliTPCcalibAlign is used.
+
+  Algorithm:
+
+  In the setup without the magnetic field the tracks are fitted using the linear track model.
+  ( in y-rphi coordinate the primary vertex is also used as a constrain)
+  Residuals (deltas0  between the track and clusters in Y and in z direction are filled in the 4 D histogram:
+  Delta: phi(180 bins): localX(53 bins): tan(phi): tan(theta)(10 bins)
+
+  Distortion map are extracted form the residual histograms as a mean value at each bin.
+  Linear fits are then performed for each pad - delta as function of theta
+  Delta Ymeas = offsetY+slopeY*tan(theta)  
+  Delta Zmeas = offsetZ+slopeZ*tan(theta)  
+  
+  Resulting residuals exported into the OCDB are:
+  DeltaY  = offsetY
+  DeltaZ  = offsetZ
+  DeltaR  = slopeZ;
+  
+  Example usage:
+
+  make calpad+ make report ps file:
+  aliroot -b -q ~/NimStyle.C ../ConfigCalibTrain.C\(119037\) $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C\(1\)
+  make only report ps file:
+  aliroot -b -q ~/NimStyle.C $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C\(3\)
+  Making fit - iterative procedure - see below:
+  
+  gROOT->Macro("~/rootlogon.C");
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libSTAT");
+  gSystem->Load("libTPCcalib");
+  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros -I$ALICE_ROOT/TPC/TPC -I$ALICE_ROOT/STAT");
+  .L $ALICE_ROOT/TPC/CalibMacros/MakeAlignCalPad.C+
+  //load OCDB 
+  gROOT->Macro("../ConfigCalibTrain.C(119037)");
+  //2
+  InitTPCalign(); 
+  MakeFits();      // this is logn proceure 30 minutes
+  UpdateOCDB(0,AliCDBRunRange::Infinity());
+  //
+   gROOT->Macro("~/NimStyle.C")
+
+*/
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include "TH1D.h"
+#include "TH2F.h"
+#include "THnSparse.h"
+#include "TVectorD.h"
+#include "TTreeStream.h"
+#include "TFile.h"
+#include "AliTPCcalibAlign.h"
+#include "AliTPCROC.h"
+#include "AliTPCCalROC.h"
+#include "AliTPCCalPad.h"
+#include "TF1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TROOT.h"
+#include "TProfile.h"
+#include "AliTPCPreprocessorOnline.h"
+#include "AliTPCcalibDB.h"
+#include "AliTPCkalmanAlign.h"
+#include "TPostScript.h"
+#include "TLegend.h"
+#include "TCut.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "AliTPCExBEffectiveSector.h"
+#include "AliTPCComposedCorrection.h"
+//
+#include "AliCDBMetaData.h"
+#include "AliCDBId.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliCDBEntry.h"
+#include "TStatToolkit.h"
+#endif
+
+
+TObject  *palign=0;
+TTree * treeP=0;
+TTree * treeM=0;
+TTree * tree0=0;
+TTree * treePZ=0;
+TTree * treeMZ=0;
+TTree * tree0Z=0;
+AliTPCROC *proc = AliTPCROC::Instance();
+AliTPCExBEffectiveSector * geffCorr=0;
+
+Double_t GetCorr(Double_t sector, Double_t localX, Double_t kZ,Int_t type);
+
+void MakeFits();
+void InitTPCalign();
+
+
+void MakeAlignCalPad(Int_t mode){
+  //
+  // Make AlignCalpad and make report
+  //
+  gSystem->Load("libANALYSIS");
+  gSystem->Load("libSTAT");
+  gSystem->Load("libTPCcalib");  
+  printf("Make report mode\t%d\n", mode);
+  if (mode<1) { 
+    InitTPCalign();
+    MakeFits();
+  }
+  printf("Make report\n");
+}
+
+
+void InitTPCalign(){
+  //
+  // read the TPC alignment
+  //
+  TFile fcalib("CalibObjects.root");
+  TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
+  if (array){
+    palign = ( AliTPCcalibAlign *)array->FindObject("alignTPC");
+  }else{
+    palign = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
+  }
+  if (!palign){
+     TFile fcalib2("TPCAlignObjects.root");
+     palign = ( AliTPCcalibAlign *)fcalib2.Get("alignTPC");
+  }
+}
+
+
+void MakeFits(){
+  AliTPCcalibAlign * align =  (AliTPCcalibAlign *)palign;
+  THnSparse * hdY = align->GetClusterDelta(0);
+  THnSparse * hdZ = align->GetClusterDelta(1);
+  AliTPCExBEffectiveSector::MakeResidualMap(hdY,"clusterDY.root");
+  AliTPCExBEffectiveSector::MakeResidualMap(hdZ,"clusterDZ.root");
+}
+
+void FitdY(TTree * tree){
+  //
+  //
+  //
+
+  AliTPCROC * roc = AliTPCROC::Instance();
+  Double_t xquadrant = roc->GetPadRowRadii(36,53); 
+  Double_t xIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
+  
+  //
+  tree->SetAlias("diffS","(sector-int(sector)-0.5)");
+  tree->SetAlias("diffQ",Form("(localX-%f)",xquadrant));
+  tree->SetAlias("diffIO",Form("(localX-%f)",xIO));
+  tree->SetAlias("iroc",Form("(localX<%f)",xIO));   // bool IROC
+  tree->SetAlias("dqLR0",Form("(localX<%f)*(-1+(diffS<0)*2)",xIO));  //sign LeftRight
+  tree->SetAlias("dqLR2",Form("(localX>%f)*(-1+(diffS<0)*2)",xquadrant));  //sign LeftRight up
+  //
+  tree->SetAlias("diffIFC","abs(localX-80)");
+  tree->SetAlias("diffOFC","abs(250-localX)");
+  tree->SetAlias("drift","(1-abs(kZ*localX)/250.)");
+  //
+  tree->SetAlias("dIFC2",Form("drift/(1+abs(diffIFC^2/(%f^2)))",10.));
+  tree->SetAlias("dIFC",Form("drift/(1+abs(diffIFC/(%f)))",5.));
+  tree->SetAlias("dOFC2",Form("drift/(1+abs(diffOFC^2/(%f^2)))",10.));
+  tree->SetAlias("dOFC",Form("drift/(1+abs(diffOFC/(%f)))",5.));
+
+
+  TStatToolkit toolkit;
+  Double_t chi2=0;
+  Int_t    npoints=0;
+  TVectorD param;
+  TMatrixD covar;
+  TString  fstringG="";                 
+  fstringG+="(iroc)++";           //1 iroc shift
+  fstringG+="(dqLR0)++";
+  fstringG+="(dqLR2)++";
+  fstringG+="(diffIO*iroc)++";     // iroc rotation
+  fstringG+="(diffIO*dqLR0)++";
+  fstringG+="(diffQ*dqLR2)++";
+  //
+  fstringG+="(dIFC+dIFC2)++";            //9
+  fstringG+="(diffS*(dIFC+dIFC2))++";
+  fstringG+="(diffS^2*(dIFC+dIFC2))++";
+  fstringG+="(dIFC-dIFC2)++";            //9
+  fstringG+="(diffS*(dIFC-dIFC2))++";
+  fstringG+="(diffS^2*(dIFC-dIFC2))++";
+  //
+  //
+  fstringG+="(dOFC+dOFC2)++";            //9
+  fstringG+="(diffS*(dOFC+dOFC2))++";
+  fstringG+="(diffS^2*(dOFC+dOFC2))++";
+  fstringG+="(dOFC-dOFC2)++";            //9
+  fstringG+="(diffS*(dOFC-dOFC2))++";
+  fstringG+="(diffS^2*(dOFC-dOFC2))++";
+  //
+  //
+  //
+  TTreeSRedirector *pcstream=new TTreeSRedirector("dyFit.root");
+  {for (Int_t isec=0; isec<18; isec++){
+      {for (Int_t iside=0; iside<=1; iside++){
+         TCut cutS=Form("abs(sector-%d.5)<0.5&&kZ*%d>0&&entries>2000",isec,(iside>0)?-1:1);
+         TString *strFitG = TStatToolkit::FitPlane(tree,"mean", fstringG.Data(),cutS, chi2,npoints,param,covar,-1,0, 10000000, kTRUE);
+         tree->SetAlias("fitG",strFitG->Data());
+         tree->Draw("mean-fitG",cutS,"");
+         //
+         printf("%d\t%d\t%f\n",iside,isec,TMath::Sqrt(chi2/npoints));
+         (*pcstream)<<"fitDy"<<
+           "iside="<<iside<<
+           "sec="<<isec<<
+           "chi2="<<chi2<<
+           "npoints="<<npoints<<
+           "p.="<<&param<<
+           "\n";
+       }
+      }
+    }
+  }
+  delete pcstream;
+}
+
+
+
+void DumpDerivative(TH3 * his){
+  //
+  //
+  //
+  Int_t nbins0=his->GetXaxis()->GetNbins();
+  Int_t nbins1=his->GetYaxis()->GetNbins();
+  Int_t nbins2=his->GetZaxis()->GetNbins();
+  TTreeSRedirector *pcstream = new TTreeSRedirector("aaa.root");
+  for (Int_t ibin0=1; ibin0<nbins0; ibin0++)
+    for (Int_t ibin1=1; ibin1<nbins1; ibin1++)
+      for (Int_t ibin2=1; ibin2<nbins2; ibin2++){
+       Float_t x= his->GetXaxis()->GetBinCenter(ibin0);
+       Float_t y= his->GetYaxis()->GetBinCenter(ibin1);
+       Float_t z= his->GetZaxis()->GetBinCenter(ibin2);
+       Float_t c= his->GetBinContent(ibin0,ibin1,ibin2);
+       Float_t c0M=his->GetBinContent(ibin0-1,ibin1,ibin2);
+       Float_t c0P=his->GetBinContent(ibin0+1,ibin1,ibin2);
+       Float_t c1M=his->GetBinContent(ibin0,ibin1-1,ibin2);
+       Float_t c1P=his->GetBinContent(ibin0,ibin1+1,ibin2);
+       Float_t c2M=his->GetBinContent(ibin0,ibin1,ibin2-1);
+       Float_t c2P=his->GetBinContent(ibin0,ibin1,ibin2+1);
+       printf("%f\t%f\t%f\t%f\n",x,y,z,c);
+       (*pcstream)<<"aaa"<<
+         "x="<<x<<
+         "y="<<y<<
+         "z="<<z<<
+         "c="<<c<<
+         "c0M="<<c0M<<
+         "c0P="<<c0P<<
+         "c1M="<<c1M<<
+         "c1P="<<c1P<<
+         "c2M="<<c2M<<
+         "c2P="<<c2P<<
+         "\n";
+      }                                                
+  delete pcstream;
+}
+
+Double_t GetCorr(Double_t sector, Double_t localX, Double_t kZ, Int_t type){
+  //
+  // calculate the correction at given position - check the geffCorr
+  //
+  Double_t phi=sector*TMath::Pi()/9.;
+  Double_t gx = localX*TMath::Cos(phi);
+  Double_t gy = localX*TMath::Sin(phi);
+  Double_t gz = localX*kZ;
+  Int_t nsector=(gz>0) ? 0:18; 
+  //
+  //
+  //
+  Float_t distPoint[3]={gx,gy,gz};
+  geffCorr->DistortPoint(distPoint, nsector);
+  Double_t r0=TMath::Sqrt(gx*gx+gy*gy);
+  Double_t r1=TMath::Sqrt(distPoint[0]*distPoint[0]+distPoint[1]*distPoint[1]);
+  Double_t phi0=TMath::ATan2(gy,gx);
+  Double_t phi1=TMath::ATan2(distPoint[1],distPoint[0]);
+  if (type==0) return r1-r0;
+  if (type==1) return (phi1-phi0)*r0;
+  return phi1-phi0;
+}
+
+
+
+void LoadDistortionTrees(){
+  //
+  // Load distortion tree
+  //
+  TFile *fp = new TFile("clusterDYPlus.root");
+  TFile *fm = new TFile("clusterDYMinus.root");
+  TFile *f0 = new TFile("clusterDY0.root");
+  TFile *fpz= new TFile("clusterDZPlus.root");
+  TFile *fmz= new TFile("clusterDZMinus.root");
+  TFile *f0z=new TFile("clusterDZ0.root");
+  treeP=(TTree*)fp->Get("delta");
+  treeM=(TTree*)fm->Get("delta");
+  tree0=(TTree*)f0->Get("delta");
+  treePZ=(TTree*)fpz->Get("delta");
+  treeMZ=(TTree*)fmz->Get("delta");
+  tree0Z=(TTree*)f0z->Get("delta");
+  tree0->AddFriend(treeP,"P");
+  tree0->AddFriend(treeM,"M");
+  tree0->AddFriend(treePZ,"PZ");
+  tree0->AddFriend(treeMZ,"MZ");
+  tree0->AddFriend(tree0Z,"Z");
+  tree0->SetMarkerStyle(25);
+  tree0->SetMarkerSize(0.4);
+}
+
+void UpdateEffSectorOCDB(){
+  //
+  // Incremeantal update ot the correction maps
+  // corrections on top of previous corrections
+  //
+  TFile fp("clusterDYPlus.root");
+  TFile fm("clusterDYMinus.root");
+  TFile f0("clusterDY0.root");
+  TFile f0z("clusterDZ0.root");
+  
+  TH3F *his3D0=(TH3F*)f0.Get("his3D");
+  TH3F *his3DP=(TH3F*)fp.Get("his3D");
+  TH3F *his3DM=(TH3F*)fm.Get("his3D");
+  TH3F *his3DZ=(TH3F*)f0z.Get("his3D");
+  TH3F *hisR=0;   //half of difference (Plus-Minus) scaled by c1
+  hisR=(TH3F*)his3DP->Clone();
+  hisR->Add(his3DM,-1);
+  hisR->Scale(0.5);         
+  hisR->Scale(1/0.3); // c1 factor to be added there //is sign correct?        
+  //
+  //
+  AliTPCExBEffectiveSector *effSector = new  AliTPCExBEffectiveSector;
+  effSector->SetName("EffSector");
+  effSector->SetTitle("EffSector");
+  effSector->fCorrectionRPhi=(TH3F*)his3D0->Clone();
+  effSector->fCorrectionR=(TH3F*)hisR->Clone();
+  effSector->fCorrectionZ=(TH3F*)his3DZ->Clone();
+  geffCorr=effSector;
+  //
+  //
+  //
+  AliTPCcalibDB * calib = AliTPCcalibDB::Instance();
+  TObjArray *arr = calib->GetTPCComposedCorrectionArray();
+  {for (Int_t i=0; i<arr->GetEntries(); i++){
+    AliTPCComposedCorrection * corr=(AliTPCComposedCorrection*)arr->At(i);
+    if (!corr) continue;
+    AliTPCExBEffectiveSector * corrSec=( AliTPCExBEffectiveSector *)corr->GetCorrections()->FindObject("EffSector");
+    if (i==0) geffCorr=(AliTPCExBEffectiveSector *)corrSec;
+    if (!corrSec) corr->GetCorrections()->Add(effSector->Clone());
+    if (corrSec){
+      if (corrSec->fCorrectionRPhi) corrSec->fCorrectionRPhi->Add(effSector->fCorrectionRPhi);
+      else corrSec->fCorrectionRPhi=(TH3F*)effSector->fCorrectionRPhi->Clone();
+      if (corrSec->fCorrectionR) corrSec->fCorrectionR->Add(effSector->fCorrectionR);
+      else corrSec->fCorrectionR=(TH3F*)effSector->fCorrectionR->Clone();
+      if (corrSec->fCorrectionZ) corrSec->fCorrectionR->Add(effSector->fCorrectionZ);
+      else corrSec->fCorrectionZ=(TH3F*)effSector->fCorrectionZ->Clone();
+    }
+    }}
+// make OCDB entries
+  TString userName=gSystem->GetFromPipe("echo $USER");
+  TString date=gSystem->GetFromPipe("date");
+  TString ocdbStorage="local:////";
+  ocdbStorage+=gSystem->GetFromPipe("pwd");
+  ocdbStorage+="/OCDB";
+  //
+  AliCDBMetaData *metaData= new AliCDBMetaData();
+  metaData->SetObjectClassName("TObjArray");
+  metaData->SetResponsible("Marian Ivanov");
+  metaData->SetBeamPeriod(1);
+  metaData->SetAliRootVersion("05-27-04"); //root version
+  metaData->SetComment(Form("Correction calibration. User: %s\n Data%s",userName.Data(),date.Data()));
+  AliCDBId* id1=NULL;
+  id1=new AliCDBId("TPC/Calib/Correction", 0, AliCDBRunRange::Infinity());
+  AliCDBStorage* gStorage = 0;
+
+  gStorage=AliCDBManager::Instance()->GetStorage((ocdbStorage+"Update").Data());
+  gStorage->Put(arr, (*id1), metaData);
+
+
+}
+
+
+void DrawDiff(){
+  TFile f0("../mergeField0/mean.root");
+  TFile fP("../mergePlus/mean.root");
+  TFile fM("../mergeMinus/mean.root");
+  //
+  TTree *itsdy0=(TTree*)f0.Get("ITSdy");
+  TTree *itsdyP=(TTree*)fP.Get("ITSdy");
+  TTree *itsdyM=(TTree*)fM.Get("ITSdy");
+  TTree *vdy0=(TTree*)f0.Get("Vertexdy");
+  TTree *vdyP=(TTree*)fP.Get("Vertexdy");
+  TTree *vdyM=(TTree*)fM.Get("Vertexdy");
+  itsdy0->SetMarkerStyle(25);
+  itsdy0->SetMarkerSize(0.3);
+  itsdy0->AddFriend(itsdyP,"P");
+  itsdy0->AddFriend(itsdyM,"M");
+  itsdy0->AddFriend(vdy0,"V");
+  itsdy0->AddFriend(vdyP,"VP");
+  itsdy0->AddFriend(vdyM,"VM");
+  itsdy0->SetMarkerStyle(25);
+  itsdy0->SetMarkerSize(0.4);
+
+}
+
+
+
+
+
+void MakePlotDeltaZ(){
+  //
+  // 
+  //
+  TCut cut="entries>500&&PZ.entries>500&&MZ.entries>500";
+  TCanvas *cA = new TCanvas("deltaZA","deltaZA",900,700);
+  TCanvas *cC = new TCanvas("deltaZC","deltaZC",900,700);
+  TCanvas *cARef = new TCanvas("deltaZARef","deltaZARef",1000,800);
+  TCanvas *cCRef = new TCanvas("deltaZCRef","deltaZCRef",1000,800);
+  //
+  TH1::AddDirectory(0);
+  TH1 * his=0;
+  cA->Divide(4,5);
+  for (Int_t isec=0; isec<18; isec++){
+    cA->cd(isec+1);
+    TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
+    tree0->Draw("Z.mean*10.:localX:abs(kZ)",cutS+"entries>4000&&P.entries>200&&localX<250&&kZ>0&&abs(kZ)<1","colz");
+    his=(TH1*)tree0->GetHistogram()->Clone();
+    his->GetXaxis()->SetTitle("r (cm)");
+    his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
+    his->GetZaxis()->SetTitle("tan($theta)");
+    his->Draw("colz");
+  }
+  cA->cd(19);
+
+  cC->Divide(4,5);
+  for (Int_t isec=0; isec<18; isec++){
+    cC->cd(isec+1);
+    TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
+    tree0->Draw("Z.mean*10.:localX:abs(kZ)",cutS+"entries>4000&&P.entries>200&&localX<250&&kZ<0&&abs(kZ)<1","colz"); 
+    his=(TH1*)tree0->GetHistogram()->Clone();
+    his->GetXaxis()->SetTitle("r (cm)");
+    his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
+    his->GetZaxis()->SetTitle("tan($theta)");
+    his->Draw("colz");
+
+  }
+
+  cARef->Divide(4,5);
+  for (Int_t isec=0; isec<18; isec++){
+    cARef->cd(isec+1);
+    TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
+    tree0->Draw("(PZ.mean-MZ.mean)*10.:localX:abs(kZ)",cutS+cut+"kZ>0&&abs(kZ)<1","colz");
+    his=(TH1*)tree0->GetHistogram()->Clone();
+    his->GetXaxis()->SetTitle("r (cm)");
+    his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
+    his->GetZaxis()->SetTitle("tan($theta)");
+    his->Draw("colz");
+
+  }
+  
+  cCRef->Divide(4,5);
+  for (Int_t isec=0; isec<18; isec++){
+    cCRef->cd(isec+1);
+    TCut cutS=Form("abs(sector-0.5-%d)<0.1",isec);
+    tree0->Draw("(PZ.mean-MZ.mean)*10.:localX:abs(kZ)",cutS+cut+"kZ0<0&&abs(kZ)<1","colz");
+    his=(TH1*)tree0->GetHistogram()->Clone();
+    his->GetXaxis()->SetTitle("r (cm)");
+    his->GetYaxis()->SetTitle("#Delta_{z} (mm)");
+    his->GetZaxis()->SetTitle("tan($theta)");
+    his->Draw("colz");
+  }
+
+  TPostScript *ps = new TPostScript("distortionZ.ps", 112);
+  ps->NewPage();
+  cA->Draw();
+  ps->NewPage();
+  cA->Draw();
+  ps->NewPage();
+  cC->Draw();
+  ps->NewPage();
+  cARef->Draw();
+  ps->NewPage();
+  cCRef->Draw();
+  ps->Close();
+}
+
+
+
+
+void MakeAlign(){
+  TFile fcalib("../mergeField0/TPCAlignObjects.root");
+  AliTPCcalibAlign * align = ( AliTPCcalibAlign *)fcalib.Get("alignTPC");
+  TFile f0("../mergeField0/mean.root");
+  TFile fP("../mergePlus/mean.root");
+  TFile fM("../mergeMinus/mean.root");
+
+  //
+  TTree *itsdy0=(TTree*)f0.Get("ITSdy");
+  TTree *itsdyP=(TTree*)fP.Get("ITSdy");
+  TTree *itsdyM=(TTree*)fM.Get("ITSdy");
+  TTree *itsdp0=(TTree*)f0.Get("ITSdsnp");
+  TTree *itsdpP=(TTree*)fP.Get("ITSdsnp");
+  TTree *itsdpM=(TTree*)fM.Get("ITSdsnp");
+  TTree *itsdz0=(TTree*)f0.Get("ITSdz");
+  TTree *itsdzP=(TTree*)fP.Get("ITSdz");
+  TTree *itsdzM=(TTree*)fM.Get("ITSdz");
+  TTree *vdy0=(TTree*)f0.Get("Vertexdy");
+  TTree *vdyP=(TTree*)fP.Get("Vertexdy");
+  TTree *vdyM=(TTree*)fM.Get("Vertexdy");
+  TTree *vds0=(TTree*)f0.Get("Vertexdsnp");
+  TTree *vdsP=(TTree*)fP.Get("Vertexdsnp");
+  TTree *vdsM=(TTree*)fM.Get("Vertexdsnp");
+  itsdy0->SetMarkerStyle(25);
+  itsdy0->SetMarkerSize(0.3);
+  itsdy0->AddFriend(itsdyP,"P");
+  itsdy0->AddFriend(itsdyM,"M");
+  itsdy0->AddFriend(vdy0,"V");
+  itsdy0->AddFriend(vdyP,"VP");
+  itsdy0->AddFriend(vdyM,"VM");
+  itsdy0->AddFriend(itsdz0,"Z");
+  itsdy0->AddFriend(itsdzP,"ZM");
+  itsdy0->AddFriend(itsdzM,"ZP");
+  itsdy0->AddFriend(itsdp0,"Snp");
+  itsdy0->AddFriend(itsdpP,"SnpM");
+  itsdy0->AddFriend(itsdpM,"SnpP");
+  itsdy0->AddFriend(vds0,"VSnp");
+  itsdy0->AddFriend(vdsP,"VSnpM");
+  itsdy0->AddFriend(vdsM,"VSnpP");
+
+  itsdy0->SetMarkerStyle(25);
+  itsdy0->SetMarkerSize(0.4);
+
+  TCut cutQ="entries>500&&abs(theta)<0.8&&abs(snp)<0.2";
+  TH1F his1("hdeltaY1","hdeltaY1",100,-0.5,0.5);
+  TMatrixD vecAlign(72,1);
+  TMatrixD covAlign(72,72);
+  AliTPCkalmanAlign::BookAlign1D(vecAlign,covAlign,0,0.0005);
+  TVectorD vecITSY(72);
+  TVectorD vecITSS(72);
+  TVectorD vecVS(72);
+  {for (Int_t isec0=0; isec0<36; isec0++){
+      Double_t phi0=(isec0%18+0.5)*TMath::Pi()/9.;
+      if (phi0>TMath::Pi()) phi0-=TMath::TwoPi();
+      Int_t iside0=(isec0%36<18)? 0:1;
+      TCut cutSector=Form("abs(%f-phi)<0.14",phi0);
+      TCut cutSide = (iside0==0)? "theta>0":"theta<0";
+      itsdy0->Draw("mean",cutQ+cutSector+cutSide);
+      Double_t meanITSY=itsdy0->GetHistogram()->GetMean()/83.6;
+      vecITSY[isec0]=meanITSY;
+      vecITSY[isec0+36]=meanITSY;
+      itsdy0->Draw("Snp.mean",cutQ+cutSector+cutSide);
+      Double_t meanITSS=itsdy0->GetHistogram()->GetMean();
+      vecITSS[isec0]=meanITSS;
+      vecITSS[isec0+36]=meanITSS;
+      itsdy0->Draw("VSnp.mean",cutQ+cutSector+cutSide);
+      Double_t meanVS=itsdy0->GetHistogram()->GetMean();
+      vecVS[isec0]=meanVS;
+      vecVS[isec0+36]=meanVS;
+    }
+  }
+  AliTPCROC * roc = AliTPCROC::Instance();
+  Double_t fXIROC = (roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(0,roc->GetNRows(0)-1))*0.5;
+  Double_t fXOROC = (roc->GetPadRowRadii(36,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+  Double_t fXmiddle   = ( roc->GetPadRowRadii(0,0)+roc->GetPadRowRadii(36,roc->GetNRows(36)-1))*0.5;
+  Double_t fXIO       = ( roc->GetPadRowRadii(0,roc->GetNRows(0)-1)+roc->GetPadRowRadii(36,0))*0.5;
+
+  TTreeSRedirector *pcstream=new TTreeSRedirector("combAlign.root");
+  
+  {
+    for (Int_t iter=0; iter<2; iter++){
+    for (Int_t isec0=0; isec0<72; isec0++){
+    for (Int_t isec1=0; isec1<72; isec1++){
+      TH1 * his = align->GetHisto(AliTPCcalibAlign::kY,isec0,isec1);
+      TH1 * hisPhi = align->GetHisto(AliTPCcalibAlign::kPhi,isec0,isec1);
+      if (!his) continue;
+      if (his->GetEntries()<100) continue;
+      Double_t xref=fXIO;
+      if (isec0<36&&isec1<36) xref=fXIROC;
+      if (isec0>=36&&isec1>=36) xref=fXOROC;
+      Double_t meanTPC=his->GetMean()/xref;
+      Double_t meanTPCPhi=hisPhi->GetMean();
+      Double_t meanITS0=vecITSY[isec0];
+      Double_t meanITS1=vecITSY[isec1];
+      Double_t meanITSS0=vecITSS[isec0];
+      Double_t meanITSS1=vecITSS[isec1];
+      Double_t meanVS0=vecVS[isec0];
+      Double_t meanVS1=vecVS[isec1];
+      AliTPCkalmanAlign::UpdateAlign1D(meanTPC, 0.001, isec0,isec1,  vecAlign,covAlign);
+      AliTPCkalmanAlign::UpdateAlign1D(meanTPCPhi, 0.001, isec0,isec1,  vecAlign,covAlign);
+      AliTPCkalmanAlign::UpdateAlign1D(meanITS1-meanITS0, 0.001, isec0,isec1,  vecAlign,covAlign);
+      AliTPCkalmanAlign::UpdateAlign1D(meanITSS1-meanITSS0, 0.001, isec0,isec1,  vecAlign,covAlign);
+      AliTPCkalmanAlign::UpdateAlign1D(meanVS1-meanVS0, 0.001, isec0,isec1,  vecAlign,covAlign);
+      //printf("isec0\t%d\tisec1\t%d\t%f\t%f\t%f\n",isec0,isec1, meanTPC, meanITS0,meanITS1);
+      Double_t kalman0= vecAlign(isec0,0);
+      Double_t kalman1= vecAlign(isec1,0);
+      if (iter>0) (*pcstream)<<"align"<<
+       "iter="<<iter<<
+       "xref="<<xref<<
+       "isec0="<<isec0<<
+       "isec1="<<isec1<<
+       "mTPC="<<meanTPC<<
+       "mTPCPhi="<<meanTPCPhi<<
+       "mITS0="<<meanITS0<<
+       "mITS1="<<meanITS1<<
+       "mITSS0="<<meanITSS0<<
+       "mITSS1="<<meanITSS1<<
+       "mVS0="<<meanVS0<<
+       "mVS1="<<meanVS1<<
+       "k0="<<kalman0<<
+       "k1="<<kalman1<<
+       "\n";
+    }          
+    }
+    }
+  }
+  delete pcstream;
+  TFile f("combAlign.root");
+  TTree * treeA = (TTree*)f.Get("align"); 
+  treeA->SetMarkerStyle(25);
+  treeA->SetMarkerSize(0.5);
+}
+
+
+