MpdTpcSectorGeo.cxx 9.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257
  1. /// \class MpdTpcSectorGeo
  2. ///
  3. /// MPD TPC sector geometry description
  4. /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
  5. #include "MpdTpcSectorGeo.h"
  6. #include "TpcGas.h"
  7. //#include "TpcGeoPar.h"
  8. //#include "FairGeoNode.h"
  9. //#include "FairRunAna.h"
  10. //#include "FairRuntimeDb.h"
  11. #include <TGeoManager.h>
  12. #include <TGeoPgon.h>
  13. #include <TGeoTube.h>
  14. #include <TMath.h>
  15. #include <TSystem.h>
  16. #include <Riostream.h>
  17. using std::cout;
  18. using std::endl;
  19. MpdTpcSectorGeo* MpdTpcSectorGeo::fgTpcSec = 0x0;
  20. //__________________________________________________________________________
  21. MpdTpcSectorGeo* MpdTpcSectorGeo::Instance()
  22. {
  23. /// Get pointer to the TPC sector singleton object
  24. if (!fgTpcSec){
  25. fgTpcSec = new MpdTpcSectorGeo;
  26. // Automatic deletion
  27. std::atexit(DestroyInstance);
  28. }
  29. return fgTpcSec;
  30. }
  31. //__________________________________________________________________________
  32. void MpdTpcSectorGeo::Init()
  33. {
  34. /// Initialization of sector parameters
  35. //FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  36. //TpcGeoPar *geoPar = (TpcGeoPar*) rtdb->getContainer("TpcGeoPar");
  37. TString volName = "tpc01sv";
  38. TGeoVolume *sv = gGeoManager->GetVolume(volName);
  39. if (!sv) {
  40. Warning("Init","No sensitive volume found !!!");
  41. return;
  42. }
  43. TGeoShape *shape = sv->GetShape();
  44. //TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  45. //cout << sensNodes->GetEntriesFast() << " " << geoPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  46. //FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->FindObject(volName));
  47. //FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->At(0));
  48. //TArrayD* params = sensVol->getParameters();
  49. fDphi = TMath::TwoPi() / fgkNsect;
  50. fPhi0 = -fDphi / 2.;
  51. cout << "\n !!!!! ***** ***** ***** !!!!!" << endl;
  52. //if (sensVol->getShape() == "TUBE") {
  53. if (TString(shape->ClassName()) == "TGeoTube") {
  54. // Old sensitive volume
  55. fNrows[0] = fNrows[1] = 50;
  56. //fNrows = 53;
  57. //Double_t rMin = params->At(0);
  58. //Double_t rMax = params->At(1);
  59. Double_t rMin = ((TGeoTube*)shape)->GetRmin();
  60. Double_t rMax = ((TGeoTube*)shape)->GetRmax();
  61. cout << " ***** TPC sensitive volume: " << sv->GetName() << " " << shape->ClassName() << " " << rMin << " " << rMax << " " << ((TGeoTube*)shape)->GetDZ() << endl;
  62. //fYsec[0] = rMin;
  63. //fYsec[1] = rMax * TMath::Cos(fPhi0);
  64. fYsec[0] = 36.; // for consistency with earlier results
  65. fYsec[1] = fYsec[2] = 96.; // for consistency with earlier results
  66. fPadH[0] = fPadH[1] = (fYsec[1] - fYsec[0]) / fNrows[0];
  67. //fYsec[0] = 40.4;
  68. //fYsec[1] = fYsec[0] + 53.0;
  69. //} else if (sensVol->getShape() == "PGON") {
  70. } else if (TString(shape->ClassName()) == "TGeoPgon") {
  71. // New sensitive volume
  72. //fYsec[0] = params->At(5);
  73. //fYsec[1] = params->At(6);
  74. fYsec[0] = ((TGeoPgon*)shape)->Rmin(0) + 0.4;
  75. //fYsec[2] = ((TGeoPgon*)shape)->Rmax(0);
  76. fPadH[0] = 1.2; // approx. pad height in inner ROC region
  77. fPadH[1] = 1.8; //1.2; // approx. pad height in outer ROC region
  78. fNrows[0] = 27; //33;
  79. fNrows[1] = 26; //33;
  80. Double_t dy = fPadH[0] * fNrows[0] + fPadH[1] * fNrows[1];
  81. fYsec[2] = fYsec[0] + dy;
  82. //Double_t scale = (fYsec[2] - fYsec[0]) / dy;
  83. //fPadH[0] *= scale;
  84. //fPadH[1] *= scale;
  85. fYsec[1] = fYsec[0] + fPadH[0] * fNrows[0];
  86. //fYsec[1] = fYsec[0] + 60;
  87. //fNrows = Int_t ((fYsec[1] - fYsec[0] + 0.1) / 1.2); // 66 lays
  88. //fNrows = Int_t ((fYsec[1] - fYsec[0] + 0.1) / 1.5); // 53 lays
  89. //fNrows = Int_t ((fYsec[1] - fYsec[0] + 0.1) / 1.0);
  90. // Starting phi
  91. fPhi0 = -((TGeoPgon*)shape)->Phi1() * TMath::DegToRad();
  92. Double_t loc[3] = {100, 0, 10}, glob[3] = {0};
  93. gGeoManager->FindNode(loc[0],loc[1],loc[2]);
  94. gGeoManager->LocalToMaster(loc,glob);
  95. fPhi0 += TMath::ATan2 (glob[1],glob[0]); // due to rotation
  96. // Extract sensitive volume Z-coordinates
  97. TGeoVolume *membr = gGeoManager->GetVolume("tpc01mb"); // membrane
  98. fZmin = ((TGeoTube*)membr->GetShape())->GetDZ();
  99. //fZmax = fZmin + 2 * ((TGeoTube*)shape)->GetDZ();
  100. fZmax = 170.;
  101. cout << " ***** TPC sensitive volume: " << sv->GetName() << ", shape: " << shape->ClassName()
  102. << ", inner radius: " << fYsec[0] << ", outer radius: " << fYsec[2]
  103. << ", Zmin, Zmax: " << fZmin << ", " << fZmax << endl;
  104. } else Fatal("MpdTpcSectorGeo::Init()"," !!! Unknown sensitive volume shape !!! ");
  105. //fPadH = (fYsec[1] - fYsec[0]) / fNrows;
  106. cout << " ***** TPC sector params - inner radius: " << fYsec[0] << ", outer radius: " << fYsec[2]
  107. << ", boundary: " << fYsec[1] << endl;
  108. cout << " Number of sectors: " << fgkNsect << ", phi0: " << fPhi0*TMath::RadToDeg() << ", numbers of padrows: "
  109. << fNrows[0] << " " << fNrows[1] << ", pad heights: " << fPadH[0] << " " << fPadH[1] << endl;
  110. cout << " !!!!! ***** ***** ***** !!!!!" << endl;
  111. fPadW[0] = fPadW[1] = 0.5; // pad widths
  112. // Numbers of pads in rows
  113. //Double_t tan = TMath::Tan(fDphi/2), dead = 1.35; // dead area on one side
  114. //Double_t tan = TMath::Tan(fDphi/2), dead = 0.15; // dead area on one side
  115. fNPadsInRows = new Int_t [NofRows()];
  116. //for (Int_t j = 0; j < fNrows[0]; ++j)
  117. // fNPadsInRows[j] = Int_t ((tan * (fYsec[0] + j * fPadH[0]) - dead) / fPadW[0]);
  118. //for (Int_t j = 0; j < fNrows[1]; ++j)
  119. // fNPadsInRows[j+fNrows[0]] = Int_t ((tan * (fYsec[1] + j * fPadH[1]) - dead) / fPadW[1]);
  120. Int_t nPiR[] = {20, 21, 21, 22, 23, 23, 24, 24, 25, 26, 26, 27, 28, 28, 29, 30, 30, 31, 32, 32, 33, 33, 34, 35, 35, 36, 37,
  121. 38, 39, 40, 41, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62};
  122. for (Int_t j = 0; j < NofRows(); ++j)
  123. fNPadsInRows[j] = nPiR[j];
  124. // Gas parameters
  125. std::string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  126. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  127. fGas = new TpcGas(tpcGasFile, 130);
  128. fTimeMax = fZmax / fGas->VDrift();
  129. fNTimeBins = 512; // max number of time bins
  130. //fTimeBin = fTimeMax / fNTimeBins;
  131. fTimeBin = 100; // 100 ns
  132. //fZ2TimeBin = fNTimeBins / fZmax;
  133. //fTimeBinMax = Int_t (fTimeMax / fTimeBin);
  134. fTimeBinMax = fTimeMax / fTimeBin;
  135. fZ2TimeBin = fTimeBinMax / fZmax;
  136. }
  137. //__________________________________________________________________________
  138. Int_t MpdTpcSectorGeo::Global2Local(const Double_t *xyzGlob, Double_t *xyzLoc, Int_t iSec0)
  139. {
  140. /// Transform global coordinates to local (sector)
  141. Double_t safety = 0.01;
  142. xyzLoc[2] = xyzGlob[2];
  143. Int_t iSec = iSec0;
  144. if (iSec0 < 0) {
  145. // Find sector No.
  146. Double_t phGlob = TMath::ATan2 (xyzGlob[1], xyzGlob[0]);
  147. if (phGlob < 0) phGlob += TMath::TwoPi();
  148. iSec = Int_t ((phGlob - fPhi0) / fDphi);
  149. if (iSec == fgkNsect) iSec = 0;
  150. } else iSec = (iSec0 >> kSectorS) & kSectorM;
  151. Double_t phSec = iSec * fDphi;
  152. Double_t cosPh = TMath::Cos(phSec);
  153. Double_t sinPh = TMath::Sin(phSec);
  154. Double_t x = xyzGlob[0] * cosPh + xyzGlob[1] * sinPh;
  155. Double_t y = -xyzGlob[0] * sinPh + xyzGlob[1] * cosPh;
  156. if (x < fYsec[0] + safety || x > fYsec[2] - safety) return -1; // outside sector in Y
  157. xyzLoc[0] = y;
  158. xyzLoc[1] = x - fYsec[0];
  159. //if (x < fYsec[0] + safety || x > fYsec[1] - safety) return -1; // outside sector in Y
  160. Int_t row = fNrows[0];
  161. if (x <= fYsec[1]) row = Int_t (xyzLoc[1] / fPadH[0]);
  162. else row += Int_t ((x - fYsec[1]) / fPadH[1]);
  163. return (iSec |= (row << kPadrowS));
  164. }
  165. //__________________________________________________________________________
  166. Int_t MpdTpcSectorGeo::Global2Local(const TVector3 &xyzGlob, TVector3 &xyzLoc, Int_t iSec0)
  167. {
  168. /// Transform global coordinates to local (sector)
  169. Double_t xyz1[3], xyz2[3];
  170. xyzGlob.GetXYZ(xyz1);
  171. Int_t iSec = Global2Local(xyz1, xyz2, iSec0);
  172. xyzLoc.SetXYZ(xyz2[0], xyz2[1], xyz2[2]);
  173. return iSec;
  174. }
  175. //__________________________________________________________________________
  176. void MpdTpcSectorGeo::Local2Global(Int_t iSec, const Double_t *xyzLoc, Double_t *xyzGlob)
  177. {
  178. /// Transform local coordinates of sector iSec to global
  179. xyzGlob[2] = xyzLoc[2];
  180. Double_t phSec = iSec * fDphi;
  181. Double_t cosPh = TMath::Cos(phSec);
  182. Double_t sinPh = TMath::Sin(phSec);
  183. Double_t x = xyzLoc[1] + fYsec[0];
  184. Double_t y = xyzLoc[0];
  185. xyzGlob[0] = x * cosPh - y * sinPh;
  186. xyzGlob[1] = x * sinPh + y * cosPh;
  187. }
  188. //__________________________________________________________________________
  189. void MpdTpcSectorGeo::Local2Global(Int_t iSec, const TVector3 &xyzLoc, TVector3 &xyzGlob)
  190. {
  191. /// Transform local coordinates of sector iSec to global
  192. Double_t xyz1[3], xyz2[3];
  193. xyzLoc.GetXYZ(xyz1);
  194. Local2Global(iSec, xyz1, xyz2);
  195. xyzGlob.SetXYZ(xyz2[0], xyz2[1], xyz2[2]);
  196. }
  197. //__________________________________________________________________________
  198. TVector2 MpdTpcSectorGeo::LocalPadPosition(Int_t padID)
  199. {
  200. /// Return local pad position for padID
  201. Int_t row = PadRow(padID);
  202. Double_t x = 0.0, y = 0.0;
  203. if (row < fNrows[0]) y = fPadH[0] * (row + 0.5);
  204. else y = fYsec[1] - fYsec[0] + fPadH[1] * (row - fNrows[0] + 0.5);
  205. return TVector2(x,y);
  206. }
  207. //__________________________________________________________________________
  208. void MpdTpcSectorGeo::PadID(Float_t xloc, Float_t yloc, UInt_t &row, UInt_t &pad, Float_t &yNext)
  209. {
  210. /// Compute row, pad and distance to the nearest row (+ or -)
  211. row = fNrows[0];
  212. Double_t lowHeight = fYsec[1] - fYsec[0];
  213. if (yloc <= lowHeight) {
  214. row = Int_t (yloc / fPadH[0]);
  215. yNext = (row + 1) * fPadH[0] - yloc + 0.000001;
  216. if (yNext > fPadH[0] / 2) yNext -= (fPadH[0] + 0.000002);
  217. pad = Int_t (xloc / fPadW[0] + fNPadsInRows[row]);
  218. } else {
  219. Int_t dRow = Int_t ((yloc - lowHeight) / fPadH[1]);
  220. row += dRow;
  221. yNext = (dRow + 1) * fPadH[1] + lowHeight - yloc + 0.000001;
  222. if (yNext > fPadH[1] / 2) yNext -= (fPadH[1] + 0.000002);
  223. pad = Int_t (xloc / fPadW[1] + fNPadsInRows[row]);
  224. }
  225. }
  226. //__________________________________________________________________________