MpdStsHitProducer.cxx 8.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247
  1. /////////////////////////////////////////////////////////////
  2. //
  3. // MpdStsHitProducer
  4. //
  5. // Filler of MpdStsHit
  6. //
  7. //
  8. //
  9. ///////////////////////////////////////////////////////////////
  10. #include "TClonesArray.h"
  11. #include "MpdStsGeoPar.h"
  12. #include "MpdStsHitProducer.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. MpdStsHitProducer::MpdStsHitProducer(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. MpdStsHitProducer::~MpdStsHitProducer() { }
  35. // -------------------------------------------------------------------------
  36. // ----- Public method Init --------------------------------------------
  37. InitStatus MpdStsHitProducer::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- MpdStsHitProducer::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- MpdStsHitProducer::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- MpdStsHitProducer: Intialization successfull" << endl;
  62. return kSUCCESS;
  63. }
  64. // -------------------------------------------------------------------------
  65. // ----- Private method SetParContainers -------------------------------
  66. void MpdStsHitProducer::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 MpdStsHitProducer::Exec(Option_t* opt) {
  77. //cout << " DIGI EXECUTION *********************" << endl;
  78. static Int_t eventCounter = 0, first = 1;
  79. // Layer Rmin, Rmax, module Z-size, layer Z-half-size
  80. static Double_t rMin[4] = {0.}, rMax[4] = {0.}, zMod[4] = {0.}, dZ[4] = {0.};
  81. static Double_t sterCos[2] = {-7.5, 7.5}, sterSin[2] = {0}; // stereo-angles
  82. const Double_t size = 6.2, pitch = 0.01; // 100 um
  83. if (first) {
  84. if ( ! fDigiArray ) Fatal("MpdStsHitProducer::Exec", "No DigiArray");
  85. // Get detector parameters
  86. first = 0;
  87. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  88. //rtdb->printParamContexts();
  89. //cout << rtdb->findContainer("TpcGeoPar") << " " << rtdb->findContainer("TpcContFact:") << endl;
  90. MpdStsGeoPar *geoPar = (MpdStsGeoPar*) rtdb->getContainer("MpdStsGeoPar");
  91. //cout << geoPar << endl;
  92. TString volName = "sts01 ";
  93. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  94. //cout << sensNodes->GetEntriesFast() << " " << geoPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  95. Int_t nLay = 4;
  96. for (Int_t i = 0; i < nLay; ++i) {
  97. volName[5] = 97 + i; // 'a', 'b', ..
  98. FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->FindObject(volName));
  99. //FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->At(0));
  100. TArrayD* params = sensVol->getParameters();
  101. rMin[i] = params->At(0);
  102. rMax[i] = params->At(1);
  103. //dR = (rMax-rMin) / 50;
  104. dZ[i] = params->At(2);
  105. Int_t nMods = Int_t (dZ[i] * 2. / size + 0.1);
  106. zMod[i] = dZ[i] * 2. / nMods;
  107. cout << " *** STS sensitive volume: " << sensVol->GetName() << " " << rMin[i] << " "
  108. << dZ[i] << " " << zMod[i] << " " << nMods << endl;
  109. }
  110. for (Int_t i = 0; i < 2; ++i) {
  111. sterCos[i] *= TMath::DegToRad(); // angle in rads
  112. sterSin[i] = TMath::Sin (sterCos[i]); // sin
  113. sterCos[i] = TMath::Cos (sterCos[i]); // cos
  114. }
  115. }
  116. cout << " StsHitProducer event " << ++eventCounter << endl;
  117. // Reset output array
  118. fDigiArray->Clear();
  119. // Declare some variables
  120. MpdStsPoint* point = NULL;
  121. map<Int_t, Float_t> fTrackEnergy;
  122. fTrackEnergy.clear();
  123. map<Int_t, Float_t>::const_iterator p;
  124. // Loop over MCPoints
  125. Int_t nPoints = fPointArray->GetEntriesFast();
  126. cout << nPoints << endl;
  127. for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
  128. point = (MpdStsPoint*) fPointArray->At(iPoint);
  129. Int_t detId = point->GetDetectorID();
  130. Int_t lay = (detId >> 1) - 1;
  131. Int_t module = Int_t ((point->GetZ() + dZ[lay]) / zMod[lay]);
  132. detId = detId | (module << 10);
  133. for (Int_t j = 0; j < 2; ++j) {
  134. // Two hits per point (for two sides)
  135. MpdStsHit *hit = AddHit(iPoint, detId+j, j);
  136. Double_t zloc = hit->GetZ() + dZ[lay];
  137. Double_t phi = TMath::ATan2 (hit->GetY(), hit->GetX());
  138. Double_t rad = (j == 0) ? rMin[lay] : rMax[lay];
  139. //Double_t x = rad * phi, xMax = TMath::TwoPi() * rad;
  140. Double_t x = rad * phi, xMax = TMath::TwoPi() * rad * sterCos[j];
  141. //Double_t xloc = x * sterCos[j] + zloc * sterSin[j]; // rotated coordinate
  142. Double_t xloc = x * sterCos[j] + hit->GetZ() * sterSin[j]; // rotated coordinate
  143. if (xloc < 0) xloc += xMax;
  144. //xloc0 = xloc0 - TMath::Floor(xloc0/xlocMx) * xlocMx; // modulo
  145. xloc = fmod (xloc, xMax);
  146. //cout << xloc0 << " " << fmod(xloc0,xlocMx) << endl;
  147. //if (TMath::Abs(xloc0-fmod(xloc0,xlocMx)) > 1.e-7) exit(0);
  148. Int_t strip = Int_t (xloc / pitch);
  149. hit->SetLocalX(xloc);
  150. hit->SetStrip(strip);
  151. //hit->SetDetID(Int_t sectorType, Int_t layer, Int_t ladder, Int_t det, Int_t side);
  152. hit->SetDetId(0, lay+1, 0, module, j);
  153. hit->SetUniqueID(1);
  154. }
  155. //cout << point->GetDetectorID() << " " << lay << " " << point->GetZ() << " " << module << " " << detId << endl;
  156. //fTrackEnergy[point->GetDetectorID()] += point->GetEnergyLoss();
  157. }
  158. // Loop to register ZdcHit
  159. /*
  160. for(p=fTrackEnergy.begin(); p!=fTrackEnergy.end(); ++p) {
  161. if ((*p).second>eneThr)
  162. AddHit(1, (*p).first, (*p).second);
  163. }
  164. */
  165. }
  166. // -------------------------------------------------------------------------
  167. // ----- Public method Create Structure --------------------------------
  168. void MpdStsHitProducer::CreateStructure() {
  169. /*
  170. TString work = getenv("VMCWORKDIR");
  171. work = work + "/geometry/" + fFileGeo;
  172. cout << "-I- <MpdStsHitProducer::CreateStructure> Zdc geometry loaded from: "
  173. << work << endl;
  174. Int_t detId = -1;
  175. MpdStsReader read(work);
  176. for(Int_t module=1; module<=read.GetMaxModules(); module++)
  177. for(Int_t row=1; row<=read.GetMaxRows(module); row++)
  178. for(Int_t crystal=1; crystal<=read.GetMaxCrystals(module,row); crystal++) {
  179. DataG4 data = read.GetData(module,row,crystal);
  180. for(Int_t copy=1; copy<=20; copy++) {
  181. detId = module*100000000 + row*1000000 + copy*10000 + crystal;
  182. emcX[detId] = data.posX; emcY[detId] = data.posY; emcZ[detId] = data.posZ;
  183. emcTheta[detId] = data.theta; emcTau[detId] = data.tau;
  184. if (module==3)
  185. emcPhi[detId] = fmod(data.phi+90.*(copy-1),360);
  186. else
  187. emcPhi[detId] = fmod(data.phi+22.5*(copy-1),360);
  188. }
  189. }
  190. */
  191. }
  192. // ----- Private method AddDigi --------------------------------------------
  193. MpdStsHit* MpdStsHitProducer::AddHit(Int_t trackID, Int_t detID, Int_t side)
  194. {
  195. // It fills the MpdStsHit category
  196. // cout << "MpdStsHitProducer: track " << trackID << " evt " << eventID << " sec " << sec << " plane " << pla << " strip " << strip << "box " << box << " tube " << tub << endl;
  197. MpdStsPoint *point = (MpdStsPoint*) fPointArray->UncheckedAt(trackID);
  198. TClonesArray& clref = *fDigiArray;
  199. Int_t size = clref.GetEntriesFast();
  200. TVector3 pos;
  201. if (side == 0) point->Position(pos); // one detector side
  202. else point->PositionOut(pos); // the other side
  203. return new(clref[size]) MpdStsHit(detID, pos, TVector3(0.,0.,0.), point->GetEnergyLoss(), trackID, 0);
  204. }
  205. // ----
  206. ClassImp(MpdStsHitProducer)