MpdEmcMatchingKI.cxx 7.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257
  1. //--------------------------------------------------------------------
  2. //
  3. // Description:
  4. // MPD TPC-EMC Matching
  5. //
  6. //
  7. // Author List:
  8. // D.Peresunko RRCKI 2019
  9. //
  10. //--------------------------------------------------------------------
  11. #include "MpdEmcMatchingKI.h"
  12. #include "MpdEmcClusterKI.h"
  13. #include "MpdEmcGeoUtils.h"
  14. #include "MpdEmcTrackExtrap.h"
  15. #include "MpdKalmanFilter.h"
  16. #include "MpdTpcHit.h"
  17. #include "MpdTpcKalmanTrack.h"
  18. #include "MpdVertex.h"
  19. #include "FairLogger.h"
  20. #include "MpdMCTrack.h"
  21. #include "FairRootManager.h"
  22. #include "FairRunAna.h"
  23. #include "MpdEmcSimParams.h"
  24. #include <TClonesArray.h>
  25. #include <TGeoManager.h>
  26. #include <TMath.h>
  27. #include <TVector3.h>
  28. #include <iterator>
  29. // ----- Default constructor -------------------------------------------
  30. MpdEmcMatchingKI::MpdEmcMatchingKI() : FairTask("TPC-EMC matching") {}
  31. // ----- Destructor ----------------------------------------------------
  32. MpdEmcMatchingKI::~MpdEmcMatchingKI() {}
  33. // -------------------------------------------------------------------------
  34. // ----- Public method Init --------------------------------------------
  35. InitStatus MpdEmcMatchingKI::Init()
  36. {
  37. // Get RootManager
  38. FairRootManager* ioman = FairRootManager::Instance();
  39. if (!ioman) {
  40. LOG(ERROR) << "MpdEmcMatchingKI::Init: "
  41. << "RootManager not instantiated!" << endl;
  42. return kFATAL;
  43. }
  44. // Get input array
  45. fClusterArray = (TObjArray*)ioman->GetObject("EmcCluster");
  46. if (!fClusterArray) {
  47. LOG(ERROR) << "MpdEmcMatchingKI::Init: "
  48. << "No EMC clusters array!" << endl;
  49. return kERROR;
  50. }
  51. fTpcTracks = (TClonesArray*)ioman->GetObject("TpcKalmanTrack");
  52. if (!fTpcTracks) {
  53. LOG(ERROR) << "MpdEmcMatchingKI::Init: "
  54. << "No TPC track array!" << endl;
  55. return kERROR;
  56. }
  57. fvtx = (TClonesArray*)ioman->GetObject("Vertex");
  58. if (!fvtx) {
  59. LOG(INFO) << "Vertex not found, assume (0,0,0)" << endl;
  60. }
  61. double rMax;
  62. MpdEmcGeoUtils* geom = MpdEmcGeoUtils::GetInstance();
  63. geom->GetECALTubeSize(frMin, rMax, fzMax);
  64. return kSUCCESS;
  65. }
  66. //__________________________________________________________________________
  67. void MpdEmcMatchingKI::Finish() {}
  68. //__________________________________________________________________________
  69. void MpdEmcMatchingKI::Exec(Option_t* opt)
  70. {
  71. // Main processing engine
  72. // Prepare list of cluster coordinates
  73. // Prepare list of coordinates of tracks extrapolated to inner ECAL radius
  74. // For each cluster find N closest tracks and extrapolate them to the 1/2 depth of ECAL - look if distance decrease
  75. // For each track find M closest clusters and extrapolate to 1/2 depth of ECAL - check if distance decrease
  76. // Calculate distance along track and perpendicular to track and (dz,dphi)
  77. Int_t nClu = fClusterArray->GetEntriesFast();
  78. Int_t ntpc = fTpcTracks->GetEntriesFast();
  79. LOG(INFO) << "MpdEmcMatchingKI::Exec started: number of clusters " << nClu << ", number of tracks " << ntpc << endl;
  80. if (ntpc == 0 || nClu == 0) {
  81. return;
  82. }
  83. CorrectClustersVtx();
  84. ExtrapolateTracks();
  85. MakeMatchToClusters();
  86. // MakeMatchToTracks() ;
  87. }
  88. //__________________________________________________________________________
  89. void MpdEmcMatchingKI::CorrectClustersVtx()
  90. {
  91. // Get vertex
  92. if (!fvtx) { // no vertex, no corrections
  93. return;
  94. }
  95. MpdEmcSimParams* simParams = MpdEmcSimParams::GetInstance();
  96. TVector3 primVert;
  97. MpdVertex* vtx = (MpdVertex*)fvtx->First();
  98. vtx->Position(primVert);
  99. if (primVert.Z() == 0. || !simParams->CorrectZ()) { // no need to correct
  100. return;
  101. }
  102. Int_t nClu = fClusterArray->GetEntriesFast();
  103. for (int i = 0; i < nClu; i++) {
  104. MpdEmcClusterKI* clu = static_cast<MpdEmcClusterKI*>(fClusterArray->UncheckedAt(i));
  105. clu->CorrectVertex(primVert.Z());
  106. }
  107. }
  108. //__________________________________________________________________________
  109. void MpdEmcMatchingKI::ExtrapolateTracks()
  110. {
  111. // Clean containers
  112. // Sort Cluster extrapolations according to grid
  113. fTrackPoints.clear();
  114. MpdKalmanFilter* pKF = MpdKalmanFilter::Instance("KF", "KF");
  115. MpdKalmanHit hEnd;
  116. hEnd.SetType(MpdKalmanHit::kFixedR);
  117. int nTracks = fTpcTracks->GetEntriesFast();
  118. for (int itrack = 0; itrack < nTracks; itrack++) {
  119. MpdTpcKalmanTrack* tr = (MpdTpcKalmanTrack*)fTpcTracks->UncheckedAt(itrack);
  120. if (TMath::Abs((*tr->GetParamAtHit())(1, 0)) > fzMax)
  121. continue;
  122. MpdTpcKalmanTrack tr1(*tr);
  123. tr1.SetParam(*tr1.GetParamAtHit());
  124. tr1.SetParamNew(*tr1.GetParamAtHit());
  125. tr1.SetWeight(*tr1.GetWeightAtHit());
  126. tr1.SetPos(tr1.GetPosAtHit());
  127. tr1.SetPosNew(tr1.GetPos());
  128. tr1.SetLength(tr1.GetLengAtHit());
  129. // Propagate to EMC inner radius
  130. hEnd.SetPos(frMin);
  131. if (!pKF->PropagateToHit(&tr1, &hEnd, kTRUE)) {
  132. continue;
  133. }
  134. double phi = tr1.GetParamNew(0) / tr1.GetPosNew();
  135. double z = tr1.GetParamNew(1);
  136. double r = tr1.GetPosNew();
  137. // Remember position at ECAL surface
  138. double surfX = r * TMath::Cos(phi);
  139. double surfY = r * TMath::Sin(phi);
  140. double surfZ = z;
  141. // Extrapolate further to the reconstruction radius of clusters
  142. double rClu = MpdEmcGeoUtils::GetInstance()->Rperp(z) ;
  143. hEnd.SetPos(rClu);
  144. bool iok = pKF->PropagateToHit(&tr1, &hEnd, kTRUE);
  145. double depthX, depthY, depthZ;
  146. if (iok) {
  147. r = tr1.GetPosNew();
  148. phi = tr1.GetParamNew(0) / r;
  149. depthX = r * TMath::Cos(phi);
  150. depthY = r * TMath::Sin(phi);
  151. depthZ = tr1.GetParamNew(1);
  152. } else { // old coordinates
  153. depthX = surfX;
  154. depthY = surfY;
  155. depthZ = surfZ;
  156. }
  157. fTrackPoints.emplace_back(itrack, surfX, surfY, surfZ, depthX, depthY, depthZ);
  158. }
  159. }
  160. //__________________________________________________________________________
  161. void MpdEmcMatchingKI::MakeMatchToClusters()
  162. {
  163. // Look over track extrapolations and find closest
  164. Int_t nClu = fClusterArray->GetEntriesFast();
  165. for (int i = 0; i < nClu; i++) {
  166. MpdEmcClusterKI* clu = static_cast<MpdEmcClusterKI*>(fClusterArray->UncheckedAt(i));
  167. std::vector<MpdEmcTrackExtrap>::iterator minTr;
  168. double minD = 9999.;
  169. std::vector<MpdEmcTrackExtrap>::iterator a = fTrackPoints.begin();
  170. while (a != fTrackPoints.end()) {
  171. double distance = a->Distance(clu);
  172. if (distance < minD) {
  173. minD = distance;
  174. minTr = a;
  175. }
  176. a++;
  177. }
  178. // Evaluate distances in (phi,z) (alongTr, perpTr)
  179. if (minD < 999) { // found some tracks
  180. double dphi, dz;
  181. minTr->DistanceDphiDz(clu, dphi, dz);
  182. clu->SetTrackIndex(minTr->GetTrackIndex());
  183. clu->SetTrackDxDz(dphi, dz);
  184. } else {
  185. // Distance to track extrapolation
  186. // Track index already set to -1 at initialization
  187. clu->SetTrackDxDz(9999., 9999.);
  188. }
  189. }
  190. }
  191. //__________________________________________________________________________
  192. void MpdEmcMatchingKI::MakeMatchToTracks()
  193. {
  194. // Look over track extrapolations and find closest
  195. Int_t nClu = fClusterArray->GetEntriesFast();
  196. std::vector<MpdEmcTrackExtrap>::iterator a = fTrackPoints.begin();
  197. while (a != fTrackPoints.end()) {
  198. double minD = 9999.;
  199. //int icluMin;
  200. for (int i = 0; i < nClu; i++) {
  201. MpdEmcClusterKI* clu = static_cast<MpdEmcClusterKI*>(fClusterArray->UncheckedAt(i));
  202. double distance = a->Distance(clu);
  203. if (distance < minD) {
  204. minD = distance;
  205. //icluMin = i;
  206. }
  207. }
  208. // Evaluate distances in (phi,z) (alongTr, perpTr)
  209. /* if(minD<999){ //found some cluster
  210. a->SetDistance(icluMin) ;
  211. }
  212. */
  213. }
  214. }
  215. ClassImp(MpdEmcMatchingKI)