]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONv1.cxx
Update of ACORDE-QA-Raw data histograms (now they go from -0.5 to 59.5)
[u/mrichter/AliRoot.git] / MUON / AliMUONv1.cxx
index e366f8340bced8b247895a3f3298e4c74f24db3f..ce0a6974e595cda09b05f8c6e02445c3eb91e4f0 100644 (file)
@@ -1,6 +1,6 @@
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *      SigmaEffect_thetadegrees                                                                  *
+ *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
  * Contributors are mentioned in the code where appropriate.              *
  *                                                                        *
 
 /* $Id$ */
 
-/////////////////////////////////////////////////////////
-//  Manager and hits classes for set:MUON version 1    //
-/////////////////////////////////////////////////////////
-
-#include <TRandom.h>
-#include <TF1.h>
-#include <TClonesArray.h>
-#include <TRandom.h> 
-#include <TVirtualMC.h>
-#include <TGeoMatrix.h>
+//-----------------------------------------------------------------------------
+// Class AliMUONv1
+// --------------------
+// AliDetector class for MUON subsystem which implements
+// functions for simulation 
+//-----------------------------------------------------------------------------
 
 #include "AliMUONv1.h"
-#include "AliConst.h" 
-#include "AliMUONChamber.h"
 #include "AliMUONConstants.h"
-#include "AliMUONFactoryV4.h"
-#include "AliMUONFactoryV3.h"
-#include "AliMUONFactoryV2.h"
+#include "AliMUONResponseFactory.h"
 #include "AliMUONHit.h"
-#include "AliMUONTriggerCircuit.h"
 #include "AliMUONGeometryBuilder.h"    
+#include "AliMUONGeometry.h"   
+#include "AliMUONGeometryTransformer.h"        
 #include "AliMUONGeometryModule.h"     
-#include "AliMUONGeometrySVMap.h"      
+#include "AliMUONStringIntMap.h"       
 #include "AliMUONGeometryDetElement.h" 
+
+#include "AliMpCDB.h"
+#include "AliMpDEManager.h"
+
+#include "AliConst.h" 
 #include "AliMagF.h"
 #include "AliRun.h"
 #include "AliMC.h"
+#include "AliTrackReference.h"
 #include "AliLog.h"
 
+#include <TClonesArray.h>
+#include <TF1.h>
+#include <TF2.h>
+#include <TGeoGlobalMagField.h>
+#include <TGeoMatrix.h>
+#include <TRandom.h>
+#include <TRandom.h> 
+#include <TVirtualMC.h>
+
 #include <string>
 
+#include "AliMUONVHitStore.h"
+
+/// \cond CLASSIMP
 ClassImp(AliMUONv1)
+/// \endcond
  
 //___________________________________________
 AliMUONv1::AliMUONv1() 
   : AliMUON(),
-    fStepManagerVersionOld(kTRUE),
-    fStepManagerVersionDE(kTRUE),
     fAngleEffect(kTRUE),
+    fMagEffect(kTRUE),
     fStepMaxInActiveGas(0.6),
     fStepSum(0x0),
     fDestepSum(0x0),
@@ -62,18 +73,18 @@ AliMUONv1::AliMUONv1()
     fElossRatio(0x0),
     fAngleEffect10(0x0),
     fAngleEffectNorma(0x0),
-    fFactory(0x0)
+    fMagAngleEffectNorma(0x0)
 {
-// Default constructor
-       AliDebug(1,Form("default (empty) ctor this=%p",this));
+/// Default constructor
+  
+  AliDebug(1,Form("default (empty) ctor this=%p",this));
 } 
 
 //___________________________________________
-AliMUONv1::AliMUONv1(const char *name, const char *title)
-  : AliMUON(name,title), 
-    fStepManagerVersionOld(kTRUE),
-    fStepManagerVersionDE(kTRUE),
+AliMUONv1::AliMUONv1(const char *name, const char* title)
+: AliMUON(name, title), 
     fAngleEffect(kTRUE),
+    fMagEffect(kTRUE),
     fStepMaxInActiveGas(0.6),
     fStepSum(0x0),
     fDestepSum(0x0),
@@ -82,12 +93,17 @@ AliMUONv1::AliMUONv1(const char *name, const char *title)
     fElossRatio(0x0),
     fAngleEffect10(0x0),
     fAngleEffectNorma(0x0),
-    fFactory(0x0)
+    fMagAngleEffectNorma(0x0)
 {
-// Standard onstructor
+/// Standard onstructor
 
-       AliDebug(1,Form("ctor this=%p",this));  
+    AliDebug(1,Form("ctor this=%p",this));     
        
+    // Load mapping
+    if ( ! AliMpCDB::LoadMpSegmentation() ) {
+      AliFatal("Could not access mapping from OCDB !");
+    }
+       
     // By default include all stations
 
     fStepSum   = new Float_t [AliMUONConstants::NCh()];
@@ -116,49 +132,31 @@ AliMUONv1::AliMUONv1(const char *name, const char *title)
     fAngleEffectNorma->SetParameter(1,-6.809e-01);
     fAngleEffectNorma->SetParameter(2,5.151e-02);
     fAngleEffectNorma->SetParameter(3,-1.490e-03);
-}
 
-//_____________________________________________________________________________
-AliMUONv1::AliMUONv1(const AliMUONv1& right) 
-  : AliMUON(right) 
-{  
-  // copy constructor (not implemented)
-
-  AliFatal("Copy constructor not provided.");
+    // Magnetic field effect: Normalisation form theta=16 degres (eq. 10 degrees B=0) to theta between -20 and 20 (Lamia Benhabib jun 2006 )  
+    // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
+    fMagAngleEffectNorma = new TF2("MagAngleEffectNorma","121.24/(([1]+[2]*abs(y))+[3]*abs(x-[0]*y)+[4]*abs((x-[0]*y)*(x-[0]*y))+[5]*abs((x-[0]*y)*(x-[0]*y)*(x-[0]*y))+[6]*abs((x-[0]*y)*(x-[0]*y)*(x-[0]*y)*(x-[0]*y)))",-20.0,20.0,-1.,1.);
+    fMagAngleEffectNorma->SetParameters(8.6995, 25.4022, 13.8822, 2.4717, 1.1551, -0.0624, 0.0012);
 }
 
 //___________________________________________
 AliMUONv1::~AliMUONv1()
 {
-// Destructor
-       AliDebug(1,Form("dtor this=%p",this));
+/// Destructor
+  
+  AliDebug(1,Form("dtor this=%p",this));
   delete [] fStepSum;
   delete [] fDestepSum;
   delete fElossRatio;
   delete fAngleEffect10;
   delete fAngleEffectNorma; 
-  delete fFactory;
+  delete fMagAngleEffectNorma;
 }
 
-//_____________________________________________________________________________
-AliMUONv1& AliMUONv1::operator=(const AliMUONv1& right)
-{
-  // assignement operator (not implemented)
-
-  // check assignement to self
-  if (this == &right) return *this;
-
-  AliFatal("Assignement operator not provided.");
-    
-  return *this;  
-}    
-
 //__________________________________________________
 void AliMUONv1::CreateGeometry()
 {
-//
-// Construct geometry using geometry builder
-//
+/// Construct geometry using geometry builder
 
   fGeometryBuilder->CreateGeometry();
 }
@@ -166,90 +164,74 @@ void AliMUONv1::CreateGeometry()
 //________________________________________________________________
 void AliMUONv1::CreateMaterials()
 {
-//
-// Construct materials using geometry builder
-//
+/// Construct materials using geometry builder
 
   fGeometryBuilder->CreateMaterials();
 }
 
