MpdTofHitProducerIdeal.cxx 7.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207
  1. //------------------------------------------------------------------------------------------------------------------------
  2. #include <assert.h>
  3. #include <iostream>
  4. #include <TClonesArray.h>
  5. #include "FairRootManager.h"
  6. #include "FairLogger.h"
  7. #include "MpdTofUtils.h"
  8. #include "MpdTofHit.h"
  9. #include "MpdTofPoint.h"
  10. #include "MpdTofHitProducerQA.h"
  11. #include "MpdTofHitProducerIdeal.h"
  12. using namespace std;
  13. ClassImp(MpdTofHitProducerIdeal)
  14. //------------------------------------------------------------------------------------------------------------------------
  15. MpdTofHitProducerIdeal::MpdTofHitProducerIdeal(const char *name, Bool_t useMCdata, Int_t verbose, Bool_t test, Bool_t merge, const char *flnm, bool IsEndcap)
  16. : FairTask(name, verbose), fDoTest(test), fDoMergeHits(merge), fUseMCData(useMCdata)
  17. {
  18. pQA = fDoTest ? new MpdTofHitProducerQA(flnm, IsEndcap) : nullptr;
  19. }
  20. //------------------------------------------------------------------------------------------------------------------------
  21. MpdTofHitProducerIdeal::~MpdTofHitProducerIdeal()
  22. {
  23. delete pQA;
  24. }
  25. //------------------------------------------------------------------------------------------------------------------------
  26. void MpdTofHitProducerIdeal::AddParameters(TString& buf)const
  27. {
  28. buf += "\n Run parameters:";
  29. buf += " fDoTest="; buf += fDoTest;
  30. buf += ", fDoMergeHits="; buf += fDoMergeHits;
  31. buf += ", fUseMCData="; buf += fUseMCData;
  32. buf += ", fOnlyPrimary="; buf += fOnlyPrimary;
  33. }
  34. //------------------------------------------------------------------------------------------------------------------------
  35. InitStatus MpdTofHitProducerIdeal::Init()
  36. {
  37. if(fUseMCData)
  38. {
  39. aMcPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("TOFPoint");
  40. aMcTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
  41. assert(aMcPoints);
  42. assert(aMcTracks);
  43. }
  44. else
  45. {
  46. aExpDigits = (TClonesArray*) FairRootManager::Instance()->GetObject("??????");// FIXME: NOW unknown name
  47. assert(aExpDigits);
  48. }
  49. // Create and register output array
  50. aTofHits = new TClonesArray("MpdTofHit");
  51. FairRootManager::Instance()->Register("TOFHit", "Tof", aTofHits, kTRUE);
  52. LOG(INFO)<<"[MpdTofHitProducerIdeal::Init] Initialization finished succesfully.";
  53. return kSUCCESS;
  54. }
  55. //------------------------------------------------------------------------------------------------------------------------
  56. void MpdTofHitProducerIdeal::Exec(Option_t* opt)
  57. {
  58. static const TVector3 XYZ_err(0., 0., 0.);
  59. aTofHits->Clear();
  60. Int_t nSingleHits = 0;
  61. TVector3 pos;
  62. if(fUseMCData)
  63. {
  64. for(Int_t pointIndex = 0, nTofPoint = aMcPoints->GetEntriesFast(); pointIndex < nTofPoint; pointIndex++ ) // cycle by TOF points
  65. {
  66. MpdTofPoint *pPoint = (MpdTofPoint*) aMcPoints->UncheckedAt(pointIndex);
  67. pPoint->Position(pos);
  68. AddHit(MpdTofPoint::ClearGap(pPoint->GetDetectorID()), pos, XYZ_err, pointIndex, pPoint->GetTrackID(), pPoint->GetTime(), MpdTofUtils::IsSingle);
  69. nSingleHits++;
  70. }
  71. }
  72. else // exp. data used
  73. {
  74. // FIXME: now not realized
  75. //AddHit(Int_t detUID, const TVector3 &posHit, const TVector3 &posHitErr, Int_t expDigitIndex, Double_t time, Int_t flag)
  76. assert(false);
  77. }
  78. int nFinally;
  79. if(fDoMergeHits)
  80. {
  81. MergeHitsOnStrip(); // save only the fastest hit in the strip
  82. nFinally = CompressHits(); // remove blank slotes
  83. }
  84. else nFinally = aTofHits->GetEntriesFast();
  85. LOG(DEBUG1)<<"[MpdTofHitProducerIdeal::Exec] single hits= "<<nSingleHits<<", final hits= "<<nFinally;
  86. }
  87. //------------------------------------------------------------------------------------------------------------------------
  88. void MpdTofHitProducerIdeal::Finish()
  89. {
  90. if(pQA) pQA->Finish();
  91. }
  92. //------------------------------------------------------------------------------------------------------------------------
  93. size_t MpdTofHitProducerIdeal::MergeHitsOnStrip(void)
  94. {
  95. size_t mergedNmb = 0;
  96. map<Int_t, MpdTofHit*> mHits; // pair<suid, MpdTofHit*> fastest hits map
  97. multiset<Int_t> msUIDs; // suid for Hits
  98. MpdTofHit *fastHit, *slowHit;
  99. for(Int_t hitIndex = 0, nHits = aTofHits->GetEntriesFast(); hitIndex < nHits; hitIndex++ ) // cycle by hits
  100. {
  101. auto pHit = (MpdTofHit*) aTofHits->UncheckedAt(hitIndex);
  102. assert(nullptr != pHit);
  103. Int_t suid = pHit->GetDetectorID();
  104. if(fDoTest) msUIDs.insert(suid);
  105. auto it = mHits.find(suid);
  106. if(it != mHits.end()) // hit for this suid already exist
  107. {
  108. mergedNmb++;
  109. if(pHit->GetTime() < it->second->GetTime()) // faster hit found
  110. {
  111. fastHit = pHit;
  112. slowHit = it->second;
  113. }
  114. else
  115. {
  116. fastHit = it->second;
  117. slowHit = pHit;
  118. }
  119. fastHit->AddLinks(slowHit->GetLinks()); // copy links
  120. aTofHits->Remove(slowHit); // remove old hit --> make blank slote !!
  121. fastHit->AddFlag(MpdTofUtils::HaveTail);
  122. it->second = fastHit; // change pair value to current suid
  123. if(pQA) pQA->FillMergedTimes(fastHit->GetTime(), slowHit->GetTime());
  124. }
  125. else mHits.insert(make_pair(suid, pHit)); // insert new suid pair
  126. } // cycle by hits
  127. if(pQA)
  128. {
  129. for(auto it = msUIDs.begin(), itEnd = msUIDs.end(); it != itEnd; it = msUIDs.upper_bound(*it)) pQA->FillOccupancy(msUIDs.count(*it));
  130. }
  131. return mergedNmb;
  132. }
  133. //------------------------------------------------------------------------------------------------------------------------
  134. void MpdTofHitProducerIdeal::AddHit(Int_t suid, const TVector3 &posHit, const TVector3 &posHitErr, Int_t expDigitIndex, Double_t time, Int_t flag)
  135. {
  136. MpdTofHit *pHit = new ((*aTofHits)[aTofHits->GetEntriesFast()]) MpdTofHit(suid, posHit, posHitErr, expDigitIndex, time, flag);
  137. pHit->AddLink(FairLink(MpdTofUtils::expDigitIndex, expDigitIndex));
  138. pHit->AddLink(FairLink(MpdTofUtils::volumeUID, suid));
  139. }
  140. //------------------------------------------------------------------------------------------------------------------------
  141. void MpdTofHitProducerIdeal::AddHit(Int_t suid, const TVector3 &posHit, const TVector3 &posHitErr, Int_t mcPointIndex, Int_t mcTrackIndex, Double_t time, Int_t flag)
  142. {
  143. MpdTofHit *pHit = new ((*aTofHits)[aTofHits->GetEntriesFast()]) MpdTofHit(suid, posHit, posHitErr, mcPointIndex, time, flag);
  144. pHit->AddLink(FairLink(MpdTofUtils::mcPointIndex, mcPointIndex));
  145. pHit->AddLink(FairLink(MpdTofUtils::mcTrackIndex, mcTrackIndex));
  146. pHit->AddLink(FairLink(MpdTofUtils::volumeUID, suid));
  147. }
  148. //------------------------------------------------------------------------------------------------------------------------
  149. Int_t MpdTofHitProducerIdeal::CompressHits()
  150. {
  151. aTofHits->Compress();
  152. return aTofHits->GetEntriesFast();
  153. }
  154. //------------------------------------------------------------------------------------------------------------------------
  155. void MpdTofHitProducerIdeal::Dump(const char* title, ostream& out) const
  156. {
  157. out<<"\n [MpdTofHitProducer::Dump] "; if(title) out<<title; out<<", size= "<<aTofHits->GetEntriesFast();
  158. MpdTofPoint *point; MpdTofHit *pHit; TVector3 hitPos, pointPos;
  159. TIterator *iter = aTofHits->MakeIterator();
  160. while( (pHit = (MpdTofHit*) iter->Next()) )
  161. {
  162. pHit->Position(hitPos);
  163. out<<"\n hit suid = "<<pHit->GetDetectorID()<<", hit pos("<<hitPos.X()<<","<<hitPos.Y()<<","<<hitPos.Z()<<"), flag ="<<pHit->GetFlag();
  164. if(aMcPoints)
  165. {
  166. point = (MpdTofPoint*) aMcPoints->UncheckedAt(pHit->GetRefIndex());
  167. point->Position(pointPos);
  168. out<<"\n point suid = "<<point->GetDetectorID()<<", point pos("<<pointPos.X()<<","<<pointPos.Y()<<","<<pointPos.Z()<<"), dev="<<(hitPos-pointPos).Mag();
  169. }
  170. }
  171. delete iter;
  172. }
  173. //------------------------------------------------------------------------------------------------------------------------