using absorber material properties (density, radiation-length, ...).
(b) All the informations about the absorber are get from the geometry rootfile
instead of from constants in AliMUONConstant.
(c) The track extrapolation to vertex has been removed from the tracking
procedure since the vertex is not known at this level.
(d) An option has been added to cancel the track extrapolation to vertex
when filling ESD in AliMUONReconstructor. ESD are then filled with track
parameters at the first chamber. Similar option have been added in
MUONmassPlot_ESD.C and MUONefficiency.C to perform extrapolation at
this level if not done before. (a detailled mail will follow this commit).
(Philippe P.)
Float_t yRec=0;
Float_t zRec=0;
- trackParam = recTrack->GetTrackParamAtVertex();
- xRec = trackParam->GetNonBendingCoor();
- yRec = trackParam->GetBendingCoor();
- zRec = trackParam->GetZ();
- points->SetPoint(iPoint,xRec,yRec,zRec);
+ // vertex unknown at the tracking level -> put it at (0,0,0)
+ points->SetPoint(iPoint,0.,0.,0.);
iPoint++;
for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
void AliMUONDisplay::PrintTrack(Int_t iRecTracks, AliMUONTrack *recTrack)
{
/// Print reconstructed track
-
+
AliMUONTrackParam *trackParam;
Float_t vertex[3], momentum[3];
Float_t pYZ, bendingSlope, nonBendingSlope, chi2dof;
Int_t charge;
- trackParam = recTrack->GetTrackParamAtVertex();
+ trackParam = recTrack->GetTrackParamAtVertex(); // meaningless since the vertex is not known at the tracking level
vertex[0] = trackParam->GetNonBendingCoor();
vertex[1] = trackParam->GetBendingCoor();
vertex[2] = trackParam->GetZ();
chi2dof = recTrack->GetFitFMin()/(2.0 * recTrack->GetNTrackHits() - 5.);
printf("===================================================\n");
+ printf("//*****************************************************************//\n");
+ printf("// meaningless since the vertex is not known at the tracking level //\n");
+ printf("//*****************************************************************//\n");
printf(" Reconstructed track # %d \n",iRecTracks);
printf(" charge: %d \n",charge);
printf(" vertex x,y,z (cm): %f %f %f \n",vertex[0],vertex[1],vertex[2]);
}
recModel->SetGhostChi2Cut(10);
recModel->SetEventNumber(evtNumber);
-
+
TTask* calibration = GetCalibrationTask();
loader->LoadRecPoints("RECREATE");
AliDebug(1,Form("Event %d",iEvent));
runLoader->GetEvent(iEvent++);
-
+
//----------------------- raw2digits & raw2trigger-------------------
// if (!loader->TreeD())
// {
clusterTimer.Start(kFALSE);
if (!loader->TreeR()) loader->MakeRecPointsContainer();
-
+
// tracking branch
fMUONData->MakeBranch("RC");
fMUONData->SetTreeAddress("RC");
loader->WriteRecPoints("OVERWRITE");
clusterTimer.Stop();
-
+
//--------------------------- Resetting branches -----------------------
fMUONData->ResetDigits();
// setting pointer for tracks, triggertracks & trackparam at vertex
AliMUONTrack* recTrack = 0;
- AliMUONTrackParam* trackParam = 0;
+ AliMUONTrackParam trackParam;
AliMUONTriggerTrack* recTriggerTrack = 0;
iEvent = runLoader->GetEventNumber();
// reading info from tracks
recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
- trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
+ trackParam = *((AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First());
- if (esdVert->GetNContributors())
- AliMUONTrackExtrap::ExtrapToVertex(trackParam, vertex[0],vertex[1],vertex[2]);
-
- bendingSlope = trackParam->GetBendingSlope();
- nonBendingSlope = trackParam->GetNonBendingSlope();
- inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
- xRec = trackParam->GetNonBendingCoor();
- yRec = trackParam->GetBendingCoor();
- zRec = trackParam->GetZ();
+ // extrapolate to the vertex if required
+ // if the vertex is not available, extrapolate to (0,0,0)
+ if (!strstr(GetOption(),"NoExtrapToVtx")) {
+ if (esdVert->GetNContributors())
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam, vertex[0],vertex[1],vertex[2]);
+ else
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0.,0.,0.);
+ }
+
+ bendingSlope = trackParam.GetBendingSlope();
+ nonBendingSlope = trackParam.GetNonBendingSlope();
+ inverseBendingMomentum = trackParam.GetInverseBendingMomentum();
+ xRec = trackParam.GetNonBendingCoor();
+ yRec = trackParam.GetBendingCoor();
+ zRec = trackParam.GetZ();
nTrackHits = recTrack->GetNTrackHits();
fitFmin = recTrack->GetFitFMin();
static const Double_t fgkMaxTrackingDistanceBending; ///< Maximum distance to the track to search for compatible hitForRec(s) in bending direction
static const Double_t fgkMaxTrackingDistanceNonBending; ///< Maximum distance to the track to search for compatible hitForRec(s) in non bending direction
- AliMUONTrackParam fTrackParamAtVertex; ///< Track parameters at vertex
+ AliMUONTrackParam fTrackParamAtVertex; //!< Track parameters at vertex
TClonesArray *fTrackParamAtHit; ///< Track parameters at hit
TClonesArray *fHitForRecAtHit; ///< Cluster parameters at hit
Int_t fNTrackHits; ///< Number of hits attached to the track
Int_t fTrackID; ///< track ID = track number in TrackRefs
- ClassDef(AliMUONTrack, 3) // Reconstructed track in ALICE dimuon spectrometer
+ ClassDef(AliMUONTrack, 4) // Reconstructed track in ALICE dimuon spectrometer
};
#endif
//
///////////////////////////////////////////////////
-#include <Riostream.h>
-#include <TMath.h>
-#include <TMatrixD.h>
-
#include "AliMUONTrackExtrap.h"
#include "AliMUONTrackParam.h"
#include "AliMUONConstants.h"
+
#include "AliMagF.h"
-#include "AliLog.h"
-#include "AliTracker.h"
+
+#include <Riostream.h>
+#include <TMath.h>
+#include <TMatrixD.h>
+#include <TGeoManager.h>
ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context
//__________________________________________________________________________
void AliMUONTrackExtrap::ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam, Double_t zVtx)
{
- /// Extrapolation to the vertex (at the z position "zVtx") without Branson and Field correction.
+ /// Extrapolation to the vertex (at the z position "zVtx") without Branson and energy loss corrections.
/// Returns the track parameters resulting from the extrapolation in the current TrackParam.
/// Include multiple Coulomb scattering effects in trackParam covariances.
+ if (trackParam->GetZ() == zVtx) return; // nothing to be done if already at vertex
+
if (trackParam->GetZ() > zVtx) { // spectro. (z<0)
cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<<trackParam->GetZ()
<<") upstream the vertex (zVtx = "<<zVtx<<")"<<endl;
exit(-1);
}
- if (zVtx < AliMUONConstants::ZAbsorberEnd()) { // spectro. (z<0)
- cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Ending Z ("<<zVtx
- <<") downstream the front absorber (zAbsorberEnd = "<<AliMUONConstants::ZAbsorberEnd()<<")"<<endl;
-
+ // Check whether the geometry is available and get absorber boundaries
+ if (!gGeoManager) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl;
+ return;
+ }
+ TGeoNode *absNode = gGeoManager->GetVolume("ALIC")->GetNode("ABSM_1");
+ if (!absNode) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: failed to get absorber node"<<endl;
+ return;
+ }
+ Double_t zAbsBeg, zAbsEnd;
+ absNode->GetVolume()->GetShape()->GetAxisRange(3,zAbsBeg,zAbsEnd);
+ const Double_t *absPos = absNode->GetMatrix()->GetTranslation();
+ zAbsBeg = absPos[2] - zAbsBeg; // spectro. (z<0)
+ zAbsEnd = absPos[2] - zAbsEnd; // spectro. (z<0)
+
+ // Check the vertex position relatively to the absorber
+ if (zVtx < zAbsBeg && zVtx > zAbsEnd) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
+ <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
+ } else if (zVtx < zAbsEnd ) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
+ <<") downstream the front absorber (zAbsorberEnd = "<<zAbsEnd<<")"<<endl;
ExtrapToZCov(trackParam,zVtx);
return;
}
- // First Extrapolates track parameters upstream to the "Z" end of the front absorber
- if (trackParam->GetZ() < AliMUONConstants::ZAbsorberEnd()) { // spectro. (z<0)
- ExtrapToZCov(trackParam,AliMUONConstants::ZAbsorberEnd());
+ // Check the track position relatively to the absorber and extrapolate track parameters to the end of the absorber if needed
+ if (trackParam->GetZ() > zAbsBeg) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
+ <<") upstream the front absorber (zAbsorberBegin = "<<zAbsBeg<<")"<<endl;
+ ExtrapToZCov(trackParam,zVtx);
+ return;
+ } else if (trackParam->GetZ() > zAbsEnd) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
+ <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
} else {
- cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<<trackParam->GetZ()
- <<") upstream or inside the front absorber (zAbsorberEnd = "<<AliMUONConstants::ZAbsorberEnd()<<")"<<endl;
+ ExtrapToZCov(trackParam,zAbsEnd);
}
- // Then go through all the absorber layers
- Double_t tan3 = TMath::Tan(3./180.*TMath::Pi());
- Double_t r0Norm, x0, z, zElement, dZ, nonBendingCoor, bendingCoor;
- for (Int_t iElement=AliMUONConstants::NAbsorberElements()-1; iElement>=0; iElement--) {
- zElement = AliMUONConstants::ZAbsorberElement(iElement);
- z = trackParam->GetZ();
- if (z > zElement) continue; // spectro. (z<0)
- nonBendingCoor = trackParam->GetNonBendingCoor();
- bendingCoor = trackParam->GetBendingCoor();
- r0Norm = nonBendingCoor * nonBendingCoor + bendingCoor * bendingCoor;
- r0Norm = TMath::Sqrt(r0Norm) / TMath::Abs(trackParam->GetZ()) / tan3;
- if (r0Norm > 1.) x0 = AliMUONConstants::X0AbsorberOut(iElement); // outer part of the absorber
- else x0 = AliMUONConstants::X0AbsorberIn(iElement); // inner part of the absorber
-
- if (zVtx > zElement) {
- ExtrapToZCov(trackParam,zElement); // extrapolate to absorber element "iElement"
- dZ = zElement - z;
- AddMCSEffectInTrackParamCov(trackParam,dZ,x0); // include MCS effect in covariances
- } else {
- ExtrapToZCov(trackParam,zVtx); // extrapolate to zVtx
- dZ = zVtx - z;
- AddMCSEffectInTrackParamCov(trackParam,dZ,x0); // include MCS effect in covariances
- break;
- }
- }
+ // Then add MCS effect in absorber to the parameters covariances
+ AliMUONTrackParam trackParamIn(*trackParam);
+ ExtrapToZ(&trackParamIn, TMath::Min(zVtx, zAbsBeg));
+ Double_t trackXYZIn[3];
+ trackXYZIn[0] = trackParamIn.GetNonBendingCoor();
+ trackXYZIn[1] = trackParamIn.GetBendingCoor();
+ trackXYZIn[2] = trackParamIn.GetZ();
+ Double_t trackXYZOut[3];
+ trackXYZOut[0] = trackParam->GetNonBendingCoor();
+ trackXYZOut[1] = trackParam->GetBendingCoor();
+ trackXYZOut[2] = trackParam->GetZ();
+ Double_t pathLength = 0.;
+ Double_t f0 = 0.;
+ Double_t f1 = 0.;
+ Double_t f2 = 0.;
+ Double_t meanRho = 0.;
+ GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho);
+ AddMCSEffectInAbsorber(trackParam,pathLength,f0,f1,f2);
// finally go to the vertex
ExtrapToZCov(trackParam,zVtx);
}
//__________________________________________________________________________
-void AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(AliMUONTrackParam *param, Double_t dZ, Double_t x0)
+void AliMUONTrackExtrap::AddMCSEffectInAbsorber(AliMUONTrackParam* param, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2)
+{
+ /// Add to the track parameter covariances the effects of multiple Coulomb scattering
+ /// at the end of the front absorber using the absorber correction parameters
+
+ // absorber related covariance parameters
+ Double_t bendingSlope = param->GetBendingSlope();
+ Double_t nonBendingSlope = param->GetNonBendingSlope();
+ Double_t inverseBendingMomentum = param->GetInverseBendingMomentum();
+ Double_t alpha2 = 0.0136 * 0.0136 * inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) /
+ (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); // velocity = 1
+ Double_t varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2);
+ Double_t covCorrSlope = alpha2 * (pathLength * f0 - f1);
+ Double_t varSlop = alpha2 * f0;
+
+ TMatrixD* paramCov = param->GetCovariances();
+ // Non bending plane
+ (*paramCov)(0,0) += varCoor; (*paramCov)(0,1) += covCorrSlope;
+ (*paramCov)(1,0) += covCorrSlope; (*paramCov)(1,1) += varSlop;
+ // Bending plane
+ (*paramCov)(2,2) += varCoor; (*paramCov)(2,3) += covCorrSlope;
+ (*paramCov)(3,2) += covCorrSlope; (*paramCov)(3,3) += varSlop;
+
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Double_t trackXYZOut[3], Double_t &pathLength,
+ Double_t &f0, Double_t &f1, Double_t &f2, Double_t &meanRho)
+{
+ /// Parameters used to correct for Multiple Coulomb Scattering and energy loss in absorber
+ /// Calculated assuming a linear propagation between track positions trackXYZIn and trackXYZOut
+ // pathLength: path length between trackXYZIn and trackXYZOut (cm)
+ // f0: 0th moment of z calculated with the inverse radiation-length distribution
+ // f1: 1st moment of z calculated with the inverse radiation-length distribution
+ // f2: 2nd moment of z calculated with the inverse radiation-length distribution
+ // meanRho: average density of crossed material (g/cm3)
+
+ // Reset absorber's parameters
+ pathLength = 0.;
+ f0 = 0.;
+ f1 = 0.;
+ f2 = 0.;
+ meanRho = 0.;
+
+ // Check whether the geometry is available
+ if (!gGeoManager) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl;
+ return;
+ }
+
+ // Initialize starting point and direction
+ pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0])*(trackXYZOut[0] - trackXYZIn[0])+
+ (trackXYZOut[1] - trackXYZIn[1])*(trackXYZOut[1] - trackXYZIn[1])+
+ (trackXYZOut[2] - trackXYZIn[2])*(trackXYZOut[2] - trackXYZIn[2]));
+ if (pathLength < TGeoShape::Tolerance()) return;
+ Double_t b[3];
+ b[0] = (trackXYZOut[0] - trackXYZIn[0]) / pathLength;
+ b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength;
+ b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength;
+ TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b);
+ if (!currentnode) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<<endl;
+ return;
+ }
+
+ // loop over absorber slices and calculate absorber's parameters
+ Double_t rho = 0.; // material density (g/cm3)
+ Double_t x0 = 0.; // radiation-length (cm-1)
+ Double_t localPathLength = 0;
+ Double_t remainingPathLength = pathLength;
+ Double_t zB = trackXYZIn[2];
+ Double_t zE, dzB, dzE;
+ do {
+ // Get material properties
+ TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
+ rho = material->GetDensity();
+ x0 = material->GetRadLen();
+ if (!material->IsMixture()) x0 /= rho; // different normalization in the modeler for mixture
+
+ // Get path length within this material
+ gGeoManager->FindNextBoundary(remainingPathLength);
+ localPathLength = gGeoManager->GetStep() + 1.e-6;
+ // Check if boundary within remaining path length. If so, make sure to cross the boundary to prepare the next step
+ if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
+ else {
+ currentnode = gGeoManager->Step();
+ if (!currentnode) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
+ f0 = f1 = f2 = meanRho = 0.;
+ return;
+ }
+ if (!gGeoManager->IsEntering()) {
+ // make another small step to try to enter in new absorber slice
+ gGeoManager->SetStep(0.001);
+ currentnode = gGeoManager->Step();
+ if (!gGeoManager->IsEntering() || !currentnode) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
+ f0 = f1 = f2 = meanRho = 0.;
+ return;
+ }
+ localPathLength += 0.001;
+ }
+ }
+
+ // calculate absorber's parameters
+ zE = b[2] * localPathLength + zB;
+ dzB = zB - trackXYZIn[2];
+ dzE = zE - trackXYZIn[2];
+ f0 += localPathLength / x0;
+ f1 += (dzE*dzE - dzB*dzB) / b[2] / b[2] / x0 / 2.;
+ f2 += (dzE*dzE*dzE - dzB*dzB*dzB) / b[2] / b[2] / b[2] / x0 / 3.;
+ meanRho += localPathLength * rho;
+
+ // prepare next step
+ zB = zE;
+ remainingPathLength -= localPathLength;
+ } while (remainingPathLength > TGeoShape::Tolerance());
+
+ meanRho /= pathLength;
+}
+
+ //__________________________________________________________________________
+void AliMUONTrackExtrap::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Double_t x0)
{
/// Add to the track parameter covariances the effects of multiple Coulomb scattering
/// through a material of thickness "dZ" and of radiation length "x0"
{
/// 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
+ /// Changes parameters according to Branson correction through the absorber and energy loss
- // Extrapolates track parameters upstream to the "Z" end of the front absorber
- ExtrapToZ(trackParam,AliMUONConstants::ZAbsorberEnd()); // !!!
- // Makes Branson correction (multiple scattering + energy loss)
- BransonCorrection(trackParam,xVtx,yVtx,zVtx);
- // Makes a simple magnetic field correction through the absorber
- FieldCorrection(trackParam,AliMUONConstants::ZAbsorberEnd());
-}
-
-
-// Keep this version for future developments
- //__________________________________________________________________________
-// void AliMUONTrackExtrap::BransonCorrection(AliMUONTrackParam* trackParam)
-// {
-// /// 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 ????
-// // and to get the Branson parameters from a function in ABSO module ????
-// // with an eventual contribution from other detectors like START ????
-// // 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;
-// Int_t iBound;
-//
-// for (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 (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 / trackParam->GetInverseBendingMomentum());
-// sign = 1;
-// if (trackParam->GetInverseBendingMomentum() < 0) sign = -1;
-// pZ = pYZ / (TMath::Sqrt(1.0 + trackParam->GetBendingSlope() * trackParam->GetBendingSlope()));
-// pX = pZ * trackParam->GetNonBendingSlope();
-// pY = pZ * trackParam->GetBendingSlope();
-// pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
-// xEndAbsorber = trackParam->GetNonBendingCoor();
-// yEndAbsorber = trackParam->GetBendingCoor();
-// radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
-//
-// if (radiusEndAbsorber2 > rLimit*rLimit) {
-// zEndAbsorber = z1[9];
-// zBP = zBP1;
-// } else {
-// zEndAbsorber = z2[3];
-// 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;
-// trackParam->SetBendingSlope(pY/pZ);
-// trackParam->SetNonBendingSlope(pX/pZ);
-//
-// pT = TMath::Sqrt(pX * pX + pY * pY);
-// theta = TMath::ATan2(pT, pZ);
-// pTotal = TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
-//
-// trackParam->SetInverseBendingMomentum((sign / pTotal) *
-// TMath::Sqrt(1.0 +
-// trackParam->GetBendingSlope() * trackParam->GetBendingSlope() +
-// trackParam->GetNonBendingSlope() * trackParam->GetNonBendingSlope()) /
-// TMath::Sqrt(1.0 + trackParam->GetBendingSlope() * trackParam->GetBendingSlope()));
-//
-// // vertex position at (0,0,0)
-// // should be taken from vertex measurement ???
-// trackParam->SetBendingCoor(0.);
-// trackParam->SetNonBendingCoor(0.);
-// trackParam->SetZ(0.);
-// }
-
-void AliMUONTrackExtrap::BransonCorrection(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx)
-{
- /// Branson correction of track parameters
- // the entry parameters have to be calculated at the end of the absorber
- // simplified version: the z positions of Branson's planes are no longer calculated
- // but are given as inputs. One can use the macros MUONTestAbso.C and DrawTestAbso.C
- // to test this correction.
- // Would it be possible to calculate all that from Geant configuration ????
- // and to get the Branson parameters from a function in ABSO module ????
- // with an eventual contribution from other detectors like START ????
- // change to take into account the vertex postition (real, reconstruct,....)
-
- Double_t zBP, xBP, yBP;
- Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
- Int_t sign;
- static Bool_t first = kTRUE;
- static Double_t zBP1, zBP2, rLimit, thetaLimit;
- // zBP1 for outer part and zBP2 for inner part (only at the first call)
- if (first) {
- first = kFALSE;
+ if (trackParam->GetZ() == zVtx) return; // nothing to be done if already at vertex
- thetaLimit = 3.0 * (TMath::Pi()) / 180.;
- rLimit = TMath::Abs(AliMUONConstants::ZAbsorberEnd()) * TMath::Tan(thetaLimit);
- zBP1 = -450; // values close to those calculated with EvalAbso.C
- zBP2 = -480;
+ if (trackParam->GetZ() > zVtx) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
+ <<") upstream the vertex (zVtx = "<<zVtx<<")"<<endl;
+ exit(-1);
}
-
- pYZ = TMath::Abs(1.0 / trackParam->GetInverseBendingMomentum());
- sign = 1;
- if (trackParam->GetInverseBendingMomentum() < 0) sign = -1;
- pZ = trackParam->Pz();
- pX = trackParam->Px();
- pY = trackParam->Py();
- pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
- xEndAbsorber = trackParam->GetNonBendingCoor();
- yEndAbsorber = trackParam->GetBendingCoor();
- radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
-
- if (radiusEndAbsorber2 > rLimit*rLimit) {
- zBP = zBP1;
+
+ // Check whether the geometry is available and get absorber boundaries
+ if (!gGeoManager) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl;
+ return;
+ }
+ TGeoNode *absNode = gGeoManager->GetVolume("ALIC")->GetNode("ABSM_1");
+ if (!absNode) {
+ cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: failed to get absorber node"<<endl;
+ return;
+ }
+ Double_t zAbsBeg, zAbsEnd;
+ absNode->GetVolume()->GetShape()->GetAxisRange(3,zAbsBeg,zAbsEnd);
+ const Double_t *absPos = absNode->GetMatrix()->GetTranslation();
+ zAbsBeg = absPos[2] - zAbsBeg; // spectro. (z<0)
+ zAbsEnd = absPos[2] - zAbsEnd; // spectro. (z<0)
+
+ // Check the vertex position relatively to the absorber
+ if (zVtx < zAbsBeg && zVtx > zAbsEnd) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
+ <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
+ } else if (zVtx < zAbsEnd ) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
+ <<") downstream the front absorber (zAbsorberEnd = "<<zAbsEnd<<")"<<endl;
+ ExtrapToZ(trackParam,zVtx);
+ return;
+ }
+
+ // Check the track position relatively to the absorber and extrapolate track parameters to the end of the absorber if needed
+ if (trackParam->GetZ() > zAbsBeg) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
+ <<") upstream the front absorber (zAbsorberBegin = "<<zAbsBeg<<")"<<endl;
+ ExtrapToZ(trackParam,zVtx);
+ return;
+ } else if (trackParam->GetZ() > zAbsEnd) { // spectro. (z<0)
+ cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
+ <<") inside the front absorber ("<<zAbsBeg<<","<<zAbsEnd<<")"<<endl;
} else {
- zBP = zBP2;
+ ExtrapToZ(trackParam,zAbsEnd);
}
-
- xBP = xEndAbsorber - (pX / pZ) * (AliMUONConstants::ZAbsorberEnd() - zBP);
- yBP = yEndAbsorber - (pY / pZ) * (AliMUONConstants::ZAbsorberEnd() - zBP);
-
- // new parameters after Branson and energy loss corrections
-// Float_t zSmear = zBP - gRandom->Gaus(0.,2.); // !!! possible smearing of Z vertex position
-
- Float_t zSmear = zBP ;
- pZ = pTotal * (zSmear-zVtx) / TMath::Sqrt((xBP-xVtx) * (xBP-xVtx) + (yBP-yVtx) * (yBP-yVtx) +( zSmear-zVtx) * (zSmear-zVtx) );
- pX = pZ * (xBP - xVtx)/ (zSmear-zVtx);
- pY = pZ * (yBP - yVtx) / (zSmear-zVtx);
- trackParam->SetBendingSlope(pY/pZ);
- trackParam->SetNonBendingSlope(pX/pZ);
-
+ // Get absorber correction parameters assuming linear propagation from vertex to the track position
+ Double_t trackXYZOut[3];
+ trackXYZOut[0] = trackParam->GetNonBendingCoor();
+ trackXYZOut[1] = trackParam->GetBendingCoor();
+ trackXYZOut[2] = trackParam->GetZ();
+ Double_t trackXYZIn[3];
+ trackXYZIn[2] = TMath::Min(zVtx, zAbsBeg); // spectro. (z<0)
+ trackXYZIn[0] = trackXYZOut[0] + (xVtx - trackXYZOut[0]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
+ trackXYZIn[1] = trackXYZOut[1] + (yVtx - trackXYZOut[1]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
+ Double_t pathLength = 0.;
+ Double_t f0 = 0.;
+ Double_t f1 = 0.;
+ Double_t f2 = 0.;
+ Double_t meanRho = 0.;
+ GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pathLength,f0,f1,f2,meanRho);
- pT = TMath::Sqrt(pX * pX + pY * pY);
- theta = TMath::ATan2(pT, TMath::Abs(pZ));
- pTotal = TotalMomentumEnergyLoss(thetaLimit, pTotal, theta);
-
- trackParam->SetInverseBendingMomentum((sign / pTotal) *
- TMath::Sqrt(1.0 +
- trackParam->GetBendingSlope() * trackParam->GetBendingSlope() +
- trackParam->GetNonBendingSlope() * trackParam->GetNonBendingSlope()) /
- TMath::Sqrt(1.0 + trackParam->GetBendingSlope() * trackParam->GetBendingSlope()));
-
- // vertex position at (0,0,0)
- // should be taken from vertex measurement ???
-
- trackParam->SetBendingCoor(xVtx);
- trackParam->SetNonBendingCoor(yVtx);
+ // Calculate energy loss
+ Double_t pTot = trackParam->P();
+ Double_t charge = TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum());
+ Double_t deltaP = TotalMomentumEnergyLoss(pTot,pathLength,meanRho);
+
+ // Correct for half of energy loss
+ pTot += 0.5 * deltaP;
+
+ // Position of the Branson plane (spectro. (z<0))
+ Double_t zB = (f1>0.) ? trackXYZIn[2] - f2/f1 : 0.;
+
+ // Get track position in the Branson plane corrected for magnetic field effect
+ ExtrapToZ(trackParam,zVtx);
+ Double_t xB = trackParam->GetNonBendingCoor() + (zB - zVtx) * trackParam->GetNonBendingSlope();
+ Double_t yB = trackParam->GetBendingCoor() + (zB - zVtx) * trackParam->GetBendingSlope();
+
+ // Get track slopes corrected for multiple scattering (spectro. (z<0))
+ Double_t nonBendingSlope = (zB<0.) ? (xB - xVtx) / (zB - zVtx) : trackParam->GetNonBendingSlope();
+ Double_t bendingSlope = (zB<0.) ? (yB - yVtx) / (zB - zVtx) : trackParam->GetBendingSlope();
+
+ // Correct for second half of energy loss
+ pTot += 0.5 * deltaP;
+
+ // Set track parameters at vertex
+ trackParam->SetNonBendingCoor(xVtx);
+ trackParam->SetBendingCoor(yVtx);
trackParam->SetZ(zVtx);
-
+ trackParam->SetNonBendingSlope(nonBendingSlope);
+ trackParam->SetBendingSlope(bendingSlope);
+ trackParam->SetInverseBendingMomentum(charge / pTot *
+ TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) /
+ TMath::Sqrt(1.0 + bendingSlope*bendingSlope));
+
}
//__________________________________________________________________________
-Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta)
+Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(Double_t pTotal, Double_t pathLength, Double_t rho)
{
- /// Returns the total momentum corrected from energy loss in the front absorber
- // One can use the macros MUONTestAbso.C and DrawTestAbso.C
- // to test this correction.
- // Momentum energy loss behaviour evaluated with the simulation of single muons (april 2002)
- Double_t deltaP, pTotalCorrected;
-
- // Parametrization to be redone according to change of absorber material ????
- // See remark in function BransonCorrection !!!!
- // The name is not so good, and there are many arguments !!!!
- if (theta < thetaLimit ) {
- if (pTotal < 20) {
- deltaP = 2.5938 + 0.0570 * pTotal - 0.001151 * pTotal * pTotal;
- } else {
- deltaP = 3.0714 + 0.011767 *pTotal;
- }
- deltaP *= 0.75; // AZ
- } else {
- if (pTotal < 20) {
- deltaP = 2.1207 + 0.05478 * pTotal - 0.00145079 * pTotal * pTotal;
- } else {
- deltaP = 2.6069 + 0.0051705 * pTotal;
- }
- deltaP *= 0.9; // AZ
- }
- pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
- return pTotalCorrected;
+ /// Returns the total momentum energy loss in the front absorber
+ Double_t muMass = 0.10566;
+ Double_t p2=pTotal*pTotal;
+ Double_t beta2=p2/(p2 + muMass*muMass);
+ Double_t dE=ApproximateBetheBloch(beta2)*pathLength*rho;
+
+ return dE;
}
//__________________________________________________________________________
-void AliMUONTrackExtrap::FieldCorrection(AliMUONTrackParam *trackParam, Double_t zEnd)
-{
- /// Correction of the effect of the magnetic field in the absorber
- // Assume a constant field along Z axis.
- Float_t b[3],x[3];
- Double_t bZ;
- Double_t pYZ,pX,pY,pZ,pT;
- Double_t pXNew,pYNew;
- Double_t c;
-
- pYZ = TMath::Abs(1.0 / trackParam->GetInverseBendingMomentum());
- c = TMath::Sign(1.0,trackParam->GetInverseBendingMomentum()); // particle charge
-
- pZ = trackParam->Pz();
- pX = trackParam->Px();
- pY = trackParam->Py();
- pT = TMath::Sqrt(pX*pX+pY*pY);
-
- if (TMath::Abs(pZ) <= 0) return;
- x[2] = zEnd/2;
- x[0] = x[2]*trackParam->GetNonBendingSlope();
- x[1] = x[2]*trackParam->GetBendingSlope();
-
- // Take magn. field value at position x.
- if (fgkField) fgkField->Field(x,b);
- else {
- cout<<"F-AliMUONTrackExtrap::FieldCorrection: fgkField = 0x0"<<endl;
- exit(-1);
- }
- bZ = b[2];
-
- // Transverse momentum rotation
- // Parameterized with the study of DeltaPhi = phiReco - phiGen as a function of pZ.
- Double_t phiShift = c*0.436*0.0003*bZ*zEnd/pZ;
- // Rotate momentum around Z axis.
- pXNew = pX*TMath::Cos(phiShift) - pY*TMath::Sin(phiShift);
- pYNew = pX*TMath::Sin(phiShift) + pY*TMath::Cos(phiShift);
-
- trackParam->SetBendingSlope(pYNew/pZ);
- trackParam->SetNonBendingSlope(pXNew/pZ);
-
- trackParam->SetInverseBendingMomentum(c/TMath::Sqrt(pYNew*pYNew+pZ*pZ));
-
+Double_t AliMUONTrackExtrap::ApproximateBetheBloch(Double_t beta2) {
+ //------------------------------------------------------------------
+ // This is an approximation of the Bethe-Bloch formula with
+ // the density effect taken into account at beta*gamma > 3.5
+ // (the approximation is reasonable only for solid materials)
+ //------------------------------------------------------------------
+ if (beta2/(1-beta2)>3.5*3.5)
+ return 0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2);
+
+ return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2);
}
//__________________________________________________________________________
return;
}
+
//__________________________________________________________________________
void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout)
{
return;
}
+
//___________________________________________________________
void AliMUONTrackExtrap::GetField(Double_t *Position, Double_t *Field)
{
return;
}
-
static void ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam, Double_t zVtx);
static void ExtrapToVertex(AliMUONTrackParam *trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx);
- static void AddMCSEffectInTrackParamCov(AliMUONTrackParam *param, Double_t dZ, Double_t x0);
+ static void AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Double_t x0);
static void ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout);
static void ConvertTrackParamForExtrap(AliMUONTrackParam* trackParam, Double_t forwardBackward, Double_t *v3);
static void RecoverTrackParam(Double_t *v3, Double_t Charge, AliMUONTrackParam* trackParam);
- static void BransonCorrection(AliMUONTrackParam *trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx);
- static Double_t TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta);
- static void FieldCorrection(AliMUONTrackParam *trackParam, Double_t Z);
+ static void AddMCSEffectInAbsorber(AliMUONTrackParam* trackParam, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2);
+ static void GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Double_t trackXYZOut[3], Double_t &pathLength,
+ Double_t &f0, Double_t &f1, Double_t &f2, Double_t &meanRho);
+
+ static Double_t TotalMomentumEnergyLoss(Double_t pTotal, Double_t pathLength, Double_t rho);
+ static Double_t ApproximateBetheBloch(Double_t beta2);
static void ExtrapOneStepHelix(Double_t charge, Double_t step, Double_t *vect, Double_t *vout);
static void ExtrapOneStepHelix3(Double_t field, Double_t step, Double_t *vect, Double_t *vout);
// Author: Alexander Zinchenko, JINR Dubna
#include "AliMUONTrackK.h"
-#include "AliMUON.h"
+#include "AliMUONData.h"
#include "AliMUONConstants.h"
#include "AliMUONTrackReconstructorK.h"
#include "AliMUONDetElement.h"
#include "AliLog.h"
-#include "AliMagF.h"
-#include "AliRunLoader.h"
#include <Riostream.h>
#include <TClonesArray.h>
/// Compute track parameters at hit positions (as for AliMUONTrack)
// 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->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);
-
- AliMUONTrackExtrap::ExtrapToVertex(&trackParam, Double_t(0.), Double_t(0.), Double_t(0.));
-
- (*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; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- printf ("%5d", hit->GetChamberNumber());
- }
- cout << endl;
- if (fgDebug > 0) {
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
- //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
- printf ("%5d", fgHitForRec->IndexOf(hit));
- }
- cout << endl;
- }
-
- // from raw clusters
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- if (hit->GetHitNumber() < 0) { // combined cluster / track finder
- Int_t index = -hit->GetHitNumber() / 100000;
- Int_t iPos = -hit->GetHitNumber() - index * 100000;
- clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
- } else {
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- }
- printf ("%5d", clus->GetTrack(1)%10000000);
- }
- cout << endl;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
- if (hit->GetHitNumber() < 0) { // combined cluster / track finder
- Int_t index = -hit->GetHitNumber() / 100000;
- Int_t iPos = -hit->GetHitNumber() - index * 100000;
- clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
- } else {
- rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
- clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
- }
- if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000);
- else printf ("%5s", " ");
- }
- cout << endl;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
- //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
- cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
- //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " ";
- }
- cout << endl;
- for (Int_t i1=0; i1<fNmbTrackHits; i1++) printf("%8.4f", (*fChi2Smooth)[i1]);
- cout << endl;
- cout << "---------------------------------------------------" << endl;
-
- // Get covariance matrix
- /* Not needed - covariance matrix is not interesting to anybody
- *fCovariance = *fWeight;
- if (fCovariance->Determinant() != 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)
{
class TObjArray;
#include <TMatrixD.h>
-#include <TObject.h>
class AliMUONTrackK : public AliMUONTrack {
void SetTrackQuality(Int_t iChi2); // compute track quality or Chi2
Bool_t KeepTrack(AliMUONTrackK* track0) const; // keep or discard track
void Kill(void); // kill track candidate
- void Branson(void); // Branson correction
- void GoToZ(Double_t zEnd); // propagate track to given Z
- void GoToVertex(Int_t iflag); // propagate track to the vertex
Bool_t Smooth(void); // apply smoother
Double_t GetChi2PerPoint(Int_t iPoint) const; // return Chi2 at point
void Print(FILE *lun) const; // print track information
hitForRec->SetHitNumber(iclus);
// Z coordinate of the raw cluster (cm)
hitForRec->SetZ(clus->GetZ(0));
- StdoutToAliDebug(3,
- cout << "Chamber " << ch <<
- " raw cluster " << iclus << " : " << endl;
- clus->Print("full");
- cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
- hitForRec->Print("full");
- );
+ if (AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 3) {
+ cout << "Chamber " << ch <<" raw cluster " << iclus << " : " << endl;
+ clus->Print("full");
+ cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
+ hitForRec->Print("full");
+ }
} // end of cluster loop
} // end of chamber loop
SortHitsForRecWithIncreasingChamber();
FollowTracks();
// Remove double tracks
RemoveDoubleTracks();
- // Propagate tracks to the vertex through absorber
- ExtrapTracksToVertex();
// Fill out the AliMUONTrack's
FillMUONTrack();
}
fNRecTracks++;
// Add MCS effects in parameter covariances
trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
- AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
+ AliMUONTrackExtrap::AddMCSEffect(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
// Printout for debuging
if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
}
// Add MCS effects in parameter covariances
trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
- AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
+ AliMUONTrackExtrap::AddMCSEffect(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
// Printout for debuging
if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
Int_t npari, nparx;
Int_t status, covStatus;
+ // Instantiate gMinuit if not already done
+ if (!gMinuit) gMinuit = new TMinuit(6);
// Clear MINUIT parameters
gMinuit->mncler();
// Give the fitted track to MINUIT
// Calculate the covariance matrix more accurately if required
if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
-
+
// get results into "invBenP", "benC", "nonBenC" ("x", "y")
gMinuit->GetParameter(0, x, errorParam);
trackParam->SetNonBendingCoor(x);
// The velocity is assumed to be 1 !!!!
varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * varMultipleScatteringAngle * varMultipleScatteringAngle;
return varMultipleScatteringAngle;
-}
-
- //__________________________________________________________________________
-void AliMUONTrackReconstructor::ExtrapTracksToVertex()
-{
- /// propagate tracks to the vertex through the absorber (Branson)
- AliMUONTrack *track;
- AliMUONTrackParam trackParamVertex;
- AliMUONHitForRec *vertex;
- Double_t vX, vY, vZ;
-
- for (Int_t iRecTrack = 0; iRecTrack <fNRecTracks; iRecTrack++) {
- track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
- trackParamVertex = *((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First()));
- vertex = track->GetVertex();
- if (vertex) { // if track as a vertex defined, use it
- vX = vertex->GetNonBendingCoor();
- vY = vertex->GetBendingCoor();
- vZ = vertex->GetZ();
- } else { // else vertex = (0.,0.,0.)
- vX = 0.;
- vY = 0.;
- vZ = 0.;
- }
- AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, vX, vY, vZ);
- track->SetTrackParamAtVertex(&trackParamVertex);
-
- if (AliLog::GetGlobalDebugLevel() > 0) {
- cout << "FollowTracks: track candidate(0..): " << iRecTrack
- << " after extrapolation to vertex" << endl;
- track->RecursiveDump();
- }
- }
-
}
//__________________________________________________________________________
/// Dump reconstructed event (track parameters at vertex and at first hit),
/// and the particle parameters
AliMUONTrack *track;
- AliMUONTrackParam *trackParam, *trackParam1;
+ AliMUONTrackParam trackParam, *trackParam1;
Double_t bendingSlope, nonBendingSlope, pYZ;
Double_t pX, pY, pZ, x, y, z, c;
Int_t trackIndex, nTrackHits;
nTrackHits = track->GetNTrackHits();
AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
// 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());
+ trackParam = (*((AliMUONTrackParam*) track->GetTrackParamAtHit()->First()));
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam,0.,0.,0.);
+ 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());
+ c = TMath::Sign(1.0, trackParam.GetInverseBendingMomentum());
AliDebug(1, Form("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));
virtual void MakeTrackCandidates(void);
virtual void FollowTracks(void);
virtual void RemoveDoubleTracks(void);
- virtual void ExtrapTracksToVertex(void);
virtual void FillMUONTrack(void);
///
////////////////////////////////////
-#include <stdlib.h>
-#include <Riostream.h>
-#include <TDirectory.h>
-#include <TFile.h>
-#include <TMatrixD.h>
-
-#include "AliMUONVTrackReconstructor.h"
#include "AliMUONTrackReconstructorK.h"
#include "AliMUONData.h"
#include "AliMUONConstants.h"
#include "AliMUONObjectPair.h"
#include "AliMUONRawCluster.h"
#include "AliMUONTrackK.h"
+
#include "AliLog.h"
+#include <Riostream.h>
+
+/// \cond CLASSIMP
ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
+ClassImp(AliMUONConstants)
+/// \endcond
//__________________________________________________________________________
AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONData* data, const Option_t* TrackMethod)
hitForRec->SetHitNumber(iclus);
// Z coordinate of the raw cluster (cm)
hitForRec->SetZ(clus->GetZ(0));
- StdoutToAliDebug(3,
- cout << "Chamber " << ch <<
- " raw cluster " << iclus << " : " << endl;
- clus->Print("full");
- cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
- hitForRec->Print("full");
- );
+ if (AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 3) {
+ cout << "Chamber " << ch <<" raw cluster " << iclus << " : " << endl;
+ clus->Print("full");
+ cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
+ hitForRec->Print("full");
+ }
} // end of cluster loop
} // end of chamber loop
SortHitsForRecWithIncreasingChamber();
FollowTracks();
// Remove double tracks
RemoveDoubleTracks();
- // Propagate tracks to the vertex through absorber
- ExtrapTracksToVertex();
// Fill AliMUONTrack data members
FillMUONTrack();
}
if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
}
- //__________________________________________________________________________
-void AliMUONTrackReconstructorK::ExtrapTracksToVertex(void)
-{
- /// Propagates track to the vertex thru absorber
- /// (using Branson correction for now)
- Double_t zVertex;
- zVertex = 0;
- for (Int_t i=0; i<fNRecTracks; i++) {
- //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
- ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
- //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
- ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber
- }
-}
-
//__________________________________________________________________________
void AliMUONTrackReconstructorK::FillMUONTrack()
{
/// MUON track reconstructor using kalman filter
////////////////////////////////////////////////
-#include <TObject.h>
#include "AliMUONVTrackReconstructor.h"
class AliMUONTrackReconstructorK : public AliMUONVTrackReconstructor {
virtual void MakeTrackCandidates(void);
virtual void FollowTracks(void);
virtual void RemoveDoubleTracks(void);
- virtual void ExtrapTracksToVertex(void);
virtual void FillMUONTrack(void);
virtual void MakeTrackCandidates(void) = 0;
virtual void FollowTracks(void) = 0;
virtual void RemoveDoubleTracks(void) = 0;
- virtual void ExtrapTracksToVertex(void) = 0;
virtual void FillMUONTrack(void) = 0;
private:
Int_t iTrackOk=0;
AliMUONTrack* track = amdi.RecTrack(iTrack);
while(track) {
- Double_t invBenMom = track->GetTrackParamAtVertex()->GetInverseBendingMomentum();
+ AliMUONTrackParam trackParam(*((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First())));
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam,0.,0.,0.);
+ Double_t invBenMom = trackParam.GetInverseBendingMomentum();
fInvBenMom->Fill(invBenMom);
fBenMom->Fill(1./invBenMom);
if (TMath::Abs(invBenMom)<=1.04) {
for(Int_t rectracki=0; rectracki < nrectracks;rectracki++)
{
rectrack = amdi.RecTrack(rectracki);
- trackparam = rectrack->GetTrackParamAtVertex();
+ trackparam = rectrack->GetTrackParamAtVertex(); // meaningless since the vertex is not known at the tracking level
x = trackparam->GetNonBendingCoor();
y = trackparam->GetBendingCoor();
z = trackparam->GetZ();
#include "AliMUONTrack.h"
#include "AliMUONRecoCheck.h"
#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
#include "AliTracker.h"
Int_t TrackCheck( Bool_t *compTrack);
p1 = trackParam->P();
// printf(" Ref. track at vertex: x,y,z: %f %f %f px,py,pz,p: %f %f %f %f \n",x1,y1,z1,pX1,pY1,pZ1,p1);
-
- trackParam = ((AliMUONTrack *)trackRecoArray->At(indexOK))->GetTrackParamAtVertex();
+ trackReco = (AliMUONTrack *)trackRecoArray->At(indexOK);
+ trackParam = new AliMUONTrackParam(*((AliMUONTrackParam*)(trackReco->GetTrackParamAtHit()->First())));
+ AliMUONTrackExtrap::ExtrapToVertex(trackParam,x1,y1,z1);
x2 = trackParam->GetNonBendingCoor();
y2 = trackParam->GetBendingCoor();
z2 = trackParam->GetZ();
pY2 = trackParam->Py();
pZ2 = trackParam->Pz();
p2 = trackParam->P();
+ delete trackParam;
// printf(" Reconst. track at vertex: x,y,z: %f %f %f px,py,pz: %f %f %f %f \n",x2,y2,z2,pX2,pY2,pZ2,p2);
hResMomVertex->Fill(p2-p1);
#include "TTree.h"
#include "TString.h"
#include <Riostream.h>
+#include <TGeoManager.h>
// STEER includes
#include "AliRun.h"
#include "AliESDMuonTrack.h"
#endif
+// Arguments:
+// ExtrapToVertex (default -1)
+// <0: no extrapolation;
+// =0: extrapolation to (0,0,0);
+// >0: extrapolation to ESDVertex if available, else to (0,0,0)
+// ResType (default 553)
+// 553 for Upsilon, anything else for J/Psi
-Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEvent = 1000000,
- char* esdFileName = "AliESDs.root", char* filename = "galice.root")
+Bool_t MUONefficiency( Int_t ExtrapToVertex = -1, Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEvent = 1000000,
+ char* geoFilename = "geometry.root", char* esdFileName = "AliESDs.root", char* filename = "galice.root")
{ // MUONefficiency starts
Double_t MUON_MASS = 0.105658369;
TLorentzVector fV1, fV2, fVtot;
+ // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
+ if (!gGeoManager) {
+ TGeoManager::Import(geoFilename);
+ if (!gGeoManager) {
+ Error("MUONmass_ESD", "getting geometry from file %s failed", filename);
+ return kFALSE;
+ }
+ }
+
// set mag field
// waiting for mag field in CDB
printf("Loading field map...\n");
// loop over all reconstructed tracks (also first track of combination)
for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
- AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
+ AliESDMuonTrack* muonTrack = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack)));
- if (!Vertex->GetNContributors()) {
- //re-extrapolate to vertex, if not kown before.
- trackParam.GetParamFrom(*muonTrack);
+ // extrapolate to vertex if required and available
+ if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
+ trackParam.GetParamFrom(*muonTrack);
AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
- trackParam.SetParamFor(*muonTrack);
+ trackParam.SetParamFor(*muonTrack); // put the new parameters in this copy of AliESDMuonTrack
+ } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
+ trackParam.GetParamFrom(*muonTrack);
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0.);
+ trackParam.SetParamFor(*muonTrack); // put the new parameters in this copy of AliESDMuonTrack
}
// Trigger
// loop over second track of combination
for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) {
- AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack2);
+ AliESDMuonTrack* muonTrack2 = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack2)));
- if (!Vertex->GetNContributors()) {
- trackParam.GetParamFrom(*muonTrack);
+ // extrapolate to vertex if required and available
+ if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
+ trackParam.GetParamFrom(*muonTrack2);
AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
- trackParam.SetParamFor(*muonTrack);
+ trackParam.SetParamFor(*muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
+ } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
+ trackParam.GetParamFrom(*muonTrack2);
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0.);
+ trackParam.SetParamFor(*muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
}
- track2Trigger = muonTrack->GetMatchTrigger();
+ track2Trigger = muonTrack2->GetMatchTrigger();
if (track2Trigger)
- track2TriggerChi2 = muonTrack->GetChi2MatchTrigger();
+ track2TriggerChi2 = muonTrack2->GetChi2MatchTrigger();
else
track2TriggerChi2 = 0. ;
- thetaX = muonTrack->GetThetaX();
- thetaY = muonTrack->GetThetaY();
+ thetaX = muonTrack2->GetThetaX();
+ thetaY = muonTrack2->GetThetaY();
- pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
+ pYZ = 1./TMath::Abs(muonTrack2->GetInverseBendingMomentum());
fPzRec2 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY));
fPxRec2 = fPzRec2 * TMath::Tan(thetaX);
fPyRec2 = fPzRec2 * TMath::Tan(thetaY);
- fCharge2 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
+ fCharge2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum()));
fE2 = TMath::Sqrt(MUON_MASS * MUON_MASS + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
- ntrackhits = muonTrack->GetNHit();
- fitfmin = muonTrack->GetChi2();
+ ntrackhits = muonTrack2->GetNHit();
+ fitfmin = muonTrack2->GetChi2();
// transverse momentum
Float_t pt2 = fV2.Pt();
hPtResonance->Fill(fVtot.Pt());
// match with trigger
- if (muonTrack->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig)) EventInMassMatch++;
+ if (muonTrack2->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig)) EventInMassMatch++;
}
} // if (fCharge1 * fCharge2) == -1)
} // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
+ delete muonTrack2;
} // for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++)
} // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
+ delete muonTrack;
} // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++)
hNumberOfTrack->Fill(nTracks);
#include "TParticle.h"
#include "TTree.h"
#include <Riostream.h>
+#include <TGeoManager.h>
// STEER includes
#include "AliRun.h"
// using Invariant Mass for rapidity.
// Arguments:
+// ExtrapToVertex (default -1)
+// <0: no extrapolation;
+// =0: extrapolation to (0,0,0);
+// >0: extrapolation to ESDVertex if available, else to (0,0,0)
// FirstEvent (default 0)
// LastEvent (default 0)
// ResType (default 553)
// Add parameters and histograms for analysis
-Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t LastEvent = 10000,
- char* esdFileName = "AliESDs.root", Int_t ResType = 553,
+Bool_t MUONmassPlot(Int_t ExtrapToVertex = -1, char* geoFilename = "geometry.root", char* filename = "galice.root",
+ Int_t FirstEvent = 0, Int_t LastEvent = 10000, char* esdFileName = "AliESDs.root", Int_t ResType = 553,
Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,
Float_t massMin = 9.17,Float_t massMax = 9.77)
{
TLorentzVector fV1, fV2, fVtot;
+ // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
+ if (!gGeoManager) {
+ TGeoManager::Import(geoFilename);
+ if (!gGeoManager) {
+ Error("MUONmass_ESD", "getting geometry from file %s failed", filename);
+ return kFALSE;
+ }
+ }
+
// set mag field
// waiting for mag field in CDB
printf("Loading field map...\n");
// open run loader and load gAlice, kinematics and header
AliRunLoader* runLoader = AliRunLoader::Open(filename);
if (!runLoader) {
- Error("MUONmass_ESD", "getting run loader from file %s failed",
- filename);
+ Error("MUONmass_ESD", "getting run loader from file %s failed", filename);
return kFALSE;
}
// get the SPD reconstructed vertex (vertexer) and fill the histogram
AliESDVertex* Vertex = (AliESDVertex*) esd->GetVertex();
-
if (Vertex->GetNContributors()) {
fZVertex = Vertex->GetZv();
fYVertex = Vertex->GetYv();
fXVertex = Vertex->GetXv();
-
}
hPrimaryVertex->Fill(fZVertex);
// loop over all reconstructed tracks (also first track of combination)
for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
- AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
+ AliESDMuonTrack* muonTrack = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack)));
- if (!Vertex->GetNContributors()) {
- //re-extrapolate to vertex, if not kown before.
- trackParam.GetParamFrom(*muonTrack);
+ // extrapolate to vertex if required and available
+ if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
+ trackParam.GetParamFrom(*muonTrack);
AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
- trackParam.SetParamFor(*muonTrack);
+ trackParam.SetParamFor(*muonTrack); // put the new parameters in this copy of AliESDMuonTrack
+ } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
+ trackParam.GetParamFrom(*muonTrack);
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0.);
+ trackParam.SetParamFor(*muonTrack); // put the new parameters in this copy of AliESDMuonTrack
}
+
thetaX = muonTrack->GetThetaX();
thetaY = muonTrack->GetThetaY();
-
+
pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
fPzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY));
fPxRec1 = fPzRec1 * TMath::Tan(thetaX);
// loop over second track of combination
for (Int_t iTrack2 = iTrack + 1; iTrack2 < nTracks; iTrack2++) {
- AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack2);
-
- if (!Vertex->GetNContributors()) {
- trackParam.GetParamFrom(*muonTrack);
+ AliESDMuonTrack* muonTrack2 = new AliESDMuonTrack(*(esd->GetMuonTrack(iTrack2)));
+
+ // extrapolate to vertex if required and available
+ if (ExtrapToVertex > 0 && Vertex->GetNContributors()) {
+ trackParam.GetParamFrom(*muonTrack2);
AliMUONTrackExtrap::ExtrapToVertex(&trackParam, fXVertex, fYVertex, fZVertex);
- trackParam.SetParamFor(*muonTrack);
+ trackParam.SetParamFor(*muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
+ } else if ((ExtrapToVertex > 0 && !Vertex->GetNContributors()) || ExtrapToVertex == 0){
+ trackParam.GetParamFrom(*muonTrack2);
+ AliMUONTrackExtrap::ExtrapToVertex(&trackParam, 0., 0., 0.);
+ trackParam.SetParamFor(*muonTrack2); // put the new parameters in this copy of AliESDMuonTrack
}
+
+ thetaX = muonTrack2->GetThetaX();
+ thetaY = muonTrack2->GetThetaY();
- thetaX = muonTrack->GetThetaX();
- thetaY = muonTrack->GetThetaY();
-
- pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
+ pYZ = 1./TMath::Abs(muonTrack2->GetInverseBendingMomentum());
fPzRec2 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaY));
fPxRec2 = fPzRec2 * TMath::Tan(thetaX);
fPyRec2 = fPzRec2 * TMath::Tan(thetaY);
- fCharge2 = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
+ fCharge2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum()));
fE2 = TMath::Sqrt(muonMass * muonMass + fPxRec2 * fPxRec2 + fPyRec2 * fPyRec2 + fPzRec2 * fPzRec2);
fV2.SetPxPyPzE(fPxRec2, fPyRec2, fPzRec2, fE2);
- ntrackhits = muonTrack->GetNHit();
- fitfmin = muonTrack->GetChi2();
+ ntrackhits = muonTrack2->GetNHit();
+ fitfmin = muonTrack2->GetChi2();
// transverse momentum
Float_t pt2 = fV2.Pt();
if (esd->GetTriggerMask() & ptTrig) NbTrigger++;
if (invMass > massMin && invMass < massMax) {
EventInMass++;
- if (muonTrack->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig))// match with trigger
+ if (muonTrack2->GetMatchTrigger() && (esd->GetTriggerMask() & ptTrig))// match with trigger
EventInMassMatch++;
hRapResonance->Fill(fVtot.Rapidity());
} // if (fCharge * fCharge2) == -1)
} // if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
+ delete muonTrack2;
} // for (Int_t iTrack2 = iTrack + 1; iTrack2 < iTrack; iTrack2++)
} // if (ch1 < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
+ delete muonTrack;
} // for (Int_t iTrack = 0; iTrack < nrectracks; iTrack++)
hNumberOfTrack->Fill(nTracks);
#include "AliMUONLocalTrigger.h"
#include "AliMUONTrack.h"
#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
#include "AliESDMuonTrack.h"
#endif
//
rectrack = (AliMUONTrack*) recTracksArray->At(irectracks);
- trackParam = rectrack->GetTrackParamAtVertex();
+ trackParam = (AliMUONTrackParam*)(rectrack->GetTrackParamAtHit()->First());
+ AliMUONTrackExtrap::ExtrapToVertex(trackParam,0.,0.,0.);
bendingSlope = trackParam->GetBendingSlope();
nonBendingSlope = trackParam->GetNonBendingSlope();
rectrack = (AliMUONTrack*) recTracksArray->At(irectracks2);
- trackParam = rectrack->GetTrackParamAtVertex();
+ trackParam = (AliMUONTrackParam*)(rectrack->GetTrackParamAtHit()->First());
+ AliMUONTrackExtrap::ExtrapToVertex(trackParam,0.,0.,0.);
bendingSlope = trackParam->GetBendingSlope();
nonBendingSlope = trackParam->GetNonBendingSlope();
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff