MpdTpcClusterFinderMlem.cxx 56 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class MpdTpcClusterFinderMlem
  7. // see MpdTpcClusterFinderMlem.h for details
  8. //
  9. // Environment:
  10. // Software developed for the MPD Detector at NICA.
  11. //
  12. // Author List:
  13. // Alexandr Zinchenko LHEP, JINR, Dubna - 11-January-2016
  14. //
  15. //-----------------------------------------------------------
  16. // This Class' Header ------------------
  17. #include "MpdTpcClusterFinderMlem.h"
  18. // Collaborating Class Headers --------
  19. #include "MpdTpc2dCluster.h"
  20. #include "MpdTpcDigit.h"
  21. #include "MpdTpcHit.h"
  22. #include "MpdTpcSectorGeo.h"
  23. #include "TpcPoint.h"
  24. #include "FairRootManager.h"
  25. #include <TClonesArray.h>
  26. #include <TFile.h>
  27. #include <TH2.h>
  28. #include <TMath.h>
  29. #include <TMatrixD.h>
  30. #include <TROOT.h>
  31. #include <TSystem.h>
  32. // C/C++ Headers ----------------------
  33. #include <iostream>
  34. #include <math.h>
  35. #include <set>
  36. #include <vector>
  37. static clock_t tStartMlem = 0;
  38. static clock_t tFinishMlem = 0;
  39. static clock_t tAllMlem = 0;
  40. using namespace std;
  41. //static Bool_t SortPix(const MpdTpcClusterFinderMlem::pixel i,
  42. // const MpdTpcClusterFinderMlem::pixel j); // sorting function
  43. //__________________________________________________________________________
  44. MpdTpcClusterFinderMlem::MpdTpcClusterFinderMlem()
  45. : FairTask("TPC Cluster finder Mlem"), fPersistence(kFALSE)
  46. {
  47. /*
  48. std::string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  49. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  50. fGas= new TpcGas(tpcGasFile, 130);
  51. std::cout<<*fGas<<std::endl;
  52. */
  53. }
  54. //__________________________________________________________________________
  55. MpdTpcClusterFinderMlem::~MpdTpcClusterFinderMlem()
  56. {
  57. //delete fGas;
  58. }
  59. //__________________________________________________________________________
  60. void MpdTpcClusterFinderMlem::FinishTask()
  61. {
  62. for (Int_t i = 0; i < fgkNsec2; ++i) {
  63. //fPrimArray[i]->Delete();
  64. delete [] fDigiSet[i];
  65. }
  66. cout << "MLEM cluster finder work time = " << ((Float_t)tAllMlem) / CLOCKS_PER_SEC << endl;
  67. }
  68. //__________________________________________________________________________
  69. InitStatus MpdTpcClusterFinderMlem::Init()
  70. {
  71. // Create containers for digits
  72. fSecGeo = MpdTpcSectorGeo::Instance();
  73. //Int_t nRows = MpdTpcSectorGeo::Instance()->NofRows();
  74. Int_t nRows = fSecGeo->NofRows();
  75. for (Int_t i = 0; i < fgkNsec2; ++i) fDigiSet[i] = new set<Int_t> [nRows];
  76. //Get ROOT Manager
  77. FairRootManager* ioman = FairRootManager::Instance();
  78. if (ioman == 0) {
  79. Error("MpdTpcClusterFinderMlem::Init","RootManager not instantiated!");
  80. return kERROR;
  81. }
  82. // Get input collection
  83. fDigiArray = (TClonesArray*) ioman->GetObject("MpdTpcDigit");
  84. if (fDigiArray == 0) {
  85. Error("TpcClusterFinderMlem::Init","Array of digits not found!");
  86. return kERROR;
  87. }
  88. // Create and register output array
  89. //SetPersistence();
  90. fClusArray = new TClonesArray("MpdTpc2dCluster");
  91. ioman->Register("TpcCluster", "Tpc", fClusArray, fPersistence);
  92. fHitArray = new TClonesArray("MpdTpcHit");
  93. ioman->Register("TpcRecPoint", "Tpc", fHitArray, fPersistence);
  94. return kSUCCESS;
  95. }
  96. //__________________________________________________________________________
  97. void MpdTpcClusterFinderMlem::Exec(Option_t* opt)
  98. {
  99. tStartMlem = clock();
  100. fClusArray->Delete();
  101. fHitArray->Delete();
  102. //const Int_t nSec = fgkNsec2 / 2; // number of TPC readout sectors
  103. // Clear digi containers
  104. Int_t nRows = fSecGeo->NofRows();
  105. for (Int_t i = 0; i < fgkNsec2; ++i) {
  106. for (Int_t j = 0; j < nRows; ++j) fDigiSet[i][j].clear();
  107. }
  108. // Fill digi containers
  109. Int_t nDigis = fDigiArray->GetEntriesFast();
  110. cout << " Total number of digits: " << nDigis << endl;
  111. for (Int_t i = 0; i < nDigis; ++i) {
  112. MpdTpcDigit *digi = (MpdTpcDigit*) fDigiArray->UncheckedAt(i);
  113. Int_t isec = digi->GetSector();
  114. Int_t irow = digi->GetRow();
  115. fDigiSet[isec][irow].insert(i);
  116. }
  117. Int_t nSum = 0;
  118. // Loop over sectors
  119. for (Int_t isec = 0; isec < fgkNsec2; ++isec) {
  120. // Loop over padrows
  121. for (Int_t irow = 0; irow < nRows; ++irow) {
  122. if (fDigiSet[isec][irow].size() == 0) continue;
  123. //cout << " Sector, row, digits: " << isec << " " << irow << " " << fDigiSet[isec][irow].size() << endl;
  124. ProcessPadrow(isec, irow);
  125. nSum += fDigiSet[isec][irow].size();
  126. }
  127. }
  128. cout << " Control sum: " << nSum << endl;
  129. // Find hits
  130. FindHits();
  131. tFinishMlem = clock();
  132. tAllMlem = tAllMlem + (tFinishMlem - tStartMlem);
  133. }
  134. //__________________________________________________________________________
  135. void MpdTpcClusterFinderMlem::ProcessPadrow(Int_t isec, Int_t irow)
  136. {
  137. // Process one padrow of a sector
  138. memset(fFlags, 0, sizeof(fFlags[0][0]) * fgkNpads * fgkNtimes);
  139. set<Int_t>::iterator it = fDigiSet[isec][irow].begin();
  140. for ( ; it != fDigiSet[isec][irow].end(); ++it) {
  141. MpdTpcDigit *digi = (MpdTpcDigit*) fDigiArray->UncheckedAt(*it);
  142. Int_t ipad = digi->GetPad(), itime = digi->GetTimeBin();
  143. if (ipad >= fgkNpads) Fatal("ProcessPadrow", "Too few pads!!! %i %i", ipad, fgkNpads);
  144. if (itime >= fgkNtimes) Fatal("ProcessPadrow", "Too few time bins!!! %i %i", itime, fgkNtimes);
  145. //AZ-110620 fCharges[ipad][itime] = digi->GetAdc();
  146. fCharges[ipad][itime] = TMath::Min (digi->GetAdc(), Float_t(fgkOvfw)); //AZ-110620
  147. //fDigis[ipad][itime] = *it;
  148. fDigis[ipad][itime] = digi->GetOrigin(); // store trackID with max charge contribution
  149. fFlags[ipad][itime] = 1;
  150. }
  151. // Find (pre)clusters
  152. Int_t nclus0 = fClusArray->GetEntriesFast(), nclus = nclus0;
  153. for (Int_t ipad = 0; ipad < fgkNpads; ++ipad) {
  154. for (Int_t itime = 0; itime < fgkNtimes; ++itime) {
  155. if (fFlags[ipad][itime] <= 0) continue;
  156. // New cluster
  157. MpdTpc2dCluster* clus = new ((*fClusArray)[nclus++]) MpdTpc2dCluster(irow, isec);
  158. //clus->Insert(fDigis[ipad][itime], irow, ipad, itime, fCharges[ipad][itime]);
  159. clus->Insert(fDigis[ipad][itime], irow, ipad+2, itime, fCharges[ipad][itime]); // extra shift
  160. clus->SetID(fDigis[ipad][itime]); // trackID
  161. clus->SetErrY(fCharges[ipad][itime]); // ADC counts
  162. fFlags[ipad][itime] = -1;
  163. for (Int_t ip = 0; ip < 2; ++ip) {
  164. for (Int_t jt = 0; jt < 2; ++jt) {
  165. if (jt == ip) continue;
  166. Int_t ip1 = ipad + ip, it1 = itime + jt;
  167. if (ip1 >= fgkNpads) continue;
  168. if (it1 >= fgkNtimes) continue;
  169. if (fFlags[ip1][it1] <= 0) continue;
  170. NextPixel(clus, ip1, it1);
  171. }
  172. }
  173. }
  174. }
  175. //cout << " Found preclusters: " << nclus - nclus0 << endl;
  176. }
  177. //__________________________________________________________________________
  178. void MpdTpcClusterFinderMlem::NextPixel(MpdTpc2dCluster* clus, Int_t ipad, Int_t itime)
  179. {
  180. // Add next pixel to the cluster
  181. //clus->Insert(fDigis[ipad][itime], clus->Row(), ipad, itime, fCharges[ipad][itime]);
  182. clus->Insert(fDigis[ipad][itime], clus->Row(), ipad+2, itime, fCharges[ipad][itime]); // extra shift
  183. clus->SetID (TMath::Min(clus->ID(),fDigis[ipad][itime])); // min trackID
  184. clus->SetErrY (TMath::Max(Double_t(clus->GetErrY()),fCharges[ipad][itime])); // max ADC counts
  185. fFlags[ipad][itime] = -1;
  186. for (Int_t ip = -1; ip < 2; ++ip) {
  187. for (Int_t it = -1; it < 2; ++it) {
  188. if (TMath::Abs(it) == TMath::Abs(ip)) continue;
  189. Int_t ip1 = ipad + ip, it1 = itime + it;
  190. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  191. if (it1 < 0 || it1 >= fgkNtimes) continue;
  192. if (fFlags[ip1][it1] <= 0) continue;
  193. NextPixel(clus, ip1, it1);
  194. }
  195. }
  196. }
  197. //__________________________________________________________________________
  198. void MpdTpcClusterFinderMlem::FindHits()
  199. {
  200. // Reconstruct hits (find local maxima)
  201. /*
  202. fCharges[ipad][itime] = digi->GetAdc();
  203. //fDigis[ipad][itime] = *it;
  204. fDigis[ipad][itime] = digi->GetOrigin(); // store trackID with max charge contribution
  205. fFlags[ipad][itime] = 1;
  206. */
  207. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  208. Int_t nclus = fClusArray->GetEntriesFast(), ihit = 0;
  209. for (Int_t iclus = 0; iclus < nclus; ++iclus) {
  210. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  211. Int_t nDigis = clus->NDigits();
  212. //if (nDigis == 1) continue; // skip too small clusters
  213. memset(fFlags, 0, sizeof(fFlags[0][0]) * fgkNpads * fgkNtimes);
  214. for (Int_t idig = 0; idig < nDigis; ++idig) {
  215. Int_t ipad = clus->Col(idig);
  216. Int_t itime = clus->Bkt(idig);
  217. fCharges[ipad][itime] = clus->Adc(idig);
  218. fFlags[ipad][itime] = 1;
  219. fDigis[ipad][itime] = idig;
  220. }
  221. // Exclude pads which are not local maxima
  222. Int_t novfw = 0;
  223. for (Int_t idig = 0; idig < nDigis; ++idig) {
  224. Int_t ipad = clus->Col(idig);
  225. Int_t itime = clus->Bkt(idig);
  226. if (clus->Adc(idig) > fgkOvfw - 1) ++novfw;
  227. for (Int_t ip = -1; ip < 2; ++ip) {
  228. for (Int_t it = -1; it < 2; ++it) {
  229. if (TMath::Abs(it) == TMath::Abs(ip)) continue; // exclude diagonals
  230. if (it == 0 && ip == 0) continue;
  231. Int_t ip1 = ipad + ip, it1 = itime + it;
  232. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  233. if (it1 < 0 || it1 >= fgkNtimes) continue;
  234. if (fFlags[ip1][it1] == 0) continue;
  235. if (clus->Adc(idig) < fCharges[ip1][it1]) { fFlags[ipad][itime] = -1; break; }
  236. else if (clus->Adc(idig) == fCharges[ip1][it1] && fFlags[ip1][it1] > 0) { fFlags[ipad][itime] = -1; break; }
  237. }
  238. }
  239. }
  240. if (novfw) clus->SetOverflows(novfw);
  241. multimap<Double_t,Int_t> localMax;
  242. for (Int_t idig = 0; idig < nDigis; ++idig) {
  243. Int_t ipad = clus->Col(idig);
  244. Int_t itime = clus->Bkt(idig);
  245. if (fFlags[ipad][itime] <= 0) continue;
  246. localMax.insert(pair<Double_t,Int_t>(clus->Adc(idig),idig));
  247. //cout << clus->Col(idig) << " " << clus->Bkt(idig) << " " << clus->Adc(idig) << endl;
  248. }
  249. //cout << " Local max: " << clus->GetSect() << " " << clus->Row() << " " << localMax.size() << endl;
  250. Int_t nLocMax0 = localMax.size();
  251. clus->SetUniqueID(nLocMax0); // for testing
  252. // Check for edge effect - maximum digit at the sector edge
  253. Int_t edge = 0;
  254. Int_t idig = localMax.rbegin()->second;
  255. Int_t ipad = clus->Col(idig);
  256. if (ipad-2 == 0 || ipad-2 == 2*fSecGeo->NPadsInRows()[clus->Row()] - 1) edge = 1; // remember extra shift !!!
  257. //if (clus->GetNumPads() > 1) edge = -edge;
  258. if (edge > 0) {
  259. // Add "virtual" digit
  260. Int_t itime = clus->Bkt(idig), ipadv = 0, ipadn = -9;
  261. if (clus->GetNumPads() == 1) ipadv = (ipad-2 == 0) ? ipad + 1 : ipad - 1;
  262. else ipadv = (ipad-2 == 0) ? ipad - 1 : ipad + 1;
  263. fFlags[ipadv][itime] = -1;
  264. //if (clus->GetNumPads() == 1) fCharges[ipadv][itime] = 20; // below threshold
  265. //AZ-161220 if (clus->GetNumPads() == 1) fCharges[ipadv][itime] = TMath::Min(29.0,fCharges[ipad][itime]*0.1); // below threshold
  266. if (clus->GetNumPads() == 1) fCharges[ipadv][itime] = TMath::Min(2.0,fCharges[ipad][itime]*0.1); // below threshold
  267. //else fCharges[ipadv][itime] = fCharges[ipad][itime] * 0.9;
  268. else {
  269. ipadn = (ipad-2 == 0) ? ipad + 1 : ipad - 1;
  270. fCharges[ipadv][itime] = fCharges[ipadn][itime];
  271. //fCharges[ipadv][itime] = fCharges[ipad][itime] * 1.5;
  272. //fCharges[ipadv][itime] = TMath::Min (fCharges[ipadv][itime], 1.0*fgkOvfw);
  273. }
  274. fDigis[ipadv][itime] = fDigis[ipad][itime];
  275. clus->Insert(fDigis[ipadv][itime], clus->Row(), ipadv, itime, fCharges[ipadv][itime]);
  276. //cout << " Edge: " << clus->Row() << " " << clus->GetId() << " " << ipad << endl;
  277. clus->SetVirtual(); // flag "virtual" pad addition
  278. }
  279. if (edge) clus->SetEdge(); // flag edge effect
  280. //AZ if (edge || localMax.size() > 1) {
  281. //if (1) {
  282. if (edge || localMax.size() > 1 || novfw) { //AZ - 030620
  283. // Run MLEM procedure
  284. //PeakAndValley(clus, localMax);
  285. //AZ clus->SetFlag(clus->Flag()|1); // flag MLEM procedure usage
  286. if (edge || localMax.size() > 1 || novfw) clus->SetFlag(clus->Flag()|1); // flag MLEM procedure usage
  287. //cout << " xxx " << edge << " " << localMax.size() << " " << novfw << " " << clus->Flag() << endl; //AZ
  288. Mlem(iclus, localMax);
  289. continue;
  290. }
  291. map<Int_t,Double_t> mapIdQ;
  292. Double_t padMean = 0, timeMean = 0, adcTot = 0, sum2t = 0, sum2p = 0;
  293. // Process simple cluster (only 1 local max)
  294. for (Int_t jdig = 0; jdig < nDigis; ++jdig) {
  295. Int_t ip = clus->Col(jdig);
  296. Int_t it = clus->Bkt(jdig);
  297. padMean += ip * fCharges[ip][it];
  298. timeMean += it * fCharges[ip][it];
  299. adcTot += fCharges[ip][it];
  300. Int_t id = clus->Sec(fDigis[ip][it]);
  301. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = fCharges[ip][it];
  302. //else mapIdQ[id] = mapIdQ[id] + fCharges[ip][it];
  303. //else mapIdQ[id] = TMath::Max (mapIdQ[id], fCharges[ip][it]);
  304. if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = 1;
  305. else mapIdQ[id] = mapIdQ[id] + 1;
  306. sum2t += it * it * fCharges[ip][it];
  307. sum2p += ip * ip * fCharges[ip][it];
  308. }
  309. //padMean /= adcTot;
  310. padMean = padMean / adcTot - 2; // correct for shift
  311. timeMean /= adcTot;
  312. Double_t xloc = fSecGeo->Pad2Xloc(padMean,clus->Row());
  313. Int_t padID = fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  314. Double_t yloc = fSecGeo->LocalPadPosition(padID).Y();
  315. Double_t zloc = fSecGeo->TimeBin2Z(timeMean);
  316. p3loc.SetXYZ(xloc, yloc, zloc);
  317. Double_t rmsZ = TMath::Sqrt (sum2t/adcTot - timeMean*timeMean);
  318. Double_t rmsX = TMath::Sqrt (sum2p/adcTot - padMean*padMean);
  319. //p3err[1] = fSecGeo->TimeBin2Z(rms); // to pass the value
  320. //p3err[1] = rms;
  321. //cout << " Result: " << nLocMax0 << " " << pixels.size() << " " << xloc << " " << zloc << endl;
  322. // Apply corrections
  323. TVector3 p3errCor(p3err);
  324. CorrectReco(p3loc, p3errCor, clus->GetNumPads(), adcTot);
  325. if (clus->GetSect() >= fSecGeo->NofSectors()) p3loc[2] = -p3loc[2];
  326. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  327. // Dip angle correction (for interaction point with Z = 0)
  328. Double_t dip = TMath::ATan2 (TMath::Abs(p3glob[2]),p3glob.Pt()) * TMath::RadToDeg();
  329. p3errCor[2] = TMath::Max (p3errCor[2], 0.116 - 0.0046 * dip + 0.00015 * dip * dip);
  330. p3errCor[2] = TMath::Min (p3errCor[2], 0.5);
  331. // Correct Z-coordinate for rows = 0 and 27
  332. //if ((clus->Row() == 0 || clus->Row() == 27) && rmsZ > 1.0) {
  333. if (clus->Row() == 0 && rmsZ > 1.0) {
  334. Double_t zcor = 0;
  335. if (clus->Row() == 0) zcor = -0.011 + 0.002 * dip;
  336. else zcor = -0.029 + 0.005 * dip + 3.886e-5 * dip * dip;
  337. zcor = TMath::Max (zcor, 0.);
  338. zcor = TMath::Min (zcor, 0.6);
  339. p3loc[2] -= TMath::Sign (zcor, p3loc[2]);
  340. p3glob[2] -= TMath::Sign (zcor, p3glob[2]);
  341. }
  342. ihit = fHitArray->GetEntriesFast();
  343. MpdTpcHit* hit = new ((*fHitArray)[ihit]) MpdTpcHit(padID, p3glob, p3errCor, iclus);
  344. hit->SetLayer(clus->Row());
  345. hit->SetLocalPosition(p3loc); // point position
  346. hit->SetEnergyLoss(adcTot);
  347. Int_t ireg = (clus->Row() < fSecGeo->NofRowsReg(0)) ? 0 : 1;
  348. Double_t step = fSecGeo->PadHeight(ireg);
  349. hit->SetStep(step);
  350. hit->SetModular(1); // modular geometry flag
  351. hit->SetPad(Int_t(padMean));
  352. hit->SetBin(Int_t(timeMean));
  353. hit->SetRMS(rmsX, 0);
  354. hit->SetRMS(rmsZ, 1);
  355. hit->SetNdigits(nDigis);
  356. hit->SetFlags(clus);
  357. // !!! Warning: FairLinks are not persistent !!!
  358. //hit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  359. vector<Int_t> vecDig;
  360. // Store maximum 3 track IDs with the highest charge contribution
  361. multimap<Double_t,Int_t> mapDig;
  362. for (map<Int_t,Double_t>::iterator mit = mapIdQ.begin(); mit != mapIdQ.end(); ++mit)
  363. mapDig.insert(pair<Double_t,Int_t>(-mit->second,mit->first));
  364. multimap<Double_t,Int_t>::iterator mit = mapDig.begin();
  365. while (mit != mapDig.end() && vecDig.size() < 9) {
  366. //while (mit != mapDig.end() && vecDig.size() < 1) {
  367. vecDig.push_back(mit->second);
  368. ++mit;
  369. }
  370. Int_t ndig = vecDig.size();
  371. /*
  372. for (Int_t jdig = 0; jdig < ndig; ++jdig)
  373. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, clus->Sec(vecDig[jdig]))); // trackID stored in Sec
  374. */
  375. for (Int_t jdig = 0; jdig < ndig; ++jdig) {
  376. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, vecDig[jdig], jdig)); // weight = jdig
  377. hit->AddID(vecDig[jdig]);
  378. }
  379. //cout << hit->GetNLinks() << " " << *(vecDig.begin()) << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(0).GetIndex() << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(hit->GetNLinks()-1).GetIndex() << endl;
  380. } // for (Int_t iclus = 0; iclus < nclus;
  381. }
  382. //__________________________________________________________________________
  383. void MpdTpcClusterFinderMlem::PeakAndValley(const MpdTpc2dCluster* clus, multimap<Double_t,Int_t> &localMax)
  384. {
  385. // Apply peak-and-valley cuts to remove some local maxima
  386. const Double_t ratio = 1.3;
  387. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin(), rit1;
  388. ++rit;
  389. for ( ; rit != localMax.rend(); ++rit) {
  390. Int_t idig = rit->second;
  391. Int_t ipad = clus->Col(idig);
  392. Int_t itime = clus->Bkt(idig);
  393. if (fFlags[ipad][itime] <= 0) continue; // merged local max
  394. // Loop over higher peaks
  395. rit1 = localMax.rbegin();
  396. for ( ; rit1 != rit; ++rit1) {
  397. Int_t idig0 = rit1->second;
  398. Int_t ipad0 = clus->Col(idig0);
  399. Int_t itime0 = clus->Bkt(idig0);
  400. if (fFlags[ipad0][itime0] <= 0) continue; // merged local max
  401. Int_t dpad = ipad - ipad0, dt = itime - itime0;
  402. Int_t i0 = itime0, i1 = itime, j0 = ipad0, j1 = ipad, intime = 1;
  403. if (TMath::Abs(dpad) > TMath::Abs(dt)) i0 = ipad0, i1 = ipad, j0 = itime0, j1 = itime, intime = 0;
  404. Int_t stepi = TMath::Sign(1,i1-i0), stepj = TMath::Sign(1,j1-j0);
  405. //Int_t valOk = 0;
  406. Int_t merge = 0;
  407. //if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) merge = 1;
  408. //else {
  409. if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) {
  410. if (fCharges[ipad-dpad][itime] > fCharges[ipad][itime] / ratio ||
  411. fCharges[ipad][itime-dt] > fCharges[ipad][itime] / ratio) merge = 1;
  412. } else {
  413. for (Int_t ii = i0 + stepi; ii != i1; ii += stepi) {
  414. merge = 0;
  415. for (Int_t jj = j0; jj != j1 + stepj; jj += stepj) {
  416. if (TMath::Abs(jj-j0) > TMath::Abs(ii-i0)) continue;
  417. if (TMath::Abs(jj-j1) > TMath::Abs(ii-i1)) continue;
  418. Int_t ip = ii, it = jj;
  419. if (intime) ip = jj, it = ii;
  420. //if (fCharges[ip][it] < fCharges[ipad][itime] / ratio) { valOk = 1; break; }
  421. if (fCharges[ip][it] > fCharges[ipad][itime] / ratio) { merge = 1; break; }
  422. }
  423. //if (valOk) break;
  424. if (!merge) break;
  425. }
  426. }
  427. //if (!valOk) fFlags[ipad][itime] = -fFlags[ipad][itime]; // not deep enough valley
  428. if (merge) fFlags[ipad][itime] = -fFlags[ipad][itime]; // not deep enough valley
  429. }
  430. }
  431. // Remove failed peaks
  432. multimap<Double_t,Int_t>::iterator it = localMax.begin();
  433. //for ( ; it != localMax.end(); ++it) cout << it->second << " ";
  434. //cout << endl;
  435. it = localMax.begin();
  436. for ( ; it != localMax.end(); ++it) {
  437. Int_t idig = it->second;
  438. Int_t ipad = clus->Col(idig);
  439. Int_t itime = clus->Bkt(idig);
  440. if (fFlags[ipad][itime] > 0) continue;
  441. //cout << " Before: " << idig << " " << itime << " " << ipad << " " << it->first << ", ";
  442. localMax.erase(it);
  443. //cout << " After: " << it->first << endl;
  444. }
  445. }
  446. //__________________________________________________________________________
  447. void MpdTpcClusterFinderMlem::CorrectReco(TVector3 &p3loc, TVector3 &p3errCor, Int_t nPads, Double_t adc)
  448. {
  449. // Correct reconstructed coordinates
  450. Double_t xsec = (p3loc.Y() + fSecGeo->GetMinY()) * TMath::Tan(fSecGeo->Dphi()/2);
  451. Double_t xloc = p3loc.X(), xloc0 = xloc, edge = 0.0; // distance to sector boundary
  452. if (xloc < 0) edge = xloc + xsec;
  453. else edge = xloc - xsec;
  454. if (TMath::Abs(edge) > 1.5 || nPads > 3) return; // no corrections
  455. if (nPads == 1) {
  456. /*
  457. xloc -= TMath::Sign(0.520,edge);
  458. Double_t corr = -0.15;
  459. if (adc < 50000) corr = 0.206 - 2.532e-5 * adc + 8.715e-10 * adc * adc
  460. - 1.610e-14 * adc * adc * adc + 1.193e-19 * adc * adc * adc * adc;
  461. p3errCor[0] = 0.116;
  462. */
  463. Double_t corr = 0.092;
  464. if (edge < 0) corr = -corr;
  465. xloc -= corr;
  466. p3errCor[0] = 0.386;
  467. } else if (nPads == 2) {
  468. /*
  469. xloc -= TMath::Sign(0.200,edge);
  470. Double_t corr = -0.07;
  471. if (adc < 200000) corr = 0.374 - 1.046e-5 * adc + 9.595e-11 * adc * adc
  472. - 4.098e-16 * adc * adc * adc + 6.811e-22 * adc * adc * adc * adc;
  473. p3errCor[0] = 0.110;
  474. */
  475. Double_t corr = 0.002;
  476. if (edge < 0) corr = -corr;
  477. xloc -= corr;
  478. p3errCor[0] = 0.03;
  479. } else if (nPads > 2) {
  480. /*
  481. xloc -= TMath::Sign(0.033,edge);
  482. Double_t corr = 0.01;
  483. if (adc < 200000) corr = 0.188 - 4.007e-6 * adc + 3.364e-11 * adc * adc
  484. - 1.265e-16 * adc * adc * adc + 1.771e-22 * adc * adc * adc * adc;
  485. p3errCor[0] = 0.073;
  486. */
  487. Double_t corr = 0.004;
  488. if (edge < 0) corr = -corr;
  489. xloc -= corr;
  490. p3errCor[0] = 0.04;
  491. } else return;
  492. if (TMath::Abs(xloc) > xsec) xloc = TMath::Sign(xsec,xloc0);
  493. p3loc.SetX(xloc);
  494. }
  495. //__________________________________________________________________________
  496. void MpdTpcClusterFinderMlem::Mlem(Int_t iclus, multimap<Double_t,Int_t> &localMax)
  497. {
  498. // Run MLEM procedure
  499. const Int_t nDigMax = 200; // max number of digits to apply pixel reduction
  500. if (gROOT->FindObject("hMlem0")) gROOT->FindObject("hMlem0")->Delete();
  501. if (gROOT->FindObject("hMlem1")) gROOT->FindObject("hMlem1")->Delete();
  502. if (gROOT->FindObject("hTimePad")) gROOT->FindObject("hTimePad")->Delete();
  503. //TCanvas *c = (TCanvas*) gROOT->FindObject("c");
  504. vector<pixel> bins, pixels[2];
  505. vector<pixel>::iterator it;
  506. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  507. Int_t nDigis = clus->NDigits();
  508. Int_t nx = clus->GetNumTimeBins();
  509. Int_t ny = clus->GetNumPads() + 2;
  510. TH2D *hXY = new TH2D("hTimePad","Pad No. vs Time bin",nx,clus->MinBkt(),clus->MaxBkt()+1,
  511. ny,clus->MinCol()-1,clus->MaxCol()+2);
  512. TH2D *hOvfw = nullptr;
  513. if (clus->Overflows() > 1 && nDigis < 50) {
  514. // Special treatment for (relatively) small clusters with many overflows
  515. if (gROOT->FindObject("hOvfw")) gROOT->FindObject("hOvfw")->Delete();
  516. hOvfw = (TH2D*) hXY->Clone("hOvfw");
  517. }
  518. // Check for edge effect - maximum digit at the edge
  519. Int_t edge = 0;
  520. Int_t idig = localMax.rbegin()->second;
  521. Int_t ipad = clus->Col(idig);
  522. if (ipad-2 == 0 || ipad-2 == 2*fSecGeo->NPadsInRows()[clus->Row()] - 1) edge = 1; // extra shift
  523. for (Int_t jdig = 0; jdig < nDigis; ++jdig) {
  524. ipad = clus->Col(jdig);
  525. Int_t itime = clus->Bkt(jdig);
  526. //fCharges[ipad][itime] = clus->Adc(jdig);
  527. //fFlags[ipad][itime] = 1;
  528. //fDigis[ipad][itime] = jdig;
  529. //cout << " digis: " << jdig << " " << ipad << " " << itime << " " << fCharges[ipad][itime] << endl;
  530. if (fCharges[ipad][itime] < 0.1) continue;
  531. pixel pix;
  532. pix.ix = hXY->GetXaxis()->FindBin(itime+0.1);
  533. pix.iy = hXY->GetYaxis()->FindBin(ipad+0.1);
  534. if (hXY->GetBinContent(pix.ix,pix.iy) > 0.1) continue; //AZ-161220
  535. pix.qq = fCharges[ipad][itime];
  536. pix.vis = 0;
  537. //if (pix.qq < fgkOvfw) bins.push_back(pix);
  538. bins.push_back(pix);
  539. pixels[1].push_back(pix);
  540. pixels[0].push_back(pix);
  541. hXY->SetBinContent(pix.ix,pix.iy,pix.qq); //AZ-161220
  542. if (hOvfw && fCharges[ipad][itime] > fgkOvfw - 1) hOvfw->Fill(itime+0.1,ipad+0.1,fCharges[ipad][itime]);
  543. // Consider edge effect
  544. //*
  545. if (edge && (ipad-2 == 0 || ipad-2 == 2*fSecGeo->NPadsInRows()[clus->Row()] - 1)) { // extra shift
  546. //cout << " Edge: " << ipad << " " << pix.iy << " " << pix.ix << " " << hXY->GetYaxis()->GetBinCenter(pix.iy) << endl;
  547. if (ipad-2 == 0) {
  548. --pix.iy;
  549. fDigis[ipad-1][itime] = fDigis[ipad][itime];
  550. } else {
  551. ++pix.iy;
  552. fDigis[ipad+1][itime] = fDigis[ipad][itime];
  553. }
  554. if (hXY->GetBinContent(pix.ix,pix.iy) > 0.1) continue;
  555. pixels[1].push_back(pix);
  556. pixels[0].push_back(pix);
  557. hXY->SetBinContent(pix.ix,pix.iy,pix.qq); //AZ-161220
  558. //cout << " Edge: " << clus->Row() << " " << clus->GetId() << " " << ipad << endl;
  559. }
  560. //*/
  561. }
  562. // Get parameters of the response function
  563. Double_t sigt, sigp, correl;
  564. GetResponse(clus, hXY, hOvfw, sigt, sigp, correl);
  565. Int_t nbins = bins.size();
  566. //sort (pixels[1].begin(), pixels[1].end(), SortPix);
  567. //if (pixels[1].size() > nbins) pixels[1].erase(pixels[1].begin()+nbins,pixels[1].end());
  568. //for (Int_t ibin = 0; ibin < nbins; ++ibin) cout << ibin << " " << pixels[1][ibin].qq << " "
  569. // << pixels[1][ibin].ix << " " << pixels[1][ibin].iy << endl;
  570. //TH2D *hMlem[2] = {nullptr, (TH2D*)hXY->Clone("hMlem1")};
  571. TH2D *hMlem[2] = {(TH2D*)hXY->Clone("hMlem0"), (TH2D*)hXY->Clone("hMlem1")};
  572. //TH2D *hExp = (TH2D*) hXY->Clone("hExp"); // for testing
  573. // Decrease pixel size
  574. //const Int_t ndiv = 5, ndiv2 = ndiv * ndiv, nloops = 2; //
  575. const Int_t ndiv = 3, ndiv2 = ndiv * ndiv, nloops = 2; //
  576. Int_t icur = 0;
  577. TMatrixD cij(2,2);
  578. vector<pair<Int_t,Int_t> > pairs;
  579. const Double_t cijMin = 1.e-5;
  580. //const Double_t cijMin = 0.001;
  581. for (Int_t iloop = 0; iloop < nloops; ++iloop) {
  582. if (iloop) {
  583. if (nDigis > nDigMax) continue;
  584. nx *= ndiv;
  585. ny *= ndiv;
  586. } else if (nDigis <= nDigMax) continue; // skip first iteration
  587. Int_t iprev = (iloop + 1) % 2;
  588. icur = iloop % 2;
  589. TString hname = "hMlem";
  590. hname += icur;
  591. if (hMlem[icur]) hMlem[icur]->Delete();
  592. hMlem[icur] = new TH2D (hname,"",nx,hXY->GetXaxis()->GetXmin(),hXY->GetXaxis()->GetXmax(),
  593. ny,hXY->GetYaxis()->GetXmin(),hXY->GetYaxis()->GetXmax());
  594. pixels[icur].clear();
  595. Int_t npix = 0;
  596. if (iloop) cij.ResizeTo(pixels[iprev].size()*ndiv2,nbins);
  597. else cij.ResizeTo(pixels[iprev].size(),nbins);
  598. Double_t visMax = 0.0;
  599. //sort (pixels[iprev].begin(), pixels[iprev].end(), SortPix);
  600. for (it = pixels[iprev].begin(); it != pixels[iprev].end(); ++it) {
  601. Double_t xc = hMlem[iprev]->GetXaxis()->GetBinCenter(it->ix);
  602. Double_t yc = hMlem[iprev]->GetYaxis()->GetBinCenter(it->iy);
  603. Double_t xmin = hMlem[iprev]->GetXaxis()->GetBinLowEdge(it->ix);
  604. Double_t ymin = hMlem[iprev]->GetYaxis()->GetBinLowEdge(it->iy);
  605. for (Int_t ij = 0; ij < ndiv2; ++ij) {
  606. if (iloop == 0 && ij) break;
  607. Double_t x0 = xc;
  608. Double_t y0 = yc;
  609. pixel pix;
  610. if (iloop == 0) {
  611. pix.ix = it->ix;
  612. pix.iy = it->iy;
  613. } else {
  614. x0 = xmin + hMlem[icur]->GetXaxis()->GetBinWidth(1) * (ij % ndiv + 0.5);
  615. y0 = ymin + hMlem[icur]->GetYaxis()->GetBinWidth(1) * (ij / ndiv + 0.5);
  616. pix.ix = hMlem[icur]->GetXaxis()->FindBin(x0);
  617. pix.iy = hMlem[icur]->GetYaxis()->FindBin(y0);
  618. //cout << xc << " " << yc << " " << x0 << " " << y0 << " " << pix.ix << " " << pix.iy << endl;
  619. }
  620. pix.qq = 1.0;
  621. //if (npix == 0) pix.qq = 9746;
  622. pix.vis = 0.0;
  623. pixels[icur].push_back(pix);
  624. // Compute couplings for this pixel
  625. Double_t vis = 0, invis = 0;
  626. for (Int_t ibin = 0; ibin < nbins; ++ibin) {
  627. Double_t x1 = hXY->GetXaxis()->GetBinCenter(bins[ibin].ix);
  628. Double_t y1 = hXY->GetYaxis()->GetBinCenter(bins[ibin].iy);
  629. cij(npix,ibin) = GetCij(clus->Row(), x0, y0, x1, y1, sigt, sigp, correl, vis);
  630. //?? pixels[icur][npix].vis += cij(npix,ibin); // visibility
  631. if (cij(npix,ibin) > cijMin) {
  632. //if (bins[ibin].qq < fgkOvfw - 1.5) pixels[icur][npix].vis += cij(npix,ibin); // visibility
  633. //else invis += cij(npix,ibin); // "invisibility"
  634. if (bins[ibin].qq > fgkOvfw - 1.5) invis += cij(npix,ibin); // invisibility
  635. pairs.push_back(pair<Int_t,Int_t>(npix,ibin));
  636. }
  637. }
  638. pixels[icur][npix].vis = vis;
  639. visMax = TMath::Max(visMax,pixels[icur][npix].vis);
  640. //AZ 04.07.20 pixels[icur][npix].vis = pixels[icur][npix].vis / (pixels[icur][npix].vis + invis); //AZ - 030620
  641. pixels[icur][npix].vis -= invis; //AZ - 050720
  642. //pixels[icur][npix].vis /= visMax; //AZ - 050720
  643. pixels[icur][npix].vis /= 1.2; //AZ - 050720
  644. //cout << npix << " " << pixels[icur][npix].vis << " " << invis << endl; //AZ
  645. ++npix;
  646. //if (npix >= nbins) break;
  647. } // for (Int_t ij = 0;
  648. //if (npix >= nbins) break;
  649. } // for (it = pixels[iprev].begin();
  650. //for (Int_t ipix = 0; ipix < npix; ++ipix) cout << ipix << " " << pixels[icur][ipix].vis << " " << pixels[icur][ipix].qq << " " << pixels[icur][ipix].ix << " " << pixels[icur][ipix].iy << endl;
  651. // MLEM procedure
  652. Int_t niter = 5; //10;
  653. vector<pair<Int_t,Int_t> >::iterator vit;
  654. for (Int_t iter = 0; iter < niter; ++iter) {
  655. // Calculate expectations for all pads
  656. for (Int_t ibin = 0; ibin < nbins; ++ibin) {
  657. bins[ibin].vis = 0.0;
  658. //AZ for (Int_t i = 0; i < npix; ++i) bins[ibin].vis += (pixels[icur][i].qq * cij(i,ibin)); // bins[ibin].vis - just storage
  659. //?? bins[ibin].vis = TMath::Min (bins[ibin].vis, fglOvfw);
  660. }
  661. for (vit = pairs.begin(); vit != pairs.end(); ++vit) {
  662. Int_t i = vit->first, ibin = vit->second;
  663. bins[ibin].vis += (pixels[icur][i].qq * cij(i,ibin)); // bins[ibin].vis - just storage
  664. }
  665. // Maximization
  666. /*
  667. for (Int_t ipix = 0; ipix < npix; ++ipix) {
  668. Double_t sum = 0.0;
  669. for (Int_t ibin = 0; ibin < nbins; ++ibin) {
  670. if (bins[ibin].vis > 1.e-6) sum += (bins[ibin].qq / bins[ibin].vis * cij(ipix,ibin));
  671. }
  672. if (pixels[icur][ipix].vis > 1.e-6) pixels[icur][ipix].qq *= (sum / pixels[icur][ipix].vis);
  673. //if (pixels[icur][ipix].vis > 1.e-6) pixels[icur][ipix].qq *= (sum / visMax);
  674. }
  675. */
  676. for (Int_t ipix = 0; ipix < npix; ++ipix) pixels[icur][ipix].sum = 0.0;
  677. for (vit = pairs.begin(); vit != pairs.end(); ++vit) {
  678. Int_t ipix = vit->first, ibin = vit->second;
  679. if (bins[ibin].vis > 1.e-6) pixels[icur][ipix].sum += (bins[ibin].qq / bins[ibin].vis * cij(ipix,ibin));
  680. //if (bins[ibin].vis > 1.e-6 && bins[ibin].qq < fgkOvfw - 1.5) pixels[icur][ipix].sum += (bins[ibin].qq / bins[ibin].vis * cij(ipix,ibin)); //AZ - 030620
  681. }
  682. for (Int_t ipix = 0; ipix < npix; ++ipix)
  683. if (pixels[icur][ipix].vis > 1.e-6) pixels[icur][ipix].qq *= (pixels[icur][ipix].sum / pixels[icur][ipix].vis);
  684. //if (pixels[icur][ipix].vis > 1.e-6) pixels[icur][ipix].qq *= (pixels[icur][ipix].sum / 1); //AZ - 030620
  685. } // for (Int_t iter = 0;
  686. // Visual control
  687. /*
  688. if (iloop == nloops - 1) {
  689. for (Int_t ipix = 0; ipix < npix; ++ipix)
  690. hMlem[icur]->SetBinContent(pixels[icur][ipix].ix,pixels[icur][ipix].iy,pixels[icur][ipix].qq);
  691. //c->cd(3);
  692. hMlem[icur]->Draw("col");
  693. }
  694. */
  695. /*
  696. hExp->Reset();
  697. for (Int_t ipix = 0; ipix < nbins; ++ipix)
  698. hExp->SetBinContent(bins[ipix].ix,bins[ipix].iy,bins[ipix].vis);
  699. */
  700. } // for (Int_t iloop = 0;
  701. //
  702. // Find local maxima
  703. //
  704. vector<Double_t> tmp;
  705. vector<vector<Double_t> > charges;
  706. tmp.assign(nx,0);
  707. charges.assign(ny,tmp);
  708. vector<Int_t> itmp;
  709. vector<vector<Int_t> > flags;
  710. itmp.assign(nx,0);
  711. flags.assign(ny,itmp);
  712. localMax.clear();
  713. Int_t npix = pixels[icur].size();
  714. /*/test
  715. Double_t qmx = -1;
  716. int ipmx, itmx, iii;
  717. for (Int_t ipix = 0; ipix < npix; ++ipix) {
  718. //cout << " pix: " << ipix << " " << pixels[icur][ipix].ix - 1 << " " << pixels[icur][ipix].iy - 1 << endl;
  719. if (pixels[icur][ipix].qq > qmx) {
  720. qmx = pixels[icur][ipix].qq;
  721. ipmx = pixels[icur][ipix].iy - 1;
  722. itmx = pixels[icur][ipix].ix - 1;
  723. iii = ipix;
  724. }
  725. }
  726. //cout << " aaaa " << qmx << " " << ipmx << " " << itmx << " " << iii << endl;
  727. //test */
  728. for (Int_t ipix = 0; ipix < npix; ++ipix) {
  729. ipad = pixels[icur][ipix].iy - 1;
  730. Int_t itime = pixels[icur][ipix].ix - 1;
  731. charges[ipad][itime] = pixels[icur][ipix].qq;
  732. flags[ipad][itime] = ipix + 1;
  733. }
  734. // Exclude pads which are not local maxima
  735. for (Int_t ipix = 0; ipix < npix; ++ipix) {
  736. ipad = pixels[icur][ipix].iy - 1;
  737. Int_t itime = pixels[icur][ipix].ix - 1;
  738. Int_t iok = 1;
  739. for (Int_t ip = -1; ip < 2; ++ip) {
  740. for (Int_t jt = -1; jt < 2; ++jt) {
  741. if (TMath::Abs(jt) == TMath::Abs(ip)) continue; // exclude diagonals
  742. if (jt == 0 && ip == 0) continue;
  743. Int_t ip1 = ipad + ip, it1 = itime + jt;
  744. if (ip1 < 0 || ip1 >= ny) continue;
  745. if (it1 < 0 || it1 >= nx) continue;
  746. if (flags[ip1][it1] == 0) continue;
  747. //if (pixels[icur][ipix].qq < charges[ip1][it1]) { flags[ipad][itime] = -1; break; }
  748. //else if (pixels[icur][ipix].qq == charges[ip1][it1] && flags[ip1][it1] > 0) { flags[ipad][itime] = -1; break; }
  749. //if (flags[ip1][it1] <= 0) continue;
  750. if (pixels[icur][ipix].qq < charges[ip1][it1]) { flags[ipad][itime] *= -1; iok = 0; break; }
  751. else if (pixels[icur][ipix].qq == charges[ip1][it1] && flags[ip1][it1] > 0) { flags[ipad][itime] *= -1; iok = 0; break; }
  752. }
  753. if (!iok) break;
  754. }
  755. }
  756. for (Int_t ipix = 0; ipix < npix; ++ipix) {
  757. ipad = pixels[icur][ipix].iy - 1;
  758. Int_t itime = pixels[icur][ipix].ix - 1;
  759. if (flags[ipad][itime] <= 0) continue;
  760. localMax.insert(pair<Double_t,Int_t>(pixels[icur][ipix].qq,ipix));
  761. }
  762. // Apply Peak-and-valley procedure in pixel domain
  763. if (localMax.size() > 1) PeakAndValley(pixels[icur], localMax, charges, flags);
  764. // Create hits
  765. Int_t nHits0 = fHitArray->GetEntriesFast();
  766. vector<multimap<Double_t,Int_t> > pixInMax;
  767. multimap<Double_t,Int_t> mtmp;
  768. pixInMax.assign(localMax.size(),mtmp);
  769. CreateHits(pixels[icur], localMax, charges, flags, iclus, pixInMax);
  770. // Get correct charge
  771. ChargeMlem(nHits0, pixels[icur], bins, pixInMax, cij, cijMin);
  772. }
  773. //__________________________________________________________________________
  774. void MpdTpcClusterFinderMlem::PeakAndValley(const vector<pixel> &pixels, multimap<Double_t,Int_t> &localMax,
  775. vector<vector<Double_t> > &charges, vector<vector<Int_t> > &flags)
  776. {
  777. // Apply peak-and-valley cuts to remove some local maxima in pixel domain
  778. const Double_t ratio = 1.2;
  779. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin(), rit1;
  780. ++rit;
  781. for ( ; rit != localMax.rend(); ++rit) {
  782. Int_t ipix = rit->second;
  783. Int_t ipad = pixels[ipix].iy - 1;
  784. Int_t itime = pixels[ipix].ix - 1;
  785. if (flags[ipad][itime] <= 0) continue; // merged local max
  786. // Loop over higher peaks
  787. rit1 = localMax.rbegin();
  788. for ( ; rit1 != rit; ++rit1) {
  789. Int_t ipix0 = rit1->second;
  790. Int_t ipad0 = pixels[ipix0].iy - 1;
  791. Int_t itime0 = pixels[ipix0].ix - 1;
  792. if (flags[ipad0][itime0] <= 0) continue; // merged local max
  793. Int_t dpad = ipad - ipad0, dt = itime - itime0;
  794. Int_t i0 = itime0, i1 = itime, j0 = ipad0, j1 = ipad, intime = 1;
  795. if (TMath::Abs(dpad) > TMath::Abs(dt)) i0 = ipad0, i1 = ipad, j0 = itime0, j1 = itime, intime = 0;
  796. Int_t stepi = TMath::Sign(1,i1-i0), stepj = TMath::Sign(1,j1-j0);
  797. //Int_t valOk = 0;
  798. Int_t merge = 0;
  799. //if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) merge = 1;
  800. //else {
  801. if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) {
  802. if (charges[ipad-dpad][itime] > charges[ipad][itime] / ratio ||
  803. charges[ipad][itime-dt] > charges[ipad][itime] / ratio) merge = 1;
  804. } else {
  805. for (Int_t ii = i0 + stepi; ii != i1; ii += stepi) {
  806. merge = 0;
  807. for (Int_t jj = j0; jj != j1 + stepj; jj += stepj) {
  808. if (TMath::Abs(jj-j0) > TMath::Abs(ii-i0)) continue;
  809. if (TMath::Abs(jj-j1) > TMath::Abs(ii-i1)) continue;
  810. Int_t ip = ii, it = jj;
  811. if (intime) ip = jj, it = ii;
  812. //if (fCharges[ip][it] < fCharges[ipad][itime] / ratio) { valOk = 1; break; }
  813. if (charges[ip][it] > charges[ipad][itime] / ratio) { merge = 1; break; }
  814. }
  815. //if (valOk) break;
  816. if (!merge) break;
  817. }
  818. }
  819. //if (!valOk) fFlags[ipad][itime] = -fFlags[ipad][itime]; // not deep enough valley
  820. if (merge) flags[ipad][itime] = -flags[ipad][itime]; // not deep enough valley
  821. }
  822. }
  823. // Remove failed peaks
  824. /*
  825. multimap<Double_t,Int_t>::iterator itt = localMax.begin();
  826. //for ( ; it != localMax.end(); ++it) cout << it->second << " ";
  827. //cout << endl;
  828. itt = localMax.begin();
  829. for ( ; itt != localMax.end(); ++itt) {
  830. Int_t ipix = itt->second;
  831. Int_t ipad = pixels[ipix].iy - 1;
  832. Int_t itime = pixels[ipix].ix - 1;
  833. cout << " Before: " << ipix << " " << itime << " " << ipad << " " << itt->first << " " << flags[ipad][itime] << endl;
  834. if (flags[ipad][itime] > 0) continue;
  835. //cout << " Before: " << idig << " " << itime << " " << ipad << " " << it->first << ", ";
  836. //localMax.erase(it1);
  837. //cout << " After: " << it->first << endl;
  838. }
  839. */
  840. // Remove failed peaks
  841. multimap<Double_t,Int_t>::iterator it = localMax.end(), it1;
  842. --it;
  843. Int_t iover = 0;
  844. while (iover == 0) {
  845. Int_t ipix = it->second;
  846. Int_t ipad = pixels[ipix].iy - 1;
  847. Int_t itime = pixels[ipix].ix - 1;
  848. if (it == localMax.begin()) iover = 1;
  849. it1 = it;
  850. --it;
  851. if (flags[ipad][itime] > 0) continue;
  852. localMax.erase(it1);
  853. }
  854. }
  855. //__________________________________________________________________________
  856. /*
  857. //static Bool_t MpdTpcClusterFinderMlem::SortPix(const pixel i, const pixel j)
  858. static Bool_t SortPix(const MpdTpcClusterFinderMlem::pixel i, const MpdTpcClusterFinderMlem::pixel j)
  859. {
  860. // Sorting function
  861. return i.qq > j.qq;
  862. }
  863. */
  864. //__________________________________________________________________________
  865. void MpdTpcClusterFinderMlem::GetResponse(const MpdTpc2dCluster* clus, TH2D *hXY, TH2D *hOvfw,
  866. Double_t &sigt, Double_t &sigp, Double_t &correl)
  867. {
  868. // Get response parameters
  869. const Double_t sigT[2][4] = {{0.853542, 0.00104733, -4.49635e-07, 0.0},
  870. {0.853845, 0.00106644, -4.84877e-07, 0.0}}; // Sigma_time vs time (in time bins)
  871. const Double_t sigP[2][4] = {{0.357956, 0.00128714, -2.07788e-06, 2.15803e-09},
  872. {0.501219, -4.32769e-05, 2.54277e-06, -3.65539e-09}}; // Sigma_pad vs time
  873. // Sigma_p vs tan(cross) (1 + A*abs(x) + B*x*x + C*x*x*x*x; A,B,C - functions of time bin)
  874. const Double_t tanCrosA[2][3] = {{0.214453, -0.00196814, 4.47426e-06},
  875. {-0.869907, 0.00733796, -1.50085e-05}};
  876. const Double_t tanCrosB[2][3] = {{1.54845, -0.00341245, 4.81089e-06},
  877. {3.12359, -0.00978931, 1.54187e-05}};
  878. const Double_t tanCrosC[2][3] = {{-0.648698, 0.00258523, -4.0762e-06},
  879. {-1.05394, 0.00179292, 1.8612e-07}};
  880. // Sigma_t vs tan(dip) (1 + A*abs(x) + B*x*x + C*x*x*x*x; A,B,C - functions of time bin)
  881. const Double_t tanDipA[2][3] = {{0.2576, -0.00243985, 5.49784e-06},
  882. {0.494572, -0.00417108, 9.28609e-06}};
  883. const Double_t tanDipB[2][3] = {{0.12751, 0.00110365, -2.57362e-06},
  884. {0.129227, 0.00267634, -6.82141e-06}};
  885. const Double_t tanDipC[2][3] = {{0.00168453, -0.000124739, 2.09058e-07},
  886. {0.0188111, -0.000514421, 1.30111e-06}};
  887. // Dip and crossing angle estimates assuming straight tracks from Z=0 for now
  888. Double_t padMean = (hXY->GetYaxis()->GetXmin() + hXY->GetYaxis()->GetXmax()) / 2;
  889. Double_t timeMean = (hXY->GetXaxis()->GetXmin() + hXY->GetXaxis()->GetXmax()) / 2;
  890. Double_t t = timeMean;
  891. Double_t xloc = fSecGeo->Pad2Xloc(padMean,clus->Row());
  892. Int_t padID = fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  893. Double_t yloc = fSecGeo->LocalPadPosition(padID).Y();
  894. Double_t zloc = fSecGeo->TimeBin2Z(timeMean);
  895. TVector3 p3loc(xloc, yloc, zloc), p3glob;
  896. if (clus->GetSect() >= fSecGeo->NofSectors()) p3loc[2] = -p3loc[2];
  897. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  898. Double_t tanDip = TMath::Abs(p3glob[2]) / p3glob.Pt();
  899. tanDip = TMath::Min (tanDip, 2.0);
  900. Double_t tanDip2 = tanDip * tanDip;
  901. Double_t tanDip4 = tanDip2 * tanDip2;
  902. //Double_t tanCros = TMath::Abs(xloc) / p3glob.Pt();
  903. //tanCros = TMath::Min (tanCros, 1.0);
  904. //Double_t tanCros2 = tanCros * tanCros;
  905. //Double_t tanCros4 = tanCros2 * tanCros2;
  906. Int_t ireg = (clus->Row() < 27) ? 0 : 1; // pad height region
  907. Double_t t2 = t * t;
  908. sigt = sigT[ireg][0] + sigT[ireg][1]*t + sigT[ireg][2]*t2;
  909. Double_t a = tanDipA[ireg][0] + tanDipA[ireg][1]*t + tanDipA[ireg][2]*t2;
  910. Double_t b = tanDipB[ireg][0] + tanDipB[ireg][1]*t + tanDipB[ireg][2]*t2;
  911. Double_t c = tanDipC[ireg][0] + tanDipC[ireg][1]*t + tanDipC[ireg][2]*t2;
  912. sigt *= (1 + a*tanDip + b*tanDip2 + c*tanDip4);
  913. sigp = sigP[ireg][0] + sigP[ireg][1]*t + sigP[ireg][2]*t2 + sigP[ireg][3]*t2*t;
  914. a = tanCrosA[ireg][0] + tanCrosA[ireg][1]*t + tanCrosA[ireg][2]*t2;
  915. b = tanCrosB[ireg][0] + tanCrosB[ireg][1]*t + tanCrosB[ireg][2]*t2;
  916. c = tanCrosC[ireg][0] + tanCrosC[ireg][1]*t + tanCrosC[ireg][2]*t2;
  917. //AZ sigp *= (1 + a*tanCros + b*tanCros2 + c*tanCros4);
  918. sigp *= (1 + 0);
  919. correl = 0.0;
  920. if (hOvfw) {
  921. correl = hOvfw->GetCorrelationFactor();
  922. //sigt = hOvfw->GetRMS(1);
  923. //sigp = hOvfw->GetRMS(2);
  924. }
  925. }
  926. /*
  927. //__________________________________________________________________________
  928. Double_t MpdTpcClusterFinderMlem::GetCij(Double_t x0, Double_t y0, Double_t x1, Double_t y1,
  929. Double_t sigt, Double_t sigp, Double_t correl)
  930. {
  931. // Compute couplings - Gauss parameter Sigma_t as a function of the dip angle
  932. Double_t dx = x1 - x0;
  933. Double_t dy = y1 - y0;
  934. dx /= sigt;
  935. dy /= sigp;
  936. Double_t cij = 0;
  937. if (TMath::Abs(correl) < 0.001) cij = TMath::Exp(-dx*dx/2 - dy*dy/2) / sigt / sigp / TMath::TwoPi();
  938. else {
  939. Double_t corr2 = 1 - correl * correl;
  940. cij = TMath::Exp((-dx*dx/2 + correl*dx*dy - dy*dy/2)/corr2) / sigt / sigp / TMath::TwoPi() / TMath::Sqrt(corr2);
  941. }
  942. cout << x1-x0 << " " << y1-y0 << " " << cij << endl;
  943. return cij;
  944. }
  945. */
  946. //*
  947. //__________________________________________________________________________
  948. Double_t MpdTpcClusterFinderMlem::GetCij(Int_t irow, Double_t x0, Double_t y0, Double_t x1, Double_t y1,
  949. Double_t sigt, Double_t sigp, Double_t correl, Double_t &vis)
  950. {
  951. // Compute couplings - Gauss parameter Sigma_t as a function of the dip angle.
  952. // Make use of cache (map) to process the same pixel.
  953. static Int_t ixpix0 = -999, iypix0 = -999;
  954. static map<Int_t,map<Int_t,Double_t> > cijCache;
  955. Int_t ixpix = x0 * 1000, iypix = y0 * 1000, ixc = 0, iyc = 0;
  956. Int_t ixbin = Int_t(x1), iybin = Int_t(y1);
  957. if (ixpix != ixpix0 || iypix != iypix0) {
  958. cijCache.clear();
  959. ixpix0 = ixpix;
  960. iypix0 = iypix;
  961. ixc = Int_t(x0);
  962. iyc = Int_t(y0);
  963. Int_t nxy = 2;
  964. vis = 0.0;
  965. // Compute Cij for the bin area around central pixel
  966. for (Int_t i = -nxy; i <= nxy; ++i) {
  967. Int_t it = ixc + i;
  968. Double_t dx = it - x0 + 0.5;
  969. for (Int_t j = -nxy; j <= nxy; ++j) {
  970. Int_t ip = iyc + j;
  971. //if (ip-2 == 0 || ip-2 == 2*fSecGeo->NPadsInRows()[irow] - 1) continue; // edge
  972. if (ip-2 < 0 || ip-2 > 2*fSecGeo->NPadsInRows()[irow] - 1) continue; // edge
  973. Double_t dy = ip - y0 + 0.5;
  974. dx /= sigt;
  975. dy /= sigp;
  976. map<Int_t,Double_t> aaa;
  977. if (cijCache.count(it) == 0) cijCache[it] = aaa;
  978. if (TMath::Abs(correl) < 0.001) cijCache[it][ip] = TMath::Exp(-dx*dx/2 - dy*dy/2) / sigt / sigp / TMath::TwoPi();
  979. else {
  980. Double_t corr2 = 1 - correl * correl;
  981. cijCache[it][ip] =
  982. TMath::Exp((-dx*dx/2 + correl*dx*dy - dy*dy/2)/corr2) / sigt / sigp / TMath::TwoPi() / TMath::Sqrt(corr2);
  983. }
  984. vis += cijCache[it][ip];
  985. //cout << it - x0 + 0.5 << " " << ip - y0 + 0.5 << " " << cijCache[it][ip] << " " << it << " " << ip << endl;
  986. }
  987. }
  988. //cout << " Visibility: " << vis << endl;
  989. } else if (cijCache.count(ixbin) == 0 || cijCache[ixbin].count(iybin) == 0) {
  990. //cout << ixbin << " " << iybin << endl;
  991. return 0.0; // outside computed area
  992. }
  993. return cijCache[ixbin][iybin];
  994. }
  995. //*/
  996. //__________________________________________________________________________
  997. void MpdTpcClusterFinderMlem::CreateHits(const vector<pixel> &pixels, multimap<Double_t,Int_t> &localMax,
  998. vector<vector<Double_t> > &charges, vector<vector<Int_t> > &flags,
  999. Int_t iclus, vector<multimap<Double_t,Int_t> > &pixInMax)
  1000. {
  1001. // Create hits from pixels
  1002. //Int_t nLocMax0 = localMax.size();
  1003. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  1004. TH2D *hXY = (TH2D*) gROOT->FindObject("hTimePad");
  1005. TH2D *hMlem = (TH2D*) gROOT->FindObject("hMlem1");
  1006. //if (hMlem == nullptr) hMlem = (TH2D*) gROOT->FindObject("hMlem1");
  1007. Double_t scale = hXY->GetXaxis()->GetBinWidth(1) / hMlem->GetXaxis()->GetBinWidth(1);
  1008. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin();
  1009. //vector<Int_t> vecDig;
  1010. map<Int_t,Double_t> mapIdQ;
  1011. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  1012. Int_t ihit = fHitArray->GetEntriesFast(), nHitsAdd = 0;
  1013. for ( ; rit != localMax.rend(); ++rit) {
  1014. Int_t ipix = rit->second;
  1015. Int_t ipad = pixels[ipix].iy - 1;
  1016. Int_t itime = pixels[ipix].ix - 1;
  1017. if (flags[ipad][itime] <= 0) continue; // merged local max
  1018. //vecDig.clear();
  1019. mapIdQ.clear();
  1020. set<pair<Int_t,Int_t> > selPix;
  1021. Double_t padMean = 0, timeMean = 0, adcTot = 0, sum2t = 0, sum2p = 0;
  1022. // Process complex cluster - start from maximum and go outward (up to 5 steps),
  1023. // adding pixels with respectively lower charges
  1024. for (Int_t idirp = -1; idirp < 2; idirp += 2) {
  1025. for (Int_t ip = 0; ip < 5; ++ip) {
  1026. //for (Int_t ip = 0; ip < 8; ++ip) {
  1027. if (idirp > 0 && ip == 0) continue;
  1028. Int_t ipsign = ip * idirp;
  1029. Int_t ip1 = ipad + ipsign;
  1030. if (ip1 < 0 || (UInt_t) ip1 >= flags.size()) break;
  1031. for (Int_t idirt = -1; idirt < 2; idirt += 2) {
  1032. for (Int_t it = 0; it < 5; ++it) {
  1033. if (idirt > 0 && it == 0) continue;
  1034. Int_t itsign = it * idirt;
  1035. Int_t it1 = itime + itsign;
  1036. if (it1 < 0 || (UInt_t) it1 >= (flags[0]).size()) break;
  1037. if (flags[ip1][it1] == 0) continue;
  1038. Int_t add = 1;
  1039. if (ip || it) {
  1040. if (flags[ip1][it1] > 0) continue; // case when 2 local max next to each other on 1 diagonal
  1041. // Check gradient
  1042. add = 0;
  1043. Int_t ipprev = ipsign, itprev = itsign;
  1044. if (it) {
  1045. itprev -= idirt;
  1046. Int_t it10 = itime + itprev;
  1047. if (it10 >= 0 && (UInt_t) it10 < (flags[0]).size() && flags[ip1][it10] != 0) {
  1048. if (selPix.find(pair<Int_t,Int_t>(ip1,it10)) != selPix.end()
  1049. && charges[ip1][it1] <= charges[ip1][it10]) add = 1;
  1050. }
  1051. }
  1052. if (add == 0 && ip) {
  1053. ipprev -= idirp;
  1054. Int_t ip10 = ipad + ipprev;
  1055. if (ip10 >= 0 && (UInt_t) ip10 < flags.size() && flags[ip10][it1] != 0) {
  1056. if (selPix.find(pair<Int_t,Int_t>(ip10,it1)) != selPix.end()
  1057. && charges[ip1][it1] <= charges[ip10][it1]) add = 1;
  1058. }
  1059. }
  1060. }
  1061. if (!add) break;
  1062. padMean += ip1 * charges[ip1][it1];
  1063. timeMean += it1 * charges[ip1][it1];
  1064. adcTot += charges[ip1][it1];
  1065. //vecDig.push_back(fDigis[ip1][it1]);
  1066. Int_t ipBin = hXY->GetYaxis()->FindBin(hMlem->GetYaxis()->GetBinCenter(ip1+1)); // original bin number
  1067. Int_t itBin = hXY->GetXaxis()->FindBin(hMlem->GetXaxis()->GetBinCenter(it1+1)); // original bin number
  1068. ipBin = TMath::Max (TMath::Nint (hXY->GetYaxis()->GetBinLowEdge(ipBin)), 0);
  1069. itBin = TMath::Nint (hXY->GetXaxis()->GetBinLowEdge(itBin));
  1070. Int_t id = clus->Sec(fDigis[ipBin][itBin]);
  1071. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = fCharges[ip1][it1];
  1072. //else mapIdQ[id] = mapIdQ[id] + fCharges[ip1][it1];
  1073. //else mapIdQ[id] = TMath::Max (mapIdQ[id], fCharges[ip1][it1]);
  1074. if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = 1;
  1075. else mapIdQ[id] = mapIdQ[id] + 1; // number of digits with the same ID
  1076. selPix.insert(pair<Int_t,Int_t>(ip1,it1));
  1077. sum2t += it1 * it1 * charges[ip1][it1];
  1078. sum2p += ip1 * ip1 * charges[ip1][it1];
  1079. pixInMax[nHitsAdd].insert(pair<Double_t,Int_t>(-charges[ip1][it1],flags[ip1][it1]));
  1080. }
  1081. }
  1082. }
  1083. }
  1084. //padMean = (padMean / adcTot + 0.5) / scale;
  1085. padMean = (padMean / adcTot + 0.5) / scale - 2; // coorrect for shift
  1086. timeMean = (timeMean / adcTot + 0.5) / scale;
  1087. Double_t xloc = fSecGeo->Pad2Xloc(padMean-0.5+hMlem->GetYaxis()->GetXmin(),clus->Row());
  1088. Int_t padID = fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  1089. Double_t yloc = fSecGeo->LocalPadPosition(padID).Y();
  1090. Double_t zloc = fSecGeo->TimeBin2Z(timeMean-0.5+hMlem->GetXaxis()->GetXmin());
  1091. p3loc.SetXYZ(xloc, yloc, zloc);
  1092. Double_t rmsZ = TMath::Sqrt (sum2t/adcTot/scale/scale - timeMean*timeMean);
  1093. Double_t rmsX = TMath::Sqrt (sum2p/adcTot/scale/scale - padMean*padMean);
  1094. //p3err[1] = fSecGeo->TimeBin2Z(rms); // to pass the value
  1095. //p3err[1] = rms;
  1096. //cout << " Result: " << nLocMax0 << " " << pixels.size() << " " << xloc << " " << zloc << endl;
  1097. // Apply corrections
  1098. TVector3 p3errCor(p3err);
  1099. CorrectRecoMlem(p3loc, p3errCor, clus, adcTot);
  1100. if (clus->GetSect() >= fSecGeo->NofSectors()) p3loc[2] = -p3loc[2];
  1101. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  1102. // Dip angle correction (for interaction point with Z = 0)
  1103. Double_t dip = TMath::ATan2 (TMath::Abs(p3glob[2]),p3glob.Pt()) * TMath::RadToDeg();
  1104. p3errCor[2] = TMath::Max (p3errCor[2], 0.116 - 0.0046 * dip + 0.00015 * dip * dip);
  1105. p3errCor[2] = TMath::Min (p3errCor[2], 0.5);
  1106. // Correct Z-coordinate for rows = 0 and 27
  1107. //if ((clus->Row() == 0 || clus->Row() == 27) && rmsZ > 1.0) {
  1108. if (clus->Row() == 0 && rmsZ > 1.0) {
  1109. Double_t zcor = 0;
  1110. if (clus->Row() == 0) zcor = -0.011 + 0.002 * dip;
  1111. else zcor = -0.029 + 0.005 * dip + 3.886e-5 * dip * dip;
  1112. zcor = TMath::Max (zcor, 0.);
  1113. zcor = TMath::Min (zcor, 0.6);
  1114. p3loc[2] -= TMath::Sign (zcor, p3loc[2]);
  1115. p3glob[2] -= TMath::Sign (zcor, p3glob[2]);
  1116. }
  1117. MpdTpcHit* hit = new ((*fHitArray)[ihit++]) MpdTpcHit(padID, p3glob, p3errCor, iclus);
  1118. hit->SetLayer(clus->Row());
  1119. hit->SetLocalPosition(p3loc); // point position
  1120. hit->SetEnergyLoss(adcTot);
  1121. Int_t ireg = (clus->Row() < fSecGeo->NofRowsReg(0)) ? 0 : 1;
  1122. Double_t step = fSecGeo->PadHeight(ireg);
  1123. hit->SetStep(step);
  1124. hit->SetModular(1); // modular geometry flag
  1125. hit->SetPad(Int_t(padMean-0.5+hMlem->GetYaxis()->GetXmin()));
  1126. hit->SetBin(Int_t(timeMean-0.5+hMlem->GetXaxis()->GetXmin()));
  1127. hit->SetRMS(rmsX, 0);
  1128. hit->SetRMS(rmsZ, 1);
  1129. hit->SetNdigits(-selPix.size()); // negative value
  1130. hit->SetFlags(clus);
  1131. ++nHitsAdd;
  1132. // !!! Warning: FairLinks are not persistent !!!
  1133. //hit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  1134. vector<Int_t> vecDig;
  1135. // Store maximum 3 track IDs with the highest charge contribution
  1136. multimap<Double_t,Int_t> mapDig;
  1137. for (map<Int_t,Double_t>::iterator mit = mapIdQ.begin(); mit != mapIdQ.end(); ++mit)
  1138. mapDig.insert(pair<Double_t,Int_t>(-mit->second,mit->first));
  1139. multimap<Double_t,Int_t>::iterator mit = mapDig.begin();
  1140. while (mit != mapDig.end() && vecDig.size() < 9) {
  1141. //while (mit != mapDig.end() && vecDig.size() < 1) {
  1142. vecDig.push_back(mit->second);
  1143. ++mit;
  1144. }
  1145. Int_t ndig = vecDig.size();
  1146. /*
  1147. for (Int_t idig = 0; idig < ndig; ++idig)
  1148. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, clus->Sec(vecDig[idig]))); // trackID stored in Sec
  1149. */
  1150. for (Int_t idig = 0; idig < ndig; ++idig) {
  1151. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, vecDig[idig], idig)); // weight = idig
  1152. hit->AddID(vecDig[idig]);
  1153. }
  1154. //cout << hit->GetNLinks() << " " << *(vecDig.begin()) << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(0).GetIndex() << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(hit->GetNLinks()-1).GetIndex() << endl;
  1155. } // for ( ; rit != localMax.rend();
  1156. }
  1157. //__________________________________________________________________________
  1158. void MpdTpcClusterFinderMlem::CorrectRecoMlem(TVector3 &p3loc, TVector3 &p3errCor, MpdTpc2dCluster *clus, Double_t adc)
  1159. {
  1160. // Correct reconstructed coordinates
  1161. Double_t xsec = (p3loc.Y() + fSecGeo->GetMinY()) * TMath::Tan(fSecGeo->Dphi()/2);
  1162. Double_t xloc = p3loc.X(), xloc0 = xloc, edge = 0.0; // distance to sector boundary
  1163. if (xloc < 0) edge = xloc + xsec;
  1164. else edge = xloc - xsec;
  1165. if (TMath::Abs(edge) > 1.1 || adc > 2000) return; // no corrections
  1166. Double_t adc2 = adc * adc;
  1167. Double_t adc3 = adc2 * adc;
  1168. if (clus->Edge()) {
  1169. // Max pad at the edge
  1170. Double_t corr = 3.886 - 0.007 * adc + 3.371e-6 * adc2 - 8.175e-11 * adc3 - 1.812e-13 * adc3 * adc;
  1171. corr /= 10; // cm
  1172. if (edge < 0) corr = -corr;
  1173. xloc -= corr;
  1174. p3errCor[0] = 0.15;
  1175. } else {
  1176. if (adc > 1000) return;
  1177. Double_t corr = -0.22 / 10; // cm
  1178. if (adc < 700) corr = 7.521 - 0.012 * adc - 9.085e-6 * adc2 + 1.365e-8 * adc3;
  1179. corr /= 10; // cm
  1180. if (edge < 0) corr = -corr;
  1181. xloc -= corr;
  1182. p3errCor[0] = 0.2;
  1183. }
  1184. if (TMath::Abs(xloc) > xsec) xloc = TMath::Sign(xsec,xloc0);
  1185. p3loc.SetX(xloc);
  1186. }
  1187. //__________________________________________________________________________
  1188. void MpdTpcClusterFinderMlem::ChargeMlem(Int_t nHits0, vector<pixel> &pixels, vector<pixel> &bins,
  1189. vector<multimap<Double_t,Int_t> > &pixInMax, const TMatrixD &cij, Double_t cijMin)
  1190. {
  1191. // Get correct charge of hits after MLEM
  1192. // (reduce number of pixels not to exceed number of time-pad bins)
  1193. Int_t nH = fHitArray->GetEntriesFast();
  1194. Double_t qTot = 0.0;
  1195. for (Int_t i = nHits0; i < nH; ++i) {
  1196. MpdTpcHit* hit = (MpdTpcHit*) fHitArray->UncheckedAt(i);
  1197. qTot += hit->GetQ();
  1198. }
  1199. Int_t nBinsTot = bins.size(), nBins = 0;//, nHits = pixInMax.size();
  1200. Double_t coef = nBinsTot / qTot;
  1201. vector<pixel> pixs;
  1202. multimap<Double_t,Int_t>::iterator it;
  1203. // Divide bins between hits proportionally to their charges
  1204. for (Int_t i = nH-1; i >= nHits0; --i) {
  1205. MpdTpcHit* hit = (MpdTpcHit*) fHitArray->UncheckedAt(i);
  1206. Int_t nb = TMath::Nint (hit->GetQ() * coef);
  1207. if (nb == 0) ++nb;
  1208. nBins += nb;
  1209. if (nBins > nBinsTot) nb -= (nBins - nBinsTot);
  1210. else if (i == nHits0 && nBins < nBinsTot) nb -= (nBins - nBinsTot);
  1211. Int_t jmax = i - nHits0, nsel = 0;
  1212. for (it = pixInMax[jmax].begin(); it != pixInMax[jmax].end(); ++it) {
  1213. ++nsel;
  1214. if (nsel > nb) break;
  1215. Int_t ipix = TMath::Abs(it->second) - 1;
  1216. //cout << nsel << " " << pixels[ipix].qq << " " << pixels[ipix].vis << endl; //AZ
  1217. pixs.push_back(pixels[ipix]);
  1218. pixs.back().qq = 1;
  1219. pixs.back().ix = ipix; // save pixel index
  1220. pixs.back().iy = i; // save hit index
  1221. pixs.back().sum = 0.0;
  1222. }
  1223. }
  1224. // Vector of active pairs (bin-pixel)
  1225. vector<pair<Int_t,Int_t> > pairs;
  1226. vector<pair<Int_t,Int_t> >::iterator vit;
  1227. Int_t npix = pixs.size();
  1228. Double_t qbins = 0; //AZ
  1229. for (Int_t ibin = 0; ibin < nBinsTot; ++ibin) {
  1230. for (Int_t i = 0; i < npix; ++i) if (cij(pixs[i].ix,ibin) > cijMin) pairs.push_back(pair<Int_t,Int_t>(i,ibin));
  1231. qbins += bins[ibin].qq; //AZ
  1232. }
  1233. // MLEM procedure
  1234. Int_t niter = 2;
  1235. for (Int_t iter = 0; iter < niter; ++iter) {
  1236. // Calculate expectations for all pads
  1237. for (Int_t ibin = 0; ibin < nBinsTot; ++ibin) {
  1238. bins[ibin].vis = 0.0;
  1239. //for (Int_t i = 0; i < npix; ++i) bins[ibin].vis += (pixs[i].qq * cij(pixs[i].ix,ibin)); // bins[ibin].vis - just storage
  1240. //?? bins[ibin].vis = TMath::Min (bins[ibin].vis, fglOvfw);
  1241. }
  1242. for (vit = pairs.begin(); vit != pairs.end(); ++vit) {
  1243. Int_t i = vit->first, ibin = vit->second;
  1244. bins[ibin].vis += (pixs[i].qq * cij(pixs[i].ix,ibin)); // bins[ibin].vis - just storage
  1245. }
  1246. // Maximization
  1247. /*
  1248. for (Int_t ipix = 0; ipix < npix; ++ipix) {
  1249. Double_t sum = 0.0;
  1250. for (Int_t ibin = 0; ibin < nBinsTot; ++ibin) {
  1251. if (bins[ibin].vis > 1.e-6) sum += (bins[ibin].qq / bins[ibin].vis * cij(pixs[ipix].ix,ibin));
  1252. }
  1253. if (pixs[ipix].vis > 1.e-6) pixs[ipix].qq *= (sum / pixs[ipix].vis);
  1254. //if (pixels[icur][ipix].vis > 1.e-6) pixels[icur][ipix].qq *= (sum / visMax);
  1255. }
  1256. */
  1257. for (Int_t ipix = 0; ipix < npix; ++ipix) pixs[ipix].sum = 0.0;
  1258. for (vit = pairs.begin(); vit != pairs.end(); ++vit) {
  1259. Int_t ipix = vit->first, ibin = vit->second;
  1260. //AZ if (bins[ibin].vis > 1.e-6) pixs[ipix].sum += (bins[ibin].qq / bins[ibin].vis * cij(pixs[ipix].ix,ibin));
  1261. if (bins[ibin].vis > 1.e-6 && bins[ibin].qq < fgkOvfw - 1.5) pixs[ipix].sum += (bins[ibin].qq / bins[ibin].vis * cij(pixs[ipix].ix,ibin)); //AZ - 030620
  1262. }
  1263. for (Int_t ipix = 0; ipix < npix; ++ipix)
  1264. if (pixs[ipix].vis > 1.e-6) pixs[ipix].qq *= (pixs[ipix].sum / pixs[ipix].vis);
  1265. //if (pixs[ipix].vis > 1.e-6) pixs[ipix].qq *= pixs[ipix].sum; //AZ 03-06-20
  1266. } // for (Int_t iter = 0;
  1267. // Compute hit charges
  1268. map<Int_t,Double_t> charges;
  1269. for (Int_t ipix = 0; ipix < npix; ++ipix) {
  1270. if (charges.find(pixs[ipix].iy) == charges.end()) charges[pixs[ipix].iy] = pixs[ipix].qq;
  1271. else charges[pixs[ipix].iy] += pixs[ipix].qq;
  1272. }
  1273. Double_t qCor = 0; //AZ
  1274. for (Int_t ih = nHits0; ih < nH; ++ih) {
  1275. MpdTpcHit* hit = (MpdTpcHit*) fHitArray->UncheckedAt(ih);
  1276. //cout << hit->GetQ() << " " << charges[ih] << " " << hit->GetNdigits() << " " << hit->GetLayer() << " " << nHits << endl;
  1277. qCor += charges[ih]; //AZ
  1278. //AZ hit->SetEnergyLoss(charges[ih]/1.089); // coefficient due to different peak with and w/out MLEM
  1279. //hit->SetEnergyLoss(charges[ih]/1.0); //AZ - 040620 coefficient due to different peak with and w/out MLEM
  1280. hit->SetEnergyLoss(charges[ih]*1.25); //AZ - 120620 coefficient due to different peak with and w/out MLEM
  1281. }
  1282. //cout << " qtot: " << qTot << " " << qCor/qTot << " " << (nH-nHits0) << " " << nBinsTot << " " << qbins << " " << pixels.size() << endl; //AZ
  1283. }
  1284. //__________________________________________________________________________
  1285. ClassImp(MpdTpcClusterFinderMlem)