MpdEmcGeoUtils.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330
  1. ////////////////////////////////////////////////////////////////
  2. // //
  3. // MpdEmcGeoUtils //
  4. // Geometrical transformations for EMC //
  5. // Author List : D.Peresunko, KI, 2019 //
  6. // //
  7. ////////////////////////////////////////////////////////////////
  8. #include "MpdEmcGeoUtils.h"
  9. #include "FairLogger.h"
  10. #include "TGeoCompositeShape.h"
  11. #include "TGeoArb8.h"
  12. #include "TGeoManager.h"
  13. #include "TGeoMatrix.h"
  14. #include "TGeoTube.h"
  15. #include "TGeoVolume.h"
  16. #include "MpdEmcSimParams.h"
  17. #include <cstdlib>
  18. // these initialisations are needed for a singleton
  19. MpdEmcGeoUtils* MpdEmcGeoUtils::sGeom = nullptr;
  20. ClassImp(MpdEmcGeoUtils)
  21. // ----- Default constructor -------------------------------------------
  22. MpdEmcGeoUtils::MpdEmcGeoUtils()
  23. : fNTowsersPerModule(16),
  24. fNTowersPerCrate(128),
  25. fNTowersPerSector(768),
  26. fNCratesPerSector(6),
  27. fNSectorsPerChamber(25),
  28. fNChambers(2),
  29. fNTowersPerChamber(19200),
  30. fEcalRmin(0.), // minimal ECAL radius
  31. fEcalRmax(0.), // maximal ECAL radius
  32. fEcalZmax(0.) // ECAL half z size
  33. {
  34. }
  35. int MpdEmcGeoUtils::AreNeighbours(int detId1, int detId2) const
  36. {
  37. // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
  38. // = 1 are neighbour
  39. // = 2 are not neighbour but do not continue searching
  40. // neighbours are defined as digits having at least a common side
  41. // The order of detId1 and detId2 is important: first (d1) should be a digit already in a cluster
  42. // which is compared to a digit (d2) not yet in a cluster
  43. int ch1, sector1, iphi1, iz1;
  44. DetIdToRelIndex(detId1, ch1, sector1, iphi1, iz1);
  45. int ch2, sector2, iphi2, iz2;
  46. DetIdToRelIndex(detId2, ch2, sector2, iphi2, iz2);
  47. if (ch1 > ch2){
  48. return -1;
  49. }
  50. if (ch1 < ch2){
  51. return 2;
  52. }
  53. if(MpdEmcSimParams::GetInstance()->AllowMultiSectorClusters()){
  54. if (sector1 == sector2) { // start looking at digits in next sector, can not create inter-sector clusters
  55. int phidiff = std::abs(iphi1 - iphi2);
  56. int zdiff = std::abs(iz1 - iz2);
  57. if (((phidiff <= 1) && (zdiff == 0)) || ((phidiff == 0) && (zdiff <= 1))) { // At least common side
  58. return 1;
  59. } else {
  60. if ((iphi2 > iphi1 + 1)) {
  61. if(iphi1!=0 && iphi1!=MaxPhiSector(sector1) && sector2>sector1){ //not at the edge and current sector alredy scanned
  62. return 2; // Difference in iphi numbers is too large to look further
  63. }
  64. else{ //current digit at the edge, continue looking at other sectors
  65. return 0 ;
  66. }
  67. }
  68. return 0; // continue searching
  69. }
  70. }
  71. //May be edge of sector
  72. if(IsPreviousSector(sector1,sector2)){ // second digit in previous sector
  73. if(iphi1==0 && iphi2==MaxPhiSector(sector2) && iz1==iz2){
  74. return 1;
  75. }
  76. else{
  77. return 0;
  78. }
  79. }
  80. else{
  81. if(IsPreviousSector(sector2,sector1)){ //second digit in next sector
  82. if(iphi1==MaxPhiSector(sector1)){
  83. if(iphi2==0 && iz1==iz2){
  84. return 1;
  85. }
  86. else{
  87. return 0;
  88. }
  89. }
  90. else{
  91. if(iphi1!=0 && sector2>sector1 ){ //not at the edge and current sector scanned
  92. return 2; // Difference in iphi numbers is too large to look further
  93. }
  94. else{ //current digit at the edge, continue looking at other sectors
  95. return 0 ;
  96. }
  97. }
  98. }
  99. else{ //not adjacent sectors
  100. if(iphi1!=0 && iphi1!=MaxPhiSector(sector1) && sector2>sector1 ){
  101. return 2; // Difference in iphi numbers is too large to look further
  102. }
  103. else{ //current digit at the edge, continue looking at other sectors
  104. return 0 ;
  105. }
  106. }
  107. }
  108. }
  109. else{ // only within sector
  110. if (sector1 > sector2) { // start looking at digits in next sector, can not create inter-sector clusters
  111. return -1;
  112. }
  113. if (sector1 < sector2) {
  114. return 2;
  115. }
  116. int phidiff = std::abs(iphi1 - iphi2);
  117. int zdiff = std::abs(iz1 - iz2);
  118. if (((phidiff <= 1) && (zdiff == 0)) || ((phidiff == 0) && (zdiff <= 1))) { // At least common side
  119. return 1;
  120. } else {
  121. if ((iphi2 > iphi1 + 1)) {
  122. return 2; // Difference in iphi numbers is too large to look further
  123. }
  124. return 0; // continue searching
  125. }
  126. }
  127. }
  128. int MpdEmcGeoUtils::AreNeighboursVertex(int detId1, int detId2) const
  129. {
  130. // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
  131. // = 1 are neighbour
  132. // In contrast to previous method, here neighbours are defined as digits having at least a common vertex
  133. // The order of detId1 and detId2 is important: first (d1) should be a digit already in a cluster
  134. // which is compared to a digit (d2) not yet in a cluster
  135. int ch1, sector1, iphi1, iz1;
  136. DetIdToRelIndex(detId1, ch1, sector1, iphi1, iz1);
  137. int ch2, sector2, iphi2, iz2;
  138. DetIdToRelIndex(detId2, ch2, sector2, iphi2, iz2);
  139. if (ch1 != ch2)
  140. return 0;
  141. if(MpdEmcSimParams::GetInstance()->AllowMultiSectorClusters()){
  142. if(sector1 == sector2){
  143. int phidiff = std::abs(iphi1 - iphi2);
  144. int zdiff = std::abs(iz1 - iz2);
  145. if (((phidiff <= 1) && (zdiff <= 1))) { // At least common vertex
  146. return 1;
  147. } else {
  148. return 0; // not a neighbour
  149. }
  150. }
  151. //May be edge of sector
  152. if(IsPreviousSector(sector1,sector2)){ //previous sector
  153. if(iphi1==0 && iphi2==MaxPhiSector(sector2) && TMath::Abs(iz1-iz2)<=1){
  154. return 1;
  155. }
  156. else{
  157. return 0;
  158. }
  159. }
  160. if(IsPreviousSector(sector2,sector1)){
  161. if(iphi1==MaxPhiSector(sector1) && iphi2==0 && TMath::Abs(iz1-iz2)<=1){
  162. return 1;
  163. }
  164. }
  165. return 0;
  166. }
  167. else{ //each sector separately
  168. if (sector1 != sector2) { // start looking at digits in next sector, can not create inter-sector clusters
  169. return 0;
  170. }
  171. int phidiff = std::abs(iphi1 - iphi2);
  172. int zdiff = std::abs(iz1 - iz2);
  173. if (((phidiff <= 1) && (zdiff <= 1))) { // At least common vertex
  174. return 1;
  175. } else {
  176. return 0; // continue searching
  177. }
  178. }
  179. }
  180. int MpdEmcGeoUtils::GeantToDetId(int chamberH, int chamber, int sector, int crate,
  181. int module, int boxA, int boxB) const // Convert Geant volume indexes to abs ID of a channel
  182. {
  183. // Convert Geant volume indexes to abs ID of a channel
  184. int num = boxB +
  185. crate * fNTowersPerCrate +
  186. sector * fNTowersPerSector +
  187. chamber* fNTowersPerChamber ;
  188. return num;
  189. }
  190. void MpdEmcGeoUtils::DetIdToRelIndex(int detId, int& chamber, int& sector, int& iphi, int& iz) const
  191. {
  192. // Convert detId to iphi,iz indexes within one sector
  193. // We use consequent numbering of sectors
  194. //int chamberH = 0;
  195. chamber = detId / fNTowersPerChamber;
  196. detId -= chamber * fNTowersPerChamber;
  197. sector = detId / fNTowersPerSector;
  198. detId-=sector * fNTowersPerSector;
  199. int crate = detId/ fNTowersPerCrate ;
  200. detId-= crate * fNTowersPerCrate ;
  201. iz = detId / 2;
  202. iphi = detId % 2 + 2 * crate;
  203. }
  204. void MpdEmcGeoUtils::DetIdToGlobalIphiIz(int detId, int& iphi, int& iz) const
  205. {
  206. int chamber,sect ;
  207. DetIdToRelIndex(detId,chamber,sect,iphi,iz) ;
  208. iphi=fNCratesPerSector*2*sect+iphi;
  209. iz =chamber*fNTowersPerCrate/2+iz ;
  210. }
  211. void MpdEmcGeoUtils::DetIdToGlobalPosition(int detId, double& x, double& y, double& z) const
  212. {
  213. // calculates senter of front surfase of tower with index detId
  214. if (!gGeoManager) {
  215. LOG(ERROR) << "Can not run without constructed geometry";
  216. x = y = z = 0;
  217. return;
  218. }
  219. int chamber = detId / (fNTowersPerChamber);
  220. detId -= chamber * fNTowersPerChamber;
  221. int sector = detId / fNTowersPerSector;
  222. detId-=sector* fNTowersPerSector ;
  223. int crate = detId / fNTowersPerCrate;
  224. detId-= crate * fNTowersPerCrate ;
  225. int module = detId/fNTowsersPerModule ;
  226. // Construct geant pass from detId:
  227. TString path(Form("cave_1/emcChamber_0/emcChH_%d/emcSector_%d/emcCrate_%d/emcModule%d_0/emc_box%d_%d/",
  228. chamber,sector,crate, module, detId/2 + 1, detId));
  229. if (gGeoManager->CheckPath(path.Data())) {
  230. gGeoManager->cd(path.Data());
  231. TGeoHMatrix* m = gGeoManager->GetCurrentMatrix();
  232. // check sizes of this tower
  233. TGeoVolume* v = gGeoManager->GetVolume(Form("emc_box%d", detId/2 + 1));
  234. double towerHalfSizeZ =0 ;
  235. TGeoCompositeShape* sh = dynamic_cast<TGeoCompositeShape*>(v->GetShape());
  236. if(sh){
  237. towerHalfSizeZ = sh->GetDZ();
  238. }
  239. else{
  240. TGeoArb8 * sh2 = dynamic_cast<TGeoArb8*>(v->GetShape());
  241. towerHalfSizeZ = sh2->GetDZ();
  242. }
  243. towerHalfSizeZ*=0.2; //V
  244. double local[3] = { 0 }, global[3] = { 0 };
  245. if (chamber == 0) { // Why this is different in Chamber0 and Chamber1? Reflection?
  246. local[2] = -towerHalfSizeZ;
  247. } else {
  248. local[2] = towerHalfSizeZ;
  249. }
  250. m->LocalToMaster(local, global);
  251. x = global[0];
  252. y = global[1];
  253. z = global[2];
  254. } else {
  255. LOG(ERROR) << "Can not find volume " << path;
  256. x = y = z = 0;
  257. }
  258. }
  259. int MpdEmcGeoUtils::MaxPhiSector(int sector)const{
  260. //returns number of crates per sector
  261. return fNCratesPerSector*2-1;
  262. }
  263. bool MpdEmcGeoUtils::IsPreviousSector(int sector1, int sector2)const{
  264. // check if sector1 goes before sector 2
  265. if(sector1==0) return sector2==24;
  266. else return sector2==sector1-1;
  267. }
  268. void MpdEmcGeoUtils::GetECALTubeSize(double& rMin, double& rMax, double& zMax)
  269. {
  270. if (fEcalRmin == 0) {
  271. TGeoVolume* v = gGeoManager->GetVolume("emcChH"); //TODO!!!!! check!!!
  272. TGeoTube* sh = dynamic_cast<TGeoTube*>(v->GetShape());
  273. fEcalRmin = sh->GetRmin();
  274. fEcalRmax = sh->GetRmax();
  275. fEcalZmax = sh->GetDz();
  276. }
  277. rMin = fEcalRmin;
  278. rMax = fEcalRmax;
  279. zMax = fEcalZmax;
  280. }
  281. double MpdEmcGeoUtils::Rperp(double z){
  282. //Parameterization of the radius in xy plane ot the surface of EMC clusters
  283. //Parameterization of singe photon simulation from V Riabov (1.12.2019)
  284. //Chi2 = 162512
  285. //NDf = 271
  286. //p0 = 188.186 +/- 6.62242e-05
  287. //p1 = 0.00882988 +/- 2.19833e-06
  288. //p2 = -0.000191223 +/- 1.9986e-08
  289. //p3 = 2.8476e-07 +/- 5.01312e-11
  290. const double p0 = 188.186;
  291. const double p1 = 0.00882988 ;
  292. const double p2 =-0.000191223 ;
  293. const double p3 = 2.8476e-07 ;
  294. return p0 + p1*TMath::Abs(z) + p2*z*z + p3*TMath::Abs(z*z*z) ;
  295. }