+//________________________________________________________________
+void AliMUONv1::UpdateInternalGeometry()
+{
+/// Update geometry after applying mis-alignment
+
+  // Load mapping
+  if ( ! AliMpCDB::LoadMpSegmentation() ) {
+    AliFatal("Could not access mapping from OCDB !");
+  }
+
+  fGeometryBuilder->UpdateInternalGeometry();
+}
+
+//________________________________________________________________
+void AliMUONv1::AddAlignableVolumes() const
+{
+/// Construct materials using geometry builder
+
+  GetGeometryTransformer()->AddAlignableVolumes();
+}
+
+
 //___________________________________________
 void AliMUONv1::Init()
 { 
+/// Initialize geometry
+
   AliDebug(1,"Start Init for version 1 - CPC chamber type");
-  Int_t i;
    
-  //
-  // Initialize geometry
-  //
   fGeometryBuilder->InitGeometry();
   AliDebug(1,"Finished Init for version 1 - CPC chamber type");   
-  
-  std::string ftype(GetTitle());
-  
-  if ( ftype == "default" || ftype == "AliMUONFactoryV2") 
-    {
-      fFactory = new AliMUONFactoryV2("New MUON Factory");
-      (static_cast<AliMUONFactoryV2*>(fFactory))->Build(this,"default");
-    }
-  else if ( ftype == "AliMUONFactoryV3" )
-    {
-      fFactory = new AliMUONFactoryV3("New MUON Factory");
-      (static_cast<AliMUONFactoryV3*>(fFactory))->Build(this,"default");
-    }
-  else if ( ftype == "AliMUONFactoryV4" )
-  {
-    fFactory = new AliMUONFactoryV4("New MUON Factory (w/ Trigger)");
-    (static_cast<AliMUONFactoryV4*>(fFactory))->Build(this,"default");
-  }
-  else
-    {
-      AliFatal(Form("Wrong factory type : %s",ftype.c_str()));      
-    }
-      
+
+  // Build segmentation
+  // using geometry parametrisation
   //
-  // Initialize segmentation
+  // Build response
   //
-  
-  for (i=0; i<AliMUONConstants::NCh(); i++) 
-    ( (AliMUONChamber*) (*fChambers)[i])->Init(2);// new segmentation
-  
-  
-  // trigger circuit
-  // cp 
-  for (i=0; i<AliMUONConstants::NTriggerCircuit(); i++) 
-  {
-    AliMUONTriggerCircuit* c = (AliMUONTriggerCircuit*)(fTriggerCircuits->At(i));
-    c->Init(i);
-//    c->Print();
-  }
+  AliMUONResponseFactory respFactory("default", fIsTailEffect);
+  respFactory.Build(this);
   
 }
 
 //__________________________________________________________________
-Int_t  AliMUONv1::GetChamberId(Int_t volId) const
+Int_t  AliMUONv1::GetGeomModuleId(Int_t volId) const
 {
-// Check if the volume with specified  volId is a sensitive volume (gas) 
-// of some chamber and returns the chamber number;
-// if not sensitive volume - return 0.
-// ---
-
-/*
-  for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
-    if (volId==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()) return i;
-*/
-  for (Int_t i = 1; i <= AliMUONConstants::NCh(); i++)
-    if ( ((AliMUONChamber*)(*fChambers)[i-1])->IsSensId(volId) ) return i;
-
-  return 0;
+/// Check if the volume with specified  volId is a sensitive volume (gas) 
+/// of some chamber and return the chamber number;
+/// if not sensitive volume - return 0.
+
+  for (Int_t i = 0; i < AliMUONConstants::NGeomModules(); i++) {
+    if ( GetGeometry()->GetModule(i)->IsSensitiveVolume(volId) )
+      return i;
+  }    
+
+  return -1;
 }
 
 //_______________________________________________________________________________
 TString  AliMUONv1::CurrentVolumePath() const
 {
-// Returns current volume path
-// (Could be removed when this function is available via gMC)
-// ---      
+/// Return current volume path
+/// (Could be removed when this function is available via gMC)
 
   TString path = "";
   TString name;
@@ -260,7 +242,7 @@ TString  AliMUONv1::CurrentVolumePath() const
     gMC->CurrentVolOffID(imother++, copyNo);
     TString add = "/";
     add += name;
-    add += ".";
+    add += "_";
     add += copyNo;
     path.Insert(0,add); 
   }
@@ -272,14 +254,7 @@ TString  AliMUONv1::CurrentVolumePath() const
 //_______________________________________________________________________________
 void AliMUONv1::StepManager()
 {
-  // Stepmanager for the chambers
-  // TBR
-
- if (fStepManagerVersionDE) {
-    StepManager2();
-    return;
-  }
-
+/// Step manager for the chambers
 
   // Only charged tracks
   if( !(gMC->TrackCharge()) ) return; 
@@ -287,192 +262,68 @@ void AliMUONv1::StepManager()
   
   // Only gas gap inside chamber
   // Tag chambers and record hits when track enters 
-  static Int_t   idvol=-1;
-  Int_t   iChamber=0;
-  Int_t   id=0;
+  static Int_t   idvol=-1, iEnter = 0;
   Int_t   copy;
   const  Float_t kBig = 1.e10;
-
+  static Double_t xyzEnter[3];
 
   //
   // Only gas gap inside chamber
   // Tag chambers and record hits when track enters 
-  id=gMC->CurrentVolID(copy);
-  iChamber = GetChamberId(id);
-  idvol = iChamber -1;
-
-  if (idvol == -1) return;
-
-  // Filling TrackRefs file for MUON. Our Track references are the active volume of the chambers
-  if ( (gMC->IsTrackEntering() || gMC->IsTrackExiting() ) )     
-    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
-  
-   if( gMC->IsTrackEntering() ) {
-     Float_t theta = fTrackMomentum.Theta();
-     if ((TMath::Pi()-theta)*kRaddeg>=15.) gMC->SetMaxStep(fStepMaxInActiveGas); // We use Pi-theta because z is negative
-  }
-
-   //  AliDebug(1,
-   //    Form("Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber,
-   //    ( (AliMUONChamber*)(*fChambers)[idvol])->Z()));
-   // Particule id and mass, 
-  Int_t     ipart = gMC->TrackPid();
-  Float_t   mass  = gMC->TrackMass();
-
-  fDestepSum[idvol]+=gMC->Edep();
-  // Get current particle id (ipart), track position (pos)  and momentum (mom)
-  if ( fStepSum[idvol]==0.0 )  gMC->TrackMomentum(fTrackMomentum);
-  fStepSum[idvol]+=gMC->TrackStep();
-  
-  //  if(AliDebugLevel()) {
-  //  AliDebug(1,
-  //        Form("iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g",
-  //             iChamber,ipart, fTrackMomentum.Theta()*kRaddeg, fTrackMomentum.Phi()*kRaddeg, 
-  //             mass, fStepSum[idvol], gMC->Edep()));
-  //  AliDebug(1,
-  //        Form("Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(), 
-  //             fTrackMomentum.Z()));
-  //  gMC->TrackPosition(fTrackPosition);
-  //  AliDebug(1,
-  //        Form("Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),
-  //     fTrackPosition.Z())) ;
-  //   }
-
-  // Track left chamber or StepSum larger than fStepMaxInActiveGas
-  if ( gMC->IsTrackExiting() || 
-       gMC->IsTrackStop() || 
-       gMC->IsTrackDisappeared()||
-       (fStepSum[idvol]>fStepMaxInActiveGas) ) {
-    
-    if   ( gMC->IsTrackExiting() || 
-           gMC->IsTrackStop() || 
-           gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
-
-    gMC->TrackPosition(fTrackPosition);
-    Float_t theta = fTrackMomentum.Theta();
-    Float_t phi   = fTrackMomentum.Phi();
-    
-    TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
-                               fStepSum[idvol]/2.*sin(theta)*sin(phi),
-                               fStepSum[idvol]/2.*cos(theta),0.0       );
-    //     AliDebug(1,Form("Exit: Track Position %f %f %f",fTrackPosition.X(),
-    //                     fTrackPosition.Y(),fTrackPosition.Z())) ;
-    //     AliDebug(1,Form("Exit: Track backToWire %f %f %f",backToWire.X(),
-    //                     backToWire.Y(),backToWire.Z()) ;
-    fTrackPosition-=backToWire;
-    
-    //-------------- Angle effect 
-    // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
-    
-    Float_t betaxGamma    = fTrackMomentum.P()/mass;//  pc/mc2
-    Float_t sigmaEffect10degrees;
-    Float_t sigmaEffectThetadegrees;
-    Float_t eLossParticleELossMip;
-    Float_t yAngleEffect=0.;
-    Float_t thetawires      =  TMath::Abs( TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) );// We use Pi-theta because z is negative
-
-
-    if (fAngleEffect){
-    if ( (betaxGamma >3.2)   &&  (thetawires*kRaddeg<=15.) ) {
-      betaxGamma=TMath::Log(betaxGamma);
-      eLossParticleELossMip = fElossRatio->Eval(betaxGamma);
-      // 10 degrees is a reference for a model (arbitrary)
-      sigmaEffect10degrees=fAngleEffect10->Eval(eLossParticleELossMip);// in micrometers
-      // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
-      sigmaEffectThetadegrees =  sigmaEffect10degrees/fAngleEffectNorma->Eval(thetawires*kRaddeg);  // For 5mm gap  
-      if ( (iChamber==1)  ||  (iChamber==2) )  
-        sigmaEffectThetadegrees/=(1.09833e+00+1.70000e-02*(thetawires*kRaddeg)); // The gap is different (4mm)
-      yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
-    }
-    }
+  Int_t id=gMC->CurrentVolID(copy);
+  Int_t iGeomModule = GetGeomModuleId(id);
+  if (iGeomModule == -1) return;
+
+  // Detection elements id
+  const AliMUONGeometryModule* kGeometryModule
+    = GetGeometry()->GetModule(iGeomModule);
+  AliMUONGeometryDetElement* detElement
+    = kGeometryModule->FindBySensitiveVolume(CurrentVolumePath());
     
-    // Detection elements ids
-    AliMUONGeometryModule* geometry
-      = Chamber(iChamber-1).GetGeometry();
-
-    AliMUONGeometryDetElement* detElement
-      = geometry->FindBySensitiveVolume(CurrentVolumePath());
-
-    Int_t detElemId = 0;
-    if (detElement) detElemId = detElement->GetUniqueID(); 
+  if (!detElement && iGeomModule < AliMUONConstants::NGeomModules()-2) {
+    iGeomModule++;
+    const AliMUONGeometryModule* kGeometryModule2
+      = GetGeometry()->GetModule(iGeomModule);
+    detElement 
+      = kGeometryModule2->FindBySensitiveVolume(CurrentVolumePath());
+  }    
+
+  Int_t detElemId = 0;
+  if (detElement) detElemId = detElement->GetUniqueID(); 
  
-    if (!detElemId) {
-      cerr << "Chamber id: "
-           << setw(3) << iChamber << "  "
-           << "Current SV: " 
-           <<  CurrentVolumePath() 
-          << "  detElemId: "
-           << setw(5) << detElemId 
-          << endl;
-      Double_t x, y, z;
-      gMC->TrackPosition(x, y, z);        
-      cerr << "   global position: "
-           << x << ", " << y << ", " << z
-           << endl;
-      AliWarning("DetElemId not identified.");
-    }  
-    
-    // One hit per chamber
-    GetMUONData()->AddHit(fIshunt, 
-                         gAlice->GetMCApp()->GetCurrentTrackNumber(), 
-                         iChamber, ipart,
-                         fTrackPosition.X(), 
-                         fTrackPosition.Y()+yAngleEffect, 
-                         fTrackPosition.Z(), 
-                         gMC->TrackTime(),
-                         fTrackMomentum.P(),
-                         theta, 
-                         phi, 
-                         fStepSum[idvol], 
-                         fDestepSum[idvol],                        
-                         fTrackPosition.X(),
-                         fTrackPosition.Y(),
-                         fTrackPosition.Z());
-
-//     if (AliDebugLevel()){
-//       AliDebug(1,Form("Exit: Particle exiting from chamber %d",iChamber));
-//       AliDebug(1,Form("Exit: StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]));
-//       AliDebug(1,Form("Exit: Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z())) ;
-//     }
-    fStepSum[idvol]  =0; // Reset for the next event
-    fDestepSum[idvol]=0; // Reset for the next event
-  }
-}
-
-//_______________________________________________________________________________
-void AliMUONv1::StepManager2()
-{
-  // Stepmanager for the chambers
-
-  // Only charged tracks
-  if( !(gMC->TrackCharge()) ) return; 
-  // Only charged tracks
+  if (!detElemId) {
+    AliErrorStream() 
+         << "Geometry module id: "
+        << setw(3) << iGeomModule << "  "
+        << "Current SV: " 
+        <<  CurrentVolumePath() 
+         << "  detElemId: "
+        << setw(5) << detElemId 
+         << endl;
+    Double_t x, y, z;
+    gMC->TrackPosition(x, y, z);        
+    AliErrorStream() 
+         << "  global position: "
+        << x << ", " << y << ", " << z
+        << endl;
+    AliErrorStream() << "DetElemId not identified." << endl;
+  } 
   
-  // Only gas gap inside chamber
-  // Tag chambers and record hits when track enters 
-  static Int_t   idvol=-1;
-  Int_t   iChamber=0;
-  Int_t   id=0;
-  Int_t   copy;
-  const  Float_t kBig = 1.e10;
-
-
-  //
-  // Only gas gap inside chamber
-  // Tag chambers and record hits when track enters 
-  id=gMC->CurrentVolID(copy);
-  iChamber = GetChamberId(id);
+  Int_t iChamber = AliMpDEManager::GetChamberId(detElemId) + 1; 
   idvol = iChamber -1;
-
-  if (idvol == -1) return;
-  
+    
   // Filling TrackRefs file for MUON. Our Track references are the active volume of the chambers
-  if ( (gMC->IsTrackEntering() || gMC->IsTrackExiting() ) )     
-    AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
+  if ( (gMC->IsTrackEntering() || gMC->IsTrackExiting() ) ) {
+    AliTrackReference* trackReference    
+      = AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kMUON);
+    trackReference->SetUserId(detElemId);
+  }  
   
-   if( gMC->IsTrackEntering() ) {
+  if( gMC->IsTrackEntering() ) {
      Float_t theta = fTrackMomentum.Theta();
      if ((TMath::Pi()-theta)*kRaddeg>=15.) gMC->SetMaxStep(fStepMaxInActiveGas); // We use Pi-theta because z is negative
+     iEnter = 1;
+     gMC->TrackPosition(xyzEnter[0], xyzEnter[1], xyzEnter[2]); // save coordinates of entrance point
   }
 
    //   AliDebug(1,
@@ -507,19 +358,56 @@ void AliMUONv1::StepManager2()
     if   ( gMC->IsTrackExiting() || 
            gMC->IsTrackStop() || 
            gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
+    if (fDestepSum[idvol] == 0) {
+      // AZ - no energy release
+      fStepSum[idvol] = 0; // Reset for the next event
+      iEnter = 0;
+      return; 
+    }
 
     gMC->TrackPosition(fTrackPosition);
     Float_t theta = fTrackMomentum.Theta();
     Float_t phi   = fTrackMomentum.Phi();
     
-    TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
-                               fStepSum[idvol]/2.*sin(theta)*sin(phi),
-                               fStepSum[idvol]/2.*cos(theta),0.0       );
-    //    AliDebug(1,
-    //      Form("Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()));
-    // AliDebug(1,
-    //      Form("Exit: Track backToWire %f %f %f",backToWire.X(),backToWire.Y(),backToWire.Z())) ;
-    fTrackPosition-=backToWire;
+    Int_t merge = 0;
+    Double_t xyz0[3]={0}, xyz1[3]={0}, tmp[3]={0};
+    if (gMC->IsTrackExiting() && iEnter != 0) {
+      // AZ - this code is to avoid artificial hit splitting at the
+      // "fake" boundary inside the same chamber. It will still produce 
+      // 2 hits but with the same coordinates (at the wire) to allow 
+      // their merging at the digitization level.
+
+      // Only for a track going from the entrance to the exit from the volume
+      // Get local coordinates
+      gMC->Gmtod(xyzEnter, xyz0, 1); // local coord. at the entrance
+
+      fTrackPosition.Vect().GetXYZ(tmp);
+      gMC->Gmtod(tmp, xyz1, 1); // local coord. at the exit
+      Float_t dx = xyz0[0] - xyz1[0];
+      Float_t dy = xyz0[1] - xyz1[1];
+      Float_t thLoc = TMath::ATan2 (TMath::Sqrt(dx*dx+dy*dy), TMath::Abs(xyz0[2]-xyz1[2]));
+      if (thLoc * TMath::RadToDeg() < 15) merge = 1; 
+    }
+
+    if (merge) {
+      Double_t dz = -0.5;
+      if (xyz1[2] != xyz0[2]) dz = xyz0[2] / (xyz1[2] - xyz0[2]);
+      tmp[0] = xyz0[0] - (xyz1[0] - xyz0[0]) * dz; // local coord. at the wire
+      tmp[1] = xyz0[1] - (xyz1[1] - xyz0[1]) * dz;
+      tmp[2] = xyz0[2] - (xyz1[2] - xyz0[2]) * dz;
+      gMC->Gdtom(tmp, xyz1, 1); // global coord. at the wire
+      fTrackPosition.SetXYZT(xyz1[0], xyz1[1], xyz1[2], fTrackPosition.T());
+    } else {
+      TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
+                                fStepSum[idvol]/2.*sin(theta)*sin(phi),
+                                fStepSum[idvol]/2.*cos(theta),0.0       );
+      fTrackPosition-=backToWire;
+      //printf(" %d %d %d %f %d \n", gMC->IsTrackExiting(), gMC->IsTrackStop(), gMC->IsTrackDisappeared(), fStepSum[idvol], iEnter);
+      //    AliDebug(1,
+      //     Form("Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()));
+      // AliDebug(1,
+      //     Form("Exit: Track backToWire %f %f %f",backToWire.X(),backToWire.Y(),backToWire.Z())) ;
+    }
     
     //-------------- Angle effect 
     // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
