From 8cde4af5f2c90b69798a53593a37b9c8b1e835f8 Mon Sep 17 00:00:00 2001 From: ivana Date: Fri, 9 Mar 2007 15:53:19 +0000 Subject: [PATCH] (a) Multiple scattering and energy loss in absorber are now calculated 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.) --- MUON/AliMUONDisplay.cxx | 14 +- MUON/AliMUONReconstructor.cxx | 36 +- MUON/AliMUONTrack.h | 4 +- MUON/AliMUONTrackExtrap.cxx | 592 ++++++++++++++-------------- MUON/AliMUONTrackExtrap.h | 11 +- MUON/AliMUONTrackK.cxx | 284 +------------ MUON/AliMUONTrackK.h | 4 - MUON/AliMUONTrackReconstructor.cxx | 76 +--- MUON/AliMUONTrackReconstructor.h | 1 - MUON/AliMUONTrackReconstructorK.cxx | 43 +- MUON/AliMUONTrackReconstructorK.h | 2 - MUON/AliMUONVTrackReconstructor.h | 1 - MUON/MUONAlignment.C | 4 +- MUON/MUONCheckDI.C | 2 +- MUON/MUONRecoCheck.C | 7 +- MUON/MUONefficiency.C | 68 +++- MUON/MUONmassPlot_ESD.C | 73 ++-- MUON/MUONmassPlot_NewIO.C | 7 +- 18 files changed, 477 insertions(+), 752 deletions(-) diff --git a/MUON/AliMUONDisplay.cxx b/MUON/AliMUONDisplay.cxx index 7522329ffa0..9205a0dffa5 100644 --- a/MUON/AliMUONDisplay.cxx +++ b/MUON/AliMUONDisplay.cxx @@ -1339,11 +1339,8 @@ void AliMUONDisplay::LoadTracks() 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++){ @@ -1365,13 +1362,13 @@ void AliMUONDisplay::LoadTracks() 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(); @@ -1385,6 +1382,9 @@ void AliMUONDisplay::PrintTrack(Int_t iRecTracks, AliMUONTrack *recTrack) 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]); diff --git a/MUON/AliMUONReconstructor.cxx b/MUON/AliMUONReconstructor.cxx index 80cbdff1991..146c49fe743 100644 --- a/MUON/AliMUONReconstructor.cxx +++ b/MUON/AliMUONReconstructor.cxx @@ -353,7 +353,7 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, } recModel->SetGhostChi2Cut(10); recModel->SetEventNumber(evtNumber); - + TTask* calibration = GetCalibrationTask(); loader->LoadRecPoints("RECREATE"); @@ -381,7 +381,7 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliDebug(1,Form("Event %d",iEvent)); runLoader->GetEvent(iEvent++); - + //----------------------- raw2digits & raw2trigger------------------- // if (!loader->TreeD()) // { @@ -411,7 +411,7 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, clusterTimer.Start(kFALSE); if (!loader->TreeR()) loader->MakeRecPointsContainer(); - + // tracking branch fMUONData->MakeBranch("RC"); fMUONData->SetTreeAddress("RC"); @@ -426,7 +426,7 @@ void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, loader->WriteRecPoints("OVERWRITE"); clusterTimer.Stop(); - + //--------------------------- Resetting branches ----------------------- fMUONData->ResetDigits(); @@ -475,7 +475,7 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const // setting pointer for tracks, triggertracks & trackparam at vertex AliMUONTrack* recTrack = 0; - AliMUONTrackParam* trackParam = 0; + AliMUONTrackParam trackParam; AliMUONTriggerTrack* recTriggerTrack = 0; iEvent = runLoader->GetEventNumber(); @@ -519,17 +519,23 @@ void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const // 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(); diff --git a/MUON/AliMUONTrack.h b/MUON/AliMUONTrack.h index f2cec41d00a..e0e50e7e05c 100644 --- a/MUON/AliMUONTrack.h +++ b/MUON/AliMUONTrack.h @@ -102,7 +102,7 @@ class AliMUONTrack : public TObject 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 @@ -120,7 +120,7 @@ class AliMUONTrack : public TObject 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 diff --git a/MUON/AliMUONTrackExtrap.cxx b/MUON/AliMUONTrackExtrap.cxx index 42fe6ffb4db..49b17a39a45 100644 --- a/MUON/AliMUONTrackExtrap.cxx +++ b/MUON/AliMUONTrackExtrap.cxx @@ -28,16 +28,16 @@ // /////////////////////////////////////////////////// -#include -#include -#include - #include "AliMUONTrackExtrap.h" #include "AliMUONTrackParam.h" #include "AliMUONConstants.h" + #include "AliMagF.h" -#include "AliLog.h" -#include "AliTracker.h" + +#include +#include +#include +#include ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context @@ -335,57 +335,76 @@ void AliMUONTrackExtrap::ExtrapToStation(AliMUONTrackParam* trackParamIn, Int_t //__________________________________________________________________________ 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 ("<GetZ() <<") upstream the vertex (zVtx = "<GetVolume("ALIC")->GetNode("ABSM_1"); + if (!absNode) { + cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: failed to get absorber node"<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 ("<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 ("<GetZ() + <<") upstream the front absorber (zAbsorberBegin = "<GetZ() > zAbsEnd) { // spectro. (z<0) + cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() + <<") inside the front absorber ("<GetZ() - <<") upstream or inside the front absorber (zAbsorberEnd = "<=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); @@ -393,7 +412,129 @@ void AliMUONTrackExtrap::ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam } //__________________________________________________________________________ -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"<InitTrack(trackXYZIn, b); + if (!currentnode) { + cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<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"<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"< 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" @@ -432,276 +573,130 @@ void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam, Double_t { /// 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 ("<GetZ() + <<") upstream the vertex (zVtx = "<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"<GetVolume("ALIC")->GetNode("ABSM_1"); + if (!absNode) { + cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: failed to get absorber node"<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 ("<GetZ() > zAbsBeg) { // spectro. (z<0) + cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() + <<") upstream the front absorber (zAbsorberBegin = "<GetZ() > zAbsEnd) { // spectro. (z<0) + cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<GetZ() + <<") inside the front absorber ("<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"<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); } //__________________________________________________________________________ @@ -896,6 +891,7 @@ void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, Doub return; } + //__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout) { @@ -1139,6 +1135,7 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, return; } + //___________________________________________________________ void AliMUONTrackExtrap::GetField(Double_t *Position, Double_t *Field) { @@ -1157,4 +1154,3 @@ void AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, return; } - diff --git a/MUON/AliMUONTrackExtrap.h b/MUON/AliMUONTrackExtrap.h index f259c7e064c..7ffdc54e87f 100644 --- a/MUON/AliMUONTrackExtrap.h +++ b/MUON/AliMUONTrackExtrap.h @@ -39,7 +39,7 @@ class AliMUONTrackExtrap : public TObject 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); @@ -61,9 +61,12 @@ class AliMUONTrackExtrap : public TObject 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); diff --git a/MUON/AliMUONTrackK.cxx b/MUON/AliMUONTrackK.cxx index 56d5ccb58c9..3d624198201 100644 --- a/MUON/AliMUONTrackK.cxx +++ b/MUON/AliMUONTrackK.cxx @@ -23,7 +23,7 @@ // Author: Alexander Zinchenko, JINR Dubna #include "AliMUONTrackK.h" -#include "AliMUON.h" +#include "AliMUONData.h" #include "AliMUONConstants.h" #include "AliMUONTrackReconstructorK.h" @@ -36,8 +36,6 @@ #include "AliMUONDetElement.h" #include "AliLog.h" -#include "AliMagF.h" -#include "AliRunLoader.h" #include #include @@ -1346,14 +1344,12 @@ void AliMUONTrackK::FillMUONTrack(void) /// 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 @@ -1381,280 +1377,6 @@ void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, trackParam->SetZ(z); } - //__________________________________________________________________________ -void AliMUONTrackK::Branson(void) -{ -/// Propagates track to the vertex thru absorber using Branson correction -/// (makes use of the AliMUONTrackParam class) - - //AliMUONTrackParam *trackParam = new AliMUONTrackParam(); - AliMUONTrackParam trackParam = AliMUONTrackParam(); - /* - trackParam->SetBendingCoor((*fTrackPar)(0,0)); - trackParam->SetNonBendingCoor((*fTrackPar)(1,0)); - trackParam->SetBendingSlope(TMath::Tan((*fTrackPar)(2,0))); - trackParam->SetNonBendingSlope(TMath::Tan((*fTrackPar)(3,0))/TMath::Cos((*fTrackPar)(2,0))); - trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0))); - trackParam->SetZ(fPosition); - */ - SetTrackParam(&trackParam, fTrackPar, fPosition); - - 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; i1GetChamberNumber()); - } - cout << endl; - if (fgDebug > 0) { - for (Int_t i1=0; i1GetHitNumber() << " "; - //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " "; - printf ("%5d", fgHitForRec->IndexOf(hit)); - } - cout << endl; - } - - // from raw clusters - for (Int_t i1=0; i1GetHitNumber() < 0) { // combined cluster / track finder - Int_t index = -hit->GetHitNumber() / 100000; - Int_t iPos = -hit->GetHitNumber() - index * 100000; - clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos); - } else { - rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber()); - clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber()); - } - printf ("%5d", clus->GetTrack(1)%10000000); - } - cout << endl; - for (Int_t i1=0; i1GetHitNumber() < 0) { // combined cluster / track finder - Int_t index = -hit->GetHitNumber() / 100000; - Int_t iPos = -hit->GetHitNumber() - index * 100000; - clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos); - } else { - rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber()); - clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber()); - } - if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000); - else printf ("%5s", " "); - } - cout << endl; - for (Int_t i1=0; i1GetHitNumber() << " "; - cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " "; - //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " "; - } - cout << endl; - for (Int_t i1=0; i1Determinant() != 0) { - Int_t ifail; - mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail); - } else { - AliWarning(" Determinant fCovariance=0:" ); - } - */ -} - //__________________________________________________________________________ void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0) { diff --git a/MUON/AliMUONTrackK.h b/MUON/AliMUONTrackK.h index 1414898b90c..c275c960730 100644 --- a/MUON/AliMUONTrackK.h +++ b/MUON/AliMUONTrackK.h @@ -23,7 +23,6 @@ class TClonesArray; class TObjArray; #include -#include class AliMUONTrackK : public AliMUONTrack { @@ -54,9 +53,6 @@ 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 diff --git a/MUON/AliMUONTrackReconstructor.cxx b/MUON/AliMUONTrackReconstructor.cxx index 0f1bec521a5..b933308a831 100644 --- a/MUON/AliMUONTrackReconstructor.cxx +++ b/MUON/AliMUONTrackReconstructor.cxx @@ -111,13 +111,12 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters() 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(); @@ -152,8 +151,6 @@ void AliMUONTrackReconstructor::MakeTracks(void) FollowTracks(); // Remove double tracks RemoveDoubleTracks(); - // Propagate tracks to the vertex through absorber - ExtrapTracksToVertex(); // Fill out the AliMUONTrack's FillMUONTrack(); } @@ -189,7 +186,7 @@ void AliMUONTrackReconstructor::MakeTrackCandidates(void) 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<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<mncler(); // Give the fitted track to MINUIT @@ -674,7 +673,7 @@ void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool // 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); @@ -933,40 +932,6 @@ Double_t MultipleScatteringAngle2(AliMUONTrackParam *param) // 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 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(); - } - } - } //__________________________________________________________________________ @@ -993,7 +958,7 @@ void AliMUONTrackReconstructor::EventDump(void) /// 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; @@ -1010,17 +975,18 @@ void AliMUONTrackReconstructor::EventDump(void) 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)); diff --git a/MUON/AliMUONTrackReconstructor.h b/MUON/AliMUONTrackReconstructor.h index eb9aacc9350..df6f017a921 100644 --- a/MUON/AliMUONTrackReconstructor.h +++ b/MUON/AliMUONTrackReconstructor.h @@ -30,7 +30,6 @@ class AliMUONTrackReconstructor : public AliMUONVTrackReconstructor { virtual void MakeTrackCandidates(void); virtual void FollowTracks(void); virtual void RemoveDoubleTracks(void); - virtual void ExtrapTracksToVertex(void); virtual void FillMUONTrack(void); diff --git a/MUON/AliMUONTrackReconstructorK.cxx b/MUON/AliMUONTrackReconstructorK.cxx index 1979e58667a..3e0a554c9d7 100644 --- a/MUON/AliMUONTrackReconstructorK.cxx +++ b/MUON/AliMUONTrackReconstructorK.cxx @@ -28,13 +28,6 @@ /// //////////////////////////////////// -#include -#include -#include -#include -#include - -#include "AliMUONVTrackReconstructor.h" #include "AliMUONTrackReconstructorK.h" #include "AliMUONData.h" #include "AliMUONConstants.h" @@ -42,9 +35,15 @@ #include "AliMUONObjectPair.h" #include "AliMUONRawCluster.h" #include "AliMUONTrackK.h" + #include "AliLog.h" +#include + +/// \cond CLASSIMP ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context +ClassImp(AliMUONConstants) +/// \endcond //__________________________________________________________________________ AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONData* data, const Option_t* TrackMethod) @@ -114,13 +113,12 @@ void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters() 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(); @@ -157,8 +155,6 @@ void AliMUONTrackReconstructorK::MakeTracks(void) FollowTracks(); // Remove double tracks RemoveDoubleTracks(); - // Propagate tracks to the vertex through absorber - ExtrapTracksToVertex(); // Fill AliMUONTrack data members FillMUONTrack(); } @@ -394,21 +390,6 @@ void AliMUONTrackReconstructorK::RemoveDoubleTracks(void) 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; iBranson(); - ((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() { diff --git a/MUON/AliMUONTrackReconstructorK.h b/MUON/AliMUONTrackReconstructorK.h index ae29f581c6f..ddb2f9f975e 100644 --- a/MUON/AliMUONTrackReconstructorK.h +++ b/MUON/AliMUONTrackReconstructorK.h @@ -13,7 +13,6 @@ /// MUON track reconstructor using kalman filter //////////////////////////////////////////////// -#include #include "AliMUONVTrackReconstructor.h" class AliMUONTrackReconstructorK : public AliMUONVTrackReconstructor { @@ -36,7 +35,6 @@ class AliMUONTrackReconstructorK : public AliMUONVTrackReconstructor { virtual void MakeTrackCandidates(void); virtual void FollowTracks(void); virtual void RemoveDoubleTracks(void); - virtual void ExtrapTracksToVertex(void); virtual void FillMUONTrack(void); diff --git a/MUON/AliMUONVTrackReconstructor.h b/MUON/AliMUONVTrackReconstructor.h index 592fd8ec86a..6b645123e21 100644 --- a/MUON/AliMUONVTrackReconstructor.h +++ b/MUON/AliMUONVTrackReconstructor.h @@ -103,7 +103,6 @@ class AliMUONVTrackReconstructor : public TObject { 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: diff --git a/MUON/MUONAlignment.C b/MUON/MUONAlignment.C index 96bd50e9f14..827e020ec8f 100644 --- a/MUON/MUONAlignment.C +++ b/MUON/MUONAlignment.C @@ -103,7 +103,9 @@ void MUONAlignment(Int_t nEvents = 100000, TString fileName = "galice.root", TSt 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) { diff --git a/MUON/MUONCheckDI.C b/MUON/MUONCheckDI.C index 5b192d5f889..0b5b5e54c3d 100644 --- a/MUON/MUONCheckDI.C +++ b/MUON/MUONCheckDI.C @@ -328,7 +328,7 @@ void MUONRecTracks (char * filename="galice.root") 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(); diff --git a/MUON/MUONRecoCheck.C b/MUON/MUONRecoCheck.C index 5a32d881745..31674e10d16 100644 --- a/MUON/MUONRecoCheck.C +++ b/MUON/MUONRecoCheck.C @@ -33,6 +33,7 @@ #include "AliMUONTrack.h" #include "AliMUONRecoCheck.h" #include "AliMUONTrackParam.h" +#include "AliMUONTrackExtrap.h" #include "AliTracker.h" Int_t TrackCheck( Bool_t *compTrack); @@ -165,8 +166,9 @@ void MUONRecoCheck (Int_t nEvent = 1, char * filename="galice.root"){ 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(); @@ -174,6 +176,7 @@ void MUONRecoCheck (Int_t nEvent = 1, char * filename="galice.root"){ 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); diff --git a/MUON/MUONefficiency.C b/MUON/MUONefficiency.C index e6adc8fcc66..74b6f421bee 100644 --- a/MUON/MUONefficiency.C +++ b/MUON/MUONefficiency.C @@ -50,6 +50,7 @@ #include "TTree.h" #include "TString.h" #include +#include // STEER includes #include "AliRun.h" @@ -67,9 +68,16 @@ #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; @@ -177,6 +185,15 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven 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"); @@ -320,13 +337,17 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven // 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 @@ -387,34 +408,39 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven // 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(); @@ -464,14 +490,16 @@ Bool_t MUONefficiency( Int_t ResType = 553, Int_t FirstEvent = 0, Int_t LastEven 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); diff --git a/MUON/MUONmassPlot_ESD.C b/MUON/MUONmassPlot_ESD.C index de460fd1824..36a0766a78a 100644 --- a/MUON/MUONmassPlot_ESD.C +++ b/MUON/MUONmassPlot_ESD.C @@ -10,6 +10,7 @@ #include "TParticle.h" #include "TTree.h" #include +#include // STEER includes #include "AliRun.h" @@ -40,6 +41,10 @@ // 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) @@ -57,8 +62,8 @@ // 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) { @@ -126,6 +131,15 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t 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"); @@ -135,8 +149,7 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t // 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; } @@ -183,12 +196,10 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t // 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); @@ -202,17 +213,22 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t // 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); @@ -258,28 +274,33 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t // 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(); @@ -310,7 +331,7 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t 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()); @@ -319,8 +340,10 @@ Bool_t MUONmassPlot(char* filename = "galice.root", Int_t FirstEvent = 0, Int_t } // 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); diff --git a/MUON/MUONmassPlot_NewIO.C b/MUON/MUONmassPlot_NewIO.C index 08332322d0d..872cadbcc43 100644 --- a/MUON/MUONmassPlot_NewIO.C +++ b/MUON/MUONmassPlot_NewIO.C @@ -26,6 +26,7 @@ #include "AliMUONLocalTrigger.h" #include "AliMUONTrack.h" #include "AliMUONTrackParam.h" +#include "AliMUONTrackExtrap.h" #include "AliESDMuonTrack.h" #endif // @@ -152,7 +153,8 @@ void MUONmassPlot(char* filename="galice.root", Int_t FirstEvent = 0, Int_t Last 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(); @@ -198,7 +200,8 @@ void MUONmassPlot(char* filename="galice.root", Int_t FirstEvent = 0, Int_t Last 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(); -- 2.43.0