MpdTpcClusterizerTask.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class MpdTpcClusterizerTask
  7. // see MpdTpcClusterizerTask.h for details
  8. //
  9. // Environment:
  10. // Software developed for the MPD Detector at NICA.
  11. //
  12. // Author List:
  13. // Alexandr Zinchenko LHEP, JINR, Dubna - adapted for MPD from PANDARoot
  14. //
  15. //-----------------------------------------------------------
  16. // This Class' Header ------------------
  17. #include "MpdTpcClusterizerTask.h"
  18. // Collaborating Class Headers --------
  19. #include "FairRootManager.h"
  20. #include "TpcPoint.h"
  21. #include "TpcGas.h"
  22. #include "TpcPrimaryCluster.h"
  23. #include "LinearInterpolPolicy.h"
  24. #include "TClonesArray.h"
  25. #include "TFile.h"
  26. #include "TMath.h"
  27. #include "TRandom.h"
  28. #include "TSystem.h"
  29. #include "TTree.h"
  30. // C/C++ Headers ----------------------
  31. #include <iostream>
  32. #include <math.h>
  33. #include <vector>
  34. using namespace std;
  35. // Class Member definitions -----------
  36. const Int_t MpdTpcClusterizerTask::fgkNsec2 = 24;
  37. //__________________________________________________________________________
  38. MpdTpcClusterizerTask::MpdTpcClusterizerTask()
  39. : FairTask("TPC Clusterizer"), fPersistence(kFALSE)
  40. {
  41. fPointBranchName = "TpcPoint";
  42. std::string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  43. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  44. fGas= new TpcGas(tpcGasFile, 130);
  45. std::cout<<*fGas<<std::endl;
  46. }
  47. //__________________________________________________________________________
  48. MpdTpcClusterizerTask::~MpdTpcClusterizerTask()
  49. {
  50. delete fGas;
  51. }
  52. //__________________________________________________________________________
  53. void MpdTpcClusterizerTask::FinishTask()
  54. {
  55. for (Int_t i = 0; i < fgkNsec2; ++i) fPrimArray[i]->Delete();
  56. }
  57. //__________________________________________________________________________
  58. InitStatus MpdTpcClusterizerTask::Init()
  59. {
  60. //Get ROOT Manager
  61. FairRootManager* ioman= FairRootManager::Instance();
  62. if(ioman==0)
  63. {
  64. Error("MpdTpcClusterizerTask::Init","RootManager not instantiated!");
  65. return kERROR;
  66. }
  67. // Get input collection
  68. fPointArray=(TClonesArray*) ioman->GetObject(fPointBranchName);
  69. if(fPointArray==0)
  70. {
  71. Error("TpcClusterizerTask::Init","Point-array not found!");
  72. return kERROR;
  73. }
  74. // create and register output array
  75. //SetPersistence();
  76. fPrimArray = new TClonesArray* [fgkNsec2];
  77. TString name0 = "TpcPrimClus";
  78. for (Int_t i = 0; i < fgkNsec2; ++i) {
  79. fPrimArray[i] = new TClonesArray("TpcPrimaryCluster");
  80. TString name = "TpcPrimClus";
  81. name += i;
  82. ioman->Register(name, "Tpc", fPrimArray[i], fPersistence);
  83. }
  84. return kSUCCESS;
  85. }
  86. //__________________________________________________________________________
  87. void MpdTpcClusterizerTask::Exec(Option_t* opt)
  88. {
  89. const Int_t nSec = fgkNsec2 / 2; // number of TPC readout sectors
  90. // Reset output Array
  91. if (fPrimArray == 0) Fatal("MpdTpcClusterizerTask::Exec)","No PrimClusterArray");
  92. for (Int_t i = 0; i < fgkNsec2; ++i) fPrimArray[i]->Delete();
  93. Int_t np = fPointArray->GetEntriesFast();
  94. if (np < 2){
  95. Warning("MpdTpcClusterizerTask::Exec","Not enough Hits in Tpc for Digitization (<2)");
  96. return;
  97. }
  98. TpcPoint *theLastPoint = (TpcPoint*)fPointArray->UncheckedAt(0), *point = theLastPoint;
  99. vector<pair<Int_t,Int_t> > sectors[fgkNsec2];
  100. // Sector definition - should be taken somewhere else
  101. const Double_t phiSec = TMath::TwoPi() / nSec, dDist = 1.2;
  102. // Find segments of tracks with the same ID inside one sector and one half
  103. // and not broken, i.e. without reentrance (short enough distance between consecutive points)
  104. Int_t curpointTrackID, ip0 = 0;
  105. curpointTrackID = point->GetTrackID();
  106. Double_t phi = TMath::ATan2 (point->GetY(), point->GetX());
  107. Int_t isec0 = (Int_t) (phi / phiSec);
  108. if (phi < 0) {
  109. --isec0;
  110. isec0 = nSec + isec0;
  111. }
  112. if (point->GetZ() < 0) isec0 += nSec;
  113. TVector3 pos0, pos;
  114. point->Position(pos0);
  115. for (Int_t ip = 1; ip < np; ++ip) {
  116. point = (TpcPoint*) fPointArray->UncheckedAt(ip);
  117. phi = TMath::ATan2 (point->GetY(), point->GetX());
  118. Int_t isec = (Int_t) (phi / phiSec);
  119. if (phi < 0) {
  120. --isec;
  121. isec = nSec + isec;
  122. }
  123. if (point->GetZ() < 0) isec += nSec;
  124. point->Position(pos);
  125. TVector3 dist = pos - pos0;
  126. if (isec != isec0 || point->GetTrackID() != curpointTrackID || dist.Mag() > dDist) {
  127. sectors[isec0].push_back(pair<Int_t,Int_t>(ip0,ip));
  128. ip0 = ip;
  129. isec0 = isec;
  130. curpointTrackID = point->GetTrackID();
  131. }
  132. pos0 = pos;
  133. }
  134. sectors[isec0].push_back(pair<Int_t,Int_t>(ip0,np));
  135. // Do clusterization sector-by-sector
  136. for (Int_t isec = 0; isec < fgkNsec2; ++isec) {
  137. //for (Int_t isec = 0; isec < 12; ++isec) {
  138. Int_t icluster = 0;
  139. // Loop over sectors
  140. Int_t nSeg = sectors[isec].size();
  141. for (Int_t iseg = 0; iseg < nSeg; ++iseg) {
  142. // Loop over track segments
  143. Int_t ibeg = sectors[isec][iseg].first, iend = sectors[isec][iseg].second;
  144. // Create "virtual" point - back-extrapolation from the first point
  145. theLastPoint = (TpcPoint*)fPointArray->UncheckedAt(ibeg);
  146. TpcPoint virt = *theLastPoint;
  147. theLastPoint->Momentum(pos0);
  148. theLastPoint->Position(pos);
  149. Double_t ppp = pos0.Mag();
  150. if (ppp > 1.e-6) {
  151. ppp = theLastPoint->GetStep() / ppp;
  152. for (Int_t i = 0; i < 3; ++i) pos[i] -= pos0[i] * ppp;
  153. } else ++ibeg;
  154. virt.SetPosition(pos);
  155. theLastPoint = &virt;
  156. for (Int_t ip = ibeg; ip < iend; ++ip) {
  157. // Loop over points
  158. point = (TpcPoint*) fPointArray->UncheckedAt(ip);
  159. //point->Print();
  160. //check if hits belong to the same track
  161. curpointTrackID = point->GetTrackID();
  162. if (curpointTrackID != theLastPoint->GetTrackID()) {
  163. cout << " -F- TpcClusterizerTask::Exec - Different track IDs " << curpointTrackID
  164. << " " << theLastPoint->GetTrackID() << endl;
  165. exit(0);
  166. }
  167. Double_t dE = point->GetEnergyLoss() * 1E9; //convert from GeV to eV
  168. //Step 0: calculate the overall ammount of charge, produced
  169. if (dE < 0) {
  170. Error("TpcClusterizerTask::Exec","Note: particle:: negative Energy loss!");
  171. continue;
  172. }
  173. Int_t qTotal = (Int_t) TMath::Abs(dE / fGas->W());
  174. Int_t qCluster = 0;
  175. Int_t nCluster = 0;
  176. //Step 1: Create Clusters
  177. Double_t pointTime = point->GetTime();
  178. assert(pointTime >= 0);
  179. // while still charge not used-up distribute charge into next cluster
  180. Int_t size = fPrimArray[isec]->GetEntriesFast();
  181. while (qTotal > 0) {
  182. //roll dice for next clustersize
  183. qCluster = fGas->GetRandomCS (gRandom->Uniform());
  184. if (qCluster > qTotal) qCluster = qTotal;
  185. qTotal -= qCluster;
  186. // create cluster
  187. new((*fPrimArray[isec])[size++]) TpcPrimaryCluster(pointTime, qCluster, TVector3(0,0,0),
  188. curpointTrackID, ip);
  189. ++nCluster;
  190. } // finish loop for cluster creation
  191. //Step 2: Distribute Clusters along track segment
  192. LinearInterpolPolicy().Interpolate(theLastPoint,point,fPrimArray[isec],icluster,nCluster);
  193. icluster += nCluster;
  194. theLastPoint = point;
  195. } // for (Int ip = ibeg; ip < iend;
  196. } // finish loop over segments
  197. std::cout<<"TpcClusterizer:: "<<fPrimArray[isec]->GetEntriesFast()<<" clusters created"<<std::endl;
  198. TString option = "";
  199. option += isec;
  200. ExecuteTasks(option);
  201. //cout << FairRootManager::Instance()->GetObject(name) << " " << FairRootManager::Instance()->GetObject(name)->ClassName() << endl;
  202. //cout << FairRootManager::Instance()->GetOutTree()->GetName() << " " << FairRootManager::Instance()->GetOutTree()->GetBranch(name) << endl;
  203. CleanTasks();
  204. // Reset data in all subtasks for selected sectors
  205. ClearData(isec);
  206. } // finish loop over sectors
  207. return;
  208. }
  209. //__________________________________________________________________________
  210. void MpdTpcClusterizerTask::ClearData(Int_t isec)
  211. {
  212. // Save and clear data structures
  213. const Int_t nSec = fgkNsec2 / 2; // number of TPC readout sectors
  214. TString option, name;
  215. TTask *task = 0x0;
  216. if (isec > 1 && isec < nSec || isec > nSec + 1 && isec < fgkNsec2) {
  217. Int_t saveSec = isec - 1;
  218. TIter next(fTasks);
  219. option = "fill";
  220. option += saveSec;
  221. while((task=(TTask*)next())) task->Clear(option);
  222. name = "TpcPrimClus";
  223. name += saveSec;
  224. if (fPersistence) FairRootManager::Instance()->GetOutTree()->GetBranch(name)->Fill();
  225. if (isec > 2 && isec < nSec - 1 || isec > nSec + 2) {
  226. fPrimArray[saveSec]->Delete();
  227. fPrimArray[saveSec]->Expand(1000); // shrink TClonesArray
  228. }
  229. option = "clea";
  230. option += saveSec;
  231. next.Reset();
  232. while((task=(TTask*)next())) task->Clear(option);
  233. }
  234. if (isec == nSec-1) {
  235. // Process remaining sectors at Z > 0 slightly differently
  236. if (fPersistence) {
  237. for (Int_t i = 0; i < nSec; i+=(nSec-1)) {
  238. name = "TpcPrimClus";
  239. name += i;
  240. FairRootManager::Instance()->GetOutTree()->GetBranch(name)->Fill();
  241. }
  242. }
  243. TIter next(fTasks);
  244. option = "fill0";
  245. next.Reset();
  246. while((task=(TTask*)next())) task->Clear(option);
  247. option = "clea0";
  248. next.Reset();
  249. while((task=(TTask*)next())) task->Clear(option);
  250. option = "fill";
  251. option += isec;
  252. next.Reset();
  253. while((task=(TTask*)next())) task->Clear(option);
  254. option = "clea";
  255. option += isec;
  256. next.Reset();
  257. while((task=(TTask*)next())) task->Clear(option);
  258. for (Int_t i = 0; i < 2; ++i) {
  259. fPrimArray[i]->Delete();
  260. fPrimArray[i]->Expand(1000); // shrink TClonesArray
  261. }
  262. for (Int_t i = nSec-2; i < nSec; ++i) {
  263. fPrimArray[i]->Delete();
  264. fPrimArray[i]->Expand(1000); // shrink TClonesArray
  265. }
  266. }
  267. if (isec == fgkNsec2 - 1) {
  268. // Process remaining sectors at Z < 0
  269. if (fPersistence) {
  270. for (Int_t i = nSec; i < fgkNsec2; i+=(nSec-1)) {
  271. name = "TpcPrimClus";
  272. name += i;
  273. FairRootManager::Instance()->GetOutTree()->GetBranch(name)->Fill();
  274. }
  275. }
  276. TIter next(fTasks);
  277. option = "fill";
  278. option += nSec;
  279. next.Reset();
  280. while((task=(TTask*)next())) task->Clear(option);
  281. option = "clea";
  282. option += nSec;
  283. next.Reset();
  284. while((task=(TTask*)next())) task->Clear(option);
  285. option = "fill";
  286. option += isec;
  287. next.Reset();
  288. while((task=(TTask*)next())) task->Clear(option);
  289. option = "clea";
  290. option += isec;
  291. next.Reset();
  292. while((task=(TTask*)next())) task->Clear(option);
  293. for (Int_t i = nSec; i < nSec+2; ++i) {
  294. fPrimArray[i]->Delete();
  295. fPrimArray[i]->Expand(1000); // shrink TClonesArray
  296. }
  297. for (Int_t i = isec-1; i <= isec; ++i) {
  298. fPrimArray[i]->Delete();
  299. fPrimArray[i]->Expand(1000); // shrink TClonesArray
  300. }
  301. }
  302. }
  303. ClassImp(MpdTpcClusterizerTask)