/*
$Log$
+Revision 1.4 2000/06/27 14:11:36 gosset
+Corrections against violations of coding conventions
+
Revision 1.3 2000/06/16 07:27:08 gosset
To remove problem in running RuleChecker, like in MUON-dev
#include "AliMUONChamber.h"
#include "AliMUONTrackHit.h"
#include "AliRun.h"
+#include "TParticle.h"
#ifndef WIN32
# define initfield initfield_
// pointer to segment
endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
// test compatibility between current segment and "extrapSegment"
+ // 4 because 4 quantities in chi2
if ((endSegment->
NormalizedChi2WithSegment(extrapSegment,
fMaxSigma2Distance)) <= 4.0) {
// pointer to HitForRec
hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
// test compatibility between current HitForRec and "extrapHitForRec"
+ // 2 because 2 quantities in chi2
if ((hit->
NormalizedChi2WithHitForRec(extrapHitForRec,
fMaxSigma2Distance)) <= 2.0) {
void AliMUONEventReconstructor::FollowTracks(void)
{
// Follow tracks in stations(1..) 3, 2 and 1
+ // too long: should be made more modular !!!!
AliMUONHitForRec *bestHit, *extrapHit, *hit;
AliMUONSegment *bestSegment, *extrapSegment, *segment;
- AliMUONTrack *track;
- AliMUONTrackParam *trackParam1, trackParam[2];
+ AliMUONTrack *track, *nextTrack;
+ AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
Int_t ch, chBestHit, iHit, iSegment, station, trackIndex;
- Double_t bestChi2, chi2, dZ1, dZ2, maxSigma2Distance, mcsFactor;
+ Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
- maxSigma2Distance = 4. * fMaxSigma2Distance; // sigma2cut increased by 4 !!!!
+ // local maxSigma2Distance, for easy increase in testing
+ maxSigma2Distance = fMaxSigma2Distance;
if (fPrintLevel >= 2)
cout << "enter FollowTracks" << endl;
// Loop over track candidates
- for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
+ track = (AliMUONTrack*) fRecTracksPtr->First();
+ trackIndex = -1;
+ while (track) {
+ // Follow function for each track candidate ????
+ trackIndex++;
+ nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
if (fPrintLevel >= 2)
cout << "FollowTracks: track candidate(0..): " << trackIndex << endl;
- // function for each track candidate ????
- track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
// Fit track candidate from vertex at X = Y = 0
+ track->SetFitMCS(0); // fit without Multiple Scattering
track->Fit(track->GetTrackParamAtVertex(), 3);
if (fPrintLevel >= 10) {
cout << "FollowTracks: track candidate(0..): " << trackIndex
track->RecursiveDump();
}
// Loop over stations(1..) 3, 2 and 1
- // something SPECIAL for stations 2 and 1 for majority coincidence ????
+ // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
+ // otherwise: majority coincidence 2 !!!!
for (station = 2; station >= 0; station--) {
// Track parameters at first track hit (smallest Z)
trackParam1 = ((AliMUONTrackHit*)
// Z difference between the two previous stations
dZ2 = (&(pMUON->Chamber(2 * station + 2)))->Z() -
(&(pMUON->Chamber(2 * station + 4)))->Z();
- extrapSegment->SetBendingCoorReso2(fBendingResolution);
- extrapSegment->SetNonBendingCoorReso2(fNonBendingResolution);
- extrapSegment->UpdateFromStationTrackParam(trackParam, mcsFactor, dZ1, dZ2);
+ // Z difference between the two chambers in the previous station
+ dZ3 = (&(pMUON->Chamber(2 * station)))->Z() -
+ (&(pMUON->Chamber(2 * station + 1)))->Z();
+ extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
+ extrapSegment->
+ SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
+ extrapSegment->UpdateFromStationTrackParam
+ (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
+ trackParam1->GetInverseBendingMomentum());
bestChi2 = 5.0;
bestSegment = NULL;
if (fPrintLevel >= 10) {
}
}
else {
- fRecTracksPtr->RemoveAt(trackIndex); // Remove candidate
+ // Remove current track candidate and update fNRecTracks
+ // To be checked: recursive delete of TrackHit's !!!!
+ // For cleaner implementation: call track->Remove()
+ // to be coded, with all cleanings,
+ // including links between HitForRec's and TrackHit's !!!!
+ fRecTracksPtr->Remove(track);
+ fNRecTracks--;
+ delete extrapSegment;
break; // stop the search for this candidate:
// exit from the loop over station
}
}
+ delete extrapSegment;
// 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();
- // Track fit from first track hit varying X and Y
+ // Track fit from first track hit varying X and Y,
+ // with multiple Coulomb scattering if all stations
+ if (station == 0) track->SetFitMCS(1);
+ else track->SetFitMCS(0);
track->Fit(trackParam1, 5);
if (fPrintLevel >= 10) {
cout << "FollowTracks: track candidate(0..): " << trackIndex
<< " after fit from station(0..): " << station << " to 4" << endl;
track->RecursiveDump();
}
- delete extrapSegment;
+ // Track extrapolation to the vertex through the absorber (Branson)
+ if (station == 0) {
+ trackParamVertex = *trackParam1;
+ (&trackParamVertex)->ExtrapToVertex();
+ track->SetTrackParamAtVertex(&trackParamVertex);
+ }
+ if (fPrintLevel >= 1) {
+ cout << "FollowTracks: track candidate(0..): " << trackIndex
+ << " after fit from station(0..): " << station << " to 4" << endl;
+ // track->RecursiveDump(); // CCC
+ }
} // for (station = 2;...
- } // for (trackIndex = 0;...
+ // go really to next track
+ track = nextTrack;
+ } // while (track)
+ // Compression of track array (necessary after Remove ????)
+ fRecTracksPtr->Compress();
+ return;
+}
+
+ //__________________________________________________________________________
+void AliMUONEventReconstructor::EventDump(void)
+{
+ // Dump reconstructed event (track parameters at vertex and at first hit),
+ // and the particle parameters
+
+ AliMUONTrack *track;
+ AliMUONTrackParam *trackParam, *trackParam1;
+ TClonesArray *particles; // pointer to the particle list
+ TParticle *p;
+ Double_t bendingSlope, nonBendingSlope, pYZ;
+ Double_t pX, pY, pZ, x, y, z, c;
+ Int_t np, trackIndex, nTrackHits;
+
+ if (fPrintLevel >= 1) cout << "****** enter EventDump ******" << endl;
+ if (fPrintLevel >= 1) {
+ cout << " Number of Reconstructed tracks :" << fNRecTracks << endl;
+ }
+ fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
+ // Loop over reconstructed tracks
+ for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
+ if (fPrintLevel >= 1)
+ cout << " track number: " << trackIndex << endl;
+ // function for each track for modularity ????
+ track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
+ nTrackHits = track->GetNTrackHits();
+ if (fPrintLevel >= 1)
+ cout << " number of track hits: " << nTrackHits << endl;
+ // track parameters at Vertex
+ trackParam = track->GetTrackParamAtVertex();
+ x = trackParam->GetNonBendingCoor();
+ y = trackParam->GetBendingCoor();
+ z = trackParam->GetZ();
+ bendingSlope = trackParam->GetBendingSlope();
+ nonBendingSlope = trackParam->GetNonBendingSlope();
+ pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
+ pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
+ pX = pZ * nonBendingSlope;
+ pY = pZ * bendingSlope;
+ c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
+ if (fPrintLevel >= 1)
+ printf(" track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
+ z, x, y, pX, pY, pZ, c);
+
+ // track parameters at first hit
+ trackParam1 = ((AliMUONTrackHit*)
+ (track->GetTrackHitsPtr()->First()))->GetTrackParam();
+ x = trackParam1->GetNonBendingCoor();
+ y = trackParam1->GetBendingCoor();
+ z = trackParam1->GetZ();
+ bendingSlope = trackParam1->GetBendingSlope();
+ nonBendingSlope = trackParam1->GetNonBendingSlope();
+ pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
+ pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
+ pX = pZ * nonBendingSlope;
+ pY = pZ * bendingSlope;
+ c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
+ if (fPrintLevel >= 1)
+ printf(" track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
+ z, x, y, pX, pY, pZ, c);
+ }
+ // informations about generated particles
+ particles = gAlice->Particles();
+ np = particles->GetEntriesFast();
+ printf(" **** number of generated particles: %d \n", np);
+
+ for (Int_t iPart = 0; iPart < np; iPart++) {
+ p = (TParticle*) particles->UncheckedAt(iPart);
+ printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
+ iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
+ }
return;
}
Double_t GetImpactParamFromBendingMomentum(Double_t BendingMomentum);
Double_t GetBendingMomentumFromImpactParam(Double_t ImpactParam);
void EventReconstruct(void);
+ void EventDump(void); // dump reconstructed event
protected:
/*
$Log$
+Revision 1.3 2000/06/25 13:06:39 hristov
+Inline functions moved from *.cxx to *.h files instead of forward declarations
+
Revision 1.2 2000/06/15 07:58:48 morsch
Code from MUON-dev joined
}
//__________________________________________________________________________
-void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
+void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum)
{
// Fill data members with values calculated from the array of track parameters
// pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
// with one chamber for the coordinate, two chambers for the angle,
// due to the arrangement in stations.
// Resolution coming from:
- // coordinate in closest station at "Dz1",
+ // coordinate in closest station at "Dz1" from current "Station",
// slope between closest stations, with "Dz2" interval between them,
- // extrapolation over "Dz" from closest station.
+ // interval "Dz3" between chambers of closest station,
+ // extrapolation over "Dz1" from closest station,
+ // "InverseMomentum".
// When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
// are assumed to be filled
// with the variance on bending and non bending coordinates.
+ // The "road" is parametrized from the old reco_muon.F
+ // with 8 cm between stations.
AliMUONTrackParam *param0;
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. };
param0 = &(TrackParam[0]);
+
+// OLD version
+// // Bending plane
+// fBendingCoor = param0->GetBendingCoor(); // coordinate
+// fBendingSlope = param0->GetBendingSlope(); // slope
+// cReso2 = fBendingCoorReso2;
+// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+// fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+// fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+// // Non bending plane
+// fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
+// fNonBendingSlope = param0->GetNonBendingSlope(); // slope
+// cReso2 = fNonBendingCoorReso2;
+// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+// fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+// fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+
+ // Coordinate and slope
// Bending plane
fBendingCoor = param0->GetBendingCoor(); // coordinate
fBendingSlope = param0->GetBendingSlope(); // slope
- cReso2 = fBendingCoorReso2;
- sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
- fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
- fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
// Non bending plane
fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
fNonBendingSlope = param0->GetNonBendingSlope(); // slope
+
+ // Resolutions
+ // cReso2 and sReso2 have to be subtracted here from the parametrization
+ // 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;
+ // Non bending plane
cReso2 = fNonBendingCoorReso2;
- sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
- fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
- fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+ sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
+ fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
+ fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
return;
}
+
+// OLD function, with roads automatically calculated instead from being parametrized
+// kept because it would be a better solution,
+// if one can really find the right values.
+// //__________________________________________________________________________
+// void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
+// {
+// // Fill data members with values calculated from the array of track parameters
+// // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
+// // of the station, respectively).
+// // Multiple Coulomb scattering is taking into account with "MCSfactor"
+// // corresponding to one chamber,
+// // with one chamber for the coordinate, two chambers for the angle,
+// // due to the arrangement in stations.
+// // Resolution coming from:
+// // coordinate in closest station at "Dz1",
+// // slope between closest stations, with "Dz2" interval between them,
+// // extrapolation over "Dz" from closest station.
+// // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
+// // are assumed to be filled
+// // with the variance on bending and non bending coordinates.
+// AliMUONTrackParam *param0;
+// Double_t cReso2, sReso2;
+// param0 = &(TrackParam[0]);
+// // Bending plane
+// fBendingCoor = param0->GetBendingCoor(); // coordinate
+// fBendingSlope = param0->GetBendingSlope(); // slope
+// cReso2 = fBendingCoorReso2;
+// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+// fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+// fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+// // Non bending plane
+// fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
+// fNonBendingSlope = param0->GetNonBendingSlope(); // slope
+// cReso2 = fNonBendingCoorReso2;
+// sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+// fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+// fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+// return;
+// }
AliMUONSegment* CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor);
Double_t NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut);
AliMUONHitForRec* CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor);
- void UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2);
+ void UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum);
// What is necessary for sorting TClonesArray's; sufficient too ????
Bool_t IsSortable() const { return kTRUE; }
/*
$Log$
+Revision 1.3 2000/06/25 13:23:28 hristov
+stdlib.h needed for non-Linux compilation
+
Revision 1.2 2000/06/15 07:58:48 morsch
Code from MUON-dev joined
#include <TClonesArray.h>
#include <TMinuit.h>
+#include <TMath.h>
+#include <TMatrix.h>
#include "AliMUONEventReconstructor.h"
#include "AliMUONHitForRec.h"
#include <stdlib.h>
+// variables to be known from minimization functions
static AliMUONTrack *trackBeingFitted;
static AliMUONTrackParam *trackParamBeingFitted;
-// Function to be minimized with Minuit
+// Functions to be minimized with Minuit
void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
+void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
+
+Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
ClassImp(AliMUONTrack) // Class implementation in ROOT context
AddSegment(EndSegment); // add hits from EndSegment
fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
SetTrackParamAtVertex(); // set track parameters at vertex
+ fFitMCS = 0;
return;
}
AddHitForRec(HitForRec); // add HitForRec
fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
SetTrackParamAtVertex(); // set track parameters at vertex
+ fFitMCS = 0;
return;
}
return *this;
}
+ //__________________________________________________________________________
+void AliMUONTrack::SetFitMCS(Int_t FitMCS)
+{
+ // Set track fit option with or without multiple Coulomb scattering
+ // from "FitMCS" argument: 0 without, 1 with
+ fFitMCS = FitMCS;
+ // better implementation with enum(with, without) ????
+ return;
+}
+
// Inline functions for Get and Set: inline removed because it does not work !!!!
AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) {
// Get pointer to fTrackParamAtVertex
TMinuit *minuit = new TMinuit(5);
trackBeingFitted = this; // for the track to be known from the function to minimize
trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ????
+ // try to use TVirtualFitter to get this feature automatically !!!!
minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ????
// how to go faster ???? choice of Minuit parameters like EDM ????
- minuit->SetFCN(TrackChi2);
+ // choice of function to be minimized according to fFitMCS
+ if (fFitMCS == 0) minuit->SetFCN(TrackChi2);
+ else minuit->SetFCN(TrackChi2MCS);
minuit->SetPrintLevel(1); // More printing !!!!
- // set first 3 parameters (try with no limits for the search: min=max=0)
+ // set first 3 parameters
+ // could be tried with no limits for the search (min=max=0) ????
minuit->mnparm(0, "InvBenP",
TrackParam->GetInverseBendingMomentum(),
- 0.003, 0.0, 0.0, error);
+ 0.003, -0.4, 0.4, error);
minuit->mnparm(1, "BenS",
TrackParam->GetBendingSlope(),
- 0.001, 0.0, 0.0, error);
+ 0.001, -0.5, 0.5, error);
minuit->mnparm(2, "NonBenS",
TrackParam->GetNonBendingSlope(),
- 0.001, 0.0, 0.0, error);
+ 0.001, -0.5, 0.5, error);
if (NParam == 5) {
- // set last 2 parameters (try with no limits for the search: min=max=0)
+ // set last 2 parameters (no limits for the search: min=max=0)
minuit->mnparm(3, "X",
TrackParam->GetNonBendingCoor(),
0.03, 0.0, 0.0, error);
// and their current values in array pointed to by "Param".
// Assumes that the track hits are sorted according to increasing Z.
// Track parameters at each TrackHit are updated accordingly.
+ // Multiple Coulomb scattering is not taken into account
AliMUONTrackHit* hit;
AliMUONTrackParam param1;
Int_t hitNumber;
dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
}
}
+
+ //__________________________________________________________________________
+void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
+{
+ // Return the "Chi2" to be minimized with Minuit for track fitting,
+ // with "NParam" parameters
+ // and their current values in array pointed to by "Param".
+ // Assumes that the track hits are sorted according to increasing Z.
+ // Track parameters at each TrackHit are updated accordingly.
+ // Multiple Coulomb scattering is taken into account with covariance matrix.
+ AliMUONTrackParam param1;
+ Chi2 = 0.0; // initialize Chi2
+ // copy of track parameters to be fitted
+ param1 = *trackParamBeingFitted;
+ // Minuit parameters to be fitted into this copy
+ param1.SetInverseBendingMomentum(Param[0]);
+ param1.SetBendingSlope(Param[1]);
+ param1.SetNonBendingSlope(Param[2]);
+ if (NParam == 5) {
+ param1.SetNonBendingCoor(Param[3]);
+ param1.SetBendingCoor(Param[4]);
+ }
+
+ AliMUONTrackHit* hit, hit1, hit2, hit3;
+ Bool_t GoodDeterminant;
+ Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3;
+ Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10];
+ Double_t hitBendingCoor[10], hitNonBendingCoor[10];
+ Double_t hitBendingReso2[10], hitNonBendingReso2[10];
+ Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
+ TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
+ TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
+
+ // Predicted coordinates and multiple scattering angles are first calculated
+ for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
+ hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
+ zHit[hitNumber] = hit->GetHitForRecPtr()->GetZ();
+ // do something special if 2 hits with same Z ????
+ // security against infinite loop ????
+ (¶m1)->ExtrapToZ(zHit[hitNumber]); // extrapolation
+ hit->SetTrackParam(¶m1);
+ paramBendingCoor[hitNumber]= (¶m1)->GetBendingCoor();
+ paramNonBendingCoor[hitNumber]= (¶m1)->GetNonBendingCoor();
+ hitBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetBendingCoor();
+ hitNonBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingCoor();
+ hitBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetBendingReso2();
+ hitNonBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingReso2();
+ ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2
+ }
+
+ // Calculates the covariance matrix
+ for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {
+ for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
+ (*covBending)(hitNumber1, hitNumber2) = 0;
+ (*covBending)(hitNumber2, hitNumber1) = 0;
+ if (hitNumber1 == hitNumber2){ // diagonal elements
+ (*covBending)(hitNumber2, hitNumber1) =
+ (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
+ }
+ // Multiple Scattering... loop on upstream chambers ??
+ for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++){
+ (*covBending)(hitNumber2, hitNumber1) =
+ (*covBending)(hitNumber2, hitNumber1) +
+ ((zHit[hitNumber1] - zHit[hitNumber3]) *
+ (zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]);
+ }
+ (*covNonBending)(hitNumber1, hitNumber2) = 0;
+ (*covNonBending)(hitNumber2, hitNumber1) =
+ (*covBending)(hitNumber2, hitNumber1);
+ if (hitNumber1 == hitNumber2) { // diagonal elements
+ (*covNonBending)(hitNumber2, hitNumber1) =
+ (*covNonBending)(hitNumber2, hitNumber1) -
+ hitBendingReso2[hitNumber1] + hitNonBendingReso2[hitNumber1] ;
+ }
+ }
+ }
+
+ // Inverts covariance matrix
+ GoodDeterminant = kTRUE;
+ if (covBending->Determinant() != 0) {
+ covBending->Invert();
+ } else {
+ GoodDeterminant = kFALSE;
+ cout << "Warning in ChiMCS Determinant Bending=0: " << endl;
+ }
+ if (covNonBending->Determinant() != 0){
+ covNonBending->Invert();
+ } else {
+ GoodDeterminant = kFALSE;
+ cout << "Warning in ChiMCS Determinant non Bending=0: " << endl;
+ }
+
+ // Calculates Chi2
+ if (GoodDeterminant) { // with Multiple Scattering if inversion correct
+ for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){
+ for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){
+ Chi2 = Chi2 +
+ ((*covBending)(hitNumber2, hitNumber1) *
+ (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
+ (hitBendingCoor[hitNumber2] - paramBendingCoor[hitNumber2]));
+ Chi2 = Chi2 +
+ ((*covNonBending)(hitNumber2, hitNumber1) *
+ (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
+ (hitNonBendingCoor[hitNumber2] - paramNonBendingCoor[hitNumber2]));
+ }
+ }
+ } else { // without Multiple Scattering if inversion impossible
+ for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++) {
+ Chi2 = Chi2 +
+ ((hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
+ (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) /
+ hitBendingReso2[hitNumber1]);
+ Chi2 = Chi2 +
+ ((hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
+ (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) /
+ hitNonBendingReso2[hitNumber1]);
+ }
+ }
+
+ delete covBending;
+ delete covNonBending;
+}
+
+Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
+{
+ // Returns square of multiple Coulomb scattering angle
+ // at TrackHit pointed to by "TrackHit"
+ Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
+ Double_t varMultipleScatteringAngle;
+ AliMUONTrackParam *param = TrackHit->GetTrackParam();
+ slopeBending = param->GetBendingSlope();
+ slopeNonBending = param->GetNonBendingSlope();
+ // thickness in radiation length for the current track,
+ // taking local angle into account
+ radiationLength =
+ trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() *
+ TMath::Sqrt(1.0 +
+ slopeBending * slopeBending + slopeNonBending * slopeNonBending);
+ inverseBendingMomentum2 =
+ param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
+ inverseTotalMomentum2 =
+ inverseBendingMomentum2 * (1.0 + slopeBending*slopeBending) /
+ (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending);
+ varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
+ varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
+ varMultipleScatteringAngle * varMultipleScatteringAngle;
+ return varMultipleScatteringAngle;
+}
AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor); // Constructor from two Segment's
AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor); // Constructor from one Segment and one HitForRec
+ AliMUONEventReconstructor* GetEventReconstructor(void) {return fEventReconstructor;};
AliMUONTrackParam* GetTrackParamAtVertex(void);
- void SetTrackParamAtVertex(void);
- AliMUONTrackParam* GetTrackParamAtFirstHit(void);
+ void SetTrackParamAtVertex(void); // Set track parameters at vertex from last stations 4 & 5
+ void SetTrackParamAtVertex(AliMUONTrackParam* TrackParam) {fTrackParamAtVertex = *TrackParam;};
+
TClonesArray* GetTrackHitsPtr(void);
Int_t GetNTrackHits(void);
+ Int_t GetFitMCS(void) {return fFitMCS;};
+ void SetFitMCS(Int_t FitMCS);
+
+ AliMUONTrackParam* GetTrackParamAtFirstHit(void);
void RecursiveDump(void); // Recursive dump (with track hits)
void Fit(AliMUONTrackParam *TrackParam, Int_t NParam); // Fit
AliMUONTrackParam fTrackParamAtVertex; // Track parameters at vertex
TClonesArray *fTrackHitsPtr; // Pointer to array of TrackHit's
Int_t fNTrackHits; // Number of TrackHit's
+ Int_t fFitMCS; // 0(1) for fit without(with) multiple Coulomb scattering
ClassDef(AliMUONTrack, 1) // Class definition in ROOT context
};
/*
$Log$
+Revision 1.2 2000/06/15 07:58:49 morsch
+Code from MUON-dev joined
+
Revision 1.1.2.3 2000/06/09 21:03:09 morsch
Make includes consistent with new file structure.
return;
}
+ //__________________________________________________________________________
+void AliMUONTrackParam::ExtrapToVertex()
+{
+ // Extrapolation to the vertex.
+ // Returns the track parameters resulting from the extrapolation,
+ // in the current TrackParam.
+ // Changes parameters according to branson correction through the absorber
+
+ Double_t zAbsorber = 503.0; // to be coherent with the Geant absorber geometry !!!!
+ // Extrapolates track parameters upstream to the "Z" end of the front absorber
+ ExtrapToZ(zAbsorber);
+ // Makes Branson correction (multiple scattering + energy loss)
+ BransonCorrection();
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackParam::BransonCorrection()
+{
+ // Branson correction of track parameters
+ // the entry parameters have to be calculated at the end of the absorber
+ Double_t zEndAbsorber, zBP, xBP, yBP;
+ Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
+ Int_t sign;
+ // Would it be possible to calculate all that from Geant configuration ????
+ // Radiation lengths outer part theta > 3 degres
+ static Double_t x01[9] = { 18.8, // C (cm)
+ 10.397, // Concrete (cm)
+ 0.56, // Plomb (cm)
+ 47.26, // Polyethylene (cm)
+ 0.56, // Plomb (cm)
+ 47.26, // Polyethylene (cm)
+ 0.56, // Plomb (cm)
+ 47.26, // Polyethylene (cm)
+ 0.56 }; // Plomb (cm)
+ // inner part theta < 3 degres
+ static Double_t x02[3] = { 18.8, // C (cm)
+ 10.397, // Concrete (cm)
+ 0.35 }; // W (cm)
+ // z positions of the materials inside the absober outer part theta > 3 degres
+ static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 };
+ // inner part theta < 3 degres
+ static Double_t z2[4] = { 90, 315, 467, 503 };
+ static Bool_t first = kTRUE;
+ static Double_t zBP1, zBP2, rLimit;
+ // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call)
+ if (first) {
+ first = kFALSE;
+ Double_t aNBP = 0.0;
+ Double_t aDBP = 0.0;
+ for (Int_t iBound = 0; iBound < 9; iBound++) {
+ aNBP = aNBP +
+ (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
+ z1[iBound] * z1[iBound] * z1[iBound] ) / x01[iBound];
+ aDBP = aDBP +
+ (z1[iBound+1] * z1[iBound+1] - z1[iBound] * z1[iBound] ) / x01[iBound];
+ }
+ zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
+ aNBP = 0.0;
+ aDBP = 0.0;
+ for (Int_t iBound = 0; iBound < 3; iBound++) {
+ aNBP = aNBP +
+ (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
+ z2[iBound] * z2[iBound ] * z2[iBound] ) / x02[iBound];
+ aDBP = aDBP +
+ (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
+ }
+ zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
+ rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
+ }
+
+ pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+ sign = 1;
+ if (fInverseBendingMomentum < 0) sign = -1;
+ pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));
+ pX = pZ * fNonBendingSlope;
+ pY = pZ * fBendingSlope;
+ pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
+ xEndAbsorber = fNonBendingCoor;
+ yEndAbsorber = fBendingCoor;
+ radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
+
+ if (radiusEndAbsorber2 > rLimit*rLimit) {
+ zEndAbsorber = z1[9];
+ zBP = zBP1;
+ } else {
+ zEndAbsorber = z2[2];
+ zBP = zBP2;
+ }
+
+ xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
+ yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
+
+ // new parameters after Branson and energy loss corrections
+ pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
+ pX = pZ * xBP / zBP;
+ pY = pZ * yBP / zBP;
+ fBendingSlope = pY / pZ;
+ fNonBendingSlope = pX / pZ;
+
+ pT = TMath::Sqrt(pX * pX + pY * pY);
+ theta = TMath::ATan2(pT, pZ);
+ pTotal =
+ TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
+
+ fInverseBendingMomentum = (sign / pTotal) *
+ TMath::Sqrt(1.0 +
+ fBendingSlope * fBendingSlope +
+ fNonBendingSlope * fNonBendingSlope) /
+ TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
+
+ // vertex position at (0,0,0)
+ // should be taken from vertex measurement ???
+ fBendingCoor = 0.0;
+ fNonBendingCoor = 0;
+ fZ= 0;
+}
+
+ //__________________________________________________________________________
+Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber)
+{
+ // Returns the total momentum corrected from energy loss in the front absorber
+ Double_t deltaP, pTotalCorrected;
+
+ Double_t radiusEndAbsorber2 =
+ xEndAbsorber *xEndAbsorber + yEndAbsorber * yEndAbsorber;
+ // Parametrization to be redone according to change of absorber material ????
+ // The name is not so good, and there are many arguments !!!!
+ if (radiusEndAbsorber2 < rLimit * rLimit) {
+ if (pTotal < 15) {
+ deltaP = 2.737 + 0.0494 * pTotal - 0.001123 * pTotal * pTotal;
+ } else {
+ deltaP = 3.0643 + 0.01346 *pTotal;
+ }
+ } else {
+ if (pTotal < 15) {
+ deltaP = 2.1380 + 0.0351 * pTotal - 0.000853 * pTotal * pTotal;
+ } else {
+ deltaP = 2.407 + 0.00702 * pTotal;
+ }
+ }
+ pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
+ return pTotalCorrected;
+}
+
void ExtrapToZ(Double_t Z);
void ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam);
+ void ExtrapToVertex(); // extrapolation to vertex through the absorber
+ void BransonCorrection(); // makes Branson correction
+ Double_t TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber); // returns total momentum after energy loss correction in the absorber
protected:
private:
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff