[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSurvey.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.                  *
 **************************************************************************/

/* $Id: $ */

// Objects of this class read txt file with survey data
// and convert the data into AliAlignObjParams of alignable EMCAL volumes.
// AliEMCALSurvey inherits TObject only to use AliLog "functions".
//
// Dummy functions originally written before EMCAL installation and
// survey are kept for backward compatibility, but now they are not
// used.
//
// Surveyed points on the EMCAL support rails were used with the CATIA
// 3D graphics program to determine the positions of the bottom
// corners of the active area for each supermodule.  These numbers are
// read in from file and converted to position of the center and roll, 
// pitch, yaw angles of each installed SM.
//
// J.L. Klay - Cal Poly
// 21-May-2010
//

#include <fstream>

#include <TClonesArray.h>
#include <TGeoManager.h>
#include <TString.h>
#include <TMath.h>

#include "AliSurveyObj.h"
#include "AliSurveyPoint.h"

#include "AliAlignObjParams.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALSurvey.h"
#include "AliLog.h"

ClassImp(AliEMCALSurvey)

//____________________________________________________________________________
AliEMCALSurvey::AliEMCALSurvey()
 : fNSuperModule(0),
   fSuperModuleData(0),
   fDataType(kSurvey)
{
  //Default constructor.
}

namespace {

  //Coordinates for each SM described in survey reports

  struct AliEMCALSuperModuleCoords {
    Double_t fX1; //x coordinate of the center of supermodule
    Double_t fY1; //y coordinate of the center of supermodule
    Double_t fZ1; //z coordinate of the center of supermodule
    Double_t fPsi;   //yaw (psi) of supermodule
    Double_t fTheta; //pitch (theta) of supermodule
    Double_t fPhi;  //roll angle (phi) of supermodule

  };

}

//____________________________________________________________________________
AliEMCALSurvey::AliEMCALSurvey(const TString &txtFileName,const SurveyDataType_t type)
  : fNSuperModule(0),
    fSuperModuleData(0),
    fDataType(type)
{
  //Get the geometry object and then attempt to
  //read survey data from a file, depending on which
  //method (kSurvey or kDummy) is selected.

  const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  if (!geom) {
    AliError("Cannot obtain AliEMCALGeometry instance.");
    return;
  }

  fNSuperModule = geom->GetNumberOfSuperModules();

  if(fDataType == kSurvey) {

    AliSurveyObj *s1 = new AliSurveyObj();
    s1->FillFromLocalFile(txtFileName);
    TObjArray* points = s1->GetData();
    InitSuperModuleData(points);

  } else {

    //Use a dummy file that stores x,y,z of the center of each SM
    //useful for testing...
    std::ifstream inputFile(txtFileName.Data());
    if (!inputFile) {
      AliError(("Cannot open the survey file " + txtFileName).Data());
      return;
    }
    Int_t dummyInt = 0;
    Double_t *xReal = new Double_t[fNSuperModule];
    Double_t *yReal = new Double_t[fNSuperModule];
    Double_t *zReal = new Double_t[fNSuperModule];
    Double_t *psiReal = new Double_t[fNSuperModule];
    Double_t *thetaReal = new Double_t[fNSuperModule];
    Double_t *phiReal = new Double_t[fNSuperModule];

    for (Int_t i = 0; i < fNSuperModule; ++i) {
      if (!inputFile) {
	AliError("Error while reading input file.");
	delete [] xReal;
	delete [] yReal;
	delete [] zReal;
	delete [] psiReal;
	delete [] thetaReal;
	delete [] phiReal;
	return;
      }
      inputFile>>dummyInt>>xReal[i]>>yReal[i]>>zReal[i]>>psiReal[i]>>thetaReal[i]>>phiReal[i];
    }

    InitSuperModuleData(xReal, yReal, zReal, psiReal, thetaReal, phiReal);

    delete [] xReal;
    delete [] yReal;
    delete [] zReal;
    delete [] psiReal;
    delete [] thetaReal;
    delete [] phiReal;

  } //kDummy way of doing it

}

//____________________________________________________________________________
AliEMCALSurvey::~AliEMCALSurvey()
{
  //destructor
  delete [] fSuperModuleData;
}

//____________________________________________________________________________
void AliEMCALSurvey::CreateAliAlignObjParams(TClonesArray &array)
{
  //Create AliAlignObjParams.
  const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
  if (!geom) {
    AliError("Cannot obtain AliEMCALGeometry instance.");
    return;
  }

  if (!gGeoManager) {
    AliWarning("Cannot create local transformations for supermodules - gGeoManager does not exist.");
    AliInfo("Null shifts and rotations will be created instead.");
    return CreateNullObjects(array, geom);
  }

  Int_t arrayInd = array.GetEntries(), iIndex = 0;
  AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
  UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);

  for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
    TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
    if(geom->GetKey110DEG() && smodnum >= 10) {
      smodName = "EMCAL/HalfSupermodule";
      smodName += (smodnum-10+1);
    }    
    AliEMCALSuperModuleDelta t(GetSuperModuleTransformation(smodnum));
    new(array[arrayInd])
      AliAlignObjParams(
			smodName.Data(), volid, 
			t.fXShift, t.fYShift, t.fZShift, 
			-t.fPsi, -t.fTheta, -t.fPhi, 
			false
			);
    ++arrayInd;

  }

}

