+ if (ExistDouble()) { delete hits; return kFALSE; }
+
+ //DropBranches(imax0, hits); // drop branches downstream the discarded hit
+ delete hits;
+
+ nRecTracks = fgTrackReconstructor->GetNRecTracks();
+ skipHit = (AliMUONHitForRec*) ((*fTrackHits)[fNmbTrackHits-1]);
+ // Remove all saved steps and smoother matrices after the skipped hit
+ RemoveMatrices(skipHit->GetZ());
+
+ //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);
+ fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
+ trackK->fRecover = 1;
+ trackK->fSkipHit = skipHit;
+ trackK->fNmbTrackHits = fNmbTrackHits;
+ delete trackK->fTrackHits; // not efficient ?
+ trackK->fTrackHits = new TObjArray(*fTrackHits);
+ if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
+ return kTRUE;
+ }
+
+ trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
+ *trackK = *this;
+ fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
+ //AZ(z->-z) trackK->fTrackDir = -1;
+ trackK->fTrackDir = 1;
+ trackK->fRecover = 1;
+ trackK->fSkipHit = skipHit;
+ Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
+ if (iD > 1) {
+ trackK->fParFilter->Last()->SetUniqueID(iD-1);
+ CreateMatrix(trackK->fParFilter);
+ }
+ *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
+ trackK->fParFilter->Last()->SetUniqueID(1);
+ *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
+ iD = trackK->fCovFilter->Last()->GetUniqueID();
+ if (iD > 1) {
+ trackK->fCovFilter->Last()->SetUniqueID(iD-1);
+ CreateMatrix(trackK->fCovFilter);
+ }
+ *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
+ trackK->fCovFilter->Last()->SetUniqueID(1);
+ *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
+ if (trackK->fWeight->Determinant() != 0) {
+ Int_t ifail;
+ mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
+ } else {
+ AliWarning(" Determinant fWeight=0:");
+ }
+ trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
+ trackK->fChi2 = 0;
+ for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
+ if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
+ return kTRUE;
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackK::AddMatrices(AliMUONTrackK *trackK, Double_t dChi2, AliMUONHitForRec *hitAdd)
+{
+/// Adds matrices for the smoother and keep Chi2 for the point
+/// Track parameters
+ //trackK->fParFilter->Last()->Print();
+ Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
+ if (iD > 1) {
+ trackK->fParFilter->Last()->SetUniqueID(iD-1);
+ CreateMatrix(trackK->fParFilter);
+ iD = 1;
+ }
+ *((TMatrixD*)(trackK->fParFilter->Last())) = *(trackK->fTrackPar);
+ trackK->fParFilter->Last()->SetUniqueID(iD);
+ if (fgDebug > 1) {
+ cout << " Add matrices" << " " << fPosition << " " << dChi2 << " " << fgHitForRec->IndexOf(hitAdd) << endl;
+ //trackK->fTrackPar->Print();
+ //trackK->fTrackParNew->Print();
+ trackK->fParFilter->Last()->Print();
+ cout << " Add matrices" << endl;
+ }
+ // Covariance
+ *(trackK->fCovariance) = *(trackK->fWeight);
+ if (trackK->fCovariance->Determinant() != 0) {
+ Int_t ifail;
+ mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
+ } else {
+ AliWarning(" Determinant fCovariance=0:");
+ }
+ iD = trackK->fCovFilter->Last()->GetUniqueID();
+ if (iD > 1) {
+ trackK->fCovFilter->Last()->SetUniqueID(iD-1);
+ CreateMatrix(trackK->fCovFilter);
+ iD = 1;
+ }
+ *((TMatrixD*)(trackK->fCovFilter->Last())) = *(trackK->fCovariance);
+ trackK->fCovFilter->Last()->SetUniqueID(iD);
+
+ // Save Chi2-increment for point
+ Int_t indx = trackK->fTrackHits->IndexOf(hitAdd);
+ if (indx < 0) indx = fNmbTrackHits;
+ if (trackK->fChi2Array->fN <= indx) trackK->fChi2Array->Set(indx+10);
+ trackK->fChi2Array->AddAt(dChi2,indx);
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackK::CreateMatrix(TObjArray *objArray) const
+{
+/// Create new matrix and add it to TObjArray
+
+ TMatrixD *matrix = (TMatrixD*) objArray->First();
+ TMatrixD *tmp = new TMatrixD(*matrix);
+ objArray->AddAtAndExpand(tmp,objArray->GetLast());
+ //cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackK::RemoveMatrices(Double_t zEnd)
+{
+/// Remove matrices (and saved steps) in the smoother part with abs(z) < abs(zEnd)
+
+ for (Int_t i=fNSteps-1; i>=0; i--) {
+ if (fgDebug > 1) printf("%10.4f %10.4f \n", (*fSteps)[i], zEnd);
+ if (TMath::Abs((*fSteps)[i]) > TMath::Abs(zEnd) - 0.01) break;
+ RemoveMatrices(this);
+ } // for (Int_t i=fNSteps-1;
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackK::RemoveMatrices(AliMUONTrackK* trackK)
+{
+/// Remove last saved matrices and steps in the smoother part
+
+ trackK->fNSteps--;
+ Int_t i = trackK->fNSteps;
+
+ Int_t id = 0;
+ // Delete only matrices not shared by several tracks
+ id = trackK->fParExtrap->Last()->GetUniqueID();
+ if (id > 1) {
+ trackK->fParExtrap->Last()->SetUniqueID(id-1);
+ trackK->fParExtrap->RemoveAt(i);
+ }
+ else ((TMatrixD*)(trackK->fParExtrap->RemoveAt(i)))->Delete();
+ id = fParFilter->Last()->GetUniqueID();
+ if (id > 1) {
+ trackK->fParFilter->Last()->SetUniqueID(id-1);
+ trackK->fParFilter->RemoveAt(i);
+ }
+ else ((TMatrixD*)(trackK->fParFilter->RemoveAt(i)))->Delete();
+ id = trackK->fCovExtrap->Last()->GetUniqueID();
+ if (id > 1) {
+ trackK->fCovExtrap->Last()->SetUniqueID(id-1);
+ trackK->fCovExtrap->RemoveAt(i);
+ }
+ else ((TMatrixD*)(trackK->fCovExtrap->RemoveAt(i)))->Delete();
+ id = trackK->fCovFilter->Last()->GetUniqueID();
+ if (id > 1) {
+ trackK->fCovFilter->Last()->SetUniqueID(id-1);
+ trackK->fCovFilter->RemoveAt(i);
+ }
+ else ((TMatrixD*)(trackK->fCovFilter->RemoveAt(i)))->Delete();
+ id = trackK->fJacob->Last()->GetUniqueID();
+ if (id > 1) {
+ trackK->fJacob->Last()->SetUniqueID(id-1);
+ trackK->fJacob->RemoveAt(i);
+ }
+ else ((TMatrixD*)(trackK->fJacob->RemoveAt(i)))->Delete();
+}
+
+ //__________________________________________________________________________
+Bool_t AliMUONTrackK::Smooth(void)
+{
+/// Apply smoother
+ Int_t ihit = fNmbTrackHits - 1;
+ AliMUONHitForRec *hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
+ fChi2Smooth = new TArrayD(fNmbTrackHits);
+ fChi2Smooth->Reset(-1);
+ fChi2 = 0;
+ fParSmooth = new TObjArray(15);
+ fParSmooth->Clear();
+
+ if (fgDebug > 0) {
+ cout << " ******** Enter Smooth " << endl;
+ cout << (*fSteps)[fNSteps-1] << " " << fPosition << " " << fPositionNew << endl;
+ /*
+ 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 << " " << (*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*)(*fTrackHits)[i])->GetZ() << " ";
+ cout << endl;
+ }
+
+ // Find last point corresponding to the last hit
+ Int_t iLast = fNSteps - 1;
+ for ( ; iLast>=0; iLast--) {
+ //AZ(z->-z) if ((*fSteps)[iLast] + 0.01 > ((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ()) break;
+ if (TMath::Abs((*fSteps)[iLast]) + 0.01 > TMath::Abs(((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ())) break;
+ }
+
+ TMatrixD parSmooth = *((TMatrixD*)(fParFilter->UncheckedAt(iLast))); // last filtered = first smoothed
+ //parSmooth.Dump();
+ TMatrixD covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(iLast))); // last filtered = first smoothed
+ TMatrixD tmp=covSmooth, weight=covSmooth, gain = covSmooth;
+ TMatrixD tmpPar = *fTrackPar;
+ //parSmooth.Print(); ((TMatrixD*)(fParExtrap->Last()))->Print();
+
+ Bool_t found;
+ 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 ifail;
+ mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
+ } else {
+ AliWarning(" Determinant weight=0:");