X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=MUON%2FAliMUONTrackK.cxx;h=3d624198201a2d1a5b102487c35feba24a4bd343;hb=184bcc16e057d9c6b9510ca1ff46bb7ef23bec45;hp=3a97361fbf94946c9b8dfdd008be62e7440e7066;hpb=54d7ba50397a146159ff9259b4de2e128aa5a7c4;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONTrackK.cxx b/MUON/AliMUONTrackK.cxx index 3a97361fbf9..3d624198201 100644 --- a/MUON/AliMUONTrackK.cxx +++ b/MUON/AliMUONTrackK.cxx @@ -18,45 +18,47 @@ // --------------------- // 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 -#include -#include -#include - #include "AliMUONTrackK.h" -#include "AliMUON.h" +#include "AliMUONData.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 +#include +#include + +/// \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() : AliMUONTrack(), @@ -95,7 +97,7 @@ AliMUONTrackK::AliMUONTrackK() } //__________________________________________________________________________ -AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClonesArray *hitForRec) +AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructorK *TrackReconstructor, TClonesArray *hitForRec) : AliMUONTrack(), fStartSegment(0x0), fPosition(0.), @@ -133,10 +135,9 @@ AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClo } //__________________________________________________________________________ -AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment) - //: AliMUONTrack(segment, segment, fgTrackReconstructor) - : AliMUONTrack(NULL, segment, fgTrackReconstructor), - fStartSegment(segment), +AliMUONTrackK::AliMUONTrackK(AliMUONObjectPair *segment) + : AliMUONTrack(NULL,NULL), + fStartSegment(new AliMUONObjectPair(*segment)), fPosition(0.), fPositionNew(0.), fChi2(0.), @@ -162,20 +163,20 @@ AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment) 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; // 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(); } // Fill array of track parameters @@ -201,32 +202,32 @@ AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment) 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); 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; + } //__________________________________________________________________________ @@ -234,6 +235,11 @@ AliMUONTrackK::~AliMUONTrackK() { /// Destructor + if (fStartSegment) { + delete fStartSegment; + fStartSegment = 0x0; + } + if (fTrackHits) { //cout << fNmbTrackHits << endl; for (Int_t i = 0; i < fNmbTrackHits; i++) { @@ -323,8 +329,11 @@ 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; @@ -442,7 +451,7 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back, } 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) { @@ -563,17 +572,15 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back, 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 @@ -670,6 +677,7 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back, } // 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; } @@ -699,7 +707,7 @@ void AliMUONTrackK::ParPropagation(Double_t zEnd) 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); @@ -719,7 +727,7 @@ void AliMUONTrackK::ParPropagation(Double_t zEnd) 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; @@ -816,7 +824,8 @@ void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth) 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; } @@ -832,6 +841,19 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int 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) { @@ -840,10 +862,15 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int 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; @@ -882,7 +909,7 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int iDhit = 1; } - + TArrayD branchChi2(20); for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) { if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem) hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]); @@ -935,7 +962,7 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int // 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())); @@ -974,7 +1001,8 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int 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(); @@ -982,7 +1010,7 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int 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; @@ -1020,12 +1048,10 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int 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 @@ -1037,6 +1063,8 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int trackK->fTrackDir = 1; trackK->fBPFlag = kTRUE; } + if (nHitsOK > branchChi2.GetSize()) branchChi2.Set(branchChi2.GetSize()+10); + branchChi2[nHitsOK-1] = dChi2; } } } @@ -1067,10 +1095,41 @@ Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int } */ } // 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) { @@ -1129,7 +1188,7 @@ void AliMUONTrackK::MSThin(Int_t sign) 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 @@ -1160,7 +1219,7 @@ void AliMUONTrackK::StartBack(void) //__________________________________________________________________________ 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; @@ -1284,26 +1343,22 @@ void AliMUONTrackK::FillMUONTrack(void) { /// 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)); } @@ -1322,288 +1377,6 @@ void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, 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; i1GetChamberNumber()); - } - cout << endl; - if (fgDebug > 0) { - for (Int_t i1=0; i1GetHitNumber() << " "; - //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " "; - printf ("%5d", fgHitForRec->IndexOf(hit)); - } - cout << endl; - } - if (fgTrackReconstructor->GetRecTrackRefHits()) { - // from track ref. hits - for (Int_t i1=0; i1GetTTRTrack() + hit->GetTrackRefSignal()*10000 << " "; - } - } else { - // from raw clusters - for (Int_t i1=0; i1GetHitNumber() < 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; i1GetHitNumber() < 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; i1GetHitNumber() << " "; - cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " "; - //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " "; - } - cout << endl; - for (Int_t i1=0; i1Determinant() != 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) { @@ -1714,7 +1487,7 @@ Bool_t AliMUONTrackK::Recover(void) 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 } @@ -1748,8 +1521,8 @@ Bool_t AliMUONTrackK::Recover(void) // 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); @@ -2067,7 +1840,7 @@ void AliMUONTrackK::Outlier() } // 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; @@ -2162,40 +1935,24 @@ void AliMUONTrackK::Print(FILE *lun) const Int_t flag = 1; AliMUONHitForRec *hit = 0; - if (fgTrackReconstructor->GetRecTrackRefHits()) { - // from track ref. hits - for (Int_t j=0; jUncheckedAt(j); - if (hit->GetTTRTrack() > 1) { flag = 0; break; } - } - for (Int_t j=0; 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; i1GetMUONData()->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; i1GetMUONData()->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; i1GetTrack(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])); @@ -2265,12 +2022,13 @@ void AliMUONTrackK::DropBranches(Int_t imax, TObjArray *hits) } //__________________________________________________________________________ -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(); @@ -2278,9 +2036,11 @@ void AliMUONTrackK::DropBranches(AliMUONSegment *segment) for (Int_t i=icand+1; ifStartSegment; 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; } @@ -2298,7 +2058,7 @@ AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void) Int_t AliMUONTrackK::GetStation0(void) { /// Return seed station number - return fStartSegment->GetHitForRec1()->GetChamberNumber() / 2; + return ((AliMUONHitForRec*) fStartSegment->First())->GetChamberNumber() / 2; } //__________________________________________________________________________ @@ -2392,3 +2152,103 @@ Bool_t AliMUONTrackK::ExistDouble(void) 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 */ +