MpdVertexZfinder.cxx 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216
  1. // -------------------------------------------------------------------------
  2. // ----- MpdVertexZfinder source file -----
  3. // ----- Created 19/08/09 by A. Zinchenko -----
  4. // -------------------------------------------------------------------------
  5. /** MpdVertexZfinder.h
  6. *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
  7. **
  8. ** Primary vertex Z-position evaluator in MPD using hits from TPC
  9. **/
  10. #include "MpdVertexZfinder.h"
  11. #include "MpdKalmanFilter.h"
  12. #include "MpdKalmanGeoScheme.h"
  13. #include "MpdKalmanHit.h"
  14. #include "MpdTpcSectorGeo.h"
  15. //#include "TpcPadPlane.h"
  16. #include <TClonesArray.h>
  17. #include <TMath.h>
  18. #include <TPad.h>
  19. #include <TVector3.h>
  20. #include <iostream>
  21. using std::cout;
  22. using std::endl;
  23. //const Double_t MpdTrackFinderIts::fgkChi2Cut = 20; //20; //100;
  24. //__________________________________________________________________________
  25. MpdVertexZfinder::MpdVertexZfinder(const char *name, Int_t iVerbose)
  26. :FairTask(name, iVerbose)
  27. {
  28. fKHits = nullptr, fhLays = nullptr;
  29. fhZ = nullptr, fUnc = nullptr;
  30. }
  31. //__________________________________________________________________________
  32. MpdVertexZfinder::~MpdVertexZfinder()
  33. {
  34. //delete fKHits;
  35. //delete fTracks;
  36. //delete fTrackCand;
  37. delete fhZ;
  38. delete fUnc;
  39. }
  40. //__________________________________________________________________________
  41. InitStatus MpdVertexZfinder::Init()
  42. {
  43. //AZ fhZ = new TH1F("hZfinder","Zv",250,-50.2,49.8);
  44. fhZ = new TH1F("hZfinder","Zv",500,-100.2,99.8);
  45. fUnc = new TF1("fUnc","[0]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])+[3]",-10,10);
  46. //return ReInit();
  47. ReInit();
  48. return kSUCCESS;
  49. }
  50. //__________________________________________________________________________
  51. InitStatus MpdVertexZfinder::ReInit()
  52. {
  53. return kSUCCESS;
  54. }
  55. //__________________________________________________________________________
  56. void MpdVertexZfinder::Reset()
  57. {
  58. ///
  59. cout << " MpdVertexZfinder::Reset " << endl;
  60. fhZ->Reset();
  61. fhZ->GetXaxis()->SetRange(10,0); // reset axis range
  62. }
  63. //__________________________________________________________________________
  64. void MpdVertexZfinder::SetParContainers()
  65. {
  66. // Get run and runtime database
  67. /*
  68. FairRunAna* run = FairRunAna::Instance();
  69. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  70. FairRuntimeDb* db = run->GetRuntimeDb();
  71. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  72. // Get STS geometry parameter container
  73. db->getContainer("MpdStsGeoPar");
  74. */
  75. }
  76. //__________________________________________________________________________
  77. void MpdVertexZfinder::Finish()
  78. {
  79. //Write();
  80. }
  81. //__________________________________________________________________________
  82. void MpdVertexZfinder::Exec(Option_t * option)
  83. {
  84. static int eventCounter = 0;
  85. cout << " - - - - \n Vertex Z-finder event " << ++eventCounter << endl;
  86. Reset();
  87. }
  88. //__________________________________________________________________________
  89. void MpdVertexZfinder::SetHits(const TClonesArray *hits, const TH1F *hLays)
  90. {
  91. /// Set hits container and histogram
  92. if (fKHits) return;
  93. fKHits = hits;
  94. fhLays = hLays;
  95. }
  96. //__________________________________________________________________________
  97. Double_t MpdVertexZfinder::FindZ(const Int_t *layPointers, Int_t &flag)
  98. {
  99. /// Evaluate vertex Z
  100. //Int_t layMax0 = ((MpdKalmanHit*)fKHits->First())->GetLayer();
  101. Int_t modular = 0;
  102. if (((MpdKalmanHit*)fKHits->First())->GetType() == MpdKalmanHit::kFixedP) modular = 1;
  103. //const TpcPadPlane *padPlane = TpcPadPlane::Instance();
  104. MpdTpcSectorGeo *secGeo = MpdTpcSectorGeo::Instance();
  105. // Estimate Z-position of vertex
  106. // Loop over layers
  107. //Int_t layBeg = layMax0, layEnd = layBeg - 2, iDir = -1, dLays = 6; //10;
  108. Int_t layBeg = ((MpdKalmanHit*)fKHits->Last())->GetLayer(), layEnd = layBeg + 2;
  109. Int_t iDir = 1, dLays = 5, isec = 0, isec1 = 0;
  110. //Int_t layMax = layMax0, dLays = 6; //10;
  111. TVector3 pos1, pos2, posLoc;
  112. Double_t rad1, phi1, rad2, phi2;
  113. for (Int_t lay = layBeg; lay != layEnd; lay+=iDir) {
  114. //for (Int_t lay = layMax; lay > layMax-2; --lay) {
  115. Int_t nHits1 = (Int_t) fhLays->GetBinContent(lay+1,0);
  116. Int_t nHits2 = (Int_t) fhLays->GetBinContent(lay+dLays*iDir+1,0);
  117. //cout << "Hits: " << nHits1 << " " << nHits2 << endl;
  118. // Loop over hits in first layer
  119. for (Int_t ihit1 = 0; ihit1 < nHits1; ++ihit1) {
  120. MpdKalmanHit *hit1 = (MpdKalmanHit*) fKHits->UncheckedAt(layPointers[lay]+ihit1);
  121. if (hit1->GetFlag() & MpdKalmanHit::kUsed) continue;
  122. if (modular) {
  123. posLoc.SetXYZ(hit1->GetMeas(0), hit1->GetDist(), hit1->GetMeas(1));
  124. isec = hit1->GetDetectorID() % 1000000;
  125. //pos1 = padPlane->fromSectorReferenceFrame(posLoc, isec);
  126. secGeo->Local2Global(isec, posLoc, pos1);
  127. rad1 = pos1.Pt();
  128. phi1 = pos1.Phi();
  129. } else {
  130. rad1 = hit1->GetPos();
  131. phi1 = hit1->GetMeas(0) / rad1;
  132. }
  133. // Loop over hits in second (closer to the beam) layer
  134. for (Int_t ihit2 = 0; ihit2 < nHits2; ++ihit2) {
  135. MpdKalmanHit *hit2 = (MpdKalmanHit*) fKHits->UncheckedAt(layPointers[lay+dLays*iDir]+ihit2);
  136. //if (h2->GetTrackID() != h1->GetTrackID()) continue; //!!! for test - exact ID matching
  137. if (hit2->GetFlag() & MpdKalmanHit::kUsed) continue;
  138. if (modular) {
  139. posLoc.SetXYZ(hit2->GetMeas(0), hit2->GetDist(), hit2->GetMeas(1));
  140. isec1 = hit2->GetDetectorID() % 1000000;
  141. //if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < padPlane->nSectors()-1) continue;
  142. //pos2 = padPlane->fromSectorReferenceFrame(posLoc, isec1);
  143. if (TMath::Abs(isec-isec1) > 1 && TMath::Abs(isec-isec1) < secGeo->NofSectors()-1) continue;
  144. secGeo->Local2Global(isec1, posLoc, pos2);
  145. rad2 = pos2.Pt();
  146. //if ((rad2 - rad1) * iDir < 0.5) continue;
  147. phi2 = pos2.Phi();
  148. } else {
  149. rad2 = hit2->GetPos();
  150. //if (TMath::Abs(rad1 - rad2) < 0.5) continue;
  151. phi2 = hit2->GetMeas(0) / rad2;
  152. }
  153. Double_t dPhi = MpdKalmanFilter::Instance()->Proxim(phi1,phi2) - phi1;
  154. if (TMath::Abs(dPhi) > TMath::PiOver4()) continue;
  155. //Double_t zvert = pos2.Z() - (pos1.Z()-pos2.Z()) / (rad1-rad2) * rad2;
  156. Double_t zvert = hit2->GetMeas(1) -
  157. (hit1->GetMeas(1)-hit2->GetMeas(1)) / (rad1-rad2) * rad2;
  158. fhZ->Fill(zvert);
  159. }
  160. }
  161. }
  162. Double_t z = fhZ->GetXaxis()->GetBinCenter(fhZ->GetMaximumBin());
  163. Double_t max = fhZ->GetMaximum();
  164. fUnc->SetParameters(max,z,5.0,1.0);
  165. fUnc->SetParLimits(0,max*0.2,max*9.0);
  166. fUnc->SetParLimits(1,z-5.,z+5.);
  167. fUnc->SetParLimits(2,0.2,10.);
  168. //fhZ->Fit(fUnc,"ww","",z-10.,z+10.);
  169. fhZ->Fit(fUnc,"wwNQ","",z-10.,z+10.);
  170. Double_t zFit = fUnc->GetParameter(1);
  171. Double_t zVert = (TMath::Abs(z-zFit) > 2.0) ? z : zFit;
  172. // Set quality flag
  173. if (max < 6 || (max < 10 && TMath::Abs(z-zFit) > 1)) { flag = -1; zVert = 0.0; }
  174. else flag = 0;
  175. cout << " MpdVertexZfinder::FindZ: " << z << " " << zFit << " " << zVert << " " << max << " " << flag << endl;
  176. /*
  177. fhZ->Draw();
  178. gPad->Update();
  179. Char_t go1[1];
  180. gets(go1);
  181. */
  182. return zVert;
  183. }
  184. ClassImp(MpdVertexZfinder);