MpdTpcHitProducer.cxx 33 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898
  1. /// \class MpdTpcHitProducer
  2. ///
  3. /// Hit producer in MPD TPC
  4. /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
  5. //---------------------------------------------------------------------------
  6. #include "MpdTpcHitProducer.h"
  7. #include "MpdTpcSectorGeo.h"
  8. #include "MpdKalmanFilter.h"
  9. #include "TpcGeoPar.h"
  10. #include "TpcPoint.h"
  11. #include "FairEventHeader.h"
  12. #include "FairGeoNode.h"
  13. #include "FairRootManager.h"
  14. #include "FairRunAna.h"
  15. #include "FairRuntimeDb.h"
  16. #include <TF1.h>
  17. #include <TGeoManager.h>
  18. #include <TGeoTube.h>
  19. #include <TSpline.h>
  20. //#include "Math/Interpolator.h"
  21. #include <iostream>
  22. using namespace std;
  23. //---------------------------------------------------------------------------
  24. MpdTpcHitProducer::MpdTpcHitProducer()
  25. : FairTask("TPC Hit Producer"),
  26. fModular(0),
  27. fPersistance(kFALSE)
  28. {
  29. fPointArray = nullptr, fHitArray = nullptr;
  30. }
  31. //---------------------------------------------------------------------------
  32. MpdTpcHitProducer::~MpdTpcHitProducer() { }
  33. //---------------------------------------------------------------------------
  34. InitStatus MpdTpcHitProducer::Init()
  35. {
  36. FairRootManager* ioman = FairRootManager::Instance();
  37. if(!ioman) {
  38. cout << "\n-E- [MpdTpcHitProducer::Init]: RootManager not instantiated!" << endl;
  39. return kFATAL;
  40. }
  41. fPointArray = (TClonesArray*) ioman->GetObject("TpcPoint");
  42. if(!fPointArray) {
  43. cout << "\n-W- [MpdTpcHitProducer::Init]: No TpcPoint array!" << endl;
  44. return kERROR;
  45. }
  46. // Create and register output array
  47. fHitArray = new TClonesArray("MpdTpcHit");
  48. if (fPersistance) ioman->Register("TpcHit", "TPC", fHitArray, kTRUE);
  49. else ioman->Register("TpcHit", "TPC", fHitArray, kFALSE);
  50. cout << "-I- MpdTpcHitProducer: Intialisation successfull." << endl;
  51. return kSUCCESS;
  52. }
  53. // ----- Private method SetParContainers -------------------------------
  54. void MpdTpcHitProducer::SetParContainers() {
  55. //return;
  56. // Get run and runtime database
  57. FairRunAna* run = FairRunAna::Instance();
  58. if ( ! run ) Fatal("SetParContainers", "No analysis run");
  59. FairRuntimeDb* db = run->GetRuntimeDb();
  60. if ( ! db ) Fatal("SetParContainers", "No runtime database");
  61. // Get TPC geometry parameter container
  62. db->getContainer("TpcGeoPar");
  63. //db->getContainer("MpdTofGeoPar");
  64. }
  65. //---------------------------------------------------------------------------
  66. void MpdTpcHitProducer::Exec(Option_t* opt)
  67. {
  68. static Int_t eventNo = 0;
  69. if (!fHitArray) Fatal("Exec", "No MpdTpcHitArray");
  70. cout << "\n-I- MpdTpcHitProducer: Event No. " << FairRun::Instance()->GetEventHeader()->GetMCEntryNumber() << " " << ++eventNo << endl;
  71. // Reset output array
  72. fHitArray->Delete();
  73. /// Merge hits if they belong to the same track and
  74. /// not too far from each other in Z: digitization emulation (interim solution).
  75. //fhLays->Reset();
  76. static Int_t first = 1;//, version3 = 0;
  77. static Double_t rMin = 99999.0, rMax = 0.0, dR;
  78. Int_t lay, layMax = 0, nPoints = fPointArray->GetEntriesFast();
  79. if (first) {
  80. first = 0;
  81. //TpcPoint *point = (TpcPoint*) fPointArray->First();
  82. //TVector3 posOut;
  83. //AZ point->PositionOut(posOut);
  84. //if (posOut != TVector3(0,0,0)) version3 = 1; // sectored and layered sensitive volume
  85. FairRuntimeDb* rtdb = FairRun::Instance()->GetRuntimeDb();
  86. //rtdb->printParamContexts();
  87. //cout << rtdb->findContainer("TpcGeoPar") << " " << rtdb->findContainer("TpcContFact:") << endl;
  88. TpcGeoPar *geoPar = (TpcGeoPar*) rtdb->getContainer("TpcGeoPar");
  89. TString volName = "tpc01sv";
  90. TObjArray* sensNodes = geoPar->GetGeoSensitiveNodes();
  91. //cout << sensNodes->GetEntriesFast() << " " << geoPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  92. //FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->FindObject(volName));
  93. Int_t nSens = sensNodes->GetEntriesFast(), nLays = 0;
  94. FairGeoNode* sensVol0 = nullptr;
  95. for (Int_t i = 0; i < nSens; ++i) {
  96. FairGeoNode* sensVol = (FairGeoNode*) (sensNodes->UncheckedAt(i));
  97. TString name = sensVol->GetName();
  98. if (!name.Contains("tpc")) continue;
  99. TArrayD* params = sensVol->getParameters();
  100. rMin = TMath::Min (params->At(0), rMin);
  101. rMax = TMath::Max (params->At(1), rMax);
  102. sensVol0 = sensVol;
  103. nLays = TMath::Max (nLays, TString(name(name.First('r')+1,15)).Atoi());
  104. }
  105. ++nLays;
  106. TObjArray* passNodes = geoPar->GetGeoPassiveNodes();
  107. if (nSens) {
  108. // Old geometry scheme
  109. if (sensVol0->getShape() == "PGON") {
  110. fModular = 1; // force using modular geometry
  111. FairGeoNode *inWall = (FairGeoNode*) passNodes->FindObject("tpc01InWall");
  112. TArrayD* params = inWall->getParameters();
  113. fZtpc = params->At(2);
  114. } else {
  115. FairGeoNode *tpc = (FairGeoNode*) passNodes->FindObject("tpcChamber1");
  116. rMax = tpc->getParameters()->At(1);
  117. dR = (rMax - rMin) / nLays;
  118. fZtpc = sensVol0->getParameters()->At(2);
  119. MpdTpcSectorGeo *geo = MpdTpcSectorGeo::Instance();
  120. geo->SetNofRows(nLays);
  121. geo->SetPadHeight(dR);
  122. geo->SetMinY(rMin);
  123. }
  124. cout << " *** TPC sensitive volume: " << sensVol0->GetName() << " "
  125. << rMin << " " << rMax << " " << fZtpc << " " << nLays << " " << dR << endl;
  126. } else {
  127. // New geometry scheme (ROOT geo)
  128. fModular = 1; // force using modular geometry
  129. TGeoVolume *inW = gGeoManager->GetVolume("tpc01InWall");
  130. TGeoTube *tube = (TGeoTube*) inW->GetShape();
  131. fZtpc = tube->GetDZ();
  132. }
  133. /*
  134. MpdTofGeoPar *tofPar = (MpdTofGeoPar*) rtdb->getContainer("MpdTofGeoPar");
  135. cout << tofPar << endl;
  136. //TString volName = "tpc01#1";
  137. //TString volName = "tpc01sv";
  138. TObjArray* sensNodesTof = tofPar->GetGeoSensitiveNodes();
  139. cout << sensNodesTof->GetEntriesFast() << " " << tofPar->GetGeoPassiveNodes()->GetEntriesFast() << endl;
  140. */
  141. //exit(0);
  142. }
  143. // !!!!!!!!! Geometry dependent values
  144. //Double_t rMin = 34.19, rMax = 99.81, dR = (rMax-rMin)/50; // 1.3124; // new version (with dead material)
  145. // !!!!!!!!!
  146. //AZ if (version3) { ExecNew(); return; } // sectored and layered sensitive volume
  147. if (1) { ExecNew(); return; } // more streamlined procedure
  148. if (fModular) { ExecModular(); return; } // emulate geometry of readout chambers
  149. TVector3 p3, p30, p3err(0.05, 0., 0.1); // R-Phi error 500 um, Z error 1 mm
  150. multiset<Int_t> layset;
  151. for (Int_t j = 0; j < nPoints; ++j ) {
  152. TpcPoint* point = (TpcPoint*) fPointArray->UncheckedAt(j);
  153. point->Position(p3);
  154. Double_t rad = p3.Pt();
  155. lay = (Int_t) ((rad - rMin) / dR);
  156. lay = TMath::Max (lay, 0);
  157. layMax = TMath::Max (lay, layMax);
  158. //MpdTpcHit *hit = new ((*fHitArray)[j]) MpdTpcHit();
  159. MpdTpcHit *hit = AddRawHit(j, point->GetDetectorID(), p3, p3err, j, point->GetTrackID());
  160. layset.insert(lay);
  161. hit->SetLayer(lay);
  162. hit->SetR(rad);
  163. hit->SetRphi(rad * p3.Phi());
  164. hit->SetLocalZ(p3.Z());
  165. hit->SetEnergyLoss(point->GetEnergyLoss());
  166. hit->SetStep(point->GetStep());
  167. //fhLays->Fill(lay+0.1);
  168. }
  169. cout << " Max layer = " << layMax << " " << fPointArray->GetEntriesFast() << endl;
  170. fHitArray->Sort(); // according to layer No. or radius (local y)
  171. Int_t ipos = 0;
  172. Int_t nHits = 0;
  173. for (Int_t i = layMax; i >= 0; --i) {
  174. ipos += nHits;
  175. nHits = layset.count(i);
  176. //nHits = TMath::Nint(fhLays->GetCellContent(i+1,0));
  177. if (nHits < 2) continue;
  178. Int_t *iHits = new Int_t [nHits];
  179. Int_t *index = new Int_t [nHits];
  180. for (Int_t j = 0; j < nHits; ++j) {
  181. Int_t k = j + ipos;
  182. iHits[j] = ((MpdTpcHit*)fHitArray->UncheckedAt(k))->GetTrackID();
  183. }
  184. TMath::Sort(nHits,iHits,index);
  185. for (Int_t j = 0; j < nHits; ++j) {
  186. // Loop over first hits
  187. MpdTpcHit *hit0 = (MpdTpcHit*) fHitArray->UncheckedAt(index[j]+ipos);
  188. if (hit0 == 0x0) continue;
  189. Int_t nMerge = 1;
  190. TpcPoint* point0 = (TpcPoint*) fPointArray->UncheckedAt(hit0->GetRefIndex());
  191. Double_t rPhi = hit0->GetRphi();
  192. Double_t z = hit0->GetLocalZ();
  193. Double_t rad = hit0->GetR();
  194. hit0->Position(p30);
  195. Double_t leng = point0->GetLength();
  196. Double_t edep = point0->GetEnergyLoss();
  197. Double_t step = point0->GetStep();
  198. for (Int_t j1 = j + 1; j1 < nHits; ++j1) {
  199. // Second hit
  200. Int_t k = index[j1] + ipos;
  201. MpdTpcHit *hit = (MpdTpcHit*) fHitArray->UncheckedAt(k);
  202. if (hit == 0x0) continue;
  203. TpcPoint* point = (TpcPoint*) fPointArray->UncheckedAt(hit->GetRefIndex());
  204. if (point->GetTrackID() != point0->GetTrackID()) break;
  205. //if (TMath::Abs(point->GetZ()-point0->GetZ()) > 5) continue;
  206. //if (TMath::Abs(hit0->GetRphi()-Proxim(hit0,hit)) > 5) continue;
  207. if (TMath::Abs(point->GetZ()-point0->GetZ()) > 2.5) continue;
  208. if (TMath::Abs(hit0->GetRphi()-Proxim(hit0,hit)) > 2.5) continue;
  209. // Merge hits
  210. hit->Position(p3);
  211. p30 += p3;
  212. rPhi += Proxim(hit0,hit);
  213. z += hit->GetLocalZ();
  214. rad += hit->GetR();
  215. leng += point->GetLength();
  216. edep += point->GetEnergyLoss();
  217. step += point->GetStep();
  218. hit0->AddLinks(hit->GetLinks()); // copy links
  219. fHitArray->RemoveAt(k);
  220. ++nMerge;
  221. }
  222. if (nMerge == 1) continue;
  223. //cout << i << " nMerge: " << nMerge << " " << hit0->GetRphi() << " " << hit0->GetZ() << " " << hit0->GetR() << endl;
  224. p30.SetMag(p30.Mag()/nMerge);
  225. hit0->SetPosition(p30);
  226. hit0->SetRphi(rPhi/nMerge);
  227. hit0->SetLocalZ(z/nMerge);
  228. hit0->SetR(rad/nMerge);
  229. //((TpcPoint*)fTpcPoints->UncheckedAt(hit0->GetRefIndex()))->SetLength(leng/nMerge); //
  230. hit0->SetLength(leng/nMerge);
  231. hit0->SetEnergyLoss(edep);
  232. hit0->SetStep(step);
  233. //cout << hit0->GetRphi() << " " << hit0->GetZ() << " " << hit0->GetR() << " " << hit0->GetTrackID() << endl;
  234. //cout << nMerge << " " << hit0->GetTrackID() << " " << hit0->GetLayer() << endl;
  235. }
  236. delete [] iHits;
  237. delete [] index;
  238. } // for (Int_t i = layMax; i >= 0;
  239. fHitArray->Compress();
  240. fHitArray->UnSort();
  241. fHitArray->Sort();
  242. /*
  243. fhLays->Reset();
  244. nHits = fHits->GetEntriesFast();
  245. cout << " Number of merged hits: " << nHits << endl;
  246. Int_t nKh = 0;
  247. for (int j = 0; j < nHits; ++j ) {
  248. TpcLheHit* hit = (TpcLheHit*) fHits->UncheckedAt(j);
  249. //if (hit->GetTrackID() != 1783) continue; // keep only one track
  250. lay = hit->GetLayer();
  251. fhLays->Fill(lay+0.1);
  252. // Add errors
  253. Double_t dRPhi = 0, dZ = 0;
  254. gRandom->Rannor(dRPhi,dZ);
  255. hit->SetZ(hit->GetZ()+dZ*hit->GetZerr()); //add error
  256. hit->SetRphi(hit->GetRphi()+dRPhi*hit->GetRphiErr()); //add error
  257. Double_t phi = hit->GetRphi() / hit->GetR();
  258. hit->SetX(hit->GetR()*TMath::Cos(phi));
  259. hit->SetY(hit->GetR()*TMath::Sin(phi));
  260. // Create Kalman hits
  261. //(Int_t detID, Int_t nDim, HitType hitType, Double_t *meas, Double_t *err, Double_t *cosSin, Double_t signal, Double_t dist, Int_t index)
  262. Double_t meas[2] = {hit->GetRphi(), hit->GetZ()};
  263. Double_t err[2] = {hit->GetRphiErr(), hit->GetZerr()};
  264. Double_t cossin[2] = {1., 0.};
  265. MpdKalmanHit *hitK = new ((*fKHits)[nKh++])
  266. MpdKalmanHit(lay*1000000+nKh-1, 2, MpdKalmanHit::kFixedR, meas, err, cossin, hit->GetEdep()/hit->GetStep(), hit->GetR(), hit->GetRefIndex());
  267. hitK->SetLength(((TpcPoint*)fTpcPoints->UncheckedAt(hit->GetRefIndex()))->GetLength());
  268. //hitK->SetDedx (hit->GetEdep()/hit->GetStep());
  269. //MpdKalmanFilter::Instance()->GetGeo()->SetGlobalPos(hitK,TVector3(hit->GetX(),hit->GetY(),hit->GetZ()));
  270. }
  271. fLayPointers = new Int_t [layMax+1];
  272. ipos = 0;
  273. for (Int_t i = layMax; i >= 0; --i) {
  274. //cout << i << " " << fhLays->GetCellContent(i+1,0) << endl;
  275. //if (ipos) cout << ((TpcLheHit*)fHits->UncheckedAt(ipos))->GetLayer() << " "
  276. // << ((TpcLheHit*)fHits->UncheckedAt(ipos-1))->GetLayer() << endl;
  277. fLayPointers[i] = ipos;
  278. ipos += (Int_t) fhLays->GetCellContent(i+1,0);
  279. }
  280. */
  281. // Event summary
  282. cout << "-I- MpdTpcHitProducer: " << nPoints << " TpcPoints, " << fHitArray->GetEntriesFast() << " Hits created." << endl;
  283. }
  284. //---------------------------------------------------------------------------
  285. void MpdTpcHitProducer::ExecModular()
  286. {
  287. // Emulate geometry of readout chambers
  288. MpdTpcSectorGeo *secGeo = MpdTpcSectorGeo::Instance();
  289. //ROOT::Math::Interpolator inter(3, ROOT::Math::Interpolation::kPOLYNOMIAL);
  290. Int_t lay, layMax = 0, nPoints = fPointArray->GetEntriesFast(), nHits = 0;
  291. TVector3 p3, p30, p3local, p3local0, pmom3, pmom3loc, p3extr, p3err(0.05, 0., 0.1); // X error 500 um, Z error 1 mm
  292. multiset<Int_t> layset;
  293. multimap<Int_t,Int_t> midIndx;
  294. cout << " MC poins: " << fPointArray->GetEntriesFast() << endl;
  295. for (Int_t j = 0; j < nPoints; ++j ) {
  296. //for (Int_t j = 0; j < 1000; ++j ) {
  297. TpcPoint* point = (TpcPoint*) fPointArray->UncheckedAt(j);
  298. //if (point->GetTrackID() != 1) continue; ///
  299. if (point->GetTrackID() < 0) continue; /// strange case - protection
  300. point->Position(p3);
  301. Int_t padID = secGeo->Global2Local(p3, p3local);
  302. //if (padID < 0) {cout << j << " " << p3.X() << " " << p3.Y() << " " << p3local.X() << " " << p3local.Y() << endl; continue; }// outside sector boundaries
  303. if (padID < 0) continue; // outside sector boundaries
  304. point->Momentum(pmom3);
  305. pmom3loc.SetXYZ(0.,0.,0.);
  306. //if (point->GetTrackID() <= 2 && point->GetTrackID() > 0) printf("%2d %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f %6d %4d %4d\n", point->GetTrackID(), p3.X(), p3.Y(), p3.Z(), p3local.X(), p3local.Y(), p3local.Z(), padID, TpcPadID::numberToPadID(padID).row(), TpcPadID::numberToPadID(padID).sector());
  307. lay = secGeo->PadRow(padID);
  308. layMax = TMath::Max (lay, layMax);
  309. //MpdTpcHit *hit = AddRawHit(j, point->GetDetectorID(), p3, p3err, j, point->GetTrackID());
  310. MpdTpcHit *hit = AddRawHit(nHits++, padID, p3, p3err, j, point->GetTrackID());
  311. hit->SetLayer(lay);
  312. hit->SetLocalPosition(p3local); // point position
  313. //cout << p3local.X() << " " << p3local.Y() << " " << p3local.Z() << " " << hit->GetLocalY() << endl;
  314. hit->SetEnergyLoss(point->GetEnergyLoss());
  315. hit->SetStep(point->GetStep());
  316. hit->SetModular(1); // modular geometry flag
  317. hit->SetLength(point->GetLength());
  318. midIndx.insert(pair<Int_t,Int_t>(point->GetTrackID(),nHits-1));
  319. } // for (Int_t j = 0; j < nPoints;
  320. // Perform track interpolation to put hits on padrow median planes
  321. //
  322. multimap<Int_t,Int_t>::iterator mit;
  323. Int_t id0 = -1, i0 = 0, nh = 0, *hindx = nullptr, *ord = nullptr, lastIndx = 0;
  324. Double_t *times = nullptr, *ys = nullptr, *zs = nullptr, *xx = nullptr, *yy = nullptr, *zz = nullptr;
  325. for (mit = midIndx.begin(); mit != midIndx.end(); ++mit) {
  326. // Loop over track hits
  327. if (mit->first != id0) {
  328. id0 = mit->first;
  329. nh = midIndx.count(id0);
  330. times = new Double_t [nh];
  331. ys = new Double_t [nh];
  332. zs = new Double_t [nh];
  333. xx = new Double_t [nh];
  334. yy = new Double_t [nh];
  335. zz = new Double_t [nh];
  336. hindx = new Int_t [nh];
  337. ord = new Int_t [nh];
  338. i0 = 0;
  339. }
  340. MpdTpcHit *h = (MpdTpcHit*) fHitArray->UncheckedAt(mit->second);
  341. TpcPoint* point = (TpcPoint*) fPointArray->UncheckedAt(h->GetRefIndex());
  342. times[i0] = point->GetTime();
  343. hindx[i0++] = mit->second;
  344. if (i0 != nh) continue;
  345. TMath::Sort (nh, times, ord, kFALSE);
  346. // Take 3 consecutive (in time) points and use parabolic interpolation
  347. // to put the hit on a padrow median plane
  348. Int_t sec0 = -99, np = 0, ibeg = 0;
  349. Double_t ypad0 = -999.0, dir = 0.0, xhit = 0.0, zhit = 0.0;
  350. Bool_t ok = kTRUE;
  351. for (Int_t i3 = 0; i3 < nh; ++i3) {
  352. h = (MpdTpcHit*) fHitArray->UncheckedAt(hindx[ord[i3]]);
  353. Int_t sec = secGeo->Sector(h->GetDetectorID());
  354. if (sec != sec0 || i3 > lastIndx) {
  355. //if (sec != sec0) {
  356. // Store coordinates of all hits in this sector
  357. //cout << " ***** New sector: " << sec << endl;
  358. sec0 = sec;
  359. np = ibeg = 0;
  360. ypad0 = -999.0;
  361. for (Int_t jj = i3; jj < nh; ++jj) {
  362. MpdTpcHit *h1 = (MpdTpcHit*) fHitArray->UncheckedAt(hindx[ord[jj]]);
  363. sec = secGeo->Sector(h1->GetDetectorID());
  364. if (sec != sec0) break; // different sector
  365. yy[np] = h1->GetLocalY();
  366. xx[np] = h1->GetLocalX();
  367. if (np > 0 && TMath::Abs(xx[np]-xx[np-1]) > 5.0) break; // broken track
  368. zz[np++] = h1->GetLocalZ();
  369. lastIndx = jj;
  370. }
  371. /*for (Int_t jj = 0; jj < 6; ++jj) cout << xx[jj] << " " << yy[jj] << " " << zz[jj] << " ";
  372. cout << endl;*/
  373. }
  374. Double_t ypad = secGeo->LocalPadPosition(h->GetDetectorID()).Y(); // padrow position
  375. if (TMath::Abs(ypad-ypad0) < 0.1) {
  376. // The same padrow
  377. if (!ok) { fHitArray->RemoveAt(hindx[ord[i3]]); continue; }
  378. times[i3] = times[i3-1];
  379. ys[i3] = ys[i3-1];
  380. zs[i3] = zs[i3-1];
  381. } else {
  382. ypad0 = ypad;
  383. //cout << "---------------------------\n";
  384. ok = Interpolate(np, ibeg, yy, xx, zz, ypad0, dir, xhit, zhit);
  385. if (!ok) { fHitArray->RemoveAt(hindx[ord[i3]]); continue; }
  386. times[i3] = xhit; // reuse array times
  387. ys[i3] = ypad0;
  388. zs[i3] = zhit;
  389. }
  390. /*cout << " np, ibeg, y, x, z: " << np << " " << ibeg << " " << yy[ibeg] << " " << xx[ibeg] << " " << zz[ibeg] << endl;
  391. cout << ys[i3] << " " << times[i3] << " " << zs[i3] << endl;*/
  392. if (ibeg < np - 2) ++ibeg;
  393. Double_t dy = ypad0 - h->GetLocalY();
  394. //h->SetLocalY(ypad0);
  395. Double_t dx = times[i3] - h->GetLocalX();
  396. Double_t dz = zs[i3] - h->GetLocalZ();
  397. Double_t dl = dx * dx + dy * dy + dz * dz;
  398. dl = TMath::Sign (TMath::Sqrt(dl), dy);
  399. if (dir > 0) h->SetLength(h->GetLength()+dl); // going outward
  400. else h->SetLength(h->GetLength()-dl); // inward
  401. } // for (Int_t i3 = 0; i3 < nh;
  402. // Change hit coordinates
  403. for (Int_t i3 = 0; i3 < nh; ++i3) {
  404. h = (MpdTpcHit*) fHitArray->UncheckedAt(hindx[ord[i3]]);
  405. if (h == nullptr) continue;
  406. Int_t padID = h->GetDetectorID();
  407. h->SetLocalX(times[i3]);
  408. h->SetLocalY(ys[i3]);
  409. h->SetLocalZ(zs[i3]);
  410. h->LocalPosition(p30);
  411. secGeo->Local2Global(secGeo->Sector(padID), p30, p3);
  412. h->SetPosition(p3);
  413. }
  414. delete [] times;
  415. delete [] ys;
  416. delete [] zs;
  417. delete [] xx;
  418. delete [] yy;
  419. delete [] zz;
  420. // Merge hits
  421. for (Int_t i3 = 0; i3 < nh; ++i3) {
  422. // Loop over first hits
  423. MpdTpcHit *hit0 = (MpdTpcHit*) fHitArray->UncheckedAt(hindx[ord[i3]]);
  424. if (hit0 == 0x0) continue;
  425. Int_t padID0 = hit0->GetDetectorID();
  426. Int_t nMerge = 1;
  427. hit0->LocalPosition(p3local0);
  428. hit0->LocalPosition(p3extr);
  429. hit0->Position(p30);
  430. Double_t leng = hit0->GetLength();
  431. Double_t edep = hit0->GetEnergyLoss();
  432. Double_t step = hit0->GetStep();
  433. for (Int_t i31 = i3 + 1; i31 < nh; ++i31) {
  434. //continue; ///
  435. // Second hit
  436. MpdTpcHit *hit = (MpdTpcHit*) fHitArray->UncheckedAt(hindx[ord[i31]]);
  437. if (hit == 0x0) continue;
  438. Int_t padID = hit->GetDetectorID();
  439. if (secGeo->PadRow(padID) != secGeo->PadRow(padID0)) break; // in different padrows
  440. if (secGeo->Sector(padID) != secGeo->Sector(padID0)) break; // in different sectors
  441. hit->LocalPosition(p3local);
  442. hit->Position(p3);
  443. if (TMath::Abs(p3extr.Z()-p3local.Z()) > 2.5) break;
  444. if (TMath::Abs(p3extr.X()-p3local.X()) > 2.5) break;
  445. // Merge hits
  446. //cout << " Merge: " << p3extr.X() << " " << p3extr.Y() << " " << p3extr.Z() << " " << p3local.X() << " " << p3local.Y() << " " << p3local.Z() << endl;
  447. p3local0 += p3local;
  448. p30 += p3;
  449. leng += hit->GetLength();
  450. edep += hit->GetEnergyLoss();
  451. step += hit->GetStep();
  452. //hit0->AddLinks(hit->GetLinks()); // copy links
  453. fHitArray->RemoveAt(hindx[ord[i31]]);
  454. ++nMerge;
  455. }
  456. if (nMerge == 1) continue;
  457. p3local0.SetMag(p3local0.Mag()/nMerge);
  458. hit0->SetLocalPosition(p3local0);
  459. p30.SetMag(p30.Mag()/nMerge);
  460. hit0->SetPosition(p30);
  461. hit0->SetLength(leng/nMerge);
  462. hit0->SetEnergyLoss(edep);
  463. hit0->SetStep(step);
  464. } // for (Int_t i3 = 0; i3 < nh;
  465. delete [] hindx;
  466. delete [] ord;
  467. } // for (mit = midIndx.begin(); mit != midIndx.end();
  468. fHitArray->Compress();
  469. fHitArray->UnSort();
  470. fHitArray->Sort();
  471. cout << " Merged hits: " << fHitArray->GetEntriesFast() << endl;
  472. // Debug
  473. /*
  474. nh = fHitArray->GetEntriesFast();
  475. multimap<Int_t,Int_t> mm;
  476. for (Int_t i = 0; i < nh; ++i) {
  477. MpdTpcHit *h = (MpdTpcHit*) fHitArray->UncheckedAt(i);
  478. //cout << i << " " << h->GetTrackID() << " " << h->GetLayer() << endl;
  479. //mm.insert(pair<Int_t,Int_t>(h->GetTrackID(),h->GetLayer()));
  480. mm.insert(pair<Int_t,Int_t>(h->GetTrackID(),i));
  481. }
  482. multimap<Int_t,Int_t>::iterator it;
  483. for (it = mm.begin(); it != mm.end(); ++it) {
  484. if (it->first != 6 && it->first != 7 && it->first != 0) continue;
  485. MpdTpcHit *h = (MpdTpcHit*) fHitArray->UncheckedAt(it->second);
  486. cout << it->first << " " << h->GetLayer() << " " << h->GetLocalY() << " " << h->GetLocalX() << " "
  487. << h->GetLocalZ() << " " << h->GetY() << " " << h->GetX() << " " << h->GetZ() << endl;
  488. }
  489. */
  490. }
  491. //---------------------------------------------------------------------------
  492. Double_t MpdTpcHitProducer::Proxim(const MpdTpcHit *hit0, const MpdTpcHit *hit)
  493. {
  494. /// Adjust hit coord. R-Phi to be "around" hit0 R-Phi - to avoid
  495. /// discontinuity around +- Pi
  496. Double_t phi0 = hit0->GetRphi() / hit0->GetR();
  497. Double_t phi = hit->GetRphi() / hit->GetR();
  498. return hit->GetR() * MpdKalmanFilter::Instance()->Proxim(phi0, phi);
  499. }
  500. //---------------------------------------------------------------------------
  501. MpdTpcHit* MpdTpcHitProducer::AddRawHit(Int_t indx, Int_t detUID, const TVector3 &posHit,
  502. const TVector3 &posHitErr, Int_t pointIndx,
  503. Int_t trackIndx)
  504. {
  505. MpdTpcHit *pHit = new ((*fHitArray)[indx]) MpdTpcHit(detUID, posHit, posHitErr, pointIndx);
  506. pHit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  507. pHit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, trackIndx));
  508. pHit->AddID(trackIndx);
  509. return pHit;
  510. }
  511. //---------------------------------------------------------------------------
  512. Bool_t MpdTpcHitProducer::Interpolate(Int_t np, Int_t& ibeg, Double_t *yp, Double_t *xp, Double_t *zp,
  513. Double_t y0, Double_t& dir, Double_t& xhit, Double_t& zhit)
  514. {
  515. // Evaluate x and z at y0 by interpolation (either parabolic (if np=3), or linear (if np=2))
  516. // (here y0 is the padrow y-coordinate)
  517. if (np == 1) {
  518. xhit = xp[0];
  519. zhit = zp[0];
  520. return kTRUE;
  521. }
  522. if (np == 2) {
  523. dir = TMath::Sign (1.0,yp[1]-yp[0]);
  524. // Linear interpolation
  525. /*
  526. Double_t expr = (y0 - yp[0]) / (yp[1] - yp[0]);
  527. xhit = xp[0] + (xp[1] - xp[0]) * expr;
  528. zhit = zp[0] + (zp[1] - zp[0]) * expr;
  529. */
  530. xhit = (xp[0] + xp[1]) / 2.0;
  531. zhit = (zp[0] + zp[1]) / 2.0;
  532. return kTRUE;
  533. } else if (np < 2) {
  534. Fatal("Interpolate", "No points !!!");
  535. }
  536. Double_t *yp3{nullptr}, *xp3{nullptr}, *zp3{nullptr};
  537. Bool_t ok = kTRUE;
  538. // Go to the right
  539. for (Int_t j = ibeg; j < np; ++j) {
  540. ok = kTRUE;
  541. yp3 = &yp[j];
  542. xp3 = &xp[j];
  543. zp3 = &zp[j];
  544. if (j > np - 3) { --yp3; --xp3; --zp3; --ibeg; }
  545. if (j > np - 2) { --yp3; --xp3; --zp3; --ibeg; }
  546. Double_t dy0 = yp3[0] - y0;
  547. Double_t dy2 = yp3[2] - y0;
  548. if (dy0 * dy2 > 0.0) {
  549. // Both points on the same side from y0
  550. if (TMath::Abs(dy2) < TMath::Abs(dy0) && j < np - 3) { ++ibeg; continue; }
  551. //ok = kFALSE;
  552. if (j < np - 3) ok = kFALSE; // go left or turning point
  553. break;
  554. } else if ((yp3[1]-yp3[0]) * (yp3[2]-yp3[1]) < 0) ok = kFALSE; // not monotone
  555. break;
  556. }
  557. // Go to the left
  558. if (!ok) {
  559. for (Int_t j = ibeg; j > -1; --j) {
  560. ok = kTRUE;
  561. yp3 = &yp[j-2];
  562. xp3 = &xp[j-2];
  563. zp3 = &zp[j-2];
  564. //if (j < 2) { ++xp3; ++yp3; ++ibeg; }
  565. //if (j < 1) { ++xp3; ++yp3; ++ibeg; }
  566. if (j < 2) { ++yp3; ++xp3; ++zp3; }
  567. if (j < 1) { ++yp3; ++xp3; ++zp3; }
  568. Double_t dy0 = yp3[0] - y0;
  569. Double_t dy2 = yp3[2] - y0;
  570. if (dy0 * dy2 > 0.0) {
  571. // Both points on the same side from y0
  572. if (TMath::Abs(dy0) < TMath::Abs(dy2) && j > 2) { --ibeg; continue; }
  573. if (j > 0) ok = kFALSE; // turning point
  574. break;
  575. } else if ((yp3[1]-yp3[0]) * (yp3[2]-yp3[1]) < 0) ok = kFALSE; // not monotone
  576. break;
  577. }
  578. }
  579. dir = TMath::Sign (1.0,yp3[1]-yp3[0]);
  580. if (!ok) {
  581. xhit = (xp3[0] + xp3[1] + xp3[2]) / 3.0;
  582. zhit = (zp3[0] + zp3[1] + zp3[2]) / 3.0;
  583. return kTRUE;
  584. }
  585. /*for (Int_t jj = 0; jj < 3; ++jj) cout << yp3[jj] << " ";
  586. cout << endl;
  587. for (Int_t jj = 0; jj < 3; ++jj) cout << xp3[jj] << " ";
  588. cout << endl;
  589. for (Int_t jj = 0; jj < 3; ++jj) cout << zp3[jj] << " ";
  590. cout << endl;*/
  591. // Parabolic interpolation
  592. Double_t dy01 = yp3[1] - yp3[0], dy02 = yp3[2] - yp3[0];
  593. Double_t y01 = yp3[1] + yp3[0], y02 = yp3[2] + yp3[0];
  594. //Double_t yp302 = yp3[0] * yp3[0];
  595. Double_t ydy = dy02 * y02 - dy02 * y01;
  596. Double_t dx01 = xp3[1] - xp3[0], dx02 = xp3[2] - xp3[0];
  597. Double_t slopex = dx01 / dy01;
  598. Double_t ax = dx02 - slopex * dy02;
  599. ax /= ydy;
  600. Double_t bx = slopex - ax * y01;
  601. Double_t cx = xp3[0] - ax * yp3[0] * yp3[0] - bx * yp3[0];
  602. xhit = ax * y0 * y0 + bx * y0 + cx;
  603. // Check if x is inside sector boundaries
  604. // (can happen when crossing sector boundaries (in some cases))
  605. MpdTpcSectorGeo *secGeo = MpdTpcSectorGeo::Instance();
  606. Double_t xEdge = (y0 + secGeo->GetMinY()) * TMath::Tan(secGeo->Dphi()/2) + 0.2; // safety margin 0.2
  607. if (TMath::Abs(xhit) > xEdge) {
  608. return kFALSE;
  609. // Find y corresponding to sector boundary
  610. /*Double_t d = bx * bx - 4 * ax * (cx - TMath::Sign(xEdge,xhit));
  611. y0 = (-bx - TMath::Sqrt(d)) / 2 / ax;
  612. xhit = ax * y0 * y0 + bx * y0 + cx;*/
  613. }
  614. Double_t dz01 = zp3[1] - zp3[0], dz02 = zp3[2] - zp3[0];
  615. Double_t slopez = dz01 / dy01;
  616. Double_t az = dz02 - slopez * dy02;
  617. az /= ydy;
  618. Double_t bz = slopez - az * y01;
  619. Double_t cz = zp3[0] - az * yp3[0] * yp3[0] - bz * yp3[0];
  620. zhit = az * y0 * y0 + bz * y0 + cz;
  621. if (TMath::Abs(zhit) > fZtpc + 1) return kFALSE;
  622. return kTRUE;
  623. }
  624. //---------------------------------------------------------------------------
  625. void MpdTpcHitProducer::ExecNew()
  626. {
  627. // More streamlined version
  628. MpdTpcSectorGeo *secGeo = MpdTpcSectorGeo::Instance();
  629. static TF1 func("func","[1]+[2]*(x-[0])+[3]*(x-[0])*(x-[0])",0,99999);
  630. Int_t nPoints = fPointArray->GetEntriesFast(), nHits = 0, lay = 0;
  631. TVector3 p3, p30, p3loc, p3loc0, pmom3, pmom3loc, p3extr, p3err(0.05, 0., 0.1); // X error 500 um, Z error 1 mm
  632. map<Int_t,map<Double_t,Int_t> > idmap;
  633. //multimap<Int_t,Int_t> midIndx;
  634. cout << " MC poins: " << nPoints << endl;
  635. // Get all points according to trackID
  636. for (Int_t j = 0; j < nPoints; ++j ) {
  637. //for (Int_t j = 0; j < 1000; ++j ) {
  638. TpcPoint* point = (TpcPoint*) fPointArray->UncheckedAt(j);
  639. Int_t id = point->GetTrackID();
  640. //if (id != 1) continue; ///
  641. if (id < 0) continue; /// strange case - protection
  642. if (idmap.count(id) == 0) {
  643. map<Double_t,Int_t> aaa;
  644. idmap[id] = aaa;
  645. }
  646. idmap[id][point->GetTime()] = j;
  647. }
  648. // Loop over trackIDs
  649. for (map<Int_t,map<Double_t,Int_t> >::iterator mit = idmap.begin(); mit != idmap.end(); ++mit) {
  650. Int_t id = mit->first;
  651. map<Double_t,Int_t>& aaa = mit->second;
  652. map<Double_t,MpdTpcHit> hitMap;
  653. // Loop over points from one track and create hit for each point
  654. for (map<Double_t,Int_t>::iterator mit1 = aaa.begin(); mit1 != aaa.end(); ++mit1) {
  655. TpcPoint *point = (TpcPoint*) fPointArray->UncheckedAt(mit1->second);
  656. point->Position(p3);
  657. Int_t padID = secGeo->Global2Local(p3, p3loc), isec = -1;
  658. if (padID < 0) continue; // outside sector boundaries
  659. lay = -1;
  660. if (padID >= 0) {
  661. lay = secGeo->PadRow(padID);
  662. isec = secGeo->Sector(padID);
  663. }
  664. MpdTpcHit hit(padID, p3, p3err, mit1->second);
  665. hit.SetLayer(lay);
  666. hit.SetLocalPosition(p3loc); // point position
  667. hit.SetEnergyLoss(point->GetEnergyLoss());
  668. hit.SetStep(point->GetStep());
  669. hit.SetModular(1); // modular geometry flag
  670. hit.SetLength(point->GetLength());
  671. // Track direction in sector frame
  672. Int_t idir = 0;
  673. if (padID >= 0) {
  674. point->Momentum(pmom3);
  675. if (pmom3.Mag() > 1.e-6) {
  676. p3extr = p3;
  677. pmom3.SetMag(0.001);
  678. p3extr += pmom3;
  679. secGeo->Global2Local(p3extr, pmom3loc, isec);
  680. pmom3loc -= p3loc;
  681. if (pmom3loc[1] < -1.e-7) idir = -1; // going inward
  682. else if (pmom3loc[1] > 1.e-7) idir = 1;
  683. }
  684. }
  685. hit.SetFlag(idir);
  686. //hit.AddLink(FairLink(MpdTpcHit::PointIndex, mit1->second));
  687. //hit.AddLink(FairLink(MpdTpcHit::MCTrackIndex, id));
  688. hitMap[mit1->first] = hit;
  689. }
  690. // Add fake hit to indicate end-of-track
  691. if (hitMap.size()) {
  692. MpdTpcHit htmp;
  693. htmp.SetDetectorID(secGeo->PadID(30,0));
  694. htmp.SetFlag(hitMap.rbegin()->second.GetFlag());
  695. hitMap[hitMap.rbegin()->first+1.0] = htmp;
  696. }
  697. // Find track segment in one sector going in one direction in sector frame
  698. Int_t isec0 = -9, idir0 = -9, layb, lay0; //laye,
  699. vector<Double_t> xyzloct[4];
  700. map<Double_t,MpdTpcHit>::iterator mitb, mite;
  701. nHits = fHitArray->GetEntriesFast();
  702. for (map<Double_t,MpdTpcHit>::iterator mit1 = hitMap.begin(); mit1 != hitMap.end(); ++mit1) {
  703. MpdTpcHit& hit = mit1->second;
  704. Int_t isec = secGeo->Sector(hit.GetDetectorID());
  705. Int_t idir = hit.GetFlag();
  706. if (isec0 < 0) {
  707. isec0 = isec;
  708. idir0 = idir;
  709. }
  710. if (isec != isec0 || idir != idir0) {
  711. // Process track segment
  712. if (xyzloct[0].size() < 2) mitb = mite = hitMap.lower_bound(xyzloct[3].front()); // Only one point
  713. else {
  714. mitb = hitMap.lower_bound(xyzloct[3].front());
  715. mite = hitMap.lower_bound(xyzloct[3].back());
  716. }
  717. layb = mitb->second.GetLayer();
  718. //laye = mite->second.GetLayer();
  719. map<Double_t,MpdTpcHit>::iterator mitt = mitb, mitend = mite, mithit;
  720. ++mitend;
  721. // Create hit on each padrow median plane
  722. Int_t spsize = TMath::Max (Int_t(xyzloct[0].size()),2);
  723. //TSpline3 tx("tx",xyzloct[3].data(),xyzloct[0].data(),spsize);
  724. //TSpline3 tz("tz",xyzloct[3].data(),xyzloct[2].data(),spsize);
  725. TSpline3* yt = nullptr, *yx = nullptr, *yz = nullptr;
  726. if (xyzloct[1].front() > xyzloct[1].back()) {
  727. // Reverse vectors (to have localY in ascending order)
  728. vector<Double_t> yrev(xyzloct[1]), trev(xyzloct[3]), xrev(xyzloct[0]), zrev(xyzloct[2]);
  729. std::reverse(yrev.begin(),yrev.end());
  730. std::reverse(trev.begin(),trev.end());
  731. std::reverse(xrev.begin(),xrev.end());
  732. std::reverse(zrev.begin(),zrev.end());
  733. yt = new TSpline3("yt",yrev.data(),trev.data(),spsize);
  734. yx = new TSpline3("yx",yrev.data(),xrev.data(),spsize);
  735. yz = new TSpline3("yz",yrev.data(),zrev.data(),spsize);
  736. }
  737. else {
  738. yt = new TSpline3("yt",xyzloct[1].data(),xyzloct[3].data(),spsize);
  739. yx = new TSpline3("yx",xyzloct[1].data(),xyzloct[0].data(),spsize);
  740. yz = new TSpline3("yz",xyzloct[1].data(),xyzloct[2].data(),spsize);
  741. }
  742. // Fit yloc vs time to find maximum and minimum of yloc achieved
  743. Double_t ylocMinMax[2] = {0, 100}, dt = TMath::Max(xyzloct[3].back()-xyzloct[3].front(),0.2);
  744. //Double_t ylocMinMax[2] = {0, 100}, dt = xyzloct[3].back()-xyzloct[3][xyzloct[3].size()-2];
  745. if (spsize > 2) {
  746. dt *= 0.2;
  747. //dt *= 5;
  748. //*
  749. TGraph gr(spsize,xyzloct[3].data(),xyzloct[1].data());
  750. func.SetParameters(0,xyzloct[1].front(),0,0);
  751. func.FixParameter(0,xyzloct[3].front());
  752. //gr.Fit("pol2","Q");
  753. gr.Fit("func","Q");
  754. ylocMinMax[0] = gr.GetFunction("func")->GetMinimum(xyzloct[3].front()-dt,xyzloct[3].back()+dt);
  755. ylocMinMax[1] = gr.GetFunction("func")->GetMaximum(xyzloct[3].front()-dt,xyzloct[3].back()+dt);
  756. //*/
  757. }
  758. Double_t time = 0.0;
  759. MpdTpcHit *hitok = nullptr;
  760. lay0 = -9;
  761. for ( ; mitt != mitend; ++mitt) {
  762. Int_t padID = mitt->second.GetDetectorID();
  763. if (mitb == mite) {
  764. // One hit
  765. for (Int_t j = 0; j < 3; ++j) p3loc[j] = xyzloct[j].front();
  766. time = xyzloct[3].front();
  767. hitok = &mitb->second;
  768. lay = layb;
  769. } else {
  770. lay = secGeo->PadRow(padID);
  771. if (lay == lay0) continue; // hit from the same padrow
  772. p3loc[1] = secGeo->LocalPadPosition(padID).Y(); // padrow position
  773. //if (p3loc[1] < ylocMinMax[0] || p3loc[1] > ylocMinMax[1]) continue; // curling track
  774. if (p3loc[1] < ylocMinMax[0]) p3loc[1] = ylocMinMax[0];
  775. else if (p3loc[1] > ylocMinMax[1]) p3loc[1] = ylocMinMax[1];
  776. time = yt->Eval(p3loc[1]);
  777. if ((idir0 == 0 && time > xyzloct[3].back()) || time < xyzloct[3].front()-dt || time > xyzloct[3].back()+dt)
  778. continue; // do not extrapolate curling track
  779. //p3loc[0] = tx.Eval(time);
  780. //p3loc[2] = tz.Eval(time);
  781. p3loc[0] = yx->Eval(p3loc[1]);
  782. p3loc[2] = yz->Eval(p3loc[1]);
  783. Double_t timh = TMath::Min (time,mite->first-0.001);
  784. timh = TMath::Max (timh,mitb->first);
  785. mithit = hitMap.lower_bound(timh);
  786. if (mithit != hitMap.end()) hitok = &mithit->second;
  787. else hitok = &hitMap.rbegin()->second;
  788. }
  789. lay0 = lay;
  790. secGeo->Local2Global(secGeo->Sector(padID),p3loc,p3);
  791. if (secGeo->Global2Local(p3, p30) < 0) continue; // cross-check
  792. MpdTpcHit *hitp = AddRawHit(nHits++, padID, p3, p3err, hitok->GetRefIndex(), id);
  793. //MpdTpcHit *hitp = new ((*fHitArray)[nHits++]) MpdTpcHit(*hitok); // copy constructor does not work
  794. hitp->SetLayer(lay);
  795. hitp->SetLocalPosition(p3loc); // point position
  796. hitp->SetPosition(p3);
  797. hitp->SetEnergyLoss(hitok->GetEnergyLoss());
  798. hitp->SetStep(hitok->GetStep());
  799. hitp->SetModular(1); // modular geometry flag
  800. hitp->SetLength(hitok->GetLength());
  801. }
  802. for (Int_t j = 0; j < 4; ++j) xyzloct[j].clear();
  803. isec0 = isec;
  804. idir0 = idir;
  805. delete yt;
  806. delete yx;
  807. delete yz;
  808. } // if (isec != isec0 || idir != idir0)
  809. if (isec == 30) break;
  810. hit.LocalPosition(p3loc);
  811. for (Int_t j = 0; j < 3; ++j) xyzloct[j].push_back(p3loc[j]);
  812. xyzloct[3].push_back(mit1->first);
  813. } // for (map<Double_t,MpdTpcHit>::iterator mit1 = hitMap.begin();
  814. } // for (map<Int_t,map<Double_t,Int_t> >::iterator mit = idmap.begin();
  815. }
  816. //___________________________________________________________________________
  817. ClassImp(MpdTpcHitProducer)