+//________________________________________________________________________________________________
+void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
+{
+ // Nov 16, 2006- float to double
+ // version for TRD1 only
+ static TVector3 vglob;
+ GetGlobal(absId, vglob);
+ eta = vglob.Eta();
+ phi = vglob.Phi();
+}
+
+//________________________________________________________________________________________________
+void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
+{
+ // Nov 16,2006 - should be discard in future
+ static TVector3 vglob;
+ GetGlobal(absId, vglob);
+ eta = float(vglob.Eta());
+ phi = float(vglob.Phi());
+}
+
+//________________________________________________________________________________________________
+Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
+{
+ // 0<= nSupMod <=11; phi in rad
+ static int i;
+ if(nSupMod<0 || nSupMod >11) return kFALSE;
+ i = nSupMod/2;
+ phiMin = fPhiBoundariesOfSM[2*i];
+ phiMax = fPhiBoundariesOfSM[2*i+1];
+ return kTRUE;
+}
+
+//________________________________________________________________________________________________
+Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const
+{
+ // 0<= nPhiSec <=4; phi in rad
+ // 0; gap boundaries between 0th&2th | 1th&3th SM
+ // 1; gap boundaries between 2th&4th | 3th&5th SM
+ // 2; gap boundaries between 4th&6th | 5th&7th SM
+ // 3; gap boundaries between 6th&8th | 7th&9th SM
+ // 4; gap boundaries between 8th&10th | 9th&11th SM
+ if(nPhiSec<0 || nPhiSec >4) return kFALSE;
+ phiMin = fPhiBoundariesOfSM[2*nPhiSec+1];
+ phiMax = fPhiBoundariesOfSM[2*nPhiSec+2];
+ return kTRUE;
+}
+
+//________________________________________________________________________________________________
+Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const
+{
+ // Return false if phi belongs a phi cracks between SM
+
+ static Int_t i;
+
+ if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE;
+
+ phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries
+ for(i=0; i<6; i++) {
+ if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) {
+ nSupMod = 2*i;
+ if(eta < 0.0) nSupMod++;
+ AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i));
+ return kTRUE;
+ }
+ }
+ return kFALSE;
+}
+
+//________________________________________________________________________________________________
+Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const
+{
+ // Nov 17,2006
+ // stay here - phi problem as usual
+ static Int_t nSupMod, i, ieta, iphi, etaShift, nphi;
+ static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc;
+ absId = nSupMod = - 1;
+ if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) {
+ // phi index first
+ phi = TVector2::Phi_0_2pi(phi);
+ phiLoc = phi - fPhiCentersOfSM[nSupMod/2];
+ nphi = fPhiCentersOfCells.GetSize();
+ if(nSupMod>=10) {
+ phiLoc = phi - 190.*TMath::DegToRad();
+ nphi /= 2;
+ }
+
+ dmin = TMath::Abs(fPhiCentersOfCells[0]-phiLoc);
+ iphi = 0;
+ for(i=1; i<nphi; i++) {
+ d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
+ if(d < dmin) {
+ dmin = d;
+ iphi = i;
+ }
+ // printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
+ }
+ // odd SM are turned with respect of even SM - reverse indexes
+ AliDebug(2,Form(" iphi %i : dmin %f (phi %f, phiLoc %f ) ", iphi, dmin, phi, phiLoc));
+ // eta index
+ absEta = TMath::Abs(eta);
+ etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
+ dmin = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
+ ieta = 0;
+ for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
+ d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
+ if(d < dmin) {
+ dmin = d;
+ ieta = i;
+ }
+ }
+ AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
+
+ if(eta<0) iphi = (nphi-1) - iphi;
+ absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
+
+ return kTRUE;
+ }
+ return kFALSE;
+}
+
+//________________________________________________________________________________________________
+AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta) const
+{
+ //This method was too long to be
+ //included in the header file - the
+ //rule checker complained about it's
+ //length, so we move it here. It returns the
+ //shishkebabmodule at a given eta index point.
+
+ static AliEMCALShishKebabTrd1Module* trd1=0;
+ if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
+ trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
+ } else trd1 = 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)
+{
+ //Browse the modules
+ if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
+}
+
+//________________________________________________________________________________________________
+Bool_t AliEMCALGeometry::IsFolder() const
+{
+ //Check if fShishKebabTrd1Modules is in folder
+ if(fShishKebabTrd1Modules) return kTRUE;
+ else return kFALSE;
+}
+
+//________________________________________________________________________________________________
+Double_t AliEMCALGeometry::GetPhiCenterOfSM(Int_t nsupmod) const
+{
+ //returns center of supermodule in phi
+ int i = nsupmod/2;
+ 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;
+}