From: shahoian Date: Thu, 22 Jan 2009 14:37:47 +0000 (+0000) Subject: Added ITS Alignment interface to MillePede2 and related supermodule class X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=7b85e477508179cc269d3c3a43d83cdd5e99cc5a;p=u%2Fmrichter%2FAliRoot.git Added ITS Alignment interface to MillePede2 and related supermodule class --- diff --git a/ITS/AliITSAlignMille2.cxx b/ITS/AliITSAlignMille2.cxx new file mode 100644 index 00000000000..be6802770fd --- /dev/null +++ b/ITS/AliITSAlignMille2.cxx @@ -0,0 +1,1842 @@ +/************************************************************************** + * 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: AliITSAlignMille.cxx 25025 2008-04-09 20:50:08Z masera $ */ + +//----------------------------------------------------------------------------- +/// \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 +#include +#include +#include + +#include "AliITSAlignMille2.h" +#include "AliITSgeomTGeo.h" +#include "AliGeomManager.h" +#include "AliMillePede2.h" +#include "AliTrackPointArray.h" +#include "AliAlignObjParams.h" +#include "AliLog.h" +#include "TSystem.h" // come si fa? +#include "AliTrackFitterRieman.h" + +/// \cond CLASSIMP +ClassImp(AliITSAlignMille2) +/// \endcond + +AliITSAlignMille2* AliITSAlignMille2::fgInstance = 0; +Int_t AliITSAlignMille2::fgInstanceID = 0; + +AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename, Bool_t initmille) +: TObject(), + fMillepede(0), + fStartFac(16.), + fResCutInitial(100.), + fResCut(100.), + fNGlobal(0), + fNLocal(4), + fNStdDev(3), + fIsMilleInit(kFALSE), + fSensorsIn(kFALSE), + fParSigTranslations(0.0100), + fParSigRotations(0.1), +// + fCurrentModule(0), + fTrack(0), + fCluster(), + fGlobalDerivatives(0), +// + fMinNPtsPerTrack(3), + fInitTrackParamsMeth(1), + fTotBadLocEqPoints(0), + fRieman(0), + // + fUseGlobalDelta(kFALSE), + fRequirePoints(kFALSE), + fTempExcludedModule(-1), + // + fGeometryFileName("geometry.root"), + fPreAlignmentFileName(""), + fGeoManager(0), + fIsConfigured(kFALSE), + fPreAlignQF(0), +// + fPrealignment(0), + fMilleModule(2), + fSuperModule(2), + fNModules(0), + fNSuperModules(0), + fUsePreAlignment(kFALSE), + fUseSortedTracks(kTRUE), + fBOn(kFALSE), + fBField(0.0), + fBug(0) +{ + /// main constructor that takes input from configuration file + // + fMillepede = new AliMillePede2(); + for (int i=3;i--;) fSigmaFactor[i] = 1.0; + // + // new RS + for (Int_t i=0; i<6; i++) { + fNReqLayUp[i]=0; + fNReqLayDown[i]=0; + fNReqLay[i]=0; + } + for (Int_t i=0; i<3; i++) { + fNReqDetUp[i]=0; + fNReqDetDown[i]=0; + fNReqDet[i]=0; + } + // + Int_t lc=LoadConfig(configFilename); + if (lc) { + AliInfo(Form("Error %d loading configuration from %s",lc,configFilename)); + } + else { + fIsConfigured=kTRUE; + if (initmille) { + 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 ... "); + } + } + // + fgInstance = this; + fgInstanceID++; + // +} + +AliITSAlignMille2::~AliITSAlignMille2() { + /// Destructor + if (fMillepede) delete fMillepede; + if (fGlobalDerivatives) delete[] fGlobalDerivatives; + if (fRieman) delete fRieman; + if (fPrealignment) delete fPrealignment; + fMilleModule.Delete(); + fSuperModule.Delete(); + if (--fgInstanceID==0) fgInstance = 0; +} + +/////////////////////////////////////////////////////////////////////// + +Int_t AliITSAlignMille2::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,f2,f3; + UShort_t voluid; + Int_t nmod=0; + // + while (fgets(st,200,pfc)) { + // + for (int i=0; iAccessPathName(st2)) { AliInfo("*** WARNING! *** geometry file not found! "); return -1;} + fGeometryFileName=st2; + InitGeometry(); + } + // + if (strstr(st,"PREALIGNMENT_FILE")) { + sscanf(st,"%s %s",tmp,st2); + if (gSystem->AccessPathName(st2)) { AliInfo("*** WARNING! *** prealignment file not found! "); return -1;} + fPreAlignmentFileName=st2; + itx=ApplyToGeometry(); + if (itx) { AliInfo(Form("*** WARNING! *** error %d reading prealignment file! ",itx)); return -6;} + } + // + if (strstr(st,"SUPERMODULE_FILE")) { + sscanf(st,"%s %s",tmp,st2); + if (gSystem->AccessPathName(st2)) { AliInfo("*** WARNING! *** supermodule file not found! "); return -1;} + if (LoadSuperModuleFile(st2)) return -1; + } + // + if (strstr(st,"SET_B_FIELD")) { + sscanf(st,"%s %f",tmp,&f1); + if (f1>0) { + fBField = f1; + fBOn = kTRUE; + fNLocal = 5; // helices + fRieman = new AliTrackFitterRieman(); + } + else { + fBField = 0.0; + fBOn = kFALSE; + fNLocal = 4; + } + } + // + if (strstr(st,"SET_MINPNT_TRA")) { + sscanf(st,"%s %d",tmp,&idx); + fMinNPtsPerTrack=idx; + } + // + 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_LOCALSIGMAFACTOR")) { + f1=f2=f3=0; + sscanf(st,"%s %f %f %f",tmp,&f1,&f2,&f3); + if (f1>0) fSigmaFactor[0] = f1; + if (f2>0) fSigmaFactor[1] = f2; else fSigmaFactor[1]=fSigmaFactor[0]; + if (f3>0) fSigmaFactor[2] = f3; else fSigmaFactor[2]=fSigmaFactor[1]; + } + // + if (strstr(st,"SET_STARTFAC")) { + sscanf(st,"%s %f",tmp,&f1); + fStartFac=f1; + } + // + // >> RS + if (strstr(st,"REQUIRE_POINT")) { + // syntax: REQUIRE_POINT where ndet updw nreqpts + // where = LAYER or DETECTOR + // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3) + // updw = 1 for Y>0, -1 for Y<0, 0 if not specified + // nreqpts = minimum number of points of that type + sscanf(st,"%s %s %d %d %d",tmp,st2,&itx,&ity,&itz); + itx--; + if (strstr(st2,"LAYER")) { + if (itx<0 || itx>5) return -7; + if (ity>0) fNReqLayUp[itx]=itz; + else if (ity<0) fNReqLayDown[itx]=itz; + else fNReqLay[itx]=itz; + fRequirePoints=kTRUE; + } + else if (strstr(st2,"DETECTOR")) { // DETECTOR + if (itx<0 || itx>2) return -7; + if (ity>0) fNReqDetUp[itx]=itz; + else if (ity<0) fNReqDetDown[itx]=itz; + else fNReqDet[itx]=itz; + fRequirePoints=kTRUE; + } + } + // << RS + + // + if (strstr(st,"MODULE_INDEX") || strstr(st,"MODULE_VOLUID")) { + f1=f2=f3=0; + sscanf(st,"%s %d %d %d %d %d %d %d %f %f %f",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips,&f1,&f2,&f3); + // + if (idx<=kMaxITSSensID) voluid=GetModuleVolumeID(idx); + else voluid = UShort_t(idx); + // + if (voluid>=kMinITSSupeModuleID) { // custom supermodule + int ism=-1; + for (int j=0; jGetVolumeID()) ism=j; + if (ism<0) return -1; // bad volid + fMilleModule.AddAtAndExpand(new AliITSAlignMille2Module(*GetSuperModule(ism)),nmod); + // >> RS +// if (f1>0) { +// for (int kk=0; kkGetNSensitiveVolumes(); kk++) { +// idx=AliITSAlignMille2Module::GetIndexFromVolumeID(GetMilleModule(nmod)->GetSensitiveVolumeVolumeID()[kk]); +// if (idx>=0) fSensVolSigmaXfactor[idx]=f1; +// } +// } +// if (f2>0) { +// for (int kk=0; kkGetNSensitiveVolumes(); kk++) { +// idx=AliITSAlignMille2Module::GetIndexFromVolumeID(GetMilleModule(nmod)->GetSensitiveVolumeVolumeID()[kk]); +// if (idx>=0) fSensVolSigmaZfactor[idx]=f2; +// } +// } + // << RS + } + else if (idx<=kMaxITSSensVID) { + fMilleModule.AddAtAndExpand(new AliITSAlignMille2Module(voluid),nmod); + AliITSAlignMille2Module* md = (AliITSAlignMille2Module*) fMilleModule[nmod]; + fSensorsIn = kTRUE; + md->SetSensorsProvided(); + } + else return -1; // bad volid + // + AliITSAlignMille2Module* mod = GetMilleModule(nmod); + mod->SetFreeDOF(kDOFTX,itx!=0); + mod->SetFreeDOF(kDOFTY,ity!=0); + mod->SetFreeDOF(kDOFTZ,itz!=0); + mod->SetFreeDOF(kDOFPH,iph!=0); + mod->SetFreeDOF(kDOFTH,ith!=0); + mod->SetFreeDOF(kDOFPS,ips!=0); + // + mod->SetUniqueID(nmod); + if (f1>0) mod->SetSigmaXFactor(f1); + if (f2>0) mod->SetSigmaYFactor(f2); else mod->SetSigmaYFactor(mod->GetSigmaXFactor()); + if (f3>0) mod->SetSigmaZFactor(f3); else mod->SetSigmaZFactor(mod->GetSigmaYFactor()); + nmod++; + } + // + } // end while + // + fNModules = nmod; + fNGlobal = fNModules*kNParCh; + // + fclose(pfc); + // + // set parent/child relationship for modules to align + printf("Setting parent/child relationships\n"); + // + for (int ipar=0;iparGetIndex()<=kMaxITSSensID) continue; // sensor cannot be a parent + // + for (int icld=0;icldBelongsTo(parent)) continue; + // + AliITSAlignMille2Module* parOld = child->GetParent(); + if (parOld && parOld->GetNSensitiveVolumes()GetNSensitiveVolumes()) continue; // parOld is closer + child->SetParent(parent); + } + // + } + // + // reorder the modules in such a way that parents come first + for (int icld=0;icldGetParent()) continue; + // + for (int icld=0;icldGetParent(); + if (!parent || parent->GetUniqueID()GetUniqueID()) continue; + // swap + fMilleModule[icld] = parent; + fMilleModule[parent->GetUniqueID()] = child; + child->SetUniqueID(parent->GetUniqueID()); + parent->SetUniqueID(icld); + } + } + // + // go over the child->parent chain and mark modules with explicitly provided sensors + for (int icld=nmod;icld--;) { + AliITSAlignMille2Module* child = GetMilleModule(icld); + AliITSAlignMille2Module* parent = child->GetParent(); + if (!parent) continue; + parent->SetSensorsProvided( child->AreSensorsProvided() ); + if (!parent->AreSensorsProvided()) continue; + // suppress unused sensors + for (int isn=0;isnGetNSensitiveVolumes();isn++) { + int snVID = parent->GetSensVolVolumeID(isn); + // check if this sensor is explicitly requested + for (int imd=nmod;imd--;) if (GetMilleModule(imd)->GetVolumeID() == snVID) {snVID = -1; break;} + // + if (snVID==-1) continue; // found this sensor, do nothing + // + // otherwise, remove this sensor from the module list + AliInfo(Form("Removing sensor %d from %s",snVID,parent->GetName())); + parent->DelSensitiveVolume(isn--); + } + } + // + fGlobalDerivatives = new Double_t[fNGlobal]; + memset(fGlobalDerivatives,0,fNGlobal*sizeof(Double_t)); + // + return 0; +} + + +void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts) +{ + // set minimum number of points in specific detector or layer + // where = LAYER or DETECTOR + // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3) + // updw = 1 for Y>0, -1 for Y<0, 0 if not specified + // nreqpts = minimum number of points of that type + ndet--; + if (strstr(where,"LAYER")) { + if (ndet<0 || ndet>5) return; + if (updw>0) fNReqLayUp[ndet]=nreqpts; + else if (updw<0) fNReqLayDown[ndet]=nreqpts; + else fNReqLay[ndet]=nreqpts; + fRequirePoints=kTRUE; + } + else if (strstr(where,"DETECTOR")) { + if (ndet<0 || ndet>2) return; + if (updw>0) fNReqDetUp[ndet]=nreqpts; + else if (updw<0) fNReqDetDown[ndet]=nreqpts; + else fNReqDet[ndet]=nreqpts; + fRequirePoints=kTRUE; + } +} + + +Int_t AliITSAlignMille2::GetModuleIndex(const Char_t *symname) { + /// index from symname + if (!symname) return -1; + for (Int_t i=0;i<=kMaxITSSensID; i++) { + if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i; + } + return -1; +} + +Int_t AliITSAlignMille2::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 && modIdGetIndex()==index) return GetSuperModule(i)->GetVolumeID(); + } + } + return 0; +} + +void AliITSAlignMille2::InitGeometry() { + /// initialize geometry + AliGeomManager::LoadGeometry(fGeometryFileName.Data()); + fGeoManager = AliGeomManager::GetGeometry(); + if (!fGeoManager) { + AliInfo("Couldn't initialize geometry"); + return; + } +} + +void AliITSAlignMille2::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; iIsFreeDOF(j)) FixParameter(i*kNParCh+j,0.0); + else { + // pepopepo: da verificare il settaggio delle sigma, ma forse va bene... + Double_t parsig=0; + if (j<3) parsig = fParSigTranslations; // translations (0.0100 cm) + else parsig = fParSigRotations; // rotations (1/10 deg) + FixParameter(i*kNParCh+j,parsig); + } + } + } + // + // Set iterations + if (fStartFac>1) fMillepede->SetIterations(fStartFac); +} + + +void AliITSAlignMille2::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 AliITSAlignMille2::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 AliITSAlignMille2::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 AliITSAlignMille2::ResetLocalEquation() +{ + /// Reset the derivative vectors + for(int i=fNLocal;i--;) fLocalDerivatives[i] = 0.0; + for(int i=fNGlobal;i--;) fGlobalDerivatives[i] = 0.0; +} + +Int_t AliITSAlignMille2::ApplyToGeometry() +{ + // apply starting realignment to ideal geometry + // + if (!fGeoManager) return -1; + TFile *pref = new TFile(fPreAlignmentFileName.Data()); + if (!pref->IsOpen()) return -2; + fPrealignment = (TClonesArray*)pref->Get("ITSAlignObjs"); + if (!fPrealignment) return -3; + Int_t nprea = fPrealignment->GetEntriesFast(); + AliInfo(Form("Array of input misalignments with %d entries",nprea)); + // + for (int ix=0; ixAt(ix); + Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(preo->GetVolUID()); + if (index>=0) { + if (index>=fPreAlignQF.GetSize()) fPreAlignQF.Set(index+10); + fPreAlignQF[index] = (int) preo->GetUniqueID()+1; + } + //TString nms = preo->GetSymName(); + //if (!nms.Contains("Ladder")) continue; //RRR + if (!preo->ApplyToGeometry()) return -4; + } + // + pref->Close(); + delete pref; + // + fUsePreAlignment = kTRUE; + return 0; +} + +Int_t AliITSAlignMille2::GetPreAlignmentQualityFactor(Int_t index) const +{ + if (!fUsePreAlignment || index<0 || index>=fPreAlignQF.GetSize()) return -1; + return fPreAlignQF[index]-1; +} + +AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *atp) { + /// create a new AliTrackPointArray keeping only defined modules + /// move points according to a given prealignment, if any + /// sort alitrackpoints w.r.t. global Y direction, if selected + + AliTrackPointArray *atps=NULL; + Int_t idx[20]; + Int_t npts=atp->GetNPoints(); + + /// checks if AliTrackPoints belong to defined modules + Int_t ngoodpts=0; + Int_t intidx[20]; + + for (int j=0; jGetVolumeID()[j]); + if (intidx[j]>=0) ngoodpts++; + } + AliDebug(3,Form("Number of points in defined modules: %d",ngoodpts)); + + // reject track if not enough points are left + if (ngoodpts> RS + AliTrackPoint p; + // check points in specific places + if (fRequirePoints) { + Int_t nlayup[6],nlaydown[6],nlay[6]; + Int_t ndetup[3],ndetdown[3],ndet[3]; + for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;} + for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;} + + for (int i=0; iGetX()[i]; + Float_t yy=atp->GetY()[i]; + Float_t r=TMath::Sqrt(xx*xx + yy*yy); + int lay=-1; + if (r<5) lay=0; + else if (r>5 && r<10) lay=1; + else if (r>10 && r<18) lay=2; + else if (r>18 && r<30) lay=3; + else if (r>30 && r<40) lay=4; + else if (r>40) lay=5; + if (lay<0) continue; + int det=lay/2; + //printf("Point %d - x=%f y=%f R=%f lay=%d det=%d\n",i,xx,yy,r,lay,det); + + if (yy>=0.0) { // UP point + nlayup[lay]++; + nlay[lay]++; + ndetup[det]++; + ndet[det]++; + } + else { + nlaydown[lay]++; + nlay[lay]++; + ndetdown[det]++; + ndet[det]++; + } + } + + // checks minimum values + Bool_t isok=kTRUE; + for (Int_t j=0; j<6; j++) { + if (nlayup[j]GetY(),idx); // sort descending... + + Int_t npto=0; + for (int i=0; iGetPoint(p,idx[i]); + + // prealign point if required + // get IDEAL matrix + TGeoHMatrix *svOrigMatrix = GetMilleModule(intidx[idx[i]])->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID()); + // get back real local coordinates: use OriginalGlobalMatrix because AliTrackPoints were written + // with idel geometry + Double_t pg[3],pl[3]; + pg[0]=p.GetX(); + pg[1]=p.GetY(); + pg[2]=p.GetZ(); + // printf("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]); + AliDebug(3,Form("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2])); + svOrigMatrix->MasterToLocal(pg,pl); + + AliDebug(3,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2])); + + // update covariance matrix + TGeoHMatrix hcov; + Double_t hcovel[9]; + hcovel[0]=double(p.GetCov()[0]); + hcovel[1]=double(p.GetCov()[1]); + hcovel[2]=double(p.GetCov()[2]); + hcovel[3]=double(p.GetCov()[1]); + hcovel[4]=double(p.GetCov()[3]); + hcovel[5]=double(p.GetCov()[4]); + hcovel[6]=double(p.GetCov()[2]); + hcovel[7]=double(p.GetCov()[4]); + hcovel[8]=double(p.GetCov()[5]); + hcov.SetRotation(hcovel); + // now rotate in local system + hcov.Multiply(svOrigMatrix); + hcov.MultiplyLeft(&svOrigMatrix->Inverse()); + // now hcov is LOCAL COVARIANCE MATRIX + + // >> RS + if (fBug==1) { + // correzione bug LAYER 5 SSD temporanea.. + int ssdidx=AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID()); + if (ssdidx>=500 && ssdidx<1248) { + int ladder=(ssdidx-500)%22; + if (ladder==18) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx+1)); + if (ladder==19) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx-1)); + } + } + + // << RS + + /// get (evenctually prealigned) matrix of sens. vol. + TGeoHMatrix *svMatrix = GetMilleModule(intidx[idx[i]])->GetSensitiveVolumeMatrix(p.GetVolumeID()); + // modify global coordinates according with pre-aligment + svMatrix->LocalToMaster(pl,pg); + // now rotate in local system + hcov.Multiply(&svMatrix->Inverse()); + hcov.MultiplyLeft(svMatrix); + // hcov is back in GLOBAL RF + Float_t pcov[6]; + pcov[0]=hcov.GetRotationMatrix()[0]; + pcov[1]=hcov.GetRotationMatrix()[1]; + pcov[2]=hcov.GetRotationMatrix()[2]; + pcov[3]=hcov.GetRotationMatrix()[4]; + pcov[4]=hcov.GetRotationMatrix()[5]; + pcov[5]=hcov.GetRotationMatrix()[8]; + + p.SetXYZ(pg[0],pg[1],pg[2],pcov); + // printf("New Gl coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]); + AliDebug(3,Form("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2])); + atps->AddPoint(npto,&p); + AliDebug(2,Form("Adding point[%d] = ( %f , %f , %f ) volid = %d",npto,atps->GetX()[npto],atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] )); + + npto++; + } + + return atps; +} + + + +AliTrackPointArray *AliITSAlignMille2::SortTrack(const AliTrackPointArray *atp) { + /// sort alitrackpoints w.r.t. global Y direction + AliTrackPointArray *atps=NULL; + Int_t idx[20]; + Int_t npts=atp->GetNPoints(); + AliTrackPoint p; + atps=new AliTrackPointArray(npts); + + TMath::Sort(npts,atp->GetY(),idx); + + for (int i=0; iGetPoint(p,idx[i]); + atps->AddPoint(i,&p); + AliDebug(2,Form("Point[%d] = ( %f , %f , %f ) volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] )); + } + return atps; +} + + +Int_t AliITSAlignMille2::InitModuleParams() { + /// initialize geometry parameters for a given detector + /// for current cluster (fCluster) + /// fGlobalInitParam[] is set as: + /// [tx,ty,tz,psi,theta,phi] + /// (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...) + /// *** At the moment: using Raffalele's angles definition *** + /// + /// return 0 if success + + if (!fGeoManager) { + AliInfo("ITS geometry not initialized!"); + return -1; + } + + // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster) + + // set the internal index (index in module list) + UShort_t voluid=fCluster.GetVolumeID(); + // + // IT IS VERY IMPORTANT: start from the end of the list, where the childs are located !!! + Int_t k=fNModules-1; + fCurrentModule = 0; + while (k>=0 && ! (fCurrentModule=GetMilleModule(k))->IsIn(voluid)) { + // VERY IMPORTANT: if the sensors were explicitly provided, don't look in the supermodules + if (fSensorsIn && fCurrentModule->GetVolumeID() > kMaxITSSensVID) {k=-1; break;} + k--; + } + if (k<0) return -3; + // fCurrentModule = GetMilleModule(k); + // + 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 (evenctually prealigned) matrix of sens. vol. + TGeoHMatrix *svMatrix = fCurrentModule->GetSensitiveVolumeMatrix(voluid); + + fMeasGlo[0] = fCluster.GetX(); + fMeasGlo[1] = fCluster.GetY(); + fMeasGlo[2] = fCluster.GetZ(); + svMatrix->MasterToLocal(fMeasGlo,fMeasLoc); + AliDebug(2,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] )); + + // 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()[2]); + hcovel[3]=double(fCluster.GetCov()[1]); + hcovel[4]=double(fCluster.GetCov()[3]); + hcovel[5]=double(fCluster.GetCov()[4]); + hcovel[6]=double(fCluster.GetCov()[2]); + hcovel[7]=double(fCluster.GetCov()[4]); + hcovel[8]=double(fCluster.GetCov()[5]); + hcov.SetRotation(hcovel); + // now rotate in local system + hcov.Multiply(svMatrix); + hcov.MultiplyLeft(&svMatrix->Inverse()); + + // set local sigmas + fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0])); + fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4])); // RS + fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8])); + + // set minimum value for SigmaLoc to 10 micron + if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010; + if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010; + + // multiply local sigmas by global and module specific factor + for (int i=3;i--;) fSigmaLoc[i] *= fSigmaFactor[i]*fCurrentModule->GetSigmaFactor(i); + // + AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] )); + + return 0; +} + +void AliITSAlignMille2::Print(Option_t*) const +{ + /// + printf("*** AliMillepede for ITS ***\n"); + printf(" Number of defined super modules: %d\n",fNModules); + printf(" Obtained parameters refer to %s Deltas\n",fUseGlobalDelta ? "GLOBAL":"LOCAL"); + // + if (fGeoManager) + printf(" geometry loaded from %s\n",fGeometryFileName.Data()); + else + printf(" geometry not loaded\n"); + // + if (fUsePreAlignment) + printf(" using prealignment from %s \n",fPreAlignmentFileName.Data()); + else + printf(" prealignment not used\n"); + // + // + if (fBOn) + printf(" B Field set to %f T - using Riemann's helices\n",fBField); + else + printf(" B Field OFF - using straight lines \n"); + // + if (fRequirePoints) printf(" Required points in tracks:\n"); + for (Int_t i=0; i<6; i++) { + if (fNReqLayUp[i]>0) printf(" Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[i]); + if (fNReqLayDown[i]>0) printf(" Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[i]); + if (fNReqLay[i]>0) printf(" Layer %d : %d points \n",i+1,fNReqLay[i]); + } + for (Int_t i=0; i<3; i++) { + if (fNReqDetUp[i]>0) printf(" Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[i]); + if (fNReqDetDown[i]>0) printf(" Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[i]); + if (fNReqDet[i]>0) printf(" Detector %d : %d points \n",i+1,fNReqDet[i]); + } + // + printf("\n 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(" mult. fact. for local sigmas : %.4f %.4f %.4f\n",fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2]); + + printf("List of defined modules:\n"); + printf(" intidx\tindex\tvoluid\tname\n"); + for (int i=0; iGetIndex(),md->GetVolumeID(),md->GetName()); + } +} + +AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleByVID(UShort_t voluid) const +{ + // return pointer to a define supermodule + // return NULL if error + Int_t i=IsDefined(voluid); + if (i<0) return NULL; + return GetMilleModule(i); +} + +Bool_t AliITSAlignMille2::InitRiemanFit() +{ + // Initialize Riemann Fitter for current track + // return kFALSE if error + + if (!fBOn) return kFALSE; + + Int_t npts=0; + AliTrackPoint ap; + npts = fTrack->GetNPoints(); + AliDebug(3,Form("Fitting track with %d points",npts)); + + fRieman->Reset(); + fRieman->SetTrackPointArray(fTrack); + + TArrayI ai(npts); + for (Int_t ipt=0; iptGetVolumeID()[ipt]; + + // fit track with 5 params in his own tracking-rotated reference system + // xc = -p[1]/p[0]; + // yc = 1/p[0]; + // R = sqrt( x0*x0 + y0*y0 - y0*p[2]); + if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) { + return kFALSE; + } + + for (int i=0; i<5; i++) + fLocalInitParam[i] = fRieman->GetParam()[i]; + + return kTRUE; +} + +Bool_t fullErr2D = kTRUE; + +void trackFit2D(Int_t &, Double_t *, double &chi2, double *par, int) +{ + const double kTiny = 1.e-14; + chi2 = 0; + static AliTrackPoint pnt; + // + enum {kAX,kAZ,kBX,kBZ}; + enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5}; + // + AliITSAlignMille2* alig = AliITSAlignMille2::GetInstance(); + AliTrackPointArray* track = alig->GetCurrentTrack(); + // + int npts = track->GetNPoints(); + for (int ip=0;ipGetPoint(pnt,ip); + const float *cov = pnt.GetCov(); + double y = pnt.GetY(); + double dx = pnt.GetX() - (par[kAX]+y*par[kBX]); + double dz = pnt.GetZ() - (par[kAZ]+y*par[kBZ]); + double xxe = cov[kXX]; + double zze = cov[kZZ]; + double xze = cov[kXZ]; + // + if (fullErr2D) { + xxe += par[kBX]*par[kBX]*cov[kYY]-2.*par[kBX]*cov[kXY]; + zze += par[kBZ]*par[kBZ]*cov[kYY]-2.*par[kBZ]*cov[kZY]; + xze += par[kBX]*par[kBZ]*cov[kYY]-cov[kYZ]*par[kBZ]-cov[kXY]*par[kBX]; + } + // + double det = xxe*zze - xze*xze; + if (det %+e\n",ip,dx,dz,xxeI,zzeI,xzeI, chi2); + } + // +} + +void AliITSAlignMille2::InitTrackParams(int meth) +{ + /// initialize local parameters with different methods + /// for current track (fTrack) + Int_t npts=0; + AliTrackPoint ap; + double sX=0,sXY=0,sZ=0,sZY=0,sY=0,sYY=0,det=0; + // 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) + npts = fTrack->GetNPoints(); + AliDebug(3,Form("*** initializing track with %d points ***",npts)); + for (int i=npts;i--;) { + sY += fTrack->GetY()[i]; + sYY += fTrack->GetY()[i]*fTrack->GetY()[i]; + sX += fTrack->GetX()[i]; + sXY += fTrack->GetX()[i]*fTrack->GetY()[i]; + sZ += fTrack->GetZ()[i]; + sZY += fTrack->GetZ()[i]*fTrack->GetY()[i]; + } + det = sYY*npts-sY*sY; + if (det==0) det = 1E-20; + fLocalInitParam[0] = (sX*sYY-sY*sXY)/det; + fLocalInitParam[2] = (sXY*npts-sY*sX)/det; + // + fLocalInitParam[1] = (sZ*sYY-sY*sZY)/det; + fLocalInitParam[3] = (sZY*npts-sY*sZ)/det; + //AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1))); + AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f ugx = %f\n",fLocalInitParam[0],fLocalInitParam[2])); + // + if (meth==1) return; + // + // perform full fit accounting for cov.matrix + static TVirtualFitter *minuit = 0; + static Double_t step[5] = {1E-3,1E-3,1E-4,1E-4,1E-5}; + static Double_t arglist[10]; + // + if (!minuit) { + minuit = TVirtualFitter::Fitter(0,4); + minuit->SetFCN(trackFit2D); + arglist[0] = 1; + minuit->ExecuteCommand("SET ERR",arglist, 1); + // + arglist[0] = -1; + minuit->ExecuteCommand("SET PRINT",arglist,1); + // + } + // + minuit->SetParameter(0, "ax", fLocalInitParam[0], step[0], 0,0); + minuit->SetParameter(1, "az", fLocalInitParam[1], step[1], 0,0); + minuit->SetParameter(2, "bx", fLocalInitParam[2], step[2], 0,0); + minuit->SetParameter(3, "bz", fLocalInitParam[3], step[3], 0,0); + // + arglist[0] = 1000; // number of function calls + arglist[1] = 0.001; // tolerance + fullErr2D = kTRUE; + minuit->ExecuteCommand("MIGRAD",arglist,2); + // + for (int i=0;i<4;i++) fLocalInitParam[i] = minuit->GetParameter(i); + for (int i=0;i<4;i++) for (int j=0;j<4;j++) fLocalInitParEr[i][j] = minuit->GetCovarianceMatrixElement(i,j); + // +} + + +Int_t AliITSAlignMille2::IsDefined(UShort_t voluid) const +{ + // checks if supermodule 'voluid' is defined and return the internal index + // return -1 if error + Int_t k = fNModules-1; + while ( k>=0 && !(voluid==GetMilleModule(k)->GetVolumeID()) ) k--; + if (k<0) return -1; + return k; +} + +Int_t AliITSAlignMille2::IsContained(UShort_t voluid) const +{ + // checks if the sensitive module 'voluid' is contained inside a supermodule and return the internal index of the last identified supermodule + // return -1 if error + if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1; + Int_t k=fNModules-1; + while (k>=0 && !(GetMilleModule(k)->IsIn(voluid)) ) k--; + if (k<0) return -1; + return k; +} + +Bool_t AliITSAlignMille2::CheckVolumeID(UShort_t voluid) const +{ + /// check if a sensitive volume is contained inside one of the defined supermodules + Int_t k=fNModules-1; + while (k>=0 && !(GetMilleModule(k)->IsIn(voluid)) ) k--; + if (k>=0) return kTRUE; + return kFALSE; +} + +Int_t AliITSAlignMille2::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++; + } + } + + if (ngoodptsGetNPoints(); + AliDebug(2,Form("*** Input track with %d points ***",npts)); + + // preprocessing of the input track: keep only points in defined volumes, + // move points if prealignment is set, sort by Yglo if required + + fTrack=PrepareTrack(track); + if (!fTrack) return -1; + + npts = fTrack->GetNPoints(); + AliDebug(2,Form("*** Processing prepared track with %d points ***",npts)); + + if (!fBOn) { // straight lines + // 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(fInitTrackParamsMeth); + } + else { + // local parms (fLocalInitParam[]) are the Riemann Fitter params + if (!InitRiemanFit()) { + AliInfo("Riemann fit failed! skipping this track..."); + delete fTrack; + fTrack=NULL; + return -5; + } + } + + Int_t nloceq=0; + static Mille2Data md[100]; + + for (Int_t ipt=0; iptGetPoint(fCluster,ipt); + AliDebug(2,Form("\n--- processing point %d --- \n",ipt)); + + // set geometry parameters for the the current module + if (InitModuleParams()) continue; + AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n", + track->GetVolumeID()[ipt], fCurrentModule->GetIndex(), + fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) )); + AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ())); + + if (!AddLocalEquation(md[nloceq])) nloceq++; + else fTotBadLocEqPoints++; + } // end loop over points + // + delete fTrack; + fTrack=NULL; + // not enough good points! + if (nloceqSaveRecordData(); // RRR + // + return 0; +} + +Int_t AliITSAlignMille2::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, ugx, ugz OR riemann fitters pars + /// 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 %g\ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],lpar[4],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5])); + AliDebug(3,Form("deltalpar = %g %g %g %g %g\n",lpar[0]-fLocalInitParam[0],lpar[1]-fLocalInitParam[1],lpar[2]-fLocalInitParam[2],lpar[3]-fLocalInitParam[3],lpar[4]-fLocalInitParam[4])); + + + // prepare the TGeoHMatrix + TGeoHMatrix *fTempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar, + !fUseGlobalDelta); + if (!fTempHMat) return -1; + + Double_t v0g[3]; // vector with straight line direction in global coord. + Double_t p0g[3]; // point of the straight line (glo) + + if (fBOn) { // B FIELD! + AliTrackPoint prf; + for (int ip=0; ip<5; ip++) + fRieman->SetParam(ip,lpar[ip]); + + if (!fRieman->GetPCA(fCluster,prf)) { + AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID())); + return -3; + } + // now determine straight line passing tangent to fit curve at prf + // ugx = dX/dY_glo = DeltaX/DeltaY_glo + // mo' P1=(X,Y,Z)_glo_prf + // => (x,y,Z)_trk_prf ruotando di alpha... + Double_t alpha=fRieman->GetAlpha(); + Double_t x1g = prf.GetX(); + Double_t y1g = prf.GetY(); + Double_t z1g = prf.GetZ(); + Double_t x1t = x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha); + Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha); + Double_t z1t = z1g; + + Double_t x2t = x1t+1.0; + Double_t y2t = y1t+fRieman->GetDYat(x1t); + Double_t z2t = z1t+fRieman->GetDZat(x1t); + Double_t x2g = x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha); + Double_t y2g = x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha); + Double_t z2g = z2t; + + AliDebug(3,Form("Riemann frame: fAlpha = %f = %f ",alpha,alpha*180./TMath::Pi())); + AliDebug(3,Form(" prf_glo=( %f , %f , %f ) prf_rf=( %f , %f , %f )\n", x1g,y1g,z1g, x1t,y1t,z1t)); + AliDebug(3,Form(" mov_glo=( %f , %f , %f ) rf=( %f , %f , %f )\n",x2g,y2g,z2g, x2t,y2t,z2t)); + + if (TMath::Abs(y2g-y1g)<1e-15) { + AliInfo("DeltaY=0! Cannot proceed..."); + return -1; + } + // ugx,1,ugz + v0g[0] = (x2g-x1g)/(y2g-y1g); + v0g[1]=1.0; + v0g[2] = (z2g-z1g)/(y2g-y1g); + + // point: just keep prf + p0g[0]=x1g; + p0g[1]=y1g; + p0g[2]=z1g; + } + else { // staight line + // vector of initial straight line direction in glob. coord + v0g[0]=lpar[2]; + v0g[1]=1.0; + v0g[2]=lpar[3]; + + // intercept in yg=0 plane in glob coord + p0g[0]=lpar[0]; + p0g[1]=0.0; + p0g[2]=lpar[1]; + } + AliDebug(3,Form("Line vector: ( %f , %f , %f ) point:( %f , %f , %f )\n",v0g[0],v0g[1],v0g[2],p0g[0],p0g[1],p0g[2])); + + // 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 AliITSAlignMille2::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 OR riemann's params + /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg) + /// return 0 if success + + // copy initial parameters + Double_t lpar[ITSMILLE2_NLOCAL]; + Double_t gpar[ITSMILLE2_NPARCH]; + for (Int_t i=0; iGetUniqueID(); + nLev++; + // + } while( (fCurrentModule=fCurrentModule->GetParent()) ); + // + // ok, can copy to m + AliDebug(2,Form("Adding local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0])); + m.measX = fMeasLoc[0]-fPintLoc0[0]; + m.sigmaX = fSigmaLoc[0]; + // + AliDebug(2,Form("Adding local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2])); + m.measZ = fMeasLoc[2]-fPintLoc0[2]; + m.sigmaZ = fSigmaLoc[2]; + // + m.levFilled = nLev; + fCurrentModule = endModule; + // + return 0; +} + +void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq) { + /// Set local equations with data stored in m + /// return 0 if success + // + for (Int_t j=0; jIncNProcessedPoints(); + int hlev = m.moduleIDX[il]*ITSMILLE2_NPARCH; // id of the supermodule + int llev = il*ITSMILLE2_NPARCH; + for (int i=ITSMILLE2_NPARCH; i--;) SetGlobalDerivative( hlev+i, m.dergloX[llev+i] ); + } + // + fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.measX, m.sigmaX); + // + // set equation for Zloc coordinate + AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",m.measZ, m.sigmaZ)); + for (int i=fNLocal; i--;) SetLocalDerivative( i, m.derlocZ[i] ); + for (int il=m.levFilled;il--;) { + int hlev = m.moduleIDX[il]*ITSMILLE2_NPARCH; // id of the supermodule + int llev = il*ITSMILLE2_NPARCH; + for (int i=ITSMILLE2_NPARCH; i--;) SetGlobalDerivative(hlev+i, m.dergloZ[llev+i] ); + } + fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.measZ, m.sigmaZ); + } +} + +Int_t AliITSAlignMille2::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 0; + } + int res = fMillepede->GlobalFit(parameters,errors,pulls); + AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed")); + return res; +} + +Double_t AliITSAlignMille2::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 AliITSAlignMille2::PrintGlobalParameters() { + /// Print global parameters + if (!fIsMilleInit) { + AliInfo("Millepede has not been initialized!"); + return; + } + fMillepede->PrintGlobalParameters(); +} + +// //_________________________________________________________________________ +Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile) +{ + // load definitions of supermodules from a root file + // return 0 if success + + TFile *smf=TFile::Open(sfile); + if (!smf->IsOpen()) { + AliInfo(Form("Cannot open supermodule file %s",sfile)); + return -1; + } + + TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules"); + if (!sma) { + AliInfo(Form("Cannot find ITSMilleSuperModules array in file")); + return -2; + } + Int_t nsma=sma->GetEntriesFast(); + AliInfo(Form("Array of SuperModules with %d entries\n",nsma)); + // + Char_t st[250]; + char symname[150]; + UShort_t volid; + TGeoHMatrix m; + // + for (Int_t i=0; iUncheckedAt(i); + volid=a->GetVolUID(); + strcpy(st,a->GetSymName()); + a->GetMatrix(m); + // + sscanf(st,"%s",symname); + // decode module list + char *stp=strstr(st,"ModuleList:"); + if (!stp) return -3; + stp += 11; + int idx[2200]; + char spp[200]; int jp=0; + char cl[20]; + strcpy(st,stp); + int l=strlen(st); + int j=0; + int n=0; + // + while (j<=l) { + if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) { + spp[jp]=0; + jp=0; + if (strlen(spp)) { + int k=strcspn(spp,"-"); + if (k%d)",idx[j],kMaxITSSensID)); + return -5; + } + } + Int_t smindex=int(2198+volid-14336); // virtual index + fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules); + // + fNSuperModules++; + } + + smf->Close(); + // + return 0; +} + +//_________________________________________________________________________ +void AliITSAlignMille2::ConstrainModuleSubUnits(Int_t idm, Double_t val, UInt_t pattern) +{ + // require that sum of modifications for the childs of this module is = val, i.e. + // the internal corrections moves the module as a whole by fixed value (0 by default). + // pattern is the bit pattern for the parameters to constrain + // + static TObjArray childs; + childs.Clear(); + // + // build list of childs for this module + int nChilds = 0; + AliITSAlignMille2Module* parent = GetMilleModule(idm); + if (!parent) return; + // + for (int i=fNModules;i--;) { + AliITSAlignMille2Module* child = GetMilleModule(i); + if (child->GetParent() == parent) childs.AddAtAndExpand(child,nChilds++); + } + if (nChilds<1) return; + // + int npc = 0; + for (int ip=0;ip>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue; + ResetLocalEquation(); + for (int ich=nChilds;ich--;) fGlobalDerivatives[childs[ich]->GetUniqueID()*kNParCh+ip] = 1.0; + AddConstraint(fGlobalDerivatives,val); + npc++; + } + // + AliInfo(Form("Constrained %d params for %d submodules of module #%d: %s",npc,nChilds,idm,parent->GetName())); + // +} + +//_________________________________________________________________________ +void AliITSAlignMille2::PostConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern) +{ + // require that median of the modifications for the childs of this module is = val, i.e. + // the internal corrections moves the module as a whole by fixed value (0 by default) + // module the outliers. + // pattern is the bit pattern for the parameters to constrain + // The difference between the mean and the median will be transfered to the parent + // + static TObjArray childs; + childs.Clear(); + // + // build list of childs for this module + int nChilds = 0; + AliITSAlignMille2Module* parent = GetMilleModule(idm); + if (!parent) return; + // + for (int i=fNModules;i--;) { + AliITSAlignMille2Module* child = GetMilleModule(i); + if (child->GetParent() == parent) childs.AddAtAndExpand(child,nChilds++); + } + if (nChilds<1) return; + // + int npc = 0; + double *deltas = fMillepede->GetDeltaPars(); + double *tmpArr = new double[nChilds]; + // + for (int ip=0;ip>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue; + // compute the median of the deltas + for (int ich=nChilds;ich--;) tmpArr[ich] = deltas[childs[ich]->GetUniqueID()*kNParCh+ip]; + for (int ic0=0;ic0tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;} + // + int kmed = nChilds/2; + double median = (tmpArr[kmed]+tmpArr[nChilds-kmed-1])/2.; + // + for (int ich=nChilds;ich--;) deltas[childs[ich]->GetUniqueID()*kNParCh+ip] -= median - val; + deltas[parent->GetUniqueID()*kNParCh+ip] += median - val; + npc++; + } + delete[] tmpArr; + // + AliInfo(Form("Applied median constraint to %d params for %d submodules of module #%d: %s",npc,nChilds,idm,parent->GetName())); + // +} + +//_________________________________________________________________________ +void AliITSAlignMille2::ConstrainOrphans(Double_t val, UInt_t pattern) +{ + // require that median of the modifications for the supermodules which have no parents is = val, i.e. + // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers. + // pattern is the bit pattern for the parameters to constrain + // + static TObjArray modSet; + modSet.Clear(); + // + // build list of childs for this module + int nModules = 0; + int *nFree = new int[kNParCh]; + for (int i=0;iGetParent()) continue; // skip this + for (int ip=0;ip>ip)&0x1) && module->IsFreeDOF(ip)) nFree[ip]++; + modSet.AddAtAndExpand(module,nModules++); + } + if (nModules<1) return; + // + int npc = 0; + for (int ip=0;ip>ip)&0x1) || !module->IsFreeDOF(ip)) continue; + fGlobalDerivatives[module->GetUniqueID()*kNParCh+ip] = 1.0; + } + AddConstraint(fGlobalDerivatives,val); + npc++; + } + // + delete[] nFree; + AliInfo(Form("Constrained %d params for %d orphan modules",npc,nModules)); + // +} +//_________________________________________________________________________ +void AliITSAlignMille2::PostConstrainOrphansMedian(Double_t val, UInt_t pattern) +{ + // require that sum of modifications for the supermodules which have no parents is = val, i.e. + // the corrections moves the whole setup by fixed value (0 by default). + // pattern is the bit pattern for the parameters to constrain + // + static TObjArray modSet; + modSet.Clear(); + // + // build list of childs for this module + int nModules = 0; + int *nFree = new int[kNParCh]; + for (int i=0;iGetParent()) continue; // skip this + for (int ip=0;ip>ip)&0x1) && module->IsFreeDOF(ip)) nFree[ip]++; + modSet.AddAtAndExpand(module,nModules++); + } + if (nModules<1) return; + // + int npc = 0; + double *deltas = fMillepede->GetDeltaPars(); + double *tmpArr = new double[nModules]; + // + for (int ip=0;ipGetUniqueID()*kNParCh+ip]; + for (int ic0=0;ic0tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}; + // + int kmed = nModules/2; + double median = (tmpArr[kmed]+tmpArr[nModules-kmed-1])/2.; + // + for (int ich=nModules;ich--;) { + AliITSAlignMille2Module* module = (AliITSAlignMille2Module*) modSet[ich]; + if ( !((pattern>>ip)&0x1) || !module->IsFreeDOF(ip)) continue; + deltas[module->GetUniqueID()*kNParCh+ip] -= median - val; + } + npc++; + } + // + delete[] nFree; + delete[] tmpArr; + AliInfo(Form("Applied median constraint to %d params for %d orphan modules",npc,nModules)); + // +} + +//_________________________________________________________________________ +void AliITSAlignMille2::ConstrainLinComb(const Int_t *vidLst, const Float_t *wghLst, Int_t nmd, Double_t val, UInt_t pattern) +{ + // require that the linear combinations of the nmd modules (refered by their volume ID) from the + // modList with the coefficients wghLst adds up to val. + // pattern is the bit pattern for the parameters to constrain. + // + static TObjArray modSet; + modSet.Clear(); + // + // build list of childs for this module + int nModules = 0; + int *nFree = new int[kNParCh]; + for (int i=0;iGetVolumeID() == vid) { + modSet.AddAtAndExpand(module,nModules++); + for (int ip=0;ip>ip)&0x1) && module->IsFreeDOF(ip)) nFree[ip]++; + break; + } + } + } + // + if (nModules != nmd) { + AliInfo(Form("Error: constraint for %d modules requested but %d are found",nmd,nModules)); + delete[] nFree; + return; + } + int npc = 0; + for (int ip=0;ip>ip)&0x1) || !module->IsFreeDOF(ip)) continue; + fGlobalDerivatives[module->GetUniqueID()*kNParCh+ip] = wghLst[ich]; + } + AddConstraint(fGlobalDerivatives,val); + npc++; + } + // + delete[] nFree; + AliInfo(Form("Constrained %d params for linerar combination of %d modules",npc,nModules)); + // +} + diff --git a/ITS/AliITSAlignMille2.h b/ITS/AliITSAlignMille2.h new file mode 100644 index 00000000000..b3017950184 --- /dev/null +++ b/ITS/AliITSAlignMille2.h @@ -0,0 +1,251 @@ +#ifndef ALIITSALIGNMILLE2_H +#define ALIITSALIGNMILLE2_H +/* Copyright(c) 2007-2009 , ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + + +/// \ingroup rec +/// \class AliITSAlignMille2 +/// \brief Class for alignment of ITS +// +// Authors: Marcello Lunardon + +#include +#include +#include "AliTrackPointArray.h" +#include "AliITSAlignMille2Module.h" + +class AliMillePede2; +class AliAlignObjParams; +class TGeoManager; +class TGeoHMatrix; +//class AliITSAlignMille2Module; +class AliTrackFitterRieman; +class TVirtualFitter; +// number of used objects + +#define ITSMILLE2_NPARCH 6 +#define ITSMILLE2_NLOCAL 5 + +struct Mille2Data { + /// structure to store data for 2 LocalEquations (X and Z) + // + Int_t moduleIDX[10]; // max 10 hierarchy levels!!! + Int_t levFilled; + // + Double_t measX; + Double_t measZ; + Double_t sigmaX; + Double_t sigmaZ; + // + Double_t derlocX[ITSMILLE2_NLOCAL]; + Double_t derlocZ[ITSMILLE2_NLOCAL]; + // + Double_t dergloX[ITSMILLE2_NPARCH*10]; + Double_t dergloZ[ITSMILLE2_NPARCH*10]; +}; + +class AliITSAlignMille2: public TObject +{ + public: + enum {kNLocal=ITSMILLE2_NLOCAL,kNParCh=ITSMILLE2_NPARCH, + kMaxITSSensID=2197,kMaxITSSensVID=14300,kMinITSSupeModuleID=14336}; + enum {kDOFTX,kDOFTY,kDOFTZ,kDOFPH,kDOFTH,kDOFPS}; + // + public: + // + AliITSAlignMille2(const Char_t *configFilename="AliITSAlignMille.conf", Bool_t initmille=kTRUE); + virtual ~AliITSAlignMille2(); + // + AliMillePede2* GetMillePede() const {return fMillepede;} + + // 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 SetParSigTranslations(double v) {fParSigTranslations = v;} + void SetParSigRotations(double v) {fParSigRotations = v;} + // + // configuration methods + void SetGeometryFileName(const Char_t* filename="geometry.root") { fGeometryFileName = filename; } + const Char_t* GetGeometryFileName() {return fGeometryFileName.Data();} + const Char_t* GetPreAlignmentFileName() {return fPreAlignmentFileName.Data();} + TClonesArray* GetPreAlignmentDeltas() const {return fPrealignment;} + void SetCurrentModule(Int_t id) {fCurrentModule = GetMilleModule(id);} + void PrintCurrentModuleInfo() const {if (fCurrentModule) fCurrentModule->Print();} + void Print(Option_t*) const; + Bool_t IsConfigured() const {return fIsConfigured;} + void SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts); + void SetUseGlobalDelta(Bool_t v=kTRUE) {fUseGlobalDelta = v;} + Bool_t GetUseGlobalDelta() const {return fUseGlobalDelta;} + // + // fitting methods + AliTrackFitterRieman *GetRiemanFitter() const {return fRieman;} + AliTrackPointArray *PrepareTrack(const AliTrackPointArray *track); + Int_t ProcessTrack(const AliTrackPointArray *track); + void InitTrackParams(int meth=1); + void SetMinNPtsPerTrack(Int_t pts=3) {fMinNPtsPerTrack=pts;} + void SetInitTrackParamsMeth(Int_t meth=1) {fInitTrackParamsMeth=meth;} + Bool_t InitRiemanFit(); + Int_t InitModuleParams(); + Int_t CheckCurrentTrack(); + Bool_t CheckVolumeID(UShort_t voluid) const; + Int_t IsDefined(UShort_t voluid) const; + Int_t IsContained(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() const {return (Double_t*)fPintLoc;} + Double_t* GetGlobalIntersectionPoint() const {return (Double_t*)fPintGlo;} + AliTrackPointArray *SortTrack(const AliTrackPointArray *atp); + void SetTemporaryExcludedModule(Int_t index) {fTempExcludedModule=index;} + Int_t GetTemporaryExcludedModule() const {return fTempExcludedModule;} + // + // Hierarchical contraints + void ConstrainModuleSubUnits(Int_t idm, Double_t val=0, UInt_t pattern=0xfffffff); + void ConstrainOrphans(Double_t val=0,UInt_t pattern=0xfffffff); + void ConstrainLinComb(const Int_t *modLst,const Float_t *wghLst, Int_t nmd, Double_t val=0, UInt_t pattern=0xfffffff); + // + void PostConstrainModuleSubUnitsMedian(Int_t idm, Double_t val=0, UInt_t pattern=0xfffffff); + void PostConstrainOrphansMedian(Double_t val=0, UInt_t pattern=0xfffffff); + // + // Double_t* GetModuleConstrDerivs(Int_t matRow,Int_t matCol, TGeoHMatrix &matLoc); + // + // 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;} + // + Int_t GlobalFit(Double_t *parameters=0,Double_t *errors=0,Double_t *pulls=0); + void PrintGlobalParameters(); + Double_t GetParError(Int_t iPar); + Int_t AddLocalEquation(Mille2Data &m); + void SetLocalEquations(const Mille2Data *marr, Int_t neq); + // + // fitting stuffs + AliTrackPointArray *GetCurrentTrack() {return fTrack;} + AliTrackPoint *GetCurrentCluster() {return &fCluster;} + void SetCurrentTrack(AliTrackPointArray *atp) {fTrack=atp;} + void SetCurrentCluster(AliTrackPoint &atp) {fCluster=atp;} + + // geometry stuffs + Int_t GetNModules() const {return fNModules;} + Int_t GetCurrentModuleIndex() const {return fCurrentModule ? fCurrentModule->GetIndex():-1;} + TGeoHMatrix *GetCurrentModuleHMatrix() const {return fCurrentModule ? fCurrentModule->GetMatrix():0;} + Double_t *GetCurrentModuleTranslation() const {return fCurrentModule ? fCurrentModule->GetMatrix()->GetTranslation():0;} + Int_t GetCurrentModuleInternalIndex() const {return fCurrentModule ? fCurrentModule->GetUniqueID():-1;} + Int_t GetTotBadLocEqPoints() const {return fTotBadLocEqPoints;} + // + AliITSAlignMille2Module* GetMilleModuleByVID(UShort_t voluid) const; // get pointer to the defined supermodule + AliITSAlignMille2Module* GetMilleModule(Int_t id) const {return (AliITSAlignMille2Module*)fMilleModule[id];} + AliITSAlignMille2Module* GetCurrentModule() const {return fCurrentModule;} + AliITSAlignMille2Module* GetSuperModule(Int_t id) const {return (AliITSAlignMille2Module*)fSuperModule[id];} + // + // debug stuffs + Double_t *GetMeasLoc() const {return (Double_t*)fMeasLoc;} + Double_t *GetSigmaLoc() const {return (Double_t*)fSigmaLoc;} + Double_t GetBField() const {return fBField;} + Double_t *GetLocalInitParam() const {return (Double_t*)fLocalInitParam;} + Double_t *GetLocalInitParEr() const {return (Double_t*)fLocalInitParEr;} + Double_t GetLocalDX() const {return fDerivativeLoc[0];} + Double_t GetLocalDY() const {return fDerivativeLoc[1];} + Double_t GetLocalDZ() const {return fDerivativeLoc[2];} + Double_t GetLocalDif(int i) const {return fDerivativeLoc[i];} + Double_t GetParSigTranslations() const {return fParSigTranslations;} + Double_t GetParSigRotations() const {return fParSigRotations;} + Int_t GetPreAlignmentQualityFactor(Int_t index) const;// if not prealign. return -1 + void SetBug(Int_t bug) {fBug=bug;} // 1:SSD inversion sens.18-19 + static AliITSAlignMille2* GetInstance() {return fgInstance;} + // + private: + // + // configuration methods + Int_t LoadConfig(const Char_t *cfile="AliITSAlignMille.conf"); + Int_t LoadSuperModuleFile(const Char_t *cfile="ITSMilleSuperModules.root"); + void ResetLocalEquation(); + void InitGeometry(); + Int_t ApplyToGeometry(); + // + // millepede methods + void Init(Int_t nGlobal, Int_t nLocal, Int_t nStdDev); + // + AliITSAlignMille2(const AliITSAlignMille2& rhs); + AliITSAlignMille2& operator=(const AliITSAlignMille2& rhs); + // + protected: + // + // millepede stuffs + AliMillePede2 *fMillepede; // Detector independent alignment class + Double_t fStartFac; // Initial value for chi2 cut + 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; // Flag for initialization + Bool_t fSensorsIn; // Are sensors explicitly provieded by the conf file? + Double_t fParSigTranslations; // init sigma for transl. params [cm] + Double_t fParSigRotations; // init sigma for rot. params [deg] + // + // fitting stuffs + AliITSAlignMille2Module *fCurrentModule; // Current SuperModule index + AliTrackPointArray *fTrack; // pointer to current track + AliTrackPoint fCluster; // current cluster + Double_t *fGlobalDerivatives; // Array of global derivatives + Double_t fLocalDerivatives[kNLocal]; // Array of local deriv. + Double_t fLocalInitParam[kNLocal]; // Array with inital values for local parameters for current track + Double_t fLocalInitParEr[kNLocal][kNLocal];// Array with inital values for local parameters for current track + Double_t fModuleInitParam[kNParCh]; // 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 + Double_t fSigmaFactor[3]; // multiplicative factor for cluster sigmaX,Y,Z + // + Double_t fDerivativeLoc[3]; // localXYZ deriv. + Int_t fMinNPtsPerTrack; // min number of points per track to accept it + Int_t fInitTrackParamsMeth; // method for track fit + Int_t fTotBadLocEqPoints; // total number of reject points because of bad EqLoc + AliTrackFitterRieman *fRieman; // riemann fitter for helices + // + // >> new members + Bool_t fUseGlobalDelta; // intetpret deltas as global + Bool_t fRequirePoints; // required points in specific layers + Int_t fNReqLayUp[6]; /// number of points required in layer[n] with Y>0 + Int_t fNReqLayDown[6]; /// number of points required in layer[n] with Y<0 + Int_t fNReqLay[6]; /// number of points required in layer[n] + Int_t fNReqDetUp[3]; /// number of points required in Detector[n] with Y>0 + Int_t fNReqDetDown[3]; /// number of points required in Detector[n] with Y<0 + Int_t fNReqDet[3]; /// number of points required in Detector[n] + Int_t fTempExcludedModule; /// single module temporary excluded from initial fit + // << new members + // + // geometry stuffs + TString fGeometryFileName; // Geometry file name + TString fPreAlignmentFileName; // file with prealigned objects + TGeoManager *fGeoManager; + Bool_t fIsConfigured; + TArrayS fPreAlignQF; + // + TClonesArray* fPrealignment; // array of prealignment global deltas + TObjArray fMilleModule; /// array of super modules to be aligned + TObjArray fSuperModule; /// array of super modules defined in supermodule file + Int_t fNModules; // number of defined modules from config file + Int_t fNSuperModules; /// number of custom supermodules in SM file + Bool_t fUsePreAlignment; // start from prealigned setup + Bool_t fUseSortedTracks; // default is kTRUE + Bool_t fBOn; // magentic field ON + Double_t fBField; // value of magnetic field + Int_t fBug; /// tag for temporary bug correction + // + static AliITSAlignMille2* fgInstance; // global pointer on itself + static Int_t fgInstanceID; // global counter of the instances + ClassDef(AliITSAlignMille2, 0) + +}; + +#endif diff --git a/ITS/AliITSAlignMille2Module.cxx b/ITS/AliITSAlignMille2Module.cxx new file mode 100644 index 00000000000..574d58276cf --- /dev/null +++ b/ITS/AliITSAlignMille2Module.cxx @@ -0,0 +1,688 @@ +/************************************************************************** + * 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 AliITSAlignMille2Module +/// Alignment class for the ALICE ITS detector +/// +/// This class is used by AliITSAlignMille to build custom supermodules +/// made of ITS sensitive modules. These supermodules are then aligned +/// +/// Custom supermodules must have VolumeID > 14335 +/// +/// \author M. Lunardon +//----------------------------------------------------------------------------- + +#include +#include + +#include "AliITSAlignMille2Module.h" +#include "AliITSgeomTGeo.h" +#include "AliGeomManager.h" +#include "AliAlignObjParams.h" +#include "AliLog.h" +#include "AliITSAlignMille2.h" + +/// \cond CLASSIMP +ClassImp(AliITSAlignMille2Module) +/// \endcond + +#define CORHW_ + +AliAlignObjParams AliITSAlignMille2Module::fgTempAlignObj; + +//------------------------------------------------------------- +AliITSAlignMille2Module::AliITSAlignMille2Module() : + TNamed(), + fNSensVol(0), + fIndex(-1), + fVolumeID(0), + fNProcPoints(0), + fSensVolIndex(0), + fSensVolVolumeID(0), + fMatrix(NULL), + fSensVolMatrix(NULL), + fSensVolModifMatrix(NULL), + fParent(NULL) +{ + /// void constructor + fMatrix = new TGeoHMatrix; + fSensVolMatrix = new TGeoHMatrix; + fSensVolModifMatrix = new TGeoHMatrix; + fSensVolIndex.Set(1); + fSensVolVolumeID.Set(1); + fSigmaFactor[0]=fSigmaFactor[1]=fSigmaFactor[2]=1.0; +} + +//------------------------------------------------------------- +AliITSAlignMille2Module::AliITSAlignMille2Module(Int_t index,UShort_t volid,char* symname,TGeoHMatrix *m,Int_t nsv,UShort_t *volidsv) : + TNamed(), + fNSensVol(0), + fIndex(-1), + fVolumeID(0), + fNProcPoints(0), + fSensVolIndex(0), + fSensVolVolumeID(0), + fMatrix(NULL), + fSensVolMatrix(NULL), + fSensVolModifMatrix(NULL), + fParent(NULL) +{ + /// void constructor + fMatrix = new TGeoHMatrix; + fSensVolMatrix = new TGeoHMatrix; + fSensVolModifMatrix = new TGeoHMatrix; + fSigmaFactor[0]=fSigmaFactor[1]=fSigmaFactor[2]=1.0; + if (Set(index,volid,symname,m,nsv,volidsv)) { + AliInfo("Error in AliITSAlignMille2Module::Set() - initializing void supermodule..."); + } +} + +//------------------------------------------------------------- +AliITSAlignMille2Module::AliITSAlignMille2Module(UShort_t volid) : + TNamed(), + fNSensVol(0), + fIndex(-1), + fVolumeID(0), + fNProcPoints(0), + fSensVolIndex(0), + fSensVolVolumeID(0), + fMatrix(NULL), + fSensVolMatrix(NULL), + fSensVolModifMatrix(NULL), + fParent(NULL) +{ + /// simple constructor building a supermodule from a single sensitive volume + fMatrix = new TGeoHMatrix; + fSensVolMatrix = new TGeoHMatrix; + fSensVolModifMatrix = new TGeoHMatrix; + // temporary align object, just use the rotation... + fSensVolIndex.Set(1); + fSensVolVolumeID.Set(1); + fSigmaFactor[0]=fSigmaFactor[1]=fSigmaFactor[2]=1.0; + // + fIndex = GetIndexFromVolumeID(volid); + if (fIndex>=0 && gGeoManager) { // good sensitive module and geometry loaded + SetName(AliGeomManager::SymName(volid)); + fVolumeID = volid; + AddSensitiveVolume(volid); + if (SensVolMatrix(volid, fMatrix)) + AliInfo("Matrix not defined"); + } + else { + AliInfo("Wrong VolumeID or Geometry not loaded - initializing void supermodule..."); + } +} + + +//_____________________________________________________________________________ +AliITSAlignMille2Module::AliITSAlignMille2Module(const AliITSAlignMille2Module &m) : + TNamed(m), + fNSensVol(m.fNSensVol), + fIndex(m.fIndex), + fVolumeID(m.fVolumeID), + fNProcPoints(0), + fSensVolIndex(m.fSensVolIndex), + fSensVolVolumeID(m.fSensVolVolumeID), + fMatrix(new TGeoHMatrix(*m.GetMatrix())), + fSensVolMatrix(new TGeoHMatrix), + fSensVolModifMatrix(new TGeoHMatrix), + fParent(m.fParent) +{ + // Copy constructor + fSensVolIndex = m.fSensVolIndex; + fSensVolVolumeID = m.fSensVolVolumeID; + for (int i=3;i--;) fSigmaFactor[i] = m.fSigmaFactor[i]; +} + +//_____________________________________________________________________________ +AliITSAlignMille2Module& AliITSAlignMille2Module::operator=(const AliITSAlignMille2Module &m) +{ + // operator = + // + if(this==&m) return *this; + ((TNamed *)this)->operator=(m); + // + fNSensVol=m.fNSensVol; + fIndex=m.fIndex; + fVolumeID=m.fVolumeID; + for (int i=3;i--;) fSigmaFactor[i] = m.fSigmaFactor[i]; + if (fMatrix) delete fMatrix; + fMatrix=new TGeoHMatrix(*m.GetMatrix()); + fSensVolIndex = m.fSensVolIndex; + fSensVolVolumeID = m.fSensVolVolumeID; + fParent = m.fParent; + return *this; +} + + +//------------------------------------------------------------- +AliITSAlignMille2Module::~AliITSAlignMille2Module() { + /// Destructor + delete fMatrix; + delete fSensVolMatrix; + delete fSensVolModifMatrix; +} + +//------------------------------------------------------------- +Int_t AliITSAlignMille2Module::Set(Int_t index, UShort_t volid, char* symname, TGeoHMatrix *m, Int_t nsv, UShort_t *volidsv) +{ + // initialize a custom supermodule + // index, volid, symname and matrix must be given + // if (volidsv) add nsv sensitive volumes to the supermodules + // return 0 if success + + if (index<2198) { + AliInfo("Index must be >= 2198"); + return -1; + } + if (volid<14336) { + AliInfo("VolumeID must be >= 14336"); + return -2; + } + + if (!symname) return -3; + for (Int_t i=0; i<2198; i++) { + if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) { + AliInfo("Symname already used by a Sensitive Volume"); + return -3; + } + } + + if (!m) return -4; + + // can initialize needed stuffs + fIndex = index; + fVolumeID = volid; + SetName(symname); + // + (*fMatrix) = (*m); + // + fSensVolIndex.Set(nsv); + fSensVolVolumeID.Set(nsv); + // add sensitive volumes + for (Int_t i=0; i6) return -1; + Int_t idx=Int_t(voluid)-2048*lay; + if (idx>=AliGeomManager::LayerSize(lay)) return -1; + for (Int_t ilay=1; ilay=parent->GetNSensitiveVolumes()) return kFALSE; + return parent->IsIn( fSensVolVolumeID[0] ); +} + +//------------------------------------------------------------- +TGeoHMatrix *AliITSAlignMille2Module::GetSensitiveVolumeModifiedMatrix(UShort_t voluid, Double_t *delta,Bool_t local) +{ + // modify the original TGeoHMatrix of the sensitive module 'voluid' according + // with a delta transform. applied to the supermodule matrix + // return NULL if error + + if (!IsIn(voluid)) return NULL; + if (!gGeoManager) return NULL; + + // prepare the TGeoHMatrix + Double_t tr[3],ang[3]; + tr[0]=delta[0]; // in centimeter + tr[1]=delta[1]; + tr[2]=delta[2]; + ang[0]=delta[3]; // psi (X) in deg + ang[1]=delta[4]; // theta (Y) + ang[2]=delta[5]; // phi (Z) + // + fgTempAlignObj.SetRotation(ang[0],ang[1],ang[2]); + fgTempAlignObj.SetTranslation(tr[0],tr[1],tr[2]); + AliDebug(3,Form("Delta angles: psi=%f theta=%f phi=%f",ang[0],ang[1],ang[2])); + TGeoHMatrix hm; + fgTempAlignObj.GetMatrix(hm); + //printf("\n0: delta matrix\n");hm.Print(); + + // 1) start setting fSensVolModif = fSensVol + if (SensVolMatrix(voluid, fSensVolModifMatrix)) return NULL; + // + if (local) { + // 2) set fSensVolModif = SensVolRel + fSensVolModifMatrix->MultiplyLeft( &fMatrix->Inverse() ); + // 3) multiply left by delta + fSensVolModifMatrix->MultiplyLeft( &hm ); + // 4) multiply left by fMatrix + fSensVolModifMatrix->MultiplyLeft( fMatrix ); + } + else fSensVolModifMatrix->MultiplyLeft( &hm ); + // + return fSensVolModifMatrix; +} + +//------------------------------------------------------------- +AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeMisalignment(UShort_t voluid, Double_t *deltalocal) +{ + // calculate misalignment of sens.vol. 'voluid' according with a displacement 'deltalocal' + // of the mother volume. The misalignment is returned as AliAlignObjParams object + + if (!IsIn(voluid)) return NULL; + if (!gGeoManager) return NULL; + + // prepare the TGeoHMatrix + Double_t tr[3],ang[3]; + tr[0]=deltalocal[0]; // in centimeter + tr[1]=deltalocal[1]; + tr[2]=deltalocal[2]; + ang[0]=deltalocal[3]; // psi (X) in deg + ang[1]=deltalocal[4]; // theta (Y) + ang[2]=deltalocal[5]; // phi (Z) + // + fgTempAlignObj.SetRotation(ang[0],ang[1],ang[2]); + fgTempAlignObj.SetTranslation(tr[0],tr[1],tr[2]); + AliDebug(3,Form("Delta angles: psi=%f theta=%f phi=%f",ang[0],ang[1],ang[2])); + // + return GetSensitiveVolumeMisalignment(voluid,&fgTempAlignObj); +} + +//------------------------------------------------------------- +AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeMisalignment(UShort_t voluid, AliAlignObjParams *a) +{ + // return the misalignment of the sens. vol. 'voluid' corresponding with + // a misalignment 'a' in the mother volume + // return NULL if error + + // Gsv = Gg * Gg-1 * Gsv -> Lsv,g = Gg-1 * Gsv + // G'sv = Gg * Dg * Lsv,g === Gsv * Dsv + // Gg * Dg * Gg-1 * Gsv = Gsv * Gsv-1 * Gg * Dg * Gg-1 * Gsv + // + // => Dsv = (Gsv-1 * Gg * Dg * Gg-1 * Gsv) + // + + if (!IsIn(voluid)) return NULL; + if (!gGeoManager) return NULL; + + //a->Print(""); + + // prepare the Delta matrix Dg + TGeoHMatrix dg; + a->GetMatrix(dg); + //dg.Print(); + + // 1) start setting fSensVolModif = Gsv + if (SensVolMatrix(voluid, fSensVolModifMatrix)) return NULL; + //printf("\n1: modif=orig del sensvol\n");fSensVolModifMatrix->Print(); + + // 2) set fSensVolModif = Gg-1 * Gsv + fSensVolModifMatrix->MultiplyLeft( &fMatrix->Inverse() ); + //printf("\n2: modif=relative del sensvol\n");fSensVolModifMatrix->Print(); + + // 3) set fSensVolModif = Dg * Gg-1 * Gsv + fSensVolModifMatrix->MultiplyLeft( &dg ); + //printf("\n3: modif= delta*relative\n");fSensVolModifMatrix->Print(); + + // 4) set fSensVolModif = Gg * Dg * Gg-1 * Gsv + fSensVolModifMatrix->MultiplyLeft( fMatrix ); + //printf("\n4: modif=quasi,manca il Gsv-1...\n");fSensVolModifMatrix->Print(); + + // 5) set fSensVolModif = Gsv-1 * Gg * Dg * Gg-1 * Gsv + if (SensVolMatrix(voluid, &dg)) return NULL; + fSensVolModifMatrix->MultiplyLeft( &dg.Inverse() ); + //printf("\n5: modif=finale\n");fSensVolModifMatrix->Print(); + // + // >> RS + // 6) mo' fSensVolModif dovrebbe essere la Dsv(loc) t.c. G'sv = Gsv*Dsv(loc) + // per trasformarla in Dsv(loc rispetto al Gsv0, non modificato) dovrebbe essere: + // Dsv(loc) -> Dpre * Dsv(loc) * Dpre-1 + //TGeoHMatrix dpre; // dpre = Gsv0-1*Gsv + //if (SensVolOrigGlobalMatrix(voluid, &dg)) return NULL; + //if (SensVolMatrix(voluid, &dpre)) return NULL; + //dpre.MultiplyLeft( &dg.Inverse() ); + //fSensVolModifMatrix->Multiply( &dpre.Inverse() ); + //fSensVolModifMatrix->MultiplyLeft( &dpre ); + // direi che NON FUNZIONA!!!! + + // << RS + + // reset align object (may not be needed...) + fgTempAlignObj.SetVolUID(0); + fgTempAlignObj.SetSymName(""); + fgTempAlignObj.SetTranslation(0,0,0); + fgTempAlignObj.SetRotation(0,0,0); + // + // >> RS +#ifdef CORHW_ + // correction for SPD y-shift + if (voluid>=2048 && voluid<4256) { + TGeoHMatrix deltay; + double dy[3]={0.,0.0081,0.}; + deltay.SetTranslation(dy); + fSensVolModifMatrix->MultiplyLeft( &deltay ); + fSensVolModifMatrix->Multiply( &deltay.Inverse() ); + } +#endif + // << RS + if (!fgTempAlignObj.SetMatrix(*fSensVolModifMatrix)) return NULL; + fgTempAlignObj.SetVolUID(voluid); + fgTempAlignObj.SetSymName(AliGeomManager::SymName(voluid)); + // + return &fgTempAlignObj; +} + +// >> RS +//------------------------------------------------------------- +AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeTotalMisalignment(UShort_t voluid, Double_t *deltalocal) +{ + // calculate misalignment of sens.vol. 'voluid' according with a displacement 'deltalocal' + // of the mother volume. The misalignment is returned as AliAlignObjParams object including + // the (evenctual) prealignment => no merging needed + + if (!IsIn(voluid)) return NULL; + if (!gGeoManager) return NULL; + + // prepare the TGeoHMatrix + Double_t tr[3],ang[3]; + tr[0]=deltalocal[0]; // in centimeter + tr[1]=deltalocal[1]; + tr[2]=deltalocal[2]; + ang[0]=deltalocal[3]; // psi (X) in deg + ang[1]=deltalocal[4]; // theta (Y) + ang[2]=deltalocal[5]; // phi (Z) + + // reset align object (may not be needed...) + fgTempAlignObj.SetVolUID(0); + fgTempAlignObj.SetSymName(""); + fgTempAlignObj.SetRotation(ang[0],ang[1],ang[2]); + fgTempAlignObj.SetTranslation(tr[0],tr[1],tr[2]); + AliDebug(3,Form("Delta angles: psi=%f theta=%f phi=%f",ang[0],ang[1],ang[2])); + + // Gsv = Gg * Gg-1 * Gsv -> Lsv,g = Gg-1 * Gsv + // G'sv = Gg * Dg * Lsv,g === DGsv * Gsv + // + // => Dsv = (G0sv-1 * Gg * Dg * Gg-1 * GMsv) // + // + + // prepare the Delta matrix Dg + TGeoHMatrix dg; + fgTempAlignObj.GetMatrix(dg); + //dg.Print(); + + // 1) start setting fSensVolModif = Gsv + if (SensVolMatrix(voluid, fSensVolModifMatrix)) return NULL; + //printf("\n1: modif=orig del sensvol\n");fSensVolModifMatrix->Print(); + + // 2) set fSensVolModif = Gg-1 * Gsv + fSensVolModifMatrix->MultiplyLeft( &fMatrix->Inverse() ); + //printf("\n2: modif=relative del sensvol\n");fSensVolModifMatrix->Print(); + + // 3) set fSensVolModif = Dg * Gg-1 * Gsv + fSensVolModifMatrix->MultiplyLeft( &dg ); + //printf("\n3: modif= delta*relative\n");fSensVolModifMatrix->Print(); + + // 4) set fSensVolModif = Gg * Dg * Gg-1 * Gsv + fSensVolModifMatrix->MultiplyLeft( fMatrix ); + //printf("\n4: modif=quasi,manca il Gsv-1...\n");fSensVolModifMatrix->Print(); + + // 5) set fSensVolModif = G0sv-1 * Gg * Dg * Gg-1 * Gsv + // qui usa l'orig anziche' la prealigned... + if (SensVolOrigGlobalMatrix(voluid, &dg)) return NULL; + fSensVolModifMatrix->MultiplyLeft( &dg.Inverse() ); + //printf("\n5: modif=finale\n");fSensVolModifMatrix->Print(); + + // reset align object (may not be needed...) + fgTempAlignObj.SetVolUID(0); + fgTempAlignObj.SetSymName(""); + fgTempAlignObj.SetTranslation(0,0,0); + fgTempAlignObj.SetRotation(0,0,0); + +#ifdef CORHW_ + // correction for SPD y-shift + if (voluid>=2048 && voluid<4256) { + TGeoHMatrix deltay; + double dy[3]={0.,0.0081,0.}; + deltay.SetTranslation(dy); + fSensVolModifMatrix->MultiplyLeft( &deltay ); + fSensVolModifMatrix->Multiply( &deltay.Inverse() ); + } +#endif + if (!fgTempAlignObj.SetMatrix(*fSensVolModifMatrix)) return NULL; + fgTempAlignObj.SetVolUID(voluid); + fgTempAlignObj.SetSymName(AliGeomManager::SymName(voluid)); + + + //fgTempAlignObj.Print(""); + + return &fgTempAlignObj; +} +//------------------------------------------------------------- + +//------------------------------------------------------------- +AliAlignObjParams *AliITSAlignMille2Module::GetSensitiveVolumeGlobalMisalignment(UShort_t voluid, Double_t *deltalocal) +{ + // calculate misalignment of sens.vol. 'voluid' according with a displacement 'deltalocal' + // of the mother volume. The misalignment is returned as AliAlignObjParams object + + if (!IsIn(voluid)) return NULL; + if (!gGeoManager) return NULL; + + // prepare the TGeoHMatrix + Double_t tr[3],ang[3]; + tr[0]=deltalocal[0]; // in centimeter + tr[1]=deltalocal[1]; + tr[2]=deltalocal[2]; + ang[0]=deltalocal[3]; // psi (X) in deg + ang[1]=deltalocal[4]; // theta (Y) + ang[2]=deltalocal[5]; // phi (Z) + + // reset align object (may not be needed...) + fgTempAlignObj.SetTranslation(0,0,0); + fgTempAlignObj.SetRotation(0,0,0); + + fgTempAlignObj.SetRotation(ang[0],ang[1],ang[2]); + fgTempAlignObj.SetTranslation(tr[0],tr[1],tr[2]); + AliDebug(3,Form("Delta angles: psi=%f theta=%f phi=%f",ang[0],ang[1],ang[2])); + + // Gsv = Gg * Gg-1 * Gsv -> Lsv,g = Gg-1 * Gsv + // G'sv = Gg * Dg * Lsv,g === DGsv * Gsv + // + // => DGsv = (Gg * Dg * Gg-1) + // + + // prepare the Delta matrix Dg + TGeoHMatrix dg; + fgTempAlignObj.GetMatrix(dg); + //dg.Print(); + + dg.MultiplyLeft( fMatrix ); + dg.Multiply( &fMatrix->Inverse() ); + + // reset align object (may not be needed...) + fgTempAlignObj.SetTranslation(0,0,0); + fgTempAlignObj.SetRotation(0,0,0); + + fgTempAlignObj.SetVolUID(voluid); + fgTempAlignObj.SetSymName(AliGeomManager::SymName(voluid)); + + if (!fgTempAlignObj.SetMatrix(dg)) return NULL; + + //fgTempAlignObj.Print(""); + + return &fgTempAlignObj; +} +// << RS + +//------------------------------------------------------------- +TGeoHMatrix *AliITSAlignMille2Module::GetSensitiveVolumeMatrix(UShort_t voluid) +{ + // return TGeoHMatrix of the sens.vol. 'voluid' in the current geometry + if (SensVolMatrix(voluid,fSensVolMatrix)) return NULL; + return fSensVolMatrix; +} + +//------------------------------------------------------------- +TGeoHMatrix *AliITSAlignMille2Module::GetSensitiveVolumeOrigGlobalMatrix(UShort_t voluid) +{ + // return original ideal position (from AliGeomManager::GetOrigGlobalMatrix()) + if (SensVolOrigGlobalMatrix(voluid,fSensVolMatrix)) return NULL; + return fSensVolMatrix; +} +//------------------------------------------------------------- +Int_t AliITSAlignMille2Module::SensVolMatrix(UShort_t volid, TGeoHMatrix *m) +{ + // set matrix for sensitive modules (SPD corrected) + // return 0 if success + Double_t rot[9]; + Int_t idx=GetIndexFromVolumeID(volid); + if (idx<0) return -1; + if (!AliITSgeomTGeo::GetRotation(idx,rot)) return -2; + m->SetRotation(rot); + Double_t oLoc[3]={0,0,0}; + Double_t oGlo[3]={0,0,0}; + if (!AliITSgeomTGeo::LocalToGlobal(idx,oLoc,oGlo)) return -3; + m->SetTranslation(oGlo); + return 0; +} + +//------------------------------------------------------------- +Int_t AliITSAlignMille2Module::SensVolOrigGlobalMatrix(UShort_t volid, TGeoHMatrix *m) +{ + // set original global matrix for sensitive modules (SPD corrected) + // return 0 if success + Int_t idx=GetIndexFromVolumeID(volid); + if (idx<0) return -1; + TGeoHMatrix mo; + if (!AliGeomManager::GetOrigGlobalMatrix(AliGeomManager::SymName(volid),mo)); + (*m)=mo; + +#ifdef CORHW_ + // SPD y-shift by 81 mu + if (volid<5000) { + Double_t oLoc[3]={0.0,0.0081,0.0}; + Double_t oGlo[3]={0,0,0}; + m->LocalToMaster(oLoc,oGlo); + m->SetTranslation(oGlo); + } +#endif + return 0; +} + +//------------------------------------------------------------- +UShort_t AliITSAlignMille2Module::GetVolumeIDFromSymname(const Char_t *symname) { + /// volume ID from symname + if (!symname) return 0; + + for (UShort_t voluid=2000; voluid<13300; voluid++) { + Int_t modId; + AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId); + if (layerId>0 && layerId<7 && modId>=0 && modId2197) return 0; + return GetVolumeIDFromSymname(AliITSgeomTGeo::GetSymName(index)); +} + +//------------------------------------------------------------- +void AliITSAlignMille2Module::Print(Option_t*) const +{ + /// + printf("*** ITS SuperModule for AliITSAlignMille ***\n"); + printf("symname : %s\n",GetName()); + printf("parent : %s\n",fParent ? fParent->GetName() : "N/A"); + printf("volumeID : %4d | index : %4d\n",fVolumeID,fIndex); + printf("Factors : X=%.2f Y=%.2f Z=%.2f | DOF: Tx:%d Ty:%d Tz:%d Phi:%d Theta:%d Psi:%d\n", + fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2], + IsFreeDOF(AliITSAlignMille2::kDOFTX),IsFreeDOF(AliITSAlignMille2::kDOFTY), + IsFreeDOF(AliITSAlignMille2::kDOFTZ),IsFreeDOF(AliITSAlignMille2::kDOFPH), + IsFreeDOF(AliITSAlignMille2::kDOFTH),IsFreeDOF(AliITSAlignMille2::kDOFPS)); + fMatrix->Print(); + printf("%4d Sensitive volumes | %6d Processed Points\n",fNSensVol,fNProcPoints); + for (Int_t i=0; iGetAlignableEntry(GetName())!=0; +} + +//------------------------------------------------------------- +void AliITSAlignMille2Module::GetLocalMatrix(TGeoHMatrix &mat) const +{ + // return the local matrix for transformation to its parent + mat = *fMatrix; + if (fParent) mat.MultiplyLeft( &fParent->GetMatrix()->Inverse() ); +} diff --git a/ITS/AliITSAlignMille2Module.h b/ITS/AliITSAlignMille2Module.h new file mode 100644 index 00000000000..7758fe8d08b --- /dev/null +++ b/ITS/AliITSAlignMille2Module.h @@ -0,0 +1,115 @@ +#ifndef ALIITSALIGNMILLE2MODULE_H +#define ALIITSALIGNMILLE2MODULE_H +/* Copyright(c) 2007-2009 , ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/// \ingroup rec +/// \class AliITSAlignMille2Module +/// \brief Class for alignment of ITS +// +// Authors: Marcello Lunardon +// +// RS Converted static arrays fSensVolVolumeID and fSensVolIndex +// to TArrays in user transparent way. +// +/* $Id$ */ +//#include +//#include +#include +#include +#include +class AliITSAlignMille2; + +class AliAlignObjParams; +class TGeoHMatrix; + +class AliITSAlignMille2Module : public TNamed +{ +public: + AliITSAlignMille2Module(); + AliITSAlignMille2Module(UShort_t volid); + AliITSAlignMille2Module(Int_t index, UShort_t volid, char* symname, TGeoHMatrix *m, Int_t nsv=0, UShort_t *volidsv=NULL); + AliITSAlignMille2Module(const AliITSAlignMille2Module& rhs); // copy constructor + AliITSAlignMille2Module& operator=(const AliITSAlignMille2Module& rhs); + // + virtual ~AliITSAlignMille2Module(); + // + // geometry methods + Int_t GetIndex() const {return fIndex;} + UShort_t GetVolumeID() const {return fVolumeID;} + Int_t GetNSensitiveVolumes() const {return fNSensVol;} + Int_t GetSensVolIndex(Int_t at) const {return fSensVolIndex[at];} + Short_t GetSensVolVolumeID(Int_t at) const {return fSensVolVolumeID[at];} + TGeoHMatrix *GetMatrix() const {return fMatrix;} + void GetLocalMatrix(TGeoHMatrix& mat) const; + UShort_t *GetSensitiveVolumeVolumeID() const {return (UShort_t*)fSensVolVolumeID.GetArray();} + Float_t GetSigmaFactor(Int_t i) const {return fSigmaFactor[i];} + Float_t GetSigmaXFactor() const {return fSigmaFactor[0];} + Float_t GetSigmaYFactor() const {return fSigmaFactor[1];} + Float_t GetSigmaZFactor() const {return fSigmaFactor[2];} + Int_t GetNProcessedPoints() const {return fNProcPoints;} + Bool_t IsFreeDOF(Int_t dof) const {return TestBit(1<