MpdStsHitProducerNew2.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289
  1. // -------------------------------------------------------------------------
  2. // ----- MpdStsHitProducerNew2 source file -----
  3. // ----- 29.04.2012 VK -----
  4. // -------------------------------------------------------------------------
  5. // Independant strip readout for each sector
  6. // Geometry: its_sec_cables.geo ( Three types of sectors consist of one, two or three sensors)
  7. #include "TClonesArray.h"
  8. #include "MpdStsGeoPar.h"
  9. #include "MpdStsHitProducerNew2.h"
  10. #include "MpdStsHit.h"
  11. #include "MpdStsPoint.h"
  12. #include "FairDetector.h"
  13. #include "FairGeoNode.h"
  14. #include "FairRootManager.h"
  15. #include "FairRun.h"
  16. #include "FairRunAna.h"
  17. #include "FairRuntimeDb.h"
  18. #include "TGeoManager.h"
  19. #include "TGeoVolume.h"
  20. #include "TGeoNode.h"
  21. #include "TGeoMatrix.h"
  22. #include "TVector3.h"
  23. // ----- Default constructor -------------------------------------------
  24. MpdStsHitProducerNew2::MpdStsHitProducerNew2(const char* fileGeo) :
  25. FairTask("Ideal STS hit Producer") {
  26. fFileGeo=fileGeo;
  27. eneThr = 0.001; // Energy threshold for STS
  28. }
  29. // -------------------------------------------------------------------------
  30. // ----- Destructor ----------------------------------------------------
  31. MpdStsHitProducerNew2::~MpdStsHitProducerNew2() { }
  32. // -------------------------------------------------------------------------
  33. // ----- Public method Init --------------------------------------------
  34. InitStatus MpdStsHitProducerNew2::Init() {
  35. cout << " INITIALIZATION *********************" << endl;
  36. //FairDetector::Initialize();
  37. //FairRun* sim = FairRun::Instance();
  38. //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
  39. // Get RootManager
  40. FairRootManager* ioman = FairRootManager::Instance();
  41. if ( ! ioman ) {
  42. cout << "-E- MpdStsHitProducerNew2::Init: "
  43. << "RootManager not instantiated!" << endl;
  44. return kFATAL;
  45. }
  46. // Get input array
  47. fPointArray = (TClonesArray*) ioman->GetObject("StsPoint");
  48. if ( ! fPointArray ) {
  49. cout << "-W- MpdStsHitProducerNew2::Init: "
  50. << "No StsPoint array!" << endl;
  51. return kERROR;
  52. }
  53. // Create and register output array
  54. fDigiArray = new TClonesArray("MpdStsHit");
  55. ioman->Register("StsHit","Sts",fDigiArray,kTRUE);
  56. // ioman->Register("StsHit","Sts",fDigiArray,kFALSE);
  57. CreateStructure();
  58. cout << "-I- MpdStsHitProducerNew2: Intialization successfull" << endl;
  59. return kSUCCESS;
  60. }
  61. // -------------------------------------------------------------------------
  62. // ----- Private method SetParContainers -------------------------------
  63. void MpdStsHitProducerNew2::SetParContainers() {
  64. // Get run and runtime database
  65. FairRunAna* run = FairRunAna::Instance();
  66. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  67. FairRuntimeDb* db = run->GetRuntimeDb();
  68. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  69. // Get STS geometry parameter container
  70. db->getContainer("MpdStsGeoPar");
  71. }
  72. // -------------------------------------------------------------------------
  73. // ----- Public method Exec --------------------------------------------
  74. void MpdStsHitProducerNew2::Exec(Option_t* opt) {
  75. cout << "-I- MpdStsHitProducerNew2: Execution *********************" << endl;
  76. fDigiArray->Delete();
  77. static Int_t eventCounter = 0, first = 1;
  78. // Detector parameters
  79. //static Double_t sterCos[2] = {-0.131, 0.131}; // stereo-angles in rad (7.5, -7.5)
  80. // static Double_t sterCos[2] = {0.131, -0.131}; // stereo-angles in rad (7.5, -7.5)
  81. static Double_t sterCos[2] = {0.131, 0.}; // stereo-angles in rad (7.5, 0)
  82. static Double_t sterSin[2] = {0.,0.};
  83. static Double_t sterTan[2] = {0.,0.};
  84. static Double_t par1[100] = {0.}, par2[100] = {0.}, par3[100] = {0.};
  85. const Double_t length_1 = 6.2; // cm, length of sector#1
  86. const Double_t length_2 = 12.4; // cm, length of sensor#2
  87. const Double_t length_3 = 18.6; // cm, length of sensor#3
  88. const Double_t width = 6.2; // cm, height of sensor
  89. const Double_t pitch = 0.00605; // cm, readout pitch
  90. const Int_t nChannels = 1024; // number of readout channels in sensor
  91. Int_t nStrips, nShortStrips;
  92. Double_t xShift;
  93. cout << " StsHitProducerNew2 event " << ++eventCounter << endl;
  94. if (first) {
  95. first = 0;
  96. for (Int_t i = 0; i < 2; ++i) {
  97. sterSin[i] = TMath::Sin (sterCos[i]); // sin
  98. sterCos[i] = TMath::Cos (sterCos[i]); // cos
  99. sterTan[i] = sterSin[i]/sterCos[i]; // tan
  100. }
  101. nStrips = Int_t (nChannels + length_1/pitch*sterTan[0]); // Total number of strips in sector#1 including short ones
  102. nShortStrips = Int_t (length_1/pitch*sterTan[0]); // Number of short strips in sector #1 not crossing readout side (extra ships)
  103. // Define sensitive volumes
  104. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  105. MpdStsGeoPar *geoPar = (MpdStsGeoPar*) rtdb->getContainer("MpdStsGeoPar");
  106. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  107. cout << "Number of sensors: "<< sensNodes->GetEntriesFast() << endl;
  108. Int_t nSensors = sensNodes->GetEntriesFast();
  109. for (Int_t i = 0; i < nSensors; ++i) {
  110. FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->At(i));
  111. TArrayD* params = sensVol->getParameters();
  112. par1[i] = params->At(0); // half length of sensor
  113. par2[i] = params->At(1); // half tickness of sensor
  114. par3[i] = params->At(2); // half width of sensor
  115. cout << " *** STS sensitive volume: " << sensVol->GetName() << " " << params->At(0) << " "
  116. << params->At(1) << " " << params->At(2) << endl;
  117. }
  118. }
  119. // Reset output array
  120. fDigiArray->Clear();
  121. // Declare some variables
  122. MpdStsPoint* point = NULL;
  123. // Loop over MCPoints
  124. Int_t nPoints = fPointArray->GetEntriesFast();
  125. cout << "********* Number of STS points: " << nPoints << endl;
  126. for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
  127. point = (MpdStsPoint*) fPointArray->UncheckedAt(iPoint);
  128. //cout << "# of current STS point: " << iPoint+1 << endl;
  129. Int_t detId = point->GetDetectorID();
  130. Int_t sector = ((detId >>1) & 31);
  131. Int_t sectorType = ((detId >> 14) & 3);
  132. Int_t ladder = ((detId >> 6) & 31);
  133. Int_t layer = ((detId >> 11) & 7);
  134. Int_t side = (detId & 1);
  135. TVector3 pos30;
  136. for (Int_t j = 0; j < 2; ++j) {
  137. Double_t pos[3] = {0};
  138. TVector3 pos3;
  139. if (j == 0) {
  140. // One detector side
  141. point->Position(pos3);
  142. pos30 = pos3;
  143. } else {
  144. // The other side
  145. point->PositionOut(pos3);
  146. TVector3 dPos = pos3 - pos30;
  147. if(dPos.Mag() > 0.) dPos.SetMag(1.); // 1 cm
  148. pos3 -= 0.001 * dPos; // small step back
  149. }
  150. pos3.GetXYZ(pos); // global coordinate
  151. TGeoNode *node = gGeoManager->FindNode(pos[0],pos[1],pos[2]);
  152. if (node == 0x0) {
  153. cout << " !!! Missing node " << pos[0] << " " << pos[1] << " " << pos[2] << endl;
  154. exit(0);
  155. }
  156. //cout << gGeoManager->GetPath() << " DetectorID = " << detId+j << endl;
  157. //cout << " Layer: "<< layer << " Ladder: " << ladder << " Sector: " << sector
  158. // << " Sector type: " << sectorType << " Side: " << side+j << " " <<endl;
  159. Double_t loc[3] = {0.};
  160. gGeoManager->MasterToLocal(pos,loc);
  161. // Two hits per point (for two sides)
  162. MpdStsHit *hit = AddHit(iPoint, detId+j, j);
  163. hit->SetUniqueID(2); // AZ - flag for sector layout
  164. hit->SetTrackID(point->GetTrackID());
  165. Double_t x = loc[0], y = loc[1], z = loc[2]; // local coordinate
  166. Double_t xReadOutSide; // X of projection point along strip to readout edge
  167. Int_t iChannel, iStripSector;
  168. if (j == 0) { // Front side
  169. if (sectorType == 1) xReadOutSide = x + (0.5*length_1 - z) * sterTan[0] + 0.5*width;
  170. if (sectorType == 2) xReadOutSide = x + (0.5*length_2 - z) * sterTan[0] + 0.5*width;
  171. if (sectorType == 3) xReadOutSide = x + (0.5*length_3 - z) * sterTan[0] + 0.5*width;
  172. }
  173. else { // Back side
  174. if (sectorType == 1) xReadOutSide = x + (0.5*length_1 - z) * sterTan[1] + 0.5*width -length_1*sterTan[1];
  175. if (sectorType == 2) xReadOutSide = x + (0.5*length_2 - z) * sterTan[1] + 0.5*width -length_2*sterTan[1];
  176. if (sectorType == 3) xReadOutSide = x + (0.5*length_3 - z) * sterTan[1] + 0.5*width -length_3*sterTan[1];
  177. }
  178. if( (sectorType == 1 && (xReadOutSide < 0. || xReadOutSide > width + length_1 * sterTan[0])) ||
  179. (sectorType == 2 && (xReadOutSide < 0. || xReadOutSide > width + length_2 * sterTan[0])) ||
  180. (sectorType == 3 && (xReadOutSide < 0. || xReadOutSide > width + length_3 * sterTan[0])) ) {
  181. cout << "************************ Point is outside sector: xReadOutSide = " << xReadOutSide << endl;
  182. cout << "Local coordinates of point in sensor: " << x << " " << y << " " << z << " " << endl;
  183. cout << gGeoManager->GetPath() << endl;
  184. xReadOutSide = -pitch;
  185. }
  186. iStripSector = Int_t (xReadOutSide/pitch + 1);
  187. iChannel = iStripSector % nChannels;
  188. Double_t xlocSensor = x * sterCos[j] + z * sterSin[j]; // rotated coordinate
  189. //hit->SetLocalX(xReadOutSide);
  190. hit->SetLocalX(xlocSensor); //AZ
  191. hit->SetStrip(iStripSector);
  192. // hit->SetStrip(iChannel);
  193. //*********************************************************************
  194. // cout << "Global coordinates of point in sensor: " << pos[0] << " " << pos[1] << " " << pos[2] << " " << endl;
  195. // cout << "Local coordinates of point in sensor: " << loc[0] << " " << loc[1] << " " << loc[2] << " " << endl;
  196. // cout << "X local along readout side: " << xReadOutSide << endl;
  197. //cout << "Strip's number in sector: " << iStripSector << " Channel's number in sector: " << iChannel<< endl;
  198. //*********************************************************************
  199. }
  200. // cout << point->GetDetectorID() << " " << lay << " " << point->GetZ() << " " << sector << " " << detId << endl;
  201. // fTrackEnergy[point->GetDetectorID()] += point->GetEnergyLoss();
  202. }
  203. }
  204. // -----------------------------------------------------------------------------
  205. // ----- Public method Create Structure ------------------------------------
  206. void MpdStsHitProducerNew2::CreateStructure() {
  207. }
  208. // -----------------------------------------------------------------------------
  209. // ----- Private method AddDigi --------------------------------------------
  210. MpdStsHit* MpdStsHitProducerNew2::AddHit(Int_t trackID, Int_t detID, Int_t side)
  211. {
  212. // It fills the MpdStsHit category
  213. MpdStsPoint *point = (MpdStsPoint*) fPointArray->UncheckedAt(trackID);
  214. TClonesArray& clref = *fDigiArray;
  215. Int_t size = clref.GetEntriesFast();
  216. TVector3 pos;
  217. if (side == 0) point->Position(pos); // one detector side
  218. else point->PositionOut(pos); // the other side
  219. return new(clref[size]) MpdStsHit(detID, pos, TVector3(0.,0.,0.), point->GetEnergyLoss(), trackID, 0);
  220. }
  221. // -----------------------------------------------------------------------------
  222. ClassImp(MpdStsHitProducerNew2)