MpdTpcDigitizerTask.cxx 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617
  1. //-----------------------------------------------------------
  2. //
  3. // Description:
  4. // Implementation of class MpdTpcDigitizerTask
  5. // see MpdTpcDigitizerTask.h for details
  6. //
  7. // Author List:
  8. // Sergey Merts
  9. //
  10. //-----------------------------------------------------------
  11. // Panda Headers ----------------------
  12. // This Class' Header ------------------
  13. #include "MpdTpcDigitizerTask.h"
  14. // C/C++ Headers ----------------------
  15. #include <math.h>
  16. #include <iostream>
  17. #include <vector>
  18. #include <algorithm>
  19. #include "MpdMultiField.h"
  20. #include "FairRunAna.h"
  21. #include "FairEventHeader.h"
  22. #include "TpcPoint.h"
  23. #include "TLorentzVector.h"
  24. #include "FairRootManager.h"
  25. #include "FairRunSim.h"
  26. #include <TGeoManager.h>
  27. #include <TRefArray.h>
  28. #include "TClonesArray.h"
  29. #include "TpcGas.h"
  30. #include "TRandom.h"
  31. #include "TMath.h"
  32. #include "MpdTpcSector.h"
  33. #include "TSystem.h"
  34. #include "TaskHelpers.h"
  35. #include "MpdTpcDigit.h"
  36. #include "TFile.h"
  37. // Class Member definitions -----------
  38. using namespace std;
  39. using namespace TMath;
  40. static Int_t nOverlapDigit;
  41. static Int_t nAllLightedDigits;
  42. static clock_t tStart = 0;
  43. static clock_t tFinish = 0;
  44. static clock_t tAll = 0;
  45. MpdTpcDigitizerTask::MpdTpcDigitizerTask() :
  46. fPersistence(kTRUE),
  47. fResponse(kTRUE),
  48. fDistribute(kTRUE),
  49. fAttach(kFALSE),
  50. fDiffuse(kFALSE),
  51. fDistort(kFALSE),
  52. fPrintDebugInfo(kFALSE),
  53. fIsHistogramsInitialized(kFALSE),
  54. fMakeQA(kFALSE),
  55. fHisto(nullptr),
  56. fPRF(nullptr),
  57. fNumOfPadsInRow(nullptr),
  58. fMCPointArray(nullptr),
  59. fMCTracksArray(nullptr),
  60. fDigits(nullptr),
  61. fSector(nullptr),
  62. fDigits4dArray(nullptr) {
  63. fInputBranchName = "TpcPoint";
  64. fOutputBranchName = "MpdTpcDigit";
  65. string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  66. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  67. fGas = new TpcGas(tpcGasFile, 130);
  68. }
  69. MpdTpcDigitizerTask::~MpdTpcDigitizerTask() {
  70. if (fIsHistogramsInitialized) {
  71. delete fHisto;
  72. }
  73. delete fGas;
  74. delete fPRF;
  75. delete fSector;
  76. }
  77. InitStatus MpdTpcDigitizerTask::Init() {
  78. //Get ROOT Manager
  79. FairRootManager* ioman = FairRootManager::Instance();
  80. fMagField = FairRunSim::Instance()->GetField();
  81. if (!ioman) {
  82. cout << "\n-E- [MpdTpcDigitizerTask::Init]: RootManager not instantiated!" << endl;
  83. return kFATAL;
  84. }
  85. fMCPointArray = (TClonesArray*) ioman->GetObject(fInputBranchName);
  86. fMCTracksArray = (TClonesArray*) ioman->GetObject("MCTrack");
  87. fSector = new TpcSector();
  88. nTimeBackets = fSector->GetNTimeBins();
  89. nSectors = fSector->GetNSectors();
  90. pwIn = fSector->GetInnerPadWidth();
  91. pwOut = fSector->GetOuterPadWidth();
  92. phIn = fSector->GetInnerPadHeight();
  93. phOut = fSector->GetOuterPadHeight();
  94. nRows = fSector->GetNumRows();
  95. nInRows = fSector->GetNumInnerRows();
  96. nOutRows = fSector->GetNumOuterRows();
  97. fSectInHeight = fSector->GetSectInnerHeight();
  98. fSectHeight = fSector->GetSectHeight();
  99. r_min = fSector->GetRmin();
  100. zCathode = fSector->GetLength(); //cm
  101. fNumOfPadsInRow = fSector->GetArrayPadsInRow();
  102. if (fPrintDebugInfo) {
  103. cout << "Number of pads in every rows is ";
  104. for (UInt_t k = 0; k < nRows; ++k)
  105. cout << fNumOfPadsInRow[k] * 2 << " ";
  106. cout << endl;
  107. }
  108. //memory allocating for output array
  109. fDigits4dArray = new DigOrigArray** [nRows];
  110. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  111. fDigits4dArray[iRow] = new DigOrigArray* [fNumOfPadsInRow[iRow] * 2];
  112. for (UInt_t iPad = 0; iPad < fNumOfPadsInRow[iRow] * 2; ++iPad) {
  113. fDigits4dArray[iRow][iPad] = new DigOrigArray [nTimeBackets];
  114. for (UInt_t iTime = 0; iTime < nTimeBackets; ++iTime) {
  115. fDigits4dArray[iRow][iPad][iTime].isOverlap = kFALSE;
  116. fDigits4dArray[iRow][iPad][iTime].origin = 0;
  117. fDigits4dArray[iRow][iPad][iTime].origins.clear();
  118. fDigits4dArray[iRow][iPad][iTime].signal = 0.0;
  119. }
  120. }
  121. }
  122. fDigits = new TClonesArray(fOutputBranchName);
  123. ioman->Register(fOutputBranchName, "TPC", fDigits, fPersistence);
  124. fNoiseThreshold = 1000.0; // electrons
  125. fGain = 5000.0; //electrons
  126. if (fResponse) {
  127. fSpread = 0.196; // cm // Value is given by TPC group
  128. k1 = 1.0 / (Sqrt(TwoPi()) * fSpread);
  129. k2 = -0.5 / fSpread / fSpread;
  130. } else {
  131. fSpread = 0.0; // cm // FOR TEST ONLY. NO RESPONSE.
  132. k1 = k2 = 1.0;
  133. }
  134. if (!fIsHistogramsInitialized && fMakeQA) {
  135. fHisto = new MpdTpcDigitizerQAHistograms();
  136. fHisto->Initialize();
  137. fIsHistogramsInitialized = kTRUE;
  138. }
  139. fPRF = padResponseFunction();
  140. nOverlapDigit = 0;
  141. nAllLightedDigits = 0;
  142. cout << "-I- MpdTpcDigitizerTask: Initialization successful." << endl;
  143. return kSUCCESS;
  144. }
  145. void MpdTpcDigitizerTask::Exec(Option_t* opt) {
  146. tStart = clock();
  147. cout << "MpdTpcDigitizer::Exec started" << endl;
  148. fDigits->Delete();
  149. Int_t nPoints = fMCPointArray->GetEntriesFast();
  150. if (nPoints < 2) {
  151. Warning("MpdTpcDigitizerTask::Exec", "Not enough Hits in TPC for Digitization (<2)");
  152. return;
  153. }
  154. if (fPrintDebugInfo) cout << "Number of MC points is " << nPoints << endl << endl;
  155. const Float_t phiStep = TwoPi() / nSectors * 2;
  156. for (UInt_t iSec = 0; iSec < nSectors; ++iSec) {
  157. //get memory for array
  158. TpcPoint* prePoint = (TpcPoint*) fMCPointArray->At(0);
  159. for (UInt_t i = 1; i < nPoints; i++) {
  160. TpcPoint* curPoint = (TpcPoint*) fMCPointArray->At(i);
  161. Float_t globPhi = ATan2(curPoint->GetY(), curPoint->GetX()); //angle in global coordinates
  162. if (globPhi < 0) globPhi += TwoPi();
  163. UInt_t curSectID = (UInt_t) (globPhi / phiStep + 0.5); //index of current sector
  164. if (curSectID == nSectors / 2) curSectID = 0;
  165. if (curPoint->GetZ() < 0.0) curSectID += (nSectors / 2);
  166. if (curSectID == iSec) { // is current point in sector
  167. TpcProcessing(prePoint, curPoint, iSec, i, nPoints);
  168. }
  169. prePoint = curPoint;
  170. }
  171. for (UInt_t iRow = 0; iRow < nRows; ++iRow) {
  172. for (UInt_t iPad = 0; iPad < fNumOfPadsInRow[iRow] * 2; ++iPad) {
  173. for (UInt_t iTime = 0; iTime < nTimeBackets; ++iTime) {
  174. if (fDigits4dArray[iRow][iPad][iTime].signal > fNoiseThreshold) {
  175. Int_t outSize = fDigits->GetEntriesFast();
  176. new((*fDigits)[outSize]) MpdTpcDigit(CalcOrigin(fDigits4dArray[iRow][iPad][iTime]), iPad, iRow, iTime, iSec, fDigits4dArray[iRow][iPad][iTime].signal);
  177. }
  178. if (fDigits4dArray[iRow][iPad][iTime].signal > 0.0) {
  179. fDigits4dArray[iRow][iPad][iTime].origins.clear();
  180. fDigits4dArray[iRow][iPad][iTime].origin = -1;
  181. fDigits4dArray[iRow][iPad][iTime].signal = 0.0;
  182. fDigits4dArray[iRow][iPad][iTime].isOverlap = kFALSE;
  183. }
  184. }
  185. }
  186. }
  187. }
  188. tFinish = clock();
  189. tAll = tAll + (tFinish - tStart);
  190. cout << "MpdTpcDigitizer::Exec finished" << endl;
  191. }
  192. Int_t MpdTpcDigitizerTask::CalcOrigin(const DigOrigArray dig) {
  193. Float_t max = 0.0;
  194. Int_t maxOrig = -1;
  195. if (dig.origins.size() > 1) {
  196. for (map<Int_t, Float_t>::const_iterator it = dig.origins.begin(); it != dig.origins.end(); ++it) {
  197. if (it->second > max) {
  198. maxOrig = it->first;
  199. max = it->second;
  200. }
  201. }
  202. } else {
  203. maxOrig = dig.origins.begin()->first;
  204. }
  205. return maxOrig;
  206. }
  207. void MpdTpcDigitizerTask::PadResponse(Float_t x, Float_t y, UInt_t timeID, Int_t origin, DigOrigArray ***arr) {
  208. vector<UInt_t> lightedPads;
  209. vector<UInt_t> lightedRows;
  210. vector<Float_t> amps;
  211. Float_t avAmp = 0.0;
  212. Float_t amplSum = 0.0;
  213. Float_t amplitude = 0.0;
  214. GetArea(x, y, fSpread * 3, lightedPads, lightedRows);
  215. for (UInt_t i = 0; i < lightedPads.size(); ++i) {
  216. amplitude = CalculatePadResponse(lightedPads.at(i), lightedRows.at(i), x, y);
  217. amps.push_back(amplitude);
  218. amplSum += amplitude;
  219. }
  220. if (amplSum > 0.0) {
  221. map<Int_t, Float_t>::iterator it;
  222. avAmp = fGain / amplSum; // Normalize amplitudes
  223. for (UInt_t i = 0; i < amps.size(); ++i) {
  224. arr[lightedRows.at(i)][lightedPads.at(i)][timeID].signal += (amps.at(i) * avAmp);
  225. it = arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.find(origin);
  226. if (it != arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.end()) {
  227. it->second += (amps.at(i) * avAmp);
  228. } else {
  229. arr[lightedRows.at(i)][lightedPads.at(i)][timeID].origins.insert(pair<Int_t, Float_t>(origin, amps.at(i) * avAmp));
  230. }
  231. }
  232. }
  233. }
  234. void MpdTpcDigitizerTask::GetArea(Float_t xEll, Float_t yEll, Float_t radius, vector<UInt_t> &padIDs, vector<UInt_t> &rowIDs) {
  235. Float_t padW = 0.0, padH = 0.0;
  236. Float_t y = 0.0, x = 0.0;
  237. UInt_t pad = 0, row = 0;
  238. Float_t delta;
  239. if (fResponse) delta = 0.0;
  240. else delta = -1000.0; //for test only!!!
  241. if (yEll > nInRows * phIn + radius) {
  242. padW = pwOut;
  243. padH = phOut;
  244. x = xEll - radius - padW;
  245. do {
  246. x += padW;
  247. y = yEll - radius - padH;
  248. do {
  249. y += padH;
  250. row = (UInt_t) ((y - fSectInHeight) / padH) + nInRows;
  251. pad = (x > 0.0) ? (fNumOfPadsInRow[row] + floor(x / padW)) : (fNumOfPadsInRow[row] - 1 + ceil(x / padW));
  252. if (row >= nRows || pad >= fNumOfPadsInRow[row] * 2) continue;
  253. padIDs.push_back(pad);
  254. rowIDs.push_back(row);
  255. } while (y < yEll + radius + delta);
  256. } while (x < xEll + radius + delta);
  257. } else if (yEll < nInRows * phIn - radius) {
  258. padW = pwIn;
  259. padH = phIn;
  260. x = xEll - radius - padW;
  261. do {
  262. x += padW;
  263. y = yEll - radius - padH;
  264. do {
  265. y += padH;
  266. row = (UInt_t) (y / padH);
  267. pad = (x > 0.0) ? (fNumOfPadsInRow[row] + floor(x / padW)) : (fNumOfPadsInRow[row] - 1 + ceil(x / padW));
  268. if (row >= nRows || pad >= fNumOfPadsInRow[row] * 2) continue;
  269. padIDs.push_back(pad);
  270. rowIDs.push_back(row);
  271. } while (y < yEll + radius + delta);
  272. } while (x < xEll + radius + delta);
  273. } else {
  274. x = xEll - radius - pwIn;
  275. do {
  276. x += pwIn;
  277. row = nInRows - 1;
  278. pad = (x > 0.0) ? (fNumOfPadsInRow[row] + floor(x / pwIn)) : (fNumOfPadsInRow[row] - 1 + ceil(x / pwIn));
  279. if (pad >= fNumOfPadsInRow[row] * 2) continue;
  280. padIDs.push_back(pad);
  281. rowIDs.push_back(row);
  282. } while (x < xEll + radius + delta);
  283. do {
  284. x += pwOut;
  285. row = nInRows;
  286. pad = (x > 0.0) ? (fNumOfPadsInRow[row] + floor(x / pwOut)) : (fNumOfPadsInRow[row] - 1 + ceil(x / pwOut));
  287. if (pad >= fNumOfPadsInRow[row] * 2) continue;
  288. padIDs.push_back(pad);
  289. rowIDs.push_back(row);
  290. } while (x < xEll + radius + delta);
  291. }
  292. }
  293. Float_t MpdTpcDigitizerTask::CalculatePadResponse(UInt_t padID, UInt_t rowID, Float_t x, Float_t y) {
  294. Float_t padX, padY;
  295. Float_t padW, padH;
  296. if (rowID < nInRows) {
  297. padW = pwIn;
  298. padH = phIn;
  299. padX = padW * ((Float_t) padID - (Float_t) fNumOfPadsInRow[rowID] + 0.5); // x-coordinate of pad center
  300. padY = padH * ((Float_t) rowID + 0.5); // y-coordinate of pad center
  301. } else {
  302. padW = pwOut;
  303. padH = phOut;
  304. padX = padW * ((Float_t) padID - (Float_t) fNumOfPadsInRow[rowID] + 0.5); // x-coordinate of pad center
  305. padY = fSectInHeight + ((Float_t) (rowID - nInRows) + 0.5) * padH; // y-coordinate of pad center
  306. }
  307. const Float_t maxX = x - (padX - padW / 2);
  308. const Float_t minX = x - (padX + padW / 2);
  309. const Float_t maxY = y - (padY - padH / 2);
  310. const Float_t minY = y - (padY + padH / 2);
  311. const Float_t i1 = (Exp(k2 * minX * minX) + Exp(k2 * maxX * maxX)) * k1 * padW / 2;
  312. const Float_t i2 = (Exp(k2 * minY * minY) + Exp(k2 * maxY * maxY)) * k1 * padH / 2;
  313. return i1 * i2;
  314. }
  315. TF1* MpdTpcDigitizerTask::padResponseFunction() {
  316. if (fPRF)
  317. return fPRF;
  318. fPRF = new TF1("Gaus PRF", "gaus", -5, 5);
  319. fPRF->SetParameter(0, 1.0 / (sqrt(2.0 * TMath::Pi()) * fSpread));
  320. fPRF->SetParameter(1, 0);
  321. fPRF->SetParameter(2, fSpread);
  322. return fPRF;
  323. }
  324. Bool_t MpdTpcDigitizerTask::isSubtrackInInwards(const TpcPoint *p1, const TpcPoint *p2) { //WHAT AM I DOING???
  325. const Float_t x1 = p1->GetX();
  326. const Float_t x2 = p2->GetX();
  327. const Float_t y1 = p1->GetY();
  328. const Float_t y2 = p2->GetY();
  329. const Float_t a = (y1 - y2) / (x1 - x2);
  330. const Float_t b = (y1 * x2 - x1 * y2) / (x2 - x1);
  331. const Float_t minR = fabs(b) / sqrt(a * a + 1);
  332. if (minR < r_min) //then check if minimal distance is between our points
  333. {
  334. const Float_t x = -a * b / (a * a + 1);
  335. const Float_t y = b / (a * a + 1);
  336. if ((x1 - x) * (x2 - x) < 0 && (y1 - y) * (y2 - y) < 0) {
  337. return kTRUE;
  338. }
  339. }
  340. return kFALSE;
  341. }
  342. void MpdTpcDigitizerTask::TpcProcessing(const TpcPoint* prePoint, const TpcPoint* curPoint, const UInt_t secID, const UInt_t iPoint, const UInt_t nPoints) {
  343. Float_t dE = 0.0; //energy loss
  344. UInt_t qTotal = 0; //sum of clusters charges (=sum of electrons between two TpcPoints)
  345. UInt_t qCluster = 0; //charge of cluster (= number of electrons)
  346. TLorentzVector curPointPos; // coordinates for current TpcPoint
  347. TLorentzVector prePointPos; // coordinates for previous TpcPoint
  348. TLorentzVector diffPointPos; // steps for clusters creation
  349. TVector3 diffuse; // vector of diffuse for every coordinates
  350. TVector3 distort; // vector of distortion for every coordinates
  351. TLorentzVector electronPos; // coordinates for created electrons
  352. TLorentzVector clustPos; // coordinates for created clusters
  353. Float_t driftl = 0.0; // length for drifting
  354. vector<UInt_t> clustArr; // vector of clusters between two TpcPoints
  355. Float_t localX = 0.0, localY = 0.0; //local coordinates of electron (sector coordinates)
  356. if ( fPrintDebugInfo && (iPoint % 1000 == 0) ) cout << UInt_t(iPoint * 1.0 / nPoints * 100.0) << " % of TPC points processed" << endl;
  357. // curPoint = (TpcPoint*) fMCPointArray->At(i);
  358. if (fOnlyPrimery == kTRUE) {
  359. MpdMCTrack* tr = (MpdMCTrack*) fMCTracksArray->At(curPoint->GetTrackID());
  360. if (tr->GetMotherId() != -1) return;
  361. }
  362. //check if hits are on the same track
  363. if (curPoint->GetTrackID() == prePoint->GetTrackID() && !isSubtrackInInwards(prePoint, curPoint)) {
  364. dE = curPoint->GetEnergyLoss() * 1E9; //convert from GeV to eV
  365. if (dE < 0) {
  366. Error("MpdTpcDigitizerTask::Exec", "Negative Energy loss!");
  367. return;
  368. }
  369. curPointPos.SetXYZT(curPoint->GetX(), curPoint->GetY(), curPoint->GetZ(), curPoint->GetTime());
  370. prePointPos.SetXYZT(prePoint->GetX(), prePoint->GetY(), prePoint->GetZ(), prePoint->GetTime());
  371. if ((curPointPos.T() < 0) || (prePointPos.T() < 0)) {
  372. Error("MpdTpcDigitizerTask::Exec", "Negative Time!");
  373. return;
  374. }
  375. diffPointPos = curPointPos - prePointPos; //differences between two points by coordinates
  376. diffPointPos *= (1 / diffPointPos.Vect().Mag()); //directional cosines //TODO! Do we need this??? Look at line #297
  377. qTotal = (UInt_t) floor(fabs(dE / fGas->W()));
  378. //while still charge not used-up distribute charge into next cluster
  379. if (fDistribute) {
  380. while (qTotal > 0) {
  381. //roll dice for next cluster
  382. qCluster = fGas->GetRandomCSUniform();
  383. if (qCluster > qTotal) qCluster = qTotal;
  384. qTotal -= qCluster;
  385. clustArr.push_back(qCluster);
  386. }// finish loop for cluster creation
  387. } else {
  388. clustArr.push_back(qTotal); // JUST FOR TEST. NO CLUSTER DISTRIBUTION!
  389. // clustArr.push_back(1); // JUST FOR TEST. NO CLUSTER DISTRIBUTION ONLY ONE ELECTRON IN CLUSTER!
  390. }
  391. diffPointPos *= (diffPointPos.Vect().Mag() / clustArr.size()); //now here are steps between clusters by coordinates TODO: correct distribution
  392. clustPos = prePointPos;
  393. for (UInt_t iClust = 0; iClust < clustArr.size(); ++iClust) {
  394. clustPos += diffPointPos;
  395. driftl = zCathode - fabs(clustPos.Z());
  396. for (UInt_t iEll = 0; iEll < clustArr.at(iClust); ++iEll) {
  397. //attachment
  398. if (fAttach)
  399. if (exp(-driftl * fGas->k()) < gRandom->Uniform()) continue; // FIXME
  400. //diffusion
  401. if (fDiffuse) {
  402. const Float_t sqrtDrift = sqrt(driftl);
  403. const Float_t sigmat = fGas->Dt() * sqrtDrift;
  404. const Float_t sigmal = fGas->Dl() * sqrtDrift;
  405. diffuse.SetXYZ(gRandom->Gaus(0, sigmat), gRandom->Gaus(0, sigmat), gRandom->Gaus(0, sigmal));
  406. }
  407. if (fDistort) {
  408. const Float_t dt = 1E-03; //time step [s]
  409. const Float_t mu = 4.23; //electron mobility [m^2 / s / V]
  410. const Float_t mu2 = mu * mu; //just square of mu
  411. 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
  412. TVector3 B(0.0, 0.0, 0.0); // vector of magnetic field components
  413. TVector3 v; // vector of current velocity components
  414. TVector3 posCur; // vector of current position components
  415. TVector3 posPre = clustPos.Vect(); // vector of previous position components
  416. TVector3 EBCross(0.0, 0.0, 0.0); // vector product of E and B vectors
  417. Bool_t inTpc = kTRUE;
  418. while (inTpc) {
  419. 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);
  420. EBCross = E.Cross(B);
  421. v = mu / (1 + mu2 * B.Mag2()) * (E - mu * EBCross + mu2 * B * (E * B));
  422. posCur = v * dt + posPre;
  423. // cout << "X = " << posCur.X() << " Y = " << posCur.Y() << " Z = " << posCur.Z() << " Vx = " << v.X() << " Vy = " << v.Y() << " Vz = " << v.Z() << endl;
  424. if ((posCur.Perp() > r_min + fSectHeight) || (posCur.Perp() < r_min) || (Abs(posCur.Z()) > zCathode)) inTpc = kFALSE;
  425. posPre = posCur;
  426. }
  427. distort.SetX(posCur.X() - clustPos.Vect().X());
  428. distort.SetY(posCur.Y() - clustPos.Vect().Y());
  429. distort.SetZ(0.0); //FIXME
  430. }
  431. electronPos.SetVect(clustPos.Vect() + diffuse + distort);
  432. electronPos.SetT(clustPos.T() + (zCathode - fabs(electronPos.Z())) / fGas->VDrift()); // Do we need to use clustPos.T() ???
  433. const Float_t phiStep = TwoPi() / nSectors * 2;
  434. const Float_t sectPhi = secID * phiStep;
  435. localY = electronPos.X() * Cos(sectPhi) + electronPos.Y() * Sin(sectPhi) - r_min; //converting from global to local coordinates
  436. localX = -electronPos.X() * Sin(sectPhi) + electronPos.Y() * Cos(sectPhi); //converting from global to local coordinates
  437. if ((localY < 0.0) || (Abs(localX) > fSector->GetMaxX()) || (localY > fSectHeight)) continue; //FIXME!!!
  438. const Float_t timeStep = (zCathode / nTimeBackets) / fGas->VDrift();
  439. const UInt_t curTimeID = (UInt_t) ((zCathode / fGas->VDrift() - electronPos.T()) / timeStep);
  440. if (curTimeID >= nTimeBackets) continue;
  441. if (fMakeQA) {
  442. fHisto->_hX_local->Fill(localX);
  443. fHisto->_hY_local->Fill(localY);
  444. fHisto->_hXY_local->Fill(localX, localY);
  445. fHisto->_hYZ_local->Fill(fabs(electronPos.Z()), localY);
  446. fHisto->_hXY_global->Fill(electronPos.X(), electronPos.Y());
  447. fHisto->_hRZ_global->Fill(electronPos.Z(), electronPos.Y());
  448. fHisto->_h3D_el->Fill(electronPos.X(), electronPos.Y(), electronPos.Z());
  449. fHisto->_hZ_local->Fill(fabs(electronPos.Z()));
  450. fHisto->_hX_global->Fill(electronPos.X());
  451. fHisto->_hY_global->Fill(electronPos.Y());
  452. fHisto->_hZ_global->Fill(electronPos.Z());
  453. fHisto->_hDiffuseXY->Fill(diffuse.X(), diffuse.Y());
  454. fHisto->_hDistortXY->Fill(distort.X(), distort.Y());
  455. }
  456. Int_t origin = prePoint->GetTrackID();
  457. PadResponse(localX, localY, curTimeID, origin, fDigits4dArray);
  458. }
  459. }
  460. }//end check for same track
  461. clustArr.clear();
  462. clustArr.resize(0);
  463. }
  464. void MpdTpcDigitizerTask::Finish() {
  465. cout << "Digitizer work time = " << ((Float_t)tAll) / CLOCKS_PER_SEC << endl;
  466. if (fMakeQA) {
  467. toDirectory("QA/TPC");
  468. Float_t digit = 0.0;
  469. UInt_t iPad_shifted = 0; //needed for correct drawing of fDigitsArray
  470. for (UInt_t iRows = 0; iRows < nRows; ++iRows) {
  471. for (UInt_t iPads = 0; iPads < fNumOfPadsInRow[iRows] * 2; ++iPads) {
  472. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[iRows];
  473. for (UInt_t iTime = 0; iTime < nTimeBackets; ++iTime) {
  474. digit = fDigits4dArray[iRows][iPads][iTime].signal;
  475. fHisto->_hXY_dig->Fill(iPad_shifted, iRows, digit);
  476. // fHisto->_hSect_dig->Fill(iSect, digit);
  477. fHisto->_hX_dig->Fill(iPad_shifted, digit);
  478. fHisto->_hY_dig->Fill(iRows, digit);
  479. fHisto->_hZ_dig->Fill(iTime, digit);
  480. fHisto->_h3D_dig->Fill(iPad_shifted, iRows, iTime, digit);
  481. if (digit > 1000.0) fHisto->_hADC_dig->Fill(digit);
  482. }
  483. }
  484. }
  485. for (UInt_t iRows = 0; iRows < nRows; ++iRows) {
  486. for (UInt_t iPads = 0; iPads < fNumOfPadsInRow[iRows] * 2; ++iPads) {
  487. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[iRows];
  488. for (UInt_t iTime = 0; iTime < nTimeBackets; ++iTime) {
  489. digit = fDigits4dArray[iRows][iPads][iTime].signal;
  490. //pad activity
  491. //if (digit > 1000.0) {
  492. // fHisto->_hXY_dig->Fill(iPad_shifted, iRows, 1.0);
  493. //}
  494. // fHisto->_hXY_dig->Fill(iPad_shifted, iRows, digit);
  495. fHisto->_h3D_dig->Fill(iPad_shifted, iRows, iTime, digit);
  496. }
  497. }
  498. }
  499. for (UInt_t iTime = 0; iTime < nTimeBackets; ++iTime) {
  500. for (UInt_t iPads = 0; iPads < fNumOfPadsInRow[1] * 2; ++iPads) {
  501. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[1];
  502. digit = fDigits4dArray[1][iPads][iTime].signal;
  503. fHisto->_hXT_dig_1->Fill(iPad_shifted, iTime, digit);
  504. }
  505. for (UInt_t iPads = 0; iPads < fNumOfPadsInRow[5] * 2; ++iPads) {
  506. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[5];
  507. digit = fDigits4dArray[5][iPads][iTime].signal;
  508. fHisto->_hXT_dig_5->Fill(iPad_shifted, iTime, digit);
  509. }
  510. for (UInt_t iPads = 0; iPads < fNumOfPadsInRow[10] * 2; ++iPads) {
  511. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[10];
  512. digit = fDigits4dArray[10][iPads][iTime].signal;
  513. fHisto->_hXT_dig_10->Fill(iPad_shifted, iTime, digit);
  514. }
  515. for (UInt_t iPads = 0; iPads < fNumOfPadsInRow[20] * 2; ++iPads) {
  516. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[20];
  517. digit = fDigits4dArray[20][iPads][iTime].signal;
  518. fHisto->_hXT_dig_20->Fill(iPad_shifted, iTime, digit);
  519. }
  520. for (UInt_t iPads = 0; iPads < fNumOfPadsInRow[40] * 2; ++iPads) {
  521. iPad_shifted = iPads + fNumOfPadsInRow[nRows - 1] - fNumOfPadsInRow[40];
  522. digit = fDigits4dArray[40][iPads][iTime].signal;
  523. fHisto->_hXT_dig_40->Fill(iPad_shifted, iTime, digit);
  524. }
  525. }
  526. fHisto->Write();
  527. gFile->cd();
  528. }
  529. }
  530. ClassImp(MpdTpcDigitizerTask)