From 483fd7d82d3965600b95cdf514b2a608f44bb736 Mon Sep 17 00:00:00 2001 From: masera Date: Thu, 8 Nov 2007 09:23:56 +0000 Subject: [PATCH] New class for ITS alignment with Millepede (M. Lunardon) --- ITS/AliITSAlignMille.cxx | 825 +++++++++++++++++++++++++++++++++++++++ ITS/AliITSAlignMille.h | 155 ++++++++ ITS/ITSrecLinkDef.h | 3 + ITS/libITSrec.pkg | 3 +- 4 files changed, 985 insertions(+), 1 deletion(-) create mode 100644 ITS/AliITSAlignMille.cxx create mode 100644 ITS/AliITSAlignMille.h diff --git a/ITS/AliITSAlignMille.cxx b/ITS/AliITSAlignMille.cxx new file mode 100644 index 00000000000..9de58c7336e --- /dev/null +++ b/ITS/AliITSAlignMille.cxx @@ -0,0 +1,825 @@ +/************************************************************************** + * Copyright(c) 2007-2009, 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$ */ + +//----------------------------------------------------------------------------- +/// \class AliITSAlignMille +/// Alignment class fro the ALICE ITS detector +/// +/// ITS specific alignment class which interface to AliMillepede. +/// For each track ProcessTrack calculates the local and global derivatives +/// at each hit and fill the corresponding local equations. Provide methods for +/// fixing or constraining detection elements for best results. +/// +/// \author M. Lunardon (thanks to J. Castillo) +//----------------------------------------------------------------------------- + +#include +#include +#include +#include + +#include "AliITSAlignMille.h" +#include "AliITSgeomTGeo.h" +#include "AliGeomManager.h" +#include "AliMillepede.h" +#include "AliTrackPointArray.h" +#include "AliAlignObjParams.h" +#include "AliLog.h" +#include "TSystem.h" // come si fa? + +/// \cond CLASSIMP +ClassImp(AliITSAlignMille) +/// \endcond + +Int_t AliITSAlignMille::fgNDetElem = ITSMILLE_NDETELEM; +Int_t AliITSAlignMille::fgNParCh = ITSMILLE_NPARCH; + +AliITSAlignMille::AliITSAlignMille(const Char_t *configFilename, Bool_t initmille) + : TObject(), + fMillepede(0), + fStartFac(16.), + fResCutInitial(100.), + fResCut(100.), + fNGlobal(ITSMILLE_NDETELEM*ITSMILLE_NPARCH), + fNLocal(ITSMILLE_NLOCAL), + fNStdDev(ITSMILLE_NSTDEV), + fIsMilleInit(kFALSE), + fParSigTranslations(0.0100), + fParSigRotations(0.1), + fTrack(NULL), + fCluster(), + fGlobalDerivatives(NULL), + fTempHMat(NULL), + fTempAlignObj(NULL), + fDerivativeXLoc(0), + fDerivativeZLoc(0), + fDeltaPar(0), + fMinNPtsPerTrack(3), + fGeometryFileName("geometry.root"), + fGeoManager(0), + fCurrentModuleIndex(0), + fCurrentModuleInternalIndex(0), + fNModules(0), + fUseLocalShifts(kTRUE), + fCurrentModuleHMatrix(NULL) +{ + /// main constructor that takes input from configuration file + + fMillepede = new AliMillepede(); + fGlobalDerivatives = new Double_t[fNGlobal]; + fTempHMat = new TGeoHMatrix; + fCurrentModuleHMatrix = new TGeoHMatrix; + + Int_t lc=LoadConfig(configFilename); + if (lc) { + AliInfo(Form("Error %d loading configuration from %s",lc,configFilename)); + } + else { + if (initmille && fNGlobal<20000) { + AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, fNLocal, fNStdDev)); + Init(fNGlobal, fNLocal, fNStdDev); + ResetLocalEquation(); + AliInfo("Parameters initialized to zero"); + } + else { + AliInfo("Millepede has not been initialized ... "); + } + } + + fDeltaPar=0.0; // not used at the moment - to be checked later... + +} + +AliITSAlignMille::~AliITSAlignMille() { + /// Destructor + if (fMillepede) delete fMillepede; + delete [] fGlobalDerivatives; + delete fCurrentModuleHMatrix; + delete fTempHMat; +} + +/////////////////////////////////////////////////////////////////////// + +Int_t AliITSAlignMille::LoadConfig(const Char_t *cfile) { + /// return 0 if success + /// 1 if error in module index or voluid + + FILE *pfc=fopen(cfile,"r"); + if (!pfc) return -1; + + Char_t st[200],st2[200]; + Char_t tmp[100]; + Int_t idx,itx,ity,itz,ith,ips,iph; + Float_t f1; + UShort_t voluid; + Int_t nmod=0; + + while (fgets(st,200,pfc)) { + + // skip comments + for (int i=0; iAccessPathName(st2)) { + AliInfo("*** WARNING! *** geometry file not found! "); + return -1; + } + fGeometryFileName=st2; + InitGeometry(); + } + + if (strstr(st,"SET_PARSIG_TRA")) { + sscanf(st,"%s %f",tmp,&f1); + fParSigTranslations=f1; + } + + if (strstr(st,"SET_PARSIG_ROT")) { + sscanf(st,"%s %f",tmp,&f1); + fParSigRotations=f1; + } + + if (strstr(st,"SET_NSTDDEV")) { + sscanf(st,"%s %d",tmp,&idx); + fNStdDev=idx; + } + + if (strstr(st,"SET_RESCUT_INIT")) { + sscanf(st,"%s %f",tmp,&f1); + fResCutInitial=f1; + } + + if (strstr(st,"SET_RESCUT_OTHER")) { + sscanf(st,"%s %f",tmp,&f1); + fResCut=f1; + } + + if (strstr(st,"SET_STARTFAC")) { + sscanf(st,"%s %f",tmp,&f1); + fStartFac=f1; + } + + if (strstr(st,"SET_LOCAL_SHIFTS")) { // only enabled mode... + fUseLocalShifts = kTRUE; + } + + if (strstr(st,"MODULE_INDEX")) { + sscanf(st,"%s %d %d %d %d %d %d %d",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips); + voluid=GetModuleVolumeID(idx); + if (!voluid) return 1; // bad index + fModuleIndex[nmod]=idx; + fModuleVolumeID[nmod]=voluid; + fFreeParam[nmod][0]=itx; + fFreeParam[nmod][1]=ity; + fFreeParam[nmod][2]=itz; + fFreeParam[nmod][3]=iph; + fFreeParam[nmod][4]=ith; + fFreeParam[nmod][5]=ips; + nmod++; + } + + if (strstr(st,"MODULE_VOLUID")) { + // to be implemented + } + + } // end while + + fNModules = nmod; + fNGlobal = fNModules*fgNParCh; + + fclose(pfc); + return 0; +} + +Int_t AliITSAlignMille::GetModuleIndex(const Char_t *symname) { + /// index from symname + if (!symname) return -1; + for (Int_t i=0; i<2198; i++) { + if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i; + } + return -1; +} + +Int_t AliITSAlignMille::GetModuleIndex(UShort_t voluid) { + /// index from volume ID + AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid); + if (lay<1|| lay>6) return -1; + Int_t idx=Int_t(voluid)-2048*lay; + if (idx>=AliGeomManager::LayerSize(lay)) return -1; + for (Int_t ilay=1; ilay0 && layerId<7 && modId>=0 && modId2197) return 0; + return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index)); +} + +void AliITSAlignMille::InitGeometry() { + /// initialize geometry + AliGeomManager::LoadGeometry(fGeometryFileName.Data()); + fGeoManager = AliGeomManager::GetGeometry(); + if (!fGeoManager) { + AliInfo("Couldn't initialize geometry"); + return; + } + // temporary align object, just use the rotation... + fTempAlignObj=new AliAlignObjParams(AliITSgeomTGeo::GetSymName(7),2055,0,0,0,0,0,0,kFALSE); +} + +void AliITSAlignMille::Init(Int_t nGlobal, /* number of global paramers */ + Int_t nLocal, /* number of local parameters */ + Int_t nStdDev /* std dev cut */ ) +{ + /// Initialization of AliMillepede. Fix parameters, define constraints ... + fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial); + fIsMilleInit = kTRUE; + + /// Fix non free parameters + for (Int_t i=0; i1) fMillepede->SetIterations(fStartFac); +} + + +void AliITSAlignMille::AddConstraint(Double_t *par, Double_t value) { + /// Constrain equation defined by par to value + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return; + } + fMillepede->SetGlobalConstraint(par, value); + AliInfo("Adding constraint"); +} + +void AliITSAlignMille::InitGlobalParameters(Double_t *par) { + /// Initialize global parameters with par array + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return; + } + fMillepede->SetGlobalParameters(par); + AliInfo("Init Global Parameters"); +} + +void AliITSAlignMille::FixParameter(Int_t iPar, Double_t value) { + /// Parameter iPar is encourage to vary in [-value;value]. + /// If value == 0, parameter is fixed + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return; + } + fMillepede->SetParSigma(iPar, value); + if (value==0) AliInfo(Form("Parameter %i Fixed", iPar)); +} + +void AliITSAlignMille::ResetLocalEquation() +{ + /// Reset the derivative vectors + for(int i=0; i=0 && !(voluid==fModuleVolumeID[k]) ) k--; + if (k<0) return -3; + fCurrentModuleInternalIndex=k; + + // set the index + Int_t index = GetModuleIndex(AliGeomManager::SymName(voluid)); + if (index<0) return -2; + fCurrentModuleIndex = index; + + fModuleInitParam[0] = 0.0; + fModuleInitParam[1] = 0.0; + fModuleInitParam[2] = 0.0; + fModuleInitParam[3] = 0.0; // psi (X) + fModuleInitParam[4] = 0.0; // theta (Y) + fModuleInitParam[5] = 0.0; // phi (Z) + + /// get global (corrected) matrix + // if (!AliITSgeomTGeo::GetOrigMatrix(index,*fCurrentModuleHMatrix)) return -3; + Double_t rott[9]; + if (!AliITSgeomTGeo::GetRotation(index,rott)) return -3; + fCurrentModuleHMatrix->SetRotation(rott); + Double_t oLoc[3]={0,0,0}; + if (!AliITSgeomTGeo::LocalToGlobal(index,oLoc,fCurrentModuleTranslation)) return -4; + fCurrentModuleHMatrix->SetTranslation(fCurrentModuleTranslation); + + /// get back local coordinates + fMeasGlo[0] = fCluster.GetX(); + fMeasGlo[1] = fCluster.GetY(); + fMeasGlo[2] = fCluster.GetZ(); + fCurrentModuleHMatrix->MasterToLocal(fMeasGlo,fMeasLoc); + + // set stdev from cluster + TGeoHMatrix hcov; + Double_t hcovel[9]; + hcovel[0]=double(fCluster.GetCov()[0]); + hcovel[1]=double(fCluster.GetCov()[1]); + hcovel[2]=double(fCluster.GetCov()[3]); + hcovel[3]=double(fCluster.GetCov()[1]); + hcovel[4]=double(fCluster.GetCov()[2]); + hcovel[5]=double(fCluster.GetCov()[4]); + hcovel[6]=double(fCluster.GetCov()[3]); + hcovel[7]=double(fCluster.GetCov()[4]); + hcovel[8]=double(fCluster.GetCov()[5]); + hcov.SetRotation(hcovel); + // now rotate in local system + hcov.MultiplyLeft(&fCurrentModuleHMatrix->Inverse()); + hcov.Multiply(fCurrentModuleHMatrix); + + // per i ruotati c'e' delle sigmaY che compaiono... prob + // e' un problema di troncamento + fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0])); + fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4])); + fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8])); + + AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%f fSigmaLocY=%f fSigmaLocZ=%f \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] )); + + return 0; +} + +void AliITSAlignMille::SetCurrentModule(Int_t index) { + /// + UShort_t voluid=GetModuleVolumeID(index); + if (voluid) { + fCluster.SetVolumeID(voluid); + fCluster.SetXYZ(0,0,0); + InitModuleParams(); + } +} + +void AliITSAlignMille::Print() { + /// + printf("*** AliMillepede for ITS ***\n"); + printf(" number of defined modules: %d\n",fNModules); + if (fGeoManager) + printf(" geometry loaded from %s\n",fGeometryFileName.Data()); + else + printf(" geometry not loaded\n"); + if (fUseLocalShifts) + printf(" Alignment shifts will be computed in LOCAL RS\n"); + else + printf(" Alignment shifts will be computed in GLOBAL RS\n"); + printf(" Millepede configuration parameters:\n"); + printf(" init parsig for translations : %.4f\n",fParSigTranslations); + printf(" init parsig for rotations : %.4f\n",fParSigRotations); + printf(" init value for chi2 cut : %.4f\n",fStartFac); + printf(" first iteration cut value : %.4f\n",fResCutInitial); + printf(" other iterations cut value : %.4f\n",fResCut); + printf(" number of stddev for chi2 cut : %d\n",fNStdDev); + + printf("List of defined modules:\n"); + printf(" intidx\tindex\tvoluid\tsymname\n"); + for (int i=0; i2197) return; + UShort_t voluid=fModuleVolumeID[fCurrentModuleInternalIndex]; + printf("Current module: index=%d voluid=%d\n",fCurrentModuleIndex,voluid); + printf(" symname:%s\n",AliGeomManager::SymName(voluid)); + printf(" TGeoHMatrix: \n"); + fCurrentModuleHMatrix->Print(); +} + + +void AliITSAlignMille::InitTrackParams(int meth) { + /// initialize local parameters with different methods + /// for current track (fTrack) + + switch (meth) { + case 1: // simple linear interpolation + // get local starting parameters (to be substituted by ESD track parms) + // local parms (fLocalInitParam[]) are: + // [0] = global x coord. of straight line intersection at y=0 plane + // [1] = global z coord. of straight line intersection at y=0 plane + // [2] = px/py + // [3] = pz/py + + // test #1: linear fit in x(y) and z(y) + Int_t npts = fTrack->GetNPoints(); + + TF1 *f1=new TF1("f1","[0]+x*[1]",-50,50); + + TGraph *g=new TGraph(npts,fTrack->GetY(),fTrack->GetX()); + g->Fit(f1,"RNQ"); + fLocalInitParam[0] = f1->GetParameter(0); + fLocalInitParam[2] = f1->GetParameter(1); + AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1))); + delete g; g=NULL; + + g=new TGraph(npts,fTrack->GetY(),fTrack->GetZ()); + g->Fit(f1,"RNQ"); + fLocalInitParam[1] = f1->GetParameter(0); + fLocalInitParam[3] = f1->GetParameter(1); + AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3])); + delete g; + delete f1; + + break; + } + +} +Bool_t AliITSAlignMille::CheckVolumeID(UShort_t voluid) const +{ + /// + Int_t k=fNModules-1; + while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--; + //printf("selected element with voluid=%d : %d\n",voluid,k); + if (k>=0) return kTRUE; + return kFALSE; +} + +Int_t AliITSAlignMille::CheckCurrentTrack() { + /// checks if AliTrackPoints belongs to defined modules + /// return number of good poins + /// return 0 if not enough points + + Int_t npts = fTrack->GetNPoints(); + Int_t ngoodpts=0; + // debug points + for (int j=0; jGetVolumeID()[j]; + if (CheckVolumeID(voluid)) { + ngoodpts++; + } + } + // pepo da controllare... + if (ngoodptsGetNPoints(); + AliDebug(2,Form("\n*** Processing track with %d points ***\n",npts)); + fTrack = track; + + // checks if there are enough good points + if (!CheckCurrentTrack()) { + AliInfo("Track with not enough good points - will not be used..."); + return -1; + } + + // set local starting parameters (to be substituted by ESD track parms) + // local parms (fLocalInitParam[]) are: + // [0] = global x coord. of straight line intersection at y=0 plane + // [1] = global z coord. of straight line intersection at y=0 plane + // [2] = px/py + // [3] = pz/py + InitTrackParams(1); + + for (Int_t ipt=0; iptGetPoint(fCluster,ipt); + if (!CheckVolumeID(fCluster.GetVolumeID())) continue; + + // set geometry parameters for the the current module + AliDebug(2,Form("\n--- processing point %d --- \n",ipt)); + if (InitModuleParams()) continue; + + AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n", track->GetVolumeID()[ipt], fCurrentModuleIndex ,fCurrentModuleInternalIndex, AliGeomManager::SymName(track->GetVolumeID()[ipt]) )); + AliDebug(2,Form(" Point = ( %f , %f , %f ) \n",track->GetX()[ipt],track->GetY()[ipt],track->GetZ()[ipt])); + + if (SetLocalEquations()) return -1; + + } // end loop over points + + return 0; +} + +Int_t AliITSAlignMille::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar) { + /// calculate track intersection point in local coordinates + /// according with a given set of parameters (local(4) and global(6)) + /// and fill fPintLoc/Glo + /// local are: pgx0, pgz0, ugx0, ugz0 + /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.) + /// return 0 if success + + AliDebug(3,Form("lpar = %g %g %g %g \ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5])); + + // vector of initial straight line direction in glob. coord + // ATTENZIONE: forse va rinormalizzato tutto... + Double_t v0g[3]; + //Double_t + v0g[0]=lpar[2]; + v0g[1]=1.0; + v0g[2]=lpar[3]; + + // intercept in yg=0 plane in glob coord + Double_t p0g[3]; + p0g[0]=lpar[0]; + p0g[1]=0.0; + p0g[2]=lpar[1]; + + + // prepare the TGeoHMatrix + Double_t tr[3],ang[3]; + //Double_t rad2deg=180./TMath::Pi(); + if (fUseLocalShifts) { // just Delta matrix + tr[0]=gpar[0]; + tr[1]=gpar[1]; + tr[2]=gpar[2]; + ang[0]=gpar[3]; // psi (X) + ang[1]=gpar[4]; // theta (Y) + ang[2]=gpar[5]; // phi (Z) + } + else { // total matrix with shifted parameter + AliInfo("global shifts not implemented yet!"); + return -1; + } + + //printf("fTempRot = 0x%x - ang = %g %g %g \n",fTempRot,gpar[5]*rad2deg,gpar[3]*rad2deg,gpar[4]*rad2deg); + + fTempAlignObj->SetRotation(ang[0],ang[1],ang[2]); + AliDebug(3,Form("Delta angles: psi=%f theta=%f phi=%f",ang[0],ang[1],ang[2])); + TGeoHMatrix hm; + fTempAlignObj->GetMatrix(hm); + fTempHMat->SetRotation(hm.GetRotationMatrix()); + fTempHMat->SetTranslation(tr); + + // in this case the gpar[] array contains only shifts + // and fInitModuleParam[] are set to 0 + // fCurrentModuleHMatrix is then modified as fCurrentHM*fTempHM + if (fUseLocalShifts) + fTempHMat->MultiplyLeft(fCurrentModuleHMatrix); + + // same in local coord. + Double_t p0l[3],v0l[3]; + fTempHMat->MasterToLocalVect(v0g,v0l); + fTempHMat->MasterToLocal(p0g,p0l); + + if (TMath::Abs(v0l[1])<1e-15) { + AliInfo("Track Y direction in local frame is zero! Cannot proceed..."); + return -1; + } + + // local intersection point + fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1]; + fPintLoc[1] = 0; + fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1]; + + // global intersection point + fTempHMat->LocalToMaster(fPintLoc,fPintGlo); + AliDebug(3,Form("Intesect. point: L( %f , %f , %f ) G( %f , %f , %f )\n",fPintLoc[0],fPintLoc[1],fPintLoc[2],fPintGlo[0],fPintGlo[1],fPintGlo[2])); + + return 0; +} + +Int_t AliITSAlignMille::CalcDerivatives(Int_t paridx, Bool_t islpar) { + /// calculate numerically (ROOT's style) the derivatives for + /// local X intersection and local Z intersection + /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0 + /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg) + /// return 0 if success + + // copy initial parameters + Double_t lpar[ITSMILLE_NLOCAL]; + Double_t gpar[ITSMILLE_NPARCH]; + for (Int_t i=0; iSetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]); + + + AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2])); + // set equation for Zloc coordinate + for (Int_t i=0; iSetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]); + + return 0; +} + + +void AliITSAlignMille::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) { + /// Call local fit for this track + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return; + } + Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit); + AliDebug(2,Form("iRes = %d",iRes)); + if (iRes && !lSingleFit) { + fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1); + } +} + +void AliITSAlignMille::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) { + /// Call global fit; Global parameters are stored in parameters + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return; + } + fMillepede->GlobalFit(parameters,errors,pulls); + AliInfo("Done fitting global parameters!"); +} + +Double_t AliITSAlignMille::GetParError(Int_t iPar) { + /// Get error of parameter iPar + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return 0; + } + Double_t lErr = fMillepede->GetParError(iPar); + return lErr; +} + +void AliITSAlignMille::PrintGlobalParameters() { + /// Print global parameters + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return; + } + fMillepede->PrintGlobalParameters(); +} + +// //_________________________________________________________________________ + diff --git a/ITS/AliITSAlignMille.h b/ITS/AliITSAlignMille.h new file mode 100644 index 00000000000..202d998e643 --- /dev/null +++ b/ITS/AliITSAlignMille.h @@ -0,0 +1,155 @@ +#ifndef ALIITSALIGNMILLE_H +#define ALIITSALIGNMILLE_H +/* Copyright(c) 2007-2009 , ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + + +/// \ingroup rec +/// \class AliITSAlignMille +/// \brief Class for alignment of ITS +// +// Authors: Marcello Lunardon + +#include +#include +#include "AliTrackPointArray.h" + +class AliMillepede; +class AliAlignObjParams; +class TGeoManager; +class TGeoHMatrix; + +// number of used objects +#define ITSMILLE_NDETELEM 2198 +#define ITSMILLE_NPARCH 6 +#define ITSMILLE_NLOCAL 4 +#define ITSMILLE_NSTDEV 3 + +class AliITSAlignMille:public TObject +{ + +public: + AliITSAlignMille(const Char_t *configFilename="AliITSAlignMille.conf", Bool_t initmille=kTRUE); + virtual ~AliITSAlignMille(); + + // geometry methods + Int_t GetModuleIndex(const Char_t *symname); + Int_t GetModuleIndex(UShort_t voluid); + UShort_t GetModuleVolumeID(const Char_t *symname); + UShort_t GetModuleVolumeID(Int_t index); + void SetCurrentModule(Int_t index); + + // configuration methods + void SetGeometryFileName(const Char_t* filename="geometry.root") + { fGeometryFileName = filename; } + const Char_t* GetGeometryFileName() {return fGeometryFileName.Data();} + void PrintCurrentModuleInfo(); + void Print(); + + // fitting methods + void SetMinNPtsPerTrack(Int_t pts=3) {fMinNPtsPerTrack=pts;} + //Bool_t CheckTrack(AliTrackPointArray *track); + Int_t ProcessTrack(AliTrackPointArray *track); + void InitTrackParams(int meth=1); + Int_t InitModuleParams(); + Int_t CheckCurrentTrack(); + Bool_t CheckVolumeID(UShort_t voluid) const ; + Int_t CalcIntersectionPoint(Double_t *lpar, Double_t *gpar); + Int_t CalcDerivatives(Int_t paridx, Bool_t islpar); + Double_t* GetLocalIntersectionPoint() {return fPintLoc;} + Double_t* GetGlobalIntersectionPoint() {return fPintGlo;} + + // millepede methods + void FixParameter(Int_t param, Double_t value); + void AddConstraint(Double_t *factor, Double_t value ); + void InitGlobalParameters(Double_t *par); + void SetLocalDerivative(Int_t index, Double_t value) + {fLocalDerivatives[index] = value;} + void SetGlobalDerivative(Int_t index, Double_t value) + {fGlobalDerivatives[index] = value;} + void LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit); + void GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls); + void PrintGlobalParameters(); + Double_t GetParError(Int_t iPar); + Int_t SetLocalEquations(); + + // fitting stuffs + AliTrackPointArray *GetCurrentTrack() {return fTrack;} + AliTrackPoint *GetCurrentCluster() {return &fCluster;} + + // geometry stuffs + Int_t GetNModules() const {return fNModules;} + Int_t GetCurrentModuleIndex() const {return fCurrentModuleIndex;} + TGeoHMatrix *GetCurrentModuleHMatrix() {return fCurrentModuleHMatrix;} + Double_t *GetCurrentModuleTranslation() {return fCurrentModuleTranslation;} + Int_t GetCurrentModuleInternalIndex() const {return fCurrentModuleInternalIndex;} + Int_t *GetModuleIndexArray() {return fModuleIndex;} + UShort_t *GetModuleVolumeIDArray() {return fModuleVolumeID;} + + private: + + // configuration methods + Int_t LoadConfig(const Char_t *cfile="AliITSAlignMille.conf"); + void ResetLocalEquation(); + void InitGeometry(); + + // millepede methods + void Init(Int_t nGlobal, Int_t nLocal, Int_t nStdDev); + + // millepede stuffs + AliMillepede *fMillepede; ///< Detector independent alignment class + static Int_t fgNParCh; ///< Number of degrees of freedom per chamber + static Int_t fgNDetElem; ///< Total number of detection elements + Double_t fStartFac; ///< Initial value for chi2 cut + ///< if > 1 Iterations in AliMil. are turned on + Double_t fResCutInitial; ///< Cut on residual for first iteration + Double_t fResCut; ///< Cut on residual for other iterations + Int_t fNGlobal; ///< Number of global parameters + Int_t fNLocal; ///< Number of local parameters + Int_t fNStdDev; ///< Number of standard deviations for chi2 cut + Bool_t fIsMilleInit; /// + Double_t fParSigTranslations; ///< init sigma for transl. params [cm] + Double_t fParSigRotations; ///< init sigma for rot. params [deg] + + // fitting stuffs + AliTrackPointArray *fTrack; ///< pointer to current track + AliTrackPoint fCluster; ///< current cluster + Double_t *fGlobalDerivatives; ///< Array of global derivatives + Double_t fLocalDerivatives[ITSMILLE_NLOCAL]; ///< Array of local deriv. + Double_t fLocalInitParam[ITSMILLE_NLOCAL]; ///< Array with inital values for local parameters for current track + Double_t fModuleInitParam[ITSMILLE_NPARCH]; ///< Array with inital values for current module parameters (init geometry) + Double_t fPintLoc[3]; /// + Double_t fPintLoc0[3]; /// + Double_t fPintGlo[3]; /// + Double_t fMeasLoc[3]; // current point local coordinates (the original ones) + Double_t fMeasGlo[3]; // current point glob. coord (AliTrackPoint) + Double_t fSigmaLoc[3]; // stdev current point + TGeoHMatrix *fTempHMat; /// + AliAlignObjParams *fTempAlignObj; /// + Double_t fDerivativeXLoc; // localX deriv. + Double_t fDerivativeZLoc; // localZ deriv. + Double_t fDeltaPar; /// + Int_t fMinNPtsPerTrack; /// + + // geometry stuffs + TString fGeometryFileName; /// + TGeoManager *fGeoManager; /// + Int_t fCurrentModuleIndex; /// + Int_t fCurrentModuleInternalIndex; /// + Double_t fCurrentModuleTranslation[3]; /// + Int_t fNModules; /// number of defined modules from config file + Int_t fModuleIndex[ITSMILLE_NDETELEM]; /// + UShort_t fModuleVolumeID[ITSMILLE_NDETELEM]; /// + Bool_t fFreeParam[ITSMILLE_NDETELEM][ITSMILLE_NPARCH]; /// + Bool_t fUseLocalShifts; /// + TGeoHMatrix *fCurrentModuleHMatrix; /// + + AliITSAlignMille(const AliITSAlignMille& rhs); + AliITSAlignMille& operator=(const AliITSAlignMille& rhs); + + +ClassDef(AliITSAlignMille, 0)}; + +#endif diff --git a/ITS/ITSrecLinkDef.h b/ITS/ITSrecLinkDef.h index 34d567fdb46..d7ae6f9a0cf 100644 --- a/ITS/ITSrecLinkDef.h +++ b/ITS/ITSrecLinkDef.h @@ -115,4 +115,7 @@ #pragma link C++ class AliITSOnlineSDDInjectors+; #pragma link C++ class AliITSPreprocessorSSD+; +// Classes for alignment +#pragma link C++ class AliITSAlignMille+; + #endif diff --git a/ITS/libITSrec.pkg b/ITS/libITSrec.pkg index c37a6521749..029eb0eec11 100644 --- a/ITS/libITSrec.pkg +++ b/ITS/libITSrec.pkg @@ -80,7 +80,8 @@ SRCS = AliITSDetTypeRec.cxx \ AliITSOnlineSDDBTP.cxx \ AliITSOnlineSDDInjectors.cxx \ AliITSOnlineSDDCMN.cxx \ - AliITSPreprocessorSSD.cxx + AliITSPreprocessorSSD.cxx \ + AliITSAlignMille.cxx \ HDRS:= $(SRCS:.cxx=.h) -- 2.39.3