MpdDchHitProducer.cxx 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424
  1. //------------------------------------------------------------------------------------------------------------------------
  2. #include <iostream>
  3. #include <assert.h>
  4. #include <TRandom2.h>
  5. #include <TClonesArray.h>
  6. #include <TGeoManager.h>
  7. #include <TGeoVolume.h>
  8. #include <TGeoNode.h>
  9. #include <TGeoMatrix.h>
  10. #include <TMath.h>
  11. #include <TRandom.h>
  12. #include <TVector3.h>
  13. #include <TH2D.h>
  14. #include "MpdDch.h"
  15. #include "MpdDchHit.h"
  16. #include "MpdDchPoint.h"
  17. #include "FairRun.h"
  18. #include "FairRuntimeDb.h"
  19. #include "FairRootManager.h"
  20. #include "FairDetector.h"
  21. #include "FairMCTrack.h"
  22. struct __ltstr
  23. {
  24. bool operator()(Double_t s1, Double_t s2) const
  25. {
  26. return s1 > s2;
  27. }
  28. };
  29. #include "MpdDchHitProducer.h"
  30. //------------------------------------------------------------------------------------------------------------------------
  31. MpdDchHitProducer::MpdDchHitProducer(const char* fileGeo, Int_t verbose, Bool_t test)
  32. : FairTask("Dch HitProducer", verbose), fDoTest(test), fRSigma(0.2000), fRPhiSigma(0.0200)
  33. {
  34. pRandom = new TRandom2;
  35. if(fDoTest)
  36. {
  37. htTime = new TH1D("htTime", "Time delay, ns", 1000, 0., 1000.); htTime->SetDirectory(0); fList.Add(htTime);
  38. htTime->SetXTitle("#Delta_{time delay}, ns"); htTime->SetYTitle("Events");
  39. htTimeA = (TH1D*) htTime->Clone("htTimeA"); htTimeA->SetDirectory(0); fList.Add(htTimeA);
  40. htGasDrift = new TH1D("htGasDrift", "", 100, 0., 0.5); htGasDrift->SetDirectory(0); fList.Add(htGasDrift);
  41. htGasDrift->SetXTitle("#Delta_{gas drift}, cm"); htGasDrift->SetYTitle("Events");
  42. htGasDriftA = (TH1D*) htGasDrift->Clone("htGasDriftA"); htGasDriftA->SetDirectory(0); fList.Add(htGasDriftA);
  43. htPerp = new TH1D("htPerp", "", 360, -140.+0.25, 140.+0.25); htPerp->SetDirectory(0); fList.Add(htPerp);
  44. htPerp->SetXTitle("wire position, cm"); htPerp->SetYTitle("Events");
  45. htPerpA = (TH1D*) htPerp->Clone("htPerpA"); htPerpA->SetDirectory(0); fList.Add(htPerpA);
  46. htOccup = new TH1D("htOccup", "Cell occupancy", 100, 0.5, 100.5);
  47. htOccup->SetDirectory(0); fList.Add(htOccup);
  48. htXYlocal = new TH2D("htXYlocal", "Local XY", 1000, 0., 0.6, 1000, -200., 200.);
  49. htXYlocal->SetDirectory(0); fList.Add(htXYlocal);
  50. htRvsR = new TH2D("htRvsR", "R global vs R local", 1000, 0., 150., 1000, 0., 150.);
  51. htRvsR->SetDirectory(0); fList.Add(htRvsR);
  52. htWireN = new TH1D("htWireN", "#wire", 16400, -200., 16000. + 200.);
  53. htWireN->SetDirectory(0); fList.Add(htWireN);
  54. htMCTime = new TH1D("htMCTime", "MC time", 100, 0., 1000.);
  55. htMCTime->SetDirectory(0); fList.Add(htMCTime);
  56. }
  57. }
  58. //------------------------------------------------------------------------------------------------------------------------
  59. MpdDchHitProducer::~MpdDchHitProducer()
  60. {
  61. delete pRandom;
  62. }
  63. //------------------------------------------------------------------------------------------------------------------------
  64. InitStatus MpdDchHitProducer::Init()
  65. {
  66. Info("Init", "Begin hit producer initialization.");
  67. FairRootManager *ioman = FairRootManager::Instance();
  68. if(ioman==0){ Error("MpdDchHitProducer::Init","FairRootManager XUINJA"); return kERROR; }
  69. pDchPoints = (TClonesArray *) ioman->GetObject("DchPoint");
  70. pMCTracks = (TClonesArray *) ioman->GetObject("MCTrack");
  71. if(!pDchPoints || !pMCTracks){ Error("MpdDchHitProducer::Init","Branch not found!"); return kERROR; }
  72. // Create and register output array
  73. pHitCollection = new TClonesArray("MpdDchHit");
  74. ioman->Register("DchHit","Dch", pHitCollection, kTRUE);
  75. Info("Init", "Initialization finished succesfully.");
  76. return kSUCCESS;
  77. }
  78. //------------------------------------------------------------------------------------------------------------------------
  79. void MpdDchHitProducer::Rotate(Int_t proj, Double_t x,Double_t y, Double_t& xRot, Double_t& yRot, Bool_t back)
  80. {
  81. // Transform to the rotated coordinate system
  82. // [0-3]==[x,y,u,v] 0, 90, -45, 45
  83. static const Double_t cosPhi_45 = TMath::Cos(-45.*TMath::DegToRad()), sinPhi_45 = TMath::Sin(-45.*TMath::DegToRad()),
  84. cosPhi45 = TMath::Cos(45.*TMath::DegToRad()), sinPhi45 = TMath::Sin(45.*TMath::DegToRad()) ;
  85. assert(proj<4);
  86. const Int_t map[]={0, 4, 3, 2};// 0<-->0, 1<-->4, 2<-->3, 3<-->2
  87. if(back) proj = map[proj];
  88. //Double_t u = -h->GetX() * cosSin[1] + h->GetY() * cosSin[0];
  89. //Double_t v = h->GetX() * cosSin[0] + h->GetY() * cosSin[1];
  90. switch(proj)
  91. {
  92. case 0: // 0 degree
  93. //xRot = x;
  94. //yRot = y;
  95. xRot = y;
  96. yRot = x;
  97. return;
  98. case 1: // 90 degree
  99. //xRot = y;
  100. //yRot = -x;
  101. xRot = -x;
  102. yRot = y;
  103. return;
  104. case 2: // -45 degree
  105. //xRot = x * cosPhi_45 + y * sinPhi_45;
  106. //yRot = -x * sinPhi_45 + y * cosPhi_45;
  107. xRot = -x * sinPhi_45 + y * cosPhi_45;
  108. yRot = x * cosPhi_45 + y * sinPhi_45;
  109. return;
  110. case 3: // 45 degree
  111. //xRot = x * cosPhi45 + y * sinPhi45;
  112. //yRot = -x * sinPhi45 + y * cosPhi45;
  113. xRot = -x * sinPhi45 + y * cosPhi45;
  114. yRot = x * cosPhi45 + y * sinPhi45;
  115. return;
  116. case 4: // -90 degree
  117. //xRot = -y;
  118. //yRot = x;
  119. xRot = x;
  120. yRot = -y;
  121. return;
  122. default: assert(false);
  123. }
  124. }
  125. //------------------------------------------------------------------------------------------------------------------------
  126. inline Double_t MpdDchHitProducer::GetPhi(Int_t proj)
  127. {
  128. static const Double_t Phi_45 = -45.*TMath::DegToRad(), Phi45 = 45.*TMath::DegToRad(), Phi90 = 90.*TMath::DegToRad();
  129. switch(proj)
  130. {
  131. case 0: return 0.; // 0 degree
  132. case 1: return Phi90; // 90 degree
  133. case 2: return Phi_45; // -45 degree
  134. case 3: return Phi45; // 45 degree
  135. default: assert(false);
  136. }
  137. }
  138. //------------------------------------------------------------------------------------------------------------------------
  139. Double_t MpdDchHitProducer::GetDriftLenght(Int_t proj, Int_t gasgap, Double_t x, Double_t& wirePos)
  140. {
  141. // ... -1 0 1 ... - first(0) gap wire position (X) [cm]
  142. // ... -0.5 0.5 1.5 ... - second(1) gap wire position (X) [cm]
  143. wirePos = (gasgap == 0) ? TMath::Nint(x) : TMath::Nint(x + 0.5) - 0.5;
  144. //AZ return TMath::Abs(x - wirePos); // [cm]
  145. return x - wirePos; // [cm]
  146. }
  147. //------------------------------------------------------------------------------------------------------------------------
  148. Bool_t MpdDchHitProducer::HitExist(Double_t delta) // [ns]
  149. {
  150. //const static Double_t _Time[4] = {0., 20., 20.01, 20.02}; // 0., 2., 10., 40.
  151. const static Double_t _Time[4] = {0., 2., 2.01, 2.02}; // 0., 2., 10., 40.
  152. const static Double_t _Eff[4] = {0., 0.01, 0.99, 1.0}; // 0., 0.2, 0.8, 1.
  153. // Efficiency
  154. //-------------------------------------
  155. // 1. x
  156. // 0.9 x
  157. //
  158. //
  159. // 0.6 x
  160. //
  161. // 0 x
  162. // 0 5 10 50 ns
  163. //-------------------------------------
  164. const static Double_t slope1 = (_Eff[1] - _Eff[0]) / (_Time[1] - _Time[0]);
  165. const static Double_t slope2 = (_Eff[2] - _Eff[1]) / (_Time[2] - _Time[1]);
  166. const static Double_t slope3 = (_Eff[3] - _Eff[2]) / (_Time[3] - _Time[2]);
  167. Double_t efficiency;
  168. if( delta > _Time[3]) return true;
  169. else if(delta > _Time[2] && delta < _Time[3]) efficiency = _Eff[2] + (delta - _Time[2]) * slope3;
  170. else if(delta > _Time[1] && delta < _Time[2]) efficiency = _Eff[1] + (delta - _Time[1]) * slope2;
  171. else if( delta < _Time[1]) efficiency = delta * slope1;
  172. ///cout<<"\n eff="<<efficiency<<" "<<delta;
  173. if(pRandom->Rndm() < efficiency) return true;
  174. return false;
  175. }
  176. //------------------------------------------------------------------------------------------------------------------------
  177. Double_t MpdDchHitProducer::GetTShift(Double_t driftLength, Double_t wirePos, Double_t R, Double_t& L)
  178. {
  179. const static Double_t gasDriftSpeed = 5.e-3; // 50 mkm/ns == 5.e-3 cm/ns
  180. const static Double_t wireDriftSpeed = 20; // 5 ns/m == 20 cm/ns
  181. const static Double_t WeelR_2 = 135*135; // cm
  182. driftLength = driftLength > 0 ? driftLength : -1.*driftLength;
  183. L = sqrt(WeelR_2 - wirePos*wirePos); // half wire length
  184. if(wirePos > -9.6 && wirePos < 9.6) L = L - TMath::Abs(R);// two wires
  185. else L = L + R; // one wire
  186. ////cout<<"\n t1="<<driftLength / gasDriftSpeed<<" ("<<driftLength<<") t2="<<(wireLength + R) / wireDriftSpeed<<" L="<<wireLength<<" Y="<<R;
  187. return driftLength / gasDriftSpeed + L / wireDriftSpeed;
  188. }
  189. //------------------------------------------------------------------------------------------------------------------------
  190. Int_t MpdDchHitProducer::WireID(Int_t uid, Double_t wirePos, Double_t R)
  191. {
  192. uid--; // uid [0,15]
  193. // wirePos [-135,135]cm
  194. // tube R 9.6 cm
  195. if(wirePos > -9.6 && wirePos < 9.6) if(R > 0)return (int)(wirePos + 1000.*uid + 500); // two wires
  196. return (int)( wirePos + 1000.*uid); // one wire
  197. }
  198. //------------------------------------------------------------------------------------------------------------------------
  199. void MpdDchHitProducer::Exec(Option_t* opt)
  200. {
  201. pHitCollection->Delete();
  202. Int_t hitID = 0;
  203. fMapOccup.clear();
  204. Int_t nDchPoint = pDchPoints->GetEntriesFast();
  205. cout<<"\n-I- [MpdDchHitProducer::Exec] "<<nDchPoint<<" points in Dch for this event."<<flush;
  206. MpdDchPoint *pPoint;
  207. MpdDchHit *pHit;
  208. Double_t R, Rphi, x, y, xRot, yRot, driftLength, wirePos, L;
  209. Double_t dRphi = 0, dR = 0;
  210. TVector3 pos, dpos;
  211. Int_t uid, wheel, gasgap, proj;
  212. for(Int_t pointIndex = 0; pointIndex < nDchPoint; pointIndex++ ) // <---Loop over the DCH points
  213. {
  214. pPoint = (MpdDchPoint*) pDchPoints->UncheckedAt(pointIndex);
  215. /////// if( ((FairMCTrack*)pMCTracks->UncheckedAt(pPoint->GetTrackID()))->GetMotherId() != -1)continue; // primary ONLY !!!
  216. uid = pPoint->GetDetectorID();
  217. wheel = MpdDch::GetWheel(uid); // [0-1] == [inner,outer]
  218. proj = MpdDch::GetProj(uid); // [0-3] == [x,y,u,v]
  219. gasgap = MpdDch::GetGasGap(uid); // [0-1] == [inner,outer]
  220. Rotate(proj, x = pPoint->GetX(), y = pPoint->GetY(), xRot, yRot); // GlobalToLocal
  221. driftLength = GetDriftLenght(proj, gasgap, xRot, wirePos); // [cm]
  222. R = yRot;
  223. Rphi = driftLength;
  224. pRandom->Rannor(dRphi,dR);
  225. if(fDoTest)
  226. {
  227. htWireN->Fill(WireID(uid, wirePos, R));
  228. htXYlocal->Fill(Rphi, R);
  229. htRvsR->Fill(sqrt(x*x + y*y), sqrt(R*R + wirePos*wirePos));
  230. htMCTime->Fill(pPoint->GetTime());
  231. }
  232. ///// xRot = pRandom->Gaus(xRot, fRSigma/10.); // [cm]
  233. ///// Rotate(proj, xRot, yRot, x, y, true); // back rotate
  234. ///// pos.SetXYZ(x, y, pPoint->GetZ());
  235. pos.SetXYZ(pPoint->GetX(), pPoint->GetY(), pPoint->GetZ());
  236. dpos.SetXYZ(0,0,0); ////////// ???????
  237. pHit = AddHit(hitID, uid, pos, dpos, pPoint->GetTrackID(), pointIndex, 0);
  238. pHit->SetPhi(GetPhi(proj));
  239. pHit->SetMeas(Rphi + dRphi * fRPhiSigma); // R-Phi
  240. pHit->SetError(fRPhiSigma);
  241. pHit->SetMeas(R + dR * fRSigma, 1); // R
  242. pHit->SetError(fRSigma, 1);
  243. pHit->SetDrift(driftLength);
  244. pHit->SetWirePosition(wirePos);
  245. pHit->SetTShift(pPoint->GetTime() + GetTShift(driftLength, wirePos, R, L));
  246. pHit->SetWireDelay(L);
  247. fMapOccup.insert(occupMap::value_type(WireID(uid, wirePos, R), hitID++)); // <Hash(wirePos,UID) == cellID, index> pair
  248. } // <---Loop over the DCH points
  249. // Double-track resolution
  250. Int_t counter;
  251. Double_t timeDelta, gasDriftDelta;
  252. MpdDchHit *hit1, *hit2, *tmp;
  253. typedef map<Double_t, MpdDchHit*, __ltstr> delayMap; delayMap mapDelay; // <time delay, index> pair, (invert sorting)
  254. for(occupIter It = fMapOccup.begin(); It != fMapOccup.end(); )
  255. {
  256. counter = fMapOccup.count(It->first); // hits in one cell
  257. if(fDoTest) htOccup->Fill(counter);
  258. if (counter == 1) ++It; // single hit (nothing to do)
  259. else if(counter == 2) // double hit
  260. {
  261. hit1 = (MpdDchHit*) pHitCollection->UncheckedAt(It->second); ++It; // slower hit
  262. hit2 = (MpdDchHit*) pHitCollection->UncheckedAt(It->second); ++It; // faster hit
  263. timeDelta = hit1->GetTShift() - hit2->GetTShift(); // [ns]
  264. if(timeDelta < 0.) // swap
  265. {
  266. timeDelta *= -1.;
  267. tmp = hit1;
  268. hit1 = hit2;
  269. hit2 = tmp;
  270. }
  271. if(fDoTest)
  272. {
  273. htTime->Fill(timeDelta);
  274. htPerp->Fill(wirePos = hit1->GetWirePosition());
  275. htGasDrift->Fill(gasDriftDelta = TMath::Abs(hit1->GetDrift() - hit2->GetDrift()));
  276. }
  277. if(!HitExist(timeDelta)) // overlap
  278. {
  279. hit2->AddLinks(hit1->GetLinks()); pHitCollection->Remove(hit1); // index2 faster
  280. // ->SetFlag(warning#=2); // FIXME: define errorFlags enum
  281. }
  282. else if(fDoTest)
  283. {
  284. htTimeA->Fill(timeDelta);
  285. htPerpA->Fill(wirePos);
  286. htGasDriftA->Fill(gasDriftDelta);
  287. }
  288. }
  289. else if(counter > 2) // multiple hit
  290. {
  291. // update map
  292. mapDelay.clear();
  293. for(Int_t i = 0; i < counter; i++)
  294. {
  295. hit1 = (MpdDchHit*) pHitCollection->UncheckedAt(It->second);
  296. mapDelay.insert(delayMap::value_type(hit1->GetTShift(), hit1));
  297. ++It;
  298. }
  299. // Cycle from biggest to smallest time delay
  300. for(delayMap::iterator it = mapDelay.begin(); ; )
  301. {
  302. hit1 = it->second; // slower hit
  303. ++it;
  304. if(it != mapDelay.end()) hit2 = it->second; // faster hit
  305. else break; // last pair
  306. timeDelta = hit1->GetTShift() - hit2->GetTShift();
  307. if(fDoTest)
  308. {
  309. htTime->Fill(timeDelta);
  310. htPerp->Fill(wirePos = hit1->GetWirePosition());
  311. htGasDrift->Fill(gasDriftDelta = TMath::Abs(hit1->GetDrift() - hit2->GetDrift()));
  312. }
  313. if(!HitExist(timeDelta)) // overlap, remove larger delay time
  314. {
  315. hit2->AddLinks(hit1->GetLinks());
  316. pHitCollection->Remove(hit1);
  317. // ->SetFlag(warning#=3); // FIXME: define errorFlags enum
  318. }
  319. else if(fDoTest)
  320. {
  321. htTimeA->Fill(timeDelta);
  322. htPerpA->Fill(wirePos);
  323. htGasDriftA->Fill(gasDriftDelta);
  324. }
  325. } // cycle by mapDelay
  326. } // multiple hit
  327. }
  328. pHitCollection->Compress();
  329. pHitCollection->Sort(); // in ascending order in abs(Z)
  330. cout<<" "<<pHitCollection->GetEntriesFast()<<"("<<hitID<<") hits created.\n";
  331. }
  332. //------------------------------------------------------------------------------------------------------------------------
  333. void MpdDchHitProducer::Finish()
  334. {
  335. if(fDoTest)
  336. {
  337. TFile file("test.MpdDchHitProducer.root", "RECREATE");
  338. TH1D *h1 = (TH1D*) htTimeA->Clone("htTimeEff"); h1->Divide(htTime); h1->SetYTitle("Efficiency"); h1->Write();
  339. h1 = (TH1D*) htGasDriftA->Clone("htGasDriftEff"); h1->Divide(htGasDrift); h1->SetYTitle("Efficiency"); h1->Write();
  340. h1 = (TH1D*) htPerpA->Clone("htPerpEff"); h1->Divide(htPerp); h1->SetYTitle("Efficiency"); h1->Write();
  341. fList.Write();
  342. file.Close();
  343. }
  344. }
  345. //------------------------------------------------------------------------------------------------------------------------
  346. MpdDchHit* MpdDchHitProducer::AddHit(Int_t index, Int_t detID, const TVector3& posHit, const TVector3& posHitErr,
  347. Int_t trackIndex, Int_t pointIndex, Int_t flag)
  348. {
  349. MpdDchHit *pHit = new ((*pHitCollection)[index]) MpdDchHit(detID, posHit, posHitErr, pointIndex, flag);
  350. pHit->AddLink(FairLink(1, pointIndex));
  351. pHit->AddLink(FairLink(2, trackIndex));
  352. return pHit;
  353. }
  354. //------------------------------------------------------------------------------------------------------------------------
  355. ClassImp(MpdDchHitProducer)