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