#include <TMatrixD.h>
#include "AliMUONTrackK.h"
-#include "AliCallf77.h"
#include "AliMUON.h"
-#include "AliMUONChamber.h"
+
#include "AliMUONTrackReconstructor.h"
#include "AliMagF.h"
#include "AliMUONSegment.h"
#include "AliMUONHitForRec.h"
#include "AliMUONRawCluster.h"
#include "AliMUONTrackParam.h"
+#include "AliMUONEventRecoCombi.h"
+#include "AliMUONDetElement.h"
#include "AliRun.h"
#include "AliLog.h"
const Int_t AliMUONTrackK::fgkSize = 5;
-const Int_t AliMUONTrackK::fgkNSigma = 6;
-const Double_t AliMUONTrackK::fgkChi2max = 25;
+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 mnvertLocalK(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
+void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); // from AliMUONTrack
ClassImp(AliMUONTrackK) // Class implementation in ROOT context
-Int_t AliMUONTrackK::fgDebug = -1;
+ Int_t AliMUONTrackK::fgDebug = -1; //-1;
Int_t AliMUONTrackK::fgNOfPoints = 0;
-AliMUON* AliMUONTrackK::fgMUON = NULL;
AliMUONTrackReconstructor* AliMUONTrackK::fgTrackReconstructor = NULL;
TClonesArray* AliMUONTrackK::fgHitForRec = NULL;
+AliMUONEventRecoCombi *AliMUONTrackK::fgCombi = NULL;
//FILE *lun1 = fopen("window.dat","w");
//__________________________________________________________________________
{
// Default constructor
- fgTrackReconstructor = NULL; // pointer to track reconstructor
- fgMUON = NULL; // pointer to Muon module
+ fgTrackReconstructor = NULL; // pointer to event reconstructor
fgHitForRec = NULL; // pointer to points
fgNOfPoints = 0; // number of points
//__________________________________________________________________________
AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClonesArray *hitForRec)
- //AZ: TObject()
: AliMUONTrack()
{
// Constructor
if (!TrackReconstructor) return;
- fgTrackReconstructor = TrackReconstructor; // pointer to track reconstructor
- fgMUON = (AliMUON*) gAlice->GetModule("MUON"); // pointer to Muon module
+ fgTrackReconstructor = TrackReconstructor; // pointer to event reconstructor
fgHitForRec = hitForRec; // pointer to points
fgNOfPoints = fgHitForRec->GetEntriesFast(); // number of points
+ fgCombi = NULL;
+ if (fgTrackReconstructor->GetTrackMethod() == 3) fgCombi = AliMUONEventRecoCombi::Instance();
fStartSegment = NULL;
fTrackHitsPtr = NULL;
//__________________________________________________________________________
AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
- //AZ: TObject()
: AliMUONTrack(segment, segment, fgTrackReconstructor)
{
// Constructor from a segment
dY = hit2->GetBendingCoor() - hit1->GetBendingCoor();
dX = hit2->GetNonBendingCoor() - hit1->GetNonBendingCoor();
(*fTrackPar)(2,0) = TMath::ATan2(dY,dZ); // alpha
- (*fTrackPar)(3,0) = TMath::ATan2(dX,dZ/TMath::Cos((*fTrackPar)(2,0))); // beta
+ 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) *= TMath::Cos((*fTrackPar)(3,0)); // 1/p
// Evaluate covariance (and weight) matrix
fParSmooth = NULL;
if (fgDebug < 0 ) return;
- cout << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
+ cout << fgTrackReconstructor->GetNRecTracks()-1 << " " << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
if (fgTrackReconstructor->GetRecTrackRefHits()) {
// from track ref. hits
cout << ((AliMUONHitForRec*)((*fTrackHitsPtr)[0]))->GetTTRTrack() << "<-->" << ((AliMUONHitForRec*)((*fTrackHitsPtr)[1]))->GetTTRTrack() << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
hit1 = (AliMUONHitForRec*) ((*fTrackHitsPtr)[i]);
rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit1->GetChamberNumber());
clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit1->GetHitNumber());
- cout << clus->GetTrack(1)-1;
- if (clus->GetTrack(2) != 0) cout << " " << clus->GetTrack(2)-1;
+ cout << clus->GetTrack(1);
+ if (clus->GetTrack(2) != -1) cout << " " << clus->GetTrack(2);
if (i == 0) cout << " <--> ";
}
cout << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
//__________________________________________________________________________
AliMUONTrackK::AliMUONTrackK (const AliMUONTrackK& source)
- //AZ: TObject(source)
: AliMUONTrack(source)
{
// Protected copy constructor
// Evaluate covariance (and weight) matrix for track candidate
Double_t sigmaB, sigmaNonB, tanA, tanB, dAdY, rad, dBdX, dBdY;
+ dZ = -dZ;
sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
(*fWeight)(1,3) = dBdX*(*fWeight)(1,1); // <xb>
(*fWeight)(3,1) = (*fWeight)(1,3);
- //(*fWeight)(4,4) = ((*fTrackPar)(4,0)*0.2)*((*fTrackPar)(4,0)*0.2); // error 20%
(*fWeight)(4,4) = ((*fTrackPar)(4,0)*0.5)*((*fTrackPar)(4,0)*0.5); // error 50%
// check whether the Invert method returns flag if matrix cannot be inverted,
// and do not calculate the Determinant in that case !!!!
if (fWeight->Determinant() != 0) {
// fWeight->Invert();
- Int_t ifailWeight;
- mnvertLocalK(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifailWeight);
+ Int_t ifail;
+ mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
AliWarning(" Determinant fWeight=0:");
}
// Follows track through detector stations
Bool_t miss, success;
Int_t ichamb, iFB, iMin, iMax, dChamb, ichambOK;
- Int_t ihit, firstIndx, lastIndx, currIndx;
+ Int_t ihit, firstIndx = -1, lastIndx = -1, currIndx = -1, iz0 = -1, iz = -1;
Double_t zEnd, dChi2;
- AliMUONHitForRec *hitAdd, *firstHit, *lastHit, *hit;
+ AliMUONHitForRec *hitAdd, *firstHit = NULL, *lastHit = NULL, *hit;
AliMUONRawCluster *clus;
TClonesArray *rawclusters;
hit = 0; clus = 0; rawclusters = 0;
miss = success = kTRUE;
Int_t endOfProp = 0;
- iFB = (ichamEnd == ichamBeg) ? fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
+ //iFB = (ichamEnd == ichamBeg) ? fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
+ iFB = (ichamEnd == ichamBeg) ? -fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
iMin = TMath::Min(ichamEnd,ichamBeg);
iMax = TMath::Max(ichamEnd,ichamBeg);
ichamb = ichamBeg;
fRecover = 0;
ichambOK = ((AliMUONHitForRec*)((*fTrackHitsPtr)[fNTrackHits-1]))->GetChamberNumber();
//if (ichambOK >= 8 && ((AliMUONHitForRec*)((*fTrackHitsPtr)[0]))->GetChamberNumber() == 6) ichambOK = 6;
- currIndx = fgHitForRec->IndexOf(fSkipHit);
+ if (fgTrackReconstructor->GetTrackMethod() == 3 &&
+ fSkipHit->GetHitNumber() < 0) {
+ iz0 = fgCombi->IZfromHit(fSkipHit);
+ currIndx = -1;
+ }
+ else currIndx = fgHitForRec->IndexOf(fSkipHit);
} else {
// outlier
fTrackHitsPtr->Compress();
// Get indices of the 1'st and last hits on the track candidate
firstHit = (AliMUONHitForRec*) fTrackHitsPtr->First();
lastHit = (AliMUONHitForRec*) fTrackHitsPtr->Last();
- firstIndx = fgHitForRec->IndexOf(firstHit);
- lastIndx = fgHitForRec->IndexOf(lastHit);
- currIndx = TMath::Abs (TMath::Max(firstIndx*iFB,lastIndx*iFB));
+ if (fgTrackReconstructor->GetTrackMethod() == 3 &&
+ lastHit->GetHitNumber() < 0) iz0 = fgCombi->IZfromHit(lastHit);
+ else {
+ firstIndx = fgHitForRec->IndexOf(firstHit);
+ lastIndx = fgHitForRec->IndexOf(lastHit);
+ currIndx = TMath::Abs (TMath::Max(firstIndx*iFB,lastIndx*iFB));
+ }
}
- while (ichamb>=iMin && ichamb<=iMax) {
+ if (iz0 < 0) iz0 = iFB;
+ while (ichamb >= iMin && ichamb <= iMax) {
// Find the closest hit in Z, not belonging to the current plane
if (Back) {
// backpropagation
break;
}
}
+
+ // Combined cluster / track finder
+ if (zEnd > 999 && iFB < 0 && fgTrackReconstructor->GetTrackMethod() == 3) {
+ currIndx = -2;
+ AliMUONHitForRec hitTmp;
+ for (iz = iz0 - iFB; iz < fgCombi->Nz(); iz++) {
+ if (TMath::Abs(fgCombi->Z(iz)-fPosition) < 0.5) continue;
+ Int_t *pDEatZ = fgCombi->DEatZ(iz);
+ //cout << iz << " " << fgCombi->Z(iz) << endl;
+ zEnd = fgCombi->Z(iz);
+ iz0 = iz;
+ AliMUONDetElement *detElem = fgCombi->DetElem(pDEatZ[0]);
+ hitAdd = &hitTmp;
+ hitAdd->SetChamberNumber(detElem->Chamber());
+ //hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
+ if (hitAdd) break;
+ }
+ }
}
- //AZ(Z->-Z) if (zEnd<-999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
if (zEnd>999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
else {
// Check if there is a chamber without hits
if (zEnd>999 || TMath::Abs(hitAdd->GetChamberNumber()-ichamb) > 1) {
if (!Back && zEnd<999) currIndx -= iFB;
ichamb += iFB;
- zEnd = (&(fgMUON->Chamber(ichamb)))->Z();
+ zEnd = AliMUONConstants::DefaultChamberZ(ichamb);
miss = kTRUE;
} else {
ichamb = hitAdd->GetChamberNumber();
Back = kFALSE;
fRecover =0;
ichambOK = ((AliMUONHitForRec*)((*fTrackHitsPtr)[fNTrackHits-1]))->GetChamberNumber();
+ //? if (fgTrackReconstructor->GetTrackMethod() != 3) currIndx = fgHitForRec->IndexOf(fSkipHit);
currIndx = fgHitForRec->IndexOf(fSkipHit);
}
TryPoint(point,pointWeight,trackParTmp,dChi2);
*fTrackPar = trackParTmp;
*fWeight += pointWeight;
- //AZ(z->-z) if (fTrackDir < 0) AddMatrices (this, dChi2, hitAdd);
if (fTrackDir > 0) AddMatrices (this, dChi2, hitAdd);
fChi2 += dChi2; // Chi2
} else *fTrackPar = *fTrackParNew; // adjust end point to the last hit
currIndx ++;
} else {
// forward propagator
- if (miss || !FindPoint(ichamb,zEnd,currIndx,iFB,hitAdd)) {
+ if (miss || !FindPoint(ichamb, zEnd, currIndx, iFB, hitAdd, iz)) {
// missing point
*fTrackPar = *fTrackParNew;
} else {
Int_t iFB, nTries;
Double_t dZ, step, distance, charge;
Double_t vGeant3[7], vGeant3New[7];
- AliMUONTrackParam* trackParam = 0;
+ AliMUONTrackParam trackParam;
+
nTries = 0;
// First step using linear extrapolation
dZ = zEnd - fPosition;
//AZ(z->-z) charge = iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
charge = -iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
SetGeantParam(vGeant3,iFB);
+ //fTrackParNew->Print();
// Check if overstep
do {
step = TMath::Abs(step);
// Propagate parameters
- trackParam->ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
- // extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
+ trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
+ //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
distance = zEnd - vGeant3New[2];
step *= dZ/(vGeant3New[2]-fPositionNew);
nTries ++;
} while (distance*iFB < 0 && TMath::Abs(distance) > fgkEpsilon);
GetFromGeantParam(vGeant3New,iFB);
+ //fTrackParNew->Print();
// Position adjustment (until within tolerance)
while (TMath::Abs(distance) > fgkEpsilon) {
do {
// binary search
// Propagate parameters
- // extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
- trackParam->ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
-
+ trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
+ //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
distance = zEnd - vGeant3New[2];
step /= 2;
nTries ++;
GetFromGeantParam(vGeant3New,iFB);
}
//cout << nTries << endl;
+ //fTrackParNew->Print();
return;
}
-/*
- //__________________________________________________________________________
-void AliMUONTrackK::WeightPropagation(void)
-{
- // Propagation of the weight matrix
- // W = DtWD, where D is Jacobian
- // !!! not implemented TMatrixD weight1(*fJacob,TMatrixD::kAtBA,*fWeight); // DtWD
- TMatrixD weight1(*fWeight,TMatrixD::kMult,*fJacob); // WD
- *fWeight = TMatrixD(*fJacob,TMatrixD::kTransposeMult,weight1); // DtWD
- return;
-}
-*/
//__________________________________________________________________________
void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
{
// Save initial and propagated parameters
TMatrixD trackPar0 = *fTrackPar;
TMatrixD trackParNew0 = *fTrackParNew;
- Double_t savePosition = fPositionNew;
// Get covariance matrix
*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) {
- // fCovariance->Invert();
- Int_t ifailCov;
- mnvertLocalK(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
+ //fCovariance->Print();
} else {
AliWarning(" Determinant fCovariance=0:");
}
- // Loop over parameters to find change of the initial vs propagated ones
- zEnd = fPosition;
- fPosition = fPositionNew;
+ // Loop over parameters to find change of the propagated vs initial ones
for (i=0; i<fgkSize; i++) {
dPar = TMath::Sqrt((*fCovariance)(i,i));
- *fTrackPar = trackParNew0;
+ *fTrackPar = trackPar0;
+ if (i == 4) dPar *= TMath::Sign(1.,-trackPar0(4,0)); // 1/p
(*fTrackPar)(i,0) += dPar;
ParPropagation(zEnd);
for (j=0; j<fgkSize; j++) {
- jacob(j,i) = ((*fTrackParNew)(j,0)-trackPar0(j,0))/dPar;
+ jacob(j,i) = ((*fTrackParNew)(j,0)-trackParNew0(j,0))/dPar;
}
}
- //jacob->Print();
//trackParNew0.Print();
//TMatrixD par1(jacob,TMatrixD::kMult,trackPar0); //
//par1.Print();
- /*
+ TMatrixD jacob0 = jacob;
if (jacob.Determinant() != 0) {
- // jacob.Invert();
+ jacob.Invert();
} else {
- cout << " ***** Warning in WeightPropagation: Determinant jacob=0:" << endl;
+ AliWarning(" Determinant jacob=0:");
}
- */
TMatrixD weight1(*fWeight,TMatrixD::kMult,jacob); // WD
*fWeight = TMatrixD(jacob,TMatrixD::kTransposeMult,weight1); // DtWD
- //fWeight->Print();
// Restore initial and propagated parameters
*fTrackPar = trackPar0;
*fTrackParNew = trackParNew0;
- fPosition = zEnd;
- fPositionNew = savePosition;
// Save for smoother
if (!smooth) return; // do not use smoother
- //AZ(z->-z) if (fTrackDir > 0) return; // only for propagation towards int. point
if (fTrackDir < 0) return; // only for propagation towards int. point
TMatrixD *tmp = new TMatrixD(*fTrackParNew); // extrapolated parameters
fParExtrap->Add(tmp);
tmp->SetUniqueID(1);
*fCovariance = *fWeight;
if (fCovariance->Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
AliWarning(" Determinant fCovariance=0:");
}
tmp = new TMatrixD(*fCovariance); // filtered covariance (if no measurement will be found)
fCovFilter->Add(tmp);
tmp->SetUniqueID(1);
- tmp = new TMatrixD(jacob); // Jacobian
- tmp->Invert();
+ tmp = new TMatrixD(jacob0); // Jacobian
fJacob->Add(tmp);
tmp->SetUniqueID(1);
if (fSteps->fN <= fNSteps) fSteps->Set(fSteps->fN+10);
}
//__________________________________________________________________________
-Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd)
+Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz)
{
// Picks up point within a window for the chamber No ichamb
// Split the track if there are more than 1 hit
TClonesArray *trackPtr;
AliMUONHitForRec *hit, *hitLoop;
AliMUONTrackK *trackK;
+ AliMUONDetElement *detElem = NULL;
Bool_t ok = kFALSE;
+
+ if (fgTrackReconstructor->GetTrackMethod() == 3 && iz >= 0) {
+ // Combined cluster / track finder
+ // Check if extrapolated track passes thru det. elems. at iz
+ Int_t *pDEatZ = fgCombi->DEatZ(iz);
+ Int_t nDetElem = pDEatZ[-1];
+ //cout << fgCombi->Z(iz) << " " << nDetElem << endl;
+ 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));
+ hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
+ ok = kTRUE;
+ break;
+ }
+ }
+ if (!ok) return ok; // outside det. elem.
+ ok = kFALSE;
+ }
+
//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) {
- // fCovariance->Invert();
- Int_t ifailCov;
- mnvertLocalK(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fCovariance=0:");
+ AliWarning(" Determinant fCovariance=0:");
}
//windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+sigmaB*sigmaB);
//windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+sigmaNonB*sigmaNonB);
TMatrixD point(fgkSize,1);
TMatrixD trackPar = point;
TMatrixD trackParTmp = point;
- Int_t nHitsOK = 0;
+ Int_t nHitsOK = 0, ihitB = currIndx, ihitE = fgNOfPoints, iDhit = iFB;
Double_t zLast;
zLast = ((AliMUONHitForRec*)fTrackHitsPtr->Last())->GetZ();
+ if (fgTrackReconstructor->GetTrackMethod() == 3 && detElem) {
+ ihitB = 0;
+ ihitE = detElem->NHitsForRec();
+ iDhit = 1;
+ }
- for (ihit=currIndx; ihit>=0 && ihit<fgNOfPoints; ihit+=iFB) {
- hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
+
+ for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) {
+ if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem)
+ hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
+ else hit = (AliMUONHitForRec*) detElem->HitsForRec()->UncheckedAt(ihit);
if (hit->GetChamberNumber() == ichamb) {
//if (TMath::Abs(hit->GetZ()-zEnd) < 0.1) {
if (TMath::Abs(hit->GetZ()-zEnd) < 0.5) {
- if (TMath::Abs(hit->GetZ()-zEnd) > 0.1) {
+ //if (TMath::Abs(hit->GetZ()-zEnd) > 0.1) {
+ if (TMath::Abs(hit->GetZ()-zEnd) > 0.05) {
// adjust position: for multiple hits in the chamber
// (mostly (only?) for GEANT hits)
+ cout << " ******* adjust " << zEnd << " " << hit->GetZ() << endl;
zEnd = hit->GetZ();
*fTrackPar = *fTrackParNew;
ParPropagation(zEnd);
// Get covariance
*fCovariance = *fWeight;
if (fCovariance->Determinant() != 0) {
- //fCovariance->Invert();
- Int_t ifailCov;
- mnvertLocalK(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fCovariance=0:" );
+ AliWarning(" Determinant fCovariance=0:" );
}
}
y = hit->GetBendingCoor();
x = hit->GetNonBendingCoor();
+ if (hit->GetBendingReso2() < 0) {
+ // Combined cluster / track finder
+ hit->SetBendingReso2(fgTrackReconstructor->GetBendingResolution()*
+ fgTrackReconstructor->GetBendingResolution());
+ hit->SetNonBendingReso2(fgTrackReconstructor->GetNonBendingResolution()*
+ fgTrackReconstructor->GetNonBendingResolution());
+ }
windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
- //fprintf(lun1,"%3d %10.4f %10.4f %10.4f \n", ichamb, windowB, windowNonB, TMath::Abs(hit->GetZ()-zLast));
- windowB = TMath::Min (windowB,5.);
+ // windowB = TMath::Min (windowB,5.);
+ if (fgkNSigma > 6) windowB = TMath::Min (windowB,5.);
/*
if (TMath::Abs(hit->GetZ()-zLast) < 50) {
windowB = TMath::Min (windowB,0.5);
// 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);
point.Zero();
point(0,0) = y;
point(1,0) = x;
pointWeight(1,1) = 1/hit->GetNonBendingReso2();
TryPoint(point,pointWeight,trackPar,dChi2);
if (TMath::Abs(1./(trackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) continue; // p < p_min - next hit
+ // if (TMath::Sign(1.,(trackPar)(4,0)*(*fTrackPar)(4,0)) < 0) continue; // change of sign
ok = kTRUE;
nHitsOK++;
//if (nHitsOK > -1) {
trackParTmp = trackPar;
pointWeightTmp = pointWeight;
hitAdd = hit;
+ if (fgDebug > 0) cout << " Added point: " << x << " " << y << " " << dChi2 << endl;
} else {
// branching: create a new track
trackPtr = fgTrackReconstructor->GetRecTracksPtr();
/*
*(trackK->fCovariance) = *(trackK->fWeight);
if (trackK->fCovariance->Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
}
cout << (*(trackK->fCovariance))(0,0) << " " << (*(trackK->fCovariance))(1,1) << " " << (*fCovariance)(0,0) << " " << (*fCovariance)(1,1) << endl;
*/
// Add filtered matrices for the smoother
- //AZ(z->-z) if (fTrackDir < 0) {
if (fTrackDir > 0) {
if (nHitsOK > 2) { // check for position adjustment
for (Int_t i=trackK->fNSteps-1; i>=0; i--) {
RemoveMatrices(trackK);
cout << " *** Position adjustment 1 " << hit->GetZ() << " "
<< (*trackK->fSteps)[i] << endl;
- AliFatal("Position adjustment 1.");
+ AliFatal(" Position adjustment 1.");
}
else break;
}
hit->SetNTrackHits(hit->GetNTrackHits()+1); // mark hit as being on track
if (ichamb == 9) {
// the last chamber
- //AZ(z->-z) trackK->fTrackDir = -1;
trackK->fTrackDir = 1;
trackK->fBPFlag = kTRUE;
}
fChi2 += dChi2Tmp; // Chi2
fPosition = savePosition;
// Add filtered matrices for the smoother
- //AZ(z->-z) if (fTrackDir < 0) {
if (fTrackDir > 0) {
for (Int_t i=fNSteps-1; i>=0; i--) {
if (TMath::Abs(fPosition-(*fSteps)[i]) > 0.1) {
//if (TMath::Abs(fPosition-fSteps[i]) > 0.1) {
RemoveMatrices(this);
- if (fgDebug > 1) cout << " *** Position adjustment 2 " << fPosition << " " << (*fSteps)[i] << endl;
+ if (fgDebug > 0) cout << " *** Position adjustment 2 " << fPosition << " " << (*fSteps)[i] << endl;
}
else break;
} // for (Int_t i=fNSteps-1;
// check whether the Invert method returns flag if matrix cannot be inverted,
// and do not calculate the Determinant in that case !!!!
if (wu.Determinant() != 0) {
- // wu.Invert();
- Int_t ifailWU;
- mnvertLocalK(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifailWU);
+ Int_t ifail;
+ mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant wu=0:");
+ AliWarning(" Determinant wu=0:");
}
trackParTmp = TMatrixD(wu,TMatrixD::kMult,right);
// check whether the Invert method returns flag if matrix cannot be inverted,
// and do not calculate the Determinant in that case !!!!
if (fWeight->Determinant() != 0) {
- //fWeight->Invert(); // covariance
-
- Int_t ifailWeight;
- mnvertLocalK(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifailWeight);
+ Int_t ifail;
+ mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fWeight=0:");
+ AliWarning(" Determinant fWeight=0:");
}
cosAlph = TMath::Cos((*fTrackParNew)(2,0));
(*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
(*fWeight)(3,3) += sign*theta0*theta0; // beta
- //fWeight->Invert(); // weight
-
- Int_t ifailWeight;
- mnvertLocalK(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifailWeight);
+ Int_t ifail;
+ mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
return;
}
VGeant3[0] = (*fTrackParNew)(1,0); // X
VGeant3[1] = (*fTrackParNew)(0,0); // Y
VGeant3[2] = fPositionNew; // Z
- //AZ(z->-z) VGeant3[3] = iFB*TMath::Sin((*fTrackParNew)(3,0)); // Px/Ptot
- //AZ)z->-z) VGeant3[4] = iFB*TMath::Cos((*fTrackParNew)(3,0))*TMath::Sin((*fTrackParNew)(2,0)); // Py/Ptot
- VGeant3[3] = -iFB*TMath::Sin((*fTrackParNew)(3,0)); // Px/Ptot
- VGeant3[4] = -iFB*TMath::Cos((*fTrackParNew)(3,0))*TMath::Sin((*fTrackParNew)(2,0)); // Py/Ptot
+ VGeant3[3] = iFB*TMath::Sin((*fTrackParNew)(3,0)); // Px/Ptot
+ VGeant3[4] = iFB*TMath::Cos((*fTrackParNew)(3,0))*TMath::Sin((*fTrackParNew)(2,0)); // Py/Ptot
VGeant3[5] = iFB*TMath::Sqrt(1.0-VGeant3[3]*VGeant3[3]-VGeant3[4]*VGeant3[4]); // Pz/Ptot
VGeant3[6] = 1/TMath::Abs((*fTrackParNew)(4,0)); // Ptot
}
fPositionNew = VGeant3[2]; // Z
(*fTrackParNew)(0,0) = VGeant3[1]; // Y
(*fTrackParNew)(1,0) = VGeant3[0]; // X
- //AZ(z->-z) (*fTrackParNew)(3,0) = TMath::ASin(iFB*VGeant3[3]); // beta
- //AZ(z->-z) (*fTrackParNew)(2,0) = TMath::ASin(iFB*VGeant3[4]/TMath::Cos((*fTrackParNew)(3,0))); // alpha
- (*fTrackParNew)(3,0) = TMath::ASin(-iFB*VGeant3[3]); // beta
- (*fTrackParNew)(2,0) = TMath::ASin(iFB*VGeant3[4]/TMath::Cos((*fTrackParNew)(3,0))) - TMath::Pi(); // alpha
+ (*fTrackParNew)(3,0) = TMath::ASin(iFB*VGeant3[3]); // beta
+ (*fTrackParNew)(2,0) = TMath::ASin(iFB*VGeant3[4]/TMath::Cos((*fTrackParNew)(3,0))); // alpha
(*fTrackParNew)(4,0) = 1/VGeant3[6]*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0)); // 1/Ptot
}
SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHitsPtr)[i]))->GetZ());
}
((AliMUONTrack*)this)->AddTrackParamAtHit(&trackParam);
+ // Fill array of HitForRec's
+ ((AliMUONTrack*)this)->AddHitForRecAtHit((AliMUONHitForRec*)fTrackHitsPtr->UncheckedAt(i));
}
}
// Propagates track to the vertex thru absorber using Branson correction
// (makes use of the AliMUONTrackParam class)
- AliMUONTrackParam *trackParam = new AliMUONTrackParam();
+ //AliMUONTrackParam *trackParam = new AliMUONTrackParam();
+ AliMUONTrackParam trackParam = AliMUONTrackParam();
/*
trackParam->SetBendingCoor((*fTrackPar)(0,0));
trackParam->SetNonBendingCoor((*fTrackPar)(1,0));
trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0)));
trackParam->SetZ(fPosition);
*/
- SetTrackParam(trackParam, fTrackPar, fPosition);
+ SetTrackParam(&trackParam, fTrackPar, fPosition);
- trackParam->ExtrapToVertex(Double_t(0.), Double_t(0.), Double_t(0.));
+ 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;
+ (*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;
cout << 1/(*fTrackPar)(4,0) << " " << fPosition << " " << (*fTrackPar)(0,0) << endl;
// Get covariance matrix
*fCovariance = *fWeight;
if (fCovariance->Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fCovariance=0:");
+ AliWarning(" Determinant fCovariance=0:");
}
}
// from raw clusters
for (Int_t i1=0; i1<fNTrackHits; i1++) {
hit = (AliMUONHitForRec*) ((*fTrackHitsPtr)[i1]);
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- printf ("%4d", clus->GetTrack(1) - 1);
+ 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 ("%4d", clus->GetTrack(1));
}
cout << endl;
for (Int_t i1=0; i1<fNTrackHits; i1++) {
hit = (AliMUONHitForRec*) ((*fTrackHitsPtr)[i1]);
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- if (clus->GetTrack(2) != 0) printf ("%4d", clus->GetTrack(2) - 1);
+ 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 ("%4d", clus->GetTrack(2));
else printf ("%4s", " ");
}
}
/* Not needed - covariance matrix is not interesting to anybody
*fCovariance = *fWeight;
if (fCovariance->Determinant() != 0) {
- // fCovariance->Invert();
- Int_t ifailCov;
- mnvertLocalK(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fCovariance=0:" );
+ AliWarning(" Determinant fCovariance=0:" );
}
*/
}
*fCovariance = *fWeight;
if (fCovariance->Determinant() != 0) {
// fCovariance->Invert();
- Int_t ifailCov;
- mnvertLocalK(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fCovariance=0:" );
+ AliWarning(" Determinant fCovariance=0:" );
}
(*fCovariance)(0,0) += dl3*(dYdT*dYdT+dYdB*dYdB); // <yy>
(*fCovariance)(0,3) += dl2*dYdB; // <yb>
(*fCovariance)(3,0) = (*fCovariance)(0,3);
- (*fCovariance)(1,2) += dl2*(dXdT*dAdT+dXdB*dAdB); // <xa>
+ (*fCovariance)(1,2) += dl2*(-dXdT*dAdT+dXdB*dAdB); // <xa>
(*fCovariance)(2,1) = (*fCovariance)(1,2);
(*fCovariance)(1,3) += dl2*dXdB; // <xb>
// Get weight matrix
*fWeight = *fCovariance;
if (fWeight->Determinant() != 0) {
- // fWeight->Invert();
- Int_t ifailWeight;
- mnvertLocalK(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifailWeight);
+ Int_t ifail;
+ mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fWeight=0:");
+ AliWarning(" Determinant fWeight=0:");
}
}
}
printf("\n");
}
- Double_t chi2_max = 0;
+ Double_t chi2max = 0;
Int_t imax = 0;
for (Int_t i=0; i<fNTrackHits; i++) {
chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
- if (chi2 < chi2_max) continue;
- chi2_max = chi2;
+ if (chi2 < chi2max) continue;
+ chi2max = chi2;
imax = i;
}
- //if (chi2_max < 10) return kFALSE; // !!!
- //if (chi2_max < 25) imax = fNTrackHits - 1;
- if (chi2_max < 15) imax = fNTrackHits - 1; // discard the last point
+ //if (chi2max < 10) return kFALSE; // !!!
+ //if (chi2max < 25) imax = fNTrackHits - 1;
+ if (chi2max < 15) imax = fNTrackHits - 1; // discard the last point
// Check if the outlier is not from the seed segment
AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHitsPtr->UncheckedAt(imax);
if (skipHit == fStartSegment->GetHitForRec1() || skipHit == fStartSegment->GetHitForRec2()) {
imax0 = imax;
// Hits after the found one will be removed
- if (GetStation0() == 3 && skipHit->GetChamberNumber() > 7) {
+ if (GetStation0() == 3 && skipHit->GetChamberNumber() >= 7) {
SortHits(1, fTrackHitsPtr); // unsort hits
imax = fTrackHitsPtr->IndexOf(skipHit);
}
*((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
*(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
if (trackK->fWeight->Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fWeight=0:");
+ AliWarning(" Determinant fWeight=0:");
}
trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
trackK->fChi2 = 0;
// Covariance
*(trackK->fCovariance) = *(trackK->fWeight);
if (trackK->fCovariance->Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fCovariance=0:");
+ AliWarning(" Determinant fCovariance=0:");
}
iD = trackK->fCovFilter->Last()->GetUniqueID();
if (iD > 1) {
TMatrixD *matrix = (TMatrixD*) objArray->First();
TMatrixD *tmp = new TMatrixD(*matrix);
- objArray->AddAt(tmp,objArray->GetLast());
+ objArray->AddAtAndExpand(tmp,objArray->GetLast());
//cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
}
else ((TMatrixD*)(trackK->fJacob->RemoveAt(i)))->Delete();
}
- //__________________________________________________________________________
-void mnvertLocalK(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 */
-
//__________________________________________________________________________
Bool_t AliMUONTrackK::Smooth(void)
{
/*
for (Int_t i=fNSteps-1; i>=0; i--) cout << (*fSteps)[i] << " ";
cout << endl;
- for (Int_t i=fNSteps-1; i>=0; i--) {cout << i; ((TMatrixD*)fParFilter->UncheckedAt(i))->Print(); }
+ for (Int_t i=fNSteps-1; i>=0; i--) {cout << i << " " << (*fSteps)[i]; ((TMatrixD*)fParFilter->UncheckedAt(i))->Print(); ((TMatrixD*)fParExtrap->UncheckedAt(i))->Print(); ((TMatrixD*)fJacob->UncheckedAt(i))->Print();}
*/
for (Int_t i=ihit; i>=0; i--) cout << ((AliMUONHitForRec*)(*fTrackHitsPtr)[i])->GetZ() << " ";
cout << endl;
//parSmooth.Print(); ((TMatrixD*)(fParExtrap->Last()))->Print();
Bool_t found;
- Double_t chi2_max = 0;
+ Double_t chi2max = 0;
for (Int_t i=iLast+1; i>0; i--) {
if (i == iLast + 1) goto L33; // temporary fix
// Smoother gain matrix
weight = *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
if (weight.Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant weight=0:");
+ AliWarning(" Determinant weight=0:");
}
// Fk'Wkk+1
tmp = TMatrixD(*((TMatrixD*)(fJacob->UncheckedAt(i))),TMatrixD::kTransposeMult,weight);
// Calculate Chi2 of smoothed residuals
if (tmp.Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant tmp=0:");
+ AliWarning(" Determinant tmp=0:");
}
TMatrixD vector(tmp,TMatrixD::kMult,tmpPar);
TMatrixD chi2(tmpPar,TMatrixD::kTransposeMult,vector);
if (fgDebug > 1) chi2.Print();
(*fChi2Smooth)[ihit] = chi2(0,0);
- if (chi2(0,0) > chi2_max) chi2_max = chi2(0,0);
+ if (chi2(0,0) > chi2max) chi2max = chi2(0,0);
fChi2 += chi2(0,0);
if (chi2(0,0) < 0) {
- chi2.Print(); cout << i << " " << iLast << endl;
- AliFatal("chi2 < 0.");
+ chi2.Print(); cout << " chi2 < 0 " << i << " " << iLast << endl;
+ //AliFatal("chi2 < 0.");
}
// Save smoothed parameters
TMatrixD *par = new TMatrixD(parSmooth);
- fParSmooth->AddAt(par, ihit);
+ fParSmooth->AddAtAndExpand(par, ihit);
} // for (Int_t i=iLast+1;
- //if (chi2_max > 16) {
- //if (chi2_max > 25) {
- //if (chi2_max > 50) {
- //if (chi2_max > 100) {
- if (chi2_max > fgkChi2max) {
+ //if (chi2max > 16) {
+ //if (chi2max > 25) {
+ //if (chi2max > 50) {
+ //if (chi2max > 100) {
+ if (chi2max > fgkChi2max) {
//if (Recover()) DropBranches();
//Recover();
Outlier();
printf("\n");
}
- Double_t chi2_max = 0;
+ Double_t chi2max = 0;
Int_t imax = 0;
for (Int_t i=0; i<fNTrackHits; i++) {
- if ((*fChi2Smooth)[i] < chi2_max) continue;
- chi2_max = (*fChi2Smooth)[i];
+ if ((*fChi2Smooth)[i] < chi2max) continue;
+ chi2max = (*fChi2Smooth)[i];
imax = i;
}
// Check if the outlier is not from the seed segment
trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
*trackK = *this;
fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
- //AZ(z->-z) trackK->fTrackDir = -1;
trackK->fTrackDir = 1;
trackK->fRecover = 2;
trackK->fSkipHit = hit;
*((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
*(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
if (trackK->fWeight->Determinant() != 0) {
- Int_t ifailCov;
- mnvertLocalK(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifailCov);
+ Int_t ifail;
+ mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
} else {
- AliWarning("Determinant fWeight=0:");
+ AliWarning(" Determinant fWeight=0:");
}
trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
trackK->fChi2 = 0;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff