MpdTpcClusterFinderAZ.cxx 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class MpdTpcClusterFinderAZ
  7. // see MpdTpcClusterFinderAZ.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 - 14-July-2015
  14. //
  15. //-----------------------------------------------------------
  16. // This Class' Header ------------------
  17. #include "MpdTpcClusterFinderAZ.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 "TpcGas.h"
  25. //#include "TpcPrimaryCluster.h"
  26. //#include "LinearInterpolPolicy.h"
  27. #include "FairRootManager.h"
  28. #include "TClonesArray.h"
  29. #include "TFile.h"
  30. #include "TMath.h"
  31. #include "TRandom.h"
  32. #include "TSystem.h"
  33. #include "TTree.h"
  34. // C/C++ Headers ----------------------
  35. #include <iostream>
  36. #include <math.h>
  37. #include <set>
  38. #include <vector>
  39. using namespace std;
  40. //__________________________________________________________________________
  41. MpdTpcClusterFinderAZ::MpdTpcClusterFinderAZ()
  42. : FairTask("TPC Cluster finder AZ"), fPersistence(kFALSE)
  43. {
  44. /*
  45. std::string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  46. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  47. fGas= new TpcGas(tpcGasFile, 130);
  48. std::cout<<*fGas<<std::endl;
  49. */
  50. }
  51. //__________________________________________________________________________
  52. MpdTpcClusterFinderAZ::~MpdTpcClusterFinderAZ()
  53. {
  54. //delete fGas;
  55. }
  56. //__________________________________________________________________________
  57. void MpdTpcClusterFinderAZ::FinishTask()
  58. {
  59. for (Int_t i = 0; i < fgkNsec2; ++i) {
  60. //fPrimArray[i]->Delete();
  61. delete [] fDigiSet[i];
  62. }
  63. }
  64. //__________________________________________________________________________
  65. InitStatus MpdTpcClusterFinderAZ::Init()
  66. {
  67. // Create containers for digits
  68. fSecGeo = MpdTpcSectorGeo::Instance();
  69. //Int_t nRows = MpdTpcSectorGeo::Instance()->NofRows();
  70. Int_t nRows = fSecGeo->NofRows();
  71. for (Int_t i = 0; i < fgkNsec2; ++i) fDigiSet[i] = new set<Int_t> [nRows];
  72. //Get ROOT Manager
  73. FairRootManager* ioman = FairRootManager::Instance();
  74. if (ioman == 0) {
  75. Error("MpdTpcClusterFinderAZ::Init","RootManager not instantiated!");
  76. return kERROR;
  77. }
  78. // Get input collection
  79. fDigiArray = (TClonesArray*) ioman->GetObject("MpdTpcDigit");
  80. if (fDigiArray == 0) {
  81. Error("TpcClusterFinderAZ::Init","Array of digits not found!");
  82. return kERROR;
  83. }
  84. // Create and register output array
  85. //SetPersistence();
  86. fClusArray = new TClonesArray("MpdTpc2dCluster");
  87. ioman->Register("TpcCluster", "Tpc", fClusArray, fPersistence);
  88. fHitArray = new TClonesArray("MpdTpcHit");
  89. ioman->Register("TpcRecPoint", "Tpc", fHitArray, fPersistence);
  90. return kSUCCESS;
  91. }
  92. //__________________________________________________________________________
  93. void MpdTpcClusterFinderAZ::Exec(Option_t* opt)
  94. {
  95. fClusArray->Delete();
  96. fHitArray->Delete();
  97. const Int_t nSec = fgkNsec2 / 2; // number of TPC readout sectors
  98. // Clear digi containers
  99. //Int_t nRows = MpdTpcSectorGeo::Instance()->NofRows();
  100. Int_t nRows = fSecGeo->NofRows();
  101. for (Int_t i = 0; i < fgkNsec2; ++i) {
  102. for (Int_t j = 0; j < nRows; ++j) fDigiSet[i][j].clear();
  103. }
  104. // Fill digi containers
  105. Int_t nDigis = fDigiArray->GetEntriesFast();
  106. cout << " Total number of digits: " << nDigis << endl;
  107. for (Int_t i = 0; i < nDigis; ++i) {
  108. MpdTpcDigit *digi = (MpdTpcDigit*) fDigiArray->UncheckedAt(i);
  109. Int_t isec = digi->GetSector();
  110. Int_t irow = digi->GetRow();
  111. fDigiSet[isec][irow].insert(i);
  112. }
  113. Int_t nSum = 0;
  114. // Loop over sectors
  115. for (Int_t isec = 0; isec < fgkNsec2; ++isec) {
  116. // Loop over padrows
  117. for (Int_t irow = 0; irow < nRows; ++irow) {
  118. if (fDigiSet[isec][irow].size() == 0) continue;
  119. //cout << " Sector, row, digits: " << isec << " " << irow << " " << fDigiSet[isec][irow].size() << endl;
  120. ProcessPadrow(isec, irow);
  121. nSum += fDigiSet[isec][irow].size();
  122. }
  123. }
  124. cout << " Control sum: " << nSum << endl;
  125. // Find hits
  126. FindHits();
  127. //FindHitsLocMax();
  128. }
  129. //__________________________________________________________________________
  130. void MpdTpcClusterFinderAZ::ProcessPadrow(Int_t isec, Int_t irow)
  131. {
  132. // Process one padrow of a sector
  133. memset(fFlags, 0, sizeof(fFlags[0][0]) * fgkNpads * fgkNtimes);
  134. set<Int_t>::iterator it = fDigiSet[isec][irow].begin();
  135. for ( ; it != fDigiSet[isec][irow].end(); ++it) {
  136. MpdTpcDigit *digi = (MpdTpcDigit*) fDigiArray->UncheckedAt(*it);
  137. Int_t ipad = digi->GetPad(), itime = digi->GetTimeBin();
  138. if (ipad >= fgkNpads) Fatal("ProcessPadrow", "Too few pads!!! %i %i", ipad, fgkNpads);
  139. if (itime >= fgkNtimes) Fatal("ProcessPadrow", "Too few time bins!!! %i %i", itime, fgkNtimes);
  140. fCharges[ipad][itime] = digi->GetAdc();
  141. //fDigis[ipad][itime] = *it;
  142. fDigis[ipad][itime] = digi->GetOrigin(); // store trackID with max charge contribution
  143. fFlags[ipad][itime] = 1;
  144. }
  145. // Find (pre)clusters
  146. Int_t nclus0 = fClusArray->GetEntriesFast(), nclus = nclus0;
  147. for (Int_t ipad = 0; ipad < fgkNpads; ++ipad) {
  148. for (Int_t itime = 0; itime < fgkNtimes; ++itime) {
  149. if (fFlags[ipad][itime] <= 0) continue;
  150. // New cluster
  151. MpdTpc2dCluster* clus = new ((*fClusArray)[nclus++]) MpdTpc2dCluster(irow, isec);
  152. clus->Insert(fDigis[ipad][itime], irow, ipad, itime, fCharges[ipad][itime]);
  153. clus->SetID(fDigis[ipad][itime]); // trackID
  154. clus->SetErrY(fCharges[ipad][itime]); // ADC counts
  155. fFlags[ipad][itime] = -1;
  156. for (Int_t ip = 0; ip < 2; ++ip) {
  157. for (Int_t itt = 0; itt < 2; ++itt) {
  158. if (itt == ip) continue;
  159. Int_t ip1 = ipad + ip, it1 = itime + itt;
  160. if (ip1 >= fgkNpads) continue;
  161. if (it1 >= fgkNtimes) continue;
  162. if (fFlags[ip1][it1] <= 0) continue;
  163. NextPixel(clus, ip1, it1);
  164. }
  165. }
  166. }
  167. }
  168. //cout << " Found preclusters: " << nclus - nclus0 << endl;
  169. }
  170. //__________________________________________________________________________
  171. void MpdTpcClusterFinderAZ::NextPixel(MpdTpc2dCluster* clus, Int_t ipad, Int_t itime)
  172. {
  173. // Add next pixel to the cluster
  174. clus->Insert(fDigis[ipad][itime], clus->Row(), ipad, itime, fCharges[ipad][itime]);
  175. clus->SetID (TMath::Min(clus->ID(),fDigis[ipad][itime])); // min trackID
  176. clus->SetErrY (TMath::Max(Double_t(clus->GetErrY()),fCharges[ipad][itime])); // max ADC counts
  177. fFlags[ipad][itime] = -1;
  178. for (Int_t ip = -1; ip < 2; ++ip) {
  179. for (Int_t it = -1; it < 2; ++it) {
  180. if (TMath::Abs(it) == TMath::Abs(ip)) continue;
  181. Int_t ip1 = ipad + ip, it1 = itime + it;
  182. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  183. if (it1 < 0 || it1 >= fgkNtimes) continue;
  184. if (fFlags[ip1][it1] <= 0) continue;
  185. NextPixel(clus, ip1, it1);
  186. }
  187. }
  188. }
  189. //__________________________________________________________________________
  190. /*
  191. void MpdTpcClusterFinderAZ::FindHits()
  192. {
  193. // Reconstruct hits (one hit per precluster)
  194. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  195. Int_t nclus = fClusArray->GetEntriesFast();
  196. for (Int_t iclus = 0; iclus < nclus; ++iclus) {
  197. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  198. Int_t nDigis = clus->NDigits();
  199. Double_t padMean = 0, timeMean = 0;
  200. for (Int_t idig = 0; idig < nDigis; ++idig) {
  201. padMean += clus->Col(idig) * clus->Adc(idig);
  202. timeMean += clus->Bkt(idig) * clus->Adc(idig);
  203. }
  204. padMean /= clus->GetADC();
  205. timeMean /= clus->GetADC();
  206. Double_t xloc = fSecGeo->Pad2Xloc(padMean,clus->Row());
  207. Int_t padID = fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  208. Double_t yloc = fSecGeo->LocalPadPosition(padID).Y();
  209. Double_t zloc = fSecGeo->TimeBin2Z(timeMean);
  210. p3loc.SetXYZ(xloc, yloc, zloc);
  211. // Apply corrections
  212. TVector3 p3errCor(p3err);
  213. CorrectReco(p3loc, p3errCor, clus->GetNumPads(), clus->GetADC());
  214. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  215. if (clus->GetSect() >= fSecGeo->NofSectors()) p3glob[2] = -p3glob[2];
  216. MpdTpcHit* hit = new ((*fHitArray)[iclus]) MpdTpcHit(padID, p3glob, p3errCor, iclus);
  217. hit->SetLayer(clus->Row());
  218. hit->SetLocalPosition(p3loc); // point position
  219. hit->SetEnergyLoss(clus->GetADC());
  220. hit->SetStep(0.0);
  221. hit->SetModular(1); // modular geometry flag
  222. hit->SetPad(Int_t(padMean));
  223. hit->SetBin(Int_t(timeMean));
  224. // !!! Warning: FairLinks are not persistent !!!
  225. //hit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  226. for (Int_t idig = 0; idig < nDigis; ++idig)
  227. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, clus->Sec(idig))); // trackID stored in Sec
  228. //cout << hit->GetNLinks();
  229. }
  230. }
  231. */
  232. //__________________________________________________________________________
  233. /*
  234. void MpdTpcClusterFinderAZ::FindHitsLocMax()
  235. {
  236. // Reconstruct hits (find local maxima)
  237. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  238. Int_t nclus = fClusArray->GetEntriesFast(), ihit = 0;
  239. for (Int_t iclus = 0; iclus < nclus; ++iclus) {
  240. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  241. Int_t nDigis = clus->NDigits();
  242. memset(fFlags, 0, sizeof(fFlags[0][0]) * fgkNpads * fgkNtimes);
  243. for (Int_t idig = 0; idig < nDigis; ++idig) {
  244. Int_t ipad = clus->Col(idig);
  245. Int_t itime = clus->Bkt(idig);
  246. fCharges[ipad][itime] = clus->Adc(idig);
  247. fFlags[ipad][itime] = 1;
  248. fDigis[ipad][itime] = idig;
  249. }
  250. // Exclude pads which are not local maxima
  251. for (Int_t idig = 0; idig < nDigis; ++idig) {
  252. Int_t ipad = clus->Col(idig);
  253. Int_t itime = clus->Bkt(idig);
  254. for (Int_t ip = -1; ip < 2; ++ip) {
  255. for (Int_t it = -1; it < 2; ++it) {
  256. //if (TMath::Abs(it) == TMath::Abs(ip)) continue; // exclude diagonals
  257. if (it == 0 && ip == 0) continue;
  258. Int_t ip1 = ipad + ip, it1 = itime + it;
  259. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  260. if (it1 < 0 || it1 >= fgkNtimes) continue;
  261. if (fFlags[ip1][it1] == 0) continue;
  262. if (clus->Adc(idig) < fCharges[ip1][it1]) { fFlags[ipad][itime] = -1; break; }
  263. }
  264. }
  265. }
  266. multimap<Double_t,Int_t> localMax;
  267. for (Int_t idig = 0; idig < nDigis; ++idig) {
  268. Int_t ipad = clus->Col(idig);
  269. Int_t itime = clus->Bkt(idig);
  270. if (fFlags[ipad][itime] <= 0) continue;
  271. localMax.insert(pair<Double_t,Int_t>(clus->Adc(idig),idig));
  272. cout << clus->Col(idig) << " " << clus->Bkt(idig) << " " << clus->Adc(idig) << endl;
  273. }
  274. cout << " Local max: " << localMax.size() << endl;
  275. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin(), rit1 = rit;
  276. vector<Int_t> vecDig;
  277. for ( ; rit != localMax.rend(); ++rit) {
  278. Int_t idig = rit->second;
  279. Int_t ipad = clus->Col(idig);
  280. Int_t itime = clus->Bkt(idig);
  281. if (fFlags[ipad][itime] <= 0) continue; // merged local max
  282. vecDig.clear();
  283. Double_t padMean = 0, timeMean = 0, adcTot = 0;
  284. for (Int_t ip = -1; ip < 2; ++ip) {
  285. for (Int_t it = -1; it < 2; ++it) {
  286. //if (TMath::Abs(it) == TMath::Abs(ip)) continue;
  287. Int_t ip1 = ipad + ip, it1 = itime + it;
  288. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  289. if (it1 < 0 || it1 >= fgkNtimes) continue;
  290. if (fFlags[ip1][it1] == 0) continue;
  291. if (fFlags[ip1][it1] > 0) fFlags[ip1][it1] = -1; // merge neighbour local max
  292. padMean += ip1 * fCharges[ip1][it1];
  293. timeMean += it1 * fCharges[ip1][it1];
  294. adcTot += fCharges[ip1][it1];
  295. vecDig.push_back(fDigis[ip1][it1]);
  296. }
  297. }
  298. padMean /= adcTot;
  299. timeMean /= adcTot;
  300. Double_t xloc = fSecGeo->Pad2Xloc(padMean,clus->Row());
  301. Int_t padID = fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  302. Double_t yloc = fSecGeo->LocalPadPosition(padID).Y();
  303. Double_t zloc = fSecGeo->TimeBin2Z(timeMean);
  304. p3loc.SetXYZ(xloc, yloc, zloc);
  305. // Apply corrections
  306. TVector3 p3errCor(p3err);
  307. CorrectReco(p3loc, p3errCor, clus->GetNumPads(), adcTot);
  308. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  309. if (clus->GetSect() >= fSecGeo->NofSectors()) p3glob[2] = -p3glob[2];
  310. MpdTpcHit* hit = new ((*fHitArray)[ihit++]) MpdTpcHit(padID, p3glob, p3errCor, iclus);
  311. hit->SetLayer(clus->Row());
  312. hit->SetLocalPosition(p3loc); // point position
  313. hit->SetEnergyLoss(adcTot);
  314. hit->SetStep(0.0);
  315. hit->SetModular(1); // modular geometry flag
  316. hit->SetPad(Int_t(padMean));
  317. hit->SetBin(Int_t(timeMean));
  318. // !!! Warning: FairLinks are not persistent !!!
  319. //hit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  320. Int_t ndig = vecDig.size();
  321. for (Int_t idig = 0; idig < ndig; ++idig)
  322. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, clus->Sec(vecDig[idig]))); // trackID stored in Sec
  323. //cout << hit->GetNLinks();
  324. } // for ( ; rit != localMax.rend();
  325. } // for (Int_t iclus = 0; iclus < nclus;
  326. }
  327. */
  328. //__________________________________________________________________________
  329. void MpdTpcClusterFinderAZ::FindHits()
  330. {
  331. // Reconstruct hits (find local maxima)
  332. /*
  333. fCharges[ipad][itime] = digi->GetAdc();
  334. //fDigis[ipad][itime] = *it;
  335. fDigis[ipad][itime] = digi->GetOrigin(); // store trackID with max charge contribution
  336. fFlags[ipad][itime] = 1;
  337. */
  338. TVector3 p3loc, p3glob, p3err(0.05,0.0,0.1);
  339. Int_t nclus = fClusArray->GetEntriesFast(), ihit = 0;
  340. for (Int_t iclus = 0; iclus < nclus; ++iclus) {
  341. MpdTpc2dCluster* clus = (MpdTpc2dCluster*) fClusArray->UncheckedAt(iclus);
  342. Int_t nDigis = clus->NDigits();
  343. memset(fFlags, 0, sizeof(fFlags[0][0]) * fgkNpads * fgkNtimes);
  344. for (Int_t idig = 0; idig < nDigis; ++idig) {
  345. Int_t ipad = clus->Col(idig);
  346. Int_t itime = clus->Bkt(idig);
  347. fCharges[ipad][itime] = clus->Adc(idig);
  348. fFlags[ipad][itime] = 1;
  349. fDigis[ipad][itime] = idig;
  350. }
  351. // Exclude pads which are not local maxima
  352. for (Int_t idig = 0; idig < nDigis; ++idig) {
  353. Int_t ipad = clus->Col(idig);
  354. Int_t itime = clus->Bkt(idig);
  355. for (Int_t ip = -1; ip < 2; ++ip) {
  356. for (Int_t it = -1; it < 2; ++it) {
  357. if (TMath::Abs(it) == TMath::Abs(ip)) continue; // exclude diagonals
  358. if (it == 0 && ip == 0) continue;
  359. Int_t ip1 = ipad + ip, it1 = itime + it;
  360. if (ip1 < 0 || ip1 >= fgkNpads) continue;
  361. if (it1 < 0 || it1 >= fgkNtimes) continue;
  362. if (fFlags[ip1][it1] == 0) continue;
  363. if (clus->Adc(idig) < fCharges[ip1][it1]) { fFlags[ipad][itime] = -1; break; }
  364. }
  365. }
  366. }
  367. multimap<Double_t,Int_t> localMax;
  368. for (Int_t idig = 0; idig < nDigis; ++idig) {
  369. Int_t ipad = clus->Col(idig);
  370. Int_t itime = clus->Bkt(idig);
  371. if (fFlags[ipad][itime] <= 0) continue;
  372. localMax.insert(pair<Double_t,Int_t>(clus->Adc(idig),idig));
  373. //cout << clus->Col(idig) << " " << clus->Bkt(idig) << " " << clus->Adc(idig) << endl;
  374. }
  375. //cout << " Local max: " << clus->GetSect() << " " << clus->Row() << " " << localMax.size() << endl;
  376. // Remove local maxima not separated by valleys - Peak and Valley
  377. Int_t nLocMax0 = localMax.size();
  378. if (localMax.size() > 1) PeakAndValley(clus, localMax);
  379. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin();
  380. //vector<Int_t> vecDig;
  381. map<Int_t,Double_t> mapIdQ;
  382. for ( ; rit != localMax.rend(); ++rit) {
  383. Int_t idig = rit->second;
  384. Int_t ipad = clus->Col(idig);
  385. Int_t itime = clus->Bkt(idig);
  386. if (fFlags[ipad][itime] <= 0) continue; // merged local max
  387. //vecDig.clear();
  388. mapIdQ.clear();
  389. set<pair<Int_t,Int_t> > pixels;
  390. Double_t padMean = 0, timeMean = 0, adcTot = 0, sum2t = 0, sum2p = 0;
  391. // Process simple cluster (only 1 local max)
  392. //if (localMax.size() == 1) {
  393. if (nLocMax0 == 1) {
  394. for (Int_t idig1 = 0; idig1 < nDigis; ++idig1) {
  395. Int_t ip = clus->Col(idig1);
  396. Int_t it = clus->Bkt(idig1);
  397. padMean += ip * fCharges[ip][it];
  398. timeMean += it * fCharges[ip][it];
  399. adcTot += fCharges[ip][it];
  400. //vecDig.push_back(fDigis[ip][it]);
  401. Int_t id = clus->Sec(fDigis[ip][it]);
  402. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = fCharges[ip][it];
  403. //else mapIdQ[id] = mapIdQ[id] + fCharges[ip][it];
  404. //else mapIdQ[id] = TMath::Max (mapIdQ[id], fCharges[ip][it]);
  405. if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = 1;
  406. else mapIdQ[id] = mapIdQ[id] + 1;
  407. sum2t += it * it * fCharges[ip][it];
  408. sum2p += ip * ip * fCharges[ip][it];
  409. }
  410. } else {
  411. // Process complex cluster - start from maximum and go outward (up to 5 steps),
  412. // adding pixels with respectively lower charges
  413. for (Int_t idirp = -1; idirp < 2; idirp += 2) {
  414. for (Int_t ip = 0; ip < 5; ++ip) {
  415. if (idirp > 0 && ip == 0) continue;
  416. Int_t ipsign = ip * idirp;
  417. Int_t ip1 = ipad + ipsign;
  418. if (ip1 < 0 || ip1 >= fgkNpads) break;
  419. for (Int_t idirt = -1; idirt < 2; idirt += 2) {
  420. for (Int_t it = 0; it < 5; ++it) {
  421. if (idirt > 0 && it == 0) continue;
  422. Int_t itsign = it * idirt;
  423. Int_t it1 = itime + itsign;
  424. if (it1 < 0 || it1 >= fgkNtimes) break;
  425. if (fFlags[ip1][it1] == 0) continue;
  426. Int_t add = 1;
  427. if (ip || it) {
  428. if (fFlags[ip1][it1] > 0) continue; // case when 2 local max next to each other on 1 diagonal
  429. // Check gradient
  430. add = 0;
  431. Int_t ipprev = ipsign, itprev = itsign;
  432. if (it) {
  433. itprev -= idirt;
  434. Int_t it10 = itime + itprev;
  435. if (it10 >= 0 && it10 < fgkNtimes && fFlags[ip1][it10] != 0) {
  436. if (pixels.find(pair<Int_t,Int_t>(ip1,it10)) != pixels.end()
  437. && fCharges[ip1][it1] <= fCharges[ip1][it10]) add = 1;
  438. }
  439. }
  440. if (add == 0 && ip) {
  441. ipprev -= idirp;
  442. Int_t ip10 = ipad + ipprev;
  443. if (ip10 >= 0 && ip10 < fgkNpads && fFlags[ip10][it1] != 0) {
  444. if (pixels.find(pair<Int_t,Int_t>(ip10,it1)) != pixels.end()
  445. && fCharges[ip1][it1] <= fCharges[ip10][it1]) add = 1;
  446. }
  447. }
  448. }
  449. if (!add) break;
  450. padMean += ip1 * fCharges[ip1][it1];
  451. timeMean += it1 * fCharges[ip1][it1];
  452. adcTot += fCharges[ip1][it1];
  453. //vecDig.push_back(fDigis[ip1][it1]);
  454. Int_t id = clus->Sec(fDigis[ip1][it1]);
  455. //if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = fCharges[ip1][it1];
  456. //else mapIdQ[id] = mapIdQ[id] + fCharges[ip1][it1];
  457. //else mapIdQ[id] = TMath::Max (mapIdQ[id], fCharges[ip1][it1]);
  458. if (mapIdQ.find(id) == mapIdQ.end()) mapIdQ[id] = 1;
  459. else mapIdQ[id] = mapIdQ[id] + 1; // number of digits with the same ID
  460. pixels.insert(pair<Int_t,Int_t>(ip1,it1));
  461. sum2t += it1 * it1 * fCharges[ip1][it1];
  462. sum2p += ip1 * ip1 * fCharges[ip1][it1];
  463. }
  464. }
  465. }
  466. }
  467. }
  468. padMean /= adcTot;
  469. timeMean /= adcTot;
  470. Double_t xloc = fSecGeo->Pad2Xloc(padMean,clus->Row());
  471. Int_t padID = fSecGeo->PadID(clus->GetSect() % fSecGeo->NofSectors(), clus->Row());
  472. Double_t yloc = fSecGeo->LocalPadPosition(padID).Y();
  473. Double_t zloc = fSecGeo->TimeBin2Z(timeMean);
  474. p3loc.SetXYZ(xloc, yloc, zloc);
  475. Double_t rmsZ = TMath::Sqrt (sum2t/adcTot - timeMean*timeMean);
  476. Double_t rmsX = TMath::Sqrt (sum2p/adcTot - padMean*padMean);
  477. //p3err[1] = fSecGeo->TimeBin2Z(rms); // to pass the value
  478. //p3err[1] = rms;
  479. //cout << " Result: " << nLocMax0 << " " << pixels.size() << " " << xloc << " " << zloc << endl;
  480. // Apply corrections
  481. TVector3 p3errCor(p3err);
  482. CorrectReco(p3loc, p3errCor, clus->GetNumPads(), adcTot);
  483. if (clus->GetSect() >= fSecGeo->NofSectors()) p3loc[2] = -p3loc[2];
  484. fSecGeo->Local2Global(clus->GetSect(), p3loc, p3glob);
  485. // Dip angle correction (for interaction point with Z = 0)
  486. Double_t dip = TMath::ATan2 (TMath::Abs(p3glob[2]),p3glob.Pt()) * TMath::RadToDeg();
  487. p3errCor[2] = TMath::Max (p3errCor[2], 0.116 - 0.0046 * dip + 0.00015 * dip * dip);
  488. p3errCor[2] = TMath::Min (p3errCor[2], 0.5);
  489. // Correct Z-coordinate for rows = 0 and 27
  490. //if ((clus->Row() == 0 || clus->Row() == 27) && rmsZ > 1.0) {
  491. if (clus->Row() == 0 && rmsZ > 1.0) {
  492. Double_t zcor = 0;
  493. if (clus->Row() == 0) zcor = -0.011 + 0.002 * dip;
  494. else zcor = -0.029 + 0.005 * dip + 3.886e-5 * dip * dip;
  495. zcor = TMath::Max (zcor, 0.);
  496. zcor = TMath::Min (zcor, 0.6);
  497. p3loc[2] -= TMath::Sign (zcor, p3loc[2]);
  498. p3glob[2] -= TMath::Sign (zcor, p3glob[2]);
  499. }
  500. MpdTpcHit* hit = new ((*fHitArray)[ihit++]) MpdTpcHit(padID, p3glob, p3errCor, iclus);
  501. hit->SetLayer(clus->Row());
  502. hit->SetLocalPosition(p3loc); // point position
  503. hit->SetEnergyLoss(adcTot);
  504. Int_t ireg = (clus->Row() < fSecGeo->NofRowsReg(0)) ? 0 : 1;
  505. Double_t step = fSecGeo->PadHeight(ireg);
  506. hit->SetStep(step);
  507. hit->SetModular(1); // modular geometry flag
  508. hit->SetPad(Int_t(padMean));
  509. hit->SetBin(Int_t(timeMean));
  510. hit->SetRMS(rmsX, 0);
  511. hit->SetRMS(rmsZ, 1);
  512. if (nLocMax0 == 1) hit->SetNdigits(nDigis);
  513. else hit->SetNdigits(pixels.size());
  514. // !!! Warning: FairLinks are not persistent !!!
  515. //hit->AddLink(FairLink(MpdTpcHit::PointIndex, pointIndx));
  516. vector<Int_t> vecDig;
  517. // Store maximum 3 track IDs with the highest charge contribution
  518. multimap<Double_t,Int_t> mapDig;
  519. for (map<Int_t,Double_t>::iterator mit = mapIdQ.begin(); mit != mapIdQ.end(); ++mit)
  520. mapDig.insert(pair<Double_t,Int_t>(-mit->second,mit->first));
  521. multimap<Double_t,Int_t>::iterator mit = mapDig.begin();
  522. while (mit != mapDig.end() && vecDig.size() < 9) {
  523. //while (mit != mapDig.end() && vecDig.size() < 1) {
  524. vecDig.push_back(mit->second);
  525. ++mit;
  526. }
  527. Int_t ndig = vecDig.size();
  528. /*
  529. for (Int_t idig = 0; idig < ndig; ++idig)
  530. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, clus->Sec(vecDig[idig]))); // trackID stored in Sec
  531. */
  532. for (Int_t idig1 = 0; idig1 < ndig; ++idig1) {
  533. hit->AddLink(FairLink(MpdTpcHit::MCTrackIndex, vecDig[idig], idig1)); // weight = idig
  534. hit->AddID(vecDig[idig1]);
  535. }
  536. //cout << hit->GetNLinks() << " " << *(vecDig.begin()) << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(0).GetIndex() << " " << hit->GetLinksWithType(MpdTpcHit::MCTrackIndex).GetLink(hit->GetNLinks()-1).GetIndex() << endl;
  537. } // for ( ; rit != localMax.rend();
  538. } // for (Int_t iclus = 0; iclus < nclus;
  539. }
  540. //__________________________________________________________________________
  541. void MpdTpcClusterFinderAZ::PeakAndValley(const MpdTpc2dCluster* clus, multimap<Double_t,Int_t> &localMax)
  542. {
  543. // Apply peak-and-valley cuts to remove some local maxima
  544. const Double_t ratio = 1.3;
  545. multimap<Double_t,Int_t>::reverse_iterator rit = localMax.rbegin(), rit1;
  546. ++rit;
  547. for ( ; rit != localMax.rend(); ++rit) {
  548. Int_t idig = rit->second;
  549. Int_t ipad = clus->Col(idig);
  550. Int_t itime = clus->Bkt(idig);
  551. if (fFlags[ipad][itime] <= 0) continue; // merged local max
  552. // Loop over higher peaks
  553. rit1 = localMax.rbegin();
  554. for ( ; rit1 != rit; ++rit1) {
  555. Int_t idig0 = rit1->second;
  556. Int_t ipad0 = clus->Col(idig0);
  557. Int_t itime0 = clus->Bkt(idig0);
  558. if (fFlags[ipad0][itime0] <= 0) continue; // merged local max
  559. Int_t dpad = ipad - ipad0, dt = itime - itime0;
  560. Int_t i0 = itime0, i1 = itime, j0 = ipad0, j1 = ipad, intime = 1;
  561. if (TMath::Abs(dpad) > TMath::Abs(dt)) i0 = ipad0, i1 = ipad, j0 = itime0, j1 = itime, intime = 0;
  562. Int_t stepi = TMath::Sign(1,i1-i0), stepj = TMath::Sign(1,j1-j0);
  563. //Int_t valOk = 0;
  564. Int_t merge = 0;
  565. //if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) merge = 1;
  566. //else {
  567. if (TMath::Abs(dpad) <= 1 && TMath::Abs(dt) <= 1) {
  568. if (fCharges[ipad-dpad][itime] > fCharges[ipad][itime] / ratio ||
  569. fCharges[ipad][itime-dt] > fCharges[ipad][itime] / ratio) merge = 1;
  570. } else {
  571. for (Int_t ii = i0 + stepi; ii != i1; ii += stepi) {
  572. merge = 0;
  573. for (Int_t jj = j0; jj != j1 + stepj; jj += stepj) {
  574. if (TMath::Abs(jj-j0) > TMath::Abs(ii-i0)) continue;
  575. if (TMath::Abs(jj-j1) > TMath::Abs(ii-i1)) continue;
  576. Int_t ip = ii, it = jj;
  577. if (intime) ip = jj, it = ii;
  578. //if (fCharges[ip][it] < fCharges[ipad][itime] / ratio) { valOk = 1; break; }
  579. if (fCharges[ip][it] > fCharges[ipad][itime] / ratio) { merge = 1; break; }
  580. }
  581. //if (valOk) break;
  582. if (!merge) break;
  583. }
  584. }
  585. //if (!valOk) fFlags[ipad][itime] = -fFlags[ipad][itime]; // not deep enough valley
  586. if (merge) fFlags[ipad][itime] = -fFlags[ipad][itime]; // not deep enough valley
  587. }
  588. }
  589. // Remove failed peaks
  590. multimap<Double_t,Int_t>::iterator it = localMax.begin();
  591. //for ( ; it != localMax.end(); ++it) cout << it->second << " ";
  592. //cout << endl;
  593. it = localMax.begin();
  594. for ( ; it != localMax.end(); ++it) {
  595. Int_t idig = it->second;
  596. Int_t ipad = clus->Col(idig);
  597. Int_t itime = clus->Bkt(idig);
  598. if (fFlags[ipad][itime] > 0) continue;
  599. //cout << " Before: " << idig << " " << itime << " " << ipad << " " << it->first << ", ";
  600. localMax.erase(it);
  601. //cout << " After: " << it->first << endl;
  602. }
  603. }
  604. //__________________________________________________________________________
  605. void MpdTpcClusterFinderAZ::CorrectReco(TVector3 &p3loc, TVector3 &p3errCor, Int_t nPads, Double_t &adc)
  606. {
  607. // Correct reconstructed coordinates and charge
  608. Double_t xsec = (p3loc.Y() + fSecGeo->GetMinY()) * TMath::Tan(fSecGeo->Dphi()/2);
  609. Double_t xloc = p3loc.X(), xloc0 = xloc, edge = 0.0; // distance to sector boundary
  610. if (xloc < 0) edge = xloc + xsec;
  611. else edge = xloc - xsec;
  612. if (TMath::Abs(edge) > 2.0) return; // no corrections
  613. Double_t adc0 = adc, coef = 1.0;
  614. // Correct charge
  615. if (nPads == 1 && adc0 < 5000) {
  616. coef = -0.0065 + 0.00015 * adc0 - 2.528e-8 * adc0 * adc0 + 1.380e-12 * adc0 * adc0 * adc0;
  617. coef = TMath::Max(coef,0.02) / 0.74 * 0.91;
  618. adc /= coef;
  619. } else if (nPads == 2 && adc0 < 10000) {
  620. coef = 0.0139 + 0.000201 * adc0 - 2.351e-8 * adc0 * adc0 + 9.750e-13 * adc0 * adc0 * adc0;
  621. coef = TMath::Max(coef,0.05) / 0.74 * 0.91;
  622. adc /= coef;
  623. } else if (nPads == 3 && adc0 < 6500) {
  624. coef = -0.068 + 0.000333 * adc0 - 4.401e-8 * adc0 * adc0 + 1.834e-12 * adc0 * adc0 * adc0;
  625. coef = TMath::Max(coef,0.02) / 0.74 * 0.91;
  626. adc /= coef;
  627. }
  628. if (TMath::Abs(edge) > 1.5 || nPads > 3) return; // no corrections
  629. // Correct coordinates
  630. if (nPads == 1) {
  631. xloc -= TMath::Sign(0.520,edge);
  632. Double_t corr = -0.15;
  633. if (adc0 < 50000) corr = 0.206 - 2.532e-5 * adc0 + 8.715e-10 * adc0 * adc0
  634. - 1.610e-14 * adc0 * adc0 * adc0 + 1.193e-19 * adc0 * adc0 * adc0 * adc0;
  635. if (edge < 0) corr = -corr;
  636. xloc -= corr;
  637. p3errCor[0] = 0.116;
  638. } else if (nPads == 2) {
  639. xloc -= TMath::Sign(0.200,edge);
  640. Double_t corr = -0.07;
  641. if (adc0 < 200000) corr = 0.374 - 1.046e-5 * adc0 + 9.595e-11 * adc0 * adc0
  642. - 4.098e-16 * adc0 * adc0 * adc0 + 6.811e-22 * adc0 * adc0 * adc0 * adc0;
  643. if (edge < 0) corr = -corr;
  644. xloc -= corr;
  645. p3errCor[0] = 0.110;
  646. } else if (nPads == 3) {
  647. xloc -= TMath::Sign(0.033,edge);
  648. Double_t corr = 0.01;
  649. if (adc0 < 200000) corr = 0.188 - 4.007e-6 * adc0 + 3.364e-11 * adc0 * adc0
  650. - 1.265e-16 * adc0 * adc0 * adc0 + 1.771e-22 * adc0 * adc0 * adc0 * adc0;
  651. if (edge < 0) corr = -corr;
  652. xloc -= corr;
  653. p3errCor[0] = 0.073;
  654. } else return;
  655. if (TMath::Abs(xloc) > xsec) xloc = TMath::Sign(xsec,xloc0);
  656. p3loc.SetX(xloc);
  657. }
  658. //__________________________________________________________________________
  659. ClassImp(MpdTpcClusterFinderAZ)