X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONAlignment.cxx;h=8c2572812775a09ee82cac2672820ab96c56301c;hb=7e7619d198263c3c576858fd0050f0cb1eb8cd0a;hp=d2077fa3d247b6145ae979a4e054cec50219855f;hpb=3b2890bef2d4aeea6b0803cc5548d2a7892f58fb;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONAlignment.cxx b/MUON/AliMUONAlignment.cxx index d2077fa3d24..8c257281277 100644 --- a/MUON/AliMUONAlignment.cxx +++ b/MUON/AliMUONAlignment.cxx @@ -16,58 +16,61 @@ /* $Id$ */ //----------------------------------------------------------------------------- -/// \class AliMUONAlignmentV5 -/// Alignment class fro the ALICE DiMuon spectrometer +/// \class AliMUONAlignment +/// Alignment class for the ALICE DiMuon spectrometer /// -/// MUON specific alignment class which interface to AliMillepede. +/// MUON 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. +/// at each cluster and fill the corresponding local equations. Provide methods +/// for fixing or constraining detection elements for best results. /// /// \author Bruce Becker, Javier Castillo //----------------------------------------------------------------------------- #include "AliMUONAlignment.h" #include "AliMUONTrack.h" -#include "AliMUONRawCluster.h" #include "AliMUONTrackParam.h" -#include "AliMUONHitForRec.h" +#include "AliMUONVCluster.h" #include "AliMUONGeometryTransformer.h" #include "AliMUONGeometryModuleTransformer.h" #include "AliMUONGeometryDetElement.h" #include "AliMUONGeometryBuilder.h" -#include "AliMUONConstants.h" #include "AliMillepede.h" #include "AliMpExMap.h" +#include "AliMpExMapIterator.h" +#include "AliAlignObjMatrix.h" #include "AliLog.h" -#include "TSystem.h" +#include "TMath.h" +#include "TMatrixDSym.h" +#include "TClonesArray.h" /// \cond CLASSIMP ClassImp(AliMUONAlignment) /// \endcond - Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2; - Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26}; - Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156}; - Int_t AliMUONAlignment::fgNParCh = 3; - Int_t AliMUONAlignment::fgNTrkMod = 16; - Int_t AliMUONAlignment::fgNCh = 10; - Int_t AliMUONAlignment::fgNSt = 5; - -AliMUONAlignment::AliMUONAlignment() +//_____________________________________________________________________ +// static variables +Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2; +Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26}; +Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156}; +Int_t AliMUONAlignment::fgNParCh = 4; +Int_t AliMUONAlignment::fgNTrkMod = 16; +Int_t AliMUONAlignment::fgNCh = 10; +Int_t AliMUONAlignment::fgNSt = 5; + +//_____________________________________________________________________ +AliMUONAlignment::AliMUONAlignment() : TObject(), fBFieldOn(kTRUE), - fStartFac(16.), - fResCutInitial(100.), + fStartFac(256.), + fResCutInitial(100.), fResCut(100.), fMillepede(0), - fTrackParamAtHit(0), - fHitForRecAtHit(0), fTrack(0), - fRecHit(0), + fCluster(0), fTrackParam(0), fNGlobal(fgNDetElem*fgNParCh), fNLocal(4), @@ -77,15 +80,18 @@ AliMUONAlignment::AliMUONAlignment() fPhi(0.), fCosPhi(1.), fSinPhi(0.), + fTrackRecord(), fTransform(0) { - /// Default constructor, setting define alignment parameters - fSigma[0] = 1.0e-1; + + fSigma[0] = 1.5e-1; fSigma[1] = 1.0e-2; - fDoF[0] = kTRUE; fDoF[1] = kTRUE; fDoF[2] = kTRUE; - fAllowVar[0] = 0.05; fAllowVar[1] = 0.05; fAllowVar[2] = 0.001; - + AliInfo(Form("fSigma[0]: %f\t fSigma[1]: %f",fSigma[0],fSigma[1])); + + fDoF[0] = kTRUE; fDoF[1] = kTRUE; fDoF[2] = kTRUE; fDoF[3] = kTRUE; + fAllowVar[0] = 0.05; fAllowVar[1] = 0.05; fAllowVar[2] = 0.001; fAllowVar[3] = 0.5; + AliInfo(Form("fAllowVar[0]: %f\t fAllowVar[1]: %f\t fPhi: %f\t fgNDetElem: %i\t fNGlobal: %i\t fNLocal: %i",fAllowVar[0],fAllowVar[1],fPhi,fgNDetElem,fNGlobal,fNLocal)); fMillepede = new AliMillepede(); @@ -97,30 +103,36 @@ AliMUONAlignment::AliMUONAlignment() } -AliMUONAlignment::~AliMUONAlignment() { - /// Destructor -} +//_____________________________________________________________________ +AliMUONAlignment::~AliMUONAlignment() +{} -void AliMUONAlignment::Init(Int_t nGlobal, /* number of global paramers */ - Int_t nLocal, /* number of local parameters */ - Int_t nStdDev /* std dev cut */ ) +//_____________________________________________________________________ +void AliMUONAlignment::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); - Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE}; - Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE}; + // Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE}; + // Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE}; + // Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE}; - AllowVariations(bStOnOff); + // AllowVariations(bChOnOff); // Fix parameters or add constraints here - for (Int_t iSt=0; iSt<5; iSt++) - if (!bStOnOff[iSt]) FixStation(iSt+1); + // for (Int_t iSt=0; iSt<5; iSt++) + // { if (!bStOnOff[iSt]) FixStation(iSt+1); } + + // for (Int_t iCh=0; iCh<10; iCh++) + // { if (!bChOnOff[iCh]) FixChamber(iCh+1); } - FixHalfSpectrometer(bStOnOff,bSpecLROnOff); +// FixHalfSpectrometer(bChOnOff,bSpecLROnOff); ResetConstraints(); - + // Define global constrains to be applied // X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE}; @@ -130,133 +142,285 @@ void AliMUONAlignment::Init(Int_t nGlobal, /* number of global paramers */ // Other possible way to add constrains bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE; bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE; -// AddConstraints(bStOnOff,bVarXYT,bDetTLBR); + // AddConstraints(bStOnOff,bVarXYT,bDetTLBR); bVarXYT[0] = kTRUE; bVarXYT[1] = kTRUE; bVarXYT[2] = kFALSE; // AddConstraints(bStOnOff,bVarXYT); - + // Set iterations - if (fStartFac>1) fMillepede->SetIterations(fStartFac); + if (fStartFac>1) fMillepede->SetIterations(fStartFac); } -void AliMUONAlignment::FixStation(Int_t iSt){ +//_____________________________________________________________________ +void AliMUONAlignment::FixStation(Int_t iSt) +{ /// Fix all detection elements of station iSt - Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0; - Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1]; - for (Int_t i = iDetElemFirst; i < iDetElemLast; i++){ + Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0; + Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1]; + for (Int_t i = iDetElemFirst; i < iDetElemLast; i++) + { FixParameter(i*fgNParCh+0, 0.0); FixParameter(i*fgNParCh+1, 0.0); FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); } + +} + +//_____________________________________________________________________ +void AliMUONAlignment::FixChamber(Int_t iCh) +{ + /// Fix all detection elements of chamber iCh + Int_t iDetElemFirst = (iCh>1) ? fgSNDetElemCh[iCh-2] : 0; + Int_t iDetElemLast = fgSNDetElemCh[iCh-1]; + for (Int_t i = iDetElemFirst; i < iDetElemLast; i++) + { + FixParameter(i*fgNParCh+0, 0.0); + FixParameter(i*fgNParCh+1, 0.0); + FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); + } +} + +//_____________________________________________________________________ +void AliMUONAlignment::FixDetElem(Int_t iDetElemId, TString sVarXYT) +{ + + /// Fix a given detection element + Int_t iDetElemNumber = iDetElemId%100; + for (int iCh=0; iCh=1 && iCh<=4){ - if ((lDetElemNumber==1 || lDetElemNumber==2) && !lSpecLROnOff[0]){ // From track crossings - FixParameter(i*fgNParCh+0, 0.0); - FixParameter(i*fgNParCh+1, 0.0); - FixParameter(i*fgNParCh+2, 0.0); - } - if ((lDetElemNumber==0 || lDetElemNumber==3) && !lSpecLROnOff[1]){ // From track crossings - FixParameter(i*fgNParCh+0, 0.0); - FixParameter(i*fgNParCh+1, 0.0); - FixParameter(i*fgNParCh+2, 0.0); - } + if (iCh>=1 && iCh<=4) + { + + if ((lDetElemNumber==1 || lDetElemNumber==2) && !lSpecLROnOff[0]) + { + // From track crossings + FixParameter(i*fgNParCh+0, 0.0); + FixParameter(i*fgNParCh+1, 0.0); + FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); + } + + if ((lDetElemNumber==0 || lDetElemNumber==3) && !lSpecLROnOff[1]) + { + // From track crossings + FixParameter(i*fgNParCh+0, 0.0); + FixParameter(i*fgNParCh+1, 0.0); + FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); + } + } - if (iCh>=5 && iCh<=6){ - if ((lDetElemNumber>=5&&lDetElemNumber<=13) && !lSpecLROnOff[0]){ - FixParameter(i*fgNParCh+0, 0.0); - FixParameter(i*fgNParCh+1, 0.0); - FixParameter(i*fgNParCh+2, 0.0); - } - if (((lDetElemNumber>=0&&lDetElemNumber<=4) || - (lDetElemNumber>=14&&lDetElemNumber<=17)) && !lSpecLROnOff[1]){ - FixParameter(i*fgNParCh+0, 0.0); - FixParameter(i*fgNParCh+1, 0.0); - FixParameter(i*fgNParCh+2, 0.0); - } + + if (iCh>=5 && iCh<=6) + { + + if ((lDetElemNumber>=5&&lDetElemNumber<=13) && !lSpecLROnOff[0]) + { + FixParameter(i*fgNParCh+0, 0.0); + FixParameter(i*fgNParCh+1, 0.0); + FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); + } + + if (((lDetElemNumber>=0&&lDetElemNumber<=4) || + (lDetElemNumber>=14&&lDetElemNumber<=17)) && !lSpecLROnOff[1]) + { + + FixParameter(i*fgNParCh+0, 0.0); + FixParameter(i*fgNParCh+1, 0.0); + FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); + } + } - if (iCh>=7 && iCh<=10){ - if ((lDetElemNumber>=7&&lDetElemNumber<=19) && !lSpecLROnOff[0]){ - FixParameter(i*fgNParCh+0, 0.0); - FixParameter(i*fgNParCh+1, 0.0); - FixParameter(i*fgNParCh+2, 0.0); - } - if (((lDetElemNumber>=0&&lDetElemNumber<=6) || - (lDetElemNumber>=20&&lDetElemNumber<=25)) && !lSpecLROnOff[1]){ - FixParameter(i*fgNParCh+0, 0.0); - FixParameter(i*fgNParCh+1, 0.0); - FixParameter(i*fgNParCh+2, 0.0); - } + + if (iCh>=7 && iCh<=10) + { + + if ((lDetElemNumber>=7&&lDetElemNumber<=19) && !lSpecLROnOff[0]) + { + FixParameter(i*fgNParCh+0, 0.0); + FixParameter(i*fgNParCh+1, 0.0); + FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); + } + + if (((lDetElemNumber>=0&&lDetElemNumber<=6) || + (lDetElemNumber>=20&&lDetElemNumber<=25)) && !lSpecLROnOff[1]) + { + FixParameter(i*fgNParCh+0, 0.0); + FixParameter(i*fgNParCh+1, 0.0); + FixParameter(i*fgNParCh+2, 0.0); + FixParameter(i*fgNParCh+3, 0.0); + } + } + } + } + } -void AliMUONAlignment::SetNonLinear(Bool_t *lStOnOff,Bool_t *lVarXYT){ - /// Set non linear parameter flag selected stations and degrees of freedom - for (Int_t i = 0; i < fgNDetElem; i++){ +//______________________________________________________________________ +void AliMUONAlignment::SetNonLinear(const Bool_t *lChOnOff, const Bool_t *lVarXYT) +{ + + /// Set non linear parameter flag selected chambers and degrees of freedom + for (Int_t i = 0; i < fgNDetElem; i++) + { + Int_t iCh=0; - for (iCh=1; iCh<=fgNCh; iCh++){ - if (iLocal2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ, - lDetElemGloX,lDetElemGloY,lDetElemGloZ); + lDetElemGloX,lDetElemGloY,lDetElemGloZ); if (iCh>=1 && iCh<=4){ - if ((lDetElemNumber==1 || lDetElemNumber==2) && lSpecLROnOff[0]){ // From track crossings - lMeanY += lDetElemGloY; - lSigmaY += lDetElemGloY*lDetElemGloY; - lMeanZ += lDetElemGloZ; - lSigmaZ += lDetElemGloZ*lDetElemGloZ; - lNDetElem++; - } - if ((lDetElemNumber==0 || lDetElemNumber==3) && lSpecLROnOff[1]){ // From track crossings - lMeanY += lDetElemGloY; - lSigmaY += lDetElemGloY*lDetElemGloY; - lMeanZ += lDetElemGloZ; - lSigmaZ += lDetElemGloZ*lDetElemGloZ; - lNDetElem++; - } + if ((lDetElemNumber==1 || lDetElemNumber==2) && lSpecLROnOff[0]){ // From track crossings + lMeanY += lDetElemGloY; + lSigmaY += lDetElemGloY*lDetElemGloY; + lMeanZ += lDetElemGloZ; + lSigmaZ += lDetElemGloZ*lDetElemGloZ; + lNDetElem++; + } + if ((lDetElemNumber==0 || lDetElemNumber==3) && lSpecLROnOff[1]){ // From track crossings + lMeanY += lDetElemGloY; + lSigmaY += lDetElemGloY*lDetElemGloY; + lMeanZ += lDetElemGloZ; + lSigmaZ += lDetElemGloZ*lDetElemGloZ; + lNDetElem++; + } } if (iCh>=5 && iCh<=6){ - if ((lDetElemNumber>=5&&lDetElemNumber<=13) && lSpecLROnOff[0]){ - lMeanY += lDetElemGloY; - lSigmaY += lDetElemGloY*lDetElemGloY; - lMeanZ += lDetElemGloZ; - lSigmaZ += lDetElemGloZ*lDetElemGloZ; - lNDetElem++; - } - if (((lDetElemNumber>=0&&lDetElemNumber<=4) || - (lDetElemNumber>=14&&lDetElemNumber<=17)) && lSpecLROnOff[1]){ - lMeanY += lDetElemGloY; - lSigmaY += lDetElemGloY*lDetElemGloY; - lMeanZ += lDetElemGloZ; - lSigmaZ += lDetElemGloZ*lDetElemGloZ; - lNDetElem++; - } + if ((lDetElemNumber>=5&&lDetElemNumber<=13) && lSpecLROnOff[0]){ + lMeanY += lDetElemGloY; + lSigmaY += lDetElemGloY*lDetElemGloY; + lMeanZ += lDetElemGloZ; + lSigmaZ += lDetElemGloZ*lDetElemGloZ; + lNDetElem++; + } + if (((lDetElemNumber>=0&&lDetElemNumber<=4) || + (lDetElemNumber>=14&&lDetElemNumber<=17)) && lSpecLROnOff[1]){ + lMeanY += lDetElemGloY; + lSigmaY += lDetElemGloY*lDetElemGloY; + lMeanZ += lDetElemGloZ; + lSigmaZ += lDetElemGloZ*lDetElemGloZ; + lNDetElem++; + } } if (iCh>=7 && iCh<=10){ - if ((lDetElemNumber>=7&&lDetElemNumber<=19) && lSpecLROnOff[0]){ - lMeanY += lDetElemGloY; - lSigmaY += lDetElemGloY*lDetElemGloY; - lMeanZ += lDetElemGloZ; - lSigmaZ += lDetElemGloZ*lDetElemGloZ; - lNDetElem++; - } - if (((lDetElemNumber>=0&&lDetElemNumber<=6) || - (lDetElemNumber>=20&&lDetElemNumber<=25)) && lSpecLROnOff[1]){ - lMeanY += lDetElemGloY; - lSigmaY += lDetElemGloY*lDetElemGloY; - lMeanZ += lDetElemGloZ; - lSigmaZ += lDetElemGloZ*lDetElemGloZ; - lNDetElem++; - } + if ((lDetElemNumber>=7&&lDetElemNumber<=19) && lSpecLROnOff[0]){ + lMeanY += lDetElemGloY; + lSigmaY += lDetElemGloY*lDetElemGloY; + lMeanZ += lDetElemGloZ; + lSigmaZ += lDetElemGloZ*lDetElemGloZ; + lNDetElem++; + } + if (((lDetElemNumber>=0&&lDetElemNumber<=6) || + (lDetElemNumber>=20&&lDetElemNumber<=25)) && lSpecLROnOff[1]){ + lMeanY += lDetElemGloY; + lSigmaY += lDetElemGloY*lDetElemGloY; + lMeanZ += lDetElemGloZ; + lSigmaZ += lDetElemGloZ*lDetElemGloZ; + lNDetElem++; + } } } } - lMeanY /= lNDetElem; - lSigmaY /= lNDetElem; - lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY); - lMeanZ /= lNDetElem; - lSigmaZ /= lNDetElem; - lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ); - AliInfo(Form("Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ)); - - for (Int_t i = 0; i < fgNDetElem; i++){ + if (lNDetElem) { + lMeanY /= lNDetElem; + lSigmaY /= lNDetElem; + lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY); + lMeanZ /= lNDetElem; + lSigmaZ /= lNDetElem; + lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ); + AliInfo(Form("Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ)); + } else { + AliError("No detection elements to constrain!!!"); + return; + } + + for (Int_t i = 0; i < fgNDetElem; i++){ Int_t iCh=0; for (iCh=1; iCh<=fgNCh; iCh++){ if (iLocal2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ, - lDetElemGloX,lDetElemGloY,lDetElemGloZ); + lDetElemGloX,lDetElemGloY,lDetElemGloZ); if (lVarXYT[0]) { // X constraint - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half } if (lVarXYT[1]) { // Y constraint - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half } if (lVarXYT[2]) { // P constraint - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half } if (lVarXYT[3]) { // X-Z shearing - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXZT,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXZL,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXZB,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXZR,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXZT,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXZL,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXZB,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXZR,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half } if (lVarXYT[4]) { // Y-Z shearing - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYZT,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYZL,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYZB,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYZR,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYZT,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYZL,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYZB,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYZR,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half } if (lVarXYT[5]) { // P-Z rotation - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPZT,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPZL,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPZB,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPZR,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPZT,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPZL,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPZB,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPZR,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half } if (lVarXYT[6]) { // X-Y shearing - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXYT,0,(lDetElemGloY-lMeanY)/lSigmaY); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXYL,0,(lDetElemGloY-lMeanY)/lSigmaY); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXYB,0,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXYR,0,(lDetElemGloY-lMeanY)/lSigmaY); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXYT,0,(lDetElemGloY-lMeanY)/lSigmaY); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXYL,0,(lDetElemGloY-lMeanY)/lSigmaY); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXYB,0,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXYR,0,(lDetElemGloY-lMeanY)/lSigmaY); // Right half } if (lVarXYT[7]) { // Y-Y scaling - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYYT,1,(lDetElemGloY-lMeanY)/lSigmaY); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYYL,1,(lDetElemGloY-lMeanY)/lSigmaY); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYYB,1,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYYR,1,(lDetElemGloY-lMeanY)/lSigmaY); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYYT,1,(lDetElemGloY-lMeanY)/lSigmaY); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYYL,1,(lDetElemGloY-lMeanY)/lSigmaY); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYYB,1,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYYR,1,(lDetElemGloY-lMeanY)/lSigmaY); // Right half } if (lVarXYT[8]) { // P-Y rotation - if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPYT,2,(lDetElemGloY-lMeanY)/lSigmaY); // Top half - if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPYL,2,(lDetElemGloY-lMeanY)/lSigmaY); // Left half - if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPYB,2,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half - if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPYR,2,(lDetElemGloY-lMeanY)/lSigmaY); // Right half + if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPYT,2,(lDetElemGloY-lMeanY)/lSigmaY); // Top half + if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPYL,2,(lDetElemGloY-lMeanY)/lSigmaY); // Left half + if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPYB,2,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half + if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPYR,2,(lDetElemGloY-lMeanY)/lSigmaY); // Right half } } } @@ -462,7 +631,7 @@ void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT, Bool_t * } } -void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar, Double_t lWeight){ +void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar, Double_t /*lWeight*/) const{ /// Set constrain equation for top half of spectrometer Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2]; if (lCh>=1 && lCh<=4){ @@ -482,7 +651,8 @@ void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstrai } } -void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar, Double_t lWeight){ +//______________________________________________________________________ +void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar, Double_t lWeight) const{ /// Set constrain equation for left half of spectrometer Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2]; if (lCh>=1 && lCh<=4){ @@ -502,7 +672,8 @@ void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstrai } } -void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar, Double_t lWeight){ +//______________________________________________________________________ +void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar, Double_t /*lWeight*/) const{ /// Set constrain equation for bottom half of spectrometer Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2]; if (lCh>=1 && lCh<=4){ @@ -511,20 +682,21 @@ void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstrai } } if (lCh>=5 && lCh<=6){ - if ((lDetElemNumber>=9&&lDetElemNumber<=17) || - (lDetElemNumber==0)){ + if ((lDetElemNumber>=9&&lDetElemNumber<=17) || + (lDetElemNumber==0)){ lConstraintB[lDetElem*fgNParCh+iVar]=1.0; } } if (lCh>=7 && lCh<=10){ - if ((lDetElemNumber>=13&&lDetElemNumber<=25) || - (lDetElemNumber==0)){ + if ((lDetElemNumber>=13&&lDetElemNumber<=25) || + (lDetElemNumber==0)){ lConstraintB[lDetElem*fgNParCh+iVar]=1.0; } } } -void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar, Double_t lWeight){ +//______________________________________________________________________ +void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar, Double_t lWeight) const{ /// Set constrain equation for right half of spectrometer Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2]; if (lCh>=1 && lCh<=4){ @@ -533,22 +705,23 @@ void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstrai } } if (lCh>=5 && lCh<=6){ - if ((lDetElemNumber>=0&&lDetElemNumber<=4) || - (lDetElemNumber>=14&&lDetElemNumber<=17)){ + if ((lDetElemNumber>=0&&lDetElemNumber<=4) || + (lDetElemNumber>=14&&lDetElemNumber<=17)){ lConstraintR[lDetElem*fgNParCh+iVar]=lWeight; } } if (lCh>=7 && lCh<=10){ - if ((lDetElemNumber>=0&&lDetElemNumber<=6) || - (lDetElemNumber>=20&&lDetElemNumber<=25)){ + if ((lDetElemNumber>=0&&lDetElemNumber<=6) || + (lDetElemNumber>=20&&lDetElemNumber<=25)){ lConstraintR[lDetElem*fgNParCh+iVar]=lWeight; } } } +//______________________________________________________________________ void AliMUONAlignment::ResetConstraints(){ /// Reset all constraint equations - for (Int_t i = 0; i < fgNDetElem; i++){ + for (Int_t i = 0; i < fgNDetElem; i++){ fConstraintX[i*fgNParCh+0]=0.0; fConstraintX[i*fgNParCh+1]=0.0; fConstraintX[i*fgNParCh+2]=0.0; @@ -672,25 +845,29 @@ void AliMUONAlignment::ResetConstraints(){ } } +//______________________________________________________________________ void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) { /// Constrain equation defined by par to value fMillepede->SetGlobalConstraint(par, value); AliInfo("Adding constraint"); } +//______________________________________________________________________ void AliMUONAlignment::InitGlobalParameters(Double_t *par) { /// Initialize global parameters with par array fMillepede->SetGlobalParameters(par); AliInfo("Init Global Parameters"); } - + +//______________________________________________________________________ void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) { - /// Parameter iPar is encourage to vary in [-value;value]. + /// Parameter iPar is encourage to vary in [-value;value]. /// If value == 0, parameter is fixed fMillepede->SetParSigma(iPar, value); - if (value==0) AliInfo(Form("Parameter %i Fixed", iPar)); + if (TMath::Abs(value)<1e-4) AliInfo(Form("Parameter %i Fixed", iPar)); } +//______________________________________________________________________ void AliMUONAlignment::ResetLocalEquation() { /// Reset the derivative vectors @@ -702,232 +879,375 @@ void AliMUONAlignment::ResetLocalEquation() } } -void AliMUONAlignment::AllowVariations(Bool_t *bStOnOff) { - /// Set allowed variation for selected stations based on fDoF and fAllowVar - for (Int_t iSt=1; iSt<=5; iSt++) { - if (bStOnOff[iSt-1]) { - Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0; - Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1]; - for (int i=0; i1) ? fgSNDetElemCh[iCh-2] : 0; + Int_t iDetElemLast = fgSNDetElemCh[iCh-1]; + for (int i=0; iSetNonLinear(iPar); AliInfo(Form("Parameter %i set to non linear", iPar)); } -void AliMUONAlignment::LocalEquationX() { - /// Define local equation for current track and hit in x coor. measurement - // set local derivatives - SetLocalDerivative(0, fCosPhi); - SetLocalDerivative(1, fCosPhi * (fTrackPos[2] - fTrackPos0[2])); - SetLocalDerivative(2, fSinPhi); - SetLocalDerivative(3, fSinPhi * (fTrackPos[2] - fTrackPos0[2])); +//______________________________________________________________________ +void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY) +{ + + /// Set expected measurement resolution + fSigma[0] = sigmaX; fSigma[1] = sigmaY; + AliInfo(Form("Using fSigma[0]=%f and fSigma[1]=%f",fSigma[0],fSigma[1])); +} + + +//______________________________________________________________________ +void AliMUONAlignment::LocalEquationX( Bool_t doAlignment ) +{ + + + // local cluster record + AliMUONAlignmentClusterRecord clusterRecord; + + // store detector and measurement + clusterRecord.SetDetElemId( fDetElemId ); + clusterRecord.SetDetElemNumber( fDetElemNumber ); + clusterRecord.SetMeas( fMeas[0] ); + clusterRecord.SetSigma( fSigma[0] ); + + // store local derivatives + clusterRecord.SetLocalDerivative( 0, fCosPhi ); + clusterRecord.SetLocalDerivative( 1, fCosPhi*(fTrackPos[2] - fTrackPos0[2]) ); + clusterRecord.SetLocalDerivative( 2, fSinPhi ); + clusterRecord.SetLocalDerivative( 3, fSinPhi*(fTrackPos[2] - fTrackPos0[2]) ); + + // store global derivatives + clusterRecord.SetGlobalDerivative( 0, -fCosPhi ); + clusterRecord.SetGlobalDerivative( 1, -fSinPhi ); + + if (fBFieldOn) + { + + clusterRecord.SetGlobalDerivative( + 2, + -fSinPhi*(fTrackPos[0]-fDetElemPos[0]) + +fCosPhi*(fTrackPos[1]-fDetElemPos[1]) ); + + } else { + + clusterRecord.SetGlobalDerivative( + 2, + -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0]) + +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]) ); - // set global derivatives - SetGlobalDerivative(fDetElemNumber*fgNParCh+0, -fCosPhi); - SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -fSinPhi); - if (fBFieldOn){ - SetGlobalDerivative(fDetElemNumber*fgNParCh+2, - -fSinPhi*(fTrackPos[0]-fDetElemPos[0]) - +fCosPhi*(fTrackPos[1]-fDetElemPos[1])); - } - else { - SetGlobalDerivative(fDetElemNumber*fgNParCh+2, - -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]* - (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0]) - +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]* - (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1])); } - fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]); + clusterRecord.SetGlobalDerivative( 3, fCosPhi*fTrackSlope0[0]+fSinPhi*fTrackSlope0[1] ); + + // append to trackRecord + fTrackRecord.AddClusterRecord( clusterRecord ); + + // store local equation + if( doAlignment ) LocalEquation( clusterRecord ); + } -void AliMUONAlignment::LocalEquationY() { - /// Define local equation for current track and hit in y coor. measurement - // set local derivatives - SetLocalDerivative(0,-fSinPhi); - SetLocalDerivative(1,-fSinPhi * (fTrackPos[2] - fTrackPos0[2])); - SetLocalDerivative(2, fCosPhi); - SetLocalDerivative(3, fCosPhi * (fTrackPos[2] - fTrackPos0[2])); +//______________________________________________________________________ +void AliMUONAlignment::LocalEquationY(Bool_t doAlignment ) +{ + + // local cluster record + AliMUONAlignmentClusterRecord clusterRecord; + + // store detector and measurement + clusterRecord.SetDetElemId( fDetElemId ); + clusterRecord.SetDetElemNumber( fDetElemNumber ); + clusterRecord.SetMeas( fMeas[1] ); + clusterRecord.SetSigma( fSigma[1] ); + + // store local derivatives + clusterRecord.SetLocalDerivative( 0, -fSinPhi ); + clusterRecord.SetLocalDerivative( 1, -fSinPhi*(fTrackPos[2] - fTrackPos0[2] ) ); + clusterRecord.SetLocalDerivative( 2, fCosPhi ); + clusterRecord.SetLocalDerivative( 3, fCosPhi*(fTrackPos[2] - fTrackPos0[2] ) ); // set global derivatives - SetGlobalDerivative(fDetElemNumber*fgNParCh+0, fSinPhi); - SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -fCosPhi); - if (fBFieldOn){ - SetGlobalDerivative(fDetElemNumber*fgNParCh+2, - -fCosPhi*(fTrackPos[0]-fDetElemPos[0]) - -fSinPhi*(fTrackPos[1]-fDetElemPos[1])); + clusterRecord.SetGlobalDerivative( 0, fSinPhi); + clusterRecord.SetGlobalDerivative( 1, -fCosPhi); + + if (fBFieldOn) + { + + clusterRecord.SetGlobalDerivative( + 2, + -fCosPhi*(fTrackPos[0]-fDetElemPos[0]) + -fSinPhi*(fTrackPos[1]-fDetElemPos[1])); + + } else { + + clusterRecord.SetGlobalDerivative( + 2, + -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0]) + -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*(fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1])); } - else { - SetGlobalDerivative(fDetElemNumber*fgNParCh+2, - -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]* - (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0]) - -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]* - (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1])); + + clusterRecord.SetGlobalDerivative( 3, -fSinPhi*fTrackSlope0[0]+fCosPhi*fTrackSlope0[1]); + + // append to trackRecord + fTrackRecord.AddClusterRecord( clusterRecord ); + + // store local equation + if( doAlignment ) LocalEquation( clusterRecord ); + +} + +//_____________________________________________________ +void AliMUONAlignment::LocalEquation( const AliMUONAlignmentClusterRecord& clusterRecord ) +{ + + // copy to local variables + for( Int_t index = 0; index < 4; ++index ) + { + SetLocalDerivative( index, clusterRecord.GetLocalDerivative( index ) ); + SetGlobalDerivative( clusterRecord.GetDetElemNumber()*fgNParCh + index, clusterRecord.GetGlobalDerivative( index ) ); } - fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]); + // pass equation parameters to millepede + fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, clusterRecord.GetMeas(), clusterRecord.GetSigma() ); + } -void AliMUONAlignment::FillRecPointData() { - /// Get information of current hit - fClustPos[0] = fRecHit->GetNonBendingCoor(); - fClustPos[1] = fRecHit->GetBendingCoor(); - fClustPos[2] = fRecHit->GetZ(); - fTransform->Global2Local(fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2], - fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]); +//______________________________________________________________________ +void AliMUONAlignment::FillRecPointData() +{ + + /// Get information of current cluster + fClustPos[0] = fCluster->GetX(); + fClustPos[1] = fCluster->GetY(); + fClustPos[2] = fCluster->GetZ(); + fTransform->Global2Local( + fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2], + fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]); + } -void AliMUONAlignment::FillTrackParamData() { - /// Get information of current track at current hit +//______________________________________________________________________ +void AliMUONAlignment::FillTrackParamData() +{ + + /// Get information of current track at current cluster fTrackPos[0] = fTrackParam->GetNonBendingCoor(); fTrackPos[1] = fTrackParam->GetBendingCoor(); fTrackPos[2] = fTrackParam->GetZ(); fTrackSlope[0] = fTrackParam->GetNonBendingSlope(); fTrackSlope[1] = fTrackParam->GetBendingSlope(); - fTransform->Global2Local(fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2], - fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]); + fTransform->Global2Local( + fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2], + fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]); + } -void AliMUONAlignment::FillDetElemData() { +//_____________________________________________________ +void AliMUONAlignment::FillDetElemData() +{ + /// Get information of current detection element Double_t lDetElemLocX = 0.; Double_t lDetElemLocY = 0.; Double_t lDetElemLocZ = 0.; - fDetElemId = fRecHit->GetDetElemId(); + fDetElemId = fCluster->GetDetElemId(); fDetElemNumber = fDetElemId%100; - for (int iCh=0; iChLocal2Global(fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ, - fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]); + for (int iCh=0; iChLocal2Global( + fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ, + fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]); + } -void AliMUONAlignment::ProcessTrack(AliMUONTrack * track) { - /// Process track; Loop over hits and set local equations +//_____________________________________________________ +AliMUONAlignmentTrackRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment ) +{ + + // store current track in running member. fTrack = track; - // get tclones arrays. - fTrackParamAtHit = fTrack->GetTrackParamAtHit(); - fHitForRecAtHit = fTrack->GetHitForRecAtHit(); - - // get size of arrays - Int_t nTrackParam = fTrackParamAtHit->GetEntries(); - Int_t nHitForRec = fHitForRecAtHit->GetEntries(); + + // clear track record + fTrackRecord.Clear(); + + // get number of tracks + Int_t nTrackParam = fTrack->GetTrackParamAtCluster()->GetEntries(); AliDebug(1,Form("Number of track param entries : %i ", nTrackParam)); - AliDebug(1,Form("Number of hit for rec entries : %i ", nHitForRec)); - for(Int_t iHit=0; iHitAt(iHit); - fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit); - if (!fRecHit || !fTrackParam) continue; - // fill local variables for this position --> one measurement - FillDetElemData(); - FillRecPointData(); - FillTrackParamData(); -// if (fDetElemId<500) continue; - fTrackPos0[0] = fTrackPos[0]; - fTrackPos0[1] = fTrackPos[1]; - fTrackPos0[2] = fTrackPos[2]; - fTrackSlope0[0] = fTrackSlope[0]; - fTrackSlope0[1] = fTrackSlope[1]; - break; - } + Bool_t first( kTRUE ); + for(Int_t iCluster=0; iClusterAt(iHit); - fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit); - if (!fRecHit || !fTrackParam) continue; + fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtCluster()->At(iCluster); + if ( ! fTrackParam ) continue; + fCluster = fTrackParam->GetClusterPtr(); + if ( ! fCluster ) continue; + // if (fDetElemId<500) continue; + // fill local variables for this position --> one measurement - FillDetElemData(); + FillDetElemData(); FillRecPointData(); FillTrackParamData(); -// if (fDetElemId<500) continue; - AliDebug(1,Form("cluster: %i", iHit)); - AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId)); - AliDebug(1,Form("fDetElemPos[0]: %f\t fDetElemPos[1]: %f\t fDetElemPos[2]: %f\t DetElemID: %i\t ", fDetElemPos[0],fDetElemPos[1],fDetElemPos[2], fDetElemId)); - AliDebug(1,Form("Track Parameter: %i", iHit)); - AliDebug(1,Form("x: %f\t y: %f\t z: %f\t slopex: %f\t slopey: %f", fTrackPos[0], fTrackPos[1], fTrackPos[2], fTrackSlope[0], fTrackSlope[1])); - AliDebug(1,Form("x0: %f\t y0: %f\t z0: %f\t slopex0: %f\t slopey0: %f", fTrackPos0[0], fTrackPos0[1], fTrackPos0[2], fTrackSlope0[0], fTrackSlope0[1])); - + if( first ) + { + + // for first valid cluster, save track position as "starting" values + first = kFALSE; + + fTrackPos0[0] = fTrackPos[0]; + fTrackPos0[1] = fTrackPos[1]; + fTrackPos0[2] = fTrackPos[2]; + fTrackSlope0[0] = fTrackSlope[0]; + fTrackSlope0[1] = fTrackSlope[1]; + + } + + // calculate measurements fCosPhi = TMath::Cos(fPhi); fSinPhi = TMath::Sin(fPhi); - if (fBFieldOn){ + if( fBFieldOn ) + { + fMeas[0] = fTrackPos[0] - fClustPos[0]; fMeas[1] = fTrackPos[1] - fClustPos[1]; - } - else { + + } else { + fMeas[0] = - fClustPos[0]; fMeas[1] = - fClustPos[1]; + } - AliDebug(1,Form("fMeas[0]: %f\t fMeas[1]: %f\t fSigma[0]: %f\t fSigma[1]: %f", fMeas[0], fMeas[1], fSigma[0], fSigma[1])); + + // soùe debugging + AliDebug(1,Form("cluster: %i", iCluster)); + AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId)); + AliDebug(1,Form("fDetElemPos[0]: %f\t fDetElemPos[1]: %f\t fDetElemPos[2]: %f\t DetElemID: %i\t ", fDetElemPos[0],fDetElemPos[1],fDetElemPos[2], fDetElemId)); + + AliDebug(1,Form("Track Parameter: %i", iCluster)); + AliDebug(1,Form("x: %f\t y: %f\t z: %f\t slopex: %f\t slopey: %f", fTrackPos[0], fTrackPos[1], fTrackPos[2], fTrackSlope[0], fTrackSlope[1])); + AliDebug(1,Form("x0: %f\t y0: %f\t z0: %f\t slopex0: %f\t slopey0: %f", fTrackPos0[0], fTrackPos0[1], fTrackPos0[2], fTrackSlope0[0], fTrackSlope0[1])); + + AliDebug(1,Form("fMeas[0]: %f\t fMeas[1]: %f\t fSigma[0]: %f\t fSigma[1]: %f", fMeas[0], fMeas[1], fSigma[0], fSigma[1])); + // Set local equations - LocalEquationX(); - LocalEquationY(); + LocalEquationX( doAlignment ); + LocalEquationY( doAlignment ); + } + + return &fTrackRecord; + } -void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) { +//______________________________________________________________________________ +void AliMUONAlignment::ProcessTrack( AliMUONAlignmentTrackRecord* track, Bool_t doAlignment ) +{ + + if( !( track && doAlignment ) ) return; + + // loop over clusters + for( Int_t index = 0; index < track->GetNRecords(); ++index ) + { if( AliMUONAlignmentClusterRecord* record = track->GetRecord( index ) ) LocalEquation( *record ); } + + return; + +} + +//_____________________________________________________ +void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) +{ + /// Call local fit for this tracks Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit); - if (iRes && !lSingleFit) { - fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1); - } + if (iRes && !lSingleFit) + { fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1); } + } -void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) { +//_____________________________________________________ +void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) +{ + /// Call global fit; Global parameters are stored in parameters fMillepede->GlobalFit(parameters,errors,pulls); AliInfo("Done fitting global parameters!"); - for (int iGlob=0; iGlobGetParError(iPar); return lErr; } -void AliMUONAlignment::PrintGlobalParameters() { +//_____________________________________________________ +void AliMUONAlignment::PrintGlobalParameters() +{ /// Print global parameters fMillepede->PrintGlobalParameters(); } //_________________________________________________________________________ -TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, double *lMisAlignment) const +TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, const double *lMisAlignment) const { /// Realign given transformation by given misalignment and return the misaligned transformation - + Double_t cartMisAlig[3] = {0,0,0}; Double_t angMisAlig[3] = {0,0,0}; // const Double_t *trans = transform.GetTranslation(); // TGeoRotation *rot; // // check if the rotation we obtain is not NULL -// if (transform.GetRotation()) { +// if (transform.GetRotation()) { // rot = transform.GetRotation(); // } -// else { +// else { // rot = new TGeoRotation("rot"); // } // default constructor. - cartMisAlig[0] = -lMisAlignment[0]; - cartMisAlig[1] = -lMisAlignment[1]; - angMisAlig[2] = -lMisAlignment[2]*180./TMath::Pi(); + cartMisAlig[0] = -TMath::Sign(1.0,transform.GetRotationMatrix()[0])*lMisAlignment[0]; + cartMisAlig[1] = -TMath::Sign(1.0,transform.GetRotationMatrix()[4])*lMisAlignment[1]; + cartMisAlig[2] = -TMath::Sign(1.0,transform.GetRotationMatrix()[8])*lMisAlignment[3]; + angMisAlig[2] = -TMath::Sign(1.0,transform.GetRotationMatrix()[0]*transform.GetRotationMatrix()[4])*lMisAlignment[2]*180./TMath::Pi(); TGeoTranslation deltaTrans(cartMisAlig[0], cartMisAlig[1], cartMisAlig[2]); TGeoRotation deltaRot; @@ -944,113 +1264,114 @@ TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, doubl //______________________________________________________________________ AliMUONGeometryTransformer * AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer, - double *misAlignments, Bool_t verbose) - + const double *misAlignments, Bool_t verbose) + { - ///////////////////////////////////////////////////////////////////// - // Takes the internal geometry module transformers, copies them + /// Returns a new AliMUONGeometryTransformer with the found misalignments + /// applied. + + // Takes the internal geometry module transformers, copies them // and gets the Detection Elements from them. // Takes misalignment parameters and applies these // to the local transform of the Detection Element // Obtains the global transform by multiplying the module transformer - // transformation with the local transformation + // transformation with the local transformation // Applies the global transform to a new detection element // Adds the new detection element to a new module transformer // Adds the new module transformer to a new geometry transformer // Returns the new geometry transformer - Double_t lModuleMisAlignment[3] = {0.,0.,0.}; - Double_t lDetElemMisAlignment[3] = {0.,0.,0.}; + Double_t lModuleMisAlignment[4] = {0.,0.,0.,0.}; + Double_t lDetElemMisAlignment[4] = {0.,0.,0.,0.}; Int_t iDetElemId = 0; Int_t iDetElemNumber = 0; AliMUONGeometryTransformer *newGeometryTransformer = - new AliMUONGeometryTransformer(kTRUE); + new AliMUONGeometryTransformer(); for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) { - // module transformers + // module transformers const AliMUONGeometryModuleTransformer *kModuleTransformer = transformer->GetModuleTransformer(iMt, true); - + AliMUONGeometryModuleTransformer *newModuleTransformer = new AliMUONGeometryModuleTransformer(iMt); newGeometryTransformer->AddModuleTransformer(newModuleTransformer); - + TGeoCombiTrans moduleTransform = TGeoCombiTrans(*kModuleTransformer->GetTransformation()); // New module transformation TGeoCombiTrans newModuleTransform = ReAlign(moduleTransform,lModuleMisAlignment); newModuleTransformer->SetTransformation(newModuleTransform); - - // Get delta transformation: + + // Get delta transformation: // Tdelta = Tnew * Told.inverse - TGeoHMatrix deltaModuleTransform = - AliMUONGeometryBuilder::Multiply(newModuleTransform, - kModuleTransformer->GetTransformation()->Inverse()); + TGeoHMatrix deltaModuleTransform = + AliMUONGeometryBuilder::Multiply(newModuleTransform, + kModuleTransformer->GetTransformation()->Inverse()); // Create module mis alignment matrix newGeometryTransformer ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform); - + AliMpExMap *detElements = kModuleTransformer->GetDetElementStore(); - + if (verbose) AliInfo(Form("%i DEs in old GeometryStore %i",detElements->GetSize(), iMt)); - for (Int_t iDe = 0; iDe < detElements->GetSize(); iDe++) { - // detection elements. - AliMUONGeometryDetElement *detElement = - (AliMUONGeometryDetElement *) detElements->GetObject(iDe); - if (!detElement) - AliFatal("Detection element not found."); - - /// make a new detection element + TIter next(detElements->CreateIterator()); + AliMUONGeometryDetElement* detElement; + Int_t iDe(-1); + while ( ( detElement = static_cast(next()) ) ) + { + ++iDe; + // make a new detection element AliMUONGeometryDetElement *newDetElement = - new AliMUONGeometryDetElement(detElement->GetId(), - detElement->GetVolumePath()); + new AliMUONGeometryDetElement(detElement->GetId(), + detElement->GetVolumePath()); TString lDetElemName(detElement->GetDEName()); lDetElemName.ReplaceAll("DE",""); iDetElemId = lDetElemName.Atoi(); iDetElemNumber = iDetElemId%100; for (int iCh=0; iChGetLocalTransformation()); + = TGeoCombiTrans(*detElement->GetLocalTransformation()); TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment); - newDetElement->SetLocalTransformation(newLocalTransform); + newDetElement->SetLocalTransformation(newLocalTransform); // global transformation TGeoHMatrix newGlobalTransform = - AliMUONGeometryBuilder::Multiply(newModuleTransform, - newLocalTransform); + AliMUONGeometryBuilder::Multiply(newModuleTransform, + newLocalTransform); newDetElement->SetGlobalTransformation(newGlobalTransform); - + // add this det element to module newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), - newDetElement); + newDetElement); // In the Alice Alignment Framework misalignment objects store // global delta transformation // Get detection "intermediate" global transformation TGeoHMatrix newOldGlobalTransform = newModuleTransform * localTransform; - // Get detection element global delta transformation: + // Get detection element global delta transformation: // Tdelta = Tnew * Told.inverse TGeoHMatrix deltaGlobalTransform - = AliMUONGeometryBuilder::Multiply(newGlobalTransform, - newOldGlobalTransform.Inverse()); - + = AliMUONGeometryBuilder::Multiply(newGlobalTransform, + newOldGlobalTransform.Inverse()); + // Create mis alignment matrix newGeometryTransformer - ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform); + ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform); } - + if (verbose) AliInfo(Form("Added module transformer %i to the transformer", iMt)); newGeometryTransformer->AddModuleTransformer(newModuleTransformer); @@ -1058,3 +1379,59 @@ AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer, return newGeometryTransformer; } +//______________________________________________________________________ +void AliMUONAlignment::SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t rChId, Double_t rChResX, Double_t rChResY, Double_t rDeResX, Double_t rDeResY) +{ + + /// Set alignment resolution to misalign objects to be stored in CDB + Int_t chIdMin = (rChId<0)? 0 : rChId; + Int_t chIdMax = (rChId<0)? 9 : rChId; + Double_t chResX = rChResX; + Double_t chResY = rChResY; + Double_t deResX = rDeResX; + Double_t deResY = rDeResY; + + TMatrixDSym mChCorrMatrix(6); + mChCorrMatrix[0][0]=chResX*chResX; + mChCorrMatrix[1][1]=chResY*chResY; + // mChCorrMatrix.Print(); + + TMatrixDSym mDECorrMatrix(6); + mDECorrMatrix[0][0]=deResX*deResX; + mDECorrMatrix[1][1]=deResY*deResY; + // mDECorrMatrix.Print(); + + AliAlignObjMatrix *alignMat = 0x0; + + for(Int_t chId=chIdMin; chId<=chIdMax; chId++) { + TString chName1; + TString chName2; + if (chId<4){ + chName1 = Form("GM%d",chId); + chName2 = Form("GM%d",chId); + } else { + chName1 = Form("GM%d",4+(chId-4)*2); + chName2 = Form("GM%d",4+(chId-4)*2+1); + } + + for (int i=0; iGetEntries(); i++) { + alignMat = (AliAlignObjMatrix*)misAlignArray->At(i); + TString volName(alignMat->GetSymName()); + if((volName.Contains(chName1)&& + ((volName.Last('/')==volName.Index(chName1)+chName1.Length())|| + (volName.Length()==volName.Index(chName1)+chName1.Length())))|| + (volName.Contains(chName2)&& + ((volName.Last('/')==volName.Index(chName2)+chName2.Length())|| + (volName.Length()==volName.Index(chName2)+chName2.Length())))){ + volName.Remove(0,volName.Last('/')+1); + if (volName.Contains("GM")) { + // alignMat->Print("NULL"); + alignMat->SetCorrMatrix(mChCorrMatrix); + } else if (volName.Contains("DE")) { + // alignMat->Print("NULL"); + alignMat->SetCorrMatrix(mDECorrMatrix); + } + } + } + } +}