MpdTpcDriftTask.cxx 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187
  1. //-----------------------------------------------------------
  2. // File and Version Information:
  3. // $Id$
  4. //
  5. // Description:
  6. // Implementation of class MpdTpcDriftTask
  7. // see MpdTpcDriftTask.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 "MpdTpcDriftTask.h"
  18. #include "MpdTpcClusterizerTask.h"
  19. // Collaborating Class Headers --------
  20. #include "TpcGas.h"
  21. #include "TpcPrimaryCluster.h"
  22. #include "TpcDriftedElectron.h"
  23. #include "FairRootManager.h"
  24. #include <TClonesArray.h>
  25. #include <TMath.h>
  26. #include <TRandom.h>
  27. #include <TSystem.h>
  28. // C/C++ Headers ----------------------
  29. #include <math.h>
  30. #include <iostream>
  31. // Class Member definitions -----------
  32. //__________________________________________________________________________
  33. MpdTpcDriftTask::MpdTpcDriftTask()
  34. : FairTask("TPC Drift"), fkNsec2(MpdTpcClusterizerTask::fgkNsec2), fPersistence(kFALSE),
  35. fAttach(kTRUE), fDiffuse(kTRUE), fDistort(kFALSE), fSec(0)
  36. {
  37. fPrimBranchName = "TpcPrimClus";
  38. std::string tpcGasFile = gSystem->Getenv("VMCWORKDIR");
  39. tpcGasFile += "/geometry/Ar-90_CH4-10.asc";
  40. fGas= new TpcGas(tpcGasFile, 130);
  41. //std::cout<<*_gas<<std::endl;
  42. }
  43. //__________________________________________________________________________
  44. MpdTpcDriftTask::~MpdTpcDriftTask()
  45. {
  46. delete fGas;
  47. }
  48. //__________________________________________________________________________
  49. InitStatus MpdTpcDriftTask::Init()
  50. {
  51. //Get ROOT Manager
  52. FairRootManager* ioman= FairRootManager::Instance();
  53. if(ioman==0)
  54. {
  55. Error("MpdTpcDriftTask::Init","RootManager not instantiated!");
  56. return kERROR;
  57. }
  58. // Get input collection
  59. //std::cout << fkNsec2 << std::endl;
  60. fPrimArray = new TClonesArray* [fkNsec2];
  61. for (Int_t i = 0; i < fkNsec2; ++i) {
  62. TString primBranchName = fPrimBranchName;
  63. primBranchName += i;
  64. fPrimArray[i] = (TClonesArray*) ioman->GetObject(primBranchName);
  65. if(fPrimArray[i] == 0x0) {
  66. Error("TpcDriftTask::Init","PrimaryElectron-array not found!");
  67. return kERROR;
  68. }
  69. }
  70. // create and register output array
  71. //SetPersistence();
  72. fDriftedArray = new TClonesArray* [fkNsec2];
  73. for (Int_t i = 0; i < fkNsec2; ++i) {
  74. TString branchName = "TpcDriftEl";
  75. branchName += i;
  76. fDriftedArray[i] = new TClonesArray("TpcDriftedElectron");
  77. ioman->Register(branchName, "Tpc", fDriftedArray[i], fPersistence);
  78. }
  79. // TODO: Get from geom!!!!
  80. fzGem = 125;
  81. return kSUCCESS;
  82. }
  83. //__________________________________________________________________________
  84. void MpdTpcDriftTask::Exec(Option_t* opt)
  85. {
  86. // Sector definition - should be taken somewhere else
  87. static const Int_t nSec = fkNsec2 / 2;
  88. static const Double_t phiSec = TMath::TwoPi() / nSec;
  89. if (TString(opt) == "") return;
  90. sscanf(opt, "%2d", &fSec);
  91. // Reset output Array
  92. if (fDriftedArray[fSec] == 0x0) Fatal("MpdTpcDriftTask::Exec)","No DriftedElectronArray");
  93. fDriftedArray[fSec]->Delete();
  94. //loop over incoming electrons
  95. Int_t nc = fPrimArray[fSec]->GetEntriesFast();
  96. for (Int_t ic = 0; ic < nc; ++ic) {
  97. TpcPrimaryCluster* pcl = (TpcPrimaryCluster*)fPrimArray[fSec]->UncheckedAt(ic);
  98. //create single electrons
  99. Int_t q = pcl->q();
  100. for(Int_t ie = 0; ie < q; ++ie) {
  101. Double_t z_length = pcl->z();
  102. Double_t driftl = 0;
  103. Double_t dx = 0, dy = 0, dt = 0;
  104. Int_t size = 0;
  105. if (z_length == 0) continue;
  106. // Calculate drift time
  107. if (z_length > 0) driftl = fzGem - z_length;
  108. else driftl = fzGem + z_length;
  109. if (driftl < 0) continue;
  110. // attachment
  111. if (fAttach) {
  112. if ( exp(-driftl * fGas->k()) < gRandom->Uniform() ) continue;
  113. }
  114. // diffusion
  115. if (fDiffuse) {
  116. Double_t sqrtDrift = sqrt(driftl);
  117. Double_t sigmat = fGas->Dt() * sqrtDrift;
  118. Double_t sigmal = fGas->Dl() * sqrtDrift;
  119. dt = (driftl+gRandom->Gaus(0,sigmal)) / fGas->VDrift();
  120. dx = gRandom->Gaus(0,sigmat);
  121. dy = gRandom->Gaus(0,sigmat);
  122. }
  123. // drift distortions
  124. if (fDistort) {
  125. // TODO: to be implemented
  126. }
  127. // Check if electron goes to another sector
  128. Double_t xend = pcl->x() + dx, yend = pcl->y() + dy;
  129. Double_t phi = TMath::ATan2 (yend, xend);
  130. Int_t isec = (Int_t) (phi / phiSec);
  131. if (phi < 0) {
  132. --isec;
  133. isec = nSec + isec;
  134. }
  135. if (fSec > nSec - 1) isec += nSec; // z < 0
  136. size = fDriftedArray[isec]->GetEntriesFast();
  137. new((*fDriftedArray[isec])[size]) TpcDriftedElectron(xend, yend,
  138. pcl->t()+dt, pcl);
  139. } // end loop over electrons
  140. } // end loop over clusters
  141. std::cout<<fDriftedArray[fSec]->GetEntriesFast()<<" electrons arriving at readout No. "
  142. << fSec << std::endl;
  143. return;
  144. }
  145. //__________________________________________________________________________
  146. void MpdTpcDriftTask::Clear(Option_t* opt)
  147. {
  148. TString option = opt;
  149. if (option == "") return;
  150. if (option.Contains("fill") && !fPersistence) return;
  151. Int_t isec = TString(option(4,2)).Atoi();
  152. //std::cout << " sec " << isec << " " << option.Data() << std::endl;
  153. TString name = "TpcDriftEl";
  154. name += isec;
  155. if (option.Contains("fill")) FairRootManager::Instance()->GetOutTree()->GetBranch(name)->Fill();
  156. else {
  157. fDriftedArray[isec]->Delete();
  158. fDriftedArray[isec]->Expand(1000); // shrink TClonesArray
  159. }
  160. }
  161. ClassImp(MpdTpcDriftTask)