// and : Aleksei Pavlinov (WSU)
//
-#include <assert.h>
+#include <cassert>
// --- Root header files ---
#include <Riostream.h>
#include <TObjString.h>
#include <TVector2.h>
#include <TVector3.h>
-
+#include <TParticle.h>
// -- ALICE Headers.
#include "AliLog.h"
// these initialisations are needed for a singleton
AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
Bool_t AliEMCALGeometry::fgInit = kFALSE;
-Char_t* AliEMCALGeometry::fgDefaultGeometryName = "EMCAL_COMPLETE";
+const Char_t* AliEMCALGeometry::fgDefaultGeometryName = "EMCAL_COMPLETE";
//
// Usage:
// You can create the AliEMCALGeometry object independently from anything.
//
// MC: If you work with MC data you have to get geometry the next way:
// == =============================
-// AliRunLoader *rl = AliRunLoader::GetRunLoader();
+// AliRunLoader *rl = AliRunLoader::Instance();
// AliEMCALGeometry *geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
// TGeoManager::Import("geometry.root");
fShellThickness(0.),fZLength(0.),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),
- fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInModule(0),fNTRUEta(0),
- fNTRUPhi(0), fNCellsInTRUEta(0), fNCellsInTRUPhi(0), fTrd1Angle(0.),f2Trd1Dx2(0.),
+ fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInModule(0),
+ // Trigger staff
+ fNTRUEta(0),
+ fNTRUPhi(0),
+ fNModulesInTRUEta(0),
+ fNModulesInTRUPhi(0),
+ fNEtaSubOfTRU(0),
+ //
+ fTrd1Angle(0.),f2Trd1Dx2(0.),
fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0),fEtaMaxOfTRD1(0),
fCentersOfCellsEtaDir(0), fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),
fEtaCentersOfCells(0),fPhiCentersOfCells(0),fShishKebabTrd1Modules(0),
fShellThickness(0.),fZLength(0.),fNZ(0),fNPhi(0),fSampling(0.),fNumberOfSuperModules(0),
fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.),
fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),
- fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInModule(0),fNTRUEta(0),
- fNTRUPhi(0), fNCellsInTRUEta(0), fNCellsInTRUPhi(0), fTrd1Angle(0.),f2Trd1Dx2(0.),
+ fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInModule(0),
+ // Trigger staff
+ fNTRUEta(0),
+ fNTRUPhi(0),
+ fNModulesInTRUEta(0),
+ fNModulesInTRUPhi(0),
+ fNEtaSubOfTRU(0),
+ //
+ fTrd1Angle(0.),f2Trd1Dx2(0.),
fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0), fEtaMaxOfTRD1(0),
fCentersOfCellsEtaDir(0),fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),
fEtaCentersOfCells(0),fPhiCentersOfCells(0),fShishKebabTrd1Modules(0),
fNCells(geom.fNCells),
fNCellsInSupMod(geom.fNCellsInSupMod),
fNCellsInModule(geom.fNCellsInModule),
+ // Trigger staff
fNTRUEta(geom.fNTRUEta),
fNTRUPhi(geom.fNTRUPhi),
- fNCellsInTRUEta(geom.fNCellsInTRUEta),
- fNCellsInTRUPhi(geom.fNCellsInTRUPhi),
+ fNModulesInTRUEta(geom.fNModulesInTRUEta),
+ fNModulesInTRUPhi(geom.fNModulesInTRUPhi),
+ fNEtaSubOfTRU(geom.fNEtaSubOfTRU),
+ //
fTrd1Angle(geom.fTrd1Angle),
f2Trd1Dx2(geom.f2Trd1Dx2),
fPhiGapForSM(geom.fPhiGapForSM),
if(fGeoName.Contains("WSUC")) fGeoName = "EMCAL_WSUC";
//check that we have a valid geometry name
- if(!(fGeoName.Contains("EMCAL_PDC06") || fGeoName.Contains("EMCAL_COMPLETE") || fGeoName.Contains("EMCAL_WSUC"))) {
+ if(!(fGeoName.Contains("EMCAL_PDC06") || fGeoName.Contains("EMCAL_COMPLETE") || fGeoName.Contains("EMCAL_WSUC") || fGeoName.Contains("EMCAL_1stYear"))) {
Fatal("Init", "%s is an undefined geometry!", fGeoName.Data()) ;
}
fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
fEtaModuleSize = fPhiModuleSize;
+ fZLength = 700.; // Z coverage (cm)
+
+
//needs to be called for each geometry and before setting geometry
//parameters which can depend on the outcome
CheckAdditionalOptions();
CheckAdditionalOptions();
}
+ if(fGeoName.Contains("1stYear")){
+ fNumberOfSuperModules = 2;
+
+ if(fGeoName.Contains("LowerEta")) {
+ fNPhiSuperModule = 1;
+ }
+ else if(fGeoName.Contains("LowerPhi_SideA")){
+ fNPhiSuperModule = 2;
+ fArm1EtaMax=0;
+ }
+ else if(fGeoName.Contains("LowerPhi_SideC")){
+ fNPhiSuperModule = 2;
+ fArm1EtaMin=0;
+ }
+
+ CheckAdditionalOptions();
+ }
+
// constant for transition absid <--> indexes
fNCellsInModule = fNPHIdiv*fNETAdiv;
fNCellsInSupMod = fNCellsInModule*fNPhi*fNZ;
if(GetKey110DEG()) fNCells -= fNCellsInSupMod;
fNPhiSuperModule = fNumberOfSuperModules/2;
- if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
+ if(fNPhiSuperModule < 1) fNPhiSuperModule = 1;
fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05
fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05
fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
if(!fGeoName.Contains("WSUC")) fShellThickness = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
-
- fZLength = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage
- fEnvelop[0] = fIPDistance; // mother volume inner radius
- fEnvelop[1] = fIPDistance + fShellThickness; // mother volume outer r.
- fEnvelop[2] = 1.00001*fZLength; // add some padding for mother volume.
+
+ //These parameters are used to create the mother volume to hold the supermodules
+ //2cm padding added to allow for misalignments - JLK 30-May-2008
+ fEnvelop[0] = fIPDistance - 1.; // mother volume inner radius
+ fEnvelop[1] = fIPDistance + fShellThickness + 1.; // mother volume outer r.
+ fEnvelop[2] = fZLength + 2.; //mother volume length
// Local coordinates
fParSM[0] = GetShellThickness()/2.;
fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
- fParSM[2] = 350./2.;
+ fParSM[2] = fZLength/4.; //divide by 4 to get half-length of SM
// SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006
fPhiBoundariesOfSM.Set(fNumberOfSuperModules);
//called after setting of scintillator and lead layer parameters
DefineSamplingFraction();
- //TRU parameters. These parameters values are not the final ones.
- fNTRUEta = 3 ;
- fNTRUPhi = 1 ;
- fNCellsInTRUEta = 16 ;
- fNCellsInTRUPhi = 24 ;
+ // TRU parameters - Apr 29,08 by PAI.
+ // These parameters values was updated at Nov 05, 2007
+ // As is on Olivier BOURRION (LPSC) ppt preasentation
+ // at ALICE trigger meeting at 13th-14th March
+ fNTRUEta = 1; // was 3
+ fNTRUPhi = 3; // was 1
+ fNModulesInTRUEta = 24; // was 8
+ fNModulesInTRUPhi = 4; // was 12
+ // Jet trigger
+ // 3*6*10 + 2*6*2 = 204 -> matrix (nphi(17), neta(12))
+ fNEtaSubOfTRU = 6;
fgInit = kTRUE;
}
+//___________________________________________________________________
void AliEMCALGeometry::PrintGeometry()
{
// Separate routine is callable from broswer; Nov 7,2006
}
//______________________________________________________________________
-void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, char *tit)
+void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, const char *tit)
{
// Service methods
Int_t nSupMod, nModule, nIphi, nIeta;
}
}
+//__________________________________________________________________
void AliEMCALGeometry::DefineSamplingFraction()
{
// Jun 05,2006
}
//______________________________________________________________________
-void AliEMCALGeometry::GetCellPhiEtaIndexInSModuleFromTRUIndex(const Int_t itru, const Int_t iphitru, const Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
+void AliEMCALGeometry::GetModulePhiEtaIndexInSModuleFromTRUIndex(Int_t itru, Int_t iphitru, Int_t ietatru, Int_t &iphiSM, Int_t &ietaSM) const
{
- // This method transforms the (eta,phi) index of cells in a
+ // This method transforms the (eta,phi) index of module in a
// TRU matrix into Super Module (eta,phi) index.
// Calculate in which row and column where the TRU are
// ordered in the SM
- Int_t col = itru/ fNTRUPhi ;
+ Int_t col = itru/ fNTRUPhi ; // indexes of TRU in SM
Int_t row = itru - col*fNTRUPhi ;
- iphiSM = fNCellsInTRUPhi*row + iphitru ;
- ietaSM = fNCellsInTRUEta*col + ietatru ;
+ iphiSM = fNModulesInTRUPhi*row + iphitru ;
+ ietaSM = fNModulesInTRUEta*col + ietatru ;
+ //printf(" GetModulePhiEtaIndexInSModuleFromTRUIndex : itru %2i iphitru %2i ietatru %2i iphiSM %2i ietaSM %2i \n",
+ // itru, iphitru, ietatru, iphiSM, ietaSM);
}
//______________________________________________________________________
AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize()));
Int_t ind=0; // this is phi index
Int_t ieta=0, nModule=0, iphiTemp;
- Double_t xr, zr, theta, phi, eta, r, x,y;
+ Double_t xr=0., zr=0., theta=0., phi=0., eta=0., r=0., x=0.,y=0.;
TVector3 vglob;
Double_t ytCenterModule=0.0, ytCenterCell=0.0;
return trd1;
}
+//________________________________________________________________________________________________
+Int_t AliEMCALGeometry::GetAbsTRUNumberFromNumberInSm(const Int_t row, const Int_t col, const Int_t sm)
+{ // Nov 6, 2007
+ Int_t itru = row + col*GetNModulesInTRUPhi() + sm*GetNTRU();
+ // printf(" GetAbsTRUNumberFromNumberInSm : row %2i col %2i sm %2i -> itru %2i\n", row, col, sm, itru);
+ return itru;
+}
+
//________________________________________________________________________________________________
void AliEMCALGeometry::Browse(TBrowser* b)
{
return fPhiCentersOfSM[i];
}
+//____________________________________________________________________________
+Bool_t AliEMCALGeometry::Impact(const TParticle * particle) const
+{
+ // Tells if a particle enters EMCAL
+ Bool_t in=kFALSE;
+ Int_t AbsID=0;
+ TVector3 vtx(particle->Vx(),particle->Vy(),particle->Vz());
+ TVector3 vimpact(0,0,0);
+ ImpactOnEmcal(vtx,particle->Theta(),particle->Phi(),AbsID,vimpact);
+ if(AbsID!=0)
+ in=kTRUE;
+ return in;
+}
+//____________________________________________________________________________
+void AliEMCALGeometry::ImpactOnEmcal(TVector3 vtx, Double_t theta, Double_t phi,
+ Int_t & absId, TVector3 & vimpact) const
+{
+ // calculates the impact coordinates on EMCAL (centre of a tower/not on EMCAL surface)
+ // of a neutral particle
+ // emitted in the vertex vtx[3] with direction theta and phi in the ALICE global coordinate system
+
+ TVector3 p(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) ;
+
+ vimpact.SetXYZ(0,0,0);
+ absId=-1;
+ if(phi==0 || theta==0) return;
+
+ TVector3 direction;
+ Double_t factor = (GetIPDistance()-vtx[1])/p[1];
+ direction = vtx + factor*p;
+
+ if (!gGeoManager){
+ AliFatal("Geo manager not initialized\n");
+ }
+ //from particle direction -> tower hitted
+ GetAbsCellIdFromEtaPhi(direction.Eta(),direction.Phi(),absId);
+
+ //tower absID hitted -> tower/module plane (evaluated at the center of the tower)
+ Int_t nSupMod, nModule, nIphi, nIeta;
+ Double_t loc[3],loc2[3],loc3[3];
+ Double_t glob[3]={},glob2[3]={},glob3[3]={};
+
+ if(!RelPosCellInSModule(absId,loc)) return;
+
+ //loc is cell center of tower
+ GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
+
+ //look at 2 neighbours-s cell using nIphi={0,1} and nIeta={0,1}
+ Int_t nIphi2,nIeta2,absId2,absId3;
+ if(nIeta==0) nIeta2=1;
+ else nIeta2=0;
+ absId2=GetAbsCellId(nSupMod,nModule,nIphi,nIeta2);
+ if(nIphi==0) nIphi2=1;
+ else nIphi2=0;
+ absId3=GetAbsCellId(nSupMod,nModule,nIphi2,nIeta);
+
+ //2nd point on emcal cell plane
+ if(!RelPosCellInSModule(absId2,loc2)) return;
+
+ //3rd point on emcal cell plane
+ if(!RelPosCellInSModule(absId3,loc3)) return;
+
+ TString volpath = "ALIC_1/XEN1_1/SMOD_";
+ volpath += (nSupMod+1);
+
+ if(GetKey110DEG() && nSupMod>=10) {
+ volpath = "ALIC_1/XEN1_1/SM10_";
+ volpath += (nSupMod-10+1);
+ }
+ if(!gGeoManager->cd(volpath.Data())){
+ AliFatal(Form("GeoManager cannot find path %s!",volpath.Data()))
+ return;
+ }
+ TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
+ if(m) {
+ m->LocalToMaster(loc, glob);
+ m->LocalToMaster(loc2, glob2);
+ m->LocalToMaster(loc3, glob3);
+ } else {
+ AliFatal("Geo matrixes are not loaded \n") ;
+ }
+
+ //Equation of Plane from glob,glob2,glob3 (Ax+By+Cz+D=0)
+ Double_t A = glob[1]*(glob2[2]-glob3[2]) + glob2[1]*(glob3[2]-glob[2]) + glob3[1]*(glob[2]-glob2[2]);
+ Double_t B = glob[2]*(glob2[0]-glob3[0]) + glob2[2]*(glob3[0]-glob[0]) + glob3[2]*(glob[0]-glob2[0]);
+ Double_t C = glob[0]*(glob2[1]-glob3[1]) + glob2[0]*(glob3[1]-glob[1]) + glob3[0]*(glob[1]-glob2[1]);
+ Double_t D = glob[0]*(glob2[1]*glob3[2]-glob3[1]*glob2[2]) + glob2[0]*(glob3[1]*glob[2]-glob[1]*glob3[2]) + glob3[0]*(glob[1]*glob2[2]-glob2[1]*glob[2]);
+ D=-D;
+
+ //shift equation of plane from tower/module center to surface along vector (A,B,C) normal to tower/module plane
+ Double_t dist = GetLongModuleSize()/2.;
+ Double_t norm = TMath::Sqrt(A*A+B*B+C*C);
+ Double_t glob4[3]={};
+ TVector3 dir(A,B,C);
+ TVector3 point(glob[0],glob[1],glob[2]);
+ if(point.Dot(dir)<0) dist*=-1;
+ glob4[0]=glob[0]-dist*A/norm;
+ glob4[1]=glob[1]-dist*B/norm;
+ glob4[2]=glob[2]-dist*C/norm;
+ D = glob4[0]*A + glob4[1]*B + glob4[2]*C ;
+ D = -D;
+
+ //Line determination (2 points for equation of line : vtx and direction)
+ //impact between line (particle) and plane (module/tower plane)
+ Double_t den = A*(vtx(0)-direction(0)) + B*(vtx(1)-direction(1)) + C*(vtx(2)-direction(2));
+ if(den==0){
+ printf("ImpactOnEmcal() No solution :\n");
+ return;
+ }
+
+ Double_t length = A*vtx(0)+B*vtx(1)+C*vtx(2)+D;
+ length /=den;
+
+ vimpact.SetXYZ(vtx(0)+length*(direction(0)-vtx(0)),vtx(1)+length*(direction(1)-vtx(1)),vtx(2)+length*(direction(2)-vtx(2)));
+
+ //shift vimpact from tower/module surface to center along vector (A,B,C) normal to tower/module plane
+ vimpact.SetXYZ(vimpact(0)+dist*A/norm,vimpact(1)+dist*B/norm,vimpact(2)+dist*C/norm);
+
+ return;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff