MUON/AliMUONv1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *      SigmaEffect_thetadegrees                                                                  *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpeateose. It is      *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////
  19 //  Manager and hits classes for set:MUON version 1    //
  20 /////////////////////////////////////////////////////////
  21
  22 #include <TRandom.h>
  23 #include <TF1.h>
  24 #include <TClonesArray.h>
  25 #include <TRandom.h>
  26 #include <TVirtualMC.h>
  27 #include <TGeoMatrix.h>
  28
  29 #include "AliMUONv1.h"
  30 #include "AliConst.h"
  31 #include "AliMUONConstants.h"
  32 #include "AliMUONSegFactory.h"
  33 #include "AliMUONResponseFactory.h"
  34 #include "AliMUONSegmentation.h"
  35 #include "AliMUONHit.h"
  36 #include "AliMUONTriggerCircuit.h"
  37 #include "AliMUONGeometryBuilder.h"
  38 #include "AliMUONGeometry.h"
  39 #include "AliMUONGeometryTransformer.h"
  40 #include "AliMUONGeometryModule.h"
  41 #include "AliMUONGeometrySVMap.h"
  42 #include "AliMUONGeometryDetElement.h"
  43 #include "AliMagF.h"
  44 #include "AliRun.h"
  45 #include "AliMC.h"
  46 #include "AliLog.h"
  47
  48 #include <string>
  49
  50 ClassImp(AliMUONv1)
  51
  52 //___________________________________________
  53 AliMUONv1::AliMUONv1()
  54   : AliMUON(),
  55     fStepManagerVersionOld(kTRUE),
  56     fStepManagerVersionDE(kTRUE),
  57     fAngleEffect(kTRUE),
  58     fStepMaxInActiveGas(0.6),
  59     fStepSum(0x0),
  60     fDestepSum(0x0),
  61     fTrackMomentum(),
  62     fTrackPosition(),
  63     fElossRatio(0x0),
  64     fAngleEffect10(0x0),
  65     fAngleEffectNorma(0x0)
  66 {
  67 // Default constructor
  68         AliDebug(1,Form("default (empty) ctor this=%p",this));
  69 }
  70
  71 //___________________________________________
  72 AliMUONv1::AliMUONv1(const char *name, const char *title,
  73                      const char* sDigitizerClassName,
  74                      const char* digitizerClassName)
  75 : AliMUON(name,title,sDigitizerClassName,digitizerClassName),
  76     fStepManagerVersionOld(kTRUE),
  77     fStepManagerVersionDE(kTRUE),
  78     fAngleEffect(kTRUE),
  79     fStepMaxInActiveGas(0.6),
  80     fStepSum(0x0),
  81     fDestepSum(0x0),
  82     fTrackMomentum(),
  83     fTrackPosition(),
  84     fElossRatio(0x0),
  85     fAngleEffect10(0x0),
  86     fAngleEffectNorma(0x0)
  87 {
  88 // Standard onstructor
  89
  90         AliDebug(1,Form("ctor this=%p",this));
  91
  92     // By default include all stations
  93
  94     fStepSum   = new Float_t [AliMUONConstants::NCh()];
  95     fDestepSum = new Float_t [AliMUONConstants::NCh()];
  96     for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {
  97       fStepSum[i] =0.0;
  98       fDestepSum[i]=0.0;
  99     }
 100     // Ratio of particle mean eloss with respect MIP's Khalil Boudjemline, sep 2003, PhD.Thesis and Particle Data Book
 101     fElossRatio = new TF1("ElossRatio","[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x",0.5,5.);
 102     fElossRatio->SetParameter(0,1.02138);
 103     fElossRatio->SetParameter(1,-9.54149e-02);
 104     fElossRatio->SetParameter(2,+7.83433e-02);
 105     fElossRatio->SetParameter(3,-9.98208e-03);
 106     fElossRatio->SetParameter(4,+3.83279e-04);
 107
 108     // Angle effect in tracking chambers at theta =10 degres as a function of ElossRatio (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis) (in micrometers)
 109     fAngleEffect10 = new TF1("AngleEffect10","[0]+[1]*x+[2]*x*x",0.5,3.0);
 110     fAngleEffect10->SetParameter(0, 1.90691e+02);
 111     fAngleEffect10->SetParameter(1,-6.62258e+01);
 112     fAngleEffect10->SetParameter(2,+1.28247e+01);
 113     // Angle effect: Normalisation form theta=10 degres to theta between 0 and 10 (Khalil BOUDJEMLINE sep 2003 Ph.D Thesis)
 114     // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
 115     fAngleEffectNorma = new TF1("AngleEffectNorma","[0]+[1]*x+[2]*x*x+[3]*x*x*x",0.0,10.0);
 116     fAngleEffectNorma->SetParameter(0,4.148);
 117     fAngleEffectNorma->SetParameter(1,-6.809e-01);
 118     fAngleEffectNorma->SetParameter(2,5.151e-02);
 119     fAngleEffectNorma->SetParameter(3,-1.490e-03);
 120 }
 121
 122 //_____________________________________________________________________________
 123 AliMUONv1::AliMUONv1(const AliMUONv1& right)
 124   : AliMUON(right)
 125 {
 126   // copy constructor (not implemented)
 127
 128   AliFatal("Copy constructor not provided.");
 129 }
 130
 131 //___________________________________________
 132 AliMUONv1::~AliMUONv1()
 133 {
 134 // Destructor
 135         AliDebug(1,Form("dtor this=%p",this));
 136   delete [] fStepSum;
 137   delete [] fDestepSum;
 138   delete fElossRatio;
 139   delete fAngleEffect10;
 140   delete fAngleEffectNorma;
 141 }
 142
 143 //_____________________________________________________________________________
 144 AliMUONv1& AliMUONv1::operator=(const AliMUONv1& right)
 145 {
 146   // assignement operator (not implemented)
 147
 148   // check assignement to self
 149   if (this == &right) return *this;
 150
 151   AliFatal("Assignement operator not provided.");
 152
 153   return *this;
 154 }
 155
 156 //__________________________________________________
 157 void AliMUONv1::CreateGeometry()
 158 {
 159 //
 160 // Construct geometry using geometry builder
 161 //
 162
 163   fGeometryBuilder->CreateGeometry();
 164 }
 165
 166 //________________________________________________________________
 167 void AliMUONv1::CreateMaterials()
 168 {
 169 //
 170 // Construct materials using geometry builder
 171 //
 172
 173   fGeometryBuilder->CreateMaterials();
 174 }
 175
 176 //___________________________________________
 177 void AliMUONv1::Init()
 178 {
 179   AliDebug(1,"Start Init for version 1 - CPC chamber type");
 180
 181   //
 182   // Initialize geometry
 183   //
 184   fGeometryBuilder->InitGeometry();
 185   AliDebug(1,"Finished Init for version 1 - CPC chamber type");
 186
 187
 188   std::string ftype(GetTitle());
 189
 190   // Build segmentation
 191   // using geometry parametrisation
 192   //
 193   AliMUONSegFactory segFactory(GetGeometryTransformer());
 194   fSegmentation = segFactory.CreateSegmentation(ftype);
 195
 196   if (!fSegmentation) {
 197     AliFatal(Form("Wrong factory type : %s",ftype.c_str()));
 198   }
 199
 200   // Build response
 201   //
 202   AliMUONResponseFactory respFactory("default");
 203   respFactory.Build(this);
 204
 205
 206   // Initialize segmentation
 207   //
 208   fSegmentation->Init();
 209
 210   // Initialize trigger circuits
 211   //
 212   for (Int_t i=0; i<AliMUONConstants::NTriggerCircuit(); i++)  {
 213     AliMUONTriggerCircuit* c = (AliMUONTriggerCircuit*)(fTriggerCircuits->At(i));
 214     c->Init(i);
 215 //    c->Print();
 216   }
 217
 218 }
 219
 220 //__________________________________________________________________
 221 Int_t  AliMUONv1::GetChamberId(Int_t volId) const
 222 {
 223 // Check if the volume with specified  volId is a sensitive volume (gas)
 224 // of some chamber and returns the chamber number;
 225 // if not sensitive volume - return 0.
 226 // ---
 227
 228   for (Int_t i = 0; i < AliMUONConstants::NCh(); i++) {
 229     if ( GetGeometry()->GetModule(i)->IsSensitiveVolume(volId) )
 230       return i+1;
 231   }
 232
 233   return 0;
 234 }
 235
 236 //_______________________________________________________________________________
 237 TString  AliMUONv1::CurrentVolumePath() const
 238 {
 239 // Returns current volume path
 240 // (Could be removed when this function is available via gMC)
 241 // ---
 242
 243   TString path = "";
 244   TString name;
 245   Int_t copyNo;
 246   Int_t imother = 0;
 247   do {
 248     name = gMC->CurrentVolOffName(imother);
 249     gMC->CurrentVolOffID(imother++, copyNo);
 250     TString add = "/";
 251     add += name;
 252     add += ".";
 253     add += copyNo;
 254     path.Insert(0,add);
 255   }
 256   while ( name != TString("ALIC") );
 257
 258   return path;
 259 }
 260
 261 //_______________________________________________________________________________
 262 void AliMUONv1::StepManager()
 263 {
 264   // Stepmanager for the chambers
 265   // TBR
 266
 267  if (fStepManagerVersionDE) {
 268     StepManager2();
 269     return;
 270   }
 271
 272
 273   // Only charged tracks
 274   if( !(gMC->TrackCharge()) ) return;
 275   // Only charged tracks
 276
 277   // Only gas gap inside chamber
 278   // Tag chambers and record hits when track enters
 279   static Int_t   idvol=-1;
 280   Int_t   iChamber=0;
 281   Int_t   id=0;
 282   Int_t   copy;
 283   const  Float_t kBig = 1.e10;
 284
 285
 286   //
 287   // Only gas gap inside chamber
 288   // Tag chambers and record hits when track enters
 289   id=gMC->CurrentVolID(copy);
 290   iChamber = GetChamberId(id);
 291   idvol = iChamber -1;
 292
 293   if (idvol == -1) return;
 294
 295   // Filling TrackRefs file for MUON. Our Track references are the active volume of the chambers
 296   if ( (gMC->IsTrackEntering() || gMC->IsTrackExiting() ) )
 297     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 298
 299    if( gMC->IsTrackEntering() ) {
 300      Float_t theta = fTrackMomentum.Theta();
 301      if ((TMath::Pi()-theta)*kRaddeg>=15.) gMC->SetMaxStep(fStepMaxInActiveGas); // We use Pi-theta because z is negative
 302   }
 303
 304    //  AliDebug(1,
 305    //    Form("Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber,
 306    //    ( (AliMUONChamber*)(*fChambers)[idvol])->Z()));
 307    // Particule id and mass,
 308   Int_t     ipart = gMC->TrackPid();
 309   Float_t   mass  = gMC->TrackMass();
 310
 311   fDestepSum[idvol]+=gMC->Edep();
 312   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 313   if ( fStepSum[idvol]==0.0 )  gMC->TrackMomentum(fTrackMomentum);
 314   fStepSum[idvol]+=gMC->TrackStep();
 315
 316   //  if(AliDebugLevel()) {
 317   //  AliDebug(1,
 318   //         Form("iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g",
 319   //              iChamber,ipart, fTrackMomentum.Theta()*kRaddeg, fTrackMomentum.Phi()*kRaddeg,
 320   //              mass, fStepSum[idvol], gMC->Edep()));
 321   //  AliDebug(1,
 322   //         Form("Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(),
 323   //              fTrackMomentum.Z()));
 324   //  gMC->TrackPosition(fTrackPosition);
 325   //  AliDebug(1,
 326   //         Form("Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),
 327   //      fTrackPosition.Z())) ;
 328   //   }
 329
 330   // Track left chamber or StepSum larger than fStepMaxInActiveGas
 331   if ( gMC->IsTrackExiting() ||
 332        gMC->IsTrackStop() ||
 333        gMC->IsTrackDisappeared()||
 334        (fStepSum[idvol]>fStepMaxInActiveGas) ) {
 335
 336     if   ( gMC->IsTrackExiting() ||
 337            gMC->IsTrackStop() ||
 338            gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
 339
 340     gMC->TrackPosition(fTrackPosition);
 341     Float_t theta = fTrackMomentum.Theta();
 342     Float_t phi   = fTrackMomentum.Phi();
 343
 344     TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
 345                                fStepSum[idvol]/2.*sin(theta)*sin(phi),
 346                                fStepSum[idvol]/2.*cos(theta),0.0       );
 347     //     AliDebug(1,Form("Exit: Track Position %f %f %f",fTrackPosition.X(),
 348     //                     fTrackPosition.Y(),fTrackPosition.Z())) ;
 349     //     AliDebug(1,Form("Exit: Track backToWire %f %f %f",backToWire.X(),
 350     //                     backToWire.Y(),backToWire.Z()) ;
 351     fTrackPosition-=backToWire;
 352
 353     //-------------- Angle effect
 354     // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
 355
 356     Float_t betaxGamma    = fTrackMomentum.P()/mass;//  pc/mc2
 357     Float_t sigmaEffect10degrees;
 358     Float_t sigmaEffectThetadegrees;
 359     Float_t eLossParticleELossMip;
 360     Float_t yAngleEffect=0.;
 361     Float_t thetawires      =  TMath::Abs( TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) );// We use Pi-theta because z is negative
 362
 363
 364     if (fAngleEffect){
 365     if ( (betaxGamma >3.2)   &&  (thetawires*kRaddeg<=15.) ) {
 366       betaxGamma=TMath::Log(betaxGamma);
 367       eLossParticleELossMip = fElossRatio->Eval(betaxGamma);
 368       // 10 degrees is a reference for a model (arbitrary)
 369       sigmaEffect10degrees=fAngleEffect10->Eval(eLossParticleELossMip);// in micrometers
 370       // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
 371       sigmaEffectThetadegrees =  sigmaEffect10degrees/fAngleEffectNorma->Eval(thetawires*kRaddeg);  // For 5mm gap
 372       if ( (iChamber==1)  ||  (iChamber==2) )
 373         sigmaEffectThetadegrees/=(1.09833e+00+1.70000e-02*(thetawires*kRaddeg)); // The gap is different (4mm)
 374       yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
 375     }
 376     }
 377
 378     // Detection elements ids
 379     const AliMUONGeometryModule* kGeometryModule
 380       = GetGeometry()->GetModule(iChamber-1);
 381
 382     AliMUONGeometryDetElement* detElement
 383       = kGeometryModule->FindBySensitiveVolume(CurrentVolumePath());
 384
 385     Int_t detElemId = 0;
 386     if (detElement) detElemId = detElement->GetUniqueID();
 387
 388     if (!detElemId) {
 389       cerr << "Chamber id: "
 390            << setw(3) << iChamber << "  "
 391            << "Current SV: "
 392            <<  CurrentVolumePath()
 393            << "  detElemId: "
 394            << setw(5) << detElemId
 395            << endl;
 396       Double_t x, y, z;
 397       gMC->TrackPosition(x, y, z);
 398       cerr << "   global position: "
 399            << x << ", " << y << ", " << z
 400            << endl;
 401       AliWarning("DetElemId not identified.");
 402     }
 403
 404     // One hit per chamber
 405     GetMUONData()->AddHit(fIshunt,
 406                           gAlice->GetMCApp()->GetCurrentTrackNumber(),
 407                           iChamber, ipart,
 408                           fTrackPosition.X(),
 409                           fTrackPosition.Y()+yAngleEffect,
 410                           fTrackPosition.Z(),
 411                           gMC->TrackTime(),
 412                           fTrackMomentum.P(),
 413                           theta,
 414                           phi,
 415                           fStepSum[idvol],
 416                           fDestepSum[idvol],
 417                           fTrackPosition.X(),
 418                           fTrackPosition.Y(),
 419                           fTrackPosition.Z());
 420
 421 //     if (AliDebugLevel()){
 422 //       AliDebug(1,Form("Exit: Particle exiting from chamber %d",iChamber));
 423 //       AliDebug(1,Form("Exit: StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]));
 424 //       AliDebug(1,Form("Exit: Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z())) ;
 425 //     }
 426     fStepSum[idvol]  =0; // Reset for the next event
 427     fDestepSum[idvol]=0; // Reset for the next event
 428   }
 429 }
 430
 431 //_______________________________________________________________________________
 432 void AliMUONv1::StepManager2()
 433 {
 434   // Stepmanager for the chambers
 435
 436   // Only charged tracks
 437   if( !(gMC->TrackCharge()) ) return;
 438   // Only charged tracks
 439
 440   // Only gas gap inside chamber
 441   // Tag chambers and record hits when track enters
 442   static Int_t   idvol=-1, iEnter = 0;
 443   Int_t   iChamber=0;
 444   Int_t   id=0;
 445   Int_t   copy;
 446   const  Float_t kBig = 1.e10;
 447   static Double_t xyzEnter[3];
 448
 449   //
 450   // Only gas gap inside chamber
 451   // Tag chambers and record hits when track enters
 452   id=gMC->CurrentVolID(copy);
 453   iChamber = GetChamberId(id);
 454   idvol = iChamber -1;
 455
 456   if (idvol == -1) return;
 457
 458   // Filling TrackRefs file for MUON. Our Track references are the active volume of the chambers
 459   if ( (gMC->IsTrackEntering() || gMC->IsTrackExiting() ) )
 460     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 461
 462    if( gMC->IsTrackEntering() ) {
 463      Float_t theta = fTrackMomentum.Theta();
 464      if ((TMath::Pi()-theta)*kRaddeg>=15.) gMC->SetMaxStep(fStepMaxInActiveGas); // We use Pi-theta because z is negative
 465      iEnter = 1;
 466      gMC->TrackPosition(xyzEnter[0], xyzEnter[1], xyzEnter[2]); // save coordinates of entrance point
 467   }
 468
 469    //   AliDebug(1,
 470    //       Form("Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber,
 471    //            ( (AliMUONChamber*)(*fChambers)[idvol])->Z())) ;
 472   // Particule id and mass,
 473   Int_t     ipart = gMC->TrackPid();
 474   Float_t   mass  = gMC->TrackMass();
 475
 476   fDestepSum[idvol]+=gMC->Edep();
 477   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 478   if ( fStepSum[idvol]==0.0 )  gMC->TrackMomentum(fTrackMomentum);
 479   fStepSum[idvol]+=gMC->TrackStep();
 480
 481   //  if (AliDebugLevel()) {
 482   //   AliDebug(1,Form("Step, iChamber %d, Particle %d, theta %f phi %f mass %f StepSum %f eloss %g",
 483   //                 iChamber,ipart, fTrackMomentum.Theta()*kRaddeg, fTrackMomentum.Phi()*kRaddeg,
 484   //         mass, fStepSum[idvol], gMC->Edep()));
 485   // AliDebug(1,Form("Step:Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(),
 486   //         fTrackMomentum.Z()));
 487   // gMC->TrackPosition(fTrackPosition);
 488   // AliDebug(1,Form("Step: Track Position %f %f %f",fTrackPosition.X(),
 489   //         fTrackPosition.Y(),fTrackPosition.Z())) ;
 490   //}
 491
 492   // Track left chamber or StepSum larger than fStepMaxInActiveGas
 493   if ( gMC->IsTrackExiting() ||
 494        gMC->IsTrackStop() ||
 495        gMC->IsTrackDisappeared()||
 496        (fStepSum[idvol]>fStepMaxInActiveGas) ) {
 497
 498     if   ( gMC->IsTrackExiting() ||
 499            gMC->IsTrackStop() ||
 500            gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig);
 501     if (fDestepSum[idvol] == 0) {
 502       // AZ - no energy release
 503       fStepSum[idvol] = 0; // Reset for the next event
 504       iEnter = 0;
 505       return;
 506     }
 507
 508     gMC->TrackPosition(fTrackPosition);
 509     Float_t theta = fTrackMomentum.Theta();
 510     Float_t phi   = fTrackMomentum.Phi();
 511
 512     if (gMC->IsTrackExiting() && iEnter != 0) {
 513       // AZ - this code is to avoid artificial hit splitting at the
 514       // "fake" boundary inside the same chamber. It will still produce
 515       // 2 hits but with the same coordinates (at the wire) to allow
 516       // their merging at the digitization level.
 517
 518       // Only for a track going from the entrance to the exit from the volume
 519       // Get local coordinates
 520       Double_t xyz0[3], xyz1[3], tmp[3];
 521       gMC->Gmtod(xyzEnter, xyz0, 1); // local coord. at the entrance
 522
 523       fTrackPosition.Vect().GetXYZ(tmp);
 524       gMC->Gmtod(tmp, xyz1, 1); // local coord. at the exit
 525
 526       Double_t dz = -0.5;
 527       if (xyz1[2] != xyz0[2]) dz = xyz0[2] / (xyz1[2] - xyz0[2]);
 528       tmp[0] = xyz0[0] - (xyz1[0] - xyz0[0]) * dz; // local coord. at the wire
 529       tmp[1] = xyz0[1] - (xyz1[1] - xyz0[1]) * dz;
 530       tmp[2] = xyz0[2] - (xyz1[2] - xyz0[2]) * dz;
 531       gMC->Gdtom(tmp, xyz1, 1); // global coord. at the wire
 532       fTrackPosition.SetXYZT(xyz1[0], xyz1[1], xyz1[2], fTrackPosition.T());
 533     } else {
 534       TLorentzVector backToWire( fStepSum[idvol]/2.*sin(theta)*cos(phi),
 535                                  fStepSum[idvol]/2.*sin(theta)*sin(phi),
 536                                  fStepSum[idvol]/2.*cos(theta),0.0       );
 537       fTrackPosition-=backToWire;
 538       //printf(" %d %d %d %f %d \n", gMC->IsTrackExiting(), gMC->IsTrackStop(), gMC->IsTrackDisappeared(), fStepSum[idvol], iEnter);
 539       //    AliDebug(1,
 540       //     Form("Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()));
 541       // AliDebug(1,
 542       //     Form("Exit: Track backToWire %f %f %f",backToWire.X(),backToWire.Y(),backToWire.Z())) ;
 543     }
 544
 545     //-------------- Angle effect
 546     // Ratio between energy loss of particle and Mip as a function of BetaGamma of particle (Energy/Mass)
 547
 548     Float_t betaxGamma    = fTrackMomentum.P()/mass;//  pc/mc2
 549     Float_t sigmaEffect10degrees;
 550     Float_t sigmaEffectThetadegrees;
 551     Float_t eLossParticleELossMip;
 552     Float_t yAngleEffect=0.;
 553     Float_t thetawires      =  TMath::Abs( TMath::ASin( TMath::Sin(TMath::Pi()-theta) * TMath::Sin(phi) ) );// We use Pi-theta because z is negative
 554
 555
 556     if (fAngleEffect){
 557     if ( (betaxGamma >3.2)   &&  (thetawires*kRaddeg<=15.) ) {
 558       betaxGamma=TMath::Log(betaxGamma);
 559       eLossParticleELossMip = fElossRatio->Eval(betaxGamma);
 560       // 10 degrees is a reference for a model (arbitrary)
 561       sigmaEffect10degrees=fAngleEffect10->Eval(eLossParticleELossMip);// in micrometers
 562       // Angle with respect to the wires assuming that chambers are perpendicular to the z axis.
 563       sigmaEffectThetadegrees =  sigmaEffect10degrees/fAngleEffectNorma->Eval(thetawires*kRaddeg);  // For 5mm gap
 564       if ( (iChamber==1)  ||  (iChamber==2) )
 565         sigmaEffectThetadegrees/=(1.09833e+00+1.70000e-02*(thetawires*kRaddeg)); // The gap is different (4mm)
 566       yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
 567     }
 568     }
 569
 570     // Detection elements ids
 571     const AliMUONGeometryModule* kGeometryModule
 572       = GetGeometry()->GetModule(iChamber-1);
 573
 574     AliMUONGeometryDetElement* detElement
 575       = kGeometryModule->FindBySensitiveVolume(CurrentVolumePath());
 576
 577     Int_t detElemId = 0;
 578     if (detElement) detElemId = detElement->GetUniqueID();
 579
 580     if (!detElemId) {
 581       cerr << "Chamber id: "
 582            << setw(3) << iChamber << "  "
 583            << "Current SV: "
 584            <<  CurrentVolumePath()
 585            << "  detElemId: "
 586            << setw(5) << detElemId
 587            << endl;
 588       Double_t x, y, z;
 589       gMC->TrackPosition(x, y, z);
 590       cerr << "   global position: "
 591            << x << ", " << y << ", " << z
 592            << endl;
 593       AliError("DetElemId not identified.");
 594     }
 595
 596     // One hit per chamber
 597     GetMUONData()->AddHit2(fIshunt,
 598                           gAlice->GetMCApp()->GetCurrentTrackNumber(),
 599                           detElemId, ipart,
 600                           fTrackPosition.X(),
 601                           fTrackPosition.Y()+yAngleEffect,
 602                           fTrackPosition.Z(),
 603                           gMC->TrackTime(),
 604                           fTrackMomentum.P(),
 605                           theta,
 606                           phi,
 607                           fStepSum[idvol],
 608                           fDestepSum[idvol],
 609                           fTrackPosition.X(),
 610                           fTrackPosition.Y(),
 611                           fTrackPosition.Z());
 612
 613     //       AliDebug(1,Form("Exit: Particle exiting from chamber %d",iChamber));
 614     //       AliDebug(1,Form("Exit: StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]));
 615     //       AliDebug(1,Form("Exit: Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
 616
 617     fStepSum[idvol]  =0; // Reset for the next event
 618     fDestepSum[idvol]=0; // Reset for the next event
 619     iEnter = 0;
 620   }
 621 }