MpdTpcDigitizerAZ.cxx 36 KB

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