CbmSttHitProducerReal.cxx 8.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262
  1. /////////////////////////////////////////////////////////////
  2. // CbmSttHitProducerReal
  3. //
  4. // Class for digitalization for STT1
  5. //
  6. // authors: Pablo Genova - Pavia University
  7. // Lia Lavezzi - Pavia University
  8. //
  9. /////////////////////////////////////////////////////////////
  10. #include "TClonesArray.h"
  11. #include "FairRootManager.h"
  12. #include "CbmSttHitProducerReal.h"
  13. #include "CbmSttHit.h"
  14. #include "CbmSttHitInfo.h"
  15. #include "CbmSttPoint.h"
  16. #include "TGeoManager.h"
  17. // #include "TGeoVolume.h"
  18. // #include "TGeoNode.h"
  19. // #include "TGeoMatrix.h"
  20. #include "TVector3.h"
  21. #include "TStraw.h"
  22. #include "TRandom.h"
  23. #include <iostream>
  24. #include "TMath.h"
  25. using std::cout;
  26. using std::cerr;
  27. using std::endl;
  28. // ----- Default constructor -------------------------------------------
  29. CbmSttHitProducerReal::CbmSttHitProducerReal() :
  30. FairTask("Ideal STT Hit Producer") { }
  31. // -------------------------------------------------------------------------
  32. // ----- Destructor ----------------------------------------------------
  33. CbmSttHitProducerReal::~CbmSttHitProducerReal() { }
  34. // -------------------------------------------------------------------------
  35. // ----- Public method Init --------------------------------------------
  36. InitStatus CbmSttHitProducerReal::Init() {
  37. // Get RootManager
  38. FairRootManager* ioman = FairRootManager::Instance();
  39. if ( ! ioman ) {
  40. cout << "-E- CbmSttHitProducerReal::Init: "
  41. << "RootManager not instantiated!" << endl;
  42. return kFATAL;
  43. }
  44. // Get input array
  45. fPointArray = (TClonesArray*) ioman->GetObject("STTPoint");
  46. if ( ! fPointArray ) {
  47. cout << "-W- CbmSttHitProducerReal::Init: "
  48. << "No STTPoint array!" << endl;
  49. return kERROR;
  50. }
  51. // Create and register output array
  52. fHitArray = new TClonesArray("CbmSttHit");
  53. ioman->Register("STTHit","STT",fHitArray,kTRUE);
  54. // Create and register output array
  55. fHitInfoArray = new TClonesArray("CbmSttHitInfo");
  56. ioman->Register("STTHitInfo", "STT", fHitInfoArray, kTRUE);
  57. // Geometry loading
  58. //TFile *tstfile=ioman->GetInFile();
  59. //TGeoManager *geoMan = (TGeoManager*) tstfile->Get("CBMGeom");
  60. TFile* f = (TFile*) ioman->GetInChain()->GetListOfFriends()->FindObject("testparams.root");
  61. if (!f)
  62. {
  63. cout<<"\n ---> ERROR: friend's list hasn't testparams.root";
  64. return kERROR;
  65. }
  66. f->Get("FairBaseParSet");
  67. if(!gGeoManager) // parameter file don't have "FairBaseParSet"
  68. {
  69. cout<<"\n ---> ERROR: parameter file don't have \"FairBaseParSet\" folder.";
  70. f->Close();
  71. return kERROR;
  72. }
  73. fVolumeArray = gGeoManager->GetListOfVolumes();
  74. cout << "-I- CbmSttHitProducerReal: Intialization successfull" << endl;
  75. return kSUCCESS;
  76. }
  77. // -------------------------------------------------------------------------
  78. // ----- Public method Exec --------------------------------------------
  79. void CbmSttHitProducerReal::Exec(Option_t* opt) {
  80. // Int_t evtn=0;
  81. // CbmSttPoint *ptemp =(CbmSttPoint*) fPointArray->At(0);
  82. // if(ptemp !=NULL) {
  83. // evtn=ptemp->GetEventID();
  84. // if(evtn%50==0)cout << "Event Number "<<evtn<<endl;
  85. // }
  86. // Reset output array
  87. if ( ! fHitArray ) Fatal("Exec", "No HitArray");
  88. fHitArray->Clear();
  89. fHitInfoArray->Clear();
  90. Int_t detID = 0; // detectorID
  91. TVector3 pos, dpos; // position and error vectors
  92. // Declare some variables
  93. CbmSttPoint* point = NULL;
  94. // Loop over SttPoints
  95. Int_t nPoints = fPointArray->GetEntriesFast();
  96. for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
  97. point = (CbmSttPoint*) fPointArray->At(iPoint);
  98. if (point == NULL) continue;
  99. double InOut[6];
  100. memset(InOut, 0, sizeof(InOut));
  101. InOut[0] = point->GetXInLocal();
  102. InOut[1] = point->GetYInLocal();
  103. InOut[2] = point->GetZInLocal();
  104. InOut[3] = point->GetXOutLocal();
  105. InOut[4] = point->GetYOutLocal();
  106. InOut[5] = point->GetZOutLocal();
  107. // single straw tube simulation -----------------------
  108. TStraw stt;
  109. //setting the single straw tube simulation constants
  110. // 3 options currently available:
  111. // TConst(tube radius (cm), gas pressure (bar), Ar%, CO2%)
  112. // stt.TConst(0.4, 1, 0.9, 0.1);
  113. // stt.TConst(0.5, 1, 0.9, 0.1);
  114. stt.TConst(0.5, 2, 0.8, 0.2);
  115. // wire positioning
  116. stt.PutWireXYZ(0., 0., -75., 0., 0., 75.);
  117. // get particle momentum
  118. TVector3 momentum(point->GetPxOut(),point->GetPyOut(),point->GetPzOut()); // GeV/c
  119. Double_t GeV=1.;
  120. // position in cm (already in cm); momentum in GeV (already in GeV); mass in GeV (already in GeV)
  121. // drift time calculation
  122. Double_t pulset = stt.PartToTime(point->GetMass()/GeV, momentum.Mag()/GeV, InOut);
  123. // simulated radius (cm)
  124. double radius = stt.TimnsToDiscm(pulset);
  125. if(radius < 0.) radius = 0.; // CHECK
  126. // true radius (cm)
  127. double true_rad = stt.TrueDist(InOut);
  128. // dE calculation
  129. // double depCharge = stt.PartToADC();
  130. // dE/dx calculation
  131. //TVector3 diff3=point1->Get
  132. //double distance =diff3.Mag(); //
  133. //double dedx = 999;
  134. //if (distance != 0) dedx = depCharge/(1000000 * distance); // in arbitrary units
  135. // stt2: detID, pos, dpos, index come from --------------
  136. // stt2 (FairHit):
  137. Double_t closestDistanceError = 0.0150; // per adesso (stessa che in Ideal:
  138. // radialResolution = 0.0150)
  139. TVector3 position(point->GetX(), point->GetY(), point->GetZ());
  140. // ----------------
  141. // stt2, ma cancellati in stt1 (controlla: in stt2 la posizione dell' hit non
  142. // corrisponde al centro del tubo (xcentro, ycentro, 35.), ma per il Real deve
  143. // essere cosi' ??perche' in stt2 non e' cosi'??
  144. // TVector3 posInLocal(point->GetXInLocal(), point->GetYInLocal(), point->GetZInLocal());
  145. // TVector3 posOutLocal(point->GetXOutLocal(), point->GetYOutLocal(), point->GetZOutLocal());
  146. // Double_t zpos = position.Z() + ((posOutLocal.Z() + posInLocal.Z()) / 2.);
  147. // Double_t zposError;
  148. // FoldZPosWithResolution(zpos, zposError, posInLocal, posOutLocal);
  149. // pos.SetXYZ(position.X(), position.Y(), zpos); // <--- stt2
  150. // ----------------
  151. pos.SetXYZ(position.X(), position.Y(), position.Z()); // <--- stt1
  152. // dpos.SetXYZ(innerStrawDiameter / 2., innerStrawDiameter / 2., GetLongitudinalResolution(position.Z()));
  153. dpos.SetXYZ(0.5, 0.5, 3.); // per adesso (stessi che in Ideal:
  154. // innerStrawDiameter/2 = 0.5,
  155. // longitudinalResolution = 3.)
  156. //----- end stt2 ------------------------------------------
  157. // wire direction from stt2 -------------------------------
  158. TVector3 wireDirection(point->GetXWireDirection(), point->GetYWireDirection(), point->GetZWireDirection());
  159. // = 0, 0, 1 if only axias tubes
  160. // --------------------------------------------------------
  161. // create hit
  162. AddHit(detID, pos, dpos, iPoint, point->GetTrackID(), pulset, radius, true_rad, closestDistanceError, wireDirection);
  163. AddHitInfo(0, 0, point->GetTrackID(), iPoint, 0, kFALSE);
  164. }// Loop over MCPoints
  165. // Event summary
  166. cout << "-I- CbmSttHitProducerReal: " << nPoints << " SttPoints, "
  167. << nPoints << " Hits created." << endl;
  168. }
  169. // -------------------------------------------------------------------------
  170. void CbmSttHitProducerReal::FoldZPosWithResolution(Double_t &zpos, Double_t &zposError,
  171. TVector3 localInPos, TVector3 localOutPos)
  172. {
  173. Double_t
  174. zPosInStrawFrame = (localOutPos.Z() - localInPos.Z()) / 2.;
  175. // zposError = gRandom->Gaus(0., GetLongitudinalResolution(zPosInStrawFrame));
  176. zposError = gRandom->Gaus(0., 3.); // per adesso (stesso che in Ideal:
  177. // longitudinalResolution = 3.)
  178. zpos += zposError;
  179. }
  180. // ----- Private method AddHit --------------------------------------------
  181. CbmSttHit* CbmSttHitProducerReal::AddHit(Int_t detID, TVector3& pos, TVector3& dpos, Int_t iPoint, Int_t trackID, Double_t p, Double_t rsim, Double_t rtrue, Double_t closestDistanceError, TVector3 wireDirection){
  182. // see CbmSttHit for hit description
  183. TClonesArray& clref = *fHitArray;
  184. Int_t size = clref.GetEntriesFast();
  185. //cout << "-I- CbmSttHitProducerReal: Adding Hit: track"<<trackID<<" event: "<<eventID<<" pulse = " << p << ", rsim = " << rsim
  186. // << ", rtrue = " << rtrue <<" name "<<nam<<"center.X "<<center.X()<< endl;
  187. return new(clref[size]) CbmSttHit(detID, pos, dpos, iPoint, trackID, p, rsim, rtrue, closestDistanceError, wireDirection);
  188. }
  189. // ----
  190. // ----- Private method AddHitInfo --------------------------------------------
  191. CbmSttHitInfo* CbmSttHitProducerReal::AddHitInfo(Int_t fileNumber, Int_t eventNumber, Int_t trackID, Int_t pointID, Int_t nMerged, Bool_t isFake){
  192. // see CbmSttHitInfo for hit description
  193. TClonesArray& clref = *fHitInfoArray;
  194. Int_t size = clref.GetEntriesFast();
  195. return new(clref[size]) CbmSttHitInfo(fileNumber, eventNumber, trackID, pointID, nMerged, isFake);
  196. }
  197. ClassImp(CbmSttHitProducerReal)