MpdTpcDigitizerAZlt.cxx 38 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027
  1. //-----------------------------------------------------------
  2. //
  3. // Description:
  4. // Implementation of class MpdTpcDigitizerAZlt
  5. // see MpdTpcDigitizerAZlt.h for details
  6. //
  7. // Author List:
  8. // Alexander Zinchenko LHEP, JINR, Dubna - 17-July-2015
  9. // (deep modernization of MpdTpcDigitizerTask.cxx)
  10. //
  11. //-----------------------------------------------------------
  12. // This Class' Header ------------------
  13. #include "MpdTpcDigitizerAZlt.h"
  14. // MPD Headers ----------------------
  15. #include "MpdTpcDigit.h"
  16. #include "MpdMultiField.h"
  17. #include "MpdTpcSector.h"
  18. #include "MpdTpcSectorGeo.h"
  19. #include "TpcGas.h"
  20. #include "TpcPoint.h"
  21. // FAIR Headers ----------------------
  22. #include "FairRunAna.h"
  23. #include "FairEventHeader.h"
  24. #include "FairRootManager.h"
  25. #include "FairRunSim.h"
  26. // ROOT Headers ----------------------
  27. #include "TClonesArray.h"
  28. #include <TGeoManager.h>
  29. #include "TLorentzVector.h"
  30. #include "TRandom.h"
  31. #include <TRefArray.h>
  32. #include "TMath.h"
  33. #include "TSystem.h"
  34. #include "TaskHelpers.h"
  35. #include <TVirtualFFT.h>
  36. #include <TFile.h>
  37. // C/C++ Headers ----------------------
  38. #include <math.h>
  39. #include <iostream>
  40. #include <vector>
  41. #include <algorithm>
  42. // Class Member definitions -----------
  43. using namespace std;
  44. using namespace TMath;
  45. static Int_t nOverlapDigit;
  46. static Int_t nAllLightedDigits;
  47. static clock_t tStart = 0;
  48. static clock_t tFinish = 0;
  49. static clock_t tAll = 0;
  50. //FILE *lunAZ = nullptr; //fopen("gasGain.dat","w");
  51. //---------------------------------------------------------------------------
  52. MpdTpcDigitizerAZlt::MpdTpcDigitizerAZlt()
  53. : FairTask("TPC digitizerAZlt"),
  54. fMCPointArray(nullptr),
  55. fMCTracksArray(nullptr),
  56. fDigits(nullptr),
  57. fDigits4dArray(nullptr),
  58. fSector(nullptr),
  59. fHisto(nullptr),
  60. fPRF(nullptr),
  61. fNumOfPadsInRow(nullptr),
  62. fIsHistogramsInitialized(kFALSE),
  63. fMakeQA(kFALSE),
  64. fOnlyPrimary(kFALSE),
  65. fPersistence(kTRUE),
  66. fAttach(kFALSE),
  67. fDiffuse(kTRUE),
  68. fDistort(kFALSE),
  69. fResponse(kTRUE),
  70. fDistribute(kTRUE),
  71. fPrintDebugInfo(kFALSE) {
  72. fInputBranchName = "TpcPoint";
  73. fOutputBranchName = "MpdTpcDigit";
  74. string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  75. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  76. fGas = new TpcGas(tpcGasFile, 130);
  77. }
  78. //---------------------------------------------------------------------------
  79. MpdTpcDigitizerAZlt::~MpdTpcDigitizerAZlt()
  80. {
  81. if (fIsHistogramsInitialized) {
  82. delete fHisto;
  83. }
  84. delete fGas;
  85. delete fPRF;
  86. delete fSector;
  87. }
  88. //---------------------------------------------------------------------------
  89. InitStatus MpdTpcDigitizerAZlt::Init()
  90. {
  91. //Get ROOT Manager
  92. FairRootManager* ioman = FairRootManager::Instance();
  93. if (FairRunSim::Instance() == nullptr)
  94. fMagField = FairRunAna::Instance()->GetField();
  95. else fMagField = FairRunSim::Instance()->GetField();
  96. //fMagField = FairRunSim::Instance()->GetField();
  97. if (!ioman) {
  98. cout << "\n-E- [MpdTpcDigitizerAZlt::Init]: RootManager not instantiated!" << endl;
  99. return kFATAL;
  100. }
  101. fMCPointArray = (TClonesArray*) ioman->GetObject(fInputBranchName);
  102. fMCTracksArray = (TClonesArray*) ioman->GetObject("MCTrack");
  103. fSector = new TpcSector();
  104. /*
  105. nTimeBackets = fSector->GetNTimeBins();
  106. nSectors = fSector->GetNSectors();
  107. pwIn = fSector->GetInnerPadWidth();
  108. pwOut = fSector->GetOuterPadWidth();
  109. phIn = fSector->GetInnerPadHeight();
  110. phOut = fSector->GetOuterPadHeight();
  111. nRows = fSector->GetNumRows();
  112. nInRows = fSector->GetNumInnerRows();
  113. nOutRows = fSector->GetNumOuterRows();
  114. fSectInHeight = fSector->GetSectInnerHeight();
  115. fSectHeight = fSector->GetSectHeight();
  116. r_min = fSector->GetRmin();
  117. zCathode = fSector->GetLength(); //cm
  118. */
  119. fSecGeo = MpdTpcSectorGeo::Instance();
  120. fNTimeBins = fSecGeo->GetNTimeBins();
  121. nSectors = fSecGeo->NofSectors() * 2;
  122. pwIn = fSecGeo->PadWidth(0);
  123. pwOut = fSecGeo->PadWidth(1);
  124. phIn = fSecGeo->PadHeight(0);
  125. phOut = fSecGeo->PadHeight(1);
  126. nRows = fSecGeo->NofRows();
  127. nInRows = fSecGeo->NofRowsReg(0);
  128. nOutRows = fSecGeo->NofRowsReg(1);
  129. fSectInHeight = fSecGeo->GetRocY(1) - fSecGeo->GetRocY(0);
  130. fSectHeight = fSecGeo->GetRocY(2) - fSecGeo->GetRocY(0);
  131. r_min = fSecGeo->GetMinY();
  132. zCathode = fSecGeo->GetZmax(); //cm
  133. //fNumOfPadsInRow = fSector->GetArrayPadsInRow();
  134. fNumOfPadsInRow = fSecGeo->NPadsInRows();
  135. //if (fPrintDebugInfo) {
  136. if (1) {
  137. cout << "Number of pads in every rows is ";
  138. for (UInt_t k = 0; k < nRows; ++k)
  139. cout << fNumOfPadsInRow[k] * 2 << " ";
  140. cout << endl;
  141. }
  142. //memory allocating for output array
  143. fDigits4dArray = new DigOrigArray** [nRows];
  144. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  145. fDigits4dArray[iRow] = new DigOrigArray* [fNumOfPadsInRow[iRow] * 2];
  146. for (UInt_t iPad = 0; iPad < (UInt_t) fNumOfPadsInRow[iRow] * 2; ++iPad) {
  147. fDigits4dArray[iRow][iPad] = new DigOrigArray [fNTimeBins];
  148. for (UInt_t iTime = 0; iTime < fNTimeBins; ++iTime) {
  149. fDigits4dArray[iRow][iPad][iTime].isOverlap = kFALSE;
  150. //AZ fDigits4dArray[iRow][iPad][iTime].origin = 0;
  151. fDigits4dArray[iRow][iPad][iTime].origin = -1;
  152. fDigits4dArray[iRow][iPad][iTime].origins.clear();
  153. fDigits4dArray[iRow][iPad][iTime].signal = 0.0;
  154. }
  155. }
  156. }
  157. fDigits = new TClonesArray(fOutputBranchName);
  158. ioman->Register(fOutputBranchName, "TPC", fDigits, fPersistence);
  159. //AZ fNoiseThreshold = 1000.0; // electrons
  160. //AZ-040720 fNoiseThreshold = 30.0; // ADC counts
  161. fNoiseThreshold = 3.0; // ADC counts
  162. //AZ fGain = 5000.0; //electrons
  163. fGain = 1000.0; //electrons
  164. if (fResponse) {
  165. fSpread = 0.196; // cm // Value is given by TPC group
  166. k1 = 1.0 / (Sqrt(TwoPi()) * fSpread);
  167. k2 = -0.5 / fSpread / fSpread;
  168. } else {
  169. fSpread = 0.0; // cm // FOR TEST ONLY. NO RESPONSE.
  170. k1 = k2 = 1.0;
  171. }
  172. if (!fIsHistogramsInitialized && fMakeQA) {
  173. fHisto = new MpdTpcDigitizerQAHistograms();
  174. fHisto->Initialize();
  175. fIsHistogramsInitialized = kTRUE;
  176. }
  177. fPRF = padResponseFunction();
  178. nOverlapDigit = 0;
  179. nAllLightedDigits = 0;
  180. cout << "-I- MpdTpcDigitizerAZlt: Initialization successful." << endl;
  181. return kSUCCESS;
  182. }
  183. //---------------------------------------------------------------------------
  184. void MpdTpcDigitizerAZlt::Exec(Option_t* opt)
  185. {
  186. tStart = clock();
  187. cout << "MpdTpcDigitizerAZlt::Exec started" << endl;
  188. fDigits->Delete();
  189. Int_t nPoints = fMCPointArray->GetEntriesFast();
  190. if (nPoints < 2) {
  191. Warning("MpdTpcDigitizerAZlt::Exec", "Not enough Hits in TPC for Digitization (<2)");
  192. return;
  193. }
  194. if (fPrintDebugInfo) cout << "Number of MC points is " << nPoints << endl << endl;
  195. const Float_t phiStep = TwoPi() / nSectors * 2;
  196. multimap<Double_t,Int_t>* pointsM = new multimap<Double_t,Int_t>[nSectors];
  197. for (Int_t ip = 0; ip < nPoints; ++ip) {
  198. TpcPoint* point = (TpcPoint*) fMCPointArray->UncheckedAt(ip);
  199. Float_t phi = ATan2(point->GetY(), point->GetX()); //angle in global coordinates
  200. if (phi < 0) phi += TMath::TwoPi();
  201. UInt_t isec = (UInt_t) (phi / phiStep + 0.5); //index of current sector
  202. if (isec == nSectors / 2) isec = 0;
  203. if (point->GetZ() < 0.0) isec += (nSectors / 2);
  204. pointsM[isec].insert(pair<Float_t,Int_t>(point->GetTime(),ip));
  205. }
  206. TpcPoint* virtPoint = new TpcPoint;
  207. TpcPoint tmpPoint;
  208. TpcPoint *ppp = &tmpPoint;
  209. for (UInt_t iSec = 0; iSec < nSectors; ++iSec) {
  210. if (pointsM[iSec].size() == 0) continue;
  211. multimap<Int_t,Int_t> pointID;
  212. multimap<Double_t,Int_t>::iterator mit = pointsM[iSec].begin();
  213. for ( ; mit != pointsM[iSec].end(); ++mit) {
  214. TpcPoint* point = (TpcPoint*) fMCPointArray->UncheckedAt(mit->second);
  215. pointID.insert(pair<Int_t,Int_t>(point->GetTrackID(),mit->second));
  216. }
  217. multimap<Int_t,Int_t>::iterator it = pointID.begin();
  218. TpcPoint* prePoint = (TpcPoint*) fMCPointArray->UncheckedAt(it->second);
  219. //cout << " prePoint: " << prePoint->GetTrackID() << " " << prePoint->GetX() << " "
  220. //<< prePoint->GetY() << " " << prePoint->GetZ() << endl;
  221. Check4Edge(iSec, prePoint, virtPoint); // check for edge-effect
  222. if (prePoint != virtPoint) ++it;
  223. for ( ; it != pointID.end(); ++it) {
  224. TpcPoint* curPoint = (TpcPoint*) fMCPointArray->UncheckedAt(it->second);
  225. if (prePoint == virtPoint) curPoint->SetEnergyLoss(virtPoint->GetEnergyLoss()); //AZ-090620 - corrected for entrance effect
  226. //if (curPoint->GetTrackID() != prePoint->GetTrackID()) Check4Edge(iSec, prePoint, virtPoint); // check for edge-effect
  227. //if (curPoint->GetTrackID() != prePoint->GetTrackID()) {
  228. if (curPoint->GetTrackID() != prePoint->GetTrackID() || curPoint->GetLength() - prePoint->GetLength() > 1.1) { // 14.11.2017
  229. TVector3 pos;
  230. curPoint->Position(pos);
  231. tmpPoint.SetPosition(pos);
  232. curPoint->Momentum(pos);
  233. tmpPoint.SetMomentum(pos);
  234. tmpPoint.SetTrackID(curPoint->GetTrackID());
  235. ppp = &tmpPoint;
  236. Check4Edge(iSec, ppp, virtPoint); // check for edge-effect
  237. prePoint = ppp;
  238. //cout << " --- " << curPoint->GetTrackID() << " " << curPoint->GetX() << " " << curPoint->GetY()
  239. // << " " << curPoint->GetZ() << endl;
  240. //cout << " --- " << prePoint->GetTrackID() << " " << prePoint->GetX() << " " << prePoint->GetY()
  241. // << " " << prePoint->GetZ() << endl;
  242. }
  243. TpcProcessing(prePoint, curPoint, iSec, it->second, nPoints);
  244. prePoint = curPoint;
  245. }
  246. //AZ Electronics response
  247. SignalShaping();
  248. UInt_t maxTimeBin = (UInt_t) MpdTpcSectorGeo::Instance()->TimeMax() / MpdTpcSectorGeo::Instance()->TimeBin() + 1;
  249. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  250. for (UInt_t iPad = 0; iPad < (UInt_t) fNumOfPadsInRow[iRow] * 2; ++iPad) {
  251. //AZ for (UInt_t iTime = 0; iTime < fNTimeBins; ++iTime) {
  252. for (UInt_t iTime = 0; iTime < maxTimeBin; ++iTime) {
  253. if (fDigits4dArray[iRow][iPad][iTime].signal > fNoiseThreshold) {
  254. Int_t outSize = fDigits->GetEntriesFast();
  255. Int_t id = CalcOrigin(fDigits4dArray[iRow][iPad][iTime]);
  256. if (id >= 0) new((*fDigits)[outSize]) MpdTpcDigit(id, iPad, iRow, iTime, iSec, fDigits4dArray[iRow][iPad][iTime].signal);
  257. }
  258. //if (fDigits4dArray[iRow][iPad][iTime].signal > 0.0) {
  259. if (CalcOrigin(fDigits4dArray[iRow][iPad][iTime]) >= 0) {
  260. fDigits4dArray[iRow][iPad][iTime].origins.clear();
  261. fDigits4dArray[iRow][iPad][iTime].origin = -1;
  262. fDigits4dArray[iRow][iPad][iTime].signal = 0.0;
  263. fDigits4dArray[iRow][iPad][iTime].isOverlap = kFALSE;
  264. }
  265. }
  266. }
  267. }
  268. } // for (UInt_t iSec = 0; iSec < nSectors;
  269. tFinish = clock();
  270. tAll = tAll + (tFinish - tStart);
  271. delete [] pointsM;
  272. delete virtPoint;
  273. cout << "MpdTpcDigitizerAZ::Exec finished" << endl;
  274. }
  275. //---------------------------------------------------------------------------
  276. void MpdTpcDigitizerAZlt::Check4Edge(UInt_t iSec, TpcPoint* &prePoint, TpcPoint* virtPoint)
  277. {
  278. // Check for edge-effect for track entering TPC from inside
  279. // (and correct for it if necessary)
  280. TVector3 posG, posL;
  281. prePoint->Position(posG);
  282. Int_t row0 = MpdTpcSectorGeo::Instance()->Global2Local(posG, posL, iSec%(nSectors/2));
  283. row0 = MpdTpcSectorGeo::Instance()->PadRow(row0);
  284. //cout << " Row: " << row0 << " " << iSec << " " << posL[1] << endl;
  285. if (row0) return;
  286. // For padrow == 0: create virtual point to correct for edge effect
  287. TVector3 mom, posL1;
  288. prePoint->Momentum(mom);
  289. if (mom.Pt() < 0.02) return; // do not adjust for very low-Pt tracks
  290. if (posL[1] < 0.01) return; // do not adjust - almost at the entrance
  291. posG += mom;
  292. MpdTpcSectorGeo::Instance()->Global2Local(posG, posL1, iSec%(nSectors/2));
  293. mom = posL1;
  294. mom -= posL; // momentum in sector frame
  295. if (mom[1] < 0.02) return; // do not adjust - going inward or parallel to sector lower edge
  296. Double_t scale = mom[1] / posL[1];
  297. mom.SetMag(mom.Mag() / scale);
  298. posL -= mom;
  299. //cout << posL[0] << " " << posL[1] << " " << posL[2] << endl;
  300. MpdTpcSectorGeo::Instance()->Local2Global(iSec%(nSectors/2), posL, posG);
  301. virtPoint->SetPosition(posG);
  302. virtPoint->SetTrackID(prePoint->GetTrackID());
  303. //AZ prePoint->SetEnergyLoss(prePoint->GetEnergyLoss()*1.3); // 29.10.16 - correct for edge-effect
  304. virtPoint->SetEnergyLoss(prePoint->GetEnergyLoss()*1.3); //AZ-090620 - correct for entrance effect
  305. prePoint = virtPoint;
  306. }
  307. //---------------------------------------------------------------------------
  308. //AZ Int_t MpdTpcDigitizerAZ::CalcOrigin(const DigOrigArray dig) {
  309. Int_t MpdTpcDigitizerAZlt::CalcOrigin(DigOrigArray& dig)
  310. {
  311. if (dig.origin >= 0) return dig.origin; // already done before
  312. if (dig.origins.size() == 0) return -1;
  313. Float_t max = 0.0;
  314. Int_t maxOrig = -1;
  315. if (dig.origins.size() > 1) {
  316. for (map<Int_t, Float_t>::const_iterator it = dig.origins.begin(); it != dig.origins.end(); ++it) {
  317. if (it->second > max) {
  318. maxOrig = it->first;
  319. max = it->second;
  320. }
  321. }
  322. } else {
  323. maxOrig = dig.origins.begin()->first;
  324. }
  325. dig.origin = maxOrig; //AZ
  326. return maxOrig;
  327. }
  328. //---------------------------------------------------------------------------
  329. void MpdTpcDigitizerAZlt::PadResponse(Float_t x, Float_t y, UInt_t timeID, Int_t origin, DigOrigArray ***arr)
  330. {
  331. vector<UInt_t> lightedPads;
  332. vector<UInt_t> lightedRows;
  333. vector<Float_t> amps;
  334. //Float_t avAmp = 0.0;
  335. Float_t amplSum = 0.0;
  336. Float_t amplitude = 0.0;
  337. GetArea(x, y, fSpread * 3, lightedPads, lightedRows);
  338. Double_t gain = fGain * Polya();
  339. //fprintf(lunAZ,"%f\n",gain/fGain);
  340. Int_t nPads = lightedPads.size();
  341. for (Int_t i = 0; i < nPads; ++i) {
  342. //AZ amplitude = CalculatePadResponse(lightedPads.at(i), lightedRows.at(i), x, y);
  343. amplitude = gain * CalculatePadResponse(i, nPads, lightedPads.at(i), lightedRows.at(i), x, y);
  344. amps.push_back(amplitude);
  345. amplSum += amplitude;
  346. }
  347. /*AZ
  348. if (amplSum > 0.0) {
  349. map<Int_t, Float_t>::iterator it;
  350. avAmp = fGain / amplSum; // Normalize amplitudes
  351. for (UInt_t i = 0; i < amps.size(); ++i) {
  352. arr[lightedRows.at(i)][lightedPads.at(i)][timeID].signal += (amps.at(i) * avAmp);
  353. it = arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.find(origin);
  354. if (it != arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.end()) {
  355. it->second += (amps.at(i) * avAmp);
  356. } else {
  357. arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.insert(pair<Int_t, Float_t>(origin, amps.at(i) * avAmp));
  358. }
  359. }
  360. }
  361. */
  362. if (amplSum > 0.0) {
  363. map<Int_t, Float_t>::iterator it;
  364. for (UInt_t i = 0; i < amps.size(); ++i) {
  365. arr[lightedRows.at(i)][lightedPads.at(i)][timeID].signal += amps.at(i);
  366. it = arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.find(origin);
  367. if (it != arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.end()) {
  368. it->second += amps.at(i);
  369. } else {
  370. arr[lightedRows[i]][lightedPads[i]][timeID].origins.insert(pair<Int_t, Float_t>(origin, amps[i]));
  371. }
  372. }
  373. }
  374. }
  375. //---------------------------------------------------------------------------
  376. Double_t MpdTpcDigitizerAZlt::Polya()
  377. {
  378. // Gas gain according to Polya-distribution with parameter \theta = 1.5
  379. static Int_t first = 0;
  380. static TH1D *hPolya;
  381. //long rseed;
  382. Double_t step = 0.01, shift = 0.005, param = 1.5, prob, lambda;
  383. if (first == 0) {
  384. hPolya = new TH1D("hPolya","Polya distribution",1000,0,10);
  385. Double_t param1 = 1 + param;
  386. Double_t coef = TMath::Exp(TMath::Log(param1)*param1) / TMath::Gamma(param1);
  387. for (Int_t i = 0; i < 1000; ++i) {
  388. lambda = i * step + shift;
  389. //prob = param / 0.8862 * TMath::Sqrt(param*lambda)*TMath::Exp(-param*lambda); // 0.8862 = Gamma(1.5)
  390. prob = coef * TMath::Exp(TMath::Log(lambda)*param) * TMath::Exp(-param1*lambda);
  391. hPolya->Fill(lambda,prob);
  392. }
  393. first = 1;
  394. }
  395. return hPolya->GetRandom();
  396. }
  397. //---------------------------------------------------------------------------
  398. void MpdTpcDigitizerAZlt::SignalShaping()
  399. {
  400. // Apply electronics response function
  401. static Int_t first = 0, nbins = 0, icent = 0, n2 = 0;
  402. static Double_t *reFilt = nullptr, *imFilt = nullptr;
  403. static TVirtualFFT *fft[2] = {nullptr,nullptr};
  404. const Double_t sigma = 190./2/TMath::Sqrt(2*TMath::Log(2)), sigma2 = sigma * sigma; // FWHM = 190 ns
  405. const Int_t maxTimeBin = MpdTpcSectorGeo::Instance()->TimeMax() / MpdTpcSectorGeo::Instance()->TimeBin() + 1;
  406. if (first == 0) {
  407. first = 1;
  408. nbins = MpdTpcSectorGeo::Instance()->GetNTimeBins();
  409. if (nbins % 2 == 0) --nbins;
  410. n2 = nbins / 2 + 1;
  411. icent = nbins / 2;
  412. reFilt = new Double_t [nbins];
  413. imFilt = new Double_t [nbins];
  414. for (Int_t i = 0; i < nbins; ++i) {
  415. Double_t t = (i - icent) * MpdTpcSectorGeo::Instance()->TimeBin();
  416. Double_t ampl = TMath::Exp (-t*t/2/sigma2);
  417. if (TMath::Abs(t) > 5*sigma) ampl = 0;
  418. reFilt[i] = ampl;
  419. }
  420. fft[0] = TVirtualFFT::FFT(1, &nbins, "R2C ES K");
  421. //fft[0] = TVirtualFFT::FFT(1, &nbins, "R2C EX K");
  422. fft[0]->SetPoints(reFilt);
  423. fft[0]->Transform();
  424. fft[0]->GetPointsComplex(reFilt,imFilt);
  425. }
  426. Double_t *reSig = new Double_t [nbins];
  427. Double_t *imSig = new Double_t [nbins];
  428. Double_t *reTot = new Double_t [nbins];
  429. Double_t *imTot = new Double_t [nbins];
  430. map<Int_t,Int_t> cumul; // cumulative active time bin counter
  431. const Double_t ScaleFactor = 0.083916084; //See ampl for details
  432. //AZ Int_t nRows = MpdTpcSectorGeo::Instance()->NofRows();
  433. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  434. for (UInt_t iPad = 0; iPad < (UInt_t) fNumOfPadsInRow[iRow] * 2; ++iPad) {
  435. memset(reSig,0,sizeof(Double_t)*nbins);
  436. UInt_t fired = 0;
  437. UInt_t ntbins = 0;
  438. //for (Int_t iTime = 0; iTime < nbins; ++iTime) {
  439. for (Int_t iTime = 0; iTime < maxTimeBin; ++iTime) {
  440. //if (fDigits4dArray[iRow][iPad][iTime].signal > 0) {
  441. if (CalcOrigin(fDigits4dArray[iRow][iPad][iTime]) >= 0) {
  442. // Fired channel
  443. fired = 1;
  444. if (ntbins == 0) cumul.clear();
  445. reSig[iTime] = fDigits4dArray[iRow][iPad][iTime].signal;
  446. cumul[iTime] = ++ntbins;
  447. }
  448. }
  449. if (!fired) continue;
  450. // !!! Formally expand each time bin backward by 3 bins (and forward if the next time bin is far enough) !!!
  451. for (map<Int_t,Int_t>::iterator mit = cumul.begin(); mit != cumul.end(); ++mit) {
  452. Int_t tbin = mit->first;
  453. Int_t orig = fDigits4dArray[iRow][iPad][tbin].origin;
  454. for (Int_t it = 1; it < 4; ++it) {
  455. Int_t iTime = tbin - it;
  456. if (iTime < 0) break;
  457. if (CalcOrigin(fDigits4dArray[iRow][iPad][iTime]) >= 0) break;
  458. fDigits4dArray[iRow][iPad][iTime].origin = orig;
  459. fDigits4dArray[iRow][iPad][iTime].origins[orig] = 1.0; // 1.0 - some amplitude
  460. }
  461. // Expand forward if needed
  462. map<Int_t,Int_t>::iterator next = cumul.upper_bound(tbin);
  463. if (next == cumul.end() || next->first - tbin > 4) {
  464. Int_t nexp = (next == cumul.end()) ? 4 : next->first - tbin - 3;
  465. nexp = TMath::Min (nexp, 4);
  466. // Expand
  467. for (Int_t it = 1; it < nexp; ++it) {
  468. Int_t iTime = tbin + it;
  469. //if (iTime >= nbins) break;
  470. if (iTime >= maxTimeBin) break;
  471. //if (CalcOrigin(fDigits4dArray[iRow][iPad][iTime]) >= 0) break;
  472. fDigits4dArray[iRow][iPad][iTime].origin = orig;
  473. fDigits4dArray[iRow][iPad][iTime].origins[orig] = 1.0; // 1.0 - some amplitude
  474. }
  475. }
  476. }
  477. // Fourier transform
  478. fft[0]->SetPoints(reSig);
  479. fft[0]->Transform();
  480. fft[0]->GetPointsComplex(reSig,imSig);
  481. // Convolution
  482. //for (Int_t i = 0; i < nbins; ++i) {
  483. for (Int_t i = 0; i < n2; ++i) {
  484. Double_t re = reSig[i] * reFilt[i] - imSig[i] * imFilt[i];
  485. Double_t im = reSig[i] * imFilt[i] + imSig[i] * reFilt[i];
  486. reTot[i] = re / nbins;
  487. imTot[i] = im / nbins;
  488. }
  489. // Inverse Fourier transform
  490. if (!fft[1]) fft[1] = TVirtualFFT::FFT(1, &nbins, "C2R ES K");
  491. //if (!fft[1]) fft[1] = TVirtualFFT::FFT(1, &nbins, "C2R EX K");
  492. fft[1]->SetPointsComplex(reTot,imTot);
  493. fft[1]->Transform();
  494. fft[1]->GetPoints(reTot);
  495. //AZ for (Int_t i = 0; i < nbins; ++i) {
  496. for (Int_t i = 0; i < maxTimeBin; ++i) {
  497. if (fDigits4dArray[iRow][iPad][i].origin < 0) continue; // !!! do not add extra time bins due to shaping !!!
  498. Int_t i1 = i;
  499. if (i1 <= icent) i1 += icent;
  500. else i1 -= (icent + 1);
  501. Double_t ampl = reTot[i1];
  502. //
  503. // Scale factor to adjust ADC counts
  504. ampl /= 30.0; //include division in scalefactor?
  505. //IR 11-MAR-2020 if (ampl > 4095.1) ampl = 4095.1; // dynamic range 12 bits
  506. //
  507. // IR03 new scale and rounding
  508. //Double_t ScaleFactor = 20./(550.*1.25); // 687.5
  509. // В старой шкале пик от mip для 15mm пэда был 687 канале,
  510. // тогда среднее (As=1.3MP) от рел-роста 1.01 в канале 550.*1.25*1.3,
  511. // а в правильной шкале должен быть в 75 канале
  512. // const Double_t ScaleFactor = 0.083916084; //Move above loop == 75./(550.*1.25*1.3); // 19-MAR-2020 Movchan
  513. // S.Movchan denies: if( iRow >= 26) ScaleFactor *=(1.2/1.8);
  514. ampl *= ScaleFactor;
  515. if( ampl > 1023.1) ampl = 1023.1;
  516. //
  517. fDigits4dArray[iRow][iPad][i].signal = ampl;
  518. }
  519. }
  520. }
  521. delete [] reSig;
  522. delete [] imSig;
  523. delete [] reTot;
  524. delete [] imTot;
  525. }
  526. //---------------------------------------------------------------------------
  527. void MpdTpcDigitizerAZlt::GetArea(Float_t xEll, Float_t yEll, Float_t radius, vector<UInt_t> &padIDs, vector<UInt_t> &rowIDs)
  528. {
  529. //Float_t padW = 0.0, padH = 0.0;
  530. //Float_t y = 0.0, x = 0.0;
  531. UInt_t pad = 0, row = 0;
  532. Float_t yNext; //delta = 0,
  533. //if (!fResponse) delta = -1000.0; //for test only!!!
  534. fSecGeo->PadID(xEll, yEll, row, pad, yNext);
  535. for (Int_t ip = -1; ip < 2; ++ip) {
  536. Int_t pad1 = pad + ip;
  537. if (pad1 < 0) continue;
  538. if (pad1 >= MpdTpcSectorGeo::Instance()->NPadsInRows()[row] * 2) break;
  539. padIDs.push_back(pad1);
  540. rowIDs.push_back(row);
  541. }
  542. // Add extra row
  543. if (TMath::Abs(yNext) < radius) {
  544. Int_t row1 = row;
  545. if (yNext < 0) {
  546. --row1;
  547. if (row1 < 0) return;
  548. if (fSecGeo->NPadsInRows()[row1] != fSecGeo->NPadsInRows()[row]) --pad; // different number of pads in rows
  549. } else if (yNext > 0) {
  550. ++row1;
  551. if ((UInt_t)row1 > nRows - 1) return;
  552. if (fSecGeo->NPadsInRows()[row1] != fSecGeo->NPadsInRows()[row]) ++pad; // different number of pads in rows
  553. }
  554. for (Int_t ip = -1; ip < 2; ++ip) {
  555. Int_t pad1 = pad + ip;
  556. if (pad1 < 0) continue;
  557. if (pad1 >= MpdTpcSectorGeo::Instance()->NPadsInRows()[row1] * 2) break;
  558. padIDs.push_back(pad1);
  559. rowIDs.push_back(row1);
  560. }
  561. }
  562. }
  563. //---------------------------------------------------------------------------
  564. Float_t MpdTpcDigitizerAZlt::CalculatePadResponse(Int_t iloop, Int_t nLoop, UInt_t padID, UInt_t rowID, Float_t x, Float_t y)
  565. {
  566. // Calculate pad response using lookup table (assumes the same pad width in both readout plane regions)
  567. const Int_t npads = 5, nposX = 100, nposX2 = nposX * 2; // 5 pads, 100 positions along padrow (step 25 um)
  568. const Int_t nposY = 100, nposY2 = nposY * 2;
  569. static Int_t first = 1, padID0, istep0, idist0;
  570. static UInt_t rowID0;
  571. static Double_t chargeX[npads][nposX], stepX = pwIn / nposX2;
  572. static Double_t chargeY[2][nposY]; // charges on padrow (inner and outer regions)
  573. static Double_t cy;
  574. Double_t padW = pwIn, padH = 0, distx, disty, maxX, minX, cx1, cx2, stepY, minY, maxY;
  575. Double_t sigma = fSpread;
  576. if (first) {
  577. // Compute lookup tables
  578. first = 0;
  579. // Along X
  580. for (Int_t i = 0; i < nposX; ++i) {
  581. distx = i * stepX; // distance to pad center
  582. maxX = padW / 2 - distx;
  583. for (Int_t j = 0; j < npads; ++j) {
  584. if (j == 0) maxX -= 2 * padW;
  585. minX = maxX - padW;
  586. if (j == 0) cx2 = TMath::Erf(minX/sigma);
  587. else cx2 = cx1;
  588. cx1 = TMath::Erf(maxX/sigma);
  589. chargeX[j][i] = TMath::Abs(cx1-cx2) / 2.;
  590. maxX += padW;
  591. }
  592. }
  593. // Along Y
  594. // 2 pad heights
  595. for (Int_t i = 0; i < 2; ++i) {
  596. padH = (i == 0) ? phIn : phOut;
  597. stepY = padH / nposY / 2;
  598. for (Int_t j = 0; j < nposY; ++j) {
  599. disty = j * stepY; // distance to pad center
  600. maxY = padH / 2 - disty;
  601. minY = maxY - padH;
  602. chargeY[i][j] = TMath::Abs (TMath::Erf(maxY/sigma) - TMath::Erf(minY/sigma)) / 2.;
  603. }
  604. }
  605. }
  606. //exit(0);
  607. Int_t izone = 0, idist;
  608. // Different padrow
  609. if (iloop == 0 || rowID != rowID0) {
  610. Double_t padY;
  611. rowID0 = rowID;
  612. if (rowID < nInRows) {
  613. //padW = pwIn;
  614. padH = phIn;
  615. padY = padH * ((Double_t)rowID + 0.5); // y-coordinate of pad center
  616. } else {
  617. //padW = pwOut;
  618. padH = phOut;
  619. padY = fSectInHeight + (((Double_t)rowID - nInRows) + 0.5) * padH; // y-coordinate of pad center
  620. izone = 1;
  621. }
  622. disty = y - padY;
  623. Int_t istep = TMath::Abs (disty / (padH / nposY2));
  624. if (istep < nposY) cy = chargeY[izone][istep];
  625. else {
  626. istep = nposY2 - istep;
  627. cy = 1.0 - chargeY[izone][istep];
  628. }
  629. Double_t padX = padW * ((Double_t)padID - fNumOfPadsInRow[rowID] + 0.5); // x-coordinate of pad center
  630. padID0 = padID;
  631. distx = x - padX;
  632. idist = TMath::Nint (distx / stepX);
  633. idist0 = idist;
  634. istep0 = TMath::Abs (idist);
  635. istep0 %= nposX2;
  636. if (istep0 > nposX) istep0 = TMath::Abs(istep0-nposX2);
  637. istep0 = TMath::Min (istep0, nposX-1);
  638. } else idist = idist0 - (padID-padID0) * nposX2;
  639. Int_t ipad = TMath::Abs (idist / nposX);
  640. if (ipad == 0) ipad = npads / 2; // central pad
  641. //else if (idist > 0) {
  642. else if (1) {
  643. if (ipad % 2 == 0) ipad = npads / 2 - ipad / 2;
  644. else ipad = npads / 2 + (ipad+1) / 2;
  645. } else {
  646. if (ipad % 2 != 0) ipad = npads / 2 + (ipad+1) / 2;
  647. else ipad = npads / 2 - ipad / 2;
  648. }
  649. cx1 = chargeX[ipad][istep0];
  650. Double_t ctot = TMath::Abs(cy) * cx1;
  651. //if (rowID == 1) cout << " Row, pad: " << rowID << " " << padID << " " << x << " " << istep0 << " " << idist << " " << ipad << " " << cx1 << " " << ctot << endl;
  652. return ctot;
  653. }
  654. //---------------------------------------------------------------------------
  655. TF1* MpdTpcDigitizerAZlt::padResponseFunction() {
  656. if (fPRF)
  657. return fPRF;
  658. fPRF = new TF1("Gaus PRF", "gaus", -5, 5);
  659. fPRF->SetParameter(0, 1.0 / (sqrt(2.0 * TMath::Pi()) * fSpread));
  660. fPRF->SetParameter(1, 0);
  661. fPRF->SetParameter(2, fSpread);
  662. return fPRF;
  663. }
  664. //---------------------------------------------------------------------------
  665. Bool_t MpdTpcDigitizerAZlt::isSubtrackInInwards(const TpcPoint *p1, const TpcPoint *p2) { //WHAT AM I DOING???
  666. const Float_t x1 = p1->GetX();
  667. const Float_t x2 = p2->GetX();
  668. const Float_t y1 = p1->GetY();
  669. const Float_t y2 = p2->GetY();
  670. const Float_t a = (y1 - y2) / (x1 - x2);
  671. const Float_t b = (y1 * x2 - x1 * y2) / (x2 - x1);
  672. const Float_t minR = fabs(b) / sqrt(a * a + 1);
  673. if (minR < r_min) //then check if minimal distance is between our points
  674. {
  675. const Float_t x = -a * b / (a * a + 1);
  676. const Float_t y = b / (a * a + 1);
  677. if ((x1 - x) * (x2 - x) < 0 && (y1 - y) * (y2 - y) < 0) {
  678. return kTRUE;
  679. }
  680. }
  681. return kFALSE;
  682. }
  683. //---------------------------------------------------------------------------
  684. void MpdTpcDigitizerAZlt::TpcProcessing(const TpcPoint* prePoint, const TpcPoint* curPoint,
  685. const UInt_t secID, const UInt_t iPoint, const UInt_t nPoints)
  686. {
  687. Float_t dE = 0.0; //energy loss
  688. UInt_t qTotal = 0; //sum of clusters charges (=sum of electrons between two TpcPoints)
  689. UInt_t qCluster = 0; //charge of cluster (= number of electrons)
  690. TLorentzVector curPointPos; // coordinates for current TpcPoint
  691. TLorentzVector prePointPos; // coordinates for previous TpcPoint
  692. TLorentzVector diffPointPos; // steps for clusters creation
  693. TVector3 diffuse; // vector of diffuse for every coordinates
  694. TVector3 distort; // vector of distortion for every coordinates
  695. TLorentzVector electronPos; // coordinates for created electrons
  696. TLorentzVector clustPos; // coordinates for created clusters
  697. Float_t driftl = 0.0; // length for drifting
  698. vector<UInt_t> clustArr; // vector of clusters between two TpcPoints
  699. Float_t localX = 0.0, localY = 0.0; //local coordinates of electron (sector coordinates)
  700. MpdTpcSectorGeo* secGeo = MpdTpcSectorGeo::Instance();
  701. if ( fPrintDebugInfo && (iPoint % 1000 == 0) ) cout << UInt_t(iPoint * 1.0 / nPoints * 100.0) << " % of TPC points processed" << endl;
  702. // curPoint = (TpcPoint*) fMCPointArray->At(i);
  703. if (fOnlyPrimary == kTRUE) {
  704. MpdMCTrack* tr = (MpdMCTrack*) fMCTracksArray->At(curPoint->GetTrackID());
  705. if (tr->GetMotherId() != -1) return;
  706. }
  707. //check if hits are on the same track
  708. if (curPoint->GetTrackID() == prePoint->GetTrackID() && !isSubtrackInInwards(prePoint, curPoint)) {
  709. dE = curPoint->GetEnergyLoss() * 1E9; //convert from GeV to eV
  710. if (dE < 0) {
  711. Error("MpdTpcDigitizerTask::Exec", "Negative Energy loss!");
  712. return;
  713. }
  714. qTotal = (UInt_t) floor(fabs(dE / fGas->W()));
  715. if (qTotal == 0) return;
  716. curPointPos.SetXYZT(curPoint->GetX(), curPoint->GetY(), curPoint->GetZ(), curPoint->GetTime());
  717. prePointPos.SetXYZT(prePoint->GetX(), prePoint->GetY(), prePoint->GetZ(), prePoint->GetTime());
  718. if ((curPointPos.T() < 0) || (prePointPos.T() < 0)) {
  719. Error("MpdTpcDigitizerTask::Exec", "Negative Time!");
  720. return;
  721. }
  722. diffPointPos = curPointPos - prePointPos; //differences between two points by coordinates
  723. //AZ diffPointPos *= (1 / diffPointPos.Vect().Mag()); //directional cosines //TODO! Do we need this??? Look at line #297
  724. //while still charge not used-up distribute charge into next cluster
  725. if (fDistribute) {
  726. while (qTotal > 0) {
  727. //roll dice for next cluster
  728. qCluster = fGas->GetRandomCSUniform();
  729. if (qCluster > qTotal) qCluster = qTotal;
  730. qTotal -= qCluster;
  731. clustArr.push_back(qCluster);
  732. }// finish loop for cluster creation
  733. } else {
  734. clustArr.push_back(qTotal); // JUST FOR TEST. NO CLUSTER DISTRIBUTION!
  735. // clustArr.push_back(1); // JUST FOR TEST. NO CLUSTER DISTRIBUTION ONLY ONE ELECTRON IN CLUSTER!
  736. }
  737. //AZ diffPointPos *= (diffPointPos.Vect().Mag() / clustArr.size()); //now here are steps between clusters by coordinates TODO: correct distribution
  738. diffPointPos *= (1./ clustArr.size());
  739. clustPos = prePointPos;
  740. for (UInt_t iClust = 0; iClust < clustArr.size(); ++iClust) {
  741. clustPos += diffPointPos;
  742. driftl = zCathode - fabs(clustPos.Z());
  743. if (driftl <= 0) continue;
  744. for (UInt_t iEll = 0; iEll < clustArr.at(iClust); ++iEll) {
  745. //attachment
  746. if (fAttach)
  747. if (exp(-driftl * fGas->k()) < gRandom->Uniform()) continue; // FIXME
  748. //diffusion
  749. if (fDiffuse) {
  750. const Float_t sqrtDrift = sqrt(driftl);
  751. const Float_t sigmat = fGas->Dt() * sqrtDrift;
  752. const Float_t sigmal = fGas->Dl() * sqrtDrift;
  753. diffuse.SetXYZ(gRandom->Gaus(0, sigmat), gRandom->Gaus(0, sigmat), gRandom->Gaus(0, sigmal));
  754. }
  755. if (fDistort) {
  756. const Float_t dt = 1E-03; //time step [s]
  757. const Float_t mu = 4.23; //electron mobility [m^2 / s / V]
  758. const Float_t mu2 = mu * mu; //just square of mu
  759. const TVector3 E(0.0, 0.0, fGas->E() * 100); // vector of electric field components (now is constant and parallel to Z axes) // 100 - convert Ez from V/cm to V/m
  760. TVector3 B(0.0, 0.0, 0.0); // vector of magnetic field components
  761. TVector3 v; // vector of current velocity components
  762. TVector3 posCur; // vector of current position components
  763. TVector3 posPre = clustPos.Vect(); // vector of previous position components
  764. TVector3 EBCross(0.0, 0.0, 0.0); // vector product of E and B vectors
  765. Bool_t inTpc = kTRUE;
  766. while (inTpc) {
  767. B.SetXYZ(fMagField->GetBx(posPre.X(), posPre.Y(), posPre.Z()) * 0.1, fMagField->GetBy(posPre.X(), posPre.Y(), posPre.Z()) * 0.1, fMagField->GetBz(posPre.X(), posPre.Y(), posPre.Z()) * 0.1);
  768. EBCross = E.Cross(B);
  769. v = mu / (1 + mu2 * B.Mag2()) * (E - mu * EBCross + mu2 * B * (E * B));
  770. posCur = v * dt + posPre;
  771. // cout << "X = " << posCur.X() << " Y = " << posCur.Y() << " Z = " << posCur.Z() << " Vx = " << v.X() << " Vy = " << v.Y() << " Vz = " << v.Z() << endl;
  772. if ((posCur.Perp() > r_min + fSectHeight) || (posCur.Perp() < r_min) || (Abs(posCur.Z()) > zCathode)) inTpc = kFALSE;
  773. posPre = posCur;
  774. }
  775. distort.SetX(posCur.X() - clustPos.Vect().X());
  776. distort.SetY(posCur.Y() - clustPos.Vect().Y());
  777. distort.SetZ(0.0); //FIXME
  778. }
  779. electronPos.SetVect(clustPos.Vect() + diffuse + distort);
  780. electronPos.SetT(clustPos.T() + (zCathode - fabs(electronPos.Z())) / fGas->VDrift()); // Do we need to use clustPos.T() ???
  781. /*AZ
  782. const Float_t phiStep = TwoPi() / nSectors * 2;
  783. const Float_t sectPhi = secID * phiStep;
  784. localY = electronPos.X() * Cos(sectPhi) + electronPos.Y() * Sin(sectPhi) - r_min; //converting from global to local coordinates
  785. localX = -electronPos.X() * Sin(sectPhi) + electronPos.Y() * Cos(sectPhi); //converting from global to local coordinates
  786. if ((localY < 0.0) || (Abs(localX) > fSector->GetMaxX()) || (localY > fSectHeight)) continue; //FIXME!!!
  787. const Float_t timeStep = (zCathode / fNTimeBins) / fGas->VDrift();
  788. const UInt_t curTimeID = (UInt_t) ((zCathode / fGas->VDrift() - electronPos.T()) / timeStep);
  789. */
  790. //AZ
  791. TVector3 xyzLoc;
  792. if (secGeo->Global2Local(electronPos.Vect(), xyzLoc, secID % secGeo->NofSectors()) < 0) continue;
  793. localX = xyzLoc[0];
  794. localY = xyzLoc[1];
  795. if (!TMath::Finite(localX) || !TMath::Finite(localY)) {
  796. // Debug
  797. cout << " !!! Not finite " << secID % secGeo->NofSectors() << endl;
  798. electronPos.Vect().Print();
  799. diffuse.Print();
  800. }
  801. UInt_t curTimeID = UInt_t (secGeo->T2TimeBin(electronPos.T()));
  802. //cout << localX << " " << xyzLoc[0] << " " << localY << " " << xyzLoc[1] << " " << curTimeID << " " << timeBin << endl;
  803. //AZ
  804. if (curTimeID >= fNTimeBins) continue;
  805. if (fMakeQA) {
  806. fHisto->_hX_local->Fill(localX);
  807. fHisto->_hY_local->Fill(localY);
  808. fHisto->_hXY_local->Fill(localX, localY);
  809. fHisto->_hYZ_local->Fill(fabs(electronPos.Z()), localY);
  810. fHisto->_hXY_global->Fill(electronPos.X(), electronPos.Y());
  811. fHisto->_hRZ_global->Fill(electronPos.Z(), electronPos.Y());
  812. fHisto->_h3D_el->Fill(electronPos.X(), electronPos.Y(), electronPos.Z());
  813. fHisto->_hZ_local->Fill(fabs(electronPos.Z()));
  814. fHisto->_hX_global->Fill(electronPos.X());
  815. fHisto->_hY_global->Fill(electronPos.Y());
  816. fHisto->_hZ_global->Fill(electronPos.Z());
  817. fHisto->_hDiffuseXY->Fill(diffuse.X(), diffuse.Y());
  818. fHisto->_hDistortXY->Fill(distort.X(), distort.Y());
  819. }
  820. Int_t origin = prePoint->GetTrackID();
  821. PadResponse(localX, localY, curTimeID, origin, fDigits4dArray);
  822. } // for (UInt_t iEll = 0; iEll < clustArr.at(iClust);
  823. } // for (UInt_t iClust = 0; iClust < clustArr.size();
  824. } // if (curPoint->GetTrackID() == prePoint->GetTrackID()
  825. clustArr.clear();
  826. clustArr.resize(0);
  827. }
  828. //---------------------------------------------------------------------------
  829. void MpdTpcDigitizerAZlt::Finish() {
  830. cout << "Digitizer work time = " << ((Float_t)tAll) / CLOCKS_PER_SEC << endl;
  831. if (fMakeQA) {
  832. toDirectory("QA/TPC");
  833. Float_t digit = 0.0;
  834. UInt_t iPad_shifted = 0; //needed for correct drawing of fDigitsArray
  835. for (UInt_t iRows = 0; iRows < nRows; ++iRows) {
  836. for (UInt_t iPads = 0; iPads < (UInt_t) fNumOfPadsInRow[iRows] * 2; ++iPads) {
  837. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[iRows];
  838. for (UInt_t iTime = 0; iTime < fNTimeBins; ++iTime) {
  839. digit = fDigits4dArray[iRows][iPads][iTime].signal;
  840. fHisto->_hXY_dig->Fill(iPad_shifted, iRows, digit);
  841. // fHisto->_hSect_dig->Fill(iSect, digit);
  842. fHisto->_hX_dig->Fill(iPad_shifted, digit);
  843. fHisto->_hY_dig->Fill(iRows, digit);
  844. fHisto->_hZ_dig->Fill(iTime, digit);
  845. fHisto->_h3D_dig->Fill(iPad_shifted, iRows, iTime, digit);
  846. if (digit > 1000.0) fHisto->_hADC_dig->Fill(digit);
  847. }
  848. }
  849. }
  850. for (UInt_t iRows = 0; iRows < nRows; ++iRows) {
  851. for (UInt_t iPads = 0; iPads < (UInt_t) fNumOfPadsInRow[iRows] * 2; ++iPads) {
  852. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[iRows];
  853. for (UInt_t iTime = 0; iTime < fNTimeBins; ++iTime) {
  854. digit = fDigits4dArray[iRows][iPads][iTime].signal;
  855. //pad activity
  856. //if (digit > 1000.0) {
  857. // fHisto->_hXY_dig->Fill(iPad_shifted, iRows, 1.0);
  858. //}
  859. // fHisto->_hXY_dig->Fill(iPad_shifted, iRows, digit);
  860. fHisto->_h3D_dig->Fill(iPad_shifted, iRows, iTime, digit);
  861. }
  862. }
  863. }
  864. for (UInt_t iTime = 0; iTime < fNTimeBins; ++iTime) {
  865. for (UInt_t iPads = 0; iPads < (UInt_t) fNumOfPadsInRow[1] * 2; ++iPads) {
  866. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[1];
  867. digit = fDigits4dArray[1][iPads][iTime].signal;
  868. fHisto->_hXT_dig_1->Fill(iPad_shifted, iTime, digit);
  869. }
  870. for (UInt_t iPads = 0; iPads < (UInt_t) fNumOfPadsInRow[5] * 2; ++iPads) {
  871. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[5];
  872. digit = fDigits4dArray[5][iPads][iTime].signal;
  873. fHisto->_hXT_dig_5->Fill(iPad_shifted, iTime, digit);
  874. }
  875. for (UInt_t iPads = 0; iPads < (UInt_t) fNumOfPadsInRow[10] * 2; ++iPads) {
  876. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[10];
  877. digit = fDigits4dArray[10][iPads][iTime].signal;
  878. fHisto->_hXT_dig_10->Fill(iPad_shifted, iTime, digit);
  879. }
  880. for (UInt_t iPads = 0; iPads < (UInt_t) fNumOfPadsInRow[20] * 2; ++iPads) {
  881. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[20];
  882. digit = fDigits4dArray[20][iPads][iTime].signal;
  883. fHisto->_hXT_dig_20->Fill(iPad_shifted, iTime, digit);
  884. }
  885. for (UInt_t iPads = 0; iPads < (UInt_t) fNumOfPadsInRow[40] * 2; ++iPads) {
  886. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[40];
  887. digit = fDigits4dArray[40][iPads][iTime].signal;
  888. fHisto->_hXT_dig_40->Fill(iPad_shifted, iTime, digit);
  889. }
  890. }
  891. fHisto->Write();
  892. gFile->cd();
  893. }
  894. }
  895. //---------------------------------------------------------------------------
  896. ClassImp(MpdTpcDigitizerAZlt)