(a) Multiple scattering and energy loss in absorber are now calculated
authorivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 9 Mar 2007 15:53:19 +0000 (15:53 +0000)
committerivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 9 Mar 2007 15:53:19 +0000 (15:53 +0000)
    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.)

18 files changed:
MUON/AliMUONDisplay.cxx
MUON/AliMUONReconstructor.cxx
MUON/AliMUONTrack.h
MUON/AliMUONTrackExtrap.cxx
MUON/AliMUONTrackExtrap.h
MUON/AliMUONTrackK.cxx
MUON/AliMUONTrackK.h
MUON/AliMUONTrackReconstructor.cxx
MUON/AliMUONTrackReconstructor.h
MUON/AliMUONTrackReconstructorK.cxx
MUON/AliMUONTrackReconstructorK.h
MUON/AliMUONVTrackReconstructor.h
MUON/MUONAlignment.C
MUON/MUONCheckDI.C
MUON/MUONRecoCheck.C
MUON/MUONefficiency.C
MUON/MUONmassPlot_ESD.C
MUON/MUONmassPlot_NewIO.C

index 7522329ffa04fb1072225938950286829b6c8a84..9205a0dffa5243335741e567114a4d4bc7283e01 100644 (file)
@@ -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]); 
index 80cbdff1991218fcd1904db92442432075fdd997..146c49fe743fa12b99078f4bd0c6c5c9e667a6cb 100644 (file)
@@ -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();
index f2cec41d00ae0c5cb567cc0ecc9a49f719b33f87..e0e50e7e05cf1df24e25d0172fed551ac662b1f5 100644 (file)
@@ -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
index 42fe6ffb4db9f55e755c34f4e88ae4ce059bafd1..49b17a39a4537306de71a35153b7f58f03930f61 100644 (file)
 //
 ///////////////////////////////////////////////////
 
-#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
 
@@ -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 ("<<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);
@@ -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"<<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"
@@ -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 ("<<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);
 }
 
  //__________________________________________________________________________
@@ -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;
 }
-
index f259c7e064c4bb0bf2721f99a71e39f49c7465d6..7ffdc54e87f1bc7f4cd9f2a838609e01db9f0e29 100644 (file)
@@ -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);
index 56d5ccb58c90d73bd98f0d1d6c8bc229a9e2f1f1..3d624198201a2d1a5b102487c35feba24a4bd343 100644 (file)
@@ -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 <Riostream.h>
 #include <TClonesArray.h>
@@ -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; 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)
 {
index 1414898b90cabfa0c6fecfcee3a2eccbf1f09e77..c275c960730ac29339779b4e30d72140e9a6db3c 100644 (file)
@@ -23,7 +23,6 @@ class TClonesArray;
 class TObjArray;
 
 #include <TMatrixD.h>
-#include <TObject.h>
 
 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
index 0f1bec521a52d4454ba55c9da10c18ba140ed27c..b933308a831afca8297d5271eab167f7fe46432b 100644 (file)
@@ -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<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
@@ -348,7 +345,7 @@ void AliMUONTrackReconstructor::FollowTracks(void)
       }
       // 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;
@@ -634,6 +631,8 @@ void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool
   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
@@ -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 <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();
-    }
-  }
-  
 }
 
   //__________________________________________________________________________
@@ -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));
 
index eb9aacc9350ad0795262de65cff130ca4412f1d0..df6f017a9216d6897a7c1e0dc1aadb82aa5dcef8 100644 (file)
@@ -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);
   
 
index 1979e58667a11d844bc9b9d87bad8c6ac34baa9e..3e0a554c9d7b5bb67eeff53c0d49b10514e521c0 100644 (file)
 ///
 ////////////////////////////////////
 
-#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)
@@ -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; 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()
 {
index ae29f581c6fc095cf21e458d44c2ad60b72cf052..ddb2f9f975e96b2c1826aef1ca58005b76f30782 100644 (file)
@@ -13,7 +13,6 @@
 /// MUON track reconstructor using kalman filter
 ////////////////////////////////////////////////
 
-#include <TObject.h>
 #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);
   
 
index 592fd8ec86a092c5dc55eeaf07719aebc42070eb..6b645123e21fa8aa30b4d301adc7f70e2644567c 100644 (file)
@@ -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:
index 96bd50e9f149cdbf09cf28045ae46fc024c3a540..827e020ec8f3baa47a53ff4e5cf9fba98262880e 100644 (file)
@@ -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) {
index 5b192d5f8890a0a46f219c479249df2b4f4e6679..0b5b5e54c3db03673f4791a3ac9b643c43d539b6 100644 (file)
@@ -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();
index 5a32d881745df12e96716ada484306c8f9152cc1..31674e10d167931e8bf0800f44a55e501b5f807f 100644 (file)
@@ -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);
index e6adc8fcc66c4577179c34657dec623f16fd3d8e..74b6f421beedfbe00b7e97f605f01de6c972da85 100644 (file)
@@ -50,6 +50,7 @@
 #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;
@@ -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);
index de460fd182489beddec64f4e36063e175ec8ccf6..36a0766a78ae9bc189e3e723f092aaf4e72381a4 100644 (file)
@@ -10,6 +10,7 @@
 #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)
@@ -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);
index 08332322d0d756a4c7cd91e657ff8d15bd6673b2..872cadbcc43f1e180aec8368ddb83c748b45d32f 100644 (file)
@@ -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();