MpdEmcClusterFinderAZ.cxx 32 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class MpdEmcClusterFinderAZ
  7. // see MpdEmcClusterFinderAZ.h for details
  8. //
  9. // Environment:
  10. // Software developed for the MPD Detector at NICA.
  11. //
  12. // Author List:
  13. // Alexander Zinchenko LHEP, JINR, Dubna - 18-May-2016
  14. // Alexander Zinchenko LHEP, JINR, Dubna - 24-June-2018 - adapted for projective geometry
  15. //
  16. //-----------------------------------------------------------
  17. // This Class' Header ------------------
  18. #include "MpdEmcClusterFinderAZ.h"
  19. // Collaborating Class Headers --------
  20. #include "MpdTpc2dCluster.h"
  21. #include "MpdEmcDigit.h"
  22. #include "MpdEmcGeoPar.h"
  23. #include "MpdEmcGeoParams.h"
  24. #include "MpdEmcPoint.h"
  25. #include "MpdTpcHit.h"
  26. #include "MpdMCTrack.h"
  27. #include "FairRootManager.h"
  28. #include "TClonesArray.h"
  29. #include "TFile.h"
  30. #include "TMath.h"
  31. #include "TRandom.h"
  32. #include "TSpline.h"
  33. #include "TSystem.h"
  34. #include "TTree.h"
  35. // C/C++ Headers ----------------------
  36. #include <iostream>
  37. #include <iterator>
  38. #include <math.h>
  39. #include <set>
  40. #include <vector>
  41. using namespace std;
  42. //__________________________________________________________________________
  43. MpdEmcClusterFinderAZ::MpdEmcClusterFinderAZ()
  44. : FairTask("EMC Cluster finder AZ"), fPersistence(kFALSE), fThresh(0.005), fTimeWindow(0)
  45. {
  46. }
  47. //__________________________________________________________________________
  48. MpdEmcClusterFinderAZ::~MpdEmcClusterFinderAZ()
  49. {
  50. }
  51. //__________________________________________________________________________
  52. void MpdEmcClusterFinderAZ::FinishTask()
  53. {
  54. }
  55. //__________________________________________________________________________
  56. InitStatus MpdEmcClusterFinderAZ::Init()
  57. {
  58. // Create containers for digits
  59. //fEmcGeo = new MpdEmcGeoPar();
  60. fEmcGeo = new MpdEmcGeoParams();
  61. Int_t nSec = fEmcGeo->GetNsec() / 2; // number of long sectors
  62. set<Int_t> aaa;
  63. fDigiSet.assign(nSec,aaa);
  64. //for (Int_t i = 0; i < fgkNsec2; ++i) fDigiSet[i] = new set<Int_t> [nRows];
  65. // Containers for all+2 sectors
  66. //const Int_t nPhi = fEmcGeo->GetNsupMod() * fEmcGeo->GetNModInSuperModByPhi() * (nSec + 2);
  67. //const Int_t nZ = fEmcGeo->GetNrows() * fEmcGeo->GetNModInSuperModByZ() * 2;
  68. const Int_t nSecRows = fEmcGeo->GetNrows() / (nSec-1); // rows per wide sector
  69. const Int_t nPhi = fEmcGeo->GetNrows() + nSecRows * 2; // extra 2 wide sectors
  70. const Int_t nZ = fEmcGeo->GetNmod();
  71. // Fill starting rows of sectors
  72. fSecRows0.push_back(0);
  73. for (Int_t isec = 1; isec <= nSec; ++isec) {
  74. if (isec == 3 || isec == 7) fSecRows0.push_back(fSecRows0[isec-1]+nSecRows/2);
  75. else fSecRows0.push_back(fSecRows0[isec-1]+nSecRows);
  76. }
  77. // Fill ending rows of sectors
  78. fSecRows1.insert(nSecRows);
  79. for (Int_t isec = 1; isec < nSec; ++isec) {
  80. if (isec == 2 || isec == 6) fSecRows1.insert(*fSecRows1.rbegin()+nSecRows/2);
  81. else fSecRows1.insert(*fSecRows1.rbegin()+nSecRows);
  82. }
  83. vector<Double_t> dVec(nZ);
  84. vector<Int_t> iVec(nZ);
  85. fCharges.assign(nPhi,dVec);
  86. fFlags.assign(nPhi,iVec);
  87. fDigis.assign(nPhi,iVec);
  88. //Get ROOT Manager
  89. FairRootManager* ioman = FairRootManager::Instance();
  90. if (ioman == 0) {
  91. Error("MpdEmcClusterFinderAZ::Init","RootManager not instantiated!");
  92. return kERROR;
  93. }
  94. // Get input collection
  95. fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigit");
  96. if (fDigiArray == 0) {
  97. Error("EmcClusterFinderAZ::Init","Array of digits not found!");
  98. return kERROR;
  99. }
  100. fMcTrArray = (TClonesArray*) ioman->GetObject("MCTrack");
  101. if (!fMcTrArray) {
  102. Error("EmcClusterFinderAZ::Init","Array of MCTracks not found!");
  103. return kERROR;
  104. }
  105. // Create and register output array
  106. SetPersistence();
  107. fClusArray = new TClonesArray("MpdTpc2dCluster");
  108. ioman->Register("EmcCluster", "Emc", fClusArray, fPersistence);
  109. fHitArray = new TClonesArray("MpdTpcHit");
  110. ioman->Register("EmcRecPoint", "Emc", fHitArray, fPersistence);
  111. cout << fEmcGeo->GetPhiRow().size() << " " << fEmcGeo->GetThetaBox().size() << endl;
  112. Double_t theBeg = fEmcGeo->GetThetaBox()[0];
  113. theBeg = TMath::Max (theBeg,180-theBeg);
  114. cout << fEmcGeo->GetPhiRow()[0] << " " << *fEmcGeo->GetPhiRow().rbegin() << " "
  115. << theBeg << " " << *fEmcGeo->GetThetaBox().rbegin() << endl;
  116. return kSUCCESS;
  117. }
  118. //__________________________________________________________________________
  119. void MpdEmcClusterFinderAZ::Exec(Option_t* opt)
  120. {
  121. fClusArray->Delete();
  122. fHitArray->Delete();
  123. static const Int_t nSec = fEmcGeo->GetNsec() / 2; // number of EMC sectors
  124. //static const Int_t nTowSec = fEmcGeo->GetNsupMod() * fEmcGeo->GetNModInSuperModByPhi();
  125. // Clear digi containers
  126. for (Int_t i = 0; i < nSec; ++i) fDigiSet[i].clear();
  127. // Fill digi containers
  128. Int_t nDigis = fDigiArray->GetEntriesFast();
  129. cout << " Total number of digits: " << nDigis << endl;
  130. for (Int_t i = 0; i < nDigis; ++i) {
  131. MpdEmcDigit *digi = (MpdEmcDigit*) fDigiArray->UncheckedAt(i);
  132. //Int_t isec = digi->GetChanPhiId() / nTowSec;
  133. Int_t isec = digi->Sector();
  134. fDigiSet[isec].insert(i);
  135. }
  136. FillEmcInfo(); // fill info (matrices)
  137. Int_t nSum = 0;
  138. // Loop over sectors
  139. for (Int_t isec = 0; isec < nSec; ++isec) {
  140. if (fDigiSet[isec].size() == 0) continue;
  141. //cout << " Sector, digits: " << isec << " " << fDigiSet[isec].size() << endl;
  142. ProcessSector(isec);
  143. nSum += fDigiSet[isec].size();
  144. }
  145. cout << " Control sum: " << nSum << endl;
  146. // Find hits
  147. FindHits();
  148. //FindHitsLocMax();
  149. }
  150. //__________________________________________________________________________
  151. void MpdEmcClusterFinderAZ::FillEmcInfo()
  152. {
  153. // Fill EMC info (matrices)
  154. static const Int_t nPhi = fFlags.size(), nZ = fFlags[0].size(), nZov2 = nZ / 2;
  155. static const Int_t nSec = fEmcGeo->GetNsec() / 2; // number of EMC sectors
  156. //static const Int_t nTowSec = fEmcGeo->GetNsupMod() * fEmcGeo->GetNModInSuperModByPhi();
  157. static const Int_t nSecRows = nPhi / (nSec + 1);
  158. static const Double_t rmin = fEmcGeo->GetRmin();
  159. for (Int_t i = 0; i < nPhi; ++i) fFlags[i].assign(nZ,0);
  160. set<Int_t>::iterator it;
  161. for (Int_t isec = 0; isec < nSec; ++isec) {
  162. for (it = fDigiSet[isec].begin(); it != fDigiSet[isec].end(); ++it) {
  163. MpdEmcDigit *digi = (MpdEmcDigit*) fDigiArray->UncheckedAt(*it);
  164. // Apply threshold
  165. if (digi->GetE() < fThresh) continue;
  166. // Apply time window
  167. //Double_t dist = TMath::Sqrt (digi->GetZcenter() * digi->GetZcenter() + rmin * rmin);
  168. Double_t dist = rmin / TMath::Sin(digi->GetZcenter());
  169. Double_t dt = digi->GetTimeStamp() - dist / 30.; // c = 30 cm/ns
  170. if (fTimeWindow && (dt < -0.5 || dt > 2.0)) continue;
  171. Int_t iphi = digi->GetChanPhiId(), iz = digi->GetChanZId();
  172. Int_t ix = iphi + nSecRows; // offset by one sector
  173. //if (digi->Side() == 1) iz = -iz - 1;
  174. if (digi->Side() == 1) iz = -iz - 0; // shift by 1 already accounted for in ChanZId
  175. iz += nZov2;
  176. fCharges[ix][iz] = digi->GetE();
  177. fFlags[ix][iz] = 1;
  178. fDigis[ix][iz] = *it; // index of digit
  179. if (isec == 0) {
  180. // Fill extra sector (for edge effect)
  181. Int_t offset = fEmcGeo->GetNrows();
  182. Int_t ix1 = ix + offset;
  183. fCharges[ix1][iz] = digi->GetE();
  184. fFlags[ix1][iz] = 1;
  185. fDigis[ix1][iz] = *it;
  186. } else if (isec == nSec - 1) {
  187. Int_t ix1 = (ix + 2 * nSecRows) % nPhi;
  188. fCharges[ix1][iz] = digi->GetE();
  189. fFlags[ix1][iz] = 1;
  190. fDigis[ix1][iz] = *it;
  191. }
  192. }
  193. }
  194. }
  195. //__________________________________________________________________________
  196. void MpdEmcClusterFinderAZ::ProcessSector(Int_t isec)
  197. {
  198. // Process one sector
  199. const Int_t nSec = fEmcGeo->GetNsec() / 2; // number of EMC sectors
  200. //const Int_t nTowSec = fEmcGeo->GetNsupMod() * fEmcGeo->GetNModInSuperModByPhi();
  201. const Int_t nPhi = fFlags.size(), nZ = fFlags[0].size();
  202. static const Int_t nSecRows = nPhi / (nSec + 1);
  203. // Find (pre)clusters in the central sector
  204. Int_t nclus0 = fClusArray->GetEntriesFast(), nclus = nclus0;
  205. //Int_t phiBeg = (isec + 1) * nTowSec, phiEnd = phiBeg + nTowSec;
  206. Int_t phiBeg = fSecRows0[isec] + nSecRows, phiEnd = fSecRows0[isec+1] + nSecRows;
  207. for (Int_t iphi = phiBeg; iphi < phiEnd; ++iphi) {
  208. for (Int_t iz = 0; iz < nZ; ++iz) {
  209. if (fFlags[iphi][iz] <= 0) continue;
  210. // New cluster
  211. MpdTpc2dCluster* clus = new ((*fClusArray)[nclus++]) MpdTpc2dCluster(0, isec);
  212. clus->Insert(fDigis[iphi][iz], 0, iphi-nSecRows, iz, fCharges[iphi][iz]);
  213. clus->SetID(fDigis[iphi][iz]); // trackID
  214. clus->SetErrY(fCharges[iphi][iz]); // ADC counts
  215. fFlags[iphi][iz] = -1;
  216. //if (isec == 0) fFlags[iphi+nPhi-2*nTowSec][iz] = -1;
  217. //else if (isec == nSec-1) fFlags[iphi-nPhi+2*nTowSec][iz] = -1;
  218. if (isec == 0) fFlags[iphi+nPhi-2*nSecRows][iz] = -1;
  219. else if (isec == nSec-1) fFlags[(iphi+2*nSecRows)%nPhi][iz] = -1;
  220. for (Int_t ip = -1; ip < 2; ++ip) {
  221. for (Int_t it = -1; it < 2; ++it) {
  222. //if (it == ip) continue;
  223. if (it == 0 && ip == 0) continue;
  224. Int_t ip1 = iphi + ip, iz1 = iz + it;
  225. if (ip1 < 0 || ip1 >= nPhi) continue;
  226. if (iz1 < 0 || iz1 >= nZ) continue;
  227. if (fFlags[ip1][iz1] <= 0) continue;
  228. NextPixel(clus, ip1, iz1);
  229. }
  230. }
  231. }
  232. }
  233. //cout << " Found preclusters: " << nclus - nclus0 << endl;
  234. }
  235. //__________________________________________________________________________
  236. void MpdEmcClusterFinderAZ::NextPixel(MpdTpc2dCluster* clus, Int_t iphi, Int_t iz)
  237. {
  238. // Add next pixel to the cluster
  239. const Int_t nPhi = fFlags.size(), nZ = fFlags[0].size();
  240. const Int_t nSec = fEmcGeo->GetNsec() / 2; // number of EMC sectors
  241. //const Int_t nTowSec = fEmcGeo->GetNsupMod() * fEmcGeo->GetNModInSuperModByPhi();
  242. const Int_t nSecRows = nPhi / (nSec + 1);
  243. clus->Insert(fDigis[iphi][iz], clus->Row(), iphi-nSecRows, iz, fCharges[iphi][iz]);
  244. clus->SetID (TMath::Min(clus->ID(),fDigis[iphi][iz])); // min trackID
  245. clus->SetErrY (TMath::Max(Double_t(clus->GetErrY()),fCharges[iphi][iz])); // max ADC counts
  246. fFlags[iphi][iz] = -1;
  247. // Edges
  248. if (iphi < nSecRows) fFlags[iphi+nPhi-2*nSecRows][iz] = -1; // last sector
  249. else if (iphi >= nPhi-nSecRows) fFlags[iphi-nPhi+2*nSecRows][iz] = -1; // first sector
  250. else if (iphi >= nPhi-2*nSecRows) fFlags[iphi-nPhi+2*nSecRows][iz] = -1; // last sector
  251. else if (iphi < 2*nSecRows) fFlags[iphi+nPhi-2*nSecRows][iz] = -1; // first sector
  252. for (Int_t ip = -1; ip < 2; ++ip) {
  253. for (Int_t it = -1; it < 2; ++it) {
  254. //if (TMath::Abs(it) == TMath::Abs(ip)) continue;
  255. Int_t ip1 = iphi + ip, it1 = iz + it;
  256. if (ip1 < 0 || ip1 >= nPhi) continue;
  257. if (it1 < 0 || it1 >= nZ) continue;
  258. if (fFlags[ip1][it1] <= 0) continue;
  259. NextPixel(clus, ip1, it1);
  260. }
  261. }
  262. }
  263. //__________________________________________________________________________
  264. /*
  265. void MpdEmcClusterFinderAZ::FindHits()
  266. {
  267. // Reconstruct hits (one hit per precluster)
  268. const Int_t nTowSec = fEmcGeo->GetNsupMod() * fEmcGeo->GetNModInSuperModByPhi();
  269. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  270. Int_t nclus = fClusArray->GetEntriesFast();
  271. for (Int_t iclus = 0; iclus < nclus; ++iclus) {
  272. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  273. Int_t nDigis = clus->NDigits();
  274. Double_t phiMean = 0, zMean = 0;
  275. for (Int_t idig = 0; idig < nDigis; ++idig) {
  276. phiMean += clus->Col(idig) * clus->Adc(idig);
  277. zMean += clus->Bkt(idig) * clus->Adc(idig);
  278. }
  279. phiMean /= clus->GetADC();
  280. phiMean -= nTowSec;
  281. zMean /= clus->GetADC();
  282. clus->SetX(phiMean);
  283. clus->SetZ(zMean);
  284. phiMean *= fEmcGeo->GetAngleOfModule() * TMath::DegToRad();
  285. //phiMean += clus->GetSect() * fEmcGeo->GetAngleOfSector() * TMath::DegToRad();
  286. clus->SetGlobX(phiMean);
  287. zMean -= (fFlags[0].size() / 2 - 0.5);
  288. zMean *= fEmcGeo->GetLengthOfModuleByZ();
  289. clus->SetGlobZ(zMean);
  290. }
  291. }
  292. */
  293. //__________________________________________________________________________
  294. /*
  295. void MpdTpcClusterFinderAZ::FindHitsLocMax()
  296. {
  297. // Reconstruct hits (find local maxima)
  298. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  299. Int_t nclus = fClusArray->GetEntriesFast(), ihit = 0;
  300. for (Int_t iclus = 0; iclus < nclus; ++iclus) {
  301. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  302. Int_t nDigis = clus->NDigits();
  303. memset(fFlags, 0, sizeof(fFlags[0][0]) * fgkNpads * fgkNtimes);
  304. for (Int_t idig = 0; idig < nDigis; ++idig) {
  305. Int_t ipad = clus->Col(idig);
  306. Int_t itime = clus->Bkt(idig);
  307. fCharges[ipad][itime] = clus->Adc(idig);
  308. fFlags[ipad][itime] = 1;
  309. fDigis[ipad][itime] = idig;
  310. }
  311. // Exclude pads which are not local maxima
  312. for (Int_t idig = 0; idig < nDigis; ++idig) {
  313. Int_t ipad = clus->Col(idig);
  314. Int_t itime = clus->Bkt(idig);
  315. for (Int_t ip = -1; ip < 2; ++ip) {
  316. for (Int_t it = -1; it < 2; ++it) {
  317. //if (TMath::Abs(it) == TMath::Abs(ip)) continue; // exclude diagonals
  318. if (it == 0 && ip == 0) continue;
  319. Int_t ip1 = ipad + ip, it1 = itime + it;
  320. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  321. if (it1 < 0 || it1 >= fgkNtimes) continue;
  322. if (fFlags[ip1][it1] == 0) continue;
  323. if (clus->Adc(idig) < fCharges[ip1][it1]) { fFlags[ipad][itime] = -1; break; }
  324. }
  325. }
  326. }
  327. multimap<Double_t,Int_t> localMax;
  328. for (Int_t idig = 0; idig < nDigis; ++idig) {
  329. Int_t ipad = clus->Col(idig);
  330. Int_t itime = clus->Bkt(idig);
  331. if (fFlags[ipad][itime] <= 0) continue;
  332. localMax.insert(pair<Double_t,Int_t>(clus->Adc(idig),idig));
  333. cout << clus->Col(idig) << " " << clus->Bkt(idig) << " " << clus->Adc(idig) << endl;
  334. }
  335. cout << " Local max: " << localMax.size() << endl;
  336. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin(), rit1 = rit;
  337. vector<Int_t> vecDig;
  338. for ( ; rit != localMax.rend(); ++rit) {
  339. Int_t idig = rit->second;
  340. Int_t ipad = clus->Col(idig);
  341. Int_t itime = clus->Bkt(idig);
  342. if (fFlags[ipad][itime] <= 0) continue; // merged local max
  343. vecDig.clear();
  344. Double_t padMean = 0, timeMean = 0, adcTot = 0;
  345. for (Int_t ip = -1; ip < 2; ++ip) {
  346. for (Int_t it = -1; it < 2; ++it) {
  347. //if (TMath::Abs(it) == TMath::Abs(ip)) continue;
  348. Int_t ip1 = ipad + ip, it1 = itime + it;
  349. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  350. if (it1 < 0 || it1 >= fgkNtimes) continue;
  351. if (fFlags[ip1][it1] == 0) continue;
  352. if (fFlags[ip1][it1] > 0) fFlags[ip1][it1] = -1; // merge neighbour local max
  353. padMean += ip1 * fCharges[ip1][it1];
  354. timeMean += it1 * fCharges[ip1][it1];
  355. adcTot += fCharges[ip1][it1];
  356. vecDig.push_back(fDigis[ip1][it1]);
  357. }
  358. }
  359. padMean /= adcTot;
  360. timeMean /= adcTot;
  361. Double_t xloc = fSecGeo->Pad2Xloc(padMean,clus->Row());
  362. Int_t padID = fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  363. Double_t yloc = fSecGeo->LocalPadPosition(padID).Y();
  364. Double_t zloc = fSecGeo->TimeBin2Z(timeMean);
  365. p3loc.SetXYZ(xloc, yloc, zloc);
  366. // Apply corrections
  367. TVector3 p3errCor(p3err);
  368. CorrectReco(p3loc, p3errCor, clus->GetNumPads(), adcTot);
  369. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  370. if (clus->GetSect() >= fSecGeo->NofSectors()) p3glob[2] = -p3glob[2];
  371. MpdTpcHit* hit = new ((*fHitArray)[ihit++]) MpdTpcHit(padID, p3glob, p3errCor, iclus);
  372. hit->SetLayer(clus->Row());
  373. hit->SetLocalPosition(p3loc); // point position
  374. hit->SetEnergyLoss(adcTot);
  375. hit->SetStep(0.0);
  376. hit->SetModular(1); // modular geometry flag
  377. hit->SetPad(Int_t(padMean));
  378. hit->SetBin(Int_t(timeMean));
  379. // !!! Warning: FairLinks are not persistent !!!
  380. //hit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  381. Int_t ndig = vecDig.size();
  382. for (Int_t idig = 0; idig < ndig; ++idig)
  383. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, clus->Sec(vecDig[idig]))); // trackID stored in Sec
  384. //cout << hit->GetNLinks();
  385. } // for ( ; rit != localMax.rend();
  386. } // for (Int_t iclus = 0; iclus < nclus;
  387. }
  388. */
  389. //__________________________________________________________________________
  390. void MpdEmcClusterFinderAZ::FindHits()
  391. {
  392. // Reconstruct hits (find local maxima)
  393. static const Int_t nPhi = fFlags.size(), nZ = fFlags[0].size();
  394. //static const Double_t dphiSec = fEmcGeo->GetAngleOfSector() * TMath::DegToRad();
  395. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  396. Int_t nclus = fClusArray->GetEntriesFast(), ihit = 0;
  397. for (Int_t iclus = 0; iclus < nclus; ++iclus) {
  398. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  399. Int_t nDigis = clus->NDigits();
  400. for (Int_t i = 0; i < nPhi; ++i) fFlags[i].assign(nZ,0);
  401. Int_t isec = clus->GetSect();
  402. Int_t ishift = 0;
  403. //if (isec == 0) ishift = nPhi / 30; // just to shift index up from 0
  404. if (isec == 0) ishift = 50; // just to shift index up from 0
  405. for (Int_t idig = 0; idig < nDigis; ++idig) {
  406. Int_t iphi = clus->Col(idig) + ishift;
  407. Int_t iz = clus->Bkt(idig);
  408. fCharges[iphi][iz] = clus->Adc(idig);
  409. fFlags[iphi][iz] = 1;
  410. //fDigis[iphi][iz] = idig;
  411. fDigis[iphi][iz] = clus->Sec(idig);
  412. }
  413. // Exclude pads which are not local maxima
  414. for (Int_t idig = 0; idig < nDigis; ++idig) {
  415. Int_t iphi = clus->Col(idig) + ishift;
  416. Int_t iz = clus->Bkt(idig);
  417. for (Int_t ip = -1; ip < 2; ++ip) {
  418. for (Int_t it = -1; it < 2; ++it) {
  419. //if (TMath::Abs(it) == TMath::Abs(ip)) continue; // exclude diagonals
  420. if (it == 0 && ip == 0) continue;
  421. Int_t ip1 = iphi + ip, it1 = iz + it;
  422. if (ip1 < 0 || ip1 >= nPhi) continue;
  423. if (it1 < 0 || it1 >= nZ) continue;
  424. if (fFlags[ip1][it1] == 0) continue;
  425. if (clus->Adc(idig) < fCharges[ip1][it1]) { fFlags[iphi][iz] = -1; break; }
  426. }
  427. }
  428. }
  429. multimap<Double_t,Int_t> localMax;
  430. for (Int_t idig = 0; idig < nDigis; ++idig) {
  431. Int_t iphi = clus->Col(idig) + ishift;
  432. Int_t iz = clus->Bkt(idig);
  433. if (fFlags[iphi][iz] <= 0) continue;
  434. localMax.insert(pair<Double_t,Int_t>(clus->Adc(idig),idig));
  435. //cout << clus->Col(idig) << " " << clus->Bkt(idig) << " " << clus->Adc(idig) << endl;
  436. }
  437. //cout << " Local max: " << clus->GetSect() << " " << clus->Row() << " " << localMax.size() << endl;
  438. // Remove local maxima not separated by valleys - Peak and Valley
  439. //Int_t nLocMax0 = localMax.size();
  440. if (localMax.size() > 1) PeakAndValley(clus, localMax, ishift);
  441. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin();
  442. //vector<Int_t> vecDig;
  443. map<Int_t,Double_t> mapIdQ;
  444. for ( ; rit != localMax.rend(); ++rit) {
  445. Int_t idig = rit->second;
  446. Int_t iphi = clus->Col(idig) + ishift;
  447. Int_t iz = clus->Bkt(idig);
  448. if (fFlags[iphi][iz] <= 0) continue; // merged local max
  449. //vecDig.clear();
  450. mapIdQ.clear();
  451. set<pair<Int_t,Int_t> > pixels;
  452. Double_t phiMean = 0, zMean = 0, adcTot = 0, sum2t = 0, sum2p = 0;
  453. // Process simple cluster (only 1 local max)
  454. if (localMax.size() == 1) {
  455. //AZ if (nLocMax0 == 1) {
  456. for (Int_t idig1 = 0; idig1 < nDigis; ++idig1) {
  457. Int_t ip = clus->Col(idig1) + ishift;
  458. Int_t it = clus->Bkt(idig1);
  459. phiMean += ip * fCharges[ip][it];
  460. zMean += it * fCharges[ip][it];
  461. adcTot += fCharges[ip][it];
  462. //vecDig.push_back(fDigis[ip][it]);
  463. //Int_t id = clus->Sec(fDigis[ip][it]);
  464. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = fCharges[ip][it];
  465. //else mapIdQ[id] = mapIdQ[id] + fCharges[ip][it];
  466. //else mapIdQ[id] = TMath::Max (mapIdQ[id], fCharges[ip][it]);
  467. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = 1;
  468. //else mapIdQ[id] = mapIdQ[id] + 1;
  469. MpdEmcDigit *dig = (MpdEmcDigit*) fDigiArray->UncheckedAt(fDigis[ip][it]);
  470. map<Int_t,Float_t> contrib = dig->GetContrib();
  471. RedoId(contrib);
  472. for (map<Int_t,Float_t>::iterator mit = contrib.begin(); mit != contrib.end(); ++mit) {
  473. if (mapIdQ.find(mit->first) == mapIdQ.end()) mapIdQ[mit->first] = mit->second;
  474. else mapIdQ[mit->first] = mapIdQ[mit->first] + mit->second;
  475. }
  476. sum2t += it * it * fCharges[ip][it];
  477. sum2p += ip * ip * fCharges[ip][it];
  478. }
  479. } else {
  480. // Process complex cluster - start from maximum and go ouward (up tp 5 steps),
  481. // adding pixels with repectively lower charges
  482. for (Int_t idirp = -1; idirp < 2; idirp += 2) {
  483. for (Int_t ip = 0; ip < 5; ++ip) {
  484. if (idirp > 0 && ip == 0) continue;
  485. Int_t ipsign = ip * idirp;
  486. Int_t ip1 = iphi + ipsign;
  487. if (ip1 < 0 || ip1 >= nPhi) break;
  488. for (Int_t idirt = -1; idirt < 2; idirt += 2) {
  489. for (Int_t it = 0; it < 5; ++it) {
  490. if (idirt > 0 && it == 0) continue;
  491. Int_t itsign = it * idirt;
  492. Int_t it1 = iz + itsign;
  493. if (it1 < 0 || it1 >= nZ) break;
  494. if (fFlags[ip1][it1] == 0) continue;
  495. Int_t add = 1;
  496. if (ip || it) {
  497. if (fFlags[ip1][it1] > 0) continue; // case when 2 local max next to each other on 1 diagonal
  498. // Check gradient
  499. add = 0;
  500. Int_t ipprev = ipsign, itprev = itsign;
  501. if (it) {
  502. itprev -= idirt;
  503. Int_t it10 = iz + itprev;
  504. if (it10 >= 0 && it10 < nZ && fFlags[ip1][it10] != 0) {
  505. if (pixels.find(pair<Int_t,Int_t>(ip1,it10)) != pixels.end()
  506. && fCharges[ip1][it1] <= fCharges[ip1][it10]) add = 1;
  507. }
  508. }
  509. if (add == 0 && ip) {
  510. ipprev -= idirp;
  511. Int_t ip10 = iphi + ipprev;
  512. if (ip10 >= 0 && ip10 < nPhi && fFlags[ip10][it1] != 0) {
  513. if (pixels.find(pair<Int_t,Int_t>(ip10,it1)) != pixels.end()
  514. && fCharges[ip1][it1] <= fCharges[ip10][it1]) add = 1;
  515. }
  516. }
  517. }
  518. if (!add) break;
  519. phiMean += ip1 * fCharges[ip1][it1];
  520. zMean += it1 * fCharges[ip1][it1];
  521. adcTot += fCharges[ip1][it1];
  522. //vecDig.push_back(fDigis[ip1][it1]);
  523. //Int_t id = clus->Sec(fDigis[ip1][it1]);
  524. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = fCharges[ip1][it1];
  525. //else mapIdQ[id] = mapIdQ[id] + fCharges[ip1][it1];
  526. //else mapIdQ[id] = TMath::Max (mapIdQ[id], fCharges[ip1][it1]);
  527. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = 1;
  528. //else mapIdQ[id] = mapIdQ[id] + 1; // number of digits with the same ID
  529. MpdEmcDigit *dig = (MpdEmcDigit*) fDigiArray->UncheckedAt(fDigis[ip1][it1]);
  530. map<Int_t,Float_t> contrib = dig->GetContrib();
  531. RedoId(contrib);
  532. for (map<Int_t,Float_t>::iterator mit = contrib.begin(); mit != contrib.end(); ++mit) {
  533. if (mapIdQ.find(mit->first) == mapIdQ.end()) mapIdQ[mit->first] = mit->second;
  534. else mapIdQ[mit->first] = mapIdQ[mit->first] + mit->second;
  535. }
  536. pixels.insert(pair<Int_t,Int_t>(ip1,it1));
  537. sum2t += it1 * it1 * fCharges[ip1][it1];
  538. sum2p += ip1 * ip1 * fCharges[ip1][it1];
  539. }
  540. }
  541. }
  542. }
  543. }
  544. phiMean /= adcTot;
  545. zMean /= adcTot;
  546. //phiMean += (0.5 - ishift); //
  547. phiMean += (0.0 - ishift); //
  548. //zMean += 0.5; //
  549. Double_t xloc = phiMean; //fSecGeo->Pad2Xloc(phiMean,clus->Row());
  550. Int_t padID = 0; //fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  551. Double_t yloc = 0; //fSecGeo->LocalPadPosition(padID).Y();
  552. Double_t zloc = zMean; //fSecGeo->TimeBin2Z(zMean);
  553. p3loc.SetXYZ(xloc, yloc, zloc);
  554. Double_t rmsZ = TMath::Sqrt (sum2t/adcTot - zMean*zMean);
  555. Double_t rmsX = TMath::Sqrt (sum2p/adcTot - phiMean*phiMean);
  556. //p3err[1] = fSecGeo->TimeBin2Z(rms); // to pass the value
  557. //p3err[1] = rms;
  558. //cout << " Result: " << nLocMax0 << " " << pixels.size() << " " << xloc << " " << zloc << endl;
  559. // Apply corrections
  560. TVector3 p3errCor(p3err);
  561. /*
  562. CorrectReco(p3loc, p3errCor, clus->GetNumPads(), adcTot);
  563. if (clus->GetSect() >= fSecGeo->NofSectors()) p3loc[2] = -p3loc[2];
  564. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  565. // Dip angle correction (for interaction point with Z = 0)
  566. Double_t dip = TMath::ATan2 (TMath::Abs(p3glob[2]),p3glob.Pt()) * TMath::RadToDeg();
  567. p3errCor[2] = TMath::Max (p3errCor[2], 0.116 - 0.0046 * dip + 0.00015 * dip * dip);
  568. p3errCor[2] = TMath::Min (p3errCor[2], 0.5);
  569. // Correct Z-coordinate for rows = 0 and 27
  570. if (clus->Row() == 0 || clus->Row() == 27) {
  571. Double_t zcor = 0;
  572. if (clus->Row() == 0) zcor = -0.011 + 0.002 * dip;
  573. else zcor = -0.029 + 0.005 * dip + 3.886e-5 * dip * dip;
  574. zcor = TMath::Max (zcor, 0.);
  575. zcor = TMath::Min (zcor, 0.6);
  576. p3loc[2] -= TMath::Sign (zcor, p3loc[2]);
  577. p3glob[2] -= TMath::Sign (zcor, p3glob[2]);
  578. }
  579. */
  580. //zMean -= (fFlags[0].size() / 2 - 0.5);
  581. //AZ zMean -= (fFlags[0].size() / 2 - 0.0);
  582. //FIXME phiMean *= fEmcGeo->GetAngleOfModule() * TMath::DegToRad();
  583. //zMean *= fEmcGeo->GetLengthOfModuleByZ(); // FIXME - Z-position calculation
  584. GetPhiTheta(phiMean, zMean); // get angular coordinates phi and theta at (near) the inner faces of towers
  585. p3glob.SetXYZ(phiMean,0,zMean);
  586. MpdTpcHit* hit = new ((*fHitArray)[ihit++]) MpdTpcHit(padID, p3glob, p3errCor, iclus);
  587. //hit->SetLayer(clus->GetSect()); // sector number
  588. set<Int_t>::iterator sit = fSecRows1.upper_bound(xloc);
  589. hit->SetLayer(std::distance(fSecRows1.begin(),sit)); // sector number
  590. hit->SetLocalPosition(p3loc); // point position
  591. hit->SetEnergyLoss(adcTot);
  592. //Int_t ireg = 0; //(clus->Row() < fSecGeo->NofRowsReg(0)) ? 0 : 1;
  593. Double_t step = 1; //fSecGeo->PadHeight(ireg);
  594. hit->SetStep(step);
  595. hit->SetModular(1); // modular geometry flag
  596. hit->SetPad(Int_t(phiMean));
  597. hit->SetBin(Int_t(zMean));
  598. hit->SetRMS(rmsX, 0);
  599. hit->SetRMS(rmsZ, 1);
  600. //if (nLocMax0 == 1) hit->SetNdigits(nDigis);
  601. if (pixels.size() == 0) hit->SetNdigits(nDigis);
  602. else hit->SetNdigits(pixels.size());
  603. // !!! Warning: FairLinks are not persistent !!!
  604. //hit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  605. vector<Int_t> vecDig;
  606. // Store maximum 5 track IDs with the highest charge contribution
  607. multimap<Double_t,Int_t> mapDig;
  608. for (map<Int_t,Double_t>::iterator mit = mapIdQ.begin(); mit != mapIdQ.end(); ++mit)
  609. mapDig.insert(pair<Double_t,Int_t>(-mit->second,mit->first));
  610. multimap<Double_t,Int_t>::iterator mit = mapDig.begin();
  611. while (mit != mapDig.end() && vecDig.size() < 5) {
  612. //while (mit != mapDig.end() && vecDig.size() < 1) {
  613. vecDig.push_back(mit->second);
  614. ++mit;
  615. }
  616. Int_t ndig = vecDig.size();
  617. for (Int_t idig1 = 0; idig1 < ndig; ++idig1) {
  618. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, vecDig[idig1], idig1)); // weight = idig
  619. hit->AddID(vecDig[idig1]);
  620. }
  621. //cout << " EMC IDs: " << hit->GetNLinks() << " " << *(vecDig.begin()) << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(0).GetIndex() << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(hit->GetNLinks()-1).GetIndex() << endl;
  622. } // for ( ; rit != localMax.rend();
  623. } // for (Int_t iclus = 0; iclus < nclus;
  624. }
  625. //*/
  626. //__________________________________________________________________________
  627. void MpdEmcClusterFinderAZ::GetPhiTheta(Double_t &phi, Double_t &theta)
  628. {
  629. // Convert COG in units of bins to angles
  630. static Int_t first = 1, offset = 0;
  631. static TSpline3 *phiS, *theS;
  632. if (first) {
  633. // Get phi and theta angles of the tower centers at their inner face
  634. first = 0;
  635. const vector<Double_t> &phis = fEmcGeo->GetPhiRow();
  636. const vector<Double_t> &thes = fEmcGeo->GetThetaBox();
  637. const vector<Double_t> &rhos = fEmcGeo->GetRhoCenterBox();
  638. const vector<Double_t> &zs = fEmcGeo->GetZCenterBox();
  639. Int_t nphi = phis.size();
  640. // Offset due to the fact that the 1'st sector starts at phi = -Phi_sec/2;
  641. offset = nphi / (fEmcGeo->GetNsec()/2 - 1) / 2;
  642. Double_t *phia = new Double_t [nphi];
  643. Double_t *ind = new Double_t [nphi];
  644. for (Int_t j = 0; j < nphi; ++j) {
  645. phia[j] = phis[j];
  646. ind[j] = j;
  647. }
  648. phiS = new TSpline3("grs",ind,phia,nphi); // phi vs ind
  649. delete [] phia;
  650. delete [] ind;
  651. Int_t nthe = thes.size();
  652. Double_t *the = new Double_t [nthe];
  653. Double_t *ind1 = new Double_t [nthe];
  654. Double_t height = fEmcGeo->GetLengthBox(); // tower half-height
  655. for (Int_t j = nthe-1; j >= 0; --j) {
  656. Double_t rho = rhos[j];
  657. Double_t z = zs[j];
  658. Double_t theta1 = thes[j];
  659. if (j < nthe-1 && thes[j] <= thes[j+1]+0.1) theta1 = 180 - theta1;
  660. Double_t costhe = TMath::Cos(theta1*TMath::DegToRad());
  661. Double_t sinthe = TMath::Sin(theta1*TMath::DegToRad());
  662. rho -= height * sinthe;
  663. z -= height * costhe;
  664. the[j] = TMath::ATan2(rho,z) * TMath::RadToDeg();
  665. ind1[j] = j; // - nthe/2;
  666. }
  667. theS = new TSpline3("grs1",ind1,the,nthe); // theta vs ind
  668. delete [] the;
  669. delete [] ind1;
  670. }
  671. phi = phiS->Eval(phi-offset);
  672. if (phi > 180) phi -= 360;
  673. theta = theS->Eval(theta);
  674. }
  675. //__________________________________________________________________________
  676. void MpdEmcClusterFinderAZ::PeakAndValley(const MpdTpc2dCluster* clus, multimap<Double_t,Int_t> &localMax, Int_t ishift)
  677. {
  678. // Apply peak-and-valley cuts to remove some local maxima
  679. //*
  680. static const Double_t ratio = 1.3;
  681. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin(), rit1;
  682. ++rit;
  683. for ( ; rit != localMax.rend(); ++rit) {
  684. Int_t idig = rit->second;
  685. Int_t ipad = clus->Col(idig) + ishift;
  686. Int_t itime = clus->Bkt(idig);
  687. if (fFlags[ipad][itime] <= 0) continue; // merged local max
  688. // Loop over higher peaks
  689. rit1 = localMax.rbegin();
  690. for ( ; rit1 != rit; ++rit1) {
  691. Int_t idig0 = rit1->second;
  692. Int_t ipad0 = clus->Col(idig0) + ishift;
  693. Int_t itime0 = clus->Bkt(idig0);
  694. if (fFlags[ipad0][itime0] <= 0) continue; // merged local max
  695. Int_t dpad = ipad - ipad0, dt = itime - itime0;
  696. Int_t i0 = itime0, i1 = itime, j0 = ipad0, j1 = ipad, intime = 1;
  697. if (TMath::Abs(dpad) > TMath::Abs(dt)) i0 = ipad0, i1 = ipad, j0 = itime0, j1 = itime, intime = 0;
  698. Int_t stepi = TMath::Sign(1,i1-i0), stepj = TMath::Sign(1,j1-j0);
  699. //Int_t valOk = 0;
  700. Int_t merge = 0;
  701. //if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) merge = 1;
  702. //else {
  703. if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) {
  704. if ( (fFlags[ipad-dpad][itime] && fCharges[ipad-dpad][itime] > fCharges[ipad][itime] / ratio) ||
  705. (fFlags[ipad][itime-dt] && fCharges[ipad][itime-dt] > fCharges[ipad][itime] / ratio) ) merge = 1;
  706. } else {
  707. for (Int_t ii = i0 + stepi; ii != i1; ii += stepi) {
  708. merge = 0;
  709. for (Int_t jj = j0; jj != j1 + stepj; jj += stepj) {
  710. if (TMath::Abs(jj-j0) > TMath::Abs(ii-i0)) continue;
  711. if (TMath::Abs(jj-j1) > TMath::Abs(ii-i1)) continue;
  712. Int_t ip = ii, it = jj;
  713. if (intime) ip = jj, it = ii;
  714. //if (fCharges[ip][it] < fCharges[ipad][itime] / ratio) { valOk = 1; break; }
  715. if (fFlags[ip][it] && fCharges[ip][it] > fCharges[ipad][itime] / ratio) { merge = 1; break; }
  716. }
  717. //if (valOk) break;
  718. if (!merge) break;
  719. }
  720. }
  721. //if (!valOk) fFlags[ipad][itime] = -fFlags[ipad][itime]; // not deep enough valley
  722. if (merge) fFlags[ipad][itime] = -fFlags[ipad][itime]; // not deep enough valley
  723. }
  724. }
  725. // Remove failed peaks
  726. multimap<Double_t,Int_t>::iterator it = localMax.begin();
  727. //for ( ; it != localMax.end(); ++it) cout << it->second << " ";
  728. //cout << endl;
  729. it = localMax.begin();
  730. for ( ; it != localMax.end(); ++it) {
  731. Int_t idig = it->second;
  732. Int_t ipad = clus->Col(idig) + ishift;
  733. Int_t itime = clus->Bkt(idig);
  734. if (fFlags[ipad][itime] > 0) continue;
  735. //cout << " Before: " << idig << " " << itime << " " << ipad << " " << it->first << ", ";
  736. localMax.erase(it);
  737. //cout << " After: " << it->first << endl;
  738. }
  739. //*/
  740. }
  741. //__________________________________________________________________________
  742. void MpdEmcClusterFinderAZ::RedoId(map<Int_t,Float_t>& contrib)
  743. {
  744. // Redo track ID numbering
  745. // Take IDs of particles produced outside EMC
  746. static const Double_t rmin = fEmcGeo->GetRmin(), rmax = fEmcGeo->GetRmax();
  747. TVector3 vert;
  748. map<Int_t,Float_t> copy(contrib);
  749. contrib.clear();
  750. for (map<Int_t,Float_t>::iterator it = copy.begin(); it != copy.end(); ++it) {
  751. Int_t id = it->first, idm = -1;
  752. MpdMCTrack *mctr = (MpdMCTrack*) fMcTrArray->UncheckedAt(id);
  753. mctr->GetStartVertex(vert);
  754. Double_t rvert = vert.Pt();
  755. while (rvert > rmin && rvert < rmax) {
  756. // Born inside EMC - find ancestor
  757. idm = mctr->GetMotherId();
  758. if (idm < 0) break;
  759. mctr = (MpdMCTrack*) fMcTrArray->UncheckedAt(idm);
  760. mctr->GetStartVertex(vert);
  761. rvert = vert.Pt();
  762. }
  763. if (idm >= 0) id = idm;
  764. if (contrib.find(id) == contrib.end()) contrib[id] = it->second;
  765. else contrib[id] += it->second;
  766. }
  767. }
  768. //__________________________________________________________________________
  769. ClassImp(MpdEmcClusterFinderAZ)