MpdEmcHitCreation.cxx 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293
  1. ////////////////////////////////////////////////////////////////
  2. // //
  3. // MpdEmcHitCreation //
  4. // Hit production for EMC, v02 //
  5. // Author List : Martemianov M., 2019 //
  6. // //
  7. ////////////////////////////////////////////////////////////////
  8. #include "TClonesArray.h"
  9. #include "FairRootManager.h"
  10. #include "TMath.h"
  11. #include "TGeoNode.h"
  12. #include "TGeoManager.h"
  13. #include "FairRunAna.h"
  14. #include "FairEventHeader.h"
  15. #include "MpdEmcHitCreation.h"
  16. #include "MpdEmcGeoParams.h"
  17. #include "MpdEmcHit.h"
  18. #include "MpdEmcPoint.h"
  19. #include "FairMCPoint.h"
  20. #include "MpdMCTrack.h"
  21. using namespace std;
  22. using namespace TMath;
  23. // ----- Default constructor -------------------------------------------
  24. MpdEmcHitCreation::MpdEmcHitCreation() :
  25. FairTask("Ideal EMC hit Creation") {
  26. }
  27. // ----- Destructor ----------------------------------------------------
  28. MpdEmcHitCreation::~MpdEmcHitCreation() {
  29. }
  30. // -------------------------------------------------------------------------
  31. // ----- Public method Init --------------------------------------------
  32. InitStatus MpdEmcHitCreation::Init() {
  33. cout << "******************* EMC INIT *********************" << endl;
  34. // Get RootManager
  35. FairRootManager* ioman = FairRootManager::Instance();
  36. if (!ioman) {
  37. cout << "-E- MpdEmcHitCreation::Init: "
  38. << "RootManager not instantiated!" << endl;
  39. return kFATAL;
  40. }
  41. // Get geometry EMC parameters
  42. fGeoPar = new MpdEmcGeoParams();
  43. // Get input array
  44. fPointArray = (TClonesArray*) ioman->GetObject("EmcPoint");
  45. if (!fPointArray) {
  46. cout << "-W- MpdEmcHitCreation::Init: " << "No EmcPoint array!" << endl;
  47. return kERROR;
  48. }
  49. fMcTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
  50. if (!fMcTrackArray) {
  51. cout << "-W- MpdEmcHitCreation::Init: " << "No MCTrack array!" << endl;
  52. return kERROR;
  53. }
  54. // Create and register output array
  55. fDigiArray = new TClonesArray("MpdEmcHit", 100);
  56. ioman->Register("MpdEmcHit", "EMC", fDigiArray, kTRUE);
  57. ioman->Register("MCTrack","EMC",fMcTrackArray, kFALSE);
  58. cout << "-I- MpdEmcHitCreation: Intialization successfull" << endl;
  59. return kSUCCESS;
  60. }
  61. void MpdEmcHitCreation::Finish() {
  62. cout << "\n-I- MpdEmcHitCreation: Finish" << endl;
  63. }
  64. void MpdEmcHitCreation::Exec(Option_t* opt) {
  65. cout << "\n-I- MpdEmcHitCreation: Event No. " << FairRun::Instance()->GetEventHeader()->GetMCEntryNumber() << endl;
  66. // Reset output Array
  67. if (!fDigiArray) Fatal("MpdEmcHitCreation::Exec)", "No array of digits");
  68. fDigiArray->Delete();
  69. //Double_t phiRow, rhoMod, zMod, thetaMod;
  70. for (UInt_t iPnt = 0; iPnt < (UInt_t) fPointArray->GetEntriesFast(); ++iPnt) {
  71. FairMCPoint* pnt = (FairMCPoint*) fPointArray->At(iPnt);
  72. Int_t trId = pnt->GetTrackID();
  73. if (trId < 0) break;
  74. MpdMCTrack* tr = (MpdMCTrack*) fMcTrackArray->At(trId);
  75. // if (tr->GetMotherId() != -1) continue;
  76. Int_t pdg = tr->GetPdgCode();
  77. Float_t x = pnt->GetX();
  78. Float_t y = pnt->GetY();
  79. Float_t z = pnt->GetZ();
  80. // Path to point
  81. //TGeoNode* currentNode = gGeoManager->FindNode(x,y,z);
  82. //TString path(gGeoManager->GetPath());
  83. // Get sector number of module
  84. Int_t sec = GetSecId(x, y);
  85. // Get row number of module
  86. Int_t row = GetRowId(x, y);
  87. // Get module number in row
  88. Double_t xTower, yTower, zTower, phiTower, thetaTower;
  89. Int_t tower = GetTowerId(x, y, z, row, xTower, yTower, zTower, phiTower, thetaTower);
  90. Int_t supMod = -1;
  91. Float_t e = pnt->GetEnergyLoss();
  92. Float_t timeEnergy = pnt->GetTime()*e;
  93. MpdEmcHit* hit = SearchHit(sec, row, tower);
  94. if (hit == NULL) {
  95. hit = new((*fDigiArray)[fDigiArray->GetEntriesFast()]) MpdEmcHit(sec, row, supMod, tower, e, timeEnergy);
  96. } else {
  97. hit->IncreaseEnergy(e);
  98. hit->IncreaseEnergyTime(timeEnergy);
  99. }
  100. hit->SetNumTracks(hit->GetNumTracks() + 1);
  101. hit->SetTrackId(trId);
  102. hit->SetPdg(pdg);
  103. hit->SetRhoCenter(sqrt(xTower*xTower+yTower*yTower));
  104. hit->SetZCenter(zTower);
  105. hit->SetPhiCenter(phiTower);
  106. hit->SetThetaCenter(thetaTower);
  107. }
  108. for (int iHit = 0; iHit < fDigiArray->GetEntriesFast(); iHit++) {
  109. MpdEmcHit* hit = (MpdEmcHit*) fDigiArray->At(iHit);
  110. hit->SetTime(hit->GetTime()/hit->GetE());
  111. if (hit->GetNumTracks() > 1) {
  112. hit->SetPdg(0);
  113. hit->SetTrackId(-1);
  114. }
  115. }
  116. }
  117. // Define sector number
  118. Int_t MpdEmcHitCreation::GetSecId(Double_t x, Double_t y) {
  119. Int_t index = -1;
  120. vector<double> diffSector;
  121. const Int_t nSec = fGeoPar->GetNsec();
  122. vector<double> phiSector = fGeoPar->GetPhiSector();
  123. Double_t ang = ATan2(y,x)*RadToDeg();
  124. const Double_t fSectorAngle = fGeoPar->GetSizeAngleSector();
  125. if (ang < 0) ang += 360.;
  126. for (int iSector = 0; iSector < nSec; iSector++)
  127. diffSector.push_back(fabs(phiSector[iSector] - ang));
  128. index = LocMin(diffSector.size(), &diffSector[0]);
  129. if ( (ang > 0) && (ang < phiSector[0] - 0.5*fSectorAngle) ) index = phiSector.size() - 1;
  130. return index;
  131. }
  132. // Define row number
  133. Int_t MpdEmcHitCreation::GetRowId(Double_t x, Double_t y) {
  134. Int_t index = -1;
  135. vector<double> diffR;
  136. const Int_t nRow = fGeoPar->GetNrows();
  137. vector<double> phiRow = fGeoPar->GetPhiRow();
  138. vector<double> xRow = fGeoPar->GetXRow();
  139. vector<double> yRow = fGeoPar->GetYRow();
  140. Float_t ang, phiAngle;
  141. for (int iRow = 0; iRow < nRow; iRow++) {
  142. ang = ATan2(y - yRow[iRow],x - xRow[iRow])*RadToDeg();
  143. phiAngle = phiRow[iRow];
  144. if ( ( phiAngle == 360.) && ( ang > 0) ) phiAngle = 0.0;
  145. if (ang < 0) ang += 360.;
  146. diffR.push_back(fabs(phiAngle - ang));
  147. }
  148. Int_t rowMod = LocMin(diffR.size(), &diffR[0]);
  149. ang = ATan2(y - yRow[rowMod],x - xRow[rowMod])*RadToDeg();
  150. phiAngle = phiRow[rowMod];
  151. if ( (phiAngle == 360.) && (ang > 0) ) phiAngle = 0.0;
  152. if (ang < 0) ang += 360.;
  153. index = 2*rowMod + 1;
  154. if (phiAngle - ang > 0) index = index - 1;
  155. return index;
  156. }
  157. Int_t MpdEmcHitCreation::GetTowerId(Double_t x, Double_t y, Double_t z, Int_t iRow, Double_t &xTower,
  158. Double_t &yTower, Double_t &zTower, Double_t &phiTower, Double_t &thetaTower) {
  159. Int_t index = -1;
  160. vector<double> diffTheta;
  161. const Int_t nRow = fGeoPar->GetNrows();
  162. const Int_t nBox = fGeoPar->GetNmod();
  163. vector<Double_t> xBox = fGeoPar->GetXBox();
  164. vector<Double_t> yBox = fGeoPar->GetYBox();
  165. vector<Double_t> zBox = fGeoPar->GetZBox();
  166. vector<Double_t> thetaBox = fGeoPar->GetThetaBox();
  167. vector<Double_t> phiBox = fGeoPar->GetPhiBox();
  168. Int_t boxStep = nBox/(2*nRow);
  169. const Double_t boxSizeLow = fGeoPar->GetSizeLowBox();
  170. const Double_t boxSizeHigh = fGeoPar->GetSizeHighBox();
  171. const Double_t boxLength = fGeoPar->GetLengthBox();
  172. Double_t ang = ATan2(y, x);
  173. if (ang < 0) ang += TwoPi();
  174. Double_t rhoMod0, zMod0, lenTotal, thetaRad, boxRho; //thetaPoint,
  175. Double_t dPolTheta = ATan2(0.5*(boxSizeHigh - boxSizeLow),boxLength)*RadToDeg();
  176. Int_t iRowStart, iRowNext = iRow*boxStep;//, iRowNum = iRow/2;
  177. if ( Odd(iRow) ) {iRowNext = (iRow-1)*boxStep + 1; /*iRowNum = (iRow - 1)/2.;*/}
  178. iRowStart = iRowNext;
  179. for (Int_t iBox = 0; iBox < boxStep; iBox++) {
  180. boxRho = sqrt(xBox[iRowNext]*xBox[iRowNext]+yBox[iRowNext]*yBox[iRowNext]);
  181. thetaRad = thetaBox[iRowNext]*DegToRad();
  182. lenTotal = boxSizeHigh/tan(dPolTheta*DegToRad()) - boxLength - boxRho/sin(thetaRad);
  183. rhoMod0 = sin(thetaRad)*lenTotal;
  184. zMod0 = zBox[iRowNext] - boxRho/tan(thetaRad) - lenTotal*cos(thetaRad);
  185. diffTheta.push_back(fabs(thetaBox[iRowNext] -
  186. ATan2(sqrt(x*x + y*y) + rhoMod0, fabs(z) - zMod0)*RadToDeg()));
  187. iRowNext += 2;
  188. }
  189. index = LocMin(diffTheta.size(), &diffTheta[0]);
  190. iRowNext = iRowStart + 2*index;
  191. xTower = xBox[iRowNext]; yTower = yBox[iRowNext]; zTower = zBox[iRowNext];
  192. phiTower = phiBox[iRowNext]; thetaTower = thetaBox[iRowNext];
  193. if (z < 0) {
  194. index = boxStep - index - 1;
  195. zTower = - zTower;
  196. }
  197. else {
  198. index = boxStep + index;
  199. }
  200. return index;
  201. }
  202. MpdEmcHit* MpdEmcHitCreation::SearchHit(UInt_t sec, UInt_t row, UInt_t mod) {
  203. MpdEmcHit* foundHit = NULL;
  204. for (Int_t i = 0; i < fDigiArray->GetEntriesFast(); ++i) {
  205. MpdEmcHit* hit = (MpdEmcHit*) fDigiArray->At(i);
  206. if ((UInt_t)hit->GetSec() != sec) continue;
  207. if ((UInt_t)hit->GetRow() != row) continue;
  208. if ((UInt_t)hit->GetMod() != mod) continue;
  209. foundHit = hit;
  210. }
  211. return foundHit;
  212. }
  213. // -------------------------------------------------------------------------
  214. ClassImp(MpdEmcHitCreation)