MpdItsHitProducer5spd.cxx 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264
  1. // -------------------------------------------------------------------------
  2. // ----- MpdItsHitProducer5spd source file -----
  3. // ----- 13.12.2016 VK -----
  4. // -------------------------------------------------------------------------
  5. // Independant pixel readout for each sector
  6. // Geometry: its_5spd.geo
  7. // spd: one type of sector consisting of one sensor (15*30 mm)
  8. #include "MpdItsHitProducer5spd.h"
  9. #include "MpdStsGeoPar.h"
  10. #include "MpdItsHit5spd.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 "TClonesArray.h"
  19. #include "TGeoBBox.h"
  20. #include "TGeoManager.h"
  21. #include "TGeoVolume.h"
  22. #include "TGeoNode.h"
  23. #include "TGeoMatrix.h"
  24. #include <TRandom.h>
  25. #include "TVector3.h"
  26. // ----- Default constructor -------------------------------------------
  27. MpdItsHitProducer5spd::MpdItsHitProducer5spd() :
  28. FairTask("Ideal ITS hit Producer5spd"),
  29. fPointArray(nullptr),
  30. fDigiArray(nullptr)
  31. {}
  32. // ----- Default constructor -------------------------------------------
  33. /*
  34. MpdItsHitProducer5spd::MpdItsHitProducer5spd(const char* fileGeo) :
  35. FairTask("Ideal ITS hit Producer_5spd") {
  36. fFileGeo=fileGeo;
  37. eneThr = 0.001; // Energy threshold for ITS
  38. fDigiArray = nullptr;
  39. }
  40. */
  41. // ----- Destructor ----------------------------------------------------
  42. //MpdItsHitProducer5spd::~MpdItsHitProducer5spd() { delete fDigiArray; fDigiArray = nullptr; }
  43. MpdItsHitProducer5spd::~MpdItsHitProducer5spd() { }
  44. // ----- Public method Init --------------------------------------------
  45. InitStatus MpdItsHitProducer5spd::Init() {
  46. cout << "-I- MpdItsHitProducer5spd: Intialization started" << endl;
  47. // Get RootManager
  48. FairRootManager* ioman = FairRootManager::Instance();
  49. if ( ! ioman ) {
  50. cout << "-E- MpdItsHitProducer5spd::Init: "
  51. << "RootManager not instantiated!" << endl;
  52. return kFATAL;
  53. }
  54. // Get input array
  55. fPointArray = (TClonesArray*) ioman->GetObject("StsPoint");
  56. if ( ! fPointArray ) {
  57. cout << "-W- MpdItsHitProducer5spd::Init: "
  58. << "No ItsPoint array!" << endl;
  59. return kERROR;
  60. }
  61. // Create and register output array
  62. fDigiArray = new TClonesArray("MpdItsHit5spd");
  63. //ioman->Register("StsHit","Sts",fDigiArray,kTRUE);
  64. ioman->Register("StsHit","Sts",fDigiArray,kFALSE);
  65. CreateStructure();
  66. cout << "-I- MpdItsHitProducer5spd: Intialization finished successfully" << endl;
  67. return kSUCCESS;
  68. }
  69. // ----- Private method SetParContainers -------------------------------
  70. void MpdItsHitProducer5spd::SetParContainers() {
  71. // Get run and runtime database
  72. FairRunAna* run = FairRunAna::Instance();
  73. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  74. FairRuntimeDb* db = run->GetRuntimeDb();
  75. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  76. // Get STS geometry parameter container
  77. db->getContainer("MpdStsGeoPar");
  78. }
  79. // ----- Public method Exec --------------------------------------------
  80. void MpdItsHitProducer5spd::Exec(Option_t* opt) {
  81. cout << "-I- MpdItsHitProducer5spd: Execution started" << endl;
  82. fDigiArray->Delete();
  83. static Int_t eventCounter = 0, first = 1;
  84. //static Double_t par1[9] = {0}, par2[9] = {0}, par3[9] = {0},
  85. static Double_t spd_length = 0, spd_width = 0; // Detector parameters
  86. Int_t detId = 0;//, sector = 0, ladder = 0, layer = 0;
  87. //------ spd parameters
  88. //const Double_t spd_length = 3.0; // cm, length of spd sensor
  89. //const Double_t spd_width = 1.5; // cm, width of spdsensor
  90. const Int_t nCol = 1024; // number of colomns in spd matrix
  91. const Int_t nRow = 512; // number of rows in spd matrix
  92. cout << " *** ItsHitProducer5spd event " << ++eventCounter << endl;
  93. /*
  94. if (first) {
  95. first = 0;
  96. //------ Define sensitive volumes
  97. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  98. MpdStsGeoPar *geoPar = (MpdStsGeoPar*) rtdb->getContainer("MpdStsGeoPar");
  99. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  100. cout << "Number of sensors: "<< sensNodes->GetEntriesFast() << endl;
  101. Int_t nSensors = sensNodes->GetEntriesFast();
  102. // !!! AZ - Crashes at root exit
  103. for (Int_t i = 0; i < nSensors; ++i) {
  104. FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->At(i));
  105. TArrayD* params = sensVol->getParameters();
  106. par1[i] = params->At(0); // half length of sensor
  107. par2[i] = params->At(1); // half tickness of sensor
  108. par3[i] = params->At(2); // half width of sensor
  109. // cout << " *** STS sensitive volume: " << sensVol->GetName() << " " << params->At(0) << " "
  110. // << params->At(1) << " " << params->At(2) << endl;
  111. }
  112. }
  113. */
  114. //------ Reset output array
  115. fDigiArray->Clear();
  116. //------ Declare some variables
  117. MpdStsPoint* point = nullptr;
  118. //------ Loop over MCPoints
  119. Int_t nPoints = fPointArray->GetEntriesFast();
  120. cout << "********* Number of STS points: " << nPoints << endl;
  121. for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
  122. point = (MpdStsPoint*) fPointArray->At(iPoint);
  123. // cout << "# of current STS point: " << iPoint+1 << endl;
  124. detId = point->GetDetectorID();
  125. //sector = (detId & 127);
  126. //ladder = ((detId >> 7) & 63);
  127. //layer = ((detId >> 13) & 7);
  128. Double_t pos[3] = {0};
  129. //AZ TVector3 pos3;
  130. TVector3 pos3, pos3out, mom3; //AZ-161020
  131. point->Position(pos3);
  132. //AZ-161020 pos3.GetXYZ(pos); // global coordinate
  133. //AZ-161020 TGeoNode *node = gGeoManager->FindNode(pos[0],pos[1],pos[2]);
  134. TGeoNode *node = gGeoManager->FindNode(pos3[0],pos3[1],pos3[2]);
  135. if (node == 0x0) {
  136. cout << " !!! Missing node " << pos3[0] << " " << pos3[1] << " " << pos3[2] << endl;
  137. exit(0);
  138. }
  139. //cout << "Node: " << node->GetName() << " Point position: " << pos[0] << " " << pos[1] << " " << pos[2] << endl;
  140. if (first) {
  141. first = 0;
  142. spd_length = ((TGeoBBox*)node->GetVolume()->GetShape())->GetDZ() * 2;
  143. spd_width = ((TGeoBBox*)node->GetVolume()->GetShape())->GetDX() * 2;
  144. }
  145. //cout << gGeoManager->GetPath() << " DetectorID = " << detId << endl;
  146. //cout << " Layer: "<< layer << " Ladder: " << ladder << " Sector: " << sector << endl;
  147. //AZ-161020
  148. point->PositionOut(pos3out);
  149. point->Momentum(mom3);
  150. mom3.SetMag(0.0002); // step 2 um
  151. mom3 *= -1;
  152. pos3out += mom3; // step back
  153. pos3 += pos3out;
  154. pos3 *= 0.5;
  155. pos3.GetXYZ(pos); // global coordinate
  156. //AZ
  157. Double_t loc[3] = {0.};
  158. gGeoManager->MasterToLocal(pos,loc);
  159. Double_t x = loc[0], y = loc[1], z = loc[2]; // local coordinate
  160. Double_t xlocSensor;
  161. Int_t iCol, iRow;
  162. if (x < -0.5*spd_width || x > 0.5*spd_width || z < -0.5*spd_length || z > 0.5*spd_length) {
  163. cout << "***** Point is outside SPD sensor! " << endl;
  164. cout << "***** Point is in volume: " << gGeoManager->GetPath() << endl;
  165. cout << "***** Local coordinates of point are: " << x << " " << y << " " << z << " " << endl;
  166. //AZ x = -spd_width*(0.5+1./nRow); // iRow=-1
  167. z = -spd_length*(0.5+1./nCol); // iCol=-1
  168. }
  169. iCol = nCol*(z/spd_length + 0.5);
  170. iRow = nRow*(x/spd_width + 0.5);
  171. xlocSensor = x;
  172. // cout << "Global coordinates of point in sensor: " << pos[0] << " " << pos[1] << " " << pos[2] << " " << endl;
  173. // cout << "Local coordinates of point in sensor: " << loc[0] << " " << loc[1] << " " << loc[2] << " " << endl;
  174. // cout << " Col= " << iCol << " Row= " << iRow << " Layer: " << layer << endl;
  175. MpdItsHit5spd *hit = AddHit(iPoint, detId); // One hit per point
  176. hit->SetUniqueID(2);
  177. hit->SetTrackID(iPoint);
  178. hit->SetLocalX(xlocSensor);
  179. hit->SetCol(iCol);
  180. hit->SetRow(iRow);
  181. } // points
  182. }
  183. //------ Public method Create Structure ------------------------------------
  184. void MpdItsHitProducer5spd::CreateStructure() {
  185. }
  186. //------ Private method AddDigi --------------------------------------------
  187. MpdItsHit5spd* MpdItsHitProducer5spd::AddHit(Int_t trackID, Int_t detID)
  188. {
  189. // It fills the MpdStsHit category
  190. MpdStsPoint *point = (MpdStsPoint*) fPointArray->UncheckedAt(trackID);
  191. TClonesArray& clref = *fDigiArray;
  192. Int_t size = clref.GetEntriesFast();
  193. TVector3 pos, posOut;
  194. point->Position(pos); // one detector side
  195. point->PositionOut(posOut); // the other side
  196. pos += posOut;
  197. pos *= 0.5;
  198. Double_t errx, errz;
  199. gRandom->Rannor(errx, errz);
  200. return new(clref[size]) MpdItsHit5spd(detID, pos, TVector3(errx,0.,errz), point->GetEnergyLoss(), trackID, 0);
  201. }
  202. //__________________________________________________________________________
  203. void MpdItsHitProducer5spd::Finish()
  204. {
  205. ///
  206. delete fDigiArray;
  207. fDigiArray = nullptr;
  208. }
  209. ClassImp(MpdItsHitProducer5spd)