MpdEmcMatching.cxx 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440
  1. //--------------------------------------------------------------------
  2. //
  3. // Description:
  4. // MPD TPC-EMC Matching
  5. //
  6. //
  7. // Author List:
  8. // Alexander Zinchenko LHEP, JINR, Dubna - 1-June-2016
  9. // Alexander Zinchenko LHEP, JINR, Dubna - 8-June-2018 - adapted for projective geometry
  10. //
  11. //--------------------------------------------------------------------
  12. #include "MpdEmcMatching.h"
  13. #include "MpdEmcGeoParams.h"
  14. //#include "MpdEmcDigit.h"
  15. #include "MpdEmcMatch.h"
  16. #include "MpdEmcPoint.h"
  17. #include "MpdKalmanFilter.h"
  18. #include "MpdTpcHit.h"
  19. #include "MpdTpcKalmanTrack.h"
  20. #include "FairMCPoint.h"
  21. #include "MpdMCTrack.h"
  22. #include "FairRootManager.h"
  23. #include "FairRunAna.h"
  24. #include <TClonesArray.h>
  25. #include <TGeoManager.h>
  26. #include <TMath.h>
  27. #include <TSpline.h>
  28. #include <iterator>
  29. using namespace std;
  30. using namespace TMath;
  31. //FILE *lun = fopen ("file.txt","w"); //AZ
  32. // ----- Default constructor -------------------------------------------
  33. MpdEmcMatching::MpdEmcMatching() :
  34. FairTask("TPC-EMC matching") {
  35. }
  36. // ----- Destructor ----------------------------------------------------
  37. MpdEmcMatching::~MpdEmcMatching() {
  38. }
  39. // -------------------------------------------------------------------------
  40. // ----- Public method Init --------------------------------------------
  41. InitStatus MpdEmcMatching::Init() {
  42. cout << "******************* EMC Matching INIT *********************" << endl;
  43. // Get RootManager
  44. FairRootManager* ioman = FairRootManager::Instance();
  45. if (!ioman) {
  46. cout << "-E- MpdEmcMatching::Init: "
  47. << "RootManager not instantiated!" << endl;
  48. return kFATAL;
  49. }
  50. fGeoPar = new MpdEmcGeoParams();
  51. // Containers for rec. points from all+2 sectors
  52. Int_t nSec = fGeoPar->GetNsec() / 2; // number of long sectors
  53. multimap<Double_t,Int_t> aaa;
  54. fRecPoints.assign(nSec+2,aaa);
  55. // Fill ending rows of sectors
  56. const Int_t nSecRows = fGeoPar->GetNrows() / (nSec-1); // rows per wide sector
  57. fSecRows0.insert(nSecRows);
  58. for (Int_t isec = 1; isec < nSec; ++isec) {
  59. if (isec == 2 || isec == 6) fSecRows0.insert(*fSecRows0.rbegin()+nSecRows/2);
  60. else fSecRows0.insert(*fSecRows0.rbegin()+nSecRows);
  61. }
  62. // Get input array
  63. fPointArray = (TClonesArray*) ioman->GetObject("EmcPoint");
  64. if (!fPointArray) {
  65. cout << "-W- MpdEmcMatching::Init: " << "No EmcPoint array!" << endl;
  66. return kERROR;
  67. }
  68. fMcTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
  69. if (!fMcTrackArray) {
  70. cout << "-W- MpdEmcMatching::Init: " << "No MCTrack array!" << endl;
  71. return kERROR;
  72. }
  73. fRecPointArray = (TClonesArray*) ioman->GetObject("EmcRecPoint");
  74. if (!fRecPointArray) {
  75. cout << "-W- MpdEmcMatching::Init: " << "No RecPoint array!" << endl;
  76. return kERROR;
  77. }
  78. fTpcTracks = (TClonesArray*) ioman->GetObject("TpcKalmanTrack");
  79. if (!fTpcTracks) {
  80. cout << "-W- MpdEmcMatching::Init: " << "No TPC track array!" << endl;
  81. return kERROR;
  82. }
  83. // Create and register output array
  84. //*
  85. fMatchArray = new TClonesArray("MpdEmcMatch", 100);
  86. ioman->Register("EmcMatch", "EMC", fMatchArray, kTRUE);
  87. //*/
  88. cout << "-I- MpdEmcMatching: Intialization successfull" << endl;
  89. return kSUCCESS;
  90. }
  91. //__________________________________________________________________________
  92. void MpdEmcMatching::Finish() {
  93. cout << "-I- MpdEmcMatching: Finish" << endl;
  94. }
  95. //__________________________________________________________________________
  96. void MpdEmcMatching::Exec(Option_t* opt)
  97. {
  98. // Main processing engine
  99. cout << "MpdEmcMatching::Exec started" << endl;
  100. static const Int_t nSec = fGeoPar->GetNsec() / 2; // number of EMC long sectors
  101. // Reset output Array
  102. //if (!fDigiArray) Fatal("MpdEmcMatching::Exec", "No array of digits");
  103. fMatchArray->Delete();
  104. // Clear rec. point containers
  105. for (Int_t i = 0; i < nSec+2; ++i) fRecPoints[i].clear();
  106. // Fill rec. point containers
  107. Int_t nRecP = fRecPointArray->GetEntriesFast();
  108. cout << " Total number of rec. points: " << nRecP << endl;
  109. for (Int_t i = 0; i < nRecP; ++i) {
  110. MpdTpcHit *recp = (MpdTpcHit*) fRecPointArray->UncheckedAt(i);
  111. Int_t isec = recp->GetLayer();
  112. //fRecPoints[isec].insert(pair<Double_t,Int_t>(recp->GetZ(),i)); // sort according to theta
  113. fRecPoints[isec].insert(pair<Double_t,Int_t>(recp->GetLocalZ(),i)); // sort according to theta bins
  114. }
  115. Int_t ntpc = fTpcTracks->GetEntriesFast();
  116. if (ntpc == 0) return;
  117. for (Int_t i = 0; i < ntpc; ++i) {
  118. //if (isec < nSec && (*tr->GetParamAtHit())(1,0) > 20 && (*tr->GetParamAtHit())(3,0) > 0) continue; // different sides
  119. //else if (isec >= nSec && (*tr->GetParamAtHit())(1,0) < -20 && (*tr->GetParamAtHit())(3,0) < 0) continue; // different sides
  120. DoMatching(i);
  121. }
  122. }
  123. //__________________________________________________________________________
  124. void MpdEmcMatching::DoMatching(Int_t itrack)
  125. {
  126. // Match track with EMC clusters
  127. static const Int_t nSec = fGeoPar->GetNsec() / 2; // number of EMC long sectors
  128. static const Double_t rMin = fGeoPar->GetRmin(); // inner radius !!!
  129. static const Double_t zMax = fGeoPar->GetLength();
  130. static Int_t first = 1;
  131. static TSpline3 *rMaxS;
  132. if (first) {
  133. // Get rmax-vs-theta dependence
  134. first = 0;
  135. const vector<Double_t> &thes = fGeoPar->GetThetaBox();
  136. const vector<Double_t> &rhos = fGeoPar->GetRhoCenterBox();
  137. const vector<Double_t> &zs = fGeoPar->GetZCenterBox();
  138. Int_t nthe = thes.size();
  139. Double_t *the = new Double_t [nthe];
  140. Double_t *rmax = new Double_t [nthe];
  141. Double_t height = fGeoPar->GetLengthBox(); // tower half-height
  142. Int_t j1 = -1;
  143. for (Int_t j = nthe-1; j >= 0; --j) {
  144. Double_t rho = rhos[j];
  145. Double_t z = zs[j];
  146. Double_t theta1 = thes[j];
  147. if (j < nthe-1 && thes[j] <= thes[j+1]+0.1) theta1 = 180 - theta1;
  148. Double_t costhe = TMath::Cos(theta1*TMath::DegToRad());
  149. Double_t sinthe = TMath::Sin(theta1*TMath::DegToRad());
  150. rho += height * sinthe;
  151. z += height * costhe;
  152. the[++j1] = TMath::ATan2(rho,z) * TMath::RadToDeg();
  153. rmax[j1] = rho;
  154. }
  155. rMaxS = new TSpline3("rmax",the,rmax,nthe); // rmax vs theta
  156. delete [] the;
  157. delete [] rmax;
  158. }
  159. MpdTpcKalmanTrack *tr = (MpdTpcKalmanTrack*) fTpcTracks->UncheckedAt(itrack);
  160. if (TMath::Abs((*tr->GetParamAtHit())(1,0)) > zMax) return;
  161. MpdTpcKalmanTrack tr1(*tr);
  162. tr1.SetParam(*tr1.GetParamAtHit());
  163. tr1.SetParamNew(*tr1.GetParamAtHit());
  164. tr1.SetWeight(*tr1.GetWeightAtHit());
  165. tr1.SetPos(tr1.GetPosAtHit());
  166. tr1.SetPosNew(tr1.GetPos());
  167. tr1.SetLength(tr1.GetLengAtHit());
  168. //Double_t eta = TMath::Abs (tr->Momentum3().Eta());
  169. Double_t theta = tr->Theta() * TMath::RadToDeg();
  170. Double_t pt = tr->Pt();
  171. // For chi2: parameters obtained for box e (pT = 0.2-1.9, eta = -1.1-1.1)
  172. //Double_t depth = 2.3 + 8.8 * eta;
  173. //Double_t sigmaL = 0.82 + 1.50 * eta + 0.91 * eta * eta;
  174. //Double_t sigmaL2 = sigmaL * sigmaL, sigmaT2 = 0.66 * 0.66;
  175. Double_t sigTel[3] = {1.004, 0.0, 1.386};
  176. sigTel[1] = 0.051 / TMath::Exp(TMath::Log(pt*pt+1.494e-4*1.494e-4)*0.616) + 0.115;
  177. Double_t sigLel[4] = {0.789, 0.0, 1.785, 0.528};
  178. Double_t depth[4] = {0,0,0,-2.0};
  179. depth[1] = 0.387 / TMath::Exp(TMath::Log(pt*pt+0.382*0.382)*0.737) + 0.210;
  180. sigLel[1] = 0.064 / TMath::Exp(TMath::Log(pt*pt+0.648*0.648)*2.426) + 0.146;
  181. depth[0] = -6.080 / TMath::Exp(TMath::Log(pt*pt+0.813*0.813)*2.588) - 2.830;
  182. // For chi2pi: parameters obtained for box pi (pT = 0.1-2.0, eta = -1.1-1.1)
  183. //Double_t sig1 = 7.8 - 10.9 * eta + 8.6 * eta * eta;
  184. //Double_t sig2 = 7.1 - 8.1 * eta + 5.3 * eta * eta;
  185. //Double_t et = TMath::Min (eta, 1.0);
  186. //Double_t sigT = 1.4 - 2.2 * et + 5.6 * et * et - 3.2 * et * et * et;
  187. Double_t sigTpi[3] = {1.222, 0.350, 1.337};
  188. Double_t sigLpi[3] = {3.21, 0.0, 3.09};
  189. Double_t depthPi = 0.611 / TMath::Exp(TMath::Log(pt*pt+0.651*0.651)*1.278) + 0.307;
  190. sigLpi[1] = 0.197 / TMath::Exp(TMath::Log(pt*pt+0.202*0.202)*0.859) + 0.126;
  191. Int_t secs[2], iok = 0;
  192. Double_t thmin = 999, thmax = -999, thetr[2], phitr[2];
  193. pair<multimap<Double_t,Int_t>::iterator,multimap<Double_t,Int_t>::iterator> pits[2]; // up to 2 sectors to check
  194. // Loop over TPC tracks
  195. MpdKalmanFilter *pKF = MpdKalmanFilter::Instance("KF","KF");
  196. MpdKalmanHit hEnd;
  197. hEnd.SetType(MpdKalmanHit::kFixedR);
  198. // Propagate to EMC min. radius
  199. hEnd.SetPos(rMin);
  200. iok = pKF->PropagateToHit(&tr1,&hEnd,kTRUE);
  201. if (!iok) return;
  202. Double_t phi = tr1.GetParamNew(0) / tr1.GetPosNew();
  203. //if (phi < 0) phi += TMath::TwoPi();
  204. //else if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
  205. //phitr[0] = phi;
  206. TVector3 pos (tr1.GetPosNew()*TMath::Cos(phi), tr1.GetPosNew()*TMath::Sin(phi), tr1.GetParamNew(1));
  207. GetTowerCoords(pos,0,phitr[0],thetr[0]);
  208. thmin = TMath::Min (thmin,thetr[0]);
  209. thmax = TMath::Max (thmax,thetr[0]);
  210. // Propagate to EMC max. radius by small steps
  211. Double_t r = rMin, rmax = rMaxS->Eval(theta), dr = 4.0, z = 0;
  212. while (r < rmax) {
  213. hEnd.SetPos(r);
  214. iok = pKF->PropagateToHit(&tr1,&hEnd,kTRUE);
  215. if (!iok) break;
  216. //if (r < rmin + 1 && TMath::Abs(theT) > zmax) break;
  217. r += dr;
  218. }
  219. if (iok) {
  220. z = tr1.GetParamNew(1);
  221. r = tr1.GetPosNew();
  222. phi = tr1.GetParamNew(0) / r;
  223. } else {
  224. z = tr1.GetParam(1);
  225. r = tr1.GetPos();
  226. phi = tr1.GetParam(0) / r;
  227. }
  228. //phitr[1] = phi;
  229. pos.SetXYZ(r*TMath::Cos(phi),r*TMath::Sin(phi),z);
  230. GetTowerCoords(pos,1,phitr[1],thetr[1]);
  231. thmin = TMath::Min (thmin,thetr[1]);
  232. thmax = TMath::Max (thmax,thetr[1]);
  233. thmin -= 2;
  234. thmax += 2;
  235. for (Int_t io = 0; io < 2; ++io) {
  236. set<Int_t>::iterator sit = fSecRows0.upper_bound(phitr[io]);
  237. secs[io] = std::distance(fSecRows0.begin(),sit);
  238. }
  239. Int_t dsecs = secs[1] - secs[0];
  240. if (TMath::Abs(dsecs) > 1 && TMath::Abs(dsecs) < nSec / 2) secs[1] = secs[0] + TMath::Sign(1,dsecs);
  241. else if (TMath::Abs(dsecs) > nSec / 2) secs[1] = secs[0] - TMath::Sign(1,dsecs);
  242. if (secs[1] < 0) secs[1] += nSec;
  243. /*
  244. for (Int_t io = 0; io < 2; ++io) {
  245. Double_t r = (io == 0) ? rMin : rMaxS->Eval(theta);
  246. hEnd.SetPos(r);
  247. iok = pKF->PropagateToHit(&tr1,&hEnd,kTRUE);
  248. if (!iok) break;
  249. ztr[io] = tr1.GetParamNew(1);
  250. Double_t phi = tr1.GetParamNew(0) / tr1.GetPosNew();
  251. if (phi < 0) phi += TMath::TwoPi();
  252. else if (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
  253. phitr[io] = phi;
  254. zminmax[0] = TMath::Min (zminmax[0],ztr[io]);
  255. zminmax[1] = TMath::Max (zminmax[1],ztr[io]);
  256. }
  257. if (!iok) return;
  258. zminmax[0] -= 2 * dzTower;
  259. zminmax[1] += 2 * dzTower;
  260. secs[0] = phitr[0] / dphiSec;
  261. secs[1] = phitr[1] / dphiSec;
  262. Int_t dsecs = secs[1] - secs[0];
  263. if (TMath::Abs(dsecs) > 1 && TMath::Abs(dsecs) < nSec / 2) secs[1] = secs[0] + TMath::Sign(1,dsecs);
  264. else if (TMath::Abs(dsecs) > nSec / 2 && TMath::Abs(dsecs) < nSec - 1) secs[1] = secs[0] - TMath::Sign(1,dsecs);
  265. if (secs[1] < 0) secs[1] += nSec;
  266. */
  267. // Select rec. points (road) to compute distance to the track
  268. for (Int_t io = 0; io < 2; ++io) {
  269. if (fRecPoints[secs[io]].count(thmin) > 1 || fRecPoints[secs[io]].count(thmax) > 1)
  270. { cout << " !!! The same coordinate !!! " << endl; exit(0); }
  271. pits[io].first = fRecPoints[secs[io]].lower_bound(thmin);
  272. pits[io].second = fRecPoints[secs[io]].upper_bound(thmax);
  273. }
  274. phitr[1] = MpdKalmanFilter::Instance()->Proxim(phitr[0],phitr[1],fGeoPar->GetNrows()/TMath::TwoPi());
  275. TVector2 track (thetr[1]-thetr[0], phitr[1]-phitr[0]);
  276. Double_t trLeng = track.Mod();
  277. Int_t imatch = fMatchArray->GetEntriesFast();
  278. for (Int_t io = 0; io < 2; ++io) {
  279. if (io && secs[io] == secs[0]) continue; // the same sector
  280. for (multimap<Double_t,Int_t>::iterator it = pits[io].first; it != pits[io].second; ++it) {
  281. MpdTpcHit *recp = (MpdTpcHit*) fRecPointArray->UncheckedAt(it->second);
  282. Double_t phiRecp = recp->GetLocalX();
  283. phiRecp = MpdKalmanFilter::Instance()->Proxim(phitr[0],phiRecp,fGeoPar->GetNrows()/TMath::TwoPi());
  284. TVector2 vect (it->first-thetr[0], phiRecp-phitr[0]);
  285. Double_t distT = TMath::Sin(vect.DeltaPhi(track)) * vect.Mod();
  286. Double_t distL = TMath::Cos(vect.DeltaPhi(track)) * vect.Mod();
  287. MpdEmcMatch *match = new ((*fMatchArray)[imatch++]) MpdEmcMatch(itrack, it->second, distT, distL,
  288. recp->GetQ(), trLeng);
  289. // Electrons
  290. Int_t ireg = 0;
  291. if (distL >= -1.0 && distL <= trLeng + 1.0) ireg = 1;
  292. else if (distL > trLeng) { distL -= (trLeng + 1.0); ireg = 2; }
  293. Double_t chi2 = distT * distT / sigTel[ireg] / sigTel[ireg];
  294. if (distL > -3.0 && distL < -1.0) ireg = 3;
  295. chi2 += (distL - depth[ireg]) * (distL - depth[ireg]) / sigLel[ireg] / sigLel[ireg];
  296. match->SetChi2(chi2);
  297. // Pions
  298. ireg = 0;
  299. if (distL >= -1.0 && distL <= trLeng + 1.0) ireg = 1;
  300. else if (distL > trLeng) { distL -= (trLeng + 1.0); ireg = 2; }
  301. chi2 = distT * distT / sigTpi[ireg] / sigTpi[ireg];
  302. if (ireg == 1) chi2 += (distL - depthPi) * (distL - depthPi) / sigLpi[ireg] / sigLpi[ireg];
  303. else chi2 += (distL * distL / sigLpi[ireg] / sigLpi[ireg]);
  304. match->SetChi2pi(chi2);
  305. }
  306. }
  307. }
  308. //__________________________________________________________________________
  309. void MpdEmcMatching::GetTowerCoords(TVector3 &pos, Int_t io, Double_t &phiT, Double_t &theT)
  310. {
  311. // Transform point coordinates to tower indices
  312. static Int_t first = 1, offset = 0, nphi = 0;
  313. static TSpline3 *phiS, *theSmin, *theSmax;
  314. if (first) {
  315. // Get phi and theta angles of the tower centers at their inner face
  316. first = 0;
  317. const vector<Double_t> &phis = fGeoPar->GetPhiRow();
  318. const vector<Double_t> &thes = fGeoPar->GetThetaBox();
  319. const vector<Double_t> &rhos = fGeoPar->GetRhoCenterBox();
  320. const vector<Double_t> &zs = fGeoPar->GetZCenterBox();
  321. nphi = phis.size();
  322. // Offset due to the fact that the 1'st sector starts at phi = -Phi_sec/2;
  323. offset = nphi / (fGeoPar->GetNsec()/2 - 1) / 2;
  324. Double_t *phia = new Double_t [nphi];
  325. Double_t *ind = new Double_t [nphi];
  326. for (Int_t j = 0; j < nphi; ++j) {
  327. phia[j] = phis[j];
  328. ind[j] = j;
  329. }
  330. phiS = new TSpline3("phis",phia,ind,nphi); // ind vs phi
  331. delete [] phia;
  332. delete [] ind;
  333. Int_t nthe = thes.size();
  334. Double_t *theMin = new Double_t [nthe];
  335. Double_t *theMax = new Double_t [nthe];
  336. ind = new Double_t [nthe];
  337. Double_t height = fGeoPar->GetLengthBox(); // tower half-height
  338. Int_t j1 = -1;
  339. for (Int_t j = nthe-1; j >= 0; --j) {
  340. Double_t rho = rhos[j];
  341. Double_t z = zs[j];
  342. Double_t theta1 = thes[j];
  343. if (j < nthe-1 && thes[j] <= thes[j+1]+0.1) theta1 = 180 - theta1;
  344. Double_t costhe = TMath::Cos(theta1*TMath::DegToRad());
  345. Double_t sinthe = TMath::Sin(theta1*TMath::DegToRad());
  346. rho -= height * sinthe;
  347. z -= height * costhe;
  348. theMin[++j1] = TMath::ATan2(rho,z) * TMath::RadToDeg(); // at inner tower face
  349. ind[j1] = j;
  350. rho += 2 * height * sinthe;
  351. z += 2 * height * costhe;
  352. theMax[j1] = TMath::ATan2(rho,z) * TMath::RadToDeg(); // at outer face
  353. }
  354. theSmin = new TSpline3("theMin",theMin,ind,nthe); // ind vs theta
  355. theSmax = new TSpline3("theMax",theMax,ind,nthe); // ind vs theta
  356. delete [] theMin;
  357. delete [] theMax;
  358. delete [] ind;
  359. }
  360. // Get Theta and Phi tower indices at inner or outer face
  361. Double_t phi = pos.Phi() * TMath::RadToDeg();
  362. Double_t theta = pos.Theta() * TMath::RadToDeg();
  363. if (phi < 0) phi += 360;
  364. phiT = phiS->Eval(phi) + offset;
  365. if (phiT > nphi) phiT -= nphi;
  366. if (io == 0) theT = theSmin->Eval(theta);
  367. else theT = theSmax->Eval(theta);
  368. }
  369. //__________________________________________________________________________
  370. ClassImp(MpdEmcMatching)