MpdStrawECTHitProducer.cxx 8.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240
  1. //--------------------------------------------------------------------------
  2. //----- MpdStrawECTHitProducer -----
  3. //----- Created 8/12/09 by A. Zinchenko -----
  4. //--------------------------------------------------------------------------
  5. /** MpdStrawECTHitProducer.cxx
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** MPD End-Cap Tracker (ECT) hit producer
  9. **/
  10. #include "MpdStrawECTHitProducer.h"
  11. #include "MpdStrawECTHit.h"
  12. #include "MpdStrawECTPoint.h"
  13. #include "FairRun.h"
  14. #include "FairRuntimeDb.h"
  15. #include "FairRootManager.h"
  16. #include "FairDetector.h"
  17. #include <TClonesArray.h>
  18. #include "TGeoManager.h"
  19. #include "TGeoVolume.h"
  20. #include "TGeoNode.h"
  21. #include "TGeoMatrix.h"
  22. #include <TMath.h>
  23. #include <TRandom.h>
  24. #include "TVector3.h"
  25. #include <iostream>
  26. // ----- Default constructor -------------------------------------------
  27. MpdStrawECTHitProducer::MpdStrawECTHitProducer(const char* fileGeo) :
  28. FairTask("Ideal StrawECT hit Producer") {
  29. fFileGeo = fileGeo;
  30. eneThr = 0.001; // Energy threshold for StrawECT
  31. }
  32. // -------------------------------------------------------------------------
  33. // ----- Destructor ----------------------------------------------------
  34. MpdStrawECTHitProducer::~MpdStrawECTHitProducer() { }
  35. // -------------------------------------------------------------------------
  36. // ----- Public method Init --------------------------------------------
  37. InitStatus MpdStrawECTHitProducer::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- MpdStrawECTHitProducer::Init: "
  46. << "RootManager not instantiated!" << endl;
  47. return kFATAL;
  48. }
  49. // Get input array
  50. fPointArray = (TClonesArray*) ioman->GetObject("STRAWPoint");
  51. if ( ! fPointArray ) {
  52. cout << "-W- MpdStrawECTHitProducer::Init: "
  53. << "No StrawECTPoint array!" << endl;
  54. return kERROR;
  55. }
  56. // Create and register output array
  57. fDigiArray = new TClonesArray("MpdStrawECTHit");
  58. ioman->Register("EctHit","EctHits",fDigiArray,kTRUE);
  59. CreateStructure();
  60. cout << "-I- MpdStrawECTHitProducer: Intialization successfull" << endl;
  61. return kSUCCESS;
  62. }
  63. // -------------------------------------------------------------------------
  64. // ----- Public method Exec --------------------------------------------
  65. void MpdStrawECTHitProducer::Exec(Option_t* opt) {
  66. //cout << " DIGI EXECUTION *********************" << endl;
  67. static int eventCounter = 0;
  68. cout << " Event " << ++eventCounter << endl;
  69. // Reset output array
  70. if ( ! fDigiArray ) Fatal("Exec", "No DigiArray");
  71. fDigiArray->Clear();
  72. Bool_t mirror = kFALSE; //kTRUE; // add mirror hits if TRUE
  73. Int_t nPoints = fPointArray->GetEntriesFast(), lay = 0, nKH = 0, itube = 0;
  74. Double_t phi, r, coord, errR = 0.2, errRphi = 0.02; // 2mm in R, 200um in R-Phi
  75. //Double_t phi, r, coord, errR = 0.01, errRphi = 0.01; // 100um in R, 100um in R-Phi
  76. Double_t firstHit[100][1000] = {{0},{0}};
  77. Int_t indxTube[100][1000];
  78. for (Int_t i = 0; i < 100; ++i) {
  79. for (Int_t j = 0; j < 1000; ++j) indxTube[i][j] = -1;
  80. }
  81. TVector3 pos, dpos(0.,0.,0.);
  82. for (Int_t ih = 0; ih < nPoints; ++ih) {
  83. MpdStrawECTPoint *h = (MpdStrawECTPoint*) fPointArray->UncheckedAt(ih);
  84. if (h->GetZ() < 0) continue; // keep only Z>0
  85. //if (TMath::Abs(TMath::ATan2(h->GetTrackY(),h->GetTrackX())*TMath::RadToDeg()-97.5) > 7.5) continue; // !!!
  86. //if (h->GetTrackID() != 1427) continue; // !!!
  87. phi = h->GetPhi(); // tube Phi
  88. lay = h->GetDetectorID() / 1000; // + 1;
  89. itube = h->GetDetectorID() % 1000; // tube number
  90. // Extrapolate track to Z = Ztube
  91. Double_t dZ = h->GetZ() - h->GetTrackZ();
  92. Double_t dt = dZ / h->GetPz();
  93. if (TMath::Abs(h->GetPz()) > 1.e-6 && h->GetPz() * h->GetZ() > 0) dt = dZ / h->GetPz();
  94. Double_t xNew = h->GetTrackX() + dt * h->GetPx();
  95. Double_t yNew = h->GetTrackY() + dt * h->GetPy();
  96. Double_t zNew = h->GetTrackZ() + dt * h->GetPz();
  97. pos.SetXYZ(xNew,yNew,zNew);
  98. //cout << dZ << " " << h->GetTrackX() << " " << xNew << " " << h->GetTrackY() << " " << yNew << " " << lay << endl;
  99. // Transform to the rotated coordinate system
  100. Double_t cosPhi = TMath::Cos(phi);
  101. Double_t sinPhi = TMath::Sin(phi);
  102. //Double_t xLoc = h->GetX() * cosPhi + h->GetY() * sinPhi; // cross-check
  103. //Double_t yLoc = -h->GetX() * sinPhi + h->GetY() * cosPhi;
  104. Double_t xRot = xNew * cosPhi + yNew * sinPhi;
  105. Double_t yRot = -xNew * sinPhi + yNew * cosPhi;
  106. //Double_t xLoc = (xNew - h->GetX()) * cosPhi + (yNew - h->GetY()) * sinPhi;
  107. //Double_t yLoc = -(xNew - h->GetX()) * sinPhi + (yNew - h->GetY()) * cosPhi;
  108. //r = xNew * xNew + yNew * yNew;
  109. //r = TMath::Sqrt (r);
  110. //r = TMath::Abs(xLoc);
  111. r = xRot;
  112. //cout << xRot << " " << yRot << " " << r << " " << h->GetPz() << endl;
  113. coord = yRot;
  114. Double_t yWire = -h->GetX() * sinPhi + h->GetY() * cosPhi;
  115. Double_t dY = TMath::Abs (yWire - coord);
  116. // Add error
  117. Double_t dRphi = 0, dR = 0;
  118. gRandom->Rannor(dRphi,dR);
  119. MpdStrawECTHit *hit = new ((*fDigiArray)[nKH]) MpdStrawECTHit(h->GetDetectorID(), pos, dpos, ih, 0);
  120. hit->SetMeas(coord+dRphi*errRphi);
  121. hit->SetError(errRphi);
  122. hit->SetPhi(phi);
  123. hit->SetIndex(ih);
  124. //MpdKalmanHitZ *hit = new ((*fKHits)[nKH]) MpdKalmanHitZ(r+dR*errR, phi, coord+dRphi*errRphi,
  125. // h->GetTrackZ(), errR, errRphi, lay, ih);
  126. //hit->SetXY(h->GetX(), h->GetY());
  127. // Add second measurement - just for test at the moment
  128. //!!!
  129. //hit->SetNofDim(2);
  130. //!!!
  131. // Check if multiple hits per tube
  132. if (1 && indxTube[lay][itube] >= 0) {
  133. // Multiple hits
  134. MpdStrawECTHit *hitOld = (MpdStrawECTHit*) fDigiArray->UncheckedAt(indxTube[lay][itube]); // previous hit
  135. if (dY < firstHit[lay][itube]) {
  136. // Closer to anode wire - remove previous hit
  137. firstHit[lay][itube] = dY;
  138. const TArrayI *inds = hitOld->Index();
  139. Int_t nOld = inds->GetSize();
  140. for (Int_t io = 0; io < nOld; ++io) hit->SetIndex((*inds)[io]);
  141. fDigiArray->RemoveAt(indxTube[lay][itube]);
  142. indxTube[lay][itube] = nKH;
  143. } else {
  144. hitOld->SetIndex(hit->GetIndex());
  145. fDigiArray->RemoveAt(nKH); // remove last hit
  146. }
  147. } else {
  148. firstHit[lay][itube] = dY;
  149. indxTube[lay][itube] = nKH;
  150. if (mirror) {
  151. // Add mirror hit
  152. // wire point transverse coordinate in rotated system:
  153. /*
  154. yRot = -h->GetX() * sinPhi + h->GetY() * cosPhi;
  155. Double_t dy = hit->GetRphi() - yRot;
  156. MpdKalmanHitZ *hitM = new ((*fDigiArray)[++nKH]) MpdKalmanHitZ(*hit);
  157. hitM->SetRphi(yRot-dy);
  158. hitM->SetMirror(kTRUE);
  159. */
  160. }
  161. }
  162. ++nKH;
  163. }
  164. fDigiArray->Compress();
  165. fDigiArray->Sort(); // in ascending order in abs(Z)
  166. cout << " Max layer = " << ((MpdStrawECTHit*)fDigiArray->Last())->GetLayer() << " "
  167. << fDigiArray->GetEntriesFast() << endl;
  168. }
  169. // -------------------------------------------------------------------------
  170. // ----- Public method Create Structure --------------------------------
  171. void MpdStrawECTHitProducer::CreateStructure() {
  172. /*
  173. TString work = getenv("VMCWORKDIR");
  174. work = work + "/geometry/" + fFileGeo;
  175. cout << "-I- <MpdStrawECTHitProducer::CreateStructure> Zdc geometry loaded from: "
  176. << work << endl;
  177. Int_t detId = -1;
  178. MpdStrawECTReader read(work);
  179. for(Int_t module=1; module<=read.GetMaxModules(); module++)
  180. for(Int_t row=1; row<=read.GetMaxRows(module); row++)
  181. for(Int_t crystal=1; crystal<=read.GetMaxCrystals(module,row); crystal++) {
  182. DataG4 data = read.GetData(module,row,crystal);
  183. for(Int_t copy=1; copy<=20; copy++) {
  184. detId = module*100000000 + row*1000000 + copy*10000 + crystal;
  185. emcX[detId] = data.posX; emcY[detId] = data.posY; emcZ[detId] = data.posZ;
  186. emcTheta[detId] = data.theta; emcTau[detId] = data.tau;
  187. if (module==3)
  188. emcPhi[detId] = fmod(data.phi+90.*(copy-1),360);
  189. else
  190. emcPhi[detId] = fmod(data.phi+22.5*(copy-1),360);
  191. }
  192. }
  193. */
  194. }
  195. // ----- Private method AddDigi --------------------------------------------
  196. MpdStrawECTHit* MpdStrawECTHitProducer::AddHit(Int_t trackID,Int_t detID, Float_t energy){
  197. // It fills the MpdStrawECTHit category
  198. // cout << "MpdStrawECTHitProducer: track " << trackID << " evt " << eventID << " sec " << sec << " plane " << pla << " strip " << strip << "box " << box << " tube " << tub << endl;
  199. TClonesArray& clref = *fDigiArray;
  200. Int_t size = clref.GetEntriesFast();
  201. return new(clref[size]) MpdStrawECTHit(); // FIXME: real hit info needed here
  202. }
  203. // ----
  204. ClassImp(MpdStrawECTHitProducer)