// ---------------------
// Class AliMUONTrackK
// ---------------------
-// Reconstructed track in the muons system based on the extended
+// Reconstructed track in the muon system based on the extended
// Kalman filter approach
// Author: Alexander Zinchenko, JINR Dubna
-#include <Riostream.h>
-#include <TClonesArray.h>
-#include <TArrayD.h>
-#include <TMatrixD.h>
-
#include "AliMUONTrackK.h"
-#include "AliMUON.h"
+#include "AliMUONRecData.h"
+#include "AliMUONConstants.h"
-#include "AliMUONTrackReconstructor.h"
-#include "AliMagF.h"
-#include "AliMUONSegment.h"
+#include "AliMUONTrackReconstructorK.h"
#include "AliMUONHitForRec.h"
+#include "AliMUONObjectPair.h"
#include "AliMUONRawCluster.h"
#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
#include "AliMUONEventRecoCombi.h"
#include "AliMUONDetElement.h"
-#include "AliRun.h"
+
#include "AliLog.h"
+#include <Riostream.h>
+#include <TClonesArray.h>
+#include <TArrayD.h>
+
+/// \cond CLASSIMP
+ClassImp(AliMUONTrackK) // Class implementation in ROOT context
+/// \endcond
+
const Int_t AliMUONTrackK::fgkSize = 5;
const Int_t AliMUONTrackK::fgkNSigma = 12;
const Double_t AliMUONTrackK::fgkChi2max = 100;
const Int_t AliMUONTrackK::fgkTriesMax = 10000;
const Double_t AliMUONTrackK::fgkEpsilon = 0.002;
-void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); // from AliMUONTrack
+void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
- Int_t AliMUONTrackK::fgDebug = -1; //-1;
+Int_t AliMUONTrackK::fgDebug = -1; //-1;
Int_t AliMUONTrackK::fgNOfPoints = 0;
-AliMUONTrackReconstructor* AliMUONTrackK::fgTrackReconstructor = NULL;
+AliMUONTrackReconstructorK* AliMUONTrackK::fgTrackReconstructor = NULL;
TClonesArray* AliMUONTrackK::fgHitForRec = NULL;
AliMUONEventRecoCombi *AliMUONTrackK::fgCombi = NULL;
-ClassImp(AliMUONTrackK) // Class implementation in ROOT context
-
//__________________________________________________________________________
AliMUONTrackK::AliMUONTrackK()
- //AZ: TObject()
- : AliMUONTrack()
+ : AliMUONTrack(),
+ fStartSegment(0x0),
+ fPosition(0.),
+ fPositionNew(0.),
+ fChi2(0.),
+ fTrackHits(0x0),
+ fNmbTrackHits(0),
+ fTrackDir(1),
+ fBPFlag(kFALSE),
+ fRecover(0),
+ fSkipHit(0x0),
+ fTrackPar(0x0),
+ fTrackParNew(0x0),
+ fCovariance(0x0),
+ fWeight(0x0),
+ fParExtrap(0x0),
+ fParFilter(0x0),
+ fParSmooth(0x0),
+ fCovExtrap(0x0),
+ fCovFilter(0x0),
+ fJacob(0x0),
+ fNSteps(0),
+ fSteps(0x0),
+ fChi2Array(0x0),
+ fChi2Smooth(0x0)
{
/// Default constructor
fgHitForRec = NULL; // pointer to points
fgNOfPoints = 0; // number of points
- fStartSegment = NULL;
- fTrackHits = NULL;
- fNmbTrackHits = 0;
- fTrackPar = fTrackParNew = NULL;
- fCovariance = fWeight = NULL;
- fSkipHit = NULL;
- fParExtrap = fParFilter = fParSmooth = NULL;
- fCovExtrap = fCovFilter = NULL;
- fJacob = NULL;
- fSteps = NULL; fNSteps = 0;
- fChi2Array = NULL;
- fChi2Smooth = NULL;
-
return;
}
//__________________________________________________________________________
-AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClonesArray *hitForRec)
- : AliMUONTrack()
+AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructorK *TrackReconstructor, TClonesArray *hitForRec)
+ : AliMUONTrack(),
+ fStartSegment(0x0),
+ fPosition(0.),
+ fPositionNew(0.),
+ fChi2(0.),
+ fTrackHits(0x0),
+ fNmbTrackHits(0),
+ fTrackDir(1),
+ fBPFlag(kFALSE),
+ fRecover(0),
+ fSkipHit(0x0),
+ fTrackPar(0x0),
+ fTrackParNew(0x0),
+ fCovariance(0x0),
+ fWeight(0x0),
+ fParExtrap(0x0),
+ fParFilter(0x0),
+ fParSmooth(0x0),
+ fCovExtrap(0x0),
+ fCovFilter(0x0),
+ fJacob(0x0),
+ fNSteps(0),
+ fSteps(0x0),
+ fChi2Array(0x0),
+ fChi2Smooth(0x0)
{
/// Constructor
fgNOfPoints = fgHitForRec->GetEntriesFast(); // number of points
fgCombi = NULL;
if (fgTrackReconstructor->GetTrackMethod() == 3) fgCombi = AliMUONEventRecoCombi::Instance();
-
- fStartSegment = NULL;
- fTrackHits = NULL;
- fNmbTrackHits = 0;
- fChi2 = 0;
- fTrackPar = fTrackParNew = NULL;
- fCovariance = fWeight = NULL;
- fSkipHit = NULL;
- fParExtrap = fParFilter = fParSmooth = NULL;
- fCovExtrap = fCovFilter = NULL;
- fJacob = NULL;
- fSteps = NULL; fNSteps = 0;
- fChi2Array = NULL;
- fChi2Smooth = NULL;
}
//__________________________________________________________________________
-AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
- //: AliMUONTrack(segment, segment, fgTrackReconstructor)
- : AliMUONTrack(NULL, segment, fgTrackReconstructor)
+AliMUONTrackK::AliMUONTrackK(AliMUONObjectPair *segment)
+ : AliMUONTrack(NULL,NULL),
+ fStartSegment(new AliMUONObjectPair(*segment)),
+ fPosition(0.),
+ fPositionNew(0.),
+ fChi2(0.),
+ fTrackHits(new TObjArray(13)),
+ fNmbTrackHits(2),
+ fTrackDir(1),
+ fBPFlag(kFALSE),
+ fRecover(0),
+ fSkipHit(0x0),
+ fTrackPar(new TMatrixD(fgkSize,1)),
+ fTrackParNew(new TMatrixD(fgkSize,1)),
+ fCovariance(new TMatrixD(fgkSize,fgkSize)),
+ fWeight(new TMatrixD(fgkSize,fgkSize)),
+ fParExtrap(new TObjArray(15)),
+ fParFilter(new TObjArray(15)),
+ fParSmooth(0x0),
+ fCovExtrap(new TObjArray(15)),
+ fCovFilter(new TObjArray(15)),
+ fJacob(new TObjArray(15)),
+ fNSteps(0),
+ fSteps(new TArrayD(15)),
+ fChi2Array(new TArrayD(13)),
+ fChi2Smooth(0x0)
{
/// Constructor from a segment
- Double_t dX, dY, dZ;
+ Double_t dX, dY, dZ, bendingSlope, bendingImpact;
AliMUONHitForRec *hit1, *hit2;
AliMUONRawCluster *clus;
TClonesArray *rawclusters;
- fStartSegment = segment;
- fRecover = 0;
- fSkipHit = NULL;
// Pointers to hits from the segment
- hit1 = segment->GetHitForRec1();
- hit2 = segment->GetHitForRec2();
+ hit1 = (AliMUONHitForRec*) segment->First();
+ hit2 = (AliMUONHitForRec*) segment->Second();
hit1->SetNTrackHits(hit1->GetNTrackHits()+1); // mark hit as being on track
hit2->SetNTrackHits(hit2->GetNTrackHits()+1); // mark hit as being on track
// check sorting in Z
if (TMath::Abs(hit1->GetZ()) > TMath::Abs(hit2->GetZ())) {
hit1 = hit2;
- hit2 = segment->GetHitForRec1();
+ hit2 = (AliMUONHitForRec*) segment->First();
}
- // memory allocation for the TObjArray of pointers to reconstructed TrackHit's
- fTrackHits = new TObjArray(13);
- fNmbTrackHits = 2;
- fChi2 = 0;
- fBPFlag = kFALSE;
- fTrackPar = new TMatrixD(fgkSize,1); // track parameters
- fTrackParNew = new TMatrixD(fgkSize,1); // track parameters
- fCovariance = new TMatrixD(fgkSize,fgkSize); // covariance matrix
- fWeight = new TMatrixD(fgkSize,fgkSize); // weight matrix (inverse of covariance)
// Fill array of track parameters
if (hit1->GetChamberNumber() > 7) {
dZ = hit2->GetZ() - hit1->GetZ();
dY = hit2->GetBendingCoor() - hit1->GetBendingCoor();
dX = hit2->GetNonBendingCoor() - hit1->GetNonBendingCoor();
+ bendingSlope = (hit2->GetBendingCoor() - hit1->GetBendingCoor()) / dZ;
+ bendingImpact = hit1->GetBendingCoor() - hit1->GetZ() * bendingSlope;
(*fTrackPar)(2,0) = TMath::ATan2(dY,dZ); // alpha
if ((*fTrackPar)(2,0) < 0.) (*fTrackPar)(2,0) += 2*TMath::Pi(); // from 0 to 2*pi
(*fTrackPar)(3,0) = TMath::ATan2(-dX,dZ/TMath::Cos((*fTrackPar)(2,0))); // beta
(*fTrackPar)(2,0) -= TMath::Pi();
- (*fTrackPar)(4,0) = 1/fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()); // 1/Pt
+ (*fTrackPar)(4,0) = 1./AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact); // 1/Pt
(*fTrackPar)(4,0) *= TMath::Cos((*fTrackPar)(3,0)); // 1/p
// Evaluate covariance (and weight) matrix
EvalCovariance(dZ);
- // For smoother
- fParExtrap = new TObjArray(15);
- fParFilter = new TObjArray(15);
- fCovExtrap = new TObjArray(15);
- fCovFilter = new TObjArray(15);
- fJacob = new TObjArray(15);
- fSteps = new TArrayD(15);
- fNSteps = 0;
- fChi2Array = new TArrayD(13);
- fChi2Smooth = NULL;
- fParSmooth = NULL;
-
if (fgDebug < 0 ) return;
- cout << fgTrackReconstructor->GetNRecTracks()-1 << " " << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
- if (fgTrackReconstructor->GetRecTrackRefHits()) {
- // from track ref. hits
- cout << ((AliMUONHitForRec*)((*fTrackHits)[0]))->GetTTRTrack() << "<-->" << ((AliMUONHitForRec*)((*fTrackHits)[1]))->GetTTRTrack() << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
- } else {
+ cout << fgTrackReconstructor->GetNRecTracks()-1 << " "
+ << AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact)
+ << " " << 1/(*fTrackPar)(4,0) << " ";
// from raw clusters
- for (Int_t i=0; i<2; i++) {
- hit1 = (AliMUONHitForRec*) ((*fTrackHits)[i]);
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit1->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit1->GetHitNumber());
- cout << clus->GetTrack(1);
- if (clus->GetTrack(2) != -1) cout << " " << clus->GetTrack(2);
- if (i == 0) cout << " <--> ";
- }
- cout << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
- }
+ for (Int_t i=0; i<2; i++) {
+ hit1 = (AliMUONHitForRec*) ((*fTrackHits)[i]);
+ rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit1->GetChamberNumber());
+ clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit1->GetHitNumber());
+ cout << clus->GetTrack(1);
+ if (clus->GetTrack(2) != -1) cout << " " << clus->GetTrack(2);
+ if (i == 0) cout << " <--> ";
+ }
+ cout << " @ " << hit1->GetChamberNumber() << endl;
+
}
//__________________________________________________________________________
{
/// Destructor
+ if (fStartSegment) {
+ delete fStartSegment;
+ fStartSegment = 0x0;
+ }
+
if (fTrackHits) {
//cout << fNmbTrackHits << endl;
for (Int_t i = 0; i < fNmbTrackHits; i++) {
delete fJacob;
}
- //__________________________________________________________________________
-AliMUONTrackK::AliMUONTrackK (const AliMUONTrackK& source)
- : AliMUONTrack(source)
-{
-/// Protected copy constructor
- AliFatal("Not implemented.");
-}
-
//__________________________________________________________________________
AliMUONTrackK & AliMUONTrackK::operator=(const AliMUONTrackK& source)
{
// base class assignement
//AZ TObject::operator=(source);
AliMUONTrack::operator=(source);
-
- fStartSegment = source.fStartSegment;
+
+ if (fStartSegment) delete fStartSegment;
+ if (source.fStartSegment) fStartSegment = new AliMUONObjectPair(*(source.fStartSegment));
+ else fStartSegment = 0x0;
+
fNmbTrackHits = source.fNmbTrackHits;
fChi2 = source.fChi2;
fPosition = source.fPosition;
} else if (fRecover) {
hit = GetHitLastOk();
currIndx = fTrackHits->IndexOf(hit);
- if (currIndx < 0) hit = fStartSegment->GetHitForRec1(); // for station 3
+ if (currIndx < 0) hit = (AliMUONHitForRec*) fStartSegment->First(); // for station 3
Back = kTRUE;
ichamb = hit->GetChamberNumber();
if (hit == fSkipHit || fRecover == 2 && currIndx >= 0) {
printf(" * %d %10.4f %10.4f %10.4f",
hit->GetChamberNumber(), hit->GetBendingCoor(),
hit->GetNonBendingCoor(), hit->GetZ());
- if (fgTrackReconstructor->GetRecTrackRefHits()) {
- // from track ref. hits
- printf(" %3d %3d \n", hit->GetTrackRefSignal(), hit->GetTTRTrack());
- } else {
- // from raw clusters
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- printf("%3d", clus->GetTrack(1)-1);
- if (clus->GetTrack(2) != 0) printf("%3d \n", clus->GetTrack(2)-1);
- else printf("\n");
- }
+ // from raw clusters
+ rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
+ clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
+ printf("%3d", clus->GetTrack(1)-1);
+ if (clus->GetTrack(2) != 0)
+ printf("%3d \n", clus->GetTrack(2)-1);
+ else
+ printf("\n");
+
}
} // if (fgDebug >= 10)
if (fNmbTrackHits>2 && fRecover==0) Recover(); // try to recover track later
} // while
if (fgDebug > 0) cout << fNmbTrackHits << " " << fChi2 << " " << 1/(*fTrackPar)(4,0) << " " << fPosition << endl;
if (1/TMath::Abs((*fTrackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) success = kFALSE; // p < p_min
+ if (GetRecover() < 0) success = kFALSE;
return success;
}
do {
step = TMath::Abs(step);
// Propagate parameters
- trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
+ AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
//extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
distance = zEnd - vGeant3New[2];
step *= dZ/(vGeant3New[2]-fPositionNew);
do {
// binary search
// Propagate parameters
- trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
+ AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
//extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
distance = zEnd - vGeant3New[2];
step /= 2;
tmp->SetUniqueID(1);
if (fSteps->fN <= fNSteps) fSteps->Set(fSteps->fN+10);
fSteps->AddAt(fPositionNew,fNSteps++);
- if (fgDebug > 0) cout << " WeightPropagation " << fNSteps << " " << fPositionNew << endl;
+ if (fgDebug > 0) printf(" WeightPropagation %d %.3f %.3f %.3f \n", fNSteps,
+ (*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPositionNew);
return;
}
AliMUONTrackK *trackK;
AliMUONDetElement *detElem = NULL;
+ //sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
+ //sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
+ *fCovariance = *fWeight;
+ // check whether the Invert method returns flag if matrix cannot be inverted,
+ // and do not calculate the Determinant in that case !!!!
+ if (fCovariance->Determinant() != 0) {
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
+ } else {
+ AliWarning(" Determinant fCovariance=0:");
+ }
+ //windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+sigmaB*sigmaB);
+ //windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+sigmaNonB*sigmaNonB);
Bool_t ok = kFALSE;
if (fgTrackReconstructor->GetTrackMethod() == 3 && iz >= 0) {
Int_t *pDEatZ = fgCombi->DEatZ(iz);
Int_t nDetElem = pDEatZ[-1];
//cout << fgCombi->Z(iz) << " " << nDetElem << endl;
+ windowB = fgkNSigma * TMath::Sqrt((*fCovariance)(0,0));
+ windowNonB = fgkNSigma * TMath::Sqrt((*fCovariance)(1,1));
+ if (fgkNSigma > 6) windowB = TMath::Min (windowB, 5.);
+ windowB = TMath::Max (windowB, 2.);
+ windowNonB = TMath::Max (windowNonB, 2.);
for (Int_t i = 0; i < nDetElem; i++) {
detElem = fgCombi->DetElem(pDEatZ[i]);
- if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition)) {
- detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0));
+ if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition, windowNonB, windowB)) {
+ detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), windowNonB, windowB);
hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
ok = kTRUE;
break;
iDhit = 1;
}
-
+ TArrayD branchChi2(20);
for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) {
if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem)
hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
// TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB &&
// hit->GetTrackRefSignal() == 1) { // just for test
// Vector of measurements and covariance matrix
- //fprintf(lun1,"%3d %3d %10.4f %10.4f \n", gAlice->GetEvNumber(), ichamb, x, y);
+ //fprintf(lun1,"%3d %3d %10.4f %10.4f \n", AliRunLoader::GetRunLoader()->GetEventNumber(), ichamb, x, y);
if (TMath::Abs(hit->GetZ()-zEnd) > 0.05) {
// Adjust position: for multiple hits in the chamber or misalignment (Z as a function of X or Y)
//AliWarning(Form(" *** adjust %f %f ", zEnd, hit->GetZ()));
trackParTmp = trackPar;
pointWeightTmp = pointWeight;
hitAdd = hit;
- if (fgDebug > 0) cout << " Added point: " << x << " " << y << " " << dChi2 << endl;
+ if (fgDebug > 0) printf(" Added point (ch, x, y, Chi2): %d %.3f %.3f %.3f\n", ichamb, x, y, dChi2);
+ branchChi2[0] = dChi2;
} else {
// branching: create a new track
trackPtr = fgTrackReconstructor->GetRecTracksPtr();
trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
*trackK = *this;
fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
- if (fgDebug > 0) cout << " ******** New track: " << ichamb << " " << hit->GetTTRTrack() << " " << 1/(trackPar)(4,0) << " " << hit->GetBendingCoor() << " " << hit->GetNonBendingCoor() << " " << fNmbTrackHits << " " << nRecTracks << endl;
+ if (fgDebug > 0) printf(" ******** New track (ch, hit, x, y, mom, Chi2, nhits, cand): %d %d %.3f %.3f %.3f %.3f %d %d\n", ichamb, hit->GetTTRTrack(), hit->GetNonBendingCoor(), hit->GetBendingCoor(), 1/(trackPar)(4,0), dChi2, fNmbTrackHits, nRecTracks);
trackK->fRecover = 0;
*(trackK->fTrackPar) = trackPar;
*(trackK->fWeight) += pointWeight;
cout << hitLoop->GetBendingCoor() << " ";
cout << hitLoop->GetNonBendingCoor() << " ";
cout << hitLoop->GetZ() << " " << " ";
- cout << hitLoop->GetTrackRefSignal() << " " << " ";
cout << hitLoop->GetTTRTrack() << endl;
- printf(" ** %d %10.4f %10.4f %10.4f %d %d \n",
+ printf(" ** %d %10.4f %10.4f %10.4f\n",
hitLoop->GetChamberNumber(), hitLoop->GetBendingCoor(),
- hitLoop->GetNonBendingCoor(), hitLoop->GetZ(),
- hitLoop->GetTrackRefSignal(), hitLoop->GetTTRTrack());
+ hitLoop->GetNonBendingCoor(), hitLoop->GetZ());
}
}
//add point
trackK->fTrackDir = 1;
trackK->fBPFlag = kTRUE;
}
+ if (nHitsOK > branchChi2.GetSize()) branchChi2.Set(branchChi2.GetSize()+10);
+ branchChi2[nHitsOK-1] = dChi2;
}
}
}
}
*/
} // if (fTrackDir > 0)
+ // Check for maximum number of branches - exclude excessive
+ if (nHitsOK > 1) CheckBranches(branchChi2, nHitsOK);
}
return ok;
}
+ //__________________________________________________________________________
+void AliMUONTrackK::CheckBranches(TArrayD &branchChi2, Int_t nBranch)
+{
+/// Check for maximum number of branches - exclude excessive
+
+ Int_t nBranchMax = 5;
+ if (nBranch <= nBranchMax) return;
+
+ Double_t *chi2 = branchChi2.GetArray();
+ Int_t *indx = new Int_t [nBranch];
+ TMath::Sort (nBranch, chi2, indx, kFALSE);
+ TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
+ Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
+ Int_t ibeg = nRecTracks - nBranch;
+
+ // Discard excessive branches with higher Chi2 contribution
+ for (Int_t i = nBranchMax; i < nBranch; ++i) {
+ if (indx[i] == 0) {
+ // Discard current track
+ SetRecover(-1);
+ continue;
+ }
+ Int_t j = ibeg + indx[i];
+ AliMUONTrackK *trackK = (AliMUONTrackK*) trackPtr->UncheckedAt(j);
+ trackK->SetRecover(-1);
+ }
+ delete [] indx;
+}
+
//__________________________________________________________________________
void AliMUONTrackK::TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatrixD &trackParTmp, Double_t &dChi2)
{
momentum = 1/(*fTrackParNew)(4,0); // particle momentum
//velo = momentum/TMath::Sqrt(momentum*momentum+muonMass*muonMass); // velocity/c for muon hypothesis
velo = 1; // relativistic
- path = TMath::Abs(fgTrackReconstructor->GetChamberThicknessInX0()/cosAlph/cosBeta); // path length
+ path = TMath::Abs(AliMUONConstants::ChamberThicknessInX0()/cosAlph/cosBeta); // path length
theta0 = 0.0136/velo/momentum*TMath::Sqrt(path)*(1+0.038*TMath::Log(path)); // projected scattering angle
(*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
//__________________________________________________________________________
void AliMUONTrackK::SortHits(Int_t iflag, TObjArray *array)
{
-/// Sort hits in Z if the seed segment in the last but one station
+/// Sort hits in Z if the seed segment is in the last but one station
/// (if iflag==0 in descending order in abs(z), if !=0 - unsort)
if (iflag && ((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetChamberNumber() == 6) return;
{
/// Compute track parameters at hit positions (as for AliMUONTrack)
- // Set number of hits per track
- SetNTrackHits(fNmbTrackHits);
// Set Chi2
+ SetTrackQuality(1); // compute Chi2
SetFitFMin(fChi2);
- // Set track parameters at vertex
+ // Set track parameters at hits
AliMUONTrackParam trackParam;
SetTrackParam(&trackParam, fTrackPar, fPosition);
- SetTrackParamAtVertex(&trackParam);
-
- // Set track parameters at hits
for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
if ((*fChi2Smooth)[i] < 0) {
// Propagate through last chambers
- trackParam.ExtrapToZ(((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
+ AliMUONTrackExtrap::ExtrapToZ(&trackParam, ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
} else {
// Take saved info
SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
}
- AddTrackParamAtHit(&trackParam);
+ AddTrackParamAtHit(&trackParam,(AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
// Fill array of HitForRec's
AddHitForRecAtHit((AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
}
trackParam->SetZ(z);
}
- //__________________________________________________________________________
-void AliMUONTrackK::Branson(void)
-{
-/// Propagates track to the vertex thru absorber using Branson correction
-/// (makes use of the AliMUONTrackParam class)
-
- //AliMUONTrackParam *trackParam = new AliMUONTrackParam();
- AliMUONTrackParam trackParam = AliMUONTrackParam();
- /*
- trackParam->SetBendingCoor((*fTrackPar)(0,0));
- trackParam->SetNonBendingCoor((*fTrackPar)(1,0));
- trackParam->SetBendingSlope(TMath::Tan((*fTrackPar)(2,0)));
- trackParam->SetNonBendingSlope(TMath::Tan((*fTrackPar)(3,0))/TMath::Cos((*fTrackPar)(2,0)));
- trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0)));
- trackParam->SetZ(fPosition);
- */
- SetTrackParam(&trackParam, fTrackPar, fPosition);
-
- trackParam.ExtrapToVertex(Double_t(0.), Double_t(0.), Double_t(0.));
- //trackParam.ExtrapToVertex();
-
- (*fTrackPar)(0,0) = trackParam.GetBendingCoor();
- (*fTrackPar)(1,0) = trackParam.GetNonBendingCoor();
- (*fTrackPar)(2,0) = TMath::ATan(trackParam.GetBendingSlope());
- (*fTrackPar)(3,0) = TMath::ATan(TMath::Cos((*fTrackPar)(2,0))*trackParam.GetNonBendingSlope());
- (*fTrackPar)(4,0) = TMath::Cos((*fTrackPar)(3,0))*trackParam.GetInverseBendingMomentum();
- fPosition = trackParam.GetZ();
- //delete trackParam;
- if (fgDebug > 0) cout << 1/(*fTrackPar)(4,0) << " " << fPosition << " " << (*fTrackPar)(0,0) << endl;
-
- // Get covariance matrix
- *fCovariance = *fWeight;
- if (fCovariance->Determinant() != 0) {
- Int_t ifail;
- mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
- } else {
- AliWarning(" Determinant fCovariance=0:");
- }
-}
-
- //__________________________________________________________________________
-void AliMUONTrackK::GoToZ(Double_t zEnd)
-{
-/// Propagates track to given Z
-
- ParPropagation(zEnd);
- MSThin(1); // multiple scattering in the chamber
- WeightPropagation(zEnd, kFALSE);
- fPosition = fPositionNew;
- *fTrackPar = *fTrackParNew;
-}
-
- //__________________________________________________________________________
-void AliMUONTrackK::GoToVertex(Int_t iflag)
-{
-/// Version 3.08
-/// Propagates track to the vertex
-/// All material constants are taken from AliRoot
-
- static Double_t x01[5] = { 24.282, // C
- 24.282, // C
- 11.274, // Concrete
- 1.758, // Fe
- 1.758}; // Fe (cm)
- // inner part theta < 3 degrees
- static Double_t x02[5] = { 30413, // Air
- 24.282, // C
- 11.274, // Concrete
- 1.758, // Fe
- 0.369}; // W (cm)
- // z positions of the materials inside the absober outer part theta > 3 degres
- static Double_t zPos[10] = {-90, -105, -315, -443, -468};
-
- Double_t dZ, r0Norm, x0, deltaP, dChi2, pTotal, pOld;
- AliMUONHitForRec *hit;
- AliMUONRawCluster *clus;
- TClonesArray *rawclusters;
-
- // First step to the rear end of the absorber
- Double_t zRear = -503;
- GoToZ(zRear);
- Double_t tan3 = TMath::Tan(3./180*TMath::Pi());
-
- // Go through absorber
- pOld = 1/(*fTrackPar)(4,0);
- Double_t r0Rear = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
- (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
- r0Rear = TMath::Sqrt(r0Rear)/TMath::Abs(fPosition)/tan3;
- r0Norm = r0Rear;
- Double_t p0, cos25, cos60;
- if (!iflag) goto vertex;
-
- for (Int_t i=4; i>=0; i--) {
- ParPropagation(zPos[i]);
- WeightPropagation(zPos[i], kFALSE);
- dZ = TMath::Abs (fPositionNew-fPosition);
- if (r0Norm > 1) x0 = x01[i];
- else x0 = x02[i];
- MSLine(dZ,x0); // multiple scattering in the medium (linear approximation)
- fPosition = fPositionNew;
- *fTrackPar = *fTrackParNew;
- r0Norm = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
- (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
- r0Norm = TMath::Sqrt(r0Norm)/TMath::Abs(fPosition)/tan3;
- }
- // Correct momentum for energy losses
- pTotal = 1/TMath::Abs((*fTrackPar)(4,0));
- p0 = pTotal;
- cos25 = TMath::Cos(2.5/180*TMath::Pi());
- cos60 = TMath::Cos(6.0/180*TMath::Pi());
- for (Int_t j=0; j<1; j++) {
- /*
- if (r0Rear > 1) {
- if (p0 < 20) {
- deltaP = 2.164 + 0.145e-1*p0 - 0.417e-3*p0*p0;
- } else {
- deltaP = 2.275 + 0.102e-2*p0 - 0.674e-6*p0*p0;
- }
- } else {
- if (p0 < 20) {
- deltaP = 2.581 + 0.188e-1*p0 - 0.398e-3*p0*p0;
- } else {
- deltaP = 2.727 + 0.356e-2*p0 + 0.242e-5*p0*p0;
- }
- }
- */
- /*
- if (r0Rear < 1) {
- //W
- if (p0<15) {
- deltaP = 2.737 + 0.0494*p0 - 0.001123*p0*p0;
- } else {
- deltaP = 3.0643 + 0.01346*p0;
- }
- deltaP *= 0.95;
- } else {
- //Pb
- if (p0<15) {
- deltaP = 2.1380 + 0.0351*p0 - 0.000853*p0*p0;
- } else {
- deltaP = 2.407 + 0.00702*p0;
- }
- deltaP *= 0.95;
- }
- */
- /*
- if (r0Rear < 1) {
- //W
- if (p0<18) {
- deltaP = 2.439 + 0.806e-1*p0 - 0.500e-2*p0*p0 + 0.106e-3*p0*p0*p0;
- } else {
- deltaP = 2.767 + 0.742e-2*p0 - 0.196e-4*p0*p0 + 0.403e-7*p0*p0*p0;
- }
- //deltaP += 0.2;
- deltaP *= cos25;
- } else {
- //Pb
- if (p0<18) {
- deltaP = 2.209 + 0.800e-2*p0;
- } else {
- deltaP = 2.285 + 0.141e-2*p0 - 0.446e-6*p0*p0;
- }
- //deltaP += 0.2;
- deltaP *= cos60;
- }
- deltaP *= 1.1;
- */
- //*
- if (r0Rear < 1) {
- if (p0 < 20) {
- deltaP = 2.5938 + 0.0570 * p0 - 0.001151 * p0 * p0;
- } else {
- deltaP = 3.0714 + 0.011767 * p0;
- }
- deltaP *= 0.75;
- } else {
- if (p0 < 20) {
- deltaP = 2.1207 + 0.05478 * p0 - 0.00145079 * p0 * p0;
- } else {
- deltaP = 2.6069 + 0.0051705 * p0;
- }
- deltaP *= 0.9;
- }
- //*/
-
- p0 = pTotal + deltaP/TMath::Abs(TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)));
- }
- (*fTrackPar)(4,0) = 1/p0*TMath::Sign((Double_t)1.,(*fTrackPar)(4,0));
-
- // Go to the vertex
-vertex:
- ParPropagation((Double_t)0.);
- WeightPropagation((Double_t)0., kFALSE);
- fPosition = fPositionNew;
- //*fTrackPar = *fTrackParNew;
- // Add vertex as a hit
- TMatrixD pointWeight(fgkSize,fgkSize);
- TMatrixD point(fgkSize,1);
- TMatrixD trackParTmp = point;
- point(0,0) = 0; // vertex coordinate - should be taken somewhere
- point(1,0) = 0; // vertex coordinate - should be taken somewhere
- pointWeight(0,0) = 1/1.e-3/1.e-3; // 10 um error
- pointWeight(1,1) = 1/1.e-3/1.e-3; // 10 um error
- TryPoint(point,pointWeight,trackParTmp,dChi2);
- *fTrackPar = trackParTmp;
- *fWeight += pointWeight;
- fChi2 += dChi2; // Chi2
- if (fgDebug < 0) return; // no output
-
- cout << pOld << " " << 1/(*fTrackPar)(4,0) << " " << dChi2 << " " << fChi2 << " " << fNmbTrackHits << endl;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- printf ("%5d", hit->GetChamberNumber());
- }
- cout << endl;
- if (fgDebug > 0) {
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
- //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
- printf ("%5d", fgHitForRec->IndexOf(hit));
- }
- cout << endl;
- }
- if (fgTrackReconstructor->GetRecTrackRefHits()) {
- // from track ref. hits
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- cout << hit->GetTTRTrack() + hit->GetTrackRefSignal()*10000 << " ";
- }
- } else {
- // from raw clusters
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- if (hit->GetHitNumber() < 0) { // combined cluster / track finder
- Int_t index = -hit->GetHitNumber() / 100000;
- Int_t iPos = -hit->GetHitNumber() - index * 100000;
- clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
- } else {
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- }
- printf ("%5d", clus->GetTrack(1)%10000000);
- }
- cout << endl;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- if (hit->GetHitNumber() < 0) { // combined cluster / track finder
- Int_t index = -hit->GetHitNumber() / 100000;
- Int_t iPos = -hit->GetHitNumber() - index * 100000;
- clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
- } else {
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- }
- if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000);
- else printf ("%5s", " ");
- }
- }
- cout << endl;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
- cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
- //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " ";
- }
- cout << endl;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) printf("%8.4f", (*fChi2Smooth)[i1]);
- cout << endl;
- cout << "---------------------------------------------------" << endl;
-
- // Get covariance matrix
- /* Not needed - covariance matrix is not interesting to anybody
- *fCovariance = *fWeight;
- if (fCovariance->Determinant() != 0) {
- Int_t ifail;
- mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
- } else {
- AliWarning(" Determinant fCovariance=0:" );
- }
- */
-}
-
//__________________________________________________________________________
void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
{
if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
// Check if the outlier is not from the seed segment
AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
- if (skipHit == fStartSegment->GetHitForRec1() || skipHit == fStartSegment->GetHitForRec2()) {
+ if (skipHit == (AliMUONHitForRec*) fStartSegment->First() || skipHit == (AliMUONHitForRec*) fStartSegment->Second()) {
//DropBranches(fStartSegment); // drop all tracks with the same seed segment
return kFALSE; // to be changed probably
}
// Remove all saved steps and smoother matrices after the skipped hit
RemoveMatrices(skipHit->GetZ());
- //AZ(z->-z) if (skipHit->GetZ() > fStartSegment->GetHitForRec2()->GetZ() || !fNSteps) {
- if (TMath::Abs(skipHit->GetZ()) > TMath::Abs(fStartSegment->GetHitForRec2()->GetZ()) || !fNSteps) {
+ //AZ(z->-z) if (skipHit->GetZ() > ((AliMUONHitForRec*) fStartSegment->Second())->GetZ() || !fNSteps) {
+ if (TMath::Abs(skipHit->GetZ()) > TMath::Abs( ((AliMUONHitForRec*) fStartSegment->Second())->GetZ()) || !fNSteps) {
// Propagation toward high Z or skipped hit next to segment -
// start track from segment
trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
}
// Check if the outlier is not from the seed segment
AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
- if (hit == fStartSegment->GetHitForRec1() || hit == fStartSegment->GetHitForRec2()) return; // to be changed probably
+ if (hit == (AliMUONHitForRec*) fStartSegment->First() || hit == (AliMUONHitForRec*) fStartSegment->Second()) return; // to be changed probably
// Check for missing station
Int_t ok = 1;
Int_t flag = 1;
AliMUONHitForRec *hit = 0;
- if (fgTrackReconstructor->GetRecTrackRefHits()) {
- // from track ref. hits
- for (Int_t j=0; j<fNmbTrackHits; j++) {
- hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(j);
- if (hit->GetTTRTrack() > 1) { flag = 0; break; }
- }
- for (Int_t j=0; j<fNmbTrackHits; j++) {
- printf("%10.4f", GetChi2PerPoint(j));
- if (GetChi2PerPoint(j) > -0.1) {
- hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(j);
- fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(j));
- fprintf(lun, "%3d %3d %3d \n", hit->GetTrackRefSignal(), hit->GetTTRTrack(), flag);
- }
- }
- printf("\n");
- } else {
// from raw clusters
- AliMUONRawCluster *clus = 0;
- TClonesArray *rawclusters = 0;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- if (TMath::Abs(clus->GetTrack(1)-1) < 2) {
- if (clus->GetTrack(2)) flag = 2;
- continue;
- }
- if (clus->GetTrack(2) && TMath::Abs(clus->GetTrack(2)-1) < 2) {
- flag = 3;
- continue;
- }
- flag = 0;
- break;
+ AliMUONRawCluster *clus = 0;
+ TClonesArray *rawclusters = 0;
+ for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
+ hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
+ rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
+ clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
+ if (TMath::Abs(clus->GetTrack(1)-1) < 2) {
+ if (clus->GetTrack(2)) flag = 2;
+ continue;
}
+ if (clus->GetTrack(2) && TMath::Abs(clus->GetTrack(2)-1) < 2) {
+ flag = 3;
+ continue;
+ }
+ flag = 0;
+ break;
+
Int_t sig[2]={1,1}, tid[2]={0};
for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
if (GetChi2PerPoint(i1) < -0.1) continue;
tid[j] = clus->GetTrack(j+1) - 1;
if (clus->GetTrack(j+1) < 0) { sig[j] = 0; tid[j] = 999; }
}
- fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
+ //fprintf(lun,"%3d %3d %10.4f", AliRunLoader::GetRunLoader()->GetEventNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
if (!(clus->GetTrack(2))) fprintf(lun, "%3d %3d", sig[0], tid[0]); // simple cluster
else { // track overlap
fprintf(lun, "%3d %3d", TMath::Max(sig[0],sig[1]), TMath::Min(tid[0],tid[1]));
}
//__________________________________________________________________________
-void AliMUONTrackK::DropBranches(AliMUONSegment *segment)
+void AliMUONTrackK::DropBranches(AliMUONObjectPair *segment)
{
/// Drop all candidates with the same seed segment
Int_t nRecTracks;
TClonesArray *trackPtr;
AliMUONTrackK *trackK;
+ AliMUONObjectPair *trackKSegment;
trackPtr = fgTrackReconstructor->GetRecTracksPtr();
nRecTracks = fgTrackReconstructor->GetNRecTracks();
for (Int_t i=icand+1; i<nRecTracks; i++) {
trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
+ trackKSegment = trackK->fStartSegment;
if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
if (trackK->GetRecover() < 0) continue;
- if (trackK->fStartSegment == segment) trackK->SetRecover(-1);
+ if (trackKSegment->First() == segment->First() &&
+ trackKSegment->Second() == segment->Second()) trackK->SetRecover(-1);
}
if (fgDebug >= 0) cout << " Drop segment " << endl;
}
Int_t AliMUONTrackK::GetStation0(void)
{
/// Return seed station number
- return fStartSegment->GetHitForRec1()->GetChamberNumber() / 2;
+ return ((AliMUONHitForRec*) fStartSegment->First())->GetChamberNumber() / 2;
}
//__________________________________________________________________________
delete hitArray;
return same;
}
+
+//______________________________________________________________________________
+ void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
+{
+///*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
+///*-* ==========================
+///*-* inverts a symmetric matrix. matrix is first scaled to
+///*-* have all ones on the diagonal (equivalent to change of units)
+///*-* but no pivoting is done since matrix is positive-definite.
+///*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
+
+ // taken from TMinuit package of Root (l>=n)
+ // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
+ // Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
+ Double_t * localVERTs = new Double_t[n];
+ Double_t * localVERTq = new Double_t[n];
+ Double_t * localVERTpp = new Double_t[n];
+ // fMaxint changed to localMaxint
+ Int_t localMaxint = n;
+
+ /* System generated locals */
+ Int_t aOffset;
+
+ /* Local variables */
+ Double_t si;
+ Int_t i, j, k, kp1, km1;
+
+ /* Parameter adjustments */
+ aOffset = l + 1;
+ a -= aOffset;
+
+ /* Function Body */
+ ifail = 0;
+ if (n < 1) goto L100;
+ if (n > localMaxint) goto L100;
+//*-*- scale matrix by sqrt of diag elements
+ for (i = 1; i <= n; ++i) {
+ si = a[i + i*l];
+ if (si <= 0) goto L100;
+ localVERTs[i-1] = 1 / TMath::Sqrt(si);
+ }
+ for (i = 1; i <= n; ++i) {
+ for (j = 1; j <= n; ++j) {
+ a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
+ }
+ }
+//*-*- . . . start main loop . . . .
+ for (i = 1; i <= n; ++i) {
+ k = i;
+//*-*- preparation for elimination step1
+ if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
+ else goto L100;
+ localVERTpp[k-1] = 1;
+ a[k + k*l] = 0;
+ kp1 = k + 1;
+ km1 = k - 1;
+ if (km1 < 0) goto L100;
+ else if (km1 == 0) goto L50;
+ else goto L40;
+L40:
+ for (j = 1; j <= km1; ++j) {
+ localVERTpp[j-1] = a[j + k*l];
+ localVERTq[j-1] = a[j + k*l]*localVERTq[k-1];
+ a[j + k*l] = 0;
+ }
+L50:
+ if (k - n < 0) goto L51;
+ else if (k - n == 0) goto L60;
+ else goto L100;
+L51:
+ for (j = kp1; j <= n; ++j) {
+ localVERTpp[j-1] = a[k + j*l];
+ localVERTq[j-1] = -a[k + j*l]*localVERTq[k-1];
+ a[k + j*l] = 0;
+ }
+//*-*- elimination proper
+L60:
+ for (j = 1; j <= n; ++j) {
+ for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
+ }
+ }
+//*-*- elements of left diagonal and unscaling
+ for (j = 1; j <= n; ++j) {
+ for (k = 1; k <= j; ++k) {
+ a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
+ a[j + k*l] = a[k + j*l];
+ }
+ }
+ delete [] localVERTs;
+ delete [] localVERTq;
+ delete [] localVERTpp;
+ return;
+//*-*- failure return
+L100:
+ delete [] localVERTs;
+ delete [] localVERTq;
+ delete [] localVERTpp;
+ ifail = 1;
+} /* mnvertLocal */
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff