/*
$Log$
+Revision 1.19 2000/11/20 11:24:10 gosset
+New package for reconstructed tracks (A. Gheata):
+* tree output in the file "tree_reco.root"
+* display events and make histograms from this file
+
Revision 1.18 2000/10/26 12:47:03 gosset
Real distance between chambers of each station taken into account
for the reconstruction parameters "fSegmentMaxDistBending[5]"
//************* Defaults parameters for reconstruction
static const Double_t kDefaultMinBendingMomentum = 3.0;
+static const Double_t kDefaultMaxBendingMomentum = 500.0;
+static const Double_t kDefaultMaxChi2 = 100.0;
static const Double_t kDefaultMaxSigma2Distance = 16.0;
static const Double_t kDefaultBendingResolution = 0.01;
static const Double_t kDefaultNonBendingResolution = 0.144;
// Set reconstruction parameters to default values
// Would be much more convenient with a structure (or class) ????
fMinBendingMomentum = kDefaultMinBendingMomentum;
+ fMaxBendingMomentum = kDefaultMaxBendingMomentum;
+ fMaxChi2 = kDefaultMaxChi2;
fMaxSigma2Distance = kDefaultMaxSigma2Distance;
AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
2.0 * TMath::Pi() / 180.0;
else fRMin[ch] = 30.0;
+ // maximum radius at 10 degrees and Z of chamber
+ fRMax[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
+ 10.0 * TMath::Pi() / 180.0;
}
- // maximum radius
- // like in TRACKF_STAT (10 degrees ????)
- fRMax[0] = fRMax[1] = 91.5;
- fRMax[2] = fRMax[3] = 122.5;
- fRMax[4] = fRMax[5] = 158.3;
- fRMax[6] = fRMax[7] = 260.0;
- fRMax[8] = fRMax[9] = 260.0;
// ******** Parameters for making segments
// should be parametrized ????
bendCoor = Hit->Y();
nonBendCoor = Hit->X();
radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor));
- if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
+ // This cut is not needed with a realistic chamber geometry !!!!
+// if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
// new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's
hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit);
fNHitsForRec++;
hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution));
hitForRec->SetNonBendingCoor(nonBendCoor
+ gRandom->Gaus(0., fNonBendingResolution));
+// // !!!! without smearing
+// hitForRec->SetBendingCoor(bendCoor);
+// hitForRec->SetNonBendingCoor(nonBendCoor);
// more information into HitForRec
// resolution: angular effect to be added here ????
hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
AliMUONSegment *segment;
Bool_t last2st;
Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
- impactParam = 0., maxImpactParam = 0.; // =0 to avoid compilation warnings.
+ impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings.
AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
if (fPrintLevel >= 1)
cout << "enter MakeSegmentsPerStation (0...) " << Station << endl;
// maximum impact parameter (cm) according to fMinBendingMomentum...
maxImpactParam =
TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
+ // minimum impact parameter (cm) according to fMaxBendingMomentum...
+ minImpactParam =
+ TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum));
}
else last2st = kFALSE;
// extrapolation factor from Z of first chamber to Z of second chamber
// Conditions "distBend" and "impactParam" correlated for these stations ????
if ((distBend < fSegmentMaxDistBending[Station]) &&
(distNonBend < fSegmentMaxDistNonBending[Station]) &&
- (!last2st || (impactParam < maxImpactParam))) {
+ (!last2st || (impactParam < maxImpactParam)) &&
+ (!last2st || (impactParam > minImpactParam))) {
// make new segment
segment = new ((*fSegmentsPtr[Station])[segmentIndex])
AliMUONSegment(hit1Ptr, hit2Ptr);
{
// Follow tracks in stations(1..) 3, 2 and 1
// too long: should be made more modular !!!!
- AliMUONHitForRec *bestHit, *extrapHit, *hit;
- AliMUONSegment *bestSegment, *extrapSegment, *segment;
+ AliMUONHitForRec *bestHit, *extrapHit, *extrapCorrHit, *hit;
+ AliMUONSegment *bestSegment, *extrapSegment, *extrapCorrSegment, *segment;
AliMUONTrack *track, *nextTrack;
AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
// -1 to avoid compilation warnings
Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex;
Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
+ Double_t bendingMomentum, chi2Norm = 0.;
AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
// local maxSigma2Distance, for easy increase in testing
maxSigma2Distance = fMaxSigma2Distance;
// extrapolation to station
trackParam1->ExtrapToStation(station, trackParam);
extrapSegment = new AliMUONSegment(); // empty segment
+ extrapCorrSegment = new AliMUONSegment(); // empty corrected segment
// multiple scattering factor corresponding to one chamber
// and momentum in bending plane (not total)
mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
extrapSegment->UpdateFromStationTrackParam
(trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
trackParam1->GetInverseBendingMomentum());
+ // same thing for corrected segment
+ // better to use copy constructor, after checking that it works properly !!!!
+ extrapCorrSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
+ extrapCorrSegment->
+ SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
+ extrapCorrSegment->UpdateFromStationTrackParam
+ (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
+ trackParam1->GetInverseBendingMomentum());
bestChi2 = 5.0;
bestSegment = NULL;
if (fPrintLevel >= 10) {
// multiple scattering ????
// separation in 2 functions: Segment and HitForRec ????
segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
- chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
+ // correction of corrected segment (fBendingCoor and fNonBendingCoor)
+ // according to real Z value of "segment" and slopes of "extrapSegment"
+ extrapCorrSegment->
+ SetBendingCoor(extrapSegment->GetBendingCoor() +
+ extrapSegment->GetBendingSlope() *
+ (segment->GetHitForRec1()->GetZ() -
+ (&(pMUON->Chamber(2 * station)))->Z()));
+ extrapCorrSegment->
+ SetNonBendingCoor(extrapSegment->GetNonBendingCoor() +
+ extrapSegment->GetNonBendingSlope() *
+ (segment->GetHitForRec1()->GetZ() -
+ (&(pMUON->Chamber(2 * station)))->Z()));
+ chi2 = segment->
+ NormalizedChi2WithSegment(extrapCorrSegment, maxSigma2Distance);
if (chi2 < bestChi2) {
// update best Chi2 and Segment if better found
bestSegment = segment;
// should consider all possibilities ????
// multiple scattering ???? do about like for extrapSegment !!!!
extrapHit = new AliMUONHitForRec(); // empty hit
+ extrapCorrHit = new AliMUONHitForRec(); // empty corrected hit
bestChi2 = 3.0;
bestHit = NULL;
if (fPrintLevel >= 10) {
// resolutions from "extrapSegment"
extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
+ // same things for corrected hit
+ // better to use copy constructor, after checking that it works properly !!!!
+ extrapCorrHit->
+ SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
+ extrapCorrHit->
+ SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
+ extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
+ extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
// Loop over hits in the chamber
ch = 2 * station + chInStation;
for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
fNHitsForRecPerChamber[ch];
iHit++) {
hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
+ // correction of corrected hit (fBendingCoor and fNonBendingCoor)
+ // according to real Z value of "hit" and slopes of right "trackParam"
+ extrapCorrHit->
+ SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor() +
+ (&(trackParam[chInStation]))->GetBendingSlope() *
+ (hit->GetZ() -
+ (&(trackParam[chInStation]))->GetZ()));
+ extrapCorrHit->
+ SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor() +
+ (&(trackParam[chInStation]))->GetNonBendingSlope() *
+ (hit->GetZ() -
+ (&(trackParam[chInStation]))->GetZ()));
// condition for hit not already in segment ????
chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
if (chi2 < bestChi2) {
// and corresponding TrackHit's, ...
track->Remove();
delete extrapSegment;
+ delete extrapCorrSegment;
+ delete extrapHit;
+ delete extrapCorrHit;
break; // stop the search for this candidate:
// exit from the loop over station
}
+ delete extrapHit;
+ delete extrapCorrHit;
}
delete extrapSegment;
+ delete extrapCorrSegment;
// Sort track hits according to increasing Z
track->GetTrackHitsPtr()->Sort();
// Update track parameters at first track hit (smallest Z)
trackParam1 = ((AliMUONTrackHit*)
(track->GetTrackHitsPtr()->First()))->GetTrackParam();
+ bendingMomentum = 0.;
+ if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.)
+ bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum()));
+ // Track removed if bendingMomentum not in window [min, max]
+ if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) {
+ track->Remove();
+ break; // stop the search for this candidate:
+ // exit from the loop over station
+ }
// Track fit
// with multiple Coulomb scattering if all stations
if (station == 0) track->SetFitMCS(1);
track->SetFitNParam(5); // with 5 parameters (momentum and position)
track->SetFitStart(1); // from parameters at first hit
track->Fit();
+ Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
+ if (numberOfDegFree > 0) {
+ chi2Norm = track->GetFitFMin() / numberOfDegFree;
+ } else {
+ chi2Norm = 1.e10;
+ }
+ // Track removed if normalized chi2 too high
+ if (chi2Norm > fMaxChi2) {
+ track->Remove();
+ break; // stop the search for this candidate:
+ // exit from the loop over station
+ }
if (fPrintLevel >= 10) {
cout << "FollowTracks: track candidate(0..): " << trackIndex
<< " after fit from station(0..): " << station << " to 4" << endl;
/*
$Log$
+Revision 1.4 2000/06/30 10:15:48 gosset
+Changes to EventReconstructor...:
+precision fit with multiple Coulomb scattering;
+extrapolation to vertex with Branson correction in absorber (JPC)
+
Revision 1.3 2000/06/25 13:06:39 hristov
Inline functions moved from *.cxx to *.h files instead of forward declarations
// The "road" is parametrized from the old reco_muon.F
// with 8 cm between stations.
AliMUONTrackParam *param0;
- Double_t cReso2, sReso2;
+// Double_t cReso2, sReso2;
// parameters to define the widths of the searching roads in station 0,1,2
// width = p0 + p1/ (momentum)^2
// station number: 0 1 2
- static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 };
- static Double_t p1BendingCoor[3] = { 986., 821., 446. };
- static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 };
- static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 };
- static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 };
- static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. };
- static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };
- static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
+// static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 };
+// static Double_t p1BendingCoor[3] = { 986., 821., 446. };
+// static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 };
+// static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 };
+// static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 };
+// static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. };
+// static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };
+// static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
+
+ static Double_t p0BendingCoor[3] = { 6.43e-2, 6.43e-2, 6.43e-2 };
+ static Double_t p1BendingCoor[3] = { 986., 986., 986. };
+ static Double_t p0BendingSlope[3] = { 3.6e-6, 3.6e-6, 3.6e-6 };
+ static Double_t p1BendingSlope[3] = { 1.1e-2, 1.1e-2, 1.1e-2 };
+ static Double_t p0NonBendingCoor[3] = { 0.049, 0.049, 0.049 };
+ static Double_t p1NonBendingCoor[3] = { 1444., 1444., 1444. };
+ static Double_t p0NonBendingSlope[3] = { 6.8e-4, 6.8e-4, 6.8e-4 };
+ static Double_t p1NonBendingSlope[3] = { 0., 0., 0. };
param0 = &(TrackParam[0]);
// OLD version
// because they are added in the functions "NormalizedChi2WithSegment"
// and "NormalizedChi2WithHitForRec"
// Bending plane
- cReso2 = fBendingCoorReso2;
- sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
- fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
- fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
+// cReso2 = fBendingCoorReso2;
+// sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
+ fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ; // - cReso2
+ fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
// Non bending plane
- cReso2 = fNonBendingCoorReso2;
- sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
- fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
- fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
+// cReso2 = fNonBendingCoorReso2;
+// sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
+ fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2;
+ fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2;
return;
}
/*
$Log$
+Revision 1.7 2000/09/19 15:50:46 gosset
+TrackChi2MCS function: covariance matrix better calculated,
+taking into account missing planes...
+
Revision 1.6 2000/07/20 12:45:27 gosset
New "EventReconstructor..." structure,
hopefully more adapted to tree/streamer.
#include <TClonesArray.h>
#include <TMath.h>
-#include <TMatrix.h>
+#include <TMatrixD.h>
#include <TObjArray.h>
#include <TVirtualFitter.h>
// choice of function to be minimized according to fFitMCS
if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
else fgFitter->SetFCN(TrackChi2MCS);
- arg[0] = 1;
+ arg[0] = -1;
fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
// Parameters according to "fFitStart"
// (should be a function to be used at every place where needed ????)
trackParam->GetNonBendingSlope(),
0.001, -0.5, 0.5);
if (fFitNParam == 5) {
- // set last 2 Minuit parameters (no limits for the search: min=max=0)
+ // set last 2 Minuit parameters
+ // mandatory limits in Bending to avoid NaN values of parameters
fgFitter->SetParameter(3, "X",
trackParam->GetNonBendingCoor(),
- 0.03, 0.0, 0.0);
+ 0.03, -500.0, 500.0);
+ // mandatory limits in non Bending to avoid NaN values of parameters
fgFitter->SetParameter(4, "Y",
trackParam->GetBendingCoor(),
- 0.10, 0.0, 0.0);
+ 0.10, -500.0, 500.0);
}
// search without gradient calculation in the function
fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
Double_t hbc1, hbc2, pbc1, pbc2;
Double_t hnbc1, hnbc2, pnbc1, pnbc2;
Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
- TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
- TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
+ TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
+ TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
Double_t *msa2 = new Double_t[numberOfHit];
// Predicted coordinates and multiple scattering angles are first calculated
}
} // for (hitNumber2 = hitNumber1;...
} // for (hitNumber1 = 0;...
-
+ // Normalization of covariance matrices
+ Double_t normCovBending2 = covBending->E2Norm();
+ Double_t normCovNonBending2 = covNonBending->E2Norm();
+ (*covBending) *= 1/normCovBending2;
+ (*covNonBending) *= 1/normCovNonBending2;
+// if (covBending->Determinant() < 1.e-33) {
+// printf(" *** covBending *** \n");
+// covBending->Print();
+// printf(" *** covNonBending *** \n");
+// covNonBending->Print();
+// cout << " number of hits " << numberOfHit << endl;
+// cout << "Momentum = " << 1/Param[0] <<endl;
+// cout << "normCovBending = " << normCovBending2 << endl;
+// cout << "normCovNonBending = " << normCovNonBending2 << endl;
+// exit(0);
+
+// }
// Inverts covariance matrix
goodDeterminant = kTRUE;
// check whether the Invert method returns flag if matrix cannot be inverted,
// Calculates Chi2
if (goodDeterminant) {
// with Multiple Scattering if inversion correct
+ // Inverse matrices without normalization
+ (*covBending) *= 1/normCovBending2;
+ (*covNonBending) *= 1/normCovNonBending2;
for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();