CbmSttHitProducerIdeal.cxx 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228
  1. // -------------------------------------------------------------------------
  2. // ----- CbmStsHitProducerIdeal source file -----
  3. // ----- Created 10/01/06 by V. Friese -----
  4. // -------------------------------------------------------------------------
  5. #include "TClonesArray.h"
  6. #include "FairRootManager.h"
  7. #include "CbmSttHitProducerIdeal.h"
  8. #include "CbmSttHit.h"
  9. #include "CbmSttHitInfo.h"
  10. #include "CbmSttPoint.h"
  11. #include "TMath.h"
  12. #include <iostream>
  13. // TODO: read this from a configuration file
  14. #define foldResolution 0 // <===== NO SMEARING
  15. #define radialResolutionPolynomialConstant1 0.0150
  16. #define radialResolutionPolynomialConstant2 0.
  17. #define radialResolutionPolynomialConstant3 0.
  18. //#define longitudinalResolutionPolynomialConstant1 3.
  19. #define longitudinalResolutionPolynomialConstant1 0.0001
  20. #define longitudinalResolutionPolynomialConstant2 0.
  21. #define longitudinalResolutionPolynomialConstant3 0.
  22. // TODO: read this from geant initialization
  23. #define innerStrawDiameter 1.
  24. // ----- Default constructor -------------------------------------------
  25. CbmSttHitProducerIdeal::CbmSttHitProducerIdeal() :
  26. FairTask("Ideal STT Hit Producer")
  27. {
  28. }
  29. // -------------------------------------------------------------------------
  30. // ----- Destructor ----------------------------------------------------
  31. CbmSttHitProducerIdeal::~CbmSttHitProducerIdeal()
  32. {
  33. }
  34. // -------------------------------------------------------------------------
  35. // ----- Public method Init --------------------------------------------
  36. InitStatus CbmSttHitProducerIdeal::Init()
  37. {
  38. // Get RootManager
  39. FairRootManager* ioman = FairRootManager::Instance();
  40. if ( ! ioman )
  41. {
  42. cout << "-E- CbmSttHitProducerIdeal::Init: "
  43. << "RootManager not instantiated!" << endl;
  44. return kFATAL;
  45. }
  46. // Get input array
  47. fPointArray = (TClonesArray*) ioman->GetObject("STTPoint");
  48. if ( ! fPointArray )
  49. {
  50. cout << "-W- CbmSttHitProducerIdeal::Init: "
  51. << "No STTPoint array!" << endl;
  52. return kERROR;
  53. }
  54. // Create and register output array
  55. fHitArray = new TClonesArray("CbmSttHit");
  56. ioman->Register("STTHit", "STT", fHitArray, kTRUE);
  57. // Create and register output array
  58. fHitInfoArray = new TClonesArray("CbmSttHitInfo");
  59. ioman->Register("STTHitInfo", "STT", fHitInfoArray, kTRUE);
  60. cout << "-I- CbmSttHitProducerIdeal: Intialisation successfull" << endl;
  61. return kSUCCESS;
  62. }
  63. // -------------------------------------------------------------------------
  64. // ----- Public method Exec --------------------------------------------
  65. void CbmSttHitProducerIdeal::Exec(Option_t* opt)
  66. {
  67. // Reset output array
  68. if ( ! fHitArray )
  69. Fatal("Exec", "No HitArray");
  70. fHitArray->Clear();
  71. fHitInfoArray->Clear();
  72. // Declare some variables
  73. CbmSttPoint
  74. *point = NULL;
  75. Int_t
  76. detID = 0, // Detector ID
  77. trackID = 0; // Track index
  78. TVector3
  79. pos, dpos; // Position and error vectors
  80. // Loop over SttPoints
  81. Int_t
  82. nPoints = fPointArray->GetEntriesFast();
  83. for (Int_t iPoint = 0; iPoint < nPoints; iPoint++)
  84. {
  85. point = (CbmSttPoint*) fPointArray->At(iPoint);
  86. if ( ! point)
  87. continue;
  88. // Detector ID
  89. detID = point->GetDetectorID();
  90. // MCTrack ID
  91. trackID = point->GetTrackID();
  92. // Determine hit position and isochrone (x,y of wire, measured z position)
  93. TVector3
  94. posInLocal(point->GetXInLocal(), point->GetYInLocal(), point->GetZInLocal()),
  95. posOutLocal(point->GetXOutLocal(), point->GetYOutLocal(), point->GetZOutLocal()),
  96. position(point->GetX(), point->GetY(), point->GetZ());
  97. Double_t
  98. closestDistance,
  99. closestDistanceError;
  100. GetClostestApproachToWire(closestDistance, closestDistanceError,
  101. posInLocal, posOutLocal);
  102. Double_t
  103. zpos = position.Z() + ((posOutLocal.Z() + posInLocal.Z()) / 2.),
  104. zposError;
  105. FoldZPosWithResolution(zpos, zposError,
  106. posInLocal, posOutLocal);
  107. TVector3
  108. wireDirection(point->GetXWireDirection(), point->GetYWireDirection(), point->GetZWireDirection());
  109. // Create new hit
  110. pos.SetXYZ(position.X(), position.Y(), zpos);
  111. dpos.SetXYZ(innerStrawDiameter / 2., innerStrawDiameter / 2., GetLongitudinalResolution(position.Z()));
  112. new ((*fHitArray)[iPoint]) CbmSttHit(detID, pos, dpos, iPoint, 0,
  113. closestDistance, closestDistanceError, wireDirection);
  114. new ((*fHitInfoArray)[iPoint]) CbmSttHitInfo(0, 0, trackID, iPoint,
  115. 0, kFALSE);
  116. } // Loop over MCPoints
  117. // Event summary
  118. cout << "-I- CbmSttHitProducerIdeal: " << nPoints << " SttPoints, "
  119. << nPoints << " Hits created." << endl;
  120. }
  121. // -------------------------------------------------------------------------
  122. // ----- Private method GetRadialResolution-------------------------------
  123. Double_t CbmSttHitProducerIdeal::GetRadialResolution(Double_t radius)
  124. {
  125. return radialResolutionPolynomialConstant1 +
  126. radius * radialResolutionPolynomialConstant2 +
  127. radius * radius * radialResolutionPolynomialConstant3;
  128. }
  129. // -------------------------------------------------------------------------
  130. // ----- Private method GetLongitudinalResolution ------------------------
  131. Double_t CbmSttHitProducerIdeal::GetLongitudinalResolution(Double_t zpos)
  132. {
  133. return longitudinalResolutionPolynomialConstant1 +
  134. zpos * longitudinalResolutionPolynomialConstant2 +
  135. zpos * zpos * longitudinalResolutionPolynomialConstant3;
  136. }
  137. // -------------------------------------------------------------------------
  138. // ----- Private method GetClostestApproachToWire ------------------------
  139. void CbmSttHitProducerIdeal::GetClostestApproachToWire(Double_t &closestDistance,
  140. Double_t &closestDistanceError,
  141. TVector3 localInPos,
  142. TVector3 localOutPos)
  143. {
  144. Double_t
  145. a = (localOutPos.X() - localInPos.X()),
  146. b = (localOutPos.Y() - localInPos.Y()),
  147. c = sqrt(a * a + b * b) / 2.;
  148. // fold with Gaussian for resolution
  149. if (c > innerStrawDiameter / 2.)
  150. {
  151. c = innerStrawDiameter / 2.;
  152. }
  153. closestDistance = sqrt((innerStrawDiameter / 2.) * (innerStrawDiameter / 2.) - c * c);
  154. closestDistanceError = 0.;
  155. if (foldResolution)
  156. {
  157. closestDistanceError = GetRadialResolution(closestDistance);
  158. closestDistance += gRandom->Gaus(0., closestDistanceError);
  159. }
  160. }
  161. // -------------------------------------------------------------------------
  162. void CbmSttHitProducerIdeal::FoldZPosWithResolution(Double_t &zpos, Double_t &zposError,
  163. TVector3 localInPos, TVector3 localOutPos)
  164. {
  165. Double_t
  166. zPosInStrawFrame = (localOutPos.Z() - localInPos.Z()) / 2.;
  167. zposError = gRandom->Gaus(0., GetLongitudinalResolution(zPosInStrawFrame));
  168. zpos += zposError;
  169. // cout << "zpos in prod: " << zpos << endl;
  170. }
  171. ClassImp(CbmSttHitProducerIdeal)