EMCAL/AliEMCALSurvey.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 /* $Id: $ */
  17
  18 // Objects of this class read txt file with survey data
  19 // and convert the data into AliAlignObjParams of alignable EMCAL volumes.
  20 // AliEMCALSurvey inherits TObject only to use AliLog "functions".
  21 //
  22 // Dummy functions originally written before EMCAL installation and
  23 // survey are kept for backward compatibility, but now they are not
  24 // used.
  25 //
  26 // Surveyed points on the EMCAL support rails were used with the CATIA
  27 // 3D graphics program to determine the positions of the bottom
  28 // corners of the active area for each supermodule.  These numbers are
  29 // read in from file and converted to position of the center and roll,
  30 // pitch, yaw angles of each installed SM.
  31 //
  32 // J.L. Klay - Cal Poly
  33 // 21-May-2010
  34 //
  35
  36 #include <fstream>
  37
  38 #include <TClonesArray.h>
  39 #include <TGeoManager.h>
  40 #include <TString.h>
  41 #include <TMath.h>
  42
  43 #include "AliSurveyObj.h"
  44 #include "AliSurveyPoint.h"
  45
  46 #include "AliAlignObjParams.h"
  47 #include "AliEMCALGeometry.h"
  48 #include "AliEMCALSurvey.h"
  49 #include "AliLog.h"
  50
  51 ClassImp(AliEMCALSurvey)
  52
  53 //____________________________________________________________________________
  54 AliEMCALSurvey::AliEMCALSurvey()
  55  : fNSuperModule(0),
  56    fSuperModuleData(0),
  57    fDataType(kSurvey)
  58 {
  59   //Default constructor.
  60 }
  61
  62 namespace {
  63
  64   //Coordinates for each SM described in survey reports
  65
  66   struct AliEMCALSuperModuleCoords {
  67     Double_t fX1; //x coordinate of the center of supermodule
  68     Double_t fY1; //y coordinate of the center of supermodule
  69     Double_t fZ1; //z coordinate of the center of supermodule
  70     Double_t fPsi;   //yaw (psi) of supermodule
  71     Double_t fTheta; //pitch (theta) of supermodule
  72     Double_t fPhi;  //roll angle (phi) of supermodule
  73
  74   };
  75
  76 }
  77
  78 //____________________________________________________________________________
  79 AliEMCALSurvey::AliEMCALSurvey(const TString &txtFileName,const SurveyDataType_t type)
  80   : fNSuperModule(0),
  81     fSuperModuleData(0),
  82     fDataType(type)
  83 {
  84   //Get the geometry object and then attempt to
  85   //read survey data from a file, depending on which
  86   //method (kSurvey or kDummy) is selected.
  87
  88   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  89   if (!geom) {
  90     AliError("Cannot obtain AliEMCALGeometry instance.");
  91     return;
  92   }
  93
  94   fNSuperModule = geom->GetNumberOfSuperModules();
  95
  96   if(fDataType == kSurvey) {
  97
  98     AliSurveyObj *s1 = new AliSurveyObj();
  99     s1->FillFromLocalFile(txtFileName);
 100     TObjArray* points = s1->GetData();
 101     InitSuperModuleData(points);
 102
 103   } else {
 104
 105     //Use a dummy file that stores x,y,z of the center of each SM
 106     //useful for testing...
 107     std::ifstream inputFile(txtFileName.Data());
 108     if (!inputFile) {
 109       AliError(("Cannot open the survey file " + txtFileName).Data());
 110       return;
 111     }
 112     Int_t dummyInt = 0;
 113     Double_t *xReal = new Double_t[fNSuperModule];
 114     Double_t *yReal = new Double_t[fNSuperModule];
 115     Double_t *zReal = new Double_t[fNSuperModule];
 116     Double_t *psiReal = new Double_t[fNSuperModule];
 117     Double_t *thetaReal = new Double_t[fNSuperModule];
 118     Double_t *phiReal = new Double_t[fNSuperModule];
 119
 120     for (Int_t i = 0; i < fNSuperModule; ++i) {
 121       if (!inputFile) {
 122         AliError("Error while reading input file.");
 123         delete [] xReal;
 124         delete [] yReal;
 125         delete [] zReal;
 126         delete [] psiReal;
 127         delete [] thetaReal;
 128         delete [] phiReal;
 129         return;
 130       }
 131       inputFile>>dummyInt>>xReal[i]>>yReal[i]>>zReal[i]>>psiReal[i]>>thetaReal[i]>>phiReal[i];
 132     }
 133
 134     InitSuperModuleData(xReal, yReal, zReal, psiReal, thetaReal, phiReal);
 135
 136     delete [] xReal;
 137     delete [] yReal;
 138     delete [] zReal;
 139     delete [] psiReal;
 140     delete [] thetaReal;
 141     delete [] phiReal;
 142
 143   } //kDummy way of doing it
 144
 145 }
 146
 147 //____________________________________________________________________________
 148 AliEMCALSurvey::~AliEMCALSurvey()
 149 {
 150   //destructor
 151   delete [] fSuperModuleData;
 152 }
 153
 154 //____________________________________________________________________________
 155 void AliEMCALSurvey::CreateAliAlignObjParams(TClonesArray &array)
 156 {
 157   //Create AliAlignObjParams.
 158   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
 159   if (!geom) {
 160     AliError("Cannot obtain AliEMCALGeometry instance.");
 161     return;
 162   }
 163
 164   if (!gGeoManager) {
 165     AliWarning("Cannot create local transformations for supermodules - gGeoManager does not exist.");
 166     AliInfo("Null shifts and rotations will be created instead.");
 167     return CreateNullObjects(array, geom);
 168   }
 169
 170   Int_t arrayInd = array.GetEntries(), iIndex = 0;
 171   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 172   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
 173
 174   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
 175     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
 176     if(geom->GetKey110DEG() && smodnum >= 10) {
 177       smodName = "EMCAL/HalfSupermodule";
 178       smodName += (smodnum-10+1);
 179     }
 180     AliEMCALSuperModuleDelta t(GetSuperModuleTransformation(smodnum));
 181     new(array[arrayInd])
 182       AliAlignObjParams(
 183                         smodName.Data(), volid,
 184                         t.fXShift, t.fYShift, t.fZShift,
 185                         -t.fPsi, -t.fTheta, -t.fPhi,
 186                         false
 187                         );
 188     ++arrayInd;
 189
 190   }
 191
 192 }
 193
 194 //____________________________________________________________________________
 195 void AliEMCALSurvey::CreateNullObjects(TClonesArray &array, const AliEMCALGeometry *geom)const
 196 {
 197   //Create null shifts and rotations.
 198   Int_t arrayInd = array.GetEntries(), iIndex = 0;
 199   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 200   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
 201
 202   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
 203     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
 204     if(geom->GetKey110DEG() && smodnum >= 10) {
 205       smodName = "EMCAL/HalfSupermodule";
 206       smodName += (smodnum-10+1);
 207     }
 208     new(array[arrayInd]) AliAlignObjParams(smodName.Data(), volid, 0., 0., 0., 0., 0., 0., true);
 209     ++arrayInd;
 210   }
 211 }
 212
 213 //____________________________________________________________________________
 214 AliEMCALSurvey::AliEMCALSuperModuleDelta AliEMCALSurvey::GetSuperModuleTransformation(Int_t supModIndex)const
 215 {
 216   //Supermodule transformation.
 217   AliEMCALSuperModuleDelta t = {0., 0., 0., 0., 0., 0.};
 218   if (!fSuperModuleData)
 219     return t;
 220
 221   return fSuperModuleData[supModIndex];
 222 }
 223
 224 //____________________________________________________________________________
 225 void AliEMCALSurvey::InitSuperModuleData(const TObjArray *svypts)
 226 {
 227   //This method uses the data points from the EMCAL survey and CATIA program to
 228   //create the alignment matrices.  Only valid for (installed)
 229   //SM, others will have null objects
 230
 231   /*--------------------------------------
 232     The bottom edges of the strip modules
 233     define the active area of the EMCAL, but
 234     in software we have a box to hold them which
 235     is longer than that.  We need to convert
 236     info about the position of the corners of the
 237     bottom of the active area to the center of
 238     the software box that contains the strip
 239     modules.
 240
 241     View from beam axis up to EMCAL
 242               Ai                Ci
 243
 244         0,1         0,0 1,0          1,1
 245     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
 246     x   x             x x              x   x
 247     x   x    % *      x x      * %     x   x
 248     x   x             x x              x   x
 249     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
 250         1,1         1,0 0,0          0,1
 251     <--> = added length                 <--> = added length
 252
 253     * represents the center of the active area
 254     % represents the center of the full box (with added length)
 255
 256          View from side of topmost SM
 257
 258               Ai                Ci
 259
 260         0,1         0,0 1,0          1,1
 261     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
 262     x   x    % *      x x      % *     x   x
 263     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
 264         1,1         1,0 0,0          0,1
 265     <--> = added length                 <--> = added length
 266
 267     * represents the center of the active area
 268     % represents the center of the full box (with added length)
 269
 270     -------------------------------------*/
 271
 272   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 273   //Center of supermodules
 274   Float_t *pars = geom->GetSuperModulesPars();
 275   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
 276   Double_t phi, phiRad, xpos, ypos, zpos;
 277
 278   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
 279   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
 280     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
 281     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
 282     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
 283     xpos = rpos * TMath::Cos(phiRad);
 284     ypos = rpos * TMath::Sin(phiRad);
 285     zpos = pars[2];
 286     if(geom->GetKey110DEG() && smodnum >= 10) {
 287       xpos += (pars[1]/2. * TMath::Sin(phiRad));
 288       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
 289     }
 290     smc.fX1 = xpos;
 291     smc.fY1 = ypos;
 292     smc.fPhi = phi; //degrees
 293     smc.fTheta = 0.; //degrees
 294     smc.fPsi = 0.; //degrees
 295     if(smodnum%2==0) {
 296       smc.fZ1 = zpos;
 297     } else {
 298       smc.fZ1 = -zpos;
 299     }
 300
 301     //printf("PHI OF IDEAL SM = %.2f\n",smc.fPhi);
 302
 303   }
 304
 305   //Real coordinates of center and rotation angles need to be computed
 306   //from the survey/CATIA points
 307
 308   char substr[100];
 309   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
 310   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
 311     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
 312     Double_t zLength = pars[2]*2.; //length of SM in z from software
 313     Double_t halfHeight = pars[0]; //half the height of the SM in y direction
 314
 315     sprintf(substr,"4096%d",smodnum);
 316     //retrieve components of four face points and determine average position of center
 317     //in x,y,z
 318
 319     std::vector<Double_t> xval;
 320     std::vector<Double_t> yval;
 321     std::vector<Double_t> zval;
 322
 323     for(Int_t i = 0; i < svypts->GetEntries(); i++) {
 324       AliSurveyPoint* pt = (AliSurveyPoint*)svypts->At(i);
 325       TString ptname = pt->GetPointName();
 326       if(ptname.Contains(substr)) {
 327         //Note: order of values is 00, 01, 10, 11
 328         xval.push_back(pt->GetX()*100.); //convert m to cm
 329         yval.push_back(pt->GetY()*100.);
 330         zval.push_back(pt->GetZ()*100.);
 331       }
 332     }
 333
 334     //compute center of active area of each SM on bottome face from survey points
 335     Double_t activeX = ((xval[0] + (xval[2] - xval[0])/2.)        //x00 and x10
 336                         +(xval[1] + (xval[3] - xval[1])/2.) ) /2.; //x01 and x11
 337
 338     Double_t activeY = ((yval[0] + (yval[2] - yval[0])/2.)
 339                         +(yval[1] + (yval[3] - yval[1])/2.) ) /2.;
 340
 341     Double_t activeZ = ((zval[0] + (zval[2] - zval[0])/2.)
 342                         +(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
 343
 344     //printf("Bottom Center of active area of SM %s: %.2f, %.2f, %.2f\n",substr,activeX,activeY,activeZ);
 345
 346     //compute angles for each SM
 347     //rotation about each axis
 348     //phi = angle in x-y plane
 349
 350     Double_t realphi = 0.;
 351     //Note: this is phi wrt y axis.  To get phi wrt to x, add pi/2
 352     if(smodnum%2 == 0) {
 353       realphi = (TMath::ATan((yval[2] - yval[0])/(xval[2] - xval[0]))
 354                  +TMath::ATan((yval[3] - yval[1])/(xval[3] - xval[1])) )/2.;
 355     } else {
 356       realphi = (TMath::ATan((yval[0] - yval[2])/(xval[0] - xval[2]))
 357                  +TMath::ATan((yval[1] - yval[3])/(xval[1] - xval[3])) )/2.;
 358     }
 359
 360     //NOTE: Psi angle is always zero because the two z values being
 361     //subtracted are exactly the same, but just in case that could change...
 362     //psi = angle in x-z plane
 363     Double_t realpsi = (TMath::ATan((zval[0] - zval[2])/(xval[2] - xval[0]))
 364                         +TMath::ATan((zval[1] - zval[3])/(xval[3] - xval[1])) )/2.;
 365
 366     //theta = angle in y-z plane
 367     Double_t realtheta = TMath::Pi()/2.
 368                        - (TMath::ATan((zval[2] - zval[3])/(yval[3] - yval[2]))
 369                           +TMath::ATan((zval[0] - zval[1])/(yval[1] - yval[0])) )/2.;
 370
 371     //printf("Old edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);
 372     //printf("Old edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
 373
 374     //Now calculate the center of the box in z with length added to the 01
 375     //and 11 corners, corrected by the theta angle
 376     Double_t activeLength = TMath::Abs(((zval[1] - zval[0]) + (zval[3] - zval[2]))/2.);
 377     //printf("ACTIVE LENGTH = %.2f\n",activeLength);
 378     if(smodnum%2 == 0) {
 379       yval[1] += (zLength - activeLength)*sin(realtheta);
 380       yval[3] += (zLength - activeLength)*sin(realtheta);
 381       zval[1] += (zLength - activeLength)*cos(realtheta);
 382       zval[3] += (zLength - activeLength)*cos(realtheta);
 383     } else {
 384       yval[1] -= (zLength - activeLength)*sin(realtheta);
 385       yval[3] -= (zLength - activeLength)*sin(realtheta);
 386       zval[1] -= (zLength - activeLength)*cos(realtheta);
 387       zval[3] -= (zLength - activeLength)*cos(realtheta);
 388     }
 389
 390     //printf("New extended edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);
 391     //printf("New extended edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
 392
 393     //Compute the center of the bottom of the box in x,y,z
 394     Double_t realX = activeX;
 395     Double_t realY = ((yval[0] + (yval[2] - yval[0])/2.)
 396                         +(yval[1] + (yval[3] - yval[1])/2.) ) /2.;
 397     Double_t realZ = ((zval[0] + (zval[2] - zval[0])/2.)
 398                         +(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
 399
 400
 401     //printf("Bottom Center of SM %s Box: %.2f, %.2f, %.2f\n",substr,realX,realY,realZ);
 402
 403     //correct the SM centers so that we have the center of the box in
 404     //x,y using the phi,theta angles
 405     realX += halfHeight*TMath::Cos(TMath::Pi()/2+realphi);
 406     realY += halfHeight*(TMath::Sin(TMath::Pi()/2+realphi) + TMath::Sin(realtheta));
 407     realZ += halfHeight*TMath::Cos(TMath::Pi()/2-realtheta);
 408
 409     //printf("Rotation angles of SM %s: %.4f, %.4f, %.4f\n",substr,realphi*TMath::RadToDeg(),realpsi*TMath::RadToDeg(),realtheta*TMath::RadToDeg());
 410     //printf("Middle of SM %s: %.2f, %.2f, %.2f\n\n",substr,realX,realY,realZ);
 411
 412     smc.fX1 = realX;
 413     smc.fY1 = realY;
 414     smc.fZ1 = realZ;
 415
 416     smc.fPhi = 90. + realphi*TMath::RadToDeg();
 417     smc.fTheta = 0. + realtheta*TMath::RadToDeg();
 418     smc.fPsi = 0. + realpsi*TMath::RadToDeg();
 419
 420   }//loop over supermodules
 421
 422   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
 423
 424   for (Int_t i = 0; i < fNSuperModule; ++i) {
 425     const AliEMCALSuperModuleCoords &real = realSM[i];
 426     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
 427     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
 428     t.fXShift = real.fX1 - ideal.fX1;
 429     t.fYShift = real.fY1 - ideal.fY1;
 430     t.fZShift = real.fZ1 - ideal.fZ1;
 431     t.fPhi = real.fPhi - ideal.fPhi;
 432     t.fTheta = real.fTheta - ideal.fTheta;
 433     t.fPsi = real.fPsi - ideal.fPsi;
 434
 435     printf("===================== SM %d =======================\n",i);
 436     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
 437     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
 438     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
 439     printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
 440     printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
 441     printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
 442     printf("===================================================\n");
 443   }
 444
 445   delete [] realSM;
 446   delete [] idealSM;
 447 }
 448
 449
 450 //____________________________________________________________________________
 451 void AliEMCALSurvey::InitSuperModuleData(const Double_t *xReal, const Double_t *yReal,
 452                                          const Double_t *zReal, const Double_t *psiReal,
 453                                          const Double_t *thetaReal, const Double_t *phiReal)
 454 {
 455   ///////////////////////////////////////
 456   //Dummy method just takes the inputted values and applies them
 457   //
 458   //Useful for testing small changes
 459   //////////////////////////////////////
 460   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 461   //Center of supermodules
 462   Float_t *pars = geom->GetSuperModulesPars();
 463   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
 464   Double_t phi, phiRad, xpos, ypos, zpos;
 465
 466   zpos = pars[2];
 467
 468   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
 469   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
 470     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
 471     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
 472     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
 473     xpos = rpos * TMath::Cos(phiRad);
 474     ypos = rpos * TMath::Sin(phiRad);
 475     if(geom->GetKey110DEG() && smodnum >= 10) {
 476       xpos += (pars[1]/2. * TMath::Sin(phiRad));
 477       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
 478     }
 479     smc.fX1 = xpos;
 480     smc.fY1 = ypos;
 481
 482     smc.fPhi = phi; //degrees
 483     smc.fTheta = 0.; //degrees
 484     smc.fPsi = 0.; //degrees
 485
 486     if(smodnum%2==0) {
 487       smc.fZ1 = zpos;
 488     } else {
 489       smc.fZ1 = -zpos;
 490     }
 491
 492   }
 493
 494   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
 495   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
 496     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
 497     smc.fX1    = xReal[smodnum];
 498     smc.fY1    = yReal[smodnum];
 499     smc.fZ1    = zReal[smodnum];
 500     smc.fTheta = thetaReal[smodnum];
 501     smc.fPsi   = psiReal[smodnum];
 502     smc.fPhi   = phiReal[smodnum];
 503   }
 504
 505   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
 506
 507   for (Int_t i = 0; i < fNSuperModule; ++i) {
 508     const AliEMCALSuperModuleCoords &real = realSM[i];
 509     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
 510     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
 511     t.fTheta = real.fTheta - ideal.fTheta;
 512     t.fPsi = 0.;
 513     t.fPhi = real.fPhi - ideal.fPhi;
 514     t.fXShift = real.fX1 - ideal.fX1;
 515     t.fYShift = real.fY1 - ideal.fY1;
 516     t.fZShift = real.fZ1 - ideal.fZ1;
 517
 518     printf("===================== SM %d =======================\n",i);
 519     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
 520     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
 521     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
 522     printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
 523     printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
 524     printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
 525     printf("===================================================\n");
 526   }
 527
 528   delete [] realSM;
 529   delete [] idealSM;
 530 }
 531