/**************************************************************************
* 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 <TGeoMatrix.h>
-#include <TVirtualMC.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 "AliMUONFactory.h"
+#include "AliMUONResponseFactory.h"
#include "AliMUONHit.h"
-#include "AliMUONTriggerCircuit.h"
-#include "AliMUONVGeometryBuilder.h"
-#include "AliMUONChamberGeometry.h"
-#include "AliMUONGeometryEnvelope.h"
-#include "AliMUONGeometryConstituent.h"
+#include "AliMUONGeometryBuilder.h"
+#include "AliMUONGeometry.h"
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONGeometryModule.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"
+
+using std::endl;
+using std::cout;
+using std::setw;
+/// \cond CLASSIMP
ClassImp(AliMUONv1)
+/// \endcond
//___________________________________________
AliMUONv1::AliMUONv1()
: AliMUON(),
- fTrackMomentum(), fTrackPosition(),fGlobalTransformation(0)
+ fAngleEffect(kTRUE),
+ fMagEffect(kTRUE),
+ fStepMaxInActiveGas(0.6),
+ fStepSum(0x0),
+ fDestepSum(0x0),
+ fTrackMomentum(),
+ fTrackPosition(),
+ fElossRatio(0x0),
+ fAngleEffect10(0x0),
+ fAngleEffectNorma(0x0),
+ fMagAngleEffectNorma(0x0)
{
-// Constructor
- fChambers = 0;
- fStepManagerVersionOld = kFALSE;
- fAngleEffect = kTRUE;
- fStepMaxInActiveGas = 0.6;
- fStepSum = 0x0;
- fDestepSum = 0x0;
- fElossRatio = 0x0;
- fAngleEffect10 = 0x0;
- fAngleEffectNorma= 0x0;
+/// Default constructor
+
+ AliDebug(1,Form("default (empty) ctor this=%p",this));
}
+
//___________________________________________
-AliMUONv1::AliMUONv1(const char *name, const char *title)
- : AliMUON(name,title), fTrackMomentum(), fTrackPosition()
+AliMUONv1::AliMUONv1(const char *name, const char* title)
+: AliMUON(name, title),
+ fAngleEffect(kTRUE),
+ fMagEffect(kTRUE),
+ fStepMaxInActiveGas(0.6),
+ fStepSum(0x0),
+ fDestepSum(0x0),
+ fTrackMomentum(),
+ fTrackPosition(),
+ fElossRatio(0x0),
+ fAngleEffect10(0x0),
+ fAngleEffectNorma(0x0),
+ fMagAngleEffectNorma(0x0)
{
-// Constructor
- // By default include all stations
- AliMUONFactory factory;
- factory.Build(this, title);
+/// Standard onstructor
- fStepManagerVersionOld = kFALSE;
- fAngleEffect = kTRUE;
- fStepMaxInActiveGas = 0.6;
+ 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()];
fDestepSum = new Float_t [AliMUONConstants::NCh()];
fAngleEffectNorma->SetParameter(2,5.151e-02);
fAngleEffectNorma->SetParameter(3,-1.490e-03);
- // Define the global transformation:
- // Transformation from the old ALICE coordinate system to a new one:
- // x->-x, z->-z
- TGeoRotation* rotGlobal
- = new TGeoRotation("rotGlobal", 90., 180., 90., 90., 180., 0.);
- fGlobalTransformation = new TGeoCombiTrans(0., 0., 0., rotGlobal);
-}
-
-//_____________________________________________________________________________
-AliMUONv1::AliMUONv1(const AliMUONv1& right)
- : AliMUON(right)
-{
- // copy constructor (not implemented)
-
- Fatal("AliMUONv1", "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
-
- delete fGlobalTransformation;
+/// 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;
-
- Fatal("operator =", "Assignement operator not provided.");
-
- return *this;
-}
-
//__________________________________________________
void AliMUONv1::CreateGeometry()
{
-//
-// Construct geometry using geometry builders.
-//
-
- for (Int_t i=0; i<fGeometryBuilders->GetEntriesFast(); i++) {
-
- // Get the builder
- AliMUONVGeometryBuilder* builder
- = (AliMUONVGeometryBuilder*)fGeometryBuilders->At(i);
-
- // Create geometry with each builder
- if (builder) {
- builder->CreateGeometry();
- builder->SetTransformations();
- }
- }
+/// Construct geometry using geometry builder
- for (Int_t j=0; j<AliMUONConstants::NCh(); j++) {
-
- AliMUONChamberGeometry* geometry = Chamber(j).GetGeometry();
-
- if (!geometry) continue;
- // Skip chambers with not defined geometry
-
- // Loop over envelopes
- const TObjArray* kEnvelopes = geometry->GetEnvelopes();
- for (Int_t k=0; k<kEnvelopes->GetEntriesFast(); k++) {
-
- // Get envelope
- AliMUONGeometryEnvelope* env = (AliMUONGeometryEnvelope*)kEnvelopes->At(k);
- const TGeoCombiTrans* kEnvTrans = env->GetTransformation();
- const char* only = "ONLY";
- if (env->IsMANY()) only = "MANY";
-
- if (env->IsVirtual() && env->GetConstituents()->GetEntriesFast() == 0 ) {
- // virtual envelope + nof constituents = 0
- // => not allowed;
- // empty virtual envelope has no sense
- Fatal("CreateGeometry", "Virtual envelope must have constituents.");
- return;
- }
-
- if (!env->IsVirtual() && env->GetConstituents()->GetEntriesFast() > 0 ) {
- // non virtual envelope + nof constituents > 0
- // => not allowed;
- // use VMC to place constituents
- Fatal("CreateGeometry", "Non virtual envelope cannot have constituents.");
- return;
- }
-
- if (!env->IsVirtual() && env->GetConstituents()->GetEntriesFast() == 0 ) {
- // non virtual envelope + nof constituents = 0
- // => place envelope in ALICE by composed transformation:
- // Tglobal * Tch * Tenv
-
- // Compound chamber transformation with the envelope one
- TGeoHMatrix total
- = (*fGlobalTransformation) *
- (*geometry->GetTransformation()) *
- (*kEnvTrans);
- PlaceVolume(env->GetName(), geometry->GetMotherVolume(),
- env->GetCopyNo(), total, 0, 0, only);
- }
-
- if (env->IsVirtual() && env->GetConstituents()->GetEntriesFast() > 0 ) {
- // virtual envelope + nof constituents > 0
- // => do not place envelope and place constituents
- // in ALICE by composed transformation:
- // Tglobal * Tch * Tenv * Tconst
-
- for (Int_t l=0; l<env->GetConstituents()->GetEntriesFast(); l++) {
- AliMUONGeometryConstituent* constituent
- = (AliMUONGeometryConstituent*)env->GetConstituents()->At(l);
-
- // Compound chamber transformation with the envelope one + the constituent one
- TGeoHMatrix total
- = (*fGlobalTransformation) *
- (*geometry->GetTransformation()) *
- (*kEnvTrans) *
- (*constituent->GetTransformation());
-
- PlaceVolume(constituent->GetName(), geometry->GetMotherVolume(),
- constituent->GetCopyNo(), total,
- constituent->GetNpar(), constituent->GetParam(), only);
- }
- }
- }
- }
-}
-//__________________________________________________________________
-Int_t AliMUONv1::GetChamberId(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;
+ fGeometryBuilder->CreateGeometry();
}
+
//________________________________________________________________
void AliMUONv1::CreateMaterials()
{
+/// Construct materials using geometry builder
- // *** DEFINITION OF AVAILABLE MUON MATERIALS ***
- //
- // Ar-CO2 gas (80%+20%)
- Float_t ag1[3] = { 39.95,12.01,16. };
- Float_t zg1[3] = { 18.,6.,8. };
- Float_t wg1[3] = { .8,.0667,.13333 };
- Float_t dg1 = .001821;
- //
- // Ar-buthane-freon gas -- trigger chambers
- Float_t atr1[4] = { 39.95,12.01,1.01,19. };
- Float_t ztr1[4] = { 18.,6.,1.,9. };
- Float_t wtr1[4] = { .56,.1262857,.2857143,.028 };
- Float_t dtr1 = .002599;
- //
- // Ar-CO2 gas
- Float_t agas[3] = { 39.95,12.01,16. };
- Float_t zgas[3] = { 18.,6.,8. };
- Float_t wgas[3] = { .74,.086684,.173316 };
- Float_t dgas = .0018327;
- //
- // Ar-Isobutane gas (80%+20%) -- tracking
- Float_t ag[3] = { 39.95,12.01,1.01 };
- Float_t zg[3] = { 18.,6.,1. };
- Float_t wg[3] = { .8,.057,.143 };
- Float_t dg = .0019596;
- //
- // Ar-Isobutane-Forane-SF6 gas (49%+7%+40%+4%) -- trigger
- Float_t atrig[5] = { 39.95,12.01,1.01,19.,32.066 };
- Float_t ztrig[5] = { 18.,6.,1.,9.,16. };
- Float_t wtrig[5] = { .49,1.08,1.5,1.84,0.04 };
- Float_t dtrig = .0031463;
- //
- // bakelite
-
- Float_t abak[3] = {12.01 , 1.01 , 16.};
- Float_t zbak[3] = {6. , 1. , 8.};
- Float_t wbak[3] = {6. , 6. , 1.};
- Float_t dbak = 1.4;
-
- Float_t epsil, stmin, deemax, tmaxfd, stemax;
-
- Int_t iSXFLD = gAlice->Field()->Integ();
- Float_t sXMGMX = gAlice->Field()->Max();
- //
- // --- Define the various materials for GEANT ---
- AliMaterial(9, "ALUMINIUM$", 26.98, 13., 2.7, 8.9, 37.2);
- AliMaterial(10, "ALUMINIUM$", 26.98, 13., 2.7, 8.9, 37.2);
- // Air
- Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
- Float_t zAir[4]={6.,7.,8.,18.};
- Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
- Float_t dAir = 1.20479E-3;
- AliMixture(15, "AIR$ ", aAir, zAir, dAir,4, wAir);
- // AliMaterial(15, "AIR$ ", 14.61, 7.3, .001205, 30423.24, 67500);
- AliMixture(19, "Bakelite$", abak, zbak, dbak, -3, wbak);
- AliMixture(20, "ArC4H10 GAS$", ag, zg, dg, 3, wg);
- AliMixture(21, "TRIG GAS$", atrig, ztrig, dtrig, -5, wtrig);
- AliMixture(22, "ArCO2 80%$", ag1, zg1, dg1, 3, wg1);
- AliMixture(23, "Ar-freon $", atr1, ztr1, dtr1, 4, wtr1);
- AliMixture(24, "ArCO2 GAS$", agas, zgas, dgas, 3, wgas);
- // materials for slat:
- // Sensitive area: gas (already defined)
- // PCB: copper
- // insulating material and frame: vetronite
- // walls: carbon, rohacell, carbon
- Float_t aglass[5]={12.01, 28.09, 16., 10.8, 23.};
- Float_t zglass[5]={ 6., 14., 8., 5., 11.};
- Float_t wglass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01};
- Float_t dglass=1.74;
-
- // rohacell: C9 H13 N1 O2
- Float_t arohac[4] = {12.01, 1.01, 14.010, 16.};
- Float_t zrohac[4] = { 6., 1., 7., 8.};
- Float_t wrohac[4] = { 9., 13., 1., 2.};
- Float_t drohac = 0.03;
-
- AliMaterial(31, "COPPER$", 63.54, 29., 8.96, 1.4, 0.);
- AliMixture(32, "Vetronite$",aglass, zglass, dglass, 5, wglass);
- AliMaterial(33, "Carbon$", 12.01, 6., 2.265, 18.8, 49.9);
- AliMixture(34, "Rohacell$", arohac, zrohac, drohac, -4, wrohac);
-
-
- epsil = .001; // Tracking precision,
- stemax = -1.; // Maximum displacement for multiple scat
- tmaxfd = -20.; // Maximum angle due to field deflection
- deemax = -.3; // Maximum fractional energy loss, DLS
- stmin = -.8;
- //
- // Air
- AliMedium(1, "AIR_CH_US ", 15, 1, iSXFLD, sXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
- //
- // Aluminum
-
- AliMedium(4, "ALU_CH_US ", 9, 0, iSXFLD, sXMGMX, tmaxfd, fMaxStepAlu,
- fMaxDestepAlu, epsil, stmin);
- AliMedium(5, "ALU_CH_US ", 10, 0, iSXFLD, sXMGMX, tmaxfd, fMaxStepAlu,
- fMaxDestepAlu, epsil, stmin);
- //
- // Ar-isoC4H10 gas
-
- AliMedium(6, "AR_CH_US ", 20, 1, iSXFLD, sXMGMX, tmaxfd, fMaxStepGas,
- fMaxDestepGas, epsil, stmin);
-//
- // Ar-Isobuthane-Forane-SF6 gas
-
- AliMedium(7, "GAS_CH_TRIGGER ", 21, 1, iSXFLD, sXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
-
- AliMedium(8, "BAKE_CH_TRIGGER ", 19, 0, iSXFLD, sXMGMX, tmaxfd, fMaxStepAlu,
- fMaxDestepAlu, epsil, stmin);
-
- AliMedium(9, "ARG_CO2 ", 22, 1, iSXFLD, sXMGMX, tmaxfd, fMaxStepGas,
- fMaxDestepAlu, epsil, stmin);
- // tracking media for slats: check the parameters!!
- AliMedium(11, "PCB_COPPER ", 31, 0, iSXFLD, sXMGMX, tmaxfd,
- fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
- AliMedium(12, "VETRONITE ", 32, 0, iSXFLD, sXMGMX, tmaxfd,
- fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
- AliMedium(13, "CARBON ", 33, 0, iSXFLD, sXMGMX, tmaxfd,
- fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
- AliMedium(14, "Rohacell ", 34, 0, iSXFLD, sXMGMX, tmaxfd,
- fMaxStepAlu, fMaxDestepAlu, epsil, stmin);
-
-
-
- //.Materials specific to stations
- // created via builders
-
- for (Int_t i=0; i<fGeometryBuilders->GetEntriesFast(); i++) {
+ fGeometryBuilder->CreateMaterials();
+}
- // Get the builder
- AliMUONVGeometryBuilder* builder
- = (AliMUONVGeometryBuilder*)fGeometryBuilders->At(i);
+//________________________________________________________________
+void AliMUONv1::UpdateInternalGeometry()
+{
+/// Update geometry after applying mis-alignment
- // Create materials with each builder
- if (builder) builder->CreateMaterials();
+ // Load mapping
+ if ( ! AliMpCDB::LoadMpSegmentation() ) {
+ AliFatal("Could not access mapping from OCDB !");
}
+
+ fGeometryBuilder->UpdateInternalGeometry();
}
-//______________________________________________________________________________
-void AliMUONv1::PlaceVolume(const TString& name, const TString& mName,
- Int_t copyNo, const TGeoHMatrix& matrix,
- Int_t npar, Double_t* param, const char* only) const
+//________________________________________________________________
+void AliMUONv1::AddAlignableVolumes() const
{
-// Place the volume specified by name with the given transformation matrix
-// ---
-
- // Do not apply global transformation
- // if mother volume == DDIP
- // (as it is applied on this volume)
- TGeoHMatrix transform(matrix);
- if (mName == TString("DDIP")) {
- transform = (*fGlobalTransformation) * transform;
- // To be changed to (*fGlobalTransformation).inverse()
- // when available in TGeo
- // To make this correct also for a general case when
- // (*fGlobalTransformation) * *fGlobalTransformation) != 1
- }
-
- // Decompose transformation
- const Double_t* xyz = transform.GetTranslation();
- const Double_t* rm = transform.GetRotationMatrix();
-
- //cout << "Got translation: "
- // << xyz[0] << " " << xyz[1] << " " << xyz[2] << endl;
-
- //cout << "Got rotation: "
- // << rm[0] << " " << rm[1] << " " << rm[2] << endl
- // << rm[3] << " " << rm[4] << " " << rm[5] << endl
- // << rm[6] << " " << rm[7] << " " << rm[8] << endl;
-
- // Check for presence of rotation
- // (will be nice to be available in TGeo)
- const Double_t kTolerance = 1e-04;
- Bool_t isRotation = true;
- if (TMath::Abs(rm[0] - 1.) < kTolerance &&
- TMath::Abs(rm[1] - 0.) < kTolerance &&
- TMath::Abs(rm[2] - 0.) < kTolerance &&
- TMath::Abs(rm[3] - 0.) < kTolerance &&
- TMath::Abs(rm[4] - 1.) < kTolerance &&
- TMath::Abs(rm[5] - 0.) < kTolerance &&
- TMath::Abs(rm[6] - 0.) < kTolerance &&
- TMath::Abs(rm[7] - 0.) < kTolerance &&
- TMath::Abs(rm[8] - 1.) < kTolerance) isRotation = false;
-
- Int_t krot = 0;
- if (isRotation) {
- TGeoRotation rot;
- rot.SetMatrix(const_cast<Double_t*>(transform.GetRotationMatrix()));
- Double_t theta1, phi1, theta2, phi2, theta3, phi3;
- rot.GetAngles(theta1, phi1, theta2, phi2, theta3, phi3);
-
- //cout << "angles: "
- // << theta1 << " " << phi1 << " "
- // << theta2 << " " << phi2 << " "
- // << theta3 << " " << phi3 << endl;
-
- AliMatrix(krot, theta1, phi1, theta2, phi2, theta3, phi3);
- }
-
- // Place the volume in ALIC
- if (npar == 0)
- gMC->Gspos(name, copyNo, mName, xyz[0], xyz[1], xyz[2] , krot, only);
- else
- gMC->Gsposp(name, copyNo, mName, xyz[0], xyz[1], xyz[2] , krot, only,
- param, npar);
+/// Construct materials using geometry builder
+
+ GetGeometryTransformer()->AddAlignableVolumes();
+}
-}
//___________________________________________
void AliMUONv1::Init()
-{
- //
- // Initialize Tracking Chambers
- //
-
- if(fDebug) printf("\n%s: Start Init for version 1 - CPC chamber type\n\n",ClassName());
- Int_t i;
- for (i=0; i<AliMUONConstants::NCh(); i++) {
- ( (AliMUONChamber*) (*fChambers)[i])->Init();
- }
+{
+/// Initialize geometry
+
+ AliDebug(1,"Start Init for version 1 - CPC chamber type");
- //
- // Set the chamber (sensitive region) GEANT identifier
- //
- for (Int_t i=0; i<fGeometryBuilders->GetEntriesFast(); i++) {
+ fGeometryBuilder->InitGeometry();
+ AliDebug(1,"Finished Init for version 1 - CPC chamber type");
+
- // Get the builder
- AliMUONVGeometryBuilder* builder
- = (AliMUONVGeometryBuilder*)fGeometryBuilders->At(i);
+ // Build segmentation
+ // using geometry parametrisation
+ //
+ // Build response
+ //
+ AliMUONResponseFactory respFactory("default", fIsTailEffect);
+ respFactory.Build(this);
+
+}
- // Set sesitive volumes with each builder
- if (builder) builder->SetSensitiveVolumes();
- }
+//__________________________________________________________________
+Int_t AliMUONv1::GetGeomModuleId(Int_t volId) const
+{
+/// 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;
+ }
-/*
- //
- // Set the chamber (sensitive region) GEANT identifier
- ((AliMUONChamber*)(*fChambers)[0])->SetGid(gMC->VolId("S01G"));
- ((AliMUONChamber*)(*fChambers)[1])->SetGid(gMC->VolId("S02G"));
-
- ((AliMUONChamber*)(*fChambers)[2])->SetGid(gMC->VolId("S03G"));
- ((AliMUONChamber*)(*fChambers)[3])->SetGid(gMC->VolId("S04G"));
-
- ((AliMUONChamber*)(*fChambers)[4])->SetGid(gMC->VolId("S05G"));
- ((AliMUONChamber*)(*fChambers)[5])->SetGid(gMC->VolId("S06G"));
-
- ((AliMUONChamber*)(*fChambers)[6])->SetGid(gMC->VolId("S07G"));
- ((AliMUONChamber*)(*fChambers)[7])->SetGid(gMC->VolId("S08G"));
-
- ((AliMUONChamber*)(*fChambers)[8])->SetGid(gMC->VolId("S09G"));
- ((AliMUONChamber*)(*fChambers)[9])->SetGid(gMC->VolId("S10G"));
-
- ((AliMUONChamber*)(*fChambers)[10])->SetGid(gMC->VolId("SG1A"));
- ((AliMUONChamber*)(*fChambers)[11])->SetGid(gMC->VolId("SG2A"));
- ((AliMUONChamber*)(*fChambers)[12])->SetGid(gMC->VolId("SG3A"));
- ((AliMUONChamber*)(*fChambers)[13])->SetGid(gMC->VolId("SG4A"));
-*/
- if(fDebug) printf("\n%s: Finished Init for version 1 - CPC chamber type\n",ClassName());
-
- //cp
- if(fDebug) printf("\n%s: Start Init for Trigger Circuits\n",ClassName());
- for (i=0; i<AliMUONConstants::NTriggerCircuit(); i++) {
- ( (AliMUONTriggerCircuit*) (*fTriggerCircuits)[i])->Init(i);
- }
- if(fDebug) printf("%s: Finished Init for Trigger Circuits\n",ClassName());
- //cp
+ return -1;
}
//_______________________________________________________________________________
-void AliMUONv1::StepManager()
+TString AliMUONv1::CurrentVolumePath() const
{
- // Stepmanager for the chambers
-
- if (fStepManagerVersionOld) {
- StepManagerOld();
- return;
+/// Return current volume path
+/// (Could be removed when this function is available via gMC)
+
+ TString path = "";
+ TString name;
+ Int_t copyNo;
+ Int_t imother = 0;
+ do {
+ name = gMC->CurrentVolOffName(imother);
+ gMC->CurrentVolOffID(imother++, copyNo);
+ TString add = "/";
+ add += name;
+ add += "_";
+ add += copyNo;
+ path.Insert(0,add);
}
+ while ( name != TString("ALIC") );
+
+ return path;
+}
+
+//_______________________________________________________________________________
+void AliMUONv1::StepManager()
+{
+/// Step manager for the chambers
// Only charged tracks
if( !(gMC->TrackCharge()) ) return;
// 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);
+ 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());
+
+ 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) {
+ 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;
+ }
+
+ 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
}
-// if (GetDebug()) {
-// Float_t z = ( (AliMUONChamber*)(*fChambers)[idvol])->Z() ;
-// Info("StepManager Step","Active volume found %d chamber %d Z chamber is %f ",idvol,iChamber, z);
-// }
+ // 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();
if ( fStepSum[idvol]==0.0 ) gMC->TrackMomentum(fTrackMomentum);
fStepSum[idvol]+=gMC->TrackStep();
-// if (GetDebug()) {
-// Info("StepManager Step","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());
-// Info("StepManager Step","Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(), fTrackMomentum.Z()) ;
-// gMC->TrackPosition(fTrackPosition);
-// Info("StepManager Step","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
-// }
+ // if (AliDebugLevel()) {
+ // AliDebug(1,Form("Step, 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("Step:Track Momentum %f %f %f", fTrackMomentum.X(), fTrackMomentum.Y(),
+ // fTrackMomentum.Z()));
+ // gMC->TrackPosition(fTrackPosition);
+ // AliDebug(1,Form("Step: Track Position %f %f %f",fTrackPosition.X(),
+ // fTrackPosition.Y(),fTrackPosition.Z())) ;
+ //}
// Track left chamber or StepSum larger than fStepMaxInActiveGas
if ( gMC->IsTrackExiting() ||
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 );
- // if (GetDebug())
- // Info("StepManager Exit","Track Position %f %f %f",fTrackPosition.X(),fTrackPosition.Y(),fTrackPosition.Z()) ;
- // if (GetDebug())
- // Info("StepManager 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)
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);
yAngleEffect=1.e-04*gRandom->Gaus(0,sigmaEffectThetadegrees); // Error due to the angle effect in cm
}
}
-
- // 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 (GetDebug()){
-// Info("StepManager Exit","Particle exiting from chamber %d",iChamber);
-// Info("StepManager Exit","StepSum %f eloss geant %g ",fStepSum[idvol],fDestepSum[idvol]);
-// Info("StepManager 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::StepManagerOld()
-{
- // Old Stepmanager for the chambers
- Int_t copy, id;
- static Int_t idvol =-1;
- static Int_t vol[2];
- Int_t ipart;
- TLorentzVector pos;
- TLorentzVector mom;
- Float_t theta,phi;
- Float_t destep, step;
-
- static Float_t sstep;
- static Float_t eloss, eloss2, xhit, yhit, zhit, tof, tlength;
- const Float_t kBig = 1.e10;
- static Float_t hits[15];
-
- TClonesArray &lhits = *fHits;
-
- //
- //
- // Only charged tracks
- if( !(gMC->TrackCharge()) ) return;
- //
- // Only gas gap inside chamber
- // Tag chambers and record hits when track enters
- id=gMC->CurrentVolID(copy);
- vol[0] = GetChamberId(id);
- idvol = vol[0] -1;
-
- if (idvol == -1) return;
-
- //
- // Get current particle id (ipart), track position (pos) and momentum (mom)
- gMC->TrackPosition(pos);
- gMC->TrackMomentum(mom);
-
- ipart = gMC->TrackPid();
+ 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
+ }
+ }
- //
- // momentum loss and steplength in last step
- destep = gMC->Edep();
- step = gMC->TrackStep();
- // cout<<"------------"<<step<<endl;
- //
- // record hits when track enters ...
- if( gMC->IsTrackEntering()) {
-
- gMC->SetMaxStep(fMaxStepGas);
- Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
- Double_t rt = TMath::Sqrt(tc);
- Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
- Double_t tx = mom[0]/pmom;
- Double_t ty = mom[1]/pmom;
- Double_t tz = mom[2]/pmom;
- Double_t s = ((AliMUONChamber*)(*fChambers)[idvol])
- ->ResponseModel()
- ->Pitch()/tz;
- theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
- phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
- hits[0] = Float_t(ipart); // Geant3 particle type
- hits[1] = pos[0]+s*tx; // X-position for hit
- hits[2] = pos[1]+s*ty; // Y-position for hit
- hits[3] = pos[2]+s*tz; // Z-position for hit
- hits[4] = theta; // theta angle of incidence
- hits[5] = phi; // phi angle of incidence
- hits[8] = 0;//PadHits does not exist anymore (Float_t) fNPadHits; // first padhit
- hits[9] = -1; // last pad hit
- hits[10] = mom[3]; // hit momentum P
- hits[11] = mom[0]; // Px
- hits[12] = mom[1]; // Py
- hits[13] = mom[2]; // Pz
- tof=gMC->TrackTime();
- hits[14] = tof; // Time of flight
- tlength = 0;
- eloss = 0;
- eloss2 = 0;
- sstep=0;
- xhit = pos[0];
- yhit = pos[1];
- zhit = pos[2];
- Chamber(idvol).ChargeCorrelationInit();
- // Only if not trigger chamber
-
-// printf("---------------------------\n");
-// printf(">>>> Y = %f \n",hits[2]);
-// printf("---------------------------\n");
+ 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());
-
-
- // if(idvol < AliMUONConstants::NTrackingCh()) {
-// //
-// // Initialize hit position (cursor) in the segmentation model
-// ((AliMUONChamber*) (*fChambers)[idvol])
-// ->SigGenInit(pos[0], pos[1], pos[2]);
-// } else {
-// //geant3->Gpcxyz();
-// //printf("In the Trigger Chamber #%d\n",idvol-9);
-// }
- }
- eloss2+=destep;
- sstep+=step;
-
- // cout<<sstep<<endl;
+ fHitStore->Add(hit);
- //
- // Calculate the charge induced on a pad (disintegration) in case
- //
- // Mip left chamber ...
- if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
- gMC->SetMaxStep(kBig);
- eloss += destep;
- tlength += step;
-
- Float_t x0,y0,z0;
- Float_t localPos[3];
- Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
- gMC->Gmtod(globalPos,localPos,1);
-
- if(idvol < AliMUONConstants::NTrackingCh()) {
-// tracking chambers
- x0 = 0.5*(xhit+pos[0]);
- y0 = 0.5*(yhit+pos[1]);
- z0 = 0.5*(zhit+pos[2]);
- } else {
-// trigger chambers
- x0 = xhit;
- y0 = yhit;
- z0 = 0.;
- }
-
-
- // if (eloss >0) MakePadHits(x0,y0,z0,eloss,tof,idvol);
-
-
- hits[6] = tlength; // track length
- hits[7] = eloss2; // de/dx energy loss
-
-
- // if (fNPadHits > (Int_t)hits[8]) {
- // hits[8] = hits[8]+1;
- // hits[9] = 0: // PadHits does not exist anymore (Float_t) fNPadHits;
- //}
-//
-// new hit
-
- new(lhits[fNhits++])
- AliMUONHit(fIshunt, gAlice->GetMCApp()->GetCurrentTrackNumber(), vol,hits);
- eloss = 0;
- //
- // Check additional signal generation conditions
- // defined by the segmentation
- // model (boundary crossing conditions)
- // only for tracking chambers
- } else if
- ((idvol < AliMUONConstants::NTrackingCh()) &&
- ((AliMUONChamber*) (*fChambers)[idvol])->SigGenCond(pos[0], pos[1], pos[2]))
- {
- ((AliMUONChamber*) (*fChambers)[idvol])
- ->SigGenInit(pos[0], pos[1], pos[2]);
-
- Float_t localPos[3];
- Float_t globalPos[3] = {pos[0], pos[1], pos[2]};
- gMC->Gmtod(globalPos,localPos,1);
-
- eloss += destep;
-
- // if (eloss > 0 && idvol < AliMUONConstants::NTrackingCh())
- // MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),pos[2],eloss,tof,idvol);
- xhit = pos[0];
- yhit = pos[1];
- zhit = pos[2];
- eloss = 0;
- tlength += step ;
- //
- // nothing special happened, add up energy loss
- } else {
- eloss += destep;
- tlength += step ;
+ fStepSum[idvol] =0; // Reset for the next event
+ fDestepSum[idvol]=0; // Reset for the next event
+ iEnter = 0;
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff