* provided "as is" without express or implied warranty. *
**************************************************************************/
-/* $Id$ */
+/* $Id: AliMUONAlignment.cxx 51000 2011-08-08 17:58:17Z ivana $ */
//-----------------------------------------------------------------------------
/// \class AliMUONAlignment
#include "AliMUONGeometryModuleTransformer.h"
#include "AliMUONGeometryDetElement.h"
#include "AliMUONGeometryBuilder.h"
-#include "AliMillepede.h"
+#include "AliMillePede2.h"
#include "AliMpExMap.h"
#include "AliMpExMapIterator.h"
#include "AliAlignObjMatrix.h"
#include "AliLog.h"
-#include "TMath.h"
-#include "TMatrixDSym.h"
-#include "TClonesArray.h"
+#include <TMath.h>
+#include <TMatrixDSym.h>
+#include <TClonesArray.h>
+#include <TGraphErrors.h>
/// \cond CLASSIMP
ClassImp(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;
+const Int_t AliMUONAlignment::fgNDetElemCh[AliMUONAlignment::fgNCh] = { 4, 4, 4, 4, 18, 18, 26, 26, 26, 26 };
+const Int_t AliMUONAlignment::fgSNDetElemCh[AliMUONAlignment::fgNCh+1] = { 0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156 };
+
+//_____________________________________________________________________
+/// self initialized array, used for adding constraints
+class Array
+{
+
+ public:
+
+ /// contructor
+ Array( void )
+ {
+ for( Int_t i=0; i < AliMUONAlignment::fNGlobal; ++i )
+ { values[i] = 0; }
+ }
+
+ /// array
+ Double_t values[AliMUONAlignment::fNGlobal];
+
+ private:
+
+ /// Not implemented
+ Array(const Array& );
+
+ /// Not implemented
+ Array& operator = (const Array& );
+
+};
//_____________________________________________________________________
AliMUONAlignment::AliMUONAlignment()
: TObject(),
- fBFieldOn(kTRUE),
- fStartFac(256.),
- fResCutInitial(100.),
- fResCut(100.),
- fMillepede(0),
- fTrack(0),
- fCluster(0),
- fTrackParam(0),
- fNGlobal(fgNDetElem*fgNParCh),
- fNLocal(4),
- fNStdDev(3),
- fDetElemId(0),
- fDetElemNumber(0),
- fPhi(0.),
- fCosPhi(1.),
- fSinPhi(0.),
+ fInitialized( kFALSE ),
+ fRunNumber( 0 ),
+ fBFieldOn( kTRUE ),
+ fStartFac( 256 ),
+ fResCutInitial( 100 ),
+ fResCut( 100 ),
+ fMillepede( 0L ),
+ fCluster( 0L ),
+ fNStdDev( 3 ),
+ fDetElemNumber( 0 ),
+ fUnbias( kFALSE ),
fTrackRecord(),
- fTransform(0)
+ fTransform( 0 ),
+ fGeoCombiTransInverse()
{
- /// Default constructor
+ /// constructor
fSigma[0] = 1.5e-1;
fSigma[1] = 1.0e-2;
- AliInfo(Form("fSigma[0]: %f\t fSigma[1]: %f",fSigma[0],fSigma[1]));
+ // default allowed variations
+ fAllowVar[0] = 0.5; // x
+ fAllowVar[1] = 0.5; // y
+ fAllowVar[2] = 0.01; // phi_z
+ fAllowVar[3] = 5; // z
- 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;
+ // initialize millepede
+ fMillepede = new AliMillePede2();
- 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));
+ // initialize degrees of freedom
+ // by default all parameters are free
+ for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
+ { fGlobalParameterStatus[iPar] = kFreeParId; }
- fMillepede = new AliMillepede();
+ // initialize local equations
+ for(int i=0; i<fNLocal; ++i )
+ { fLocalDerivatives[i] = 0.0; }
- Init(fNGlobal, fNLocal, fNStdDev);
-
- ResetLocalEquation();
- AliInfo("Parameters initialized to zero");
+ for(int i=0; i<fNGlobal; ++i )
+ { fGlobalDerivatives[i] = 0.0; }
}
//_____________________________________________________________________
-AliMUONAlignment::AliMUONAlignment(TRootIOCtor* /*dummy*/)
-: TObject(),
-fBFieldOn(kFALSE),
-fStartFac(0.),
-fResCutInitial(0.),
-fResCut(0.),
-fMillepede(0),
-fTrack(0),
-fCluster(0),
-fTrackParam(0),
-fNGlobal(0),
-fNLocal(0),
-fNStdDev(0),
-fDetElemId(0),
-fDetElemNumber(0),
-fPhi(0.),
-fCosPhi(0.),
-fSinPhi(0.),
-fTrackRecord(),
-fTransform(0)
+AliMUONAlignment::~AliMUONAlignment()
{
- /// Root IO constructor
- ResetConstraints();
- for (Int_t iCh=0; iCh<10; iCh++) {
- fChOnOff[iCh] = kFALSE;
- }
- fSpecLROnOff[0] = kFALSE; fSpecLROnOff[1] = kFALSE;
- for (Int_t iDoF=0; iDoF<4; iDoF++) {
- fDoF[iDoF] = kFALSE;
- fAllowVar[iDoF] = kFALSE;
- }
- for (Int_t i=0; i<3; i++) {
- fClustPos[i] = 0;
- fClustPosLoc[i] = 0;
- fTrackPos0[i] = 0;
- fTrackPos[i] = 0;
- fTrackPosLoc[i] = 0;
- fDetElemPos[i] = 0;
- }
- fTrackSlope0[0]=0; fTrackSlope0[1]=0;
- fTrackSlope[0]=0; fTrackSlope[1]=0;
- fMeas[0]=0; fMeas[1]=0;
- fSigma[0]=0; fSigma[1]=0;
- for (Int_t iLPar=0; iLPar<4; iLPar++) {
- fLocalDerivatives[iLPar]=0;
- }
- for (Int_t iGPar=0; iGPar<624; iGPar++) {
- fGlobalDerivatives[iGPar]=0;
- }
+ /// 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( void )
{
- /// 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 bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
- // Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
+ /// initialize
+ /**
+ initialize millipede
+ must be called after necessary detectors have been fixed,
+ but before constrains are added and before global parameters initial value are set
+ */
+ if( fInitialized )
+ { AliFatal( "Millepede already initialized" ); }
+
+ // assign proper groupID to free parameters
+ Int_t nGlobal = 0;
+ for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
+ {
- // AllowVariations(bChOnOff);
+ if( fGlobalParameterStatus[iPar] == kFixedParId )
+ {
+ // fixed parameters are left unchanged
+ continue;
- // Fix parameters or add constraints here
- // for (Int_t iSt=0; iSt<5; iSt++)
- // { if (!bStOnOff[iSt]) FixStation(iSt+1); }
+ } else if( fGlobalParameterStatus[iPar] == kFreeParId || fGlobalParameterStatus[iPar] == kGroupBaseId ) {
- // for (Int_t iCh=0; iCh<10; iCh++)
- // { if (!bChOnOff[iCh]) FixChamber(iCh+1); }
+ // free parameters or first element of group are assigned a new group id
+ fGlobalParameterStatus[iPar] = nGlobal++;
+ continue;
-// FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
+ } else if( fGlobalParameterStatus[iPar] < kGroupBaseId ) {
- ResetConstraints();
+ // get detector element id from status, get chamber parameter id
+ const Int_t iDeBase( kGroupBaseId - 1 - fGlobalParameterStatus[iPar] );
+ const Int_t iParBase = iPar%fgNParCh;
- // 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};
- Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
- // AddConstraints(bStOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
+ // check
+ if( iDeBase < 0 || iDeBase >= iPar/fgNParCh )
+ { AliFatal( Form( "Group for parameter index %i has wrong base detector element: %i", iPar, iDeBase ) ); }
- // 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);
+ // assign identical group id to current
+ fGlobalParameterStatus[iPar] = fGlobalParameterStatus[iDeBase*fgNParCh + iParBase];
+ AliInfo( Form( "Parameter %i grouped to detector %i (%s)", iPar, iDeBase, GetParameterMaskString( 1<<iParBase ).Data() ) );
- bVarXYT[0] = kTRUE; bVarXYT[1] = kTRUE; bVarXYT[2] = kFALSE;
- // AddConstraints(bStOnOff,bVarXYT);
+ } else AliFatal( Form( "Unrecognized parameter status for index %i: %i", iPar, fGlobalParameterStatus[iPar] ) );
+
+ }
+
+ AliInfo( Form( "Free Parameters: %i out of %i", nGlobal, fNGlobal ) );
+
+ // initialize millepede
+ fMillepede->InitMille( fNGlobal, fNLocal, fNStdDev, fResCut, fResCutInitial, fGlobalParameterStatus );
+ fInitialized = kTRUE;
+
+ // some debug output
+ for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
+ { AliInfo( Form( "fAllowVar[%i]= %f", iPar, fAllowVar[iPar] ) ); }
+
+ // set allowed variations for all parameters
+ for( Int_t iDet = 0; iDet < fgNDetElem; ++iDet )
+ {
+ for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
+ { fMillepede->SetParSigma( iDet*fgNParCh + iPar, fAllowVar[iPar] ); }
+ }
// Set iterations
if (fStartFac>1) fMillepede->SetIterations(fStartFac);
+
}
-//_____________________________________________________________________
-void AliMUONAlignment::FixStation(Int_t iSt)
+//_____________________________________________________
+AliMillePedeRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment, Double_t weight )
{
- /// 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++)
+ /// process track for alignment minimization
+ /**
+ returns the alignment records for this track.
+ They can be stored in some output for later reprocessing.
+ */
+
+ // reset track records
+ fTrackRecord.Reset();
+ if( fMillepede->GetRecord() ) fMillepede->GetRecord()->Reset();
+
+ // get number of track parameters
+ Int_t nTrackParam = track->GetTrackParamAtCluster()->GetEntries();
+
+ Bool_t first( kTRUE );
+ for( Int_t iTrackParam = 0; iTrackParam < nTrackParam; ++iTrackParam )
{
- FixParameter(i*fgNParCh+0, 0.0);
- FixParameter(i*fgNParCh+1, 0.0);
- FixParameter(i*fgNParCh+2, 0.0);
- FixParameter(i*fgNParCh+3, 0.0);
+
+ // get new pointers
+ AliMUONTrackParam* trackParam( (AliMUONTrackParam *) track->GetTrackParamAtCluster()->At(iTrackParam) );
+ if( !trackParam ) continue;
+
+ AliMUONVCluster* cluster = trackParam->GetClusterPtr();
+ if( !cluster ) continue;
+
+ // fill local variables for this position --> one measurement
+ FillDetElemData( cluster );
+ FillRecPointData( cluster );
+
+ // unbias and store track parameters
+ if( fUnbias && !UnbiasTrackParamData( trackParam ) ) continue;
+ FillTrackParamData( trackParam );
+
+ 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];
+
+ }
+
+ // 'inverse' (GlobalToLocal) rotation matrix
+ const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
+
+ // calculate measurements
+ if( fBFieldOn )
+ {
+
+ // use residuals (cluster - track) for measurement
+ fMeas[0] = r[0]*(fClustPos[0] - fTrackPos[0]) + r[1]*(fClustPos[1] - fTrackPos[1]);
+ fMeas[1] = r[3]*(fClustPos[0] - fTrackPos[0]) + r[4]*(fClustPos[1] - fTrackPos[1]);
+
+ } else {
+
+ // use cluster position for measurement
+ fMeas[0] = ( r[0]*fClustPos[0] + r[1]*fClustPos[1] );
+ fMeas[1] = ( r[3]*fClustPos[0] + r[4]*fClustPos[1] );
+
+ }
+
+ // Set local equations
+ LocalEquationX();
+ LocalEquationY();
+
+ }
+
+ // copy track record
+ fMillepede->SetRecordRun(fRunNumber);
+ fMillepede->SetRecordWeight(weight);
+ fTrackRecord = *fMillepede->GetRecord();
+
+ // save record data
+ if( doAlignment )
+ {
+ fMillepede->SaveRecordData();
}
+ // return record
+ return &fTrackRecord;
+
+}
+
+//______________________________________________________________________________
+void AliMUONAlignment::ProcessTrack( AliMillePedeRecord* trackRecord )
+{
+ /// process track record
+ if( !trackRecord ) return;
+
+ // make sure record storage is initialized
+ if( !fMillepede->GetRecord() ) fMillepede->InitDataRecStorage();
+
+ // copy content
+ *fMillepede->GetRecord() = *trackRecord;
+
+ // save record
+ fMillepede->SaveRecordData();
+
+ return;
+
}
//_____________________________________________________________________
-void AliMUONAlignment::FixChamber(Int_t iCh)
+void AliMUONAlignment::FixAll( UInt_t mask )
{
- /// 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++)
+ /// fix parameters matching mask, for all chambers
+ AliInfo( Form( "Fixing %s for all detector elements", GetParameterMaskString( mask ).Data() ) );
+
+ // fix all stations
+ for( Int_t i = 0; i < fgNDetElem; ++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);
+ if( mask & ParX ) FixParameter(i, 0);
+ if( mask & ParY ) FixParameter(i, 1);
+ if( mask & ParTZ ) FixParameter(i, 2);
+ if( mask & ParZ ) FixParameter(i, 3);
}
+
}
//_____________________________________________________________________
-void AliMUONAlignment::FixDetElem(Int_t iDetElemId, TString sVarXYT)
+void AliMUONAlignment::FixChamber( Int_t iCh, UInt_t mask )
{
+ /// fix parameters matching mask, for all detector elements in a given chamber, counting from 1
- /// Fix a given detection element
- Int_t iDetElemNumber = iDetElemId%100;
- for (int iCh=0; iCh<iDetElemId/100-1; iCh++)
- {
- iDetElemNumber += fgNDetElemCh[iCh];
- }
+ // check boundaries
+ if( iCh < 1 || iCh > 10 )
+ { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
- if (sVarXYT.Contains("X"))
+ // get first and last element
+ const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
+ const Int_t iDetElemLast = fgSNDetElemCh[iCh];
+ for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
{
- // X constraint
- FixParameter(iDetElemNumber*fgNParCh+0, 0.0);
- }
- if (sVarXYT.Contains("Y"))
- {
- // Y constraint
- FixParameter(iDetElemNumber*fgNParCh+1, 0.0);
- }
+ AliInfo( Form( "Fixing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
- if (sVarXYT.Contains("T"))
- {
- // T constraint
- FixParameter(iDetElemNumber*fgNParCh+2, 0.0);
- }
+ if( mask & ParX ) FixParameter(i, 0);
+ if( mask & ParY ) FixParameter(i, 1);
+ if( mask & ParTZ ) FixParameter(i, 2);
+ if( mask & ParZ ) FixParameter(i, 3);
- if (sVarXYT.Contains("Z"))
- {
- // T constraint
- FixParameter(iDetElemNumber*fgNParCh+3, 0.0);
}
+}
+
+//_____________________________________________________________________
+void AliMUONAlignment::FixDetElem( Int_t iDetElemId, UInt_t mask )
+{
+ /// fix parameters matching mask, for a given detector element, counting from 0
+ const Int_t iDet( GetDetElemNumber( iDetElemId ) );
+ if ( mask & ParX ) FixParameter(iDet, 0);
+ if ( mask & ParY ) FixParameter(iDet, 1);
+ if ( mask & ParTZ ) FixParameter(iDet, 2);
+ if ( mask & ParZ ) FixParameter(iDet, 3);
}
//_____________________________________________________________________
-void AliMUONAlignment::FixHalfSpectrometer(const Bool_t *lChOnOff, const Bool_t *lSpecLROnOff)
+void AliMUONAlignment::FixHalfSpectrometer( const Bool_t *lChOnOff, UInt_t sidesMask, UInt_t mask )
{
- /// Fix left or right detector
- for (Int_t i = 0; i < fgNDetElem; i++)
+ /// Fix parameters matching mask for all detectors in selected chambers and selected sides of the spectrometer
+ for( Int_t i = 0; i < fgNDetElem; ++i )
{
- Int_t iCh=0;
- for (iCh=1; iCh<=fgNCh; iCh++)
- { if (i<fgSNDetElemCh[iCh-1]) break; }
+ // get chamber matching detector
+ const Int_t iCh( GetChamberId(i) );
+ if( !lChOnOff[iCh-1] ) continue;
- if (lChOnOff[iCh-1])
+ // get detector element in chamber
+ Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
+
+ // skip detector if its side is off
+ // stations 1 and 2
+ if( iCh>=1 && iCh<=4 )
{
+ if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
+ if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
+ if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
+ if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
+ }
- Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
- if (iCh>=1 && iCh<=4)
- {
+ // station 3
+ if (iCh>=5 && iCh<=6)
+ {
+ if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
+ if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
+ if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
+ if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
+ }
- 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);
- }
+ // stations 4 and 5
+ if (iCh>=7 && iCh<=10)
+ {
+ if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
+ if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
+ if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
+ if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
+ }
- 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);
- }
+ // detector is accepted, fix it
+ FixDetElem( i, mask );
- }
+ }
- 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);
- }
+//______________________________________________________________________
+void AliMUONAlignment::FixParameter( Int_t iPar )
+{
- if (((lDetElemNumber>=0&&lDetElemNumber<=4) ||
- (lDetElemNumber>=14&&lDetElemNumber<=17)) && !lSpecLROnOff[1])
- {
+ /// fix a given parameter, counting from 0
+ if( fInitialized )
+ { AliFatal( "Millepede already initialized" ); }
- FixParameter(i*fgNParCh+0, 0.0);
- FixParameter(i*fgNParCh+1, 0.0);
- FixParameter(i*fgNParCh+2, 0.0);
- FixParameter(i*fgNParCh+3, 0.0);
- }
+ fGlobalParameterStatus[iPar] = kFixedParId;
- }
+}
- 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);
- }
+//_____________________________________________________________________
+void AliMUONAlignment::ReleaseChamber( Int_t iCh, UInt_t mask )
+{
+ /// release parameters matching mask, for all detector elements in a given chamber, counting from 1
- 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);
- }
+ // check boundaries
+ if( iCh < 1 || iCh > 10 )
+ { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
- }
+ // get first and last element
+ const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
+ const Int_t iDetElemLast = fgSNDetElemCh[iCh];
+ for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
+ {
- }
+ AliInfo( Form( "Releasing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
+
+ if( mask & ParX ) ReleaseParameter(i, 0);
+ if( mask & ParY ) ReleaseParameter(i, 1);
+ if( mask & ParTZ ) ReleaseParameter(i, 2);
+ if( mask & ParZ ) ReleaseParameter(i, 3);
}
+}
+
+//_____________________________________________________________________
+void AliMUONAlignment::ReleaseDetElem( Int_t iDetElemId, UInt_t mask )
+{
+ /// release parameters matching mask, for a given detector element, counting from 0
+ const Int_t iDet( GetDetElemNumber( iDetElemId ) );
+ if ( mask & ParX ) ReleaseParameter(iDet, 0);
+ if ( mask & ParY ) ReleaseParameter(iDet, 1);
+ if ( mask & ParTZ ) ReleaseParameter(iDet, 2);
+ if ( mask & ParZ ) ReleaseParameter(iDet, 3);
}
//______________________________________________________________________
-void AliMUONAlignment::SetNonLinear(const Bool_t *lChOnOff, const Bool_t *lVarXYT)
+void AliMUONAlignment::ReleaseParameter( Int_t iPar )
{
- /// Set non linear parameter flag selected chambers and degrees of freedom
- for (Int_t i = 0; i < fgNDetElem; i++)
- {
+ /// release a given parameter, counting from 0
+ if( fInitialized )
+ { AliFatal( "Millepede already initialized" ); }
- Int_t iCh=0;
- for (iCh=1; iCh<=fgNCh; iCh++)
- { if (i<fgSNDetElemCh[iCh-1]) break; }
+ fGlobalParameterStatus[iPar] = kFreeParId;
- if (lChOnOff[iCh-1])
- {
+}
- if (lVarXYT[0])
- {
- // X constraint
- SetNonLinear(i*fgNParCh+0);
- }
+//_____________________________________________________________________
+void AliMUONAlignment::GroupChamber( Int_t iCh, UInt_t mask )
+{
+ /// group parameters matching mask for all detector elements in a given chamber, counting from 1
+ if( iCh < 1 || iCh > 10 )
+ { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
- if (lVarXYT[1])
- {
- // Y constraint
- SetNonLinear(i*fgNParCh+1);
- }
+ const Int_t detElemMin = 100*iCh;
+ const Int_t detElemMax = 100*iCh + fgNDetElemCh[iCh]-1;
+ GroupDetElems( detElemMin, detElemMax, mask );
- if (lVarXYT[2])
- {
- // T constraint
- SetNonLinear(i*fgNParCh+2);
- }
+}
- if (lVarXYT[3])
- {
- // Z constraint
- SetNonLinear(i*fgNParCh+3);
- }
+//_____________________________________________________________________
+void AliMUONAlignment::GroupDetElems( Int_t detElemMin, Int_t detElemMax, UInt_t mask )
+{
+ /// group parameters matching mask for all detector elements between min and max
+ // check number of detector elements
+ const Int_t nDetElem = detElemMax - detElemMin + 1;
+ if( nDetElem<2 )
+ { AliFatal( Form( "Requested group of DEs %d-%d contains less than 2 DE's", detElemMin, detElemMax ) ); }
- }
+ // create list
+ Int_t* detElemList = new int[nDetElem];
+ for( Int_t i = 0; i < nDetElem; ++i )
+ { detElemList[i] = detElemMin+i; }
+
+ // group
+ GroupDetElems( detElemList, nDetElem, mask );
+ delete[] detElemList;
+
+}
+
+//_____________________________________________________________________
+void AliMUONAlignment::GroupDetElems( Int_t* detElemList, Int_t nDetElem, UInt_t mask )
+{
+ /// group parameters matching mask for all detector elements in list
+ if( fInitialized )
+ { AliFatal( "Millepede already initialized" ); }
+ const Int_t iDeBase( GetDetElemNumber( detElemList[0] ) );
+ for( Int_t i = 0; i < nDetElem; ++i )
+ {
+ const Int_t iDeCurrent( GetDetElemNumber( detElemList[i] ) );
+ if( mask & ParX ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 0] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
+ if( mask & ParY ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 1] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
+ if( mask & ParTZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 2] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
+ if( mask & ParZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 3] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
+
+ if( i== 0 ) AliInfo( Form( "Creating new group for detector %i and variable %s", detElemList[i], GetParameterMaskString( mask ).Data() ) );
+ else AliInfo( Form( "Adding detector element %i to current group", detElemList[i] ) );
}
}
//______________________________________________________________________
-void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT)
+void AliMUONAlignment::SetChamberNonLinear( Int_t iCh, UInt_t mask )
{
-
- /// Add constraint equations for selected chambers and degrees of freedom
- for (Int_t i = 0; i < fgNDetElem; i++)
+ /// Set parameters matching mask as non linear, for all detector elements in a given chamber, counting from 1
+ const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
+ const Int_t iDetElemLast = fgSNDetElemCh[iCh];
+ for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
{
- Int_t iCh=0;
- for (iCh=1; iCh<=fgNCh; iCh++)
- {
- if (i<fgSNDetElemCh[iCh-1]) break;
- }
+ if( mask & ParX ) SetParameterNonLinear(i, 0);
+ if( mask & ParY ) SetParameterNonLinear(i, 1);
+ if( mask & ParTZ ) SetParameterNonLinear(i, 2);
+ if( mask & ParZ ) SetParameterNonLinear(i, 3);
- if (lChOnOff[iCh-1])
- {
- if (lVarXYT[0])
- {
- // X constraint
- fConstraintX[i*fgNParCh+0]=1.0;
- }
+ }
- if (lVarXYT[1])
- {
- // Y constraint
- fConstraintY[i*fgNParCh+1]=1.0;
- }
+}
- if (lVarXYT[2])
- {
- // T constraint
- fConstraintP[i*fgNParCh+2]=1.0;
- }
- // if (lVarXYT[3]) { // Z constraint
- // fConstraintP[i*fgNParCh+3]=1.0;
- // }
+//_____________________________________________________________________
+void AliMUONAlignment::SetDetElemNonLinear( Int_t iDetElemId, UInt_t mask )
+{
+ /// Set parameters matching mask as non linear, for a given detector element, counting from 0
+ const Int_t iDet( GetDetElemNumber( iDetElemId ) );
+ if ( mask & ParX ) SetParameterNonLinear(iDet, 0);
+ if ( mask & ParY ) SetParameterNonLinear(iDet, 1);
+ if ( mask & ParTZ ) SetParameterNonLinear(iDet, 2);
+ if ( mask & ParZ ) SetParameterNonLinear(iDet, 3);
- }
- }
+}
- if (lVarXYT[0])
- {
- // X constraint
- AddConstraint(fConstraintX,0.0);
- }
+//______________________________________________________________________
+void AliMUONAlignment::SetParameterNonLinear( Int_t iPar )
+{
+ /// Set nonlinear flag for parameter iPar
+ if( !fInitialized )
+ { AliFatal( "Millepede not initialized" ); }
- if (lVarXYT[1])
- {
- // Y constraint
- AddConstraint(fConstraintY,0.0);
- }
+ fMillepede->SetNonLinear( iPar );
+ AliInfo( Form( "Parameter %i set to non linear", iPar ) );
+}
+
+//______________________________________________________________________
+void AliMUONAlignment::AddConstraints( const Bool_t *lChOnOff, UInt_t mask )
+{
+ /// Add constraint equations for selected chambers and degrees of freedom
+
+ Array fConstraintX;
+ Array fConstraintY;
+ Array fConstraintTZ;
+ Array fConstraintZ;
- if (lVarXYT[2])
+ for( Int_t i = 0; i < fgNDetElem; ++i )
{
- // T constraint
- AddConstraint(fConstraintP,0.0);
+
+ // get chamber matching detector
+ const Int_t iCh( GetChamberId(i) );
+ if (lChOnOff[iCh-1])
+ {
+
+ if( mask & ParX ) fConstraintX.values[i*fgNParCh+0]=1.0;
+ if( mask & ParY ) fConstraintY.values[i*fgNParCh+1]=1.0;
+ if( mask & ParTZ ) fConstraintTZ.values[i*fgNParCh+2]=1.0;
+ if( mask & ParZ ) fConstraintTZ.values[i*fgNParCh+3]=1.0;
+
+ }
}
-// if (lVarXYT[3]) { // Z constraint
-// AddConstraint(fConstraintP,0.0);
-// }
+ if( mask & ParX ) AddConstraint(fConstraintX.values,0.0);
+ if( mask & ParY ) AddConstraint(fConstraintY.values,0.0);
+ if( mask & ParTZ ) AddConstraint(fConstraintTZ.values,0.0);
+ if( mask & ParZ ) AddConstraint(fConstraintZ.values,0.0);
+
}
//______________________________________________________________________
-void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, const Bool_t *lDetTLBR, const Bool_t *lSpecLROnOff)
+void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, UInt_t sidesMask )
{
+ /*
+ questions:
+ - is there not redundancy/inconsistency between lDetTLBR and lSpecLROnOff ? shouldn't we use only lDetTLBR ?
+ - why is weight ignored for ConstrainT and ConstrainB
+ - why is there no constrain on z
+ */
+
/// Add constraint equations for selected chambers, degrees of freedom and detector half
- Double_t lDetElemLocX = 0.;
- Double_t lDetElemLocY = 0.;
- Double_t lDetElemLocZ = 0.;
- Double_t lDetElemGloX = 0.;
- Double_t lDetElemGloY = 0.;
- Double_t lDetElemGloZ = 0.;
Double_t lMeanY = 0.;
Double_t lSigmaY = 0.;
Double_t lMeanZ = 0.;
Double_t lSigmaZ = 0.;
Int_t lNDetElem = 0;
- for (Int_t i = 0; i < fgNDetElem; i++)
+
+ for( Int_t i = 0; i < fgNDetElem; ++i )
{
- Int_t iCh=0;
- for (iCh=1; iCh<=fgNCh; iCh++){
- if (i<fgSNDetElemCh[iCh-1]) break;
+ // get chamber matching detector
+ const Int_t iCh( GetChamberId(i) );
+
+ // skip detector if chamber is off
+ if( lChOnOff[iCh-1] ) continue;
+
+ // get detector element id from detector element number
+ const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
+ const Int_t lDetElemId = iCh*100+lDetElemNumber;
+
+ // skip detector if its side is off
+ // stations 1 and 2
+ if( iCh>=1 && iCh<=4 )
+ {
+ if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
+ if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
+ if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
+ if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
}
- if (lChOnOff[iCh-1]){
- Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
- Int_t lDetElemId = iCh*100+lDetElemNumber;
- fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
- 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 (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 (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]){
+
+ // station 3
+ if (iCh>=5 && iCh<=6)
+ {
+ if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
+ if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
+ if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
+ if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
+ }
+
+ // stations 4 and 5
+ if (iCh>=7 && iCh<=10)
+ {
+ if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
+ if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
+ if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
+ if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
+ }
+
+ // get global x, y and z position
+ Double_t lDetElemGloX = 0.;
+ Double_t lDetElemGloY = 0.;
+ Double_t lDetElemGloZ = 0.;
+ fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
+
+ // increment mean Y, mean Z, sigmas and number of accepted detectors
lMeanY += lDetElemGloY;
lSigmaY += lDetElemGloY*lDetElemGloY;
lMeanZ += lDetElemGloZ;
lSigmaZ += lDetElemGloZ*lDetElemGloZ;
lNDetElem++;
+
}
+
+ // calculate mean values
+ 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 ) );
+
+ // create all possible arrays
+ Array fConstraintX[4]; //Array for constraint equation X
+ Array fConstraintY[4]; //Array for constraint equation Y
+ Array fConstraintP[4]; //Array for constraint equation P
+ Array fConstraintXZ[4]; //Array for constraint equation X vs Z
+ Array fConstraintYZ[4]; //Array for constraint equation Y vs Z
+ Array fConstraintPZ[4]; //Array for constraint equation P vs Z
+
+ // do we really need these ?
+ Array fConstraintXY[4]; //Array for constraint equation X vs Y
+ Array fConstraintYY[4]; //Array for constraint equation Y vs Y
+ Array fConstraintPY[4]; //Array for constraint equation P vs Y
+
+ // fill Bool_t sides array based on masks, for convenience
+ Bool_t lDetTLBR[4];
+ lDetTLBR[0] = sidesMask & SideTop;
+ lDetTLBR[1] = sidesMask & SideLeft;
+ lDetTLBR[2] = sidesMask & SideBottom;
+ lDetTLBR[3] = sidesMask & SideRight;
+
+ for( Int_t i = 0; i < fgNDetElem; ++i )
+ {
+
+ // get chamber matching detector
+ const Int_t iCh( GetChamberId(i) );
+
+ // skip detector if chamber is off
+ if( !lChOnOff[iCh-1] ) continue;
+
+ // get detector element id from detector element number
+ const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
+ const Int_t lDetElemId = iCh*100+lDetElemNumber;
+
+ // get global x, y and z position
+ Double_t lDetElemGloX = 0.;
+ Double_t lDetElemGloY = 0.;
+ Double_t lDetElemGloZ = 0.;
+ fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
+
+ // loop over sides
+ for( Int_t iSide = 0; iSide < 4; iSide++ )
+ {
+
+ // skip if side is not selected
+ if( !lDetTLBR[iSide] ) continue;
+
+ // skip detector if it is not in the selected side
+ // stations 1 and 2
+ if( iCh>=1 && iCh<=4 )
+ {
+ if( lDetElemNumber == 0 && !(iSide == 0 || iSide == 3) ) continue; // top-right
+ if( lDetElemNumber == 1 && !(iSide == 0 || iSide == 1) ) continue; // top-left
+ if( lDetElemNumber == 2 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
+ if( lDetElemNumber == 3 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
}
- }
- }
- 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 (i<fgSNDetElemCh[iCh-1]) break;
- }
- if (lChOnOff[iCh-1]){
- Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
- Int_t lDetElemId = iCh*100+lDetElemNumber;
- fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
- 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 (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 (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 (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 (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 (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 (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 (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
+
+ // station 3
+ if (iCh>=5 && iCh<=6)
+ {
+ if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(iSide == 0 || iSide == 3) ) continue; // top-right
+ if( lDetElemNumber >= 5 && lDetElemNumber <= 9 && !(iSide == 0 || iSide == 1) ) continue; // top-left
+ if( lDetElemNumber >= 10 && lDetElemNumber <= 13 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
+ if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
}
- 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
+
+ // stations 4 and 5
+ if (iCh>=7 && iCh<=10)
+ {
+ if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(iSide == 0 || iSide == 3) ) continue; // top-right
+ if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(iSide == 0 || iSide == 1) ) continue; // top-left
+ if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
+ if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
}
- }
- }
- if (lVarXYT[0]) { // X constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
- }
- if (lVarXYT[1]) { // Y constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
- }
- if (lVarXYT[2]) { // T constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
- }
- if (lVarXYT[3]) { // X-Z constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintXZT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintXZL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintXZB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintXZR,0.0); // Right half
- }
- if (lVarXYT[4]) { // Y-Z constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintYZT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintYZL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintYZB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintYZR,0.0); // Right half
- }
- if (lVarXYT[5]) { // P-Z constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintPZT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintPZL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintPZB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintPZR,0.0); // Right half
- }
- if (lVarXYT[6]) { // X-Y constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintXYT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintXYL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintXYB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintXYR,0.0); // Right half
- }
- if (lVarXYT[7]) { // Y-Y constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintYYT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintYYL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintYYB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintYYR,0.0); // Right half
- }
- if (lVarXYT[8]) { // P-Y constraint
- if (lDetTLBR[0]) AddConstraint(fConstraintPYT,0.0); // Top half
- if (lDetTLBR[1]) AddConstraint(fConstraintPYL,0.0); // Left half
- if (lDetTLBR[2]) AddConstraint(fConstraintPYB,0.0); // Bottom half
- if (lDetTLBR[3]) AddConstraint(fConstraintPYR,0.0); // Right half
- }
-}
-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){
- if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
- lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
- }
- }
- if (lCh>=5 && lCh<=6){
- if (lDetElemNumber>=0&&lDetElemNumber<=9){
- lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
+ // constrain x
+ if( lVarXYT[0] ) fConstraintX[iSide].values[i*fgNParCh+0] = 1;
+
+ // constrain y
+ if( lVarXYT[1] ) fConstraintY[iSide].values[i*fgNParCh+1] = 1;
+
+ // constrain phi (rotation around z)
+ if( lVarXYT[2] ) fConstraintP[iSide].values[i*fgNParCh+2] = 1;
+
+ // x-z shearing
+ if( lVarXYT[3] ) fConstraintXZ[iSide].values[i*fgNParCh+0] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
+
+ // y-z shearing
+ if( lVarXYT[4] ) fConstraintYZ[iSide].values[i*fgNParCh+1] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
+
+ // phi-z shearing
+ if( lVarXYT[5] ) fConstraintPZ[iSide].values[i*fgNParCh+2] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
+
+ // x-y shearing
+ if( lVarXYT[6] ) fConstraintXY[iSide].values[i*fgNParCh+0] = (lDetElemGloY-lMeanY)/lSigmaY;
+
+ // y-y shearing
+ if( lVarXYT[7] ) fConstraintYY[iSide].values[i*fgNParCh+1] = (lDetElemGloY-lMeanY)/lSigmaY;
+
+ // phi-y shearing
+ if( lVarXYT[8] ) fConstraintPY[iSide].values[i*fgNParCh+2] = (lDetElemGloY-lMeanY)/lSigmaY;
+
}
+
}
- if (lCh>=7 && lCh<=10){
- if (lDetElemNumber>=0&&lDetElemNumber<=13){
- lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
- }
+
+ // pass constraints to millepede
+ for( Int_t iSide = 0; iSide < 4; iSide++ )
+ {
+ // skip if side is not selected
+ if( !lDetTLBR[iSide] ) continue;
+
+ if( lVarXYT[0] ) AddConstraint(fConstraintX[iSide].values,0.0);
+ if( lVarXYT[1] ) AddConstraint(fConstraintY[iSide].values,0.0);
+ if( lVarXYT[2] ) AddConstraint(fConstraintP[iSide].values,0.0);
+ if( lVarXYT[3] ) AddConstraint(fConstraintXZ[iSide].values,0.0);
+ if( lVarXYT[4] ) AddConstraint(fConstraintYZ[iSide].values,0.0);
+ if( lVarXYT[5] ) AddConstraint(fConstraintPZ[iSide].values,0.0);
+ if( lVarXYT[6] ) AddConstraint(fConstraintXY[iSide].values,0.0);
+ if( lVarXYT[7] ) AddConstraint(fConstraintYY[iSide].values,0.0);
+ if( lVarXYT[8] ) AddConstraint(fConstraintPY[iSide].values,0.0);
}
+
}
//______________________________________________________________________
-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){
- if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
- lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
- }
- }
- if (lCh>=5 && lCh<=6){
- if (lDetElemNumber>=5&&lDetElemNumber<=13){
- lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
- }
- }
- if (lCh>=7 && lCh<=10){
- if (lDetElemNumber>=7&&lDetElemNumber<=19){
- lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
- }
- }
+void AliMUONAlignment::InitGlobalParameters(Double_t *par)
+{
+ /// Initialize global parameters with par array
+ if( !fInitialized )
+ { AliFatal( "Millepede is not initialized" ); }
+
+ fMillepede->SetGlobalParameters(par);
}
//______________________________________________________________________
-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){
- if (lDetElemNumber==2 || lDetElemNumber==3){ // From track crossings
- lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
- }
- }
- if (lCh>=5 && lCh<=6){
- 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)){
- lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
- }
- }
+void AliMUONAlignment::SetAllowedVariation( Int_t iPar, Double_t value )
+{
+ /// "Encouraged" variation for degrees of freedom
+ // check initialization
+ if( fInitialized )
+ { AliFatal( "Millepede already initialized" ); }
+
+ // check initialization
+ if( !(iPar >= 0 && iPar < fgNParCh ) )
+ { AliFatal( Form( "Invalid index: %i", iPar ) ); }
+
+ fAllowVar[iPar] = value;
}
//______________________________________________________________________
-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){
- if (lDetElemNumber==0 || lDetElemNumber==3){ // From track crossings
- lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
- }
- }
- if (lCh>=5 && lCh<=6){
- 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)){
- lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
- }
- }
+void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
+{
+
+ /// Set expected measurement resolution
+ fSigma[0] = sigmaX;
+ fSigma[1] = sigmaY;
+
+ // print
+ for( Int_t i=0; i<2; ++i )
+ { AliInfo( Form( "fSigma[%i]=%f", i, fSigma[i] ) ); }
+
}
-//______________________________________________________________________
-void AliMUONAlignment::ResetConstraints(){
- /// Reset all constraint equations
- for (Int_t i = 0; i < fgNDetElem; i++){
- fConstraintX3[i*fgNParCh+0]=0.0;
- fConstraintY3[i*fgNParCh+0]=0.0;
- fConstraintX4[i*fgNParCh+0]=0.0;
- fConstraintY4[i*fgNParCh+0]=0.0;
- fConstraintP4[i*fgNParCh+0]=0.0;
- fConstraintX5[i*fgNParCh+0]=0.0;
- fConstraintY5[i*fgNParCh+0]=0.0;
- fConstraintX[i*fgNParCh+0]=0.0;
- fConstraintX[i*fgNParCh+1]=0.0;
- fConstraintX[i*fgNParCh+2]=0.0;
- fConstraintY[i*fgNParCh+0]=0.0;
- fConstraintY[i*fgNParCh+1]=0.0;
- fConstraintY[i*fgNParCh+2]=0.0;
- fConstraintP[i*fgNParCh+0]=0.0;
- fConstraintP[i*fgNParCh+1]=0.0;
- fConstraintP[i*fgNParCh+2]=0.0;
- fConstraintXT[i*fgNParCh+0]=0.0;
- fConstraintXT[i*fgNParCh+1]=0.0;
- fConstraintXT[i*fgNParCh+2]=0.0;
- fConstraintYT[i*fgNParCh+0]=0.0;
- fConstraintYT[i*fgNParCh+1]=0.0;
- fConstraintYT[i*fgNParCh+2]=0.0;
- fConstraintPT[i*fgNParCh+0]=0.0;
- fConstraintPT[i*fgNParCh+1]=0.0;
- fConstraintPT[i*fgNParCh+2]=0.0;
- fConstraintXZT[i*fgNParCh+0]=0.0;
- fConstraintXZT[i*fgNParCh+1]=0.0;
- fConstraintXZT[i*fgNParCh+2]=0.0;
- fConstraintYZT[i*fgNParCh+0]=0.0;
- fConstraintYZT[i*fgNParCh+1]=0.0;
- fConstraintYZT[i*fgNParCh+2]=0.0;
- fConstraintPZT[i*fgNParCh+0]=0.0;
- fConstraintPZT[i*fgNParCh+1]=0.0;
- fConstraintPZT[i*fgNParCh+2]=0.0;
- fConstraintXYT[i*fgNParCh+0]=0.0;
- fConstraintXYT[i*fgNParCh+1]=0.0;
- fConstraintXYT[i*fgNParCh+2]=0.0;
- fConstraintYYT[i*fgNParCh+0]=0.0;
- fConstraintYYT[i*fgNParCh+1]=0.0;
- fConstraintYYT[i*fgNParCh+2]=0.0;
- fConstraintPYT[i*fgNParCh+0]=0.0;
- fConstraintPYT[i*fgNParCh+1]=0.0;
- fConstraintPYT[i*fgNParCh+2]=0.0;
- fConstraintXL[i*fgNParCh+0]=0.0;
- fConstraintXL[i*fgNParCh+1]=0.0;
- fConstraintXL[i*fgNParCh+2]=0.0;
- fConstraintYL[i*fgNParCh+0]=0.0;
- fConstraintYL[i*fgNParCh+1]=0.0;
- fConstraintYL[i*fgNParCh+2]=0.0;
- fConstraintPL[i*fgNParCh+0]=0.0;
- fConstraintPL[i*fgNParCh+1]=0.0;
- fConstraintPL[i*fgNParCh+2]=0.0;
- fConstraintXZL[i*fgNParCh+0]=0.0;
- fConstraintXZL[i*fgNParCh+1]=0.0;
- fConstraintXZL[i*fgNParCh+2]=0.0;
- fConstraintYZL[i*fgNParCh+0]=0.0;
- fConstraintYZL[i*fgNParCh+1]=0.0;
- fConstraintYZL[i*fgNParCh+2]=0.0;
- fConstraintPZL[i*fgNParCh+0]=0.0;
- fConstraintPZL[i*fgNParCh+1]=0.0;
- fConstraintPZL[i*fgNParCh+2]=0.0;
- fConstraintXYL[i*fgNParCh+0]=0.0;
- fConstraintXYL[i*fgNParCh+1]=0.0;
- fConstraintXYL[i*fgNParCh+2]=0.0;
- fConstraintYYL[i*fgNParCh+0]=0.0;
- fConstraintYYL[i*fgNParCh+1]=0.0;
- fConstraintYYL[i*fgNParCh+2]=0.0;
- fConstraintPYL[i*fgNParCh+0]=0.0;
- fConstraintPYL[i*fgNParCh+1]=0.0;
- fConstraintPYL[i*fgNParCh+2]=0.0;
- fConstraintXB[i*fgNParCh+0]=0.0;
- fConstraintXB[i*fgNParCh+1]=0.0;
- fConstraintXB[i*fgNParCh+2]=0.0;
- fConstraintYB[i*fgNParCh+0]=0.0;
- fConstraintYB[i*fgNParCh+1]=0.0;
- fConstraintYB[i*fgNParCh+2]=0.0;
- fConstraintPB[i*fgNParCh+0]=0.0;
- fConstraintPB[i*fgNParCh+1]=0.0;
- fConstraintPB[i*fgNParCh+2]=0.0;
- fConstraintXZB[i*fgNParCh+0]=0.0;
- fConstraintXZB[i*fgNParCh+1]=0.0;
- fConstraintXZB[i*fgNParCh+2]=0.0;
- fConstraintYZB[i*fgNParCh+0]=0.0;
- fConstraintYZB[i*fgNParCh+1]=0.0;
- fConstraintYZB[i*fgNParCh+2]=0.0;
- fConstraintPZB[i*fgNParCh+0]=0.0;
- fConstraintPZB[i*fgNParCh+1]=0.0;
- fConstraintPZB[i*fgNParCh+2]=0.0;
- fConstraintXYB[i*fgNParCh+0]=0.0;
- fConstraintXYB[i*fgNParCh+1]=0.0;
- fConstraintXYB[i*fgNParCh+2]=0.0;
- fConstraintYYB[i*fgNParCh+0]=0.0;
- fConstraintYYB[i*fgNParCh+1]=0.0;
- fConstraintYYB[i*fgNParCh+2]=0.0;
- fConstraintPYB[i*fgNParCh+0]=0.0;
- fConstraintPYB[i*fgNParCh+1]=0.0;
- fConstraintPYB[i*fgNParCh+2]=0.0;
- fConstraintXR[i*fgNParCh+0]=0.0;
- fConstraintXR[i*fgNParCh+1]=0.0;
- fConstraintXR[i*fgNParCh+2]=0.0;
- fConstraintYR[i*fgNParCh+0]=0.0;
- fConstraintYR[i*fgNParCh+1]=0.0;
- fConstraintYR[i*fgNParCh+2]=0.0;
- fConstraintPR[i*fgNParCh+0]=0.0;
- fConstraintPR[i*fgNParCh+1]=0.0;
- fConstraintPR[i*fgNParCh+2]=0.0;
- fConstraintXZR[i*fgNParCh+0]=0.0;
- fConstraintXZR[i*fgNParCh+1]=0.0;
- fConstraintXZR[i*fgNParCh+2]=0.0;
- fConstraintYZR[i*fgNParCh+0]=0.0;
- fConstraintYZR[i*fgNParCh+1]=0.0;
- fConstraintYZR[i*fgNParCh+2]=0.0;
- fConstraintPZR[i*fgNParCh+0]=0.0;
- fConstraintPZR[i*fgNParCh+1]=0.0;
- fConstraintPZR[i*fgNParCh+2]=0.0;
- fConstraintPZR[i*fgNParCh+0]=0.0;
- fConstraintPZR[i*fgNParCh+1]=0.0;
- fConstraintPZR[i*fgNParCh+2]=0.0;
- fConstraintXYR[i*fgNParCh+0]=0.0;
- fConstraintXYR[i*fgNParCh+1]=0.0;
- fConstraintXYR[i*fgNParCh+2]=0.0;
- fConstraintYYR[i*fgNParCh+0]=0.0;
- fConstraintYYR[i*fgNParCh+1]=0.0;
- fConstraintYYR[i*fgNParCh+2]=0.0;
- fConstraintPYR[i*fgNParCh+0]=0.0;
- fConstraintPYR[i*fgNParCh+1]=0.0;
- fConstraintPYR[i*fgNParCh+2]=0.0;
+//_____________________________________________________
+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 iDet=0; iDet<fgNDetElem; ++iDet )
+ {
+ AliInfo( Form( "%d\t %f\t %f\t %f\t %f",
+ iDet,
+ parameters[iDet*fgNParCh+0],parameters[iDet*fgNParCh+1],
+ parameters[iDet*fgNParCh+3],parameters[iDet*fgNParCh+2]
+ ) );
}
+
}
-//______________________________________________________________________
-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::PrintGlobalParameters() const
+{ fMillepede->PrintGlobalParameters(); }
+
+//_____________________________________________________
+Double_t AliMUONAlignment::GetParError(Int_t iPar) const
+{ return fMillepede->GetParError(iPar); }
+
+//______________________________________________________________________
+AliMUONGeometryTransformer* AliMUONAlignment::ReAlign(
+ const AliMUONGeometryTransformer * transformer,
+ const double *misAlignments, Bool_t )
+{
+
+ /// 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
+ // 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[fgNParCh] = {0};
+ Double_t lDetElemMisAlignment[fgNParCh] = {0};
+ const TClonesArray* oldMisAlignArray( transformer->GetMisAlignmentData() );
-//______________________________________________________________________
-void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
- /// Initialize global parameters with par array
- fMillepede->SetGlobalParameters(par);
- AliInfo("Init Global Parameters");
-}
+ AliMUONGeometryTransformer *newGeometryTransformer = new AliMUONGeometryTransformer();
+ for( Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); ++iMt )
+ {
-//______________________________________________________________________
-void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
- /// Parameter iPar is encourage to vary in [-value;value].
- /// If value == 0, parameter is fixed
- fMillepede->SetParSigma(iPar, value);
- if (TMath::Abs(value)<1e-4) AliInfo(Form("Parameter %i Fixed", iPar));
-}
+ // module transformers
+ const AliMUONGeometryModuleTransformer *kModuleTransformer = transformer->GetModuleTransformer(iMt, kTRUE);
-//______________________________________________________________________
-void AliMUONAlignment::ResetLocalEquation()
-{
- /// Reset the derivative vectors
- for(int i=0; i<fNLocal; i++) {
- fLocalDerivatives[i] = 0.0;
- }
- for(int i=0; i<fNGlobal; i++) {
- fGlobalDerivatives[i] = 0.0;
- }
-}
+ AliMUONGeometryModuleTransformer *newModuleTransformer = new AliMUONGeometryModuleTransformer(iMt);
+ newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
-//______________________________________________________________________
-void AliMUONAlignment::AllowVariations(const Bool_t *bChOnOff)
-{
+ // get transformation
+ TGeoHMatrix deltaModuleTransform( DeltaTransform( lModuleMisAlignment ) );
- /// Set allowed variation for selected chambers based on fDoF and fAllowVar
- for (Int_t iCh=1; iCh<=10; iCh++)
- {
- if (bChOnOff[iCh-1])
+ // update module
+ TGeoHMatrix moduleTransform( *kModuleTransformer->GetTransformation() );
+ TGeoHMatrix newModuleTransform( AliMUONGeometryBuilder::Multiply( deltaModuleTransform, moduleTransform ) );
+ newModuleTransformer->SetTransformation(newModuleTransform);
+
+ // Get matching old alignment and update current matrix accordingly
+ if( oldMisAlignArray )
{
- Int_t iDetElemFirst = (iCh>1) ? fgSNDetElemCh[iCh-2] : 0;
- Int_t iDetElemLast = fgSNDetElemCh[iCh-1];
- for (int i=0; i<fgNParCh; i++)
+ const AliAlignObjMatrix* oldAlignObj(0);
+ const Int_t moduleId( kModuleTransformer->GetModuleId() );
+ const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, moduleId );
+ for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
{
- AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));
- if (fDoF[i])
+
+ const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
+ if( localAlignObj && localAlignObj->GetVolUID() == volId )
{
- for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){
- FixParameter(j*fgNParCh+i, fAllowVar[i]);
- }
+ oldAlignObj = localAlignObj;
+ break;
}
}
+ // multiply
+ if( oldAlignObj )
+ {
+
+ TGeoHMatrix oldMatrix;
+ oldAlignObj->GetMatrix( oldMatrix );
+ deltaModuleTransform.Multiply( &oldMatrix );
+
+ }
+
}
- }
+ // Create module mis alignment matrix
+ newGeometryTransformer ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
-}
+ AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
-//______________________________________________________________________
-void AliMUONAlignment::SetNonLinear(Int_t iPar /* set non linear flag */ )
-{
- /// Set nonlinear flag for parameter iPar
- fMillepede->SetNonLinear(iPar);
- AliInfo(Form("Parameter %i set to non linear", iPar));
-}
+ TIter next(detElements->CreateIterator());
+ AliMUONGeometryDetElement* detElement;
+ Int_t iDe(-1);
+ while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
+ {
+ ++iDe;
+ // make a new detection element
+ AliMUONGeometryDetElement *newDetElement = new AliMUONGeometryDetElement(detElement->GetId(), detElement->GetVolumePath());
+ TString lDetElemName(detElement->GetDEName());
+ lDetElemName.ReplaceAll("DE","");
-//______________________________________________________________________
-void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
-{
+ // store detector element id and number
+ const Int_t iDetElemId = lDetElemName.Atoi();
+ if( !DetElemIsValid( iDetElemId ) )
+ {
+ AliInfo( Form( "Skipping invalid detector element %i", iDetElemId ) );
+ continue;
+ }
- /// Set expected measurement resolution
- fSigma[0] = sigmaX; fSigma[1] = sigmaY;
- AliInfo(Form("Using fSigma[0]=%f and fSigma[1]=%f",fSigma[0],fSigma[1]));
-}
+ const Int_t iDetElemNumber( GetDetElemNumber( iDetElemId ) );
+ for( int i=0; i<fgNParCh; ++i )
+ {
+ lDetElemMisAlignment[i] = 0.0;
+ if( iMt<fgNTrkMod ) { lDetElemMisAlignment[i] = misAlignments[iDetElemNumber*fgNParCh+i]; }
+ }
-//______________________________________________________________________
-void AliMUONAlignment::LocalEquationX( Bool_t doAlignment )
-{
+ // get transformation
+ TGeoHMatrix deltaGlobalTransform( DeltaTransform( lDetElemMisAlignment ) );
+ // update module
+ TGeoHMatrix globalTransform( *detElement->GetGlobalTransformation() );
+ TGeoHMatrix newGlobalTransform( AliMUONGeometryBuilder::Multiply( deltaGlobalTransform, globalTransform ) );
+ newDetElement->SetGlobalTransformation( newGlobalTransform );
+ newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
- // local cluster record
- AliMUONAlignmentClusterRecord clusterRecord;
+ // Get matching old alignment and update current matrix accordingly
+ if( oldMisAlignArray )
+ {
- // store detector and measurement
- clusterRecord.SetDetElemId( fDetElemId );
- clusterRecord.SetDetElemNumber( fDetElemNumber );
- clusterRecord.SetMeas( fMeas[0] );
- clusterRecord.SetSigma( fSigma[0] );
+ const AliAlignObjMatrix* oldAlignObj(0);
+ const int detElemId( detElement->GetId() );
+ const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId );
+ for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
+ {
- // 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]) );
+ const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
+ if( localAlignObj && localAlignObj->GetVolUID() == volId )
+ {
+ oldAlignObj = localAlignObj;
+ break;
+ }
- // store global derivatives
- clusterRecord.SetGlobalDerivative( 0, -fCosPhi );
- clusterRecord.SetGlobalDerivative( 1, -fSinPhi );
+ }
- if (fBFieldOn)
- {
+ // multiply
+ if( oldAlignObj )
+ {
- clusterRecord.SetGlobalDerivative(
- 2,
- -fSinPhi*(fTrackPos[0]-fDetElemPos[0])
- +fCosPhi*(fTrackPos[1]-fDetElemPos[1]) );
+ TGeoHMatrix oldMatrix;
+ oldAlignObj->GetMatrix( oldMatrix );
+ deltaGlobalTransform.Multiply( &oldMatrix );
- } 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]) );
+ }
- }
+ // Create misalignment matrix
+ newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
- clusterRecord.SetGlobalDerivative( 3, fCosPhi*fTrackSlope0[0]+fSinPhi*fTrackSlope0[1] );
+ }
- // append to trackRecord
- fTrackRecord.AddClusterRecord( clusterRecord );
+ newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
+ }
- // store local equation
- if( doAlignment ) LocalEquation( clusterRecord );
+ return newGeometryTransformer;
}
//______________________________________________________________________
-void AliMUONAlignment::LocalEquationY(Bool_t doAlignment )
+void AliMUONAlignment::SetAlignmentResolution( const TClonesArray* misAlignArray, Int_t rChId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY )
{
- // local cluster record
- AliMUONAlignmentClusterRecord clusterRecord;
-
- // store detector and measurement
- clusterRecord.SetDetElemId( fDetElemId );
- clusterRecord.SetDetElemNumber( fDetElemNumber );
- clusterRecord.SetMeas( fMeas[1] );
- clusterRecord.SetSigma( fSigma[1] );
+ /// Set alignment resolution to misalign objects to be stored in CDB
+ /// if rChId is > 0 set parameters for this chamber only, counting from 1
+ TMatrixDSym mChCorrMatrix(6);
+ mChCorrMatrix[0][0]=chResX*chResX;
+ mChCorrMatrix[1][1]=chResY*chResY;
- // 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] ) );
+ TMatrixDSym mDECorrMatrix(6);
+ mDECorrMatrix[0][0]=deResX*deResX;
+ mDECorrMatrix[1][1]=deResY*deResY;
- // set global derivatives
- clusterRecord.SetGlobalDerivative( 0, fSinPhi);
- clusterRecord.SetGlobalDerivative( 1, -fCosPhi);
+ AliAlignObjMatrix *alignMat = 0x0;
- if (fBFieldOn)
+ for( Int_t chId = 0; chId <= 9; ++chId )
{
- clusterRecord.SetGlobalDerivative(
- 2,
- -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
- -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
+ // skip chamber if selection is valid, and does not match
+ if( rChId > 0 && chId+1 != rChId ) continue;
- } else {
+ TString chName1;
+ TString chName2;
+ if (chId<4)
+ {
- 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]));
- }
+ chName1 = Form("GM%d",chId);
+ chName2 = Form("GM%d",chId);
- clusterRecord.SetGlobalDerivative( 3, -fSinPhi*fTrackSlope0[0]+fCosPhi*fTrackSlope0[1]);
+ } else {
- // append to trackRecord
- fTrackRecord.AddClusterRecord( clusterRecord );
+ chName1 = Form("GM%d",4+(chId-4)*2);
+ chName2 = Form("GM%d",4+(chId-4)*2+1);
- // store local equation
- if( doAlignment ) LocalEquation( clusterRecord );
+ }
+
+ for( int i=0; i<misAlignArray->GetEntries(); ++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->SetCorrMatrix(mChCorrMatrix);
+ else if (volName.Contains("DE")) alignMat->SetCorrMatrix(mDECorrMatrix);
+
+ }
+
+ }
+
+ }
}
+
//_____________________________________________________
-void AliMUONAlignment::LocalEquation( const AliMUONAlignmentClusterRecord& clusterRecord )
+void AliMUONAlignment::FillDetElemData( AliMUONVCluster* cluster )
{
- // 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 ) );
- }
+ /// Get information of current detection element
+ // get detector element number from Alice ID
+ const Int_t detElemId = cluster->GetDetElemId();
+ fDetElemNumber = GetDetElemNumber( detElemId );
+
+ // get detector element
+ const AliMUONGeometryDetElement* detElement = fTransform->GetDetElement( detElemId );
- // pass equation parameters to millepede
- fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, clusterRecord.GetMeas(), clusterRecord.GetSigma() );
+ /*
+ get the global transformation matrix and store its inverse, in order to manually perform
+ the global to Local transformations needed to calculate the derivatives
+ */
+ fGeoCombiTransInverse = detElement->GetGlobalTransformation()->Inverse();
}
//______________________________________________________________________
-void AliMUONAlignment::FillRecPointData()
+void AliMUONAlignment::FillRecPointData( AliMUONVCluster* cluster )
{
/// 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]);
+ fClustPos[0] = cluster->GetX();
+ fClustPos[1] = cluster->GetY();
+ fClustPos[2] = cluster->GetZ();
}
//______________________________________________________________________
-void AliMUONAlignment::FillTrackParamData()
+void AliMUONAlignment::FillTrackParamData( AliMUONTrackParam* trackParam )
{
/// 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]);
+ fTrackPos[0] = trackParam->GetNonBendingCoor();
+ fTrackPos[1] = trackParam->GetBendingCoor();
+ fTrackPos[2] = trackParam->GetZ();
+ fTrackSlope[0] = trackParam->GetNonBendingSlope();
+ fTrackSlope[1] = trackParam->GetBendingSlope();
}
-//_____________________________________________________
-void AliMUONAlignment::FillDetElemData()
+//______________________________________________________________________
+Bool_t AliMUONAlignment::UnbiasTrackParamData( AliMUONTrackParam* trackParam ) const
{
- /// Get information of current detection element
- Double_t lDetElemLocX = 0.;
- Double_t lDetElemLocY = 0.;
- Double_t lDetElemLocZ = 0.;
- fDetElemId = fCluster->GetDetElemId();
- fDetElemNumber = fDetElemId%100;
- for (int iCh=0; iCh<fDetElemId/100-1; iCh++)
- { fDetElemNumber += fgNDetElemCh[iCh]; }
-
- fTransform->Local2Global(
- fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
- fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
+ /**
+ calculate unbiased track parameters at a given detector, that is,
+ after taking out the contribution of the detector's cluster from the track
+ */
+
+ // check track parameters
+ if( !trackParam ) return kFALSE;
+
+ // Remove cluster contibution from smoothed track param
+ TMatrixD smoothParameters( trackParam->GetSmoothParameters() );
+ TMatrixD smoothCovariances( trackParam->GetSmoothCovariances() );
+
+ AliMUONVCluster* cluster = trackParam->GetClusterPtr();
+ // p' = p + K(m - H*p)
+ // K = C H (-V + H C H^t)^-1
+ // C' = C - K H C
+ // where p,C are smoothed param,cov at cluster position
+ // H is the matrix converting the state vector to measurement
+ // V is the measurement cov.matrix
+ // m is the measurement vector
+ static TMatrixD H(2,5);
+ H.Zero();
+ H(0,0)=H(1,2) = 1.;
+
+ // (-Vk+H_k C^n_k H_k^T)^-1
+ static TMatrixD df(2,2);
+ df.Zero();
+ df(0,0) = smoothCovariances(0,0) - cluster->GetErrX2();
+ df(1,1) = smoothCovariances(2,2) - cluster->GetErrY2();
+ df(0,1) = smoothCovariances(0,2);
+ df(1,0) = smoothCovariances(2,0);
+
+ if (df.Determinant() != 0) df.Invert();
+ else {
+ AliInfo( "Determinant = 0\n" );
+ return kFALSE;
+ }
-}
+ // gain matrix
+ TMatrixD kTmp( smoothCovariances, TMatrixD::kMultTranspose, H );
+ TMatrixD K(kTmp, TMatrixD::kMult, df);
-//_____________________________________________________
-AliMUONAlignmentTrackRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment )
-{
- /// Process track and set Local Equations
- // store current track in running member.
- fTrack = track;
+ TMatrixD dfc(2,1);
+ dfc.Zero();
+ dfc(0,0) = cluster->GetX()-smoothParameters(0,0);
+ dfc(1,0) = cluster->GetY()-smoothParameters(2,0);
+ TMatrixD tmp0(K,TMatrixD::kMult, dfc);
+ smoothParameters += tmp0;
- // clear track record
- fTrackRecord.Clear();
+ TMatrixD tmp1(K, TMatrixD::kMult, H);
+ TMatrixD tmp2(tmp1,TMatrixD::kMult, smoothCovariances);
+ smoothCovariances -= tmp2;
- // get number of tracks
- Int_t nTrackParam = fTrack->GetTrackParamAtCluster()->GetEntries();
- AliDebug(1,Form("Number of track param entries : %i ", nTrackParam));
+ // update track parameters
+ trackParam->SetParameters( smoothParameters );
+ trackParam->SetCovariances( smoothCovariances );
- Bool_t first( kTRUE );
- for(Int_t iCluster=0; iCluster<nTrackParam; iCluster++)
- {
+ return kTRUE;
- // and get new pointers
- 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();
- FillRecPointData();
- FillTrackParamData();
+//______________________________________________________________________
+void AliMUONAlignment::LocalEquationX( void )
+{
+ /// local equation along X
+
+ // 'inverse' (GlobalToLocal) rotation matrix
+ const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
+
+ // local derivatives
+ SetLocalDerivative( 0, r[0] );
+ SetLocalDerivative( 1, r[0]*(fTrackPos[2] - fTrackPos0[2]) );
+ SetLocalDerivative( 2, r[1] );
+ SetLocalDerivative( 3, r[1]*(fTrackPos[2] - fTrackPos0[2]) );
+
+ // global derivatives
+ /*
+ alignment parameters are
+ 0: delta_x
+ 1: delta_y
+ 2: delta_phiz
+ 3: delta_z
+ */
+
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[0] );
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[1] );
+
+ if( fBFieldOn )
+ {
- if( first )
- {
+ // use local position for derivatives vs 'delta_phi_z'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*fTrackPos[0] + r[0]*fTrackPos[1] );
- // for first valid cluster, save track position as "starting" values
- first = kFALSE;
+ // use local slopes for derivatives vs 'delta_z'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope[0] + r[1]*fTrackSlope[1] );
- fTrackPos0[0] = fTrackPos[0];
- fTrackPos0[1] = fTrackPos[1];
- fTrackPos0[2] = fTrackPos[2];
- fTrackSlope0[0] = fTrackSlope[0];
- fTrackSlope0[1] = fTrackSlope[1];
+ } else {
- }
+ // local copy of extrapolated track positions
+ const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
+ const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
- // calculate measurements
- fCosPhi = TMath::Cos(fPhi);
- fSinPhi = TMath::Sin(fPhi);
- if( fBFieldOn )
- {
+ // use properly extrapolated position for derivatives vs 'delta_phi_z'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*trackPosX + r[0]*trackPosY );
- fMeas[0] = fTrackPos[0] - fClustPos[0];
- fMeas[1] = fTrackPos[1] - fClustPos[1];
+ // use slopes at origin for derivatives vs 'delta_z'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope0[0] + r[1]*fTrackSlope0[1] );
- } else {
+ }
- fMeas[0] = - fClustPos[0];
- fMeas[1] = - fClustPos[1];
+ // store local equation
+ fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0] );
- }
+}
- // 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));
+//______________________________________________________________________
+void AliMUONAlignment::LocalEquationY( void )
+{
+ /// local equation along Y
- 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]));
+ // 'inverse' (GlobalToLocal) rotation matrix
+ const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
- 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]));
+ // store local derivatives
+ SetLocalDerivative( 0, r[3] );
+ SetLocalDerivative( 1, r[3]*(fTrackPos[2] - fTrackPos0[2] ) );
+ SetLocalDerivative( 2, r[4] );
+ SetLocalDerivative( 3, r[4]*(fTrackPos[2] - fTrackPos0[2] ) );
- // Set local equations
- LocalEquationX( doAlignment );
- LocalEquationY( doAlignment );
+ // set global derivatives
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[3]);
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[4]);
- }
+ if( fBFieldOn )
+ {
- return &fTrackRecord;
+ // use local position for derivatives vs 'delta_phi'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*fTrackPos[0] + r[3]*fTrackPos[1]);
-}
+ // use local slopes for derivatives vs 'delta_z'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope[0]+r[4]*fTrackSlope[1] );
-//______________________________________________________________________________
-void AliMUONAlignment::ProcessTrack( AliMUONAlignmentTrackRecord* track, Bool_t doAlignment )
-{
- /// Process track (from record) and set Local Equations
- if( !( track && doAlignment ) ) return;
+ } else {
- // loop over clusters
- for( Int_t index = 0; index < track->GetNRecords(); ++index )
- { if( AliMUONAlignmentClusterRecord* record = track->GetRecord( index ) ) LocalEquation( *record ); }
+ // local copy of extrapolated track positions
+ const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
+ const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
- return;
+ // use properly extrapolated position for derivatives vs 'delta_phi'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*trackPosX + r[3]*trackPosY );
-}
+ // use slopes at origin for derivatives vs 'delta_z'
+ SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope0[0]+r[4]*fTrackSlope0[1] );
-//_____________________________________________________
-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); }
+ // store local equation
+ fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1] );
}
-//_____________________________________________________
-void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls)
+//_________________________________________________________________________
+TGeoCombiTrans AliMUONAlignment::DeltaTransform( const double *lMisAlignment) const
{
+ /// Get Delta Transformation, based on alignment parameters
- /// Call global fit; Global parameters are stored in parameters
- fMillepede->GlobalFit(parameters,errors,pulls);
+ // translation
+ const TGeoTranslation deltaTrans( lMisAlignment[0], lMisAlignment[1], lMisAlignment[3]);
- AliInfo("Done fitting global parameters!");
- for (int iGlob=0; iGlob<fgNDetElem; iGlob++)
- { printf("%d\t %f\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+3],parameters[iGlob*fgNParCh+2]); }
+ // rotation
+ TGeoRotation deltaRot;
+ deltaRot.RotateZ(lMisAlignment[2]*180./TMath::Pi());
+ // combined rotation and translation.
+ return TGeoCombiTrans(deltaTrans,deltaRot);
}
-//_____________________________________________________
-Double_t AliMUONAlignment::GetParError(Int_t iPar)
+//______________________________________________________________________
+void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value)
{
- /// Get error of parameter iPar
- Double_t lErr = fMillepede->GetParError(iPar);
- return lErr;
-}
+ /// Constrain equation defined by par to value
+ if( !fInitialized )
+ { AliFatal( "Millepede is not initialized" ); }
-//_____________________________________________________
-void AliMUONAlignment::PrintGlobalParameters()
-{
- /// Print global parameters
- fMillepede->PrintGlobalParameters();
+ fMillepede->SetGlobalConstraint(par, value);
}
-//_________________________________________________________________________
-TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, const double *lMisAlignment) const
+//______________________________________________________________________
+Bool_t AliMUONAlignment::DetElemIsValid( Int_t iDetElemId ) 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()) {
-// rot = transform.GetRotation();
-// }
-// else {
-// rot = new TGeoRotation("rot");
-// } // default constructor.
-
- 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;
- deltaRot.RotateX(angMisAlig[0]);
- deltaRot.RotateY(angMisAlig[1]);
- deltaRot.RotateZ(angMisAlig[2]);
-
- TGeoCombiTrans deltaTransf(deltaTrans,deltaRot);
- TGeoHMatrix newTransfMat = transform * deltaTransf;
-
- return TGeoCombiTrans(newTransfMat);
+ /// return true if given detector element is valid (and belongs to muon tracker)
+ const Int_t iCh = iDetElemId/100;
+ const Int_t iDet = iDetElemId%100;
+ return ( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] );
}
//______________________________________________________________________
-AliMUONGeometryTransformer *
-AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
- const double *misAlignments, Bool_t verbose)
-
+Int_t AliMUONAlignment::GetDetElemNumber( Int_t iDetElemId ) const
{
- /// 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
- // 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[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();
- for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) {
- // 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:
- // Tdelta = Tnew * Told.inverse
- TGeoHMatrix deltaModuleTransform =
- AliMUONGeometryBuilder::Multiply(newModuleTransform,
- kModuleTransformer->GetTransformation()->Inverse());
- // Create module mis alignment matrix
- newGeometryTransformer
- ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
-
- AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
+ /// get det element number from ID
+ // get chamber and element number in chamber
+ const Int_t iCh = iDetElemId/100;
+ const Int_t iDet = iDetElemId%100;
- if (verbose)
- AliInfo(Form("%i DEs in old GeometryStore %i",detElements->GetSize(), iMt));
+ // make sure detector index is valid
+ if( !( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] ) )
+ { AliFatal( Form( "Invalid detector element id: %i", iDetElemId ) ); }
- TIter next(detElements->CreateIterator());
- AliMUONGeometryDetElement* detElement;
- Int_t iDe(-1);
- while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
- {
- ++iDe;
- // make a new detection element
- AliMUONGeometryDetElement *newDetElement =
- new AliMUONGeometryDetElement(detElement->GetId(),
- detElement->GetVolumePath());
- TString lDetElemName(detElement->GetDEName());
- lDetElemName.ReplaceAll("DE","");
- iDetElemId = lDetElemName.Atoi();
- iDetElemNumber = iDetElemId%100;
- for (int iCh=0; iCh<iDetElemId/100-1; iCh++){
- iDetElemNumber += fgNDetElemCh[iCh];
- }
- for (int i=0; i<fgNParCh; i++) {
- lDetElemMisAlignment[i] = 0.0;
- if (iMt<fgNTrkMod) {
- AliInfo(Form("iMt %i, iCh %i, iDe %i, iDeId %i, iDeNb %i, iPar %i",iMt, iDetElemId/100, iDe, iDetElemId, iDetElemNumber, iDetElemNumber*fgNParCh+i));
- lDetElemMisAlignment[i] = misAlignments[iDetElemNumber*fgNParCh+i];
- }
- }
- // local transformation of this detection element.
- TGeoCombiTrans localTransform
- = TGeoCombiTrans(*detElement->GetLocalTransformation());
- TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
- newDetElement->SetLocalTransformation(newLocalTransform);
-
- // global transformation
- TGeoHMatrix newGlobalTransform =
- AliMUONGeometryBuilder::Multiply(newModuleTransform,
- newLocalTransform);
- newDetElement->SetGlobalTransformation(newGlobalTransform);
-
- // add this det element to module
- newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
- 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:
- // Tdelta = Tnew * Told.inverse
- TGeoHMatrix deltaGlobalTransform
- = AliMUONGeometryBuilder::Multiply(newGlobalTransform,
- newOldGlobalTransform.Inverse());
-
- // Create mis alignment matrix
- newGeometryTransformer
- ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
- }
+ // add number of detectors up to this chamber
+ return iDet + fgSNDetElemCh[iCh-1];
- if (verbose)
- AliInfo(Form("Added module transformer %i to the transformer", iMt));
- newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
- }
- return newGeometryTransformer;
}
//______________________________________________________________________
-void AliMUONAlignment::SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t rChId, Double_t rChResX, Double_t rChResY, Double_t rDeResX, Double_t rDeResY)
+Int_t AliMUONAlignment::GetChamberId( Int_t iDetElemNumber ) const
{
+ /// get chamber (counting from 1) matching a given detector element id
+ Int_t iCh( 0 );
+ for( iCh=0; iCh<fgNCh; iCh++ )
+ { if( iDetElemNumber < fgSNDetElemCh[iCh] ) break; }
- /// 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;
+ return iCh;
+}
- 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);
- }
+//______________________________________________________________________
+TString AliMUONAlignment::GetParameterMaskString( UInt_t mask ) const
+{
+ TString out;
+ if( mask & ParX ) out += "X";
+ if( mask & ParY ) out += "Y";
+ if( mask & ParZ ) out += "Z";
+ if( mask & ParTZ ) out += "T";
+ return out;
+}
- for (int i=0; i<misAlignArray->GetEntries(); 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);
- }
- }
- }
- }
+//______________________________________________________________________
+TString AliMUONAlignment::GetSidesMaskString( UInt_t mask ) const
+{
+ TString out;
+ if( mask & SideTop ) out += "T";
+ if( mask & SideLeft ) out += "L";
+ if( mask & SideBottom ) out += "B";
+ if( mask & SideRight ) out += "R";
+ return out;
}
/// \class AliMUONAlignmentTask
/// AliAnalysisTask to align the MUON spectrometer.
/// The Task reads as input ESDs and feeds the MUONTracks to AliMUONAlignment.
-/// The alignment itself is performed by AliMillepede.
+/// The alignment itself is performed by AliMillePede2.
/// A OCDB entry is written with the alignment parameters.
-/// A list of graph are written to monitor the alignment parameters.
///
/// \author Javier Castillo, CEA/Saclay - Irfu/SPhN
+/// \author Hugo Pereira Da Costa, CEA/Saclay - Irfu/SPhN
//-----------------------------------------------------------------------------
// this class header must always come first
// local header must come before system headers
#include "AliAnalysisManager.h"
+#include "AliAlignObjMatrix.h"
#include "AliAODInputHandler.h"
#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliAODHeader.h"
+#include "AliCDBEntry.h"
#include "AliESDInputHandler.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
-#include "AliMagF.h"
#include "AliCDBStorage.h"
#include "AliCDBManager.h"
#include "AliGRPManager.h"
#include "AliCDBId.h"
#include "AliGeomManager.h"
+#include "AliMagF.h"
+#include "AliMillePedeRecord.h"
#include "AliMpCDB.h"
#include "AliMUONCDB.h"
#include "AliMUONAlignment.h"
+#include "AliMUONConstants.h"
+#include "AliMUONRecoParam.h"
#include "AliMUONTrack.h"
#include "AliMUONTrackExtrap.h"
#include "AliMUONTrackParam.h"
#include "AliMUONESDInterface.h"
// system headers
+#include <cassert>
#include <fstream>
+
+// root headers
#include <TString.h>
#include <TError.h>
+#include <TFile.h>
#include <TGraphErrors.h>
#include <TTree.h>
#include <TChain.h>
///\endcond
//________________________________________________________________________
-AliMUONAlignmentTask::AliMUONAlignmentTask(const char *name, const char *newalignocdb, const char *oldalignocdb, const char *defaultocdb, const char *geofilename):
+AliMUONAlignmentTask::AliMUONAlignmentTask( const char *name ):
AliAnalysisTaskSE(name),
fReadRecords( kFALSE ),
fWriteRecords( kFALSE ),
fDoAlignment( kTRUE ),
+ fMergeAlignmentCDBs( kTRUE ),
+ fForceBField( kFALSE ),
+ fBFieldOn( kFALSE ),
+ fUnbias( kFALSE ),
fAlign(0x0),
- fGeoFilename(geofilename),
- fDefaultStorage(defaultocdb),
- fOldAlignStorage(oldalignocdb),
- fNewAlignStorage(newalignocdb),
- fOldGeoTransformer(NULL),
- fNewGeoTransformer(NULL),
+ fNewAlignStorage( "local://ReAlignOCDB" ),
+ fOldGeoTransformer(0x0),
+ fNewGeoTransformer(0x0),
fLoadOCDBOnce(kFALSE),
fOCDBLoaded(kFALSE),
+ fEvent(0),
fTrackTot(0),
fTrackOk(0),
- fLastRunNumber(-1),
- fMSDEx(0x0),
- fMSDEy(0x0),
- fMSDEz(0x0),
- fMSDEp(0x0),
- fList(0x0),
+ fRunNumberMin(0),
+ fRunNumberMax(0),
fRecords(0x0),
fRecordCount(0)
{
// Input slot #0 works with a TChain
DefineInput(0, TChain::Class());
- // Output slot #0 writes NTuple/histos into a TList
- DefineOutput(1, TList::Class());
-
// initialize parameters ...
- for(Int_t k=0;k<4*156;k++)
+ for(Int_t k=0;k<AliMUONAlignment::fNGlobal;k++)
{
fParameters[k]=0.;
fErrors[k]=0.;
fPulls[k]=0.;
}
- fOldGeoTransformer = new AliMUONGeometryTransformer();
-
-}
-
-//________________________________________________________________________
-AliMUONAlignmentTask::AliMUONAlignmentTask(const AliMUONAlignmentTask& other):
- AliAnalysisTaskSE(other),
- fReadRecords( other.fReadRecords ),
- fWriteRecords( other.fWriteRecords ),
- fDoAlignment( other.fDoAlignment ),
- fAlign( other.fAlign ),
- fGeoFilename( other.fGeoFilename ),
- fDefaultStorage( other.fDefaultStorage ),
- fOldAlignStorage( other.fOldAlignStorage ),
- fNewAlignStorage( other.fNewAlignStorage ),
- fOldGeoTransformer( other.fOldGeoTransformer ),
- fNewGeoTransformer( other.fNewGeoTransformer ),
- fLoadOCDBOnce( other.fLoadOCDBOnce ),
- fOCDBLoaded( other.fOCDBLoaded ),
- fTrackTot( other.fTrackTot ),
- fTrackOk( other.fTrackOk ),
- fLastRunNumber( other.fLastRunNumber ),
- fMSDEx( other.fMSDEx ),
- fMSDEy( other.fMSDEy ),
- fMSDEz( other.fMSDEz ),
- fMSDEp( other.fMSDEp ),
- fList( other.fList ),
- fRecords( other.fRecords ),
- fRecordCount( other.fRecordCount )
-{
-
- // initialize parameters
- for(Int_t k=0;k<4*156;k++)
- {
-
- fParameters[k]=other.fParameters[k];
- fErrors[k]=other.fErrors[k];
- fPulls[k]=other.fPulls[k];
-
- }
-
-}
-
-//________________________________________________________________________
-AliMUONAlignmentTask& AliMUONAlignmentTask::operator=(const AliMUONAlignmentTask& other)
-{
- /// Assignment
- if(&other == this) return *this;
- AliAnalysisTaskSE::operator=(other);
-
- fReadRecords = other.fReadRecords;
- fWriteRecords = other.fWriteRecords;
- fDoAlignment = other.fDoAlignment;
-
- // this breaks in destructor
- fAlign = other.fAlign;
-
- fGeoFilename = other.fGeoFilename;
- fDefaultStorage = other.fDefaultStorage;
- fOldAlignStorage = other.fOldAlignStorage;
- fNewAlignStorage = other.fNewAlignStorage;
-
- // this breaks in destructor
- fOldGeoTransformer = other.fOldGeoTransformer;
- fNewGeoTransformer = other.fNewGeoTransformer;
-
- fLoadOCDBOnce = other.fLoadOCDBOnce;
- fOCDBLoaded = other.fOCDBLoaded;
- fTrackTot = other.fTrackTot;
- fTrackOk = other.fTrackOk;
- fLastRunNumber = other.fLastRunNumber;
- fMSDEx = other.fMSDEx;
- fMSDEy = other.fMSDEy;
- fMSDEz = other.fMSDEz;
- fMSDEp = other.fMSDEp;
- fList = other.fList;
- fRecords = other.fRecords;
- fRecordCount = other.fRecordCount;
-
- // initialize parameters ...
- for( Int_t k=0; k<4*156; ++k)
- {
-
- fParameters[k]=other.fParameters[k];
- fErrors[k]=other.fErrors[k];
- fPulls[k]=other.fPulls[k];
+ // create alignment object
+ fAlign = new AliMUONAlignment();
- }
+ // create old geometry transformer
+ fOldGeoTransformer = new AliMUONGeometryTransformer();
- return *this;
}
//________________________________________________________________________
AliMUONAlignmentTask::~AliMUONAlignmentTask()
{
- /*
- it is safe to delete NULL pointers, so
- there is no need to test their validity here.
- However, it crashes here, because of incorrect assignment operator
- and copy constructor, resulting in double deletion of objects.
- Would require deep copy instead.
- */
+ /// destructor
delete fAlign;
delete fOldGeoTransformer;
delete fNewGeoTransformer;
//________________________________________________________________________
void AliMUONAlignmentTask::LocalInit()
{
- /// Local initialization, called once per task on the client machine
- /// where the analysis train is assembled
+ /**
+ Local initialization, called once per task on the client machine
+ where the analysis train is assembled
+ **/
+
+ /* must run alignment when reading records */
+ if( fReadRecords && !fDoAlignment )
+ {
+
+ AliInfo( "Automatically setting fDoAlignment to kTRUE because fReadRecords is kTRUE" );
+ SetDoAlignment( kTRUE );
+
+ }
// print configuration
+ if( fRunNumberMin > 0 || fRunNumberMax > 0 )
+ { AliInfo( Form( "run range: %i - %i", fRunNumberMin, fRunNumberMax ) ); }
+
AliInfo( Form( "fReadRecords: %s", (fReadRecords ? "kTRUE":"kFALSE" ) ) );
AliInfo( Form( "fWriteRecords: %s", (fWriteRecords ? "kTRUE":"kFALSE" ) ) );
AliInfo( Form( "fDoAlignment: %s", (fDoAlignment ? "kTRUE":"kFALSE" ) ) );
+ if( fDoAlignment )
+ {
+ // merge alignment DB flag is irrelevant if not actually performing the alignemnt
+ AliInfo( Form( "fMergeAlignmentCDBs: %s", (fMergeAlignmentCDBs ? "kTRUE":"kFALSE" ) ) );
+ }
+
+ // storage elements
+ if( !fDefaultStorage.IsNull() ) AliInfo( Form( "fDefaultStorage: %s", fDefaultStorage.Data() ) );
+ if( !fOldAlignStorage.IsNull() ) AliInfo( Form( "fOldAlignStorage: %s", fOldAlignStorage.Data() ) );
+
+ if( fDoAlignment )
+ {
+ // new alignment storage is irrelevant if not actually performing the alignemnt
+ if( !fNewAlignStorage.IsNull() ) AliInfo( Form( "fNewAlignStorage: %s", fNewAlignStorage.Data() ) );
+ else AliFatal( "Invalid new alignment storage path" );
+ }
+
+ if( !fReadRecords )
+ {
+ // following flags are only relevant if not reading records
+ if( fForceBField ) AliInfo( Form( "fBFieldOn: %s", (fBFieldOn ? "kTRUE":"kFALSE" ) ) );
+ else AliInfo( "fBFieldOn: from GRP" );
+ AliInfo( Form( "fUnbias: %s", (fUnbias ? "kTRUE":"kFALSE" ) ) );
+ }
+
// consistency checks between flags
/* need at least one of the flags set to true */
if( !( fReadRecords || fWriteRecords || fDoAlignment ) )
if( fReadRecords && fWriteRecords )
{ AliFatal( "Cannot set both fReadRecords and fWriteRecords to kTRUE at the same time" ); }
- /* must run alignment when reading records */
- if( fReadRecords && !fDoAlignment )
+ /*
+ fix detectors
+ warning, counting chambers from 1.
+ this must be done before calling the Init method
+ */
+ if( fDoAlignment )
{
- AliInfo( "Automatically setting fDoAlignment to kTRUE because fReadRecords is kTRUE" );
- SetDoAlignment( kTRUE );
+ fAlign->FixChamber(1);
+ fAlign->FixChamber(10);
- }
+ } else {
- // Set initial values here, good guess may help convergence
- // St 1
- // Int_t iPar = 0;
- // fParameters[iPar++] = 0.010300 ; fParameters[iPar++] = 0.010600 ; fParameters[iPar++] = 0.000396 ;
+ AliInfo( "Not fixing detector elements, since alignment is not required" );
- fAlign = new AliMUONAlignment();
- fAlign->InitGlobalParameters(fParameters);
-
-// AliCDBManager::Instance()->Print();
-//
-// fAlign->SetGeometryTransformer(fOldGeoTransformer);
+ }
- // Do alignment with magnetic field off
- fAlign->SetBFieldOn(kFALSE);
+ // initialize
+ fAlign->Init();
- // Set tracking station to use
- // Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
- Bool_t bChOnOff[10] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
+ // use unbiased residuals
+ fAlign->SetUnbias( fUnbias );
- // Set degrees of freedom to align (see AliMUONAlignment)
- fAlign->AllowVariations(bChOnOff);
+ // Do alignment with magnetic field off
+ fAlign->SetBFieldOn( fBFieldOn );
// Set expected resolution (see AliMUONAlignment)
fAlign->SetSigmaXY(0.15,0.10);
- // Fix parameters or add constraints here
- // for (Int_t iSt=0; iSt<5; iSt++)
- // if (!bStOnOff[iSt]) fAlign->FixStation(iSt+1);
- for (Int_t iCh=0; iCh<10; iCh++)
- if (!bChOnOff[iCh]) fAlign->FixChamber(iCh+1);
-
- // Left and right sides of the detector are independent, one can choose to align
- // only one side
- Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
- fAlign->FixHalfSpectrometer(bChOnOff,bSpecLROnOff);
-
- fAlign->SetChOnOff(bChOnOff);
- fAlign->SetSpecLROnOff(bChOnOff);
-
- // Here we can fix some detection elements
- // fAlign->FixDetElem(908);
- // fAlign->FixDetElem(1012);
- fAlign->FixDetElem(608);
-
- // Set predifined global constrains: X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
-// Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
-// Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
- // fAlign->AddConstraints(bChOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
+ // initialize global parameters to provided values
+ fAlign->InitGlobalParameters(fParameters);
}
void AliMUONAlignmentTask::UserCreateOutputObjects()
{
- if( fDoAlignment )
- {
-
- // Creating graphs
- fMSDEx = new TGraphErrors(156);
- fMSDEy = new TGraphErrors(156);
- fMSDEz = new TGraphErrors(156);
- fMSDEp = new TGraphErrors(156);
-
- // Add Ntuples to the list
- fList = new TList();
- fList->Add(fMSDEx);
- fList->Add(fMSDEy);
- fList->Add(fMSDEz);
- fList->Add(fMSDEp);
-
- fList->SetOwner(kTRUE);
- PostData(1, fList);
- }
-
// connect AOD output
if( fWriteRecords )
{
{
// get AOD output handler and add Branch
- fRecords = new TClonesArray( "AliMUONAlignmentTrackRecord", 0 );
+ fRecords = new TClonesArray( "AliMillePedeRecord", 0 );
fRecords->SetName( "records" );
handler->AddBranch("TClonesArray", &fRecords);
fRecordCount = 0;
- } else AliInfo( "Error: invalid output event handler" );
+ } else AliFatal( "Error: invalid output event handler" );
}
//________________________________________________________________________
void AliMUONAlignmentTask::UserExec(Option_t *)
{
- // print configuration
-// AliInfo( Form( "fReadRecords: %s", (fReadRecords ? "kTRUE":"kFALSE" ) ) );
-// AliInfo( Form( "fWriteRecords: %s", (fWriteRecords ? "kTRUE":"kFALSE" ) ) );
-// AliInfo( Form( "fDoAlignment: %s", (fDoAlignment ? "kTRUE":"kFALSE" ) ) );
+
+ // do nothing if run number not in range
+ if( fRunNumberMin > 0 && fCurrentRunNumber < fRunNumberMin ) return;
+ if( fRunNumberMax > 0 && fCurrentRunNumber > fRunNumberMax ) return;
+
+ // increase event count
+ ++fEvent;
+
// clear array
if( fWriteRecords && fRecords )
{
fRecordCount = 0;
}
- // local track parameters
- Double_t trackParams[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
-
- if( fReadRecords ) {
+ if( fReadRecords )
+ {
AliAODEvent* lAOD( dynamic_cast<AliAODEvent*>(InputEvent() ) );
for( Int_t index = 0; index < lRecordsCount; ++index )
{
- if( AliMUONAlignmentTrackRecord* record = dynamic_cast<AliMUONAlignmentTrackRecord*>( records->UncheckedAt(index) ) )
+ if( AliMillePedeRecord* record = dynamic_cast<AliMillePedeRecord*>( records->UncheckedAt(index) ) )
{
- fAlign->ProcessTrack( record, fDoAlignment );
- if( fDoAlignment )
- { fAlign->LocalFit( fTrackOk++, trackParams, 0 ); }
+ fAlign->ProcessTrack( record );
+ ++fTrackOk;
+
+ if(!(fTrackTot%100))
+ { AliInfo( Form( "Processed %i Tracks and %i were fitted.", fTrackTot, fTrackOk ) ); }
} else AliInfo( Form( "Invalid record at %i", index ) );
return;
}
-// HUGO: Comment out check on run number, to be able to run on MC
-// if (!lESD->GetRunNumber())
-// {
-// AliInfo( Form( "Current Run Number: %i", lESD->GetRunNumber() ) );
-// return;
-// }
-
- // if (lESD->GetRunNumber()!=fLastRunNumber){
- // fLastRunNumber = lESD->GetRunNumber();
- // Prepare(fGeoFilename.Data(),fDefaultOCDB.Data(),fMisAlignOCDB.Data());
- // }
-
Int_t nTracks = Int_t(lESD->GetNumberOfMuonTracks());
-// cout << " there are " << nTracks << " tracks" << endl;
for( Int_t iTrack = 0; iTrack < nTracks; iTrack++ )
{
if (!esdTrack->ContainTriggerData()) continue;
Double_t invBenMom = esdTrack->GetInverseBendingMomentum();
- // fInvBenMom->Fill(invBenMom);
- // fBenMom->Fill(1./invBenMom);
if (TMath::Abs(invBenMom)<=1.04)
{
AliMUONESDInterface::ESDToMUON(*esdTrack, track);
// process track and retrieve corresponding records, for storage
- const AliMUONAlignmentTrackRecord* lRecords( fAlign->ProcessTrack( &track, fDoAlignment ) );
-
- // do the fit, if required
- if( fDoAlignment ) fAlign->LocalFit(fTrackOk++,trackParams,0);
- else fTrackOk++;
+ const AliMillePedeRecord* lRecords( fAlign->ProcessTrack( &track, fDoAlignment ) );
+ ++fTrackOk;
// store in array
if( fWriteRecords && fRecords )
- {
- new((*fRecords)[fRecordCount]) AliMUONAlignmentTrackRecord( *lRecords );
+ {
+ new((*fRecords)[fRecordCount]) AliMillePedeRecord( *lRecords );
++fRecordCount;
}
- }
-
+ }
+
++fTrackTot;
- if(!(fTrackTot%1000))
+ if(!(fTrackTot%100))
{ AliInfo( Form( "Processed %i Tracks and %i were fitted.", fTrackTot, fTrackOk ) ); }
-// cout << "Processed " << fTrackTot << " Tracks." << endl;
// Post final data. Write histo list to a file with option "RECREATE"
- PostData(1,fList);
+ // PostData(1,fList);
}
// save AOD
- if( fWriteRecords && fRecordCount > 0 ) {
+ if( fWriteRecords && fRecordCount > 0 )
+ {
AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
if( handler )
- {
-// printf("handler: %p\n",handler);
- AliAODEvent* aod = handler->GetAOD();
-// printf("aod: %p\n",aod);
- AliAODHeader* header = aod->GetHeader();
-// printf("header: %p\n",header);
- header->SetRunNumber(lESD->GetRunNumber());
-// printf("RunNumber: %d\n",lESD->GetRunNumber());
- AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
- } else AliInfo( "Error: invalid output event handler" );
+ {
+ AliAODEvent* aod = handler->GetAOD();
+ AliAODHeader* header = aod->GetHeader();
+ header->SetRunNumber(lESD->GetRunNumber());
+ AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
+
+ } else AliInfo( "Error: invalid output event handler" );
+
}
}
}
-//________________________________________________________________________
-void AliMUONAlignmentTask::Terminate(const Option_t*)
-{ return; }
-
//________________________________________________________________________
void AliMUONAlignmentTask::FinishTaskOutput()
{
// Perform global fit
fAlign->GlobalFit(fParameters,fErrors,fPulls);
-
-// // Update pointers reading them from the output slot
-// fList = (TList*)GetOutputData(0);
-// fMSDEx = (TGraphErrors*)fList->At(0);
-// fMSDEy = (TGraphErrors*)fList->At(1);
-// fMSDEz = (TGraphErrors*)fList->At(2);
-// fMSDEp = (TGraphErrors*)fList->At(3);
-
- // Store results
- Double_t DEid[156] = {0};
- Double_t MSDEx[156] = {0};
- Double_t MSDEy[156] = {0};
- Double_t MSDEz[156] = {0};
- Double_t MSDEp[156] = {0};
- Double_t DEidErr[156] = {0};
- Double_t MSDExErr[156] = {0};
- Double_t MSDEyErr[156] = {0};
- Double_t MSDEzErr[156] = {0};
- Double_t MSDEpErr[156] = {0};
- Int_t lNDetElem = 4*2+4*2+18*2+26*2+26*2;
- Int_t lNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
- // Int_t lSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
Int_t idOffset = 0; // 400
Int_t lSDetElemCh = 0;
- for(Int_t iDE=0; iDE<lNDetElem; iDE++)
+
+ for( Int_t iDE=0; iDE<AliMUONAlignment::fgNDetElem; iDE++ )
{
- DEidErr[iDE] = 0.;
- DEid[iDE] = idOffset+100;
- DEid[iDE] += iDE;
+ Int_t detElemId = idOffset+100;
+ detElemId += iDE;
lSDetElemCh = 0;
for(Int_t iCh=0; iCh<9; iCh++)
{
- lSDetElemCh += lNDetElemCh[iCh];
+ lSDetElemCh += AliMUONAlignment::fgNDetElemCh[iCh];
if (iDE>=lSDetElemCh)
{
- DEid[iDE] += 100;
- DEid[iDE] -= lNDetElemCh[iCh];
+ detElemId += 100;
+ detElemId -= AliMUONAlignment::fgNDetElemCh[iCh];
}
}
- MSDEx[iDE]=fParameters[4*iDE+0];
- MSDEy[iDE]=fParameters[4*iDE+1];
- MSDEz[iDE]=fParameters[4*iDE+3];
- MSDEp[iDE]=fParameters[4*iDE+2];
- MSDExErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+0);
- MSDEyErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+1);
- MSDEzErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+3);
- MSDEpErr[iDE]=(Double_t)fAlign->GetParError(4*iDE+2);
- fMSDEx->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+0]);
- fMSDEx->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+0));
- fMSDEy->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+1]);
- fMSDEy->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+1));
- fMSDEz->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+3]);
- fMSDEz->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+3));
- fMSDEp->SetPoint(iDE,DEid[iDE],fParameters[4*iDE+2]);
- fMSDEp->SetPointError(iDE,DEidErr[iDE],(Double_t)fAlign->GetParError(4*iDE+2));
-
}
// Post final data. Write histo list to a file with option "RECREATE"
- PostData(1,fList);
+ // PostData(1,fList);
- // HUGO: stop here to test reproducibility
- // return;
+ // store misalignments from OCDB into old transformers
+ if( fMergeAlignmentCDBs )
+ { SaveMisAlignmentData( fOldGeoTransformer ); }
// Re Align
- fNewGeoTransformer = fAlign->ReAlign(fOldGeoTransformer,fParameters,true);
- // newTransform->WriteTransformations("transform2ReAlign.dat");
+ fNewGeoTransformer = fAlign->ReAlign( fOldGeoTransformer, fParameters, true );
// Generate realigned data in local cdb
const TClonesArray* array = fNewGeoTransformer->GetMisAlignmentData();
// 100 mum residual resolution for chamber misalignments?
- fAlign->SetAlignmentResolution(array,-1,0.01,0.01,0.004,0.003);
+ fAlign->SetAlignmentResolution( array, -1, 0.01, 0.01, 0.004, 0.003 );
// CDB manager
AliLog::SetGlobalDebugLevel(2);
- AliCDBManager* cdbm = AliCDBManager::Instance();
- // cdbManager->SetDefaultStorage(fDefaultOCDB.Data());
- // recover default storage full name (raw:// cannot be used to set specific storage)
- TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
- if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
- else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
+ // recover default storage full name (raw:// cannot be used to set specific storage)
+ AliCDBManager* cdbManager = AliCDBManager::Instance();
+
+ // unload old alignment path
+ if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
+ { cdbManager->UnloadFromCache("MUON/Align/Data"); }
- if (fOldAlignStorage != "none") cdbm->UnloadFromCache("MUON/Align/Data");
- if (!fNewAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
- // else cdbm->SetSpecificStorage("MUON/Align/Data",fDefaultStorage.Data());
- else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+ // load new alignment path
+ if( !fNewAlignStorage.IsNull() ) cdbManager->SetSpecificStorage("MUON/Align/Data",fNewAlignStorage.Data());
+ else {
+ TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
+ if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+ else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+ cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+
+ }
+ // create new DB entry
AliCDBMetaData* cdbData = new AliCDBMetaData();
cdbData->SetResponsible("Dimuon Offline project");
cdbData->SetComment("MUON alignment objects with residual misalignment");
AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
- cdbm->Put(const_cast<TClonesArray*>(array), id, cdbData);
-
- gSystem->Exec("cp -a ReAlignOCDB/MUON/Align/Data/Run0_999999999_v0_s0.root Run0_999999999_v0_s0.root");
- gSystem->Exec("ls -l");
+ cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
}
//________________________________________________________________________
void AliMUONAlignmentTask::NotifyRun()
{
-
- if (fOCDBLoaded && fLoadOCDBOnce) { AliError(Form("OCDB already loaded %d %d",fOCDBLoaded,fLoadOCDBOnce));
+
+ /// run number (re)initialization
+
+ // propagate run number to fAlign
+ if( fAlign ) fAlign->SetRunNumber( fCurrentRunNumber );
+
+ // update ocdb
+ if (fOCDBLoaded && fLoadOCDBOnce)
+ {
+ AliError(Form("OCDB already loaded %d %d",fOCDBLoaded,fLoadOCDBOnce));
return;
}
- AliCDBManager* cdbm = AliCDBManager::Instance();
- cdbm->SetDefaultStorage(fDefaultStorage.Data());
-// cdbm->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+ AliCDBManager* cdbManager = AliCDBManager::Instance();
+
+ // Initialize default storage
+ if( !( cdbManager->IsDefaultStorageSet() || fDefaultStorage.IsNull() ) )
+ {
+
+ AliInfo( Form( "Initializing default storage: %s", fDefaultStorage.Data() ) );
+ cdbManager->SetDefaultStorage(fDefaultStorage.Data());
+
+ } else if( !fDefaultStorage.IsNull() ) {
+
+ AliInfo( "Default storage already set. Ignoring fDefaultStorage" );
+
+ } else {
+
+ AliInfo( "Default storage already set" );
+
+ }
+
+ // Initialize run number
+ if( cdbManager->GetRun() < 0 )
+ {
+ AliInfo( Form( "Setting run number: %i", fCurrentRunNumber ) );
+ cdbManager->SetRun(fCurrentRunNumber);
+ }
+
+ // following initialization is not needed when reading records
+ if( !fReadRecords )
+ {
+
+ // load magnetic field if needed
+ if( !TGeoGlobalMagField::Instance()->IsLocked() )
+ {
- // HUGO: add to comment out the specific settings below in order to be able to run.
- //cdbm->SetSpecificStorage("GRP/GRP/Data",fDefaultStorage.Data());
- //cdbm->SetSpecificStorage("GRP/Geometry/data",fDefaultStorage.Data());
- //cdbm->SetSpecificStorage("MUON/Calib/MappingData",fDefaultStorage.Data());
- cdbm->SetRun(fCurrentRunNumber);
+ AliInfo( "Loading magnetic field" );
+ if( !AliMUONCDB::LoadField() )
+ {
+ AliError( "Failed to load magnetic field" );
+ return;
+ }
- if (!AliMUONCDB::LoadField()) { AliError("Problem field"); return;}
+ } else { AliInfo( "Magnetic field already locked" ); }
- // set the magnetic field for track extrapolations
- AliMUONTrackExtrap::SetField();
+ // checking magnitic field
+ if( !fForceBField )
+ {
+ AliInfo( "Reading magnetic field setup from GRP" );
+
+ // decide bFieldOn value base on dipole current
+ // propagete to Alignment class
+ // and printout
+ AliMagF* magF = dynamic_cast<AliMagF*>( TGeoGlobalMagField::Instance()->GetField() );
+ fBFieldOn = TMath::Abs( magF->GetFactorDip() ) > 1e-5;
+ fAlign->SetBFieldOn( fBFieldOn );
+ AliInfo( Form( "Dipole magnetic field factor: %.2f", magF->GetFactorDip() ) );
+ AliInfo( Form( "fBFieldOn = %s", (fBFieldOn ? "kTRUE":"kFALSE") ) );
+ }
- if (!AliMUONCDB::LoadMapping(kTRUE)) { AliError("Problem mapping"); return;}
+ AliInfo( "Loading muon mapping" );
+ if( !AliMUONCDB::LoadMapping(kFALSE) )
+ {
+ AliError( "Failed to load muon mapping" );
+ return;
+ }
- // recover default storage full name (raw:// cannot be used to set specific storage)
- TString defaultStorage(cdbm->GetDefaultStorage()->GetType());
- if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
- else defaultStorage += Form("://%s", cdbm->GetDefaultStorage()->GetBaseFolder().Data());
+ AliInfo( "Assigning field to Track extrapolator" );
+ AliMUONTrackExtrap::SetField();
- // reset existing geometry/alignment if any
- if (cdbm->GetEntryCache()->Contains("GRP/Geometry/Data")) cdbm->UnloadFromCache("GRP/Geometry/Data");
- if (cdbm->GetEntryCache()->Contains("MUON/Align/Data")) cdbm->UnloadFromCache("MUON/Align/Data");
- if (AliGeomManager::GetGeometry()) AliGeomManager::GetGeometry()->UnlockGeometry();
+ }
- // get original geometry transformer
- AliGeomManager::LoadGeometry();
- if (!AliGeomManager::GetGeometry()) { AliError("Problem geometry"); return;}
- if (fOldAlignStorage != "none")
+ // load geometry if needed
+ if( !AliGeomManager::GetGeometry() )
{
- if (!fOldAlignStorage.IsNull()) cdbm->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
- else cdbm->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
- // else cdbm->SetSpecificStorage("MUON/Align/Data",fDefaultStorage.Data());
-
+ // reset existing geometry/alignment if any
+ if( cdbManager->GetEntryCache()->Contains("GRP/Geometry/Data") )
+ {
+ AliInfo( "Unloading GRP/Geometry/Data" );
+ cdbManager->UnloadFromCache("GRP/Geometry/Data");
+ }
+
+ if( cdbManager->GetEntryCache()->Contains("MUON/Align/Data") )
+ {
+ AliInfo( Form( "Unloading MUON/Align/Data from %s", cdbManager->GetSpecificStorage( "MUON/Align/Data" )->GetBaseFolder().Data() ) );
+ cdbManager->UnloadFromCache("MUON/Align/Data");
+ }
+
+ // get original geometry transformer
+ AliInfo( "Loading geometry" );
+ AliGeomManager::LoadGeometry();
+ if (!AliGeomManager::GetGeometry())
+ {
+ AliError("Failed to load geometry");
+ return;
+ }
+
+ if (!fOldAlignStorage.IsNull())
+ {
+
+ AliInfo( Form( "Initializing MUON/Align/Data using: %s", fOldAlignStorage.Data() ) );
+ cdbManager->SetSpecificStorage("MUON/Align/Data",fOldAlignStorage.Data());
+
+ } else {
+
+ // recover default storage full name (raw:// cannot be used to set specific storage)
+ TString defaultStorage(cdbManager->GetDefaultStorage()->GetType());
+ if (defaultStorage == "alien") defaultStorage += Form("://folder=%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+ else defaultStorage += Form("://%s", cdbManager->GetDefaultStorage()->GetBaseFolder().Data());
+
+ AliInfo( Form( "Re-initializing MUON/Align/Data using: %s", defaultStorage.Data() ) );
+ cdbManager->SetSpecificStorage("MUON/Align/Data",defaultStorage.Data());
+
+ }
+
+ AliInfo("Loading muon Alignment objects");
AliGeomManager::ApplyAlignObjsFromCDB("MUON");
+ } else { AliInfo( "Geometry already initialized - fOldAlignStorage ignored" ); }
+
+ // update geometry transformer and pass to Alignment object
+ AliInfo("Loading muon geometry in transformer");
+ fOldGeoTransformer->LoadGeometryData();
+ fAlign->SetGeometryTransformer(fOldGeoTransformer);
+
+ fOCDBLoaded = kTRUE;
+
+}
+
+//_____________________________________________________________________________________
+void AliMUONAlignmentTask::SaveMisAlignmentData( AliMUONGeometryTransformer* transformer ) const
+{
+
+ // clear transformer
+ transformer->ClearMisAlignmentData();
+
+ // load MUON/Align/Data from OCDB
+ AliCDBManager* cdbManager = AliCDBManager::Instance();
+ AliCDBEntry* cdbEntry = cdbManager->Get( "MUON/Align/Data" );
+ if (!cdbEntry)
+ {
+
+ AliError( "unable to load entry for path MUON/Align/Data" );
+ return;
+
+ }
+
+ // get TClonesArray and check
+ TClonesArray* misArray = (TClonesArray*)cdbEntry->GetObject();
+ if( !misArray )
+ {
+
+ AliError( "unable to load old misalignment array" );
+ return;
+
}
- // fOldGeoTransformer = new AliMUONGeometryTransformer();
- fOldGeoTransformer->LoadGeometryData();
- fAlign->SetGeometryTransformer(fOldGeoTransformer);
+ // loop over stored entries
+ for (Int_t index=0; index<misArray->GetEntriesFast(); ++index )
+ {
+
+ // load matrix and check
+ AliAlignObjMatrix* matrix = (AliAlignObjMatrix*) misArray->At( index );
+ if( !matrix )
+ {
+ AliError( Form( "unable to load matrix for index %i", index ) );
+ continue;
+ }
+
+ // get volume ID
+ AliGeomManager::ELayerID layerId;
+ Int_t moduleId;
+ matrix->GetVolUID( layerId, moduleId);
+
+ // make sure ELayerID is correct
+ assert( layerId == AliGeomManager::kMUON );
- fOCDBLoaded = kTRUE;
+ // printout
+ // AliInfo( Form( "Found matrix for %s %i", matrix->GetSymName(), moduleId ) );
+ // get matrix
+ TGeoHMatrix misMatrix;
+ matrix->GetMatrix(misMatrix);
+
+ // add to geometry transformer
+ // need the detector element
+ // "detElement->GetId()"
+ // see fOldGeoTransformer->GetMisAlignmentData( ... )
+
+ if( TString( matrix->GetSymName() ).Contains( "DE" ) ) transformer->AddMisAlignDetElement( moduleId, misMatrix);
+ else transformer->AddMisAlignModule( moduleId, misMatrix);
+ }
+
+ return;
}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff