MpdTpcClusterFinderTask.cxx 50 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class TpcClusterFinderTask
  7. //
  8. // Environment:
  9. // Software developed for the MPD at NICA.
  10. //
  11. // Author List:
  12. // Sergey Merts
  13. // Artem Basalaev
  14. //
  15. //
  16. //-----------------------------------------------------------
  17. // Panda Headers ----------------------
  18. // This Class' Header ------------------
  19. #include "MpdTpcClusterFinderTask.h"
  20. #include "MpdTpcHit.h"
  21. #include "TpcCluster.h"
  22. #include "TpcPoint.h"
  23. #include "MpdTpcPeak.h"
  24. #include "FairRootManager.h"
  25. #include "MpdMCTrack.h"
  26. #include "TClonesArray.h"
  27. #include <TSystem.h>
  28. #include <iostream>
  29. #include "TaskHelpers.h"
  30. #include <TGraph.h>
  31. #include "MpdTpcSector.h"
  32. #include "TMath.h"
  33. #include "TVector3.h"
  34. #include "TFile.h"
  35. #include "sys/timeb.h"
  36. #include <sys/time.h>
  37. using namespace std;
  38. using namespace TMath;
  39. ClassImp(MpdTpcClusterFinderTask);
  40. static Float_t kOneOverSqrt12 = 1 / Sqrt(12);
  41. static Float_t kErrorCorrection = 0.115;
  42. static clock_t tStart = 0;
  43. static clock_t tFinish = 0;
  44. static clock_t tAll = 0;
  45. MpdTpcClusterFinderTask::MpdTpcClusterFinderTask() :
  46. fDigits(nullptr),
  47. fHitsArray(nullptr),
  48. fMCPointArray(nullptr),
  49. fMCTracks(nullptr),
  50. fDigitsArray(nullptr),
  51. fHisto(nullptr),
  52. fPersistence(kTRUE),
  53. fPrintDebugInfo(kFALSE),
  54. isHistogramsInitialized(kFALSE),
  55. fMakeQA(kFALSE),
  56. fCalcResiduals(kFALSE),
  57. fitGamma(kFALSE),
  58. fNumOfPadsInRow(nullptr) {
  59. inputBranchName = "MpdTpcDigit";
  60. outputBranchName = "MpdTpcHit";
  61. string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  62. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  63. fGas = new TpcGas(tpcGasFile, 130);
  64. }
  65. MpdTpcClusterFinderTask::~MpdTpcClusterFinderTask() {
  66. if (isHistogramsInitialized) delete fHisto;
  67. delete fHitsArray;
  68. delete fGas;
  69. }
  70. InitStatus MpdTpcClusterFinderTask::Init() {
  71. FairRootManager* ioman = FairRootManager::Instance();
  72. if (!ioman) {
  73. cout << "\n-E- [MpdTpcClusterFinderTask::Init]: RootManager not instantiated!" << endl;
  74. return kFATAL;
  75. }
  76. if (fCalcResiduals) {
  77. fMCPointArray = (TClonesArray*) ioman->GetObject("TpcPoint");
  78. fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
  79. }
  80. fDigits = (TClonesArray*) ioman->GetObject(inputBranchName);
  81. fHitsArray = new TClonesArray(outputBranchName);
  82. ioman->Register(outputBranchName, "TPC", fHitsArray, fPersistence);
  83. TpcSector* sector = new TpcSector();
  84. zDrift = sector->GetLength(); //cm
  85. nSect = sector->GetNSectors();
  86. nRows = sector->GetNumRows();
  87. fNumOfPadsInRow = (UInt_t*)sector->GetArrayPadsInRow();
  88. nTimeBins = sector->GetNTimeBins();
  89. pwIn = sector->GetInnerPadWidth();
  90. pwOut = sector->GetOuterPadWidth();
  91. phIn = sector->GetInnerPadHeight();
  92. phOut = sector->GetOuterPadHeight();
  93. nInRows = sector->GetNumInnerRows();
  94. nOutRows = sector->GetNumOuterRows();
  95. fSectInHeight = sector->GetSectInnerHeight();
  96. r_min = sector->GetRmin();
  97. delete sector;
  98. if (!isHistogramsInitialized && fMakeQA) {
  99. fHisto = new MpdTpcClusterFinderQAHistograms();
  100. fHisto->Initialize();
  101. isHistogramsInitialized = kTRUE;
  102. }
  103. fNoiseThreshold = 1000.0; //electrons
  104. fSpread = 0.196; // cm // Value is given by TPC group
  105. fEventId = 0;
  106. OnePeakCntr = 0;
  107. MoreThenOnePeakCntr = 0;
  108. OneDigitHitCntr = 0;
  109. AllCntr = 0;
  110. ThreePadsCntr = 0;
  111. //memory allocating for input array
  112. fDigitsArray = new DigOrigArray** [nRows];
  113. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  114. fDigitsArray[iRow] = new DigOrigArray* [fNumOfPadsInRow[iRow] * 2];
  115. for (UInt_t iPad = 0; iPad < fNumOfPadsInRow[iRow] * 2; ++iPad) {
  116. fDigitsArray[iRow][iPad] = new DigOrigArray [nTimeBins];
  117. for (UInt_t iTime = 0; iTime < nTimeBins; ++iTime) {
  118. fDigitsArray[iRow][iPad][iTime].isOverlap = kFALSE;
  119. fDigitsArray[iRow][iPad][iTime].origin = 0;
  120. fDigitsArray[iRow][iPad][iTime].origins.clear();
  121. fDigitsArray[iRow][iPad][iTime].signal = 0.0;
  122. }
  123. }
  124. }
  125. cout << "-I- MpdTpcClusterFinderTask: Initialization successful." << endl;
  126. return kSUCCESS;
  127. }
  128. void MpdTpcClusterFinderTask::Exec(Option_t* opt) {
  129. cout << "MpdTpcClusterFinder::Exec started! Event #" << ++fEventId << endl;
  130. tStart = clock();
  131. // Reset output Array
  132. if (!fHitsArray) Fatal("MpdTpcClusterFinder::Exec)", "No FoundClustersArray");
  133. fHitsArray->Delete();
  134. // parallel processing for different sectors
  135. #pragma omp parallel //num_threads(8)
  136. {
  137. #pragma omp for
  138. for (UInt_t iSec = 0; iSec < nSect; ++iSec) {
  139. for (Int_t i = 0; i < fDigits->GetEntriesFast(); ++i) {
  140. MpdTpcDigit* fDigit = (MpdTpcDigit*) fDigits->At(i);
  141. if ( (UInt_t) fDigit->GetSector() == iSec) {
  142. fDigitsArray[fDigit->GetRow()][fDigit->GetPad()][fDigit->GetTimeBin()].origin = fDigit->GetOrigin();
  143. fDigitsArray[fDigit->GetRow()][fDigit->GetPad()][fDigit->GetTimeBin()].signal = fDigit->GetAdc();
  144. }
  145. }
  146. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  147. vector<MpdTpcFoundHit*> tpcHitList; // output vector of TpcHits
  148. vector<MpdTpc2dCluster*> extClusters; // vector of extended clusters
  149. vector<MpdTpcPeak*> peakList; // vector of all peaks in cluster (for each cluster))
  150. Find2DClusters(&extClusters, fDigitsArray[iRow], iRow, iSec);
  151. for (UInt_t iClust = 0; iClust < extClusters.size(); ++iClust) {
  152. FindPeaksInCluster(extClusters.at(iClust), &peakList, fDigitsArray[iRow]);
  153. //mybeg
  154. //Peaks drawing
  155. if(fMakeQA) {
  156. Int_t numSect = fHisto->NumSector_hist;
  157. if (numSect == -1) {
  158. if(iRow == fHisto->NumRow_hist) {
  159. for(UInt_t iPeak = 0; iPeak < peakList.size(); iPeak++) {
  160. MpdTpcPeak *peak = peakList.at(iPeak);
  161. fHisto->_hXT_peak_row->Fill(peak->Col(), peak->MaxAdcBkt(), /*peak->SumADC()*/0.9);
  162. }
  163. }
  164. } else {
  165. if(iRow == fHisto->NumRow_hist && iSec == (UInt_t) numSect) {
  166. for(UInt_t iPeak = 0; iPeak < peakList.size(); iPeak++) {
  167. MpdTpcPeak *peak = peakList.at(iPeak);
  168. fHisto->_hXT_peak_row->Fill(peak->Col(), peak->MaxAdcBkt(), /*peak->SumADC()*/0.9);
  169. }
  170. }
  171. }
  172. }
  173. Float_t iSizeCell = 0.1; //for collection painting (cell of some coolection has its size)
  174. //myend
  175. vector<MpdTpcPeak*> collectedPeakList; //for each cluster
  176. if (peakList.empty()) {
  177. CreateHit(peakList, extClusters.at(iClust), &tpcHitList, fDigitsArray[iRow]);
  178. } else {
  179. while (!peakList.empty()) {
  180. CollectPeaks(peakList, extClusters.at(iClust), &collectedPeakList);
  181. //mybeg
  182. //Collected Peaks drawing
  183. if(fMakeQA) {
  184. if (fHisto->NumSector_hist == -1) {
  185. if(iRow == fHisto->NumRow_hist) {
  186. for(UInt_t iPeak = 0; iPeak < collectedPeakList.size(); iPeak++) {
  187. MpdTpcPeak *cpeak = collectedPeakList.at(iPeak);
  188. if(cpeak->Chi2() > 0.0) {
  189. fHisto->_hXT_collected_peak_row->Fill(cpeak->Col(), cpeak->PeakTime(), /*peak->SumADC()*/iSizeCell);
  190. }
  191. else {
  192. fHisto->_hXT_collected_peak_row->Fill(cpeak->Col(), cpeak->Mean()+cpeak->BktOff(), /*peak->SumADC()*/iSizeCell);
  193. }
  194. }
  195. }
  196. }
  197. else {
  198. if(iRow == fHisto->NumRow_hist && iSec == (UInt_t) fHisto->NumSector_hist) {
  199. for(UInt_t iPeak = 0; iPeak < collectedPeakList.size(); iPeak++) {
  200. MpdTpcPeak *cpeak = collectedPeakList.at(iPeak);
  201. if(cpeak->Chi2() > 0.0) {
  202. fHisto->_hXT_collected_peak_row->Fill(cpeak->Col(), cpeak->PeakTime(), /*peak->SumADC()*/iSizeCell);
  203. }
  204. else {
  205. fHisto->_hXT_collected_peak_row->Fill(cpeak->Col(), cpeak->Mean()+cpeak->BktOff(), /*peak->SumADC()*/iSizeCell);
  206. }
  207. }
  208. }
  209. }
  210. iSizeCell += 0.1; //increse in size for the next collection
  211. if(iSizeCell > 0.7) iSizeCell = 0.1;
  212. }
  213. //myend
  214. CreateHit(collectedPeakList, extClusters.at(iClust), &tpcHitList, fDigitsArray[iRow]);
  215. for (UInt_t i = 0; i < collectedPeakList.size(); ++i) {
  216. vector<MpdTpcPeak*>::iterator pkItr = find(peakList.begin(), peakList.end(), collectedPeakList[i]);
  217. delete *pkItr;
  218. peakList.erase(pkItr);
  219. }
  220. ///fixed bug with creating multiple identical clusters when cluster is oneHitCluster AND peakList.size()>1
  221. if (extClusters.at(iClust)->GetNumPads() == 1 || extClusters.at(iClust)->GetNumTimeBins() == 1) break;
  222. }
  223. }
  224. collectedPeakList.clear();
  225. peakList.clear();
  226. }
  227. for (UInt_t iHit = 0; iHit < tpcHitList.size(); ++iHit) {
  228. MpdTpcFoundHit* hit = tpcHitList.at(iHit);
  229. const UInt_t Sect = hit->GetSect();
  230. const Float_t sectPhi = Sect * TwoPi() / 12;
  231. const Float_t lX = hit->GetLocalX(); // local X coordinate
  232. const Float_t lY = hit->GetLocalY(); // local Y coordinate
  233. const Float_t lZ = hit->GetLocalZ(); // local Z coordinate
  234. const Float_t gX = (lY + r_min) * Cos(sectPhi) - lX * Sin(sectPhi); // global X coordinate
  235. const Float_t gY = (lY + r_min) * Sin(sectPhi) + lX * Cos(sectPhi); // global Y coordinate
  236. const Float_t gZ = (Sect < nSect / 2) ? lZ : lZ * (-1); // global Z coordinate
  237. hit->SetGlobalX(gX);
  238. hit->SetGlobalY(gY);
  239. hit->SetGlobalZ(gZ);
  240. const Float_t Q = hit->QADC(); //charge of hit
  241. if (fCalcResiduals) CalcResiduals(hit, sectPhi);
  242. Int_t row = hit->Cluster()->Row();
  243. UInt_t sec12 = (iSec < (nSect / 2)) ? iSec : (iSec - (nSect / 2));
  244. UInt_t padID = (sec12 | (row << 5));
  245. MpdTpcHit* outHit;
  246. #pragma omp critical
  247. {
  248. outHit = new((*fHitsArray)[fHitsArray->GetEntriesFast()]) MpdTpcHit(padID, TVector3(gX, gY, gZ), TVector3(hit->errX(), hit->errY(), hit->errZ()), hit->GetOrigin());
  249. }
  250. outHit->SetQ(Q);
  251. outHit->SetLocalXYZ(lX, lY, lZ);
  252. outHit->SetLayer(row);
  253. outHit->SetModular(1);
  254. if (fMakeQA) {
  255. #pragma omp critical
  256. {
  257. fHisto->_hXY->Fill(lX, lY, Q);
  258. fHisto->_hYZ_local->Fill(lZ, lY, Q);
  259. fHisto->_hRZ_global->Fill(gZ, gY, Q);
  260. fHisto->_hX->Fill(lX, Q);
  261. fHisto->_hY->Fill(lY, Q);
  262. fHisto->_hZ->Fill(lZ, Q);
  263. fHisto->_hX_global->Fill(gX, Q);
  264. fHisto->_hY_global->Fill(gY, Q);
  265. fHisto->_hZ_global->Fill(gZ, Q);
  266. fHisto->_hXY_global->Fill(gX, gY, Q);
  267. fHisto->_hPeak->Fill(Q);
  268. if ((UInt_t)hit->PadRow() < nInRows) {
  269. fHisto->_hErrX_inner->Fill(hit->errX());
  270. fHisto->_hErrZ_inner->Fill(hit->errZ());
  271. } else {
  272. fHisto->_hErrX_outer->Fill(hit->errX());
  273. fHisto->_hErrZ_outer->Fill(hit->errZ());
  274. }
  275. fHisto->_hErrY->Fill(hit->errY());
  276. fHisto->_h3D->Fill(gX, gY, gZ, Q);
  277. if (hit->GetSect() == 3 && lY <= 1.0) fHisto->_hXT_clust_row1->Fill(lX, lZ, Q); // bmy
  278. //mybeg
  279. //Clusters and hits drawing (together - on one histogram and separately - on different histograms)
  280. if (fHisto->NumSector_hist == -1) {
  281. if( (UInt_t) hit->Cluster()->Row() == fHisto->NumRow_hist) {
  282. for(UInt_t iDig = 0; iDig < hit->Cluster()->GetNumDigits(); ++iDig) {
  283. fHisto->_hXT_clust_hit_row->Fill(hit->Cluster()->Col(iDig),hit->Cluster()->Bkt(iDig), /*hit->Cluster()->Adc(iDig)*/0.001);
  284. fHisto->_hXT_clust_row->Fill(hit->Cluster()->Col(iDig),hit->Cluster()->Bkt(iDig), hit->Cluster()->Adc(iDig));
  285. }
  286. fHisto->_hXT_clust_hit_row->Fill(hit->PadCol(), hit->TimeBkt(), /*hit->QADC()*/1000.0);
  287. fHisto->_hXT_hit_row->Fill(hit->PadCol(), hit->TimeBkt(), hit->QADC());
  288. }
  289. } else {
  290. if( (UInt_t) hit->Cluster()->Row() == fHisto->NumRow_hist && hit->GetSect() == fHisto->NumSector_hist) {
  291. for(UInt_t iDig = 0; iDig < hit->Cluster()->GetNumDigits(); ++iDig) {
  292. fHisto->_hXT_clust_hit_row->Fill(hit->Cluster()->Col(iDig),hit->Cluster()->Bkt(iDig), /*hit->Cluster()->Adc(iDig)*/0.001);
  293. fHisto->_hXT_clust_row->Fill(hit->Cluster()->Col(iDig),hit->Cluster()->Bkt(iDig), hit->Cluster()->Adc(iDig));
  294. }
  295. fHisto->_hXT_clust_hit_row->Fill(hit->PadCol(), hit->TimeBkt(), /*hit->QADC()*/1000.0);
  296. fHisto->_hXT_hit_row->Fill(hit->PadCol(), hit->TimeBkt(), hit->QADC());
  297. }
  298. }
  299. //myend
  300. fHisto->_hNumOfDigitsInCluster->Fill(hit->Cluster()->GetNumDigits());
  301. fHisto->_hSect->Fill(hit->Cluster()->GetSect(), hit->Cluster()->GetADC());
  302. fHisto->_hNumOfPadsInCluster->Fill(hit->Cluster()->GetNumPads());
  303. fHisto->_hNumOfTimeBinsInCluster->Fill(hit->Cluster()->GetNumTimeBins());
  304. }
  305. }
  306. }
  307. if (fMakeQA && tpcHitList.size() != 0) fHisto->_hHitDistr->Fill(tpcHitList.size());
  308. for (UInt_t i = 0; i < extClusters.size(); ++i)
  309. delete extClusters.at(i);
  310. extClusters.clear();
  311. for (UInt_t i = 0; i < tpcHitList.size(); ++i)
  312. delete tpcHitList.at(i);
  313. tpcHitList.clear();
  314. for (UInt_t iPad = 0; iPad < fNumOfPadsInRow[iRow] * 2; ++iPad) {
  315. for (UInt_t iTime = 0; iTime < nTimeBins; ++iTime) {
  316. if (fDigitsArray[iRow][iPad][iTime].signal > 0.0) {
  317. fDigitsArray[iRow][iPad][iTime].origin = -1;
  318. fDigitsArray[iRow][iPad][iTime].origins.clear();
  319. fDigitsArray[iRow][iPad][iTime].signal = 0.0;
  320. fDigitsArray[iRow][iPad][iTime].isOverlap = kFALSE;
  321. }
  322. }
  323. }
  324. }
  325. }//for (UInt_t iSec = 0; iSec < nSect; ++iSec)
  326. }
  327. tFinish = clock();
  328. tAll = tAll + (tFinish - tStart);
  329. if (fPrintDebugInfo) cout << "Number of Found Clusters = " << fHitsArray->GetEntriesFast() << endl;
  330. cout << "MpdTpcClusterFinder::Exec finished" << endl;
  331. }
  332. Bool_t MpdTpcClusterFinderTask::Find2DClusters(vector<MpdTpc2dCluster*> *extClusters, DigOrigArray **f2dArray, UInt_t row, UInt_t sec) {
  333. UInt_t curDigit[3];
  334. Bool_t result = kFALSE;
  335. UInt_t cluscount = 0;
  336. Bool_t** fADCMarks = new Bool_t* [fNumOfPadsInRow[row] * 2];
  337. for (UInt_t i = 0; i < fNumOfPadsInRow[row] * 2; ++i) {
  338. fADCMarks[i] = new Bool_t[nTimeBins];
  339. for (UInt_t j = 0; j < nTimeBins; ++j)
  340. fADCMarks[i][j] = kFALSE;
  341. }
  342. for (UInt_t pad = 0; pad < fNumOfPadsInRow[row] * 2; ++pad) {
  343. for (UInt_t tBin = 0; tBin < nTimeBins; ++tBin) {
  344. if (f2dArray[pad][tBin].signal < fNoiseThreshold) continue;
  345. curDigit[0] = row;
  346. curDigit[1] = pad;
  347. curDigit[2] = tBin;
  348. if (!fADCMarks[pad][tBin]) { //if current ADC has no mark yet
  349. MpdTpc2dCluster* cluster = new MpdTpc2dCluster(row, sec);
  350. fADCMarks[pad][tBin] = kTRUE;
  351. cluster->SetID(cluscount++);
  352. if (!cluster->Insert(row, pad, tBin, f2dArray[pad][tBin].signal)) {
  353. #pragma omp critical
  354. {
  355. cout << "Failed to insert digit into cluster..." << endl;
  356. }
  357. return kFALSE;
  358. }
  359. // now that we have the first digit in the 2d cluster, the following
  360. // will find the rest of the connected digits
  361. result = GetNextDigit(curDigit, cluster, f2dArray, fADCMarks);
  362. if (!result) return kFALSE;
  363. extClusters->push_back(cluster);
  364. }
  365. }
  366. }
  367. for (UInt_t iPad = 0; iPad < fNumOfPadsInRow[row] * 2; ++iPad)
  368. delete [] fADCMarks[iPad];
  369. delete [] fADCMarks;
  370. return kTRUE;
  371. }
  372. Bool_t MpdTpcClusterFinderTask::GetNextDigit(UInt_t* currdig, MpdTpc2dCluster* Clus2d, DigOrigArray **fDigitsArray1, Bool_t **fADCMarks) {
  373. UInt_t thisCol, thisRow, thisBkt;
  374. UInt_t nextCol, nextBkt;
  375. thisRow = currdig[0];
  376. thisCol = currdig[1];
  377. thisBkt = currdig[2];
  378. for (int i = -1; i < 2; ++i) {
  379. // check left-right bounds...
  380. if (thisCol == 0 && i < 0) continue;
  381. if (thisCol == fNumOfPadsInRow[thisRow] * 2 - 1 && i >= 0) continue;
  382. nextCol = thisCol + i; // look right, left, check if we're at the edge
  383. for (int j = -1; j < 2; ++j) {
  384. if (i == 0 && j == 0) continue;
  385. // check up-down bounds...
  386. if (thisBkt == 0 && j < 0) continue;
  387. if (thisBkt == nTimeBins - 1 && j >= 0) continue;
  388. nextBkt = thisBkt + j; // look down, up, check if we're at the edge
  389. if (fDigitsArray1[nextCol][nextBkt].signal < fNoiseThreshold) continue;
  390. UInt_t nextdig[3] = {thisRow, nextCol, nextBkt};
  391. if (!fADCMarks[nextCol][nextBkt]) {
  392. if (!Clus2d->Insert(thisRow, nextCol, nextBkt, fDigitsArray1[nextCol][nextBkt].signal)) {
  393. #pragma omp critical
  394. {
  395. cout << "Failed to insert digit into cluster..." << endl;
  396. }
  397. return kFALSE;
  398. }
  399. fADCMarks[nextCol][nextBkt] = kTRUE;
  400. GetNextDigit(nextdig, Clus2d, fDigitsArray1, fADCMarks);
  401. }
  402. }
  403. }
  404. return kTRUE;
  405. }
  406. void MpdTpcClusterFinderTask::FindPeaksInCluster(MpdTpc2dCluster* clust, vector<MpdTpcPeak*> *PeakList, DigOrigArray** fOriginsArray) {
  407. const UInt_t fPeakValleyRatio = 2; // maximum of peak / minimum of peak
  408. const Int_t mincol = clust->MinCol();
  409. const Int_t minbkt = clust->MinBkt();
  410. const Int_t nbkt = clust->GetNumTimeBins();
  411. const Int_t ncol = clust->GetNumPads();
  412. const Int_t ndig = clust->GetNumDigits();
  413. Float_t xyADC[ncol][nbkt]; //box around 2D cluster
  414. for (Int_t i = 0; i < ncol; ++i)
  415. for (Int_t j = 0; j < nbkt; ++j)
  416. xyADC[i][j] = 0.0;
  417. for (Int_t i = 0; i < ndig; ++i) // data extraction block.
  418. xyADC[clust->Col(i) - mincol][clust->Bkt(i) - minbkt] = clust->Adc(i);
  419. // now loop through each pad and extrapolate individual peaks
  420. /// TODO!!! Set proper values for our TPC
  421. ///minthresh is a minimum ADC required for peak to be created
  422. const Float_t adcmin = 1000.0;
  423. const Float_t minthresh = adcmin + 2 * sqrt(adcmin); //???
  424. for (Int_t i = 0; i < ncol; ++i) {
  425. Float_t thresh = 0.;
  426. Float_t max = 0.;
  427. Float_t min = 1.e9;
  428. //Int_t imax = -1;
  429. //Int_t imin = -1;
  430. Int_t npeak = 0;
  431. MpdTpcPeak* peak = new MpdTpcPeak();
  432. peak->SetCol(i + mincol);
  433. peak->AttachCluster(clust);
  434. Bool_t FoundPeak = kFALSE;
  435. for (Int_t j = 0; j < nbkt; ++j) {
  436. Float_t adc = xyADC[i][j];
  437. if (adc > minthresh) {
  438. if (peak->NSamples() == 0) peak->SetBktOff(minbkt + j);
  439. peak->Insert(adc, i + mincol, j + minbkt);
  440. if (adc > max) {
  441. max = adc;
  442. //imax = i;
  443. thresh = max - 2 * sqrt(max); // ???
  444. }
  445. }
  446. if (FoundPeak && adc < min) {
  447. min = adc;
  448. //imin = i;
  449. }
  450. if (adc < thresh) FoundPeak = kTRUE;
  451. if ((FoundPeak && adc < minthresh) ||
  452. // ^-- here we're below our min. threshold, so we save
  453. // the peak info. and create a new peak
  454. (FoundPeak && adc > min)) {
  455. // ^-- here we've found a peak but the slope is rising... could be
  456. // entering a new peak, so save the peak info. and create a new one
  457. if (adc > min) // remove current ADC hit from peak
  458. peak->Remove(peak->NSamples() - 1);
  459. FoundPeak = kFALSE;
  460. if (peak->NSamples() > 1 && (peak->Max() / peak->Min() > fPeakValleyRatio)) {
  461. peak->SetOrigin(fOriginsArray[peak->MaxAdcCol()][peak->MaxAdcBkt()].origin);
  462. PeakList->push_back(peak);
  463. ++npeak;
  464. peak = new MpdTpcPeak();
  465. } else {
  466. peak->Clear();
  467. }
  468. peak->SetCol(i + mincol);
  469. peak->SetBktOff(minbkt + j);
  470. peak->AttachCluster(clust);
  471. if (adc > minthresh) peak->Insert(adc, i + mincol, j + minbkt);
  472. thresh = 0.;
  473. max = 0.;
  474. min = 1.e9;
  475. //imax = -1;
  476. //imin = -1;
  477. } // end check for 2nd peak (or clear end of 1st peak)
  478. } // end loop through buckets in this column
  479. // need to add check that _last_ peak was not already added to list!
  480. if (FoundPeak && peak->NSamples() > 1 && (peak->Max() / peak->Min() > fPeakValleyRatio)) {
  481. // check that we haven't already added this peak, just to be safe...
  482. vector<MpdTpcPeak*>::iterator pkItr = find(PeakList->begin(), PeakList->end(), peak);
  483. if (pkItr == PeakList->end()) {
  484. peak->SetOrigin(fOriginsArray[peak->MaxAdcCol()][peak->MaxAdcBkt()].origin);
  485. PeakList->push_back(peak);
  486. ++npeak;
  487. }
  488. }
  489. if (npeak == 0) { // no peaks were found, so we'll store the ADC info
  490. // as a single "peak" anyway...
  491. peak->Clear();
  492. peak->SetCol(i + mincol);
  493. peak->AttachCluster(clust);
  494. for (Int_t j = 0; j < nbkt; ++j) {
  495. if (xyADC[i][j] > minthresh) {
  496. if (peak->NSamples() == 0) peak->SetBktOff(minbkt + j);
  497. peak->Insert(xyADC[i][j], i + mincol, j + minbkt);
  498. }
  499. }
  500. if (peak->NSamples() > 1) {
  501. peak->SetOrigin(fOriginsArray[peak->MaxAdcCol()][peak->MaxAdcBkt()].origin);
  502. PeakList->push_back(peak);
  503. }
  504. }
  505. // finally, check that peak has not been added to the peak list...
  506. // if not, then delete it!
  507. vector<MpdTpcPeak*>::iterator pkItr = find(PeakList->begin(), PeakList->end(), peak);
  508. if (pkItr == PeakList->end()) delete peak;
  509. } // end loop thru pads in this column
  510. }
  511. void MpdTpcClusterFinderTask::CollectPeaks(vector<MpdTpcPeak*> peakList, MpdTpc2dCluster* clust, vector<MpdTpcPeak*>* collectedPeakList) {
  512. const Int_t mincol = clust->MinCol();
  513. const Int_t ncol = clust->GetNumPads();
  514. Int_t npeaks[ncol];
  515. for (Int_t i = 0; i < ncol; ++i) npeaks[i] = 0;
  516. if (fitGamma)
  517. FitPeaks(peakList, npeaks, mincol);
  518. else {
  519. for (UInt_t i = 0; i < peakList.size(); ++i) {
  520. MpdTpcPeak* peak = peakList[i];
  521. peak->SetChi2(0.0); // TMP
  522. }
  523. }
  524. collectedPeakList->clear();
  525. collectedPeakList->push_back(peakList[0]);
  526. // collect all peaks that line-up with each other
  527. Float_t tprev, tnext, tdiff; //time of current and previous peaks and time shift between two peaks
  528. Int_t cprev, cnext, cdiff; //pad of current and previous peaks and pad shift between two peaks
  529. if (peakList[0]->Chi2() > 0.)
  530. tprev = peakList[0]->PeakTime();
  531. else {
  532. //Float_t x1 = peakList[0]->NSamples();
  533. //Float_t dt1 = par[0]*(1 - exp(-(par[1] - x1)*(par[1] - x1) / par[2]));
  534. tprev = peakList[0]->Mean() + peakList[0]->BktOff(); // - dt1;
  535. }
  536. cprev = peakList[0]->Col();
  537. for (UInt_t i = 1; i < peakList.size(); ++i) {
  538. if (peakList[i]->Chi2() > 0.)
  539. tnext = peakList[i]->PeakTime();
  540. else {
  541. //Float_t x1 = peakList[i]->NSamples();
  542. //Float_t dt1 = par[0]*(1 - exp(-(par[1] - x1)*(par[1] - x1) / par[2]));
  543. tnext = peakList[i]->Mean() + peakList[i]->BktOff(); // - dt1;
  544. }
  545. tdiff = tnext - tprev;
  546. cnext = peakList[i]->Col();
  547. cdiff = cnext - cprev;
  548. if (fabs(tdiff) < 2. && abs(cdiff) == 1) {
  549. collectedPeakList->push_back(peakList[i]);
  550. cprev = cnext;
  551. }
  552. }
  553. }
  554. void MpdTpcClusterFinderTask::CreateHit(vector<MpdTpcPeak*> collectedPeakList, MpdTpc2dCluster* clust, vector<MpdTpcFoundHit*> *hitList, DigOrigArray** fOriginsArray) {
  555. AllCntr++;
  556. if (clust->GetNumPads() == 3) ThreePadsCntr++;
  557. Bool_t oneHitCluster = kFALSE;
  558. if (clust->GetNumPads() == 1 || clust->GetNumTimeBins() == 1) oneHitCluster = kTRUE;
  559. if (!oneHitCluster && !collectedPeakList.empty()) {
  560. // TF1* fitfcn = nullptr;
  561. Float_t sigmaX = -1.0;
  562. MpdTpcFoundHit* hit = new MpdTpcFoundHit(clust);
  563. // now calculate <t>, sig<t>, <x> and sig<x>:
  564. Float_t avgt = 0.0, avgx = 0.0;
  565. Float_t sumq = 0.0;
  566. Float_t sigt = 0.0, sigx = 0.0;
  567. Float_t w = 0.0, dw = 0.;
  568. Bool_t GoodHit = kFALSE;
  569. if (collectedPeakList.size() == 1) { // handle the case where a single _clear_
  570. // peak is found, but no neighboring peaks ...
  571. OnePeakCntr++;
  572. MpdTpcPeak* peak = collectedPeakList[0];
  573. Float_t pTime = peak->PeakTime();
  574. Float_t pInt = peak->Integral();
  575. Int_t pCol = peak->Col();
  576. //if (peak->Chi2() > 0. && peak->Chi2() < 30. && peak->Max() > 100.) {
  577. hit->SetQFit(pInt);
  578. hit->SetSigQFit(peak->IntegSig());
  579. hit->SetType(MpdTpcFoundHit::kFitPeak);
  580. // now estimate sigx and sigt...
  581. // draw a 3x3 box around peak
  582. Float_t xyADC[3][5];
  583. Float_t xADC[3];
  584. memset(xyADC, 0, sizeof (xyADC));
  585. memset(xADC, 0, sizeof (xADC));
  586. int ncol = 1;
  587. for (int ih = 0; ih < clust->NDigits(); ++ih) {
  588. int ic = (Int_t) (clust->Col(ih) - pCol);
  589. int ib = (Int_t) (clust->Bkt(ih) - pTime);
  590. if (abs(ic) <= 1 && abs(ib) <= 2) {
  591. xyADC[ic + 1][ib + 2] = clust->Adc(ih);
  592. xADC[ic + 1] += xyADC[ic + 1][ib + 2];
  593. }
  594. }
  595. if (xADC[0] > 0) ++ncol;
  596. if (xADC[2] > 0) ++ncol;
  597. if (ncol > 1) {
  598. // first deal with x-position
  599. avgx = 0.33 * xADC[0] * (pCol - 1) + 0.33 * xADC[2] * (pCol + 1) + pInt * pCol;
  600. avgx /= (0.33 * (xADC[0] + xADC[2]) + pInt);
  601. sigx = 0.33 * xADC[0] * (pCol - 1 - avgx) * (pCol - 1 - avgx) + 0.33 * xADC[2] * (pCol + 1 - avgx) * (pCol + 1 - avgx) + pInt * (pCol - avgx) * (pCol - avgx);
  602. sigx /= (0.33 * xADC[0] + pInt + 0.33 * xADC[2]) * (ncol - 1);
  603. if (sigx < 0.0001) // this is completely unrealistic, so we assume something really bad happened here...
  604. sigx = kOneOverSqrt12;
  605. else
  606. sigx = sqrt(sigx) * kErrorCorrection;
  607. // now deal with y-position
  608. Float_t y1 = 0;
  609. for (int ib = -2; ib <= 2; ++ib)
  610. y1 += (Int_t) (pTime + ib) * xyADC[0][ib + 2];
  611. if (xADC[0] > 0) y1 /= xADC[0];
  612. Float_t y2 = 0;
  613. for (int ib = -2; ib <= 2; ++ib)
  614. y2 += (Int_t) (pTime + ib) * xyADC[2][ib + 2];
  615. if (xADC[2] > 0) y2 /= xADC[2];
  616. avgt = 0.33 * xADC[0] * y1 + 0.33 * xADC[2] * y2 + pTime * pInt;
  617. avgt /= (0.33 * (xADC[0] + xADC[2]) + pInt);
  618. if (fabs(pTime - avgt) < 1.) {
  619. // looks reasonable, so continue
  620. sigt = 0.33 * xADC[0] * (y1 - avgt) * (y1 - avgt) + 0.33 * xADC[2] * (y2 - avgt) * (y2 - avgt) + pInt * (pTime - avgt) * (pTime - avgt);
  621. sigt /= (0.33 * xADC[0] + pInt + 0.33 * xADC[2]) * (ncol - 1);
  622. if (sigt < 0.001) sigt = peak->SigMean() / peak->SumADC();
  623. else sigt = sqrt(sigt);
  624. GoodHit = kTRUE;
  625. }
  626. }
  627. //}
  628. } else { // more than one peak found in cluster...
  629. MoreThenOnePeakCntr++;
  630. GoodHit = kTRUE;
  631. // TH1F* h = new TH1F("h", "h", 50, 0, 50);
  632. for (UInt_t i = 0; i < collectedPeakList.size(); ++i) {
  633. MpdTpcPeak* peak = collectedPeakList[i];
  634. if (peak->Chi2() > 0.) { // this means the peak was successfully fit.
  635. w = peak->Integral();
  636. dw = peak->IntegSig();
  637. hit->SetQFit(hit->QFit() + w);
  638. hit->SetSigQFit(hit->QFitSig() + dw * dw);
  639. hit->SetType(hit->Type() | MpdTpcFoundHit::kFitPeak);
  640. Float_t myt = peak->PeakTime();
  641. /// Correction for TPC MPD may be added here
  642. // if (fUseTCorr) myt -= fTCorrTable->TCorr(peak->Row(),peak->Col());
  643. avgt += myt*w;
  644. sigt += myt * myt*w;
  645. } else { // fit was not successfully fit
  646. w = peak->SumADC();
  647. hit->SetQADC(hit->QADC() + peak->SumADC());
  648. hit->SetType(hit->Type() | MpdTpcFoundHit::kWMPeak);
  649. hit->AddOrigin(peak->GetOrigin());
  650. //Float_t x1 = peak->NSamples();
  651. //Float_t dt1 = par[0]*(1 - exp(-(par[1] - x1)*(par[1] - x1) / par[2]));
  652. Float_t myt = peak->Mean() + peak->BktOff(); // - dt1;
  653. /// Correction for TPC MPD may be added here
  654. // if (fUseTCorr) myt -= fTCorrTable->TCorr(peak->Row(),peak->Col());
  655. avgt += myt * w;
  656. sigt += myt * myt * w;
  657. }
  658. // h->Fill(peak->Col(), w);
  659. avgx += peak->Col() * w;
  660. sigx += peak->Col() * peak->Col() * w;
  661. sumq += w;
  662. }
  663. // if (h->GetEntries() > 3) {
  664. // h->Fit("gaus", "wQ");
  665. // fitfcn = h->GetFunction("gaus");
  666. // if (fitfcn) sigmaX = fitfcn->GetParameter(2);
  667. // }
  668. // delete h;
  669. hit->SetSigQFit(Sqrt(hit->QFitSig()));
  670. avgt /= sumq;
  671. avgx /= sumq;
  672. sigt = (sigt / sumq - avgt * avgt);
  673. sigx = (sigx / sumq - avgx * avgx);
  674. if (sigt <= 0.0) sigt = kOneOverSqrt12; //(hitPeaks.size()) * kOneOverSqrt12;
  675. else sigt = Sqrt(sigt);
  676. if (sigx <= 0.0) sigx = kOneOverSqrt12; //sigx = (hitPeaks.size()) * kOneOverSqrt12;
  677. else sigx = Sqrt(sigx) * kErrorCorrection;
  678. }
  679. Float_t pos[3], dpos[3];
  680. if (GoodHit) {
  681. Float_t padW, padH;
  682. UInt_t Row = (UInt_t) clust->Row();
  683. if (Row < nInRows) {
  684. padW = pwIn;
  685. padH = phIn;
  686. pos[1] = padH * ((Float_t) Row + 0.5); // y-coordinate of pad center
  687. } else {
  688. padW = pwOut;
  689. padH = phOut;
  690. pos[1] = fSectInHeight + ((Float_t) (Row - nInRows) + 0.5) * padH; // y-coordinate of pad center
  691. }
  692. pos[0] = padW * ((Float_t) avgx - (Float_t) fNumOfPadsInRow[Row] + 0.5);
  693. pos[2] = ((Float_t) avgt + 0.5) * (zDrift / nTimeBins);
  694. hit->SetPadCol(avgx);
  695. hit->SetTimeBkt(avgt);
  696. //adding errors
  697. const Float_t timeBktSize = zDrift / nTimeBins;
  698. //cout << sigx * padW << " " << sigmaX * padW << endl;
  699. if (sigmaX > 0.0) dpos[0] = sigmaX * padW; //dx
  700. else dpos[0] = sigx * padW; //dx
  701. dpos[1] = padH * kOneOverSqrt12; //dy
  702. dpos[2] = sigt * timeBktSize; //dz
  703. hit->SetPos(pos, dpos);
  704. hitList->push_back(hit);
  705. } else
  706. delete hit;
  707. } else { //oneHitCluster is true or collectedPeakList is empty
  708. Int_t iPad, iRow, iTime;
  709. Float_t y;
  710. // Initialize variables to 0
  711. Float_t avgx = 0.0, avgz = 0.0;
  712. Float_t rmsx = 0.0, rmsz = 0.0;
  713. Float_t adc = 0.0, sumadc = 0.0, invsum = 0.0;
  714. UInt_t ndig = (UInt_t) clust->GetNumDigits();
  715. Int_t Origin{0};
  716. for (UInt_t i = 0; i < ndig; ++i) {
  717. iPad = clust->Col(i);
  718. iTime = clust->Bkt(i);
  719. iRow = clust->Row(i);
  720. adc = clust->Adc(i);
  721. Int_t CurrDigOrigin = fOriginsArray[iPad][iTime].origin;
  722. Origin = CurrDigOrigin;
  723. // y-coordinate of pad center
  724. y = ((UInt_t)iRow < nInRows) ? (phIn * ((Float_t) iRow + 0.5)) : (fSectInHeight + ((Float_t) (iRow - nInRows) + 0.5) * phOut);
  725. // Compute weighted sums of x, t, x^2, and t^2
  726. avgx += adc * iPad;
  727. avgz += adc * iTime;
  728. rmsx += adc * iPad * iPad;
  729. rmsz += adc * iTime * iTime;
  730. sumadc += adc;
  731. }
  732. // Make sure that we had at least one digit with charge
  733. if (sumadc > 0.0) {
  734. invsum = 1.0 / sumadc;
  735. avgx *= invsum;
  736. avgz *= invsum;
  737. rmsx = rmsx * invsum - avgx * avgx;
  738. rmsz = rmsz * invsum - avgz * avgz;
  739. }
  740. if (clust->GetNumPads() == 1 || rmsx <= 0.0) {
  741. Float_t xErrOnePad = (y < nInRows * phIn) ? (pwIn * kOneOverSqrt12) : (pwOut * kOneOverSqrt12);
  742. clust->SetErrX(xErrOnePad);
  743. } else {
  744. Float_t xErr = (y < nInRows * phIn) ? (pwIn * Sqrt(rmsx)) : (pwOut * Sqrt(rmsx));
  745. clust->SetErrX(xErr);
  746. }
  747. if (rmsz <= 0.0) clust->SetErrZ(zDrift / nTimeBins * kOneOverSqrt12);
  748. else clust->SetErrZ(Sqrt(rmsz) * zDrift / nTimeBins);
  749. ///mistake fixed!
  750. //if (clust->GetSect() > nSect / 2) avgz *= -1;
  751. clust->SetX(avgx);
  752. clust->SetY(y);
  753. clust->SetZ(avgz);
  754. Float_t yErrOneRow = (y < nInRows * phIn) ? (phIn * kOneOverSqrt12) : (phOut * kOneOverSqrt12);
  755. clust->SetErrY(yErrOneRow);
  756. clust->SetADC(sumadc);
  757. Float_t pos[3];
  758. Float_t hitErr[3];
  759. const Float_t padW = ( (UInt_t) clust->Row() < nInRows) ? pwIn : pwOut;
  760. pos[0] = padW * ((Float_t) clust->GetX() - (Float_t) fNumOfPadsInRow[clust->Row()] + 0.5);
  761. pos[1] = y;
  762. pos[2] = ((Float_t) clust->GetZ()) * (zDrift / nTimeBins);
  763. //adding errors
  764. hitErr[0] = clust->GetErrX(); //dx
  765. hitErr[1] = clust->GetErrY(); //dy
  766. hitErr[2] = clust->GetErrZ(); //dz
  767. MpdTpcFoundHit* tmphit = new MpdTpcFoundHit(pos, hitErr);
  768. tmphit->AttachCluster(clust);
  769. tmphit->SetTimeBkt(clust->GetZ());
  770. tmphit->SetPadCol(clust->GetX());
  771. tmphit->SetNumHits(1);
  772. tmphit->SetType(MpdTpcFoundHit::kWMPeak);
  773. tmphit->SetQADC(clust->GetADC());
  774. tmphit->AddOrigin(Origin);
  775. hitList->push_back(tmphit);
  776. OneDigitHitCntr++;
  777. }
  778. }
  779. void MpdTpcClusterFinderTask::FitPeaks(vector<MpdTpcPeak*> PeakList, Int_t npeaks[], Int_t mincol) {
  780. Float_t par[3], epar[3], chi2;
  781. Float_t lnq[100], w[100], dt[100], lndt[100];
  782. for (UInt_t i = 0; i < PeakList.size(); ++i) {
  783. MpdTpcPeak* peak = PeakList[i];
  784. Int_t ny = peak->NSamples();
  785. // skip the fitting if this peak has too few bins...
  786. //Int_t fMinNBktGamma = 5; ///because this is the minimum number of bins in hit, must be re-checked!
  787. if (ny < 5) continue;
  788. Float_t yADC[ny];
  789. npeaks[peak->Col() - mincol] += 1;
  790. for (Int_t j = 0; j < ny; ++j) {
  791. yADC[j] = peak->Sample(j);
  792. lnq[j] = log(yADC[j]);
  793. w[j] = yADC[j] * yADC[j];
  794. }
  795. // search for good guess for t0
  796. Float_t t0, t0_save;
  797. t0_save = -100.;
  798. Float_t minchi2 = 1.e30;
  799. for (Float_t k = -2.; k < 2.1; k += 0.5) {
  800. t0 = peak->BktOff() + k;
  801. for (Int_t j = 0; j < ny; ++j) {
  802. dt[j] = peak->BktOff() + (Float_t) j - t0;
  803. if (dt[j] < 1.e-3) dt[j] = 1.e-10;
  804. lndt[j] = Log(dt[j]);
  805. }
  806. LinGammaFit(lnq, w, dt, lndt, ny, par, epar, chi2);
  807. if (fPrintDebugInfo) {
  808. #pragma omp critical
  809. {
  810. cout << "t0 = " << t0 << endl;
  811. for (int j = 0; j < 3; ++j) cout << "\tpar[" << j << "] = " << par[j] << endl;
  812. cout << "\t chi2 = " << chi2 << endl;
  813. }
  814. }
  815. if (chi2 > 0 && par[0] == par[0] && par[1] == par[1] && par[2] == par[2] && chi2 < minchi2) {
  816. minchi2 = chi2;
  817. t0_save = t0;
  818. }
  819. }
  820. if (fPrintDebugInfo)
  821. #pragma omp critical
  822. {
  823. cout << "t0_save = " << t0_save << endl;
  824. }
  825. // check that we found a good guess for t0; if not, skip this peak
  826. if (t0_save < peak->BktOff() - 2.) continue;
  827. Float_t t0_start = t0_save - 0.5;
  828. Float_t t0_stop = t0_save + 0.6;
  829. t0_save = -100.;
  830. minchi2 = 1.e30;
  831. Float_t par_save[3] = {0., 0., 0.};
  832. Float_t epar_save[3] = {0., 0., 0.};
  833. Float_t chi2_save = 0.;
  834. // now do a finer-binned search for the best t0...
  835. for (t0 = t0_start; t0 < t0_stop; t0 += 0.1) {
  836. for (Int_t j = 0; j < ny; ++j) {
  837. dt[j] = peak->BktOff() + (Float_t) j - t0;
  838. if (dt[j] < 1.e-3) dt[j] = 1.e-10;
  839. lndt[j] = log(dt[j]);
  840. }
  841. LinGammaFit(lnq, w, dt, lndt, ny, par, epar, chi2);
  842. if (fPrintDebugInfo) {
  843. #pragma omp critical
  844. {
  845. cout << "t0 = " << t0 << endl;
  846. for (Int_t j = 0; j < 3; ++j) cout << "\tpar[" << j << "] = " << par[j] << endl;
  847. cout << "\t chi2 = " << chi2 << endl;
  848. }
  849. }
  850. if (chi2 > 0 && par[0] == par[0] && par[1] == par[1] && par[2] == par[2] && chi2 < minchi2) {
  851. minchi2 = chi2;
  852. t0_save = t0;
  853. for (int m = 0; m < 3; ++m) {
  854. par_save[m] = par[m];
  855. epar_save[m] = epar[m];
  856. chi2_save = chi2;
  857. }
  858. }
  859. }
  860. if (fPrintDebugInfo) {
  861. #pragma omp critical
  862. {
  863. cout << "Final fit result:" << endl;
  864. cout << "\t t0 = " << t0_save << endl;
  865. for (int m = 0; m < 3; ++m) cout << "\t par[" << m << "] = " << par_save[m] << " +- " << epar_save[m] << endl;
  866. cout << "\t chi^2/ndf = " << chi2_save << endl;
  867. }
  868. }
  869. // now estimate the uncertainty on the integrated charge
  870. Float_t gam1 = Gamma(par_save[2] + 1);
  871. Float_t dgam1 = (Gamma(par_save[2] + 1 + 1.e-5) - Gamma(par_save[2] + 1 - 1.e-5)) / 2.e-5;
  872. Float_t Q = par_save[0] * par_save[1] * gam1;
  873. Float_t dQ = epar_save[0] * epar_save[0] / (par_save[0] * par_save[0]) + epar_save[1] * epar_save[1] / (par_save[1] * par_save[1]) + dgam1 * dgam1 * epar_save[2] * epar_save[2] / (gam1 * gam1);
  874. dQ = sqrt(fabs(dQ));
  875. dQ *= Q;
  876. peak->SetPeakTime(t0_save + par_save[1] * par_save[2]);
  877. // peak->SetPeakTime(WeightedAverage(peak));
  878. peak->SetIntegral(Q);
  879. peak->SetIntegSig(dQ);
  880. peak->SetChi2(chi2_save);
  881. }
  882. }
  883. //.......................................................................
  884. // FitCluster finds "peaks" in each pad of a cluster, then tries to
  885. // fit these peaks to a Gamma distribution. We then collect all peaks
  886. // that are within 2 time buckets of each other to form a hit. The hit
  887. // time position is taken as the weighted mean of the peak positions, and
  888. // the x-position is taken as the weighted pad number, where the weights
  889. // in both cases are the total charge. For peaks that were fit to Gamma,
  890. // the charge is the integral of the Gamma Dist. For peaks that were not
  891. // fit to Gamma, the charge is set to the summed ADC values.
  892. //.......................................................................
  893. void MpdTpcClusterFinderTask::LinGammaFit(Float_t lnq[], Float_t w[], Float_t dt[], Float_t lndt[], Int_t nbin, Float_t par[], Float_t epar[], Float_t& chi2) {
  894. Float_t a[3], b[3];
  895. Float_t C[3][3]; // C is symmetric
  896. Float_t CInv[3][3];
  897. for (Int_t i = 0; i < 3; ++i) {
  898. a[i] = b[i] = 0.;
  899. for (Int_t j = 0; j < 3; ++j) C[i][j] = CInv[i][j] = 0.;
  900. }
  901. Float_t sumlnq2 = 0;
  902. for (int i = 0; i < nbin; ++i) {
  903. sumlnq2 += w[i] * lnq[i] * lnq[i];
  904. b[0] += w[i] * lnq[i];
  905. b[1] += w[i] * dt[i] * lnq[i];
  906. b[2] += w[i] * lndt[i] * lnq[i];
  907. C[0][0] += w[i];
  908. C[0][1] += w[i] * dt[i];
  909. C[0][2] += w[i] * lndt[i];
  910. C[1][1] += w[i] * dt[i] * dt[i];
  911. C[1][2] += w[i] * dt[i] * lndt[i];
  912. C[2][2] += w[i] * lndt[i] * lndt[i];
  913. }
  914. C[1][0] = C[0][1];
  915. C[2][0] = C[0][2];
  916. C[2][1] = C[1][2];
  917. Float_t Cdet = C[0][0] * C[1][1] * C[2][2] + C[0][1] * C[1][2] * C[2][0]
  918. + C[0][2] * C[1][0] * C[2][1] - C[0][2] * C[1][1] * C[2][0]
  919. - C[1][2] * C[2][1] * C[0][0] - C[2][2] * C[0][1] * C[1][0];
  920. CInv[0][0] = (C[1][1] * C[2][2] - C[1][2] * C[2][1]) / Cdet;
  921. CInv[0][1] = (C[0][2] * C[2][1] - C[0][1] * C[2][2]) / Cdet;
  922. CInv[0][2] = (C[0][1] * C[1][2] - C[0][2] * C[1][1]) / Cdet;
  923. CInv[1][0] = (C[1][2] * C[2][0] - C[1][0] * C[2][2]) / Cdet;
  924. CInv[1][1] = (C[0][0] * C[2][2] - C[0][2] * C[2][0]) / Cdet;
  925. CInv[1][2] = (C[0][2] * C[1][0] - C[0][0] * C[1][2]) / Cdet;
  926. CInv[2][0] = (C[1][0] * C[2][1] - C[1][1] * C[2][0]) / Cdet;
  927. CInv[2][1] = (C[0][1] * C[2][0] - C[0][0] * C[2][1]) / Cdet;
  928. CInv[2][2] = (C[0][0] * C[1][1] - C[0][1] * C[1][0]) / Cdet;
  929. for (Int_t i = 0; i < 3; ++i)
  930. for (Int_t j = 0; j < 3; ++j)
  931. a[i] += CInv[i][j] * b[j];
  932. par[2] = a[2];
  933. par[1] = -1. / a[1];
  934. par[0] = pow(par[1], par[2]) * exp(a[0]);
  935. chi2 = sumlnq2 - 2. * (a[0] * b[0] + a[1] * b[1] + a[2] * b[2]) +
  936. a[0] * a[0] * C[0][0] + a[1] * a[1] * C[1][1] + a[2] * a[2] * C[2][2] +
  937. 2. * (a[0] * a[1] * C[0][1] + a[0] * a[2] * C[0][2] + a[1] * a[2] * C[1][2]);
  938. chi2 /= nbin - 2;
  939. Float_t A = par[0];
  940. Float_t K = a[2];
  941. Float_t tau = par[1];
  942. Float_t tau2 = tau*tau;
  943. Float_t tau3 = tau*tau2;
  944. Float_t lnT = log(par[1]);
  945. epar[0] = sqrt(fabs(A * A / (C[0][0]*(1 - a[0]) + b[0] - K * C[0][2] + C[0][1] / tau)));
  946. epar[1] = C[0][0]*((K / tau2)*(K + a[0])) -
  947. C[0][1]*((2 / (tau3))*(a[0] + 3 * K / 2)) +
  948. C[0][2] * K * K / tau2 + C[1][1] / (tau2 * tau2) -
  949. 2 * K * C[1][2] / (tau * tau * tau) + 2 * C[2][2] / (tau2 * tau2) -
  950. K * b[0] / tau2 + 2 * b[1] / tau3;
  951. epar[1] = sqrt(fabs(1. / epar[1]));
  952. epar[2] = sqrt(fabs(1. / (lnT * lnT * C[0][0] - lnT * C[0][2] + C[1][1])));
  953. }
  954. void MpdTpcClusterFinderTask::CalcResiduals(MpdTpcFoundHit* hit, Float_t sectPhi) {
  955. Int_t mcPntIdxZ = 0;
  956. Int_t mcPntIdxX = 0;
  957. TGraph* gResX = new TGraph();
  958. TGraph* gResZ = new TGraph();
  959. Int_t j, mctrackid;
  960. TpcPoint *point;
  961. const Float_t gX = hit->GetGlobalX();
  962. const Float_t gY = hit->GetGlobalY();
  963. const Float_t gZ = hit->GetGlobalZ();
  964. const Float_t lX = hit->GetLocalX();
  965. const Float_t lY = hit->GetLocalY();
  966. const Float_t lZ = hit->GetLocalZ();
  967. const Bool_t isInnerPad = Sqrt(gX * gX + gY * gY) < fSectInHeight + r_min ? kTRUE : kFALSE;
  968. const Float_t cosPhi = Cos(sectPhi);
  969. const Float_t sinPhi = Sin(sectPhi);
  970. for (j = 0; j < fMCPointArray->GetEntriesFast(); ++j) {
  971. point = (TpcPoint*) fMCPointArray->At(j);
  972. mctrackid = point->GetTrackID();
  973. if (mctrackid == hit->GetOrigin()) {
  974. gResX->SetPoint(mcPntIdxX++, -point->GetX() * sinPhi + point->GetY() * cosPhi, point->GetX() * cosPhi + point->GetY() * sinPhi - r_min);
  975. gResZ->SetPoint(mcPntIdxZ++, point->GetZ(), Sqrt(Power(point->GetX(), 2) + Power(point->GetY(), 2)));
  976. }
  977. }
  978. if (gResX->GetN() < 2 || gResZ->GetN() < 2) return;
  979. gResX->Fit("pol1", "Q");
  980. TF1* fitfcn_x = gResX->GetFunction("pol1");
  981. const Float_t bx = fitfcn_x->GetParameter(0);
  982. const Float_t ax = fitfcn_x->GetParameter(1);
  983. const Float_t xRes = (lY - bx) / ax - lX;
  984. gResZ->Fit("pol1", "Q");
  985. TF1* fitfcn_z = gResZ->GetFunction("pol1");
  986. const Float_t bz = fitfcn_z->GetParameter(0);
  987. const Float_t az = fitfcn_z->GetParameter(1);
  988. const Float_t zRes = (Sqrt(gX * gX + gY * gY) - bz) / az - gZ;
  989. fstream f;
  990. f.open(fNameResFile, fstream::out | fstream::app);
  991. f << isInnerPad << " " << xRes << " " << zRes << " " << hit->errX() << " " << hit->errZ() << " "
  992. << gX << " " << gY << " " << gZ << " " << lX << " " << lY << " " << lZ << " " << hit->QADC() << " " << hit->GetOrigin() << endl;
  993. f.close();
  994. delete gResX;
  995. delete gResZ;
  996. }
  997. void MpdTpcClusterFinderTask::CalcNewErrors(MpdTpcFoundHit* hit) {
  998. Float_t sig_preamp2 = 0.0 * 0.0; //FIXME!!!
  999. Float_t sig_prf2 = fSpread * fSpread;
  1000. Float_t Ldrift = zDrift - hit->GetGlobalZ();
  1001. Float_t Dt2 = fGas->Dt() * fGas->Dt();
  1002. Float_t Dl2 = fGas->Dl() * fGas->Dl();
  1003. Float_t tanA2 = 0.33; // tg(30) * tg(30) ---> supremum
  1004. Float_t tanB2 = 1.42; // tg(55) * tg(55) ---> supremum???
  1005. Float_t Lpad2 = 0.4 * 0.4; // TODO: get from geometry
  1006. Float_t SigT = Sqrt(Dt2 * Ldrift + sig_preamp2 + tanA2 * Lpad2 / 12); //sigma time
  1007. Float_t SigP = Sqrt(Dl2 * Ldrift + sig_prf2 + tanB2 * Lpad2 / 12); //sigma pad
  1008. hit->SetXerr(SigP * 0.4);
  1009. hit->SetZerr(SigT * 0.4);
  1010. }
  1011. void MpdTpcClusterFinderTask::Finish() {
  1012. cout << "ClusterFinder work time = " << ((Float_t)tAll) / CLOCKS_PER_SEC << endl;
  1013. if (fMakeQA) {
  1014. toDirectory("QA/TPC");
  1015. fHisto->Write();
  1016. gFile->cd();
  1017. }
  1018. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  1019. for (UInt_t iPad = 0; iPad < fNumOfPadsInRow[iRow] * 2; ++iPad)
  1020. delete [] fDigitsArray[iRow][iPad];
  1021. delete [] fDigitsArray[iRow];
  1022. }
  1023. delete [] fDigitsArray;
  1024. delete [] fNumOfPadsInRow;
  1025. //mybeg
  1026. /*
  1027. TH2F *_thist = new TH2F("thist","thist_title",300,-150,150,300,-150,150);
  1028. for(int i=0; i<fHitsArray->GetEntriesFast(); i++) {
  1029. MpdTpcHit *fHit = (MpdTpcHit*) fHitsArray->At(i);
  1030. _thist->Fill(fHit->GetX(), fHit->GetY(), fHit->GetQ());
  1031. tfile2 <<i<<") "<<"XYZ="<<fHit->GetLocalX()<<", "<<fHit->GetLocalY()<<", "<<fHit->GetLocalZ()<<"; Layer= "<<fHit->GetLayer()<< "; Pad="<<fHit->GetPad()<<"; Bin"<<fHit->GetBin()<<"\n";
  1032. }
  1033. _thist->Write("THIST", kOverwrite);
  1034. */
  1035. //myend
  1036. }