@@ -529,10 +417,13 @@ void AliMUONv1::StepManager2()
     Float_t sigmaEffectThetadegrees;
     Float_t eLossParticleELossMip;
     Float_t yAngleEffect=0.;
-    Float_t thetawires      =  TMath::Abs( TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) );// We use Pi-theta because z is negative
-
+    Float_t thetawires      =  TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) ;// We use Pi-theta because z is negative
+    Double_t bField[3] = {0};
+    fTrackPosition.Vect().GetXYZ(tmp);
+    TGeoGlobalMagField::Instance()->Field(tmp,bField);
 
-    if (fAngleEffect){
+    if (fAngleEffect && !fMagEffect){
+      thetawires = TMath::Abs(thetawires);
     if ( (betaxGamma >3.2)   &&  (thetawires*kRaddeg<=15.) ) {
       betaxGamma=TMath::Log(betaxGamma);
       eLossParticleELossMip = fElossRatio->Eval(betaxGamma);
@@ -545,55 +436,40 @@ void AliMUONv1::StepManager2()
       yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
     }
     }
-    
-    // Detection elements ids
-    AliMUONGeometryModule* geometry
-      = Chamber(iChamber-1).GetGeometry();
-
-    AliMUONGeometryDetElement* detElement
-      = geometry->FindBySensitiveVolume(CurrentVolumePath());
+    else if (fAngleEffect && fMagEffect) {
+      if ( (betaxGamma >3.2)   &&  (TMath::Abs(thetawires*kRaddeg)<=15.) ) {
+       betaxGamma=TMath::Log(betaxGamma);
+       eLossParticleELossMip = fElossRatio->Eval(betaxGamma);
+       // 10 degrees is a reference for a model (arbitrary)
+       sigmaEffect10degrees=fAngleEffect10->Eval(eLossParticleELossMip);// in micrometers
+       // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
+       sigmaEffectThetadegrees =  sigmaEffect10degrees/fMagAngleEffectNorma->Eval(thetawires*kRaddeg,bField[0]/10.);  // For 5mm gap  
+      if ( (iChamber==1)  ||  (iChamber==2) )  
+        sigmaEffectThetadegrees/=(1.09833e+00+1.70000e-02*(thetawires*kRaddeg)); // The gap is different (4mm)
+       yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
+      }
+    }
 
-    Int_t detElemId = 0;
-    if (detElement) detElemId = detElement->GetUniqueID(); 
-    if (!detElemId) {
-      cerr << "Chamber id: "
-           << setw(3) << iChamber << "  "
-           << "Current SV: " 
-           <<  CurrentVolumePath() 
-          << "  detElemId: "
-           << setw(5) << detElemId 
-          << endl;
-      Double_t x, y, z;
-      gMC->TrackPosition(x, y, z);        
-      cerr << "   global position: "
-           << x << ", " << y << ", " << z
-           << endl;
-      AliError("DetElemId not identified.");
-    }  
+    AliMUONHit hit(fIshunt, 
+                  gAlice->GetMCApp()->GetCurrentTrackNumber(), 
+                  detElemId, ipart,
+                  fTrackPosition.X(), 
+                  fTrackPosition.Y()+yAngleEffect, 
+                  fTrackPosition.Z(), 
+                  gMC->TrackTime(),
+                  fTrackMomentum.P(),
+                  theta, 
+                  phi, 
+                  fStepSum[idvol], 
+                  fDestepSum[idvol],                        
+                  fTrackPosition.X(),
+                  fTrackPosition.Y(),
+                  fTrackPosition.Z());
     
-    // One hit per chamber
-    GetMUONData()->AddHit2(fIshunt, 
-                         gAlice->GetMCApp()->GetCurrentTrackNumber(), 
-                         detElemId, ipart,
-                         fTrackPosition.X(), 
-                         fTrackPosition.Y()+yAngleEffect, 
-                         fTrackPosition.Z(), 
-                         gMC->TrackTime(),
-                         fTrackMomentum.P(),
-                         theta, 
-                         phi, 
-                         fStepSum[idvol], 
-                         fDestepSum[idvol],                        
-                         fTrackPosition.X(),
-                         fTrackPosition.Y(),
-                         fTrackPosition.Z());
-
-    //       AliDebug(1,Form("Exit: Particle exiting from chamber %d",iChamber));
-    //       AliDebug(1,Form("Exit: StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]));
-    //       AliDebug(1,Form("Exit: Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
+    fHitStore->Add(hit);
 
     fStepSum[idvol]  =0; // Reset for the next event
     fDestepSum[idvol]=0; // Reset for the next event
+    iEnter = 0;
   }
 }