MpdStsHitProducerV1.cxx 7.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222
  1. /////////////////////////////////////////////////////////////
  2. //
  3. // MpdStsHitProducerV1
  4. //
  5. // Filler of MpdStsHit
  6. //
  7. //
  8. //
  9. ///////////////////////////////////////////////////////////////
  10. #include "TClonesArray.h"
  11. #include "MpdStsGeoPar.h"
  12. #include "MpdStsHitProducerV1.h"
  13. #include "MpdStsHit.h"
  14. #include "MpdStsPoint.h"
  15. #include "FairDetector.h"
  16. #include "FairGeoNode.h"
  17. #include "FairRootManager.h"
  18. #include "FairRun.h"
  19. #include "FairRunAna.h"
  20. #include "FairRuntimeDb.h"
  21. #include "TGeoManager.h"
  22. #include "TGeoVolume.h"
  23. #include "TGeoNode.h"
  24. #include "TGeoMatrix.h"
  25. #include "TVector3.h"
  26. // ----- Default constructor -------------------------------------------
  27. MpdStsHitProducerV1::MpdStsHitProducerV1(const char* fileGeo) :
  28. FairTask("Ideal STS hit Producer") {
  29. fFileGeo=fileGeo;
  30. eneThr = 0.001; // Energy threshold for STS
  31. }
  32. // -------------------------------------------------------------------------
  33. // ----- Destructor ----------------------------------------------------
  34. MpdStsHitProducerV1::~MpdStsHitProducerV1() { }
  35. // -------------------------------------------------------------------------
  36. // ----- Public method Init --------------------------------------------
  37. InitStatus MpdStsHitProducerV1::Init() {
  38. cout << " INITIALIZATION *********************" << endl;
  39. //FairDetector::Initialize();
  40. //FairRun* sim = FairRun::Instance();
  41. //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
  42. // Get RootManager
  43. FairRootManager* ioman = FairRootManager::Instance();
  44. if ( ! ioman ) {
  45. cout << "-E- MpdStsHitProducerV1::Init: "
  46. << "RootManager not instantiated!" << endl;
  47. return kFATAL;
  48. }
  49. // Get input array
  50. fPointArray = (TClonesArray*) ioman->GetObject("StsPoint");
  51. if ( ! fPointArray ) {
  52. cout << "-W- MpdStsHitProducerV1::Init: "
  53. << "No StsPoint array!" << endl;
  54. return kERROR;
  55. }
  56. // Create and register output array
  57. fDigiArray = new TClonesArray("MpdStsHit");
  58. ioman->Register("StsHit","Sts",fDigiArray,kTRUE);
  59. //ioman->Register("StsHit","Sts",fDigiArray,kFALSE);
  60. CreateStructure();
  61. cout << "-I- MpdStsHitProducerV1: Intialization successfull" << endl;
  62. return kSUCCESS;
  63. }
  64. // -------------------------------------------------------------------------
  65. // ----- Private method SetParContainers -------------------------------
  66. void MpdStsHitProducerV1::SetParContainers() {
  67. // Get run and runtime database
  68. FairRunAna* run = FairRunAna::Instance();
  69. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  70. FairRuntimeDb* db = run->GetRuntimeDb();
  71. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  72. // Get STS geometry parameter container
  73. db->getContainer("MpdStsGeoPar");
  74. }
  75. // ----- Public method Exec --------------------------------------------
  76. void MpdStsHitProducerV1::Exec(Option_t* opt) {
  77. //cout << " DIGI EXECUTION *********************" << endl;
  78. static Int_t eventCounter = 0, first = 1;
  79. // Detector parameters
  80. //static Double_t sterCos[2] = {-0.0075, 0.0275}, sterSin[2] = {0}; // stereo-angles in rad
  81. static Double_t sterCos[2] = {-0.131, 0.131}, sterSin[2] = {0}; // stereo-angles in rad
  82. //static Double_t sterCos[2] = {-0., 0.}, sterSin[2] = {0}; // stereo-angles in rad
  83. const Double_t length = 7.3; // cm // length of sensor
  84. const Double_t height = 4.0; // cm // height of sensor
  85. const Double_t pitch = 0.0095; // cm //
  86. const Int_t nstrip = 768; // number of strips in sensor
  87. if (first) {
  88. first = 0;
  89. for (Int_t i = 0; i < 2; ++i) {
  90. sterSin[i] = TMath::Sin (sterCos[i]); // sin
  91. sterCos[i] = TMath::Cos (sterCos[i]); // cos
  92. }
  93. }
  94. cout << " StsHitProducerV1 event " << ++eventCounter << endl;
  95. // Reset output array
  96. fDigiArray->Clear();
  97. // Declare some variables
  98. MpdStsPoint* point = NULL;
  99. map<Int_t, Float_t> fTrackEnergy;
  100. fTrackEnergy.clear();
  101. map<Int_t, Float_t>::const_iterator p;
  102. // Loop over MCPoints
  103. Int_t nPoints = fPointArray->GetEntriesFast();
  104. cout << nPoints << endl;
  105. for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
  106. point = (MpdStsPoint*) fPointArray->At(iPoint);
  107. Int_t detId = point->GetDetectorID();
  108. Int_t lay = ((detId >>1) & 15);
  109. Int_t lad = ((detId >> 5) & 31);
  110. Int_t det = (detId >> 10);
  111. TVector3 pos30;
  112. for (Int_t j = 0; j < 2; ++j) {
  113. Double_t pos[3] = {0};
  114. TVector3 pos3;
  115. if (j == 0) {
  116. // One detector side
  117. point->Position(pos3);
  118. pos30 = pos3;
  119. } else {
  120. // The other side
  121. point->PositionOut(pos3);
  122. TVector3 dPos = pos3 - pos30;
  123. Double_t mag = dPos.Mag();
  124. if (mag > 1.e-7) dPos *= (1. / dPos.Mag());
  125. pos3 -= 0.001 * dPos; // small step back
  126. }
  127. pos3.GetXYZ(pos); // global coordinate
  128. TGeoNode *node = gGeoManager->FindNode(pos[0],pos[1],pos[2]);
  129. if (node == 0x0) {
  130. cout << " !!! Missing node " << pos[0] << " " << pos[1] << " " << pos[2] << endl;
  131. exit(0);
  132. }
  133. cout << gGeoManager->GetPath() << " " << detId+j << endl;
  134. Double_t loc[3] = {0.};
  135. gGeoManager->MasterToLocal(pos,loc);
  136. // Two hits per point (for two sides)
  137. MpdStsHit *hit = AddHit(iPoint, detId+j, j);
  138. Double_t x = loc[0], y = loc[1], z = loc[2]; // local coordinate
  139. Double_t xloc = x * sterCos[j] + z * sterSin[j]; // rotated coordinate
  140. Double_t xpitch = xloc / pitch;
  141. Int_t strip = Int_t (xpitch + (nstrip/2) + 1);
  142. if(strip > nstrip) strip = nstrip;
  143. if(strip < 1) strip = 1;
  144. hit->SetLocalX(xloc);
  145. hit->SetStrip(strip);
  146. //*********************************************************************
  147. cout << pos[0] << " " << pos[1] << " " << pos[2] << " " << endl;
  148. cout << loc[0] << " " << loc[1] << " " << loc[2] << " " << endl;
  149. cout << "xloc=" << x << " xrot=" << xloc << " xpitch=" << xpitch
  150. << " strip=" << strip << endl;
  151. //*********************************************************************
  152. }
  153. //cout << point->GetDetectorID() << " " << lay << " " << point->GetZ() << " " << module << " " << detId << endl;
  154. //fTrackEnergy[point->GetDetectorID()] += point->GetEnergyLoss();
  155. }
  156. }
  157. // -------------------------------------------------------------------------
  158. // ----- Public method Create Structure --------------------------------
  159. void MpdStsHitProducerV1::CreateStructure() {
  160. }
  161. // ----- Private method AddDigi --------------------------------------------
  162. MpdStsHit* MpdStsHitProducerV1::AddHit(Int_t trackID, Int_t detID, Int_t side)
  163. {
  164. // It fills the MpdStsHit category
  165. // cout << "MpdStsHitProducer: track " << trackID << " evt " << eventID << " sec " << sec << " plane " << pla << " strip " << strip << "box " << box << " tube " << tub << endl;
  166. MpdStsPoint *point = (MpdStsPoint*) fPointArray->UncheckedAt(trackID);
  167. TClonesArray& clref = *fDigiArray;
  168. Int_t size = clref.GetEntriesFast();
  169. TVector3 pos;
  170. if (side == 0) point->Position(pos); // one detector side
  171. else point->PositionOut(pos); // the other side
  172. return new(clref[size]) MpdStsHit(detID, pos, TVector3(0.,0.,0.), point->GetEnergyLoss(), trackID, 0);
  173. }
  174. // ----
  175. ClassImp(MpdStsHitProducerV1)