MpdTgemHitProducer.cxx 8.6 KB

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