TPC/AliTPCCalibGlobalMisalignment.cxx

   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 ////////////////////////////////////////////////////////////////////////////
  18 //                                                                        //
  19 // AliTPCCalibGlobalMisalignment class                                        //
  20 // The class calculates the space point distortions due to simple         //
  21 // misalignments like shifts in caresian coordinates or a rotation        //
  22 // of the TPC read out planes (A and C side)                              //
  23 // Optionaly possible to use it for visualization of the alignemnt form the Alignment OCDB //
  24 // fUseGeoManager has to be set to kTRUE to enable this option
  25 //                                                                        //
  26 // date: 06/05/2010                                                       //
  27 // Authors: Stefan Rossegger, Jim Thomas, Magnus Mager                    //
  28 ////////////////////////////////////////////////////////////////////////////
  29
  30 #include "AliTPCCalibGlobalMisalignment.h"
  31 #include "TMath.h"
  32 #include "TGeoMatrix.h"
  33 #include "AliTPCROC.h"
  34 #include "AliTPCcalibDB.h"
  35 #include "AliTPCParam.h"
  36 #include <TGeoPhysicalNode.h>
  37 //
  38 #include "AliAlignObjParams.h"
  39 #include "AliGeomManager.h"
  40 #include "AliCDBManager.h"
  41 #include "AliCDBEntry.h"
  42
  43 AliTPCCalibGlobalMisalignment::AliTPCCalibGlobalMisalignment()
  44   : AliTPCCorrection("mialign","Misalignment"),
  45     fXShift(0.),fYShift(0.),fZShift(0.),
  46     fRotPhiA(0.),fRotPhiC(0.),
  47     fdRPhiOffsetA(0.),
  48     fdRPhiOffsetC(0.),
  49     fQuadrantQ0(0),   //OROC medium pads -delta ly+ - ly - shift
  50     fQuadrantRQ0(0),  //OROC medium pads -delta ly+ - ly - rotation
  51     fQuadrantQ1(0),   //OROC long   pads -delta ly+ - ly - shift
  52     fQuadrantQ2(0),   //OROC long   pads -shift
  53     fQuadrantRQ1(0),  //OROC long   pads -delta ly+ - ly - rotation
  54     fQuadrantRQ2(0),  //OROC long   pads -rotation
  55     fMatrixGlobal(0), // global Alignment common
  56     fArraySector(0)   // fArraySector
  57 {
  58   //
  59   // default constructor
  60   //
  61 }
  62
  63 AliTPCCalibGlobalMisalignment::~AliTPCCalibGlobalMisalignment() {
  64   //
  65   // destructor
  66   //
  67   delete fQuadrantQ0;   //OROC medium pads -delta ly+ - ly - shift
  68   delete fQuadrantRQ0;  //OROC medium pads -delta ly+ - ly - rotation
  69   delete fQuadrantQ1;   //OROC long   pads -delta ly+ - ly - shift
  70   delete fQuadrantQ2;   //OROC long   pads -shift
  71   delete fQuadrantRQ1;  //OROC long   pads -delta ly+ - ly - rotation
  72   delete fQuadrantRQ2;  //OROC long   pads -rotation
  73   delete fMatrixGlobal; // global matrix
  74   delete fArraySector;  // sector matrices
  75 }
  76
  77 void AliTPCCalibGlobalMisalignment::SetQuadranAlign(const TVectorD *quadrantQ0, const TVectorD *quadrantRQ0, const TVectorD *quadrantQ1,const TVectorD *quadrantRQ1,  const TVectorD *quadrantQ2,  const TVectorD *quadrantRQ2){
  78   //
  79   // Set quadrant alignment
  80   // 6 vectors for 36 (super) sectors
  81   //
  82   if (quadrantQ0) fQuadrantQ0   = new TVectorD(*quadrantQ0);
  83   if (quadrantRQ0) fQuadrantRQ0 = new TVectorD(*quadrantRQ0);
  84   //
  85   if (quadrantQ1) fQuadrantQ1   = new TVectorD(*quadrantQ1);
  86   if (quadrantQ1) fQuadrantRQ1  = new TVectorD(*quadrantRQ1);
  87   if (quadrantQ2) fQuadrantQ2   = new TVectorD(*quadrantQ2);
  88   if (quadrantQ2) fQuadrantRQ2  = new TVectorD(*quadrantRQ2);
  89 }
  90
  91
  92
  93 void AliTPCCalibGlobalMisalignment::SetAlignGlobal(const TGeoMatrix * matrixGlobal){
  94   //
  95   // Set global misalignment
  96   // Object is OWNER
  97   //
  98   if (fMatrixGlobal) delete fMatrixGlobal;
  99   fMatrixGlobal=0;
 100   if (matrixGlobal) fMatrixGlobal = new TGeoHMatrix(*matrixGlobal);
 101 }
 102
 103 void AliTPCCalibGlobalMisalignment::SetAlignSectors(const TObjArray *arraySector){
 104   //
 105   // Set misalignment TObjArray of TGeoMatrices  - for each sector
 106   // Object is OWNER
 107   //
 108   if (fArraySector) delete fArraySector;
 109   fArraySector=0;
 110   if (arraySector) fArraySector = (TObjArray*)arraySector->Clone();
 111 }
 112
 113
 114 //void AliTPCCalibGlobalMisalignment::Init() {
 115 //  //
 116 // // Initialization funtion
 117 //  //
 118
 119 //  // nothing to be initialized, results of this calibration class will go to the global aligment structure
 120
 121 //}
 122
 123 //void AliTPCCalibGlobalMisalignment::Update(const TTimeStamp &/*timeStamp*/) {
 124 //  //
 125 //  // Update function
 126 //  //
 127 //
 128 //  // nothing to be updated, results of this calibration class will go to the global aligment structure
 129 //
 130 //}
 131
 132
 133
 134 void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 135   //
 136   // Calculates the simple correction due to a shift (in x,y,z) or an rotation of the TPC (around z)
 137   //
 138   static AliTPCROC *tpcRoc =AliTPCROC::Instance();
 139   Double_t xref  = ( tpcRoc->GetPadRowRadii(0,0)+tpcRoc->GetPadRowRadii(36,tpcRoc->GetNRows(36)-1))*0.5;
 140
 141   Double_t r=0, phi=0;
 142   r   = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] );
 143   phi = TMath::ATan2(x[1],x[0]);
 144   // Getsector number
 145   Double_t sec=TMath::Nint(-0.5+(phi*9./TMath::Pi()));
 146   if (sec<0) sec+=18;
 147   Int_t isec = TMath::Nint(sec);
 148   if  (roc%36>=18) isec+=18;
 149   //
 150   // Get the point on the local coordiante frame
 151   //
 152   Double_t alpha=(sec+0.5)*TMath::Pi()/9;
 153   Double_t pos[3]={0,0,x[2]};
 154   pos[0]=  TMath::Cos(alpha)*x[0]+TMath::Sin(alpha)*x[1];
 155   pos[1]= -TMath::Sin(alpha)*x[0]+TMath::Cos(alpha)*x[1];
 156   if (pos[0]>tpcRoc->GetPadRowRadiiUp(0))  isec+=36;
 157
 158   //
 159   // apply quadrant alignment if available - in local coordinate frame
 160   //
 161   //
 162   Double_t posQG[3]={x[0],x[1],x[2]};
 163   {
 164     Double_t dly=0;
 165     Bool_t isQ0 = TMath::Abs(pos[0]-161)<28;
 166     Bool_t isQ1 = (pos[0]>189);
 167     Double_t  sign = (pos[1]>0)? 1.: -1.;
 168     if (isQ0){
 169       if (fQuadrantQ0)  dly+=sign*(*fQuadrantQ0)[isec%36];  // shift in cm
 170       if (fQuadrantRQ0) dly+=sign*(*fQuadrantRQ0)[isec%36]*(pos[0]-xref);
 171     }
 172     if (isQ1){
 173       if (fQuadrantQ1)  dly+=sign*(*fQuadrantQ1)[isec%36];  // shift in cm
 174       if (fQuadrantRQ1) dly+=sign*(*fQuadrantRQ1)[isec%36]*(pos[0]-xref);
 175       if (fQuadrantQ2)  dly+=(*fQuadrantQ2)[isec%36];  // shift in cm
 176       if (fQuadrantRQ2) dly+=(*fQuadrantRQ2)[isec%36]*(pos[0]-xref);
 177     }
 178     // Tranform the corrected point to the global frame
 179     posQG[0]=  TMath::Cos(alpha)*pos[0]-TMath::Sin(alpha)*(pos[1]+dly);
 180     posQG[1]=  TMath::Sin(alpha)*pos[0]+TMath::Cos(alpha)*(pos[1]+dly);
 181   }
 182   //
 183   // rotation of the read-out planes
 184   if  (roc%36<18) // A side
 185     phi += fRotPhiA;
 186   else         // C side
 187     phi += fRotPhiC;
 188
 189   // Simply adding a constant dRPHi residual. PURELY FOR CALIBRATION PURPOSES
 190   if  (roc%36<18) // A side
 191     phi += fdRPhiOffsetA/r;
 192   else         // C side
 193     phi += fdRPhiOffsetC/r;
 194
 195   dx[0] = r * TMath::Cos(phi) - x[0];
 196   dx[1] = r * TMath::Sin(phi) - x[1];
 197   dx[2] = 0.;
 198
 199   // Simple shifts
 200   dx[0] -= fXShift;
 201   dx[1] -= fYShift;
 202   dx[2] -= fZShift;
 203   // quadrant shifts
 204   dx[0] += (posQG[0]-x[0]);
 205   dx[1] += (posQG[1]-x[1]);
 206   //
 207   //
 208   if (fMatrixGlobal){
 209     // apply global alignment matrix
 210     Double_t ppos[3]={x[0],x[1],x[2]};
 211     Double_t pposC[3]={x[0],x[1],x[2]};
 212     fMatrixGlobal->LocalToMaster(ppos,pposC);
 213     dx[0]+=pposC[0]-ppos[0];
 214     dx[1]+=pposC[1]-ppos[1];
 215     dx[2]+=pposC[2]-ppos[2];
 216   }
 217
 218   if (fArraySector){
 219     // apply global alignment matrix
 220     TGeoMatrix  *mat = (TGeoMatrix*)fArraySector->At(isec);
 221     if (mat){
 222       Double_t ppos[3]={x[0],x[1],x[2]};
 223       Double_t pposC[3]={x[0],x[1],x[2]};
 224       mat->LocalToMaster(ppos,pposC);
 225       dx[0]+=pposC[0]-ppos[0];
 226       dx[1]+=pposC[1]-ppos[1];
 227       dx[2]+=pposC[2]-ppos[2];
 228     }
 229   }
 230 }
 231
 232 void AliTPCCalibGlobalMisalignment::Print(Option_t* /*option*/ ) const {
 233   //
 234   // Print function to check the settings
 235   //
 236   printf("%s",GetTitle());
 237   printf(" - Trivial Misalignments for calibration purposes: \n");
 238   printf(" - X-Shift: %1.3f cm, Y-Shift: %1.3f cm, Z-Shift: %1.3f cm \n",fXShift,fYShift,fZShift);
 239   printf(" - Phi-Rotations: A side: %1.5f rad, C side: %1.5f rad\n",fRotPhiA,fRotPhiC);
 240   printf(" - dRPhi offsets: A side: %1.5f cm, C side: %1.5f cm\n",fdRPhiOffsetA,fdRPhiOffsetC);
 241
 242
 243 }
 244
 245 void AliTPCCalibGlobalMisalignment::AddAlign(const  AliTPCCalibGlobalMisalignment & add){
 246   //
 247   // Add the alignmnet to current object
 248   //
 249   fXShift+=add.fXShift;               // Shift in global X [cm]
 250   fYShift+=add.fYShift;               // Shift in global Y [cm]
 251   fZShift+=add.fZShift;               // Shift in global Z [cm]
 252
 253   fRotPhiA+=add.fRotPhiA;      // simple rotation of A side read-out plane around the Z axis [rad]
 254   fRotPhiC+=add.fRotPhiC;      // simple rotation of C side read-out plane around the Z axis [rad]
 255   fdRPhiOffsetA+=add.fdRPhiOffsetA;  // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes!
 256   fdRPhiOffsetC+=add.fdRPhiOffsetC;  // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes!
 257   //
 258   // Quadrant alignment
 259   //
 260   if (add.fQuadrantQ0) {
 261     if (fQuadrantQ0) fQuadrantQ0->Add(*(add.fQuadrantQ0));
 262     if (!fQuadrantQ0) fQuadrantQ0 = (TVectorD*)(add.fQuadrantQ0->Clone());
 263   }
 264   if (add.fQuadrantRQ0) {
 265     if (fQuadrantRQ0) fQuadrantRQ0->Add(*(add.fQuadrantRQ0));
 266     if (!fQuadrantRQ0) fQuadrantRQ0 = (TVectorD*)(add.fQuadrantRQ0->Clone());
 267   }
 268   //
 269   if (add.fQuadrantQ1) {
 270     if (fQuadrantQ1) fQuadrantQ1->Add(*(add.fQuadrantQ1));
 271     if (!fQuadrantQ1) fQuadrantQ1 = (TVectorD*)(add.fQuadrantQ1->Clone());
 272   }
 273   if (add.fQuadrantRQ1) {
 274     if (fQuadrantRQ1) fQuadrantRQ1->Add(*(add.fQuadrantRQ1));
 275     if (!fQuadrantRQ1) fQuadrantRQ1 = (TVectorD*)(add.fQuadrantRQ1->Clone());
 276   }
 277   //
 278   if (add.fQuadrantQ2) {
 279     if (fQuadrantQ2) fQuadrantQ2->Add(*(add.fQuadrantQ2));
 280     if (!fQuadrantQ2) fQuadrantQ2 = (TVectorD*)(add.fQuadrantQ2->Clone());
 281   }
 282   if (add.fQuadrantRQ2) {
 283     if (fQuadrantRQ2) fQuadrantRQ2->Add(*(add.fQuadrantRQ2));
 284     if (!fQuadrantRQ2) fQuadrantRQ2 = (TVectorD*)(add.fQuadrantRQ2->Clone());
 285   }
 286   //
 287   // Global alignment - use native ROOT representation
 288   //
 289   if (add.fMatrixGlobal){
 290     if (!fMatrixGlobal)  fMatrixGlobal = new TGeoHMatrix(*(add.fMatrixGlobal));
 291     if (fMatrixGlobal)   ((TGeoHMatrix*)fMatrixGlobal)->Multiply(add.fMatrixGlobal);
 292   }
 293   if (add.fArraySector){
 294     if (!fArraySector) {SetAlignSectors(add.fArraySector);
 295     }else{
 296       for (Int_t isec=0; isec<72; isec++){
 297         TGeoHMatrix *mat0= (TGeoHMatrix*)fArraySector->At(isec);
 298         TGeoHMatrix *mat1= (TGeoHMatrix*)add.fArraySector->At(isec);
 299         if (mat0&&mat1) mat0->Multiply(mat1);
 300       }
 301     }
 302   }
 303 }
 304
 305
 306 AliTPCCalibGlobalMisalignment *  AliTPCCalibGlobalMisalignment::CreateOCDBAlign(){
 307   //
 308   // Create  AliTPCCalibGlobalMisalignment from OCDB Alignment entry
 309   // OCDB has to be initialized before in user code
 310   // All storages (defualt and specific)  and run number
 311   //
 312   AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
 313   if (!entry){
 314     printf("Missing alignmnet entry. OCDB not initialized?\n");
 315     return 0;
 316   }
 317   TClonesArray * array = (TClonesArray*)entry->GetObject();
 318   Int_t entries = array->GetEntries();
 319   TGeoHMatrix matrixGlobal;
 320   TObjArray *alignArrayOCDB= new TObjArray(73);  // sector misalignment + global misalignment
 321   //                                            // global is number 72
 322   //
 323   { for (Int_t i=0;i<entries; i++){
 324       //
 325       //
 326       TGeoHMatrix matrix;
 327       AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
 328       alignP->GetMatrix(matrix);
 329       Int_t imod;
 330       AliGeomManager::ELayerID ilayer;
 331       alignP->GetVolUID(ilayer, imod);
 332       if (ilayer==AliGeomManager::kInvalidLayer) {
 333         alignArrayOCDB->AddAt(matrix.Clone(),72);
 334         alignP->GetMatrix(matrixGlobal);
 335       }else{
 336         Int_t sector=imod;
 337         if (ilayer==AliGeomManager::kTPC2) sector+=36;
 338         alignArrayOCDB->AddAt(matrix.Clone(),sector);
 339       }
 340     }
 341   }
 342   AliTPCCalibGlobalMisalignment *align = new  AliTPCCalibGlobalMisalignment;
 343   align->SetAlignGlobal(&matrixGlobal);
 344   align->SetAlignSectors(alignArrayOCDB);
 345   return align;
 346 }
 347
 348
 349 AliTPCCalibGlobalMisalignment *  AliTPCCalibGlobalMisalignment::CreateMeanAlign(const AliTPCCalibGlobalMisalignment *alignIn){
 350   //
 351   // Create new object, disantangle common mean alignmnet and sector alignment
 352   //
 353   // 1. Try to get mean alignment
 354   // 2. Remove mean alignment from sector alignment
 355   // 3. Create new object
 356   //
 357   TObjArray * array = alignIn->GetAlignSectors();
 358   TObjArray * arrayNew = new TObjArray(72);
 359   //
 360   //Get mean transformation
 361   TGeoHMatrix matrix;
 362   {for (Int_t isec=0; isec<72; isec++){
 363       const TGeoMatrix* cmatrix=(TGeoMatrix*)array->At(isec);
 364       if (!cmatrix) continue;
 365       matrix.Multiply(cmatrix);
 366     }}
 367   TGeoHMatrix matrixMean(matrix);
 368   matrixMean.SetDx(matrix.GetTranslation()[0]/72.);
 369   matrixMean.SetDy(matrix.GetTranslation()[1]/72.);
 370   matrixMean.SetDz(matrix.GetTranslation()[2]/72.);
 371   Double_t rotation[12];
 372   {for (Int_t i=0; i<12; i++) {
 373       rotation[i]=1.0;
 374       if (TMath::Abs(matrix.GetRotationMatrix()[i]-1.)>0.1){
 375       rotation[i]=matrix.GetRotationMatrix()[i]/72.;
 376       }
 377     }}
 378   matrixMean.SetRotation(rotation);
 379   TGeoHMatrix matrixInv = matrixMean.Inverse();
 380   //
 381   {for (Int_t isec=0; isec<72; isec++){
 382       TGeoHMatrix* amatrix=(TGeoHMatrix*)(array->At(isec)->Clone());
 383       if (!amatrix) continue;
 384       amatrix->Multiply(&matrixInv);
 385       arrayNew->AddAt(amatrix,isec);
 386     }}
 387   if (alignIn->GetAlignGlobal()) matrixMean.Multiply((alignIn->GetAlignGlobal()));
 388   AliTPCCalibGlobalMisalignment *alignOut = new  AliTPCCalibGlobalMisalignment;
 389   alignOut->SetAlignGlobal(&matrixMean);
 390   alignOut->SetAlignSectors(arrayNew);
 391   /*
 392     Checks transformation:
 393     AliTPCCalibGlobalMisalignment *  alignIn =  AliTPCCalibGlobalMisalignment::CreateOCDBAlign()
 394     AliTPCCalibGlobalMisalignment * alignOut =  AliTPCCalibGlobalMisalignment::CreateMeanAlign(alignIn)
 395     alignOutM= (AliTPCCalibGlobalMisalignment*)alignOut->Clone();
 396     alignOutS= (AliTPCCalibGlobalMisalignment*)alignOut->Clone();
 397     alignOutS->SetAlignGlobal(0);
 398     alignOutM->SetAlignSectors(0);
 399     //
 400     AliTPCCorrection::AddVisualCorrection(alignOut,0);
 401     AliTPCCorrection::AddVisualCorrection(alignOutM,1);
 402     AliTPCCorrection::AddVisualCorrection(alignOutS,2);
 403     AliTPCCorrection::AddVisualCorrection(alignIn,3);
 404
 405     TF1 f0("f0","AliTPCCorrection::GetCorrSector(x,85,0.9,1,0)",0,18);
 406     TF1 f1("f1","AliTPCCorrection::GetCorrSector(x,85,0.9,1,1)",0,18);
 407     TF1 f2("f2","AliTPCCorrection::GetCorrSector(x,85,0.9,1,2)",0,18);
 408     TF1 f3("f3","AliTPCCorrection::GetCorrSector(x,85,0.9,1,3)",0,18);
 409     f0->SetLineColor(1);
 410     f1->SetLineColor(2);
 411     f2->SetLineColor(3);
 412     f3->SetLineColor(4);
 413     f0->Draw();
 414     f1->Draw("same");
 415     f2->Draw("same");
 416     f3->Draw("same");
 417
 418     TF2 f2D("f2D","AliTPCCorrection::GetCorrSector(x,y,0.9,1,0)-AliTPCCorrection::GetCorrSector(x,y,0.9,1,3)",0,18,85,245);
 419   */
 420   return alignOut;
 421 }
 422
 423
 424 void AliTPCCalibGlobalMisalignment::DumpAlignment( AliTPCCalibGlobalMisalignment* align, TTreeSRedirector *pcstream, const char *name){
 425   //
 426   // Dump alignment per sector into tree
 427   //
 428   TObjArray * array = align->GetAlignSectors();
 429   if (!array) return;
 430   //
 431   //Get mean transformation
 432   TGeoHMatrix matrix;
 433   {for (Int_t isec=0; isec<72; isec++){
 434       TGeoHMatrix* cmatrix=(TGeoHMatrix*)array->At(isec);
 435       TGeoHMatrix* cmatrixDown=(TGeoHMatrix*)array->At(isec%36);
 436       TGeoHMatrix* cmatrixUp=(TGeoHMatrix*)array->At(isec%36+36);
 437       TGeoHMatrix diff(*cmatrixDown);
 438       diff.Multiply(&(cmatrixUp->Inverse()));
 439       (*pcstream)<<name<<
 440         "isec="<<isec<<
 441         "m0.="<<cmatrix<<
 442         "diff.="<<&diff<<
 443         "\n";
 444     }
 445   }
 446
 447 }