MpdStrawendcapHitProducer.cxx 9.4 KB

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