// Geometry class for EMCAL : singleton
// EMCAL consists of layers of scintillator and lead
// Places the the Barrel Geometry of The EMCAL at Midrapidity
-// between 0 and 120 degrees of Phi and
+// between 80 and 180(or 190) degrees of Phi and
// -0.7 to 0.7 in eta
// Number of Modules and Layers may be controlled by
// the name of the instance defined
// and : Yves Schutz (SUBATECH)
// and : Jennifer Klay (LBL)
// SHASHLYK : Aleksei Pavlinov (WSU)
+// SuperModules -> module(or tower) -> cell
// --- AliRoot header files ---
+#include <assert.h>
#include <TMath.h>
#include <TVector3.h>
+#include <TRegexp.h>
+#include <TObjArray.h>
+#include <TObjString.h>
+#include <assert.h>
+#include <TMatrixD.h>
+#include <TClonesArray.h>
// -- ALICE Headers.
//#include "AliConst.h"
// --- EMCAL headers
#include "AliEMCALGeometry.h"
+#include "AliEMCALDigit.h"
ClassImp(AliEMCALGeometry)
Bool_t AliEMCALGeometry::fgInit = kFALSE;
TString name; // contains name of geometry
+char *additionalOpts[]={"nl=", // number of sampling layers
+ "pbTh=", // cm, Thickness of the Pb
+ "scTh=" // cm, Thickness of the Sc
+};
+int nAdditionalOpts = sizeof(additionalOpts) / sizeof(char*);
+
//______________________________________________________________________
AliEMCALGeometry::~AliEMCALGeometry(void){
// dtor
//______________________________________________________________________
Bool_t AliEMCALGeometry::AreInSameTower(Int_t id1, Int_t id2) const {
- // Find out whether two hits are in the same tower
+ // Find out whether two hits are in the same tower - have to be change
Int_t idmax = TMath::Max(id1, id2) ;
Int_t idmin = TMath::Min(id1, id2) ;
if ( ((idmax - GetNZ() * GetNPhi()) == idmin ) ||
// New geometry: EMCAL_55_25
// 24-aug-04 for shish-kebab
// SHISH_25 or SHISH_62
+ // 11-oct-05 - correction for pre final design
+ // Feb 06,2006 - decrease the weight of EMCAL
fgInit = kFALSE; // Assume failed until proven otherwise.
name = GetName();
- name.ToUpper();
+ name.ToUpper();
+ fKey110DEG = 0;
+ if(name.Contains("110DEG")) fKey110DEG = 1; // for GetAbsCellId
fNZ = 114; // granularity along Z (eta)
fNPhi = 168; // granularity in phi (azimuth)
fArm1EtaMin = -0.7; // pseudorapidity, Starting EMCAL Eta position
fArm1EtaMax = +0.7; // pseudorapidity, Ending EMCAL Eta position
fIPDistance = 454.0; // cm, Radial distance to inner surface of EMCAL
+ fPhiGapForSM = 0.; // cm, only for final TRD1 geometry
// geometry
- if (name == "EMCAL_55_25") {
- fECPbRadThickness = 0.5; // cm, Thickness of the Pb radiators
- fECScintThick = 0.5; // cm, Thickness of the scintillator
- fNECLayers = 25; // number of scintillator layers
-
- fSampling = 13.1; // calculated with Birk's law implementation
-
- fAlFrontThick = 3.5; // cm, Thickness of front Al layer
- fGap2Active = 1.0; // cm, Gap between Al and 1st Scintillator
- }
- else if( name == "G56_2_55_19" || name == "EMCAL_5655_21" || name == "G56_2_55_19_104_14"|| name == "G65_2_64_19" || name == "EMCAL_6564_21"){
- Fatal("Init", "%s is an old geometry! Please update your Config file", name.Data()) ;
- }
- else if(name.Contains("SHISH")){
- fNumberOfSuperModules = 12; // 12 = 6 * 2 (6 in phi, 2 in Z);
+ if(name.Contains("SHISH")){ // Only shahslyk now
+ // 7-sep-05; integration issue
+ fArm1PhiMin = 80.0; // 60 -> 80
+ fArm1PhiMax = 180.0; // 180 -> 190
+
+ fNumberOfSuperModules = 10; // 12 = 6 * 2 (6 in phi, 2 in Z);
fSteelFrontThick = 2.54; // 9-sep-04
fIPDistance = 460.0;
fFrontSteelStrip = fPassiveScintThick = 0.0; // 13-may-05
// http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
if(name.Contains("TRD1")) { // 30-jan-05
// for final design
- if(name.Contains("MAY05") || name.Contains("WSUC")){
+ fPhiGapForSM = 2.; // cm, only for final TRD1 geometry
+ if(name.Contains("MAY05") || name.Contains("WSUC") || name.Contains("FINAL")){
fNumberOfSuperModules = 12; // 20-may-05
if(name.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05
fNECLayers = 77; // (13-may-05 from V.Petrov)
fPassiveScintThick = 0.8; // 0.8cm = 8mm (13-may-05 from V.Petrov)
fNZ = 24;
fTrd1Angle = 1.5; // 1.3 or 1.5
+
+ if(name.Contains("FINAL")) { // 9-sep-05
+ fNumberOfSuperModules = 10;
+ if(name.Contains("110DEG")) {
+ fNumberOfSuperModules = 12;// last two modules have size 10 degree in phi (180<phi<190)
+ fArm1PhiMax = 200.0; // for XEN1 and turn angle of super modules
+ }
+ fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
+ fEtaModuleSize = fPhiModuleSize;
+ if(name.Contains("HUGE")) fNECLayers *= 3; // 28-oct-05 for analysing leakage
+ }
}
} else if(name.Contains("TRD2")) { // 30-jan-05
fSteelFrontThick = 0.0; // 11-mar-05
fShellThickness = 30.; // should be change
fNPhi = fNZ = 4;
}
+
+ CheckAditionalOptions();
+
// constant for transition absid <--> indexes
fNCellsInTower = fNPHIdiv*fNETAdiv;
fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ;
fNCells = fNCellsInSupMod*fNumberOfSuperModules;
+ if(name.Contains("110DEG")) fNCells -= fNCellsInSupMod;
fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
if(name.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick);
f2Trd2Dy2 = fPhiModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
}
}
- }
- else
- Fatal("Init", "%s is an undefined geometry!", name.Data()) ;
-
+ } else Fatal("Init", "%s is an undefined geometry!", name.Data()) ;
fNPhiSuperModule = fNumberOfSuperModules/2;
if(fNPhiSuperModule<1) fNPhiSuperModule = 1;
if (kTRUE) {
printf("Init: geometry of EMCAL named %s is as follows:\n", name.Data());
- printf( " ECAL : %d x (%f mm Pb, %f mm Sc) \n", GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
+ printf( " ECAL : %d x (%f cm Pb, %f cm Sc) \n", GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ;
if(name.Contains("SHISH")){
printf(" fIPDistance %6.3f cm \n", fIPDistance);
if(fSteelFrontThick>0.)
printf(" fSteelFrontThick %6.3f cm \n", fSteelFrontThick);
printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ);
+ printf(" fNCellsInTower %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInTower, fNCellsInSupMod, fNCells);
if(name.Contains("MAY05")){
printf(" fFrontSteelStrip %6.4f cm (thickness of front steel strip)\n",
fFrontSteelStrip);
printf(" fPassiveScintThick %6.4f cm (thickness of front passive Sc tile)\n",
fPassiveScintThick);
}
- printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
- printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
- printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize);
+ printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize);
+ printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize);
+ printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers);
+ printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize);
printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
}
if(name.Contains("TRD")) {
if(name.Contains("TRD2")) {
printf(" fTrd2AngleY %7.4f\n", fTrd2AngleY);
printf(" f2Trd2Dy2 %7.4f\n", f2Trd2Dy2);
- printf(" fTubsR %7.2f\n", fTubsR);
+ printf(" fTubsR %7.2f cm\n", fTubsR);
printf(" fTubsTurnAngle %7.4f\n", fTubsTurnAngle);
- printf(" fEmptySpace %7.4f\n", fEmptySpace);
+ printf(" fEmptySpace %7.4f cm\n", fEmptySpace);
+ } else if(name.Contains("TRD1") && name.Contains("FINAL")){
+ printf(" fPhiGapForSM %7.4f cm \n", fPhiGapForSM);
+ if(name.Contains("110DEG"))printf(" Last two modules have size 10 degree in phi (180<phi<190)\n");
}
}
printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f\n",
GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
}
+ //TRU parameters. These parameters values are not the final ones.
+ fNTRU = 3 ;
+ fNTRUEta = 3 ;
+ fNTRUPhi = 1 ;
}
+//______________________________________________________________________
+
+void AliEMCALGeometry::CheckAditionalOptions()
+{ // Feb 06,2006
+ fArrayOpts = new TObjArray;
+ Int_t nopt = ParseString(name, *fArrayOpts);
+ if(nopt==1) { // no aditional option(s)
+ fArrayOpts->Delete();
+ delete fArrayOpts;
+ fArrayOpts = 0;
+ return;
+ }
+ for(Int_t i=1; i<nopt; i++){
+ TObjString *o = (TObjString*)fArrayOpts->At(i);
+
+ TString addOpt = o->String();
+ Int_t indj=-1;
+ for(Int_t j=0; j<nAdditionalOpts; j++) {
+ TString opt = additionalOpts[j];
+ if(addOpt.Contains(opt,TString::kIgnoreCase)) {
+ indj = j;
+ break;
+ }
+ }
+ if(indj<0) {
+ printf("<E> option |%s| unavailable : ** look to the file AliEMCALGeometry.h **\n",
+ addOpt.Data());
+ assert(0);
+ } else {
+ printf("<I> option |%s| is valid : number %i : |%s|\n",
+ addOpt.Data(), indj, additionalOpts[indj]);
+ if (addOpt.Contains("NL=",TString::kIgnoreCase)) {// number of sampling layers
+ sscanf(addOpt.Data(),"NL=%i", &fNECLayers);
+ printf(" fNECLayers %i (new) \n", fNECLayers);
+ } else if(addOpt.Contains("PBTH=",TString::kIgnoreCase)) {//Thickness of the Pb
+ sscanf(addOpt.Data(),"PBTH=%f", &fECPbRadThickness);
+ } else if(addOpt.Contains("SCTH=",TString::kIgnoreCase)) {//Thickness of the Sc
+ sscanf(addOpt.Data(),"SCTH=%f", &fECScintThick);
+ }
+ }
+ }
+}
+
+//____________________________________________________________________________
+TClonesArray * AliEMCALGeometry::FillTRU(const TClonesArray * digits) {
+
+
+ //Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
+ //Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
+ //TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta. For the
+ //moment the TRU of the 2 smaller supermodules are considered to be equal
+ //to the rest.
+
+ //Check data members
+
+ if(fNTRUEta*fNTRUPhi != fNTRU)
+ Error("FillTRU"," Wrong number of TRUS per Eta or Phi");
+
+ //Initilize variables
+ //List of TRU matrices initialized to 0.
+ Int_t nCellsPhi = fNPhi*2/fNTRUPhi;
+ Int_t nCellsEta = fNZ*2/fNTRUEta;
+ TClonesArray * matrix = new TClonesArray("TMatrixD",1000);
+
+ for(Int_t k = 0; k < fNTRU*fNumberOfSuperModules; k++){
+ TMatrixD * trus = new TMatrixD(nCellsPhi,nCellsEta) ;
+ for(Int_t i = 0; i < nCellsPhi; i++)
+ for(Int_t j = 0; j < nCellsEta; j++)
+ (*trus)(i,j) = 0.0;
+
+ new((*matrix)[k]) TMatrixD(*trus) ;
+ }
+
+ AliEMCALDigit * dig ;
+
+ //Declare variables
+ Int_t id = -1;
+ Float_t amp = -1;
+ Int_t iSupMod = -1;
+ Int_t nTower = -1;
+ Int_t nIphi = -1;
+ Int_t nIeta = -1;
+ Int_t iphi = -1;
+ Int_t ieta = -1;
+
+ //Digits loop to fill TRU matrices with amplitudes.
+
+ for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
+
+ dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
+ amp = dig->GetAmp() ; //Energy of the digit (arbitrary units)
+ id = dig->GetId() ; //Id label of the cell
+ //cout<<"idig "<<idig<<" Amp "<<amp<<" Id "<<id<<endl;
+
+ //Get eta and phi cell position in supermodule
+ Bool_t bCell = GetCellIndex(id, iSupMod, nTower, nIphi, nIeta) ;
+ if(!bCell)
+ Error("FillTRU","Wrong cell id number") ;
+
+ GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta);
+
+ //Check to which TRU in the supermodule belongs the cell.
+ //Supermodules are divided in a TRU matrix of dimension
+ //(fNTRUPhi,fNTRUEta).
+ //Each TRU is a cell matrix of dimension (nCellsPhi,nCellsEta)
+
+ //First calculate the row and column in the supermodule
+ //of the TRU to which the cell belongs.
+
+ Int_t col = (ieta-1)/nCellsEta+1;
+ Int_t row = (iphi-1)/nCellsPhi+1;
+ Int_t itru = col*row + (iSupMod-1)*fNTRU - 1; //Label number of the TRU
+// Info("FillTRU","SM %d, cell: phi %d, eta %d",iSupMod,iphi,ieta);
+// Info("FillTRU","SM TRU: SMrow %d, SMcol %d, SMtru %d,",row,col,itru);
+
+
+ //Fill TRU matrix with cell values
+
+ TMatrixD * trus = dynamic_cast<TMatrixD *>(matrix->At(itru)) ;
+
+ //Calculate row and column of the cell inside the TRU with number itru
+
+ Int_t irow = (iphi-1) - (row-1) * nCellsPhi;
+ Int_t icol = (ieta-1) - (col-1) * nCellsEta;
+
+ (*trus)(irow,icol) = amp ;
+
+
+ // Info("FillTRU","TRU: row %d, col %d",irow,icol);
+
+ }
+ return matrix;
+}
+
+
//______________________________________________________________________
AliEMCALGeometry * AliEMCALGeometry::GetInstance(){
// Returns the pointer of the unique instance
return rv;
}
+// These methods are obsolete but use in AliEMCALRecPoint - keep it now
//______________________________________________________________________
Int_t AliEMCALGeometry::TowerIndex(Int_t ieta,Int_t iphi) const {
// Returns the tower index number from the based on the Z and Phi
}
return 0;
}
+// ==
//
// == Shish-kebab cases ==
//
-Int_t AliEMCALGeometry::GetAbsCellId(const int nSupMod, const int nTower, const int nIphi, const int nIeta)
-{ // 27-aug-04; corr. 21-sep-04
- static Int_t id; // have to change from 1 to fNCells
- id = fNCellsInSupMod*(nSupMod-1);
+Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta)
+{ // 27-aug-04;
+ // corr. 21-sep-04;
+ // 13-oct-05; 110 degree case
+ // 1 <= nSupMod <= fNumberOfSuperModules
+ // 1 <= nTower <= fNPHI * fNZ ( fNPHI * fNZ/2 for fKey110DEG=1)
+ // 1 <= nIphi <= fNPHIdiv
+ // 1 <= nIeta <= fNETAdiv
+ // 1 <= absid <= fNCells
+ static Int_t id=0; // have to change from 1 to fNCells
+ if(fKey110DEG == 1 && nSupMod > 10) { // 110 degree case; last two supermodules
+ id = fNCellsInSupMod*10 + (fNCellsInSupMod/2)*(nSupMod-11);
+ } else {
+ id = fNCellsInSupMod*(nSupMod-1);
+ }
id += fNCellsInTower *(nTower-1);
id += fNPHIdiv *(nIphi-1);
id += nIeta;
if(id<=0 || id > fNCells) {
- printf(" wrong numerations !!\n");
- printf(" id %6i(will be force to -1)\n", id);
- printf(" fNCells %6i\n", fNCells);
- printf(" nSupMod %6i\n", nSupMod);
- printf(" nTower %6i\n", nTower);
- printf(" nIphi %6i\n", nIphi);
- printf(" nIeta %6i\n", nIeta);
- id = -1;
+// printf(" wrong numerations !!\n");
+// printf(" id %6i(will be force to -1)\n", id);
+// printf(" fNCells %6i\n", fNCells);
+// printf(" nSupMod %6i\n", nSupMod);
+// printf(" nTower %6i\n", nTower);
+// printf(" nIphi %6i\n", nIphi);
+// printf(" nIeta %6i\n", nIeta);
+ id = -TMath::Abs(id);
}
return id;
}
} else return IsInECA(ind);
}
-Bool_t AliEMCALGeometry::GetCellIndex(const Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta)
+Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta)
{ // 21-sep-04
- static Int_t tmp=0;
+ // 19-oct-05;
+ static Int_t tmp=0, sm10=0;
if(absId<=0 || absId>fNCells) {
- Info("GetCellIndex"," wrong abs Id %i !! \n", absId);
+// Info("GetCellIndex"," wrong abs Id %i !! \n", absId);
return kFALSE;
}
- nSupMod = (absId-1) / fNCellsInSupMod + 1;
- tmp = (absId-1) % fNCellsInSupMod;
+ sm10 = fNCellsInSupMod*10;
+ if(fKey110DEG == 1 && absId > sm10) { // 110 degree case; last two supermodules
+ nSupMod = (absId-1-sm10) / (fNCellsInSupMod/2) + 11;
+ tmp = (absId-1-sm10) % (fNCellsInSupMod/2);
+ } else {
+ nSupMod = (absId-1) / fNCellsInSupMod + 1;
+ tmp = (absId-1) % fNCellsInSupMod;
+ }
nTower = tmp / fNCellsInTower + 1;
tmp = tmp % fNCellsInTower;
-
- nIphi = tmp / fNPHIdiv + 1;
- nIeta = tmp % fNPHIdiv + 1;
+ nIphi = tmp / fNPHIdiv + 1;
+ nIeta = tmp % fNPHIdiv + 1;
return kTRUE;
}
-void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(const int nTower, const int nIphi, const int nIeta,
+void AliEMCALGeometry::GetTowerPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, int &iphit, int &ietat)
+{ // added nSupMod; have to check - 19-oct-05 !
+ static Int_t nphi;
+
+ if(fKey110DEG == 1 && nSupMod>=11) nphi = fNPhi/2;
+ else nphi = fNPhi;
+
+ ietat = (nTower-1)/nphi + 1; // have to change from 1 to fNZ
+ iphit = (nTower-1)%nphi + 1; // have to change from 1 to fNPhi
+}
+
+void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta,
int &iphi, int &ieta)
-{ // don't check validity of nTower, nIphi and nIeta index
- // have to change - 1-nov-04 ??
+{ // added nSupMod; Nov 25, 05
static Int_t iphit, ietat;
- ietat = (nTower-1)/fNPhi;
- ieta = ietat*fNETAdiv + nIeta; // change from 1 to fNZ*fNETAdiv
-
- iphit = (nTower-1)%fNPhi;
- iphi = iphit*fNPHIdiv + nIphi; // change from 1 to fNPhi*fNPHIdiv
+ GetTowerPhiEtaIndexInSModule(nSupMod,nTower, iphit, ietat);
+ // have to change from 1 to fNZ*fNETAdiv
+ ieta = (ietat-1)*fNETAdiv + (3-nIeta); // x(module) = -z(SM)
+ // iphi - have to change from 1 to fNPhi*fNPHIdiv
+ iphi = (iphit-1)*fNPHIdiv + nIphi; // y(module) = y(SM)
+}
+// Service routine
+int AliEMCALGeometry::ParseString(const TString &topt, TObjArray &Opt)
+{ // Feb 06, 2006
+ Ssiz_t begin, index, end, end2;
+ begin = index = end = end2 = 0;
+ TRegexp separator("[^ ;,\\t\\s/]+");
+ while ( (begin < topt.Length()) && (index != kNPOS) ) {
+ // loop over given options
+ index = topt.Index(separator,&end,begin);
+ if (index >= 0 && end >= 1) {
+ TString substring(topt(index,end));
+ Opt.Add(new TObjString(substring.Data()));
+ }
+ begin += end+1;
+ }
+ return Opt.GetEntries();
}