X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONv1.cxx;h=d25249e171405fa90f5cc8d0ea6be370f8fad68e;hb=7227b364280147c60cc3e9eccb6339b394e94011;hp=724125a4366b259cdbc78afbced51feb816f3539;hpb=5a0e88a7b299211ae37c9f14aaf5e1404173863e;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONv1.cxx b/MUON/AliMUONv1.cxx index 724125a4366..d25249e1714 100644 --- a/MUON/AliMUONv1.cxx +++ b/MUON/AliMUONv1.cxx @@ -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. * * * @@ -15,41 +15,59 @@ /* $Id$ */ -///////////////////////////////////////////////////////// -// Manager and hits classes for set:MUON version 1 // -///////////////////////////////////////////////////////// - -#include -#include -#include -#include -#include -#include +//----------------------------------------------------------------------------- +// 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 "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 +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "AliMUONVHitStore.h" + +using std::endl; +using std::cout; +using std::setw; +/// \cond CLASSIMP ClassImp(AliMUONv1) +/// \endcond //___________________________________________ AliMUONv1::AliMUONv1() : AliMUON(), - fStepManagerVersionOld(kFALSE), - fStepManagerVersionDE(kFALSE), fAngleEffect(kTRUE), + fMagEffect(kTRUE), fStepMaxInActiveGas(0.6), fStepSum(0x0), fDestepSum(0x0), @@ -57,17 +75,19 @@ AliMUONv1::AliMUONv1() fTrackPosition(), fElossRatio(0x0), fAngleEffect10(0x0), - fAngleEffectNorma(0x0) + fAngleEffectNorma(0x0), + fMagAngleEffectNorma(0x0) { -// Default constructor +/// Default constructor + + AliDebug(1,Form("default (empty) ctor this=%p",this)); } //___________________________________________ -AliMUONv1::AliMUONv1(const char *name, const char *title) - : AliMUON(name,title), - fStepManagerVersionOld(kFALSE), - fStepManagerVersionDE(kFALSE), +AliMUONv1::AliMUONv1(const char *name, const char* title) +: AliMUON(name, title), fAngleEffect(kTRUE), + fMagEffect(kTRUE), fStepMaxInActiveGas(0.6), fStepSum(0x0), fDestepSum(0x0), @@ -75,10 +95,18 @@ AliMUONv1::AliMUONv1(const char *name, const char *title) fTrackPosition(), fElossRatio(0x0), fAngleEffect10(0x0), - fAngleEffectNorma(0x0) + fAngleEffectNorma(0x0), + fMagAngleEffectNorma(0x0) { -// Standard onstructor +/// Standard onstructor + 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()]; @@ -107,48 +135,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 - +/// Destructor + + AliDebug(1,Form("dtor this=%p",this)); delete [] fStepSum; delete [] fDestepSum; delete fElossRatio; delete fAngleEffect10; delete fAngleEffectNorma; + 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(); } @@ -156,84 +167,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() { - AliDebug(1,"Start Init for version 1 - CPC chamber type"); - Int_t i; +/// Initialize geometry - - // - // Initialize geometry - // - fGeometryBuilder->InitGeometry(); - AliDebug(1,"Finished Init for version 1 - CPC chamber type"); - - if (fSegmentationType == 1) { - AliFatal("Old Segmentation no more supported"); - } - - if (fSegmentationType == 2) { - fFactory = new AliMUONFactoryV2("New MUON Factory"); - AliInfo("New Segmentation"); - } - - fFactory->Build(this, "default"); - - // - // Initialize segmentation - // - if (!fSegmentationType) { - AliFatal("No Segmentation Type defined."); - return; - } - - for (i=0; iInit(fSegmentationType);// new segmentation + AliDebug(1,"Start Init for version 1 - CPC chamber type"); + fGeometryBuilder->InitGeometry(); + AliDebug(1,"Finished Init for version 1 - CPC chamber type"); - // trigger circuit - // cp - for (i=0; iInit(i); - - - - + // Build segmentation + // using geometry parametrisation + // + // Build response + // + 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; @@ -244,7 +245,7 @@ TString AliMUONv1::CurrentVolumePath() const gMC->CurrentVolOffID(imother++, copyNo); TString add = "/"; add += name; - add += "."; + add += "_"; add += copyNo; path.Insert(0,add); } @@ -256,14 +257,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; @@ -271,192 +265,70 @@ 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 + if ( fIsMaxStep && (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, @@ -488,22 +360,60 @@ void AliMUONv1::StepManager2() gMC->IsTrackDisappeared()|| (fStepSum[idvol]>fStepMaxInActiveGas) ) { - if ( gMC->IsTrackExiting() || - gMC->IsTrackStop() || - gMC->IsTrackDisappeared() ) gMC->SetMaxStep(kBig); + if ( fIsMaxStep && + ( 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) @@ -513,10 +423,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); @@ -529,55 +442,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; } }