//____________________________________________________________________________
void AliEMCALSurvey::CreateNullObjects(TClonesArray &array, const AliEMCALGeometry *geom)const
{
  //Create null shifts and rotations.
  Int_t arrayInd = array.GetEntries(), iIndex = 0;
  AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
  UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);

  for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
    TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
    if(geom->GetKey110DEG() && smodnum >= 10) {
      smodName = "EMCAL/HalfSupermodule";
      smodName += (smodnum-10+1);
    }
    new(array[arrayInd]) AliAlignObjParams(smodName.Data(), volid, 0., 0., 0., 0., 0., 0., true);
    ++arrayInd;
  }
}

//____________________________________________________________________________
AliEMCALSurvey::AliEMCALSuperModuleDelta AliEMCALSurvey::GetSuperModuleTransformation(Int_t supModIndex)const
{
  //Supermodule transformation.
  AliEMCALSuperModuleDelta t = {0., 0., 0., 0., 0., 0.};
  if (!fSuperModuleData)
    return t;

  return fSuperModuleData[supModIndex];
}

//____________________________________________________________________________
void AliEMCALSurvey::InitSuperModuleData(const TObjArray *svypts)
{
  //This method uses the data points from the EMCAL survey and CATIA program to
  //create the alignment matrices.  Only valid for (installed)
  //SM, others will have null objects

  /*--------------------------------------
    The bottom edges of the strip modules
    define the active area of the EMCAL, but
    in software we have a box to hold them which
    is longer than that.  We need to convert
    info about the position of the corners of the
    bottom of the active area to the center of
    the software box that contains the strip
    modules.

    View from beam axis up to EMCAL
              Ai                Ci

        0,1         0,0 1,0          1,1
    xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
    x   x             x x              x   x
    x   x    % *      x x      * %     x   x
    x   x             x x              x   x
    xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
        1,1         1,0 0,0          0,1
    <--> = added length                 <--> = added length

    * represents the center of the active area
    % represents the center of the full box (with added length)

         View from side of topmost SM

              Ai                Ci

        0,1         0,0 1,0          1,1
    xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
    x   x    % *      x x      % *     x   x
    xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
        1,1         1,0 0,0          0,1
    <--> = added length                 <--> = added length

    * represents the center of the active area
    % represents the center of the full box (with added length)

    -------------------------------------*/

  AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  //Center of supermodules
  Float_t *pars = geom->GetSuperModulesPars();
  Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
  Double_t phi, phiRad, xpos, ypos, zpos;

  AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
  for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
    AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
    phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
    phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
    xpos = rpos * TMath::Cos(phiRad);
    ypos = rpos * TMath::Sin(phiRad);
    zpos = pars[2];
    if(geom->GetKey110DEG() && smodnum >= 10) {
      xpos += (pars[1]/2. * TMath::Sin(phiRad));
      ypos -= (pars[1]/2. * TMath::Cos(phiRad));
    }
    smc.fX1 = xpos;
    smc.fY1 = ypos;
    smc.fPhi = phi; //degrees
    smc.fTheta = 0.; //degrees
    smc.fPsi = 0.; //degrees
    if(smodnum%2==0) {
      smc.fZ1 = zpos;
    } else {
      smc.fZ1 = -zpos;
    }

    //printf("PHI OF IDEAL SM = %.2f\n",smc.fPhi);

  }

  //Real coordinates of center and rotation angles need to be computed
  //from the survey/CATIA points

  char substr[100];
  AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
  for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
    AliEMCALSuperModuleCoords &smc = realSM[smodnum];
    Double_t zLength = pars[2]*2.; //length of SM in z from software
    Double_t halfHeight = pars[0]; //half the height of the SM in y direction

    sprintf(substr,"4096%d",smodnum);
    //retrieve components of four face points and determine average position of center
    //in x,y,z

    std::vector<Double_t> xval;
    std::vector<Double_t> yval;
    std::vector<Double_t> zval;
    
    for(Int_t i = 0; i < svypts->GetEntries(); i++) {
      AliSurveyPoint* pt = (AliSurveyPoint*)svypts->At(i);
      TString ptname = pt->GetPointName();
      if(ptname.Contains(substr)) {
	//Note: order of values is 00, 01, 10, 11
	xval.push_back(pt->GetX()*100.); //convert m to cm
	yval.push_back(pt->GetY()*100.); 
	zval.push_back(pt->GetZ()*100.); 
      }
    }

    //compute center of active area of each SM on bottome face from survey points
    Double_t activeX = ((xval[0] + (xval[2] - xval[0])/2.)        //x00 and x10
			+(xval[1] + (xval[3] - xval[1])/2.) ) /2.; //x01 and x11
    
    Double_t activeY = ((yval[0] + (yval[2] - yval[0])/2.)
			+(yval[1] + (yval[3] - yval[1])/2.) ) /2.;
    
    Double_t activeZ = ((zval[0] + (zval[2] - zval[0])/2.)
			+(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
    
    //printf("Bottom Center of active area of SM %s: %.2f, %.2f, %.2f\n",substr,activeX,activeY,activeZ);
    
    //compute angles for each SM
    //rotation about each axis
    //phi = angle in x-y plane
    
    Double_t realphi = 0.;
    //Note: this is phi wrt y axis.  To get phi wrt to x, add pi/2
    if(smodnum%2 == 0) {
      realphi = (TMath::ATan((yval[2] - yval[0])/(xval[2] - xval[0])) 
		 +TMath::ATan((yval[3] - yval[1])/(xval[3] - xval[1])) )/2.;
    } else {
      realphi = (TMath::ATan((yval[0] - yval[2])/(xval[0] - xval[2]))
		 +TMath::ATan((yval[1] - yval[3])/(xval[1] - xval[3])) )/2.;
    }

    //NOTE: Psi angle is always zero because the two z values being
    //subtracted are exactly the same, but just in case that could change...
    //psi = angle in x-z plane
    Double_t realpsi = (TMath::ATan((zval[0] - zval[2])/(xval[2] - xval[0]))
			+TMath::ATan((zval[1] - zval[3])/(xval[3] - xval[1])) )/2.;
    
    //theta = angle in y-z plane
    Double_t realtheta = TMath::Pi()/2. 
                       - (TMath::ATan((zval[2] - zval[3])/(yval[3] - yval[2]))
			  +TMath::ATan((zval[0] - zval[1])/(yval[1] - yval[0])) )/2.;

    //printf("Old edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);
    //printf("Old edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);

    //Now calculate the center of the box in z with length added to the 01
    //and 11 corners, corrected by the theta angle
    Double_t activeLength = TMath::Abs(((zval[1] - zval[0]) + (zval[3] - zval[2]))/2.);
    //printf("ACTIVE LENGTH = %.2f\n",activeLength);
    if(smodnum%2 == 0) {
      yval[1] += (zLength - activeLength)*sin(realtheta);
      yval[3] += (zLength - activeLength)*sin(realtheta);
      zval[1] += (zLength - activeLength)*cos(realtheta);
      zval[3] += (zLength - activeLength)*cos(realtheta);
    } else {
      yval[1] -= (zLength - activeLength)*sin(realtheta);
      yval[3] -= (zLength - activeLength)*sin(realtheta);
      zval[1] -= (zLength - activeLength)*cos(realtheta);
      zval[3] -= (zLength - activeLength)*cos(realtheta);
    }

    //printf("New extended edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);            
    //printf("New extended edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);

    //Compute the center of the bottom of the box in x,y,z
    Double_t realX = activeX;    
    Double_t realY = ((yval[0] + (yval[2] - yval[0])/2.)
			+(yval[1] + (yval[3] - yval[1])/2.) ) /2.;    
    Double_t realZ = ((zval[0] + (zval[2] - zval[0])/2.)
			+(zval[1] + (zval[3] - zval[1])/2.) ) /2.;


    //printf("Bottom Center of SM %s Box: %.2f, %.2f, %.2f\n",substr,realX,realY,realZ);

    //correct the SM centers so that we have the center of the box in
    //x,y using the phi,theta angles                   
    realX += halfHeight*TMath::Cos(TMath::Pi()/2+realphi);
    realY += halfHeight*(TMath::Sin(TMath::Pi()/2+realphi) + TMath::Sin(realtheta));
    realZ += halfHeight*TMath::Cos(TMath::Pi()/2-realtheta);

    //printf("Rotation angles of SM %s: %.4f, %.4f, %.4f\n",substr,realphi*TMath::RadToDeg(),realpsi*TMath::RadToDeg(),realtheta*TMath::RadToDeg());
    //printf("Middle of SM %s: %.2f, %.2f, %.2f\n\n",substr,realX,realY,realZ);

    smc.fX1 = realX;
    smc.fY1 = realY;
    smc.fZ1 = realZ;

    smc.fPhi = 90. + realphi*TMath::RadToDeg();
    smc.fTheta = 0. + realtheta*TMath::RadToDeg();
    smc.fPsi = 0. + realpsi*TMath::RadToDeg();

  }//loop over supermodules
  
  fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
  
  for (Int_t i = 0; i < fNSuperModule; ++i) {
    const AliEMCALSuperModuleCoords &real = realSM[i];
    const AliEMCALSuperModuleCoords &ideal = idealSM[i];
    AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
    t.fXShift = real.fX1 - ideal.fX1;
    t.fYShift = real.fY1 - ideal.fY1;
    t.fZShift = real.fZ1 - ideal.fZ1;
    t.fPhi = real.fPhi - ideal.fPhi;
    t.fTheta = real.fTheta - ideal.fTheta;
    t.fPsi = real.fPsi - ideal.fPsi;

    printf("===================== SM %d =======================\n",i);
    printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
    printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
    printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
    printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
    printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
    printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
    printf("===================================================\n");    
  }

  delete [] realSM;
  delete [] idealSM;
}


//____________________________________________________________________________
void AliEMCALSurvey::InitSuperModuleData(const Double_t *xReal, const Double_t *yReal, 
					 const Double_t *zReal, const Double_t *psiReal,
					 const Double_t *thetaReal, const Double_t *phiReal)
{
  ///////////////////////////////////////
  //Dummy method just takes the inputted values and applies them
  //
  //Useful for testing small changes
  //////////////////////////////////////
  AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
  //Center of supermodules
  Float_t *pars = geom->GetSuperModulesPars();
  Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
  Double_t phi, phiRad, xpos, ypos, zpos;

  zpos = pars[2];

  AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
  for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
    AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
    phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
    phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
    xpos = rpos * TMath::Cos(phiRad);
    ypos = rpos * TMath::Sin(phiRad);
    if(geom->GetKey110DEG() && smodnum >= 10) {
      xpos += (pars[1]/2. * TMath::Sin(phiRad));
      ypos -= (pars[1]/2. * TMath::Cos(phiRad));
    }
    smc.fX1 = xpos;
    smc.fY1 = ypos;

    smc.fPhi = phi; //degrees
    smc.fTheta = 0.; //degrees
    smc.fPsi = 0.; //degrees

    if(smodnum%2==0) {
      smc.fZ1 = zpos;
    } else {
      smc.fZ1 = -zpos;
    }

  }

  AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
  for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
    AliEMCALSuperModuleCoords &smc = realSM[smodnum];
    smc.fX1    = xReal[smodnum];  
    smc.fY1    = yReal[smodnum];  
    smc.fZ1    = zReal[smodnum];  
    smc.fTheta = thetaReal[smodnum];
    smc.fPsi   = psiReal[smodnum];
    smc.fPhi   = phiReal[smodnum];
  }
  
  fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
  
  for (Int_t i = 0; i < fNSuperModule; ++i) {
    const AliEMCALSuperModuleCoords &real = realSM[i];
    const AliEMCALSuperModuleCoords &ideal = idealSM[i];
    AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
    t.fTheta = real.fTheta - ideal.fTheta;
    t.fPsi = 0.;
    t.fPhi = real.fPhi - ideal.fPhi;
    t.fXShift = real.fX1 - ideal.fX1;
    t.fYShift = real.fY1 - ideal.fY1;
    t.fZShift = real.fZ1 - ideal.fZ1;

    printf("===================== SM %d =======================\n",i);
    printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
    printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
    printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
    printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
    printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
    printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
    printf("===================================================\n");    
  }

  delete [] realSM;
  delete [] idealSM;
}
Commit	Line	Data
289cc8a7	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	//
68711871	22	// Dummy functions originally written before EMCAL installation and
	23	// survey are kept for backward compatibility, but now they are not
	24	// used.
54a4c3a6	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.
289cc8a7	31	//
289cc8a7	32	// J.L. Klay - Cal Poly
54a4c3a6	33	// 21-May-2010
289cc8a7	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"
68711871	44	#include "AliSurveyPoint.h"
289cc8a7	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()
54a4c3a6	55	: fNSuperModule(0),
	56	fSuperModuleData(0),
	57	fDataType(kSurvey)
289cc8a7	58	{
	59	//Default constructor.
	60	}
	61
	62	namespace {
	63
68711871	64	//Coordinates for each SM described in survey reports
289cc8a7	65
289cc8a7	66	struct AliEMCALSuperModuleCoords {
68711871	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
68711871	70	Double_t fPsi; //yaw (psi) of supermodule
54a4c3a6	71	Double_t fTheta; //pitch (theta) of supermodule
	72	Double_t fPhi; //roll angle (phi) of supermodule
	73
289cc8a7	74	};
	75
	76	}
	77
	78	//____________________________________________________________________________
54a4c3a6	79	AliEMCALSurvey::AliEMCALSurvey(const TString &txtFileName,const SurveyDataType_t type)
	80	: fNSuperModule(0),
	81	fSuperModuleData(0),
	82	fDataType(type)
289cc8a7	83	{
54a4c3a6	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
289cc8a7	88	const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
	89	if (!geom) {
	90	AliError("Cannot obtain AliEMCALGeometry instance.");
	91	return;
	92	}
	93
54a4c3a6	94	fNSuperModule = geom->GetNumberOfSuperModules();
289cc8a7	95
54a4c3a6	96	if(fDataType == kSurvey) {
289cc8a7	97
54a4c3a6	98	AliSurveyObj *s1 = new AliSurveyObj();
	99	s1->FillFromLocalFile(txtFileName);
	100	TObjArray* points = s1->GetData();
	101	InitSuperModuleData(points);
289cc8a7	102
54a4c3a6	103	} else {
68711871	104
54a4c3a6	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
289cc8a7	144
289cc8a7	145	}
	146
	147	//____________________________________________________________________________
	148	AliEMCALSurvey::~AliEMCALSurvey()
	149	{
68711871	150	//destructor
289cc8a7	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;
68711871	189
289cc8a7	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	//____________________________________________________________________________
68711871	225	void AliEMCALSurvey::InitSuperModuleData(const TObjArray *svypts)
289cc8a7	226	{
54a4c3a6	227	//This method uses the data points from the EMCAL survey and CATIA program to
54a4c3a6	228	//create the alignment matrices. Only valid for (installed)
68711871	229	//SM, others will have null objects
54a4c3a6	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
68711871	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
54a4c3a6	301	//printf("PHI OF IDEAL SM = %.2f\n",smc.fPhi);
68711871	302
	303	}
	304
	305	//Real coordinates of center and rotation angles need to be computed
54a4c3a6	306	//from the survey/CATIA points
68711871	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];
54a4c3a6	312	Double_t zLength = pars[2]*2.; //length of SM in z from software
54a4c3a6	313	Double_t halfHeight = pars[0]; //half the height of the SM in y direction
68711871	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)) {
54a4c3a6	327	//Note: order of values is 00, 01, 10, 11
68711871	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
54a4c3a6	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.;
68711871	358	}
68711871	359
54a4c3a6	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);
68711871	388	}
68711871	389
54a4c3a6	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);
68711871	402
54a4c3a6	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);
68711871	408
54a4c3a6	409	//printf("Rotation angles of SM %s: %.4f, %.4f, %.4f\n",substr,realphiTMath::RadToDeg(),realpsiTMath::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();
68711871	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	}
289cc8a7	444
68711871	445	delete [] realSM;
	446	delete [] idealSM;
	447	}
	448
	449
	450	//____________________________________________________________________________
54a4c3a6	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)
68711871	454	{
68711871	455	///////////////////////////////////////
54a4c3a6	456	//Dummy method just takes the inputted values and applies them
	457	//
	458	//Useful for testing small changes
68711871	459	//////////////////////////////////////
289cc8a7	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
54a4c3a6	466	zpos = pars[2];
54a4c3a6	467
289cc8a7	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);
289cc8a7	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;
54a4c3a6	481
	482	smc.fPhi = phi; //degrees
	483	smc.fTheta = 0.; //degrees
	484	smc.fPsi = 0.; //degrees
	485
289cc8a7	486	if(smodnum%2==0) {
	487	smc.fZ1 = zpos;
	488	} else {
	489	smc.fZ1 = -zpos;
	490	}
	491
	492	}
3ddecb84	493
289cc8a7	494	AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
	495	for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
	496	AliEMCALSuperModuleCoords &smc = realSM[smodnum];
54a4c3a6	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];
289cc8a7	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];
54a4c3a6	511	t.fTheta = real.fTheta - ideal.fTheta;
289cc8a7	512	t.fPsi = 0.;
54a4c3a6	513	t.fPhi = real.fPhi - ideal.fPhi;
289cc8a7	514	t.fXShift = real.fX1 - ideal.fX1;
	515	t.fYShift = real.fY1 - ideal.fY1;
	516	t.fZShift = real.fZ1 - ideal.fZ1;
	517
68711871	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);
54a4c3a6	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);
68711871	525	printf("===================================================\n");
289cc8a7	526	}
	527
	528	delete [] realSM;
	529	delete [] idealSM;
	530	}
68711871	531