MpdZdc.cxx 26 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636
  1. /*************************************************************************************
  2. *
  3. * Class MpdZdc
  4. *
  5. * Adopted for MPD by: Elena Litvinenko (EL)
  6. * e-mail: litvin@nf.jinr.ru
  7. * Version: 8-Apr-2008
  8. * Last update: 22-Feb-2012 (EL)
  9. *
  10. ************************************************************************************/
  11. #include <iostream>
  12. #include "TClonesArray.h"
  13. #include "TGeoMCGeometry.h"
  14. #include "TGeoManager.h"
  15. #include "TLorentzVector.h"
  16. #include "TParticle.h"
  17. #include "TVirtualMC.h"
  18. #include "TGeoArb8.h"
  19. #include "FairGeoInterface.h"
  20. #include "FairGeoLoader.h"
  21. #include "FairGeoNode.h"
  22. #include "MpdZdcGeo.h"
  23. #include "FairGeoRootBuilder.h"
  24. #include "MpdStack.h"
  25. #include "MpdZdc.h"
  26. #include "MpdZdcPoint.h"
  27. #include "FairRootManager.h"
  28. #include "FairVolume.h"
  29. // add on for debug
  30. //#include "FairGeoG3Builder.h"
  31. #include "FairRuntimeDb.h"
  32. #include "TObjArray.h"
  33. #include "FairRun.h"
  34. #include "TParticlePDG.h"
  35. // ----- Default constructor -------------------------------------------
  36. MpdZdc::MpdZdc() {
  37. fZdcCollection = new TClonesArray("MpdZdcPoint");
  38. volDetector = 0;
  39. fPosIndex = 0;
  40. // fpreflag = 0;
  41. //fpostflag = 0;
  42. fEventID=-1;
  43. fHitNb=0;
  44. fVSCVolId=0;
  45. fVSCHVolId=0;
  46. fVerboseLevel = 1;
  47. }
  48. // -------------------------------------------------------------------------
  49. // ----- Standard constructor ------------------------------------------
  50. MpdZdc::MpdZdc(const char* name, Bool_t active)
  51. : FairDetector(name, active) {
  52. fZdcCollection = new TClonesArray("MpdZdcPoint");
  53. fPosIndex = 0;
  54. volDetector = 0;
  55. //fpreflag = 0;
  56. //fpostflag = 0;
  57. fEventID=-1;
  58. fHitNb=0;
  59. fVSCVolId=0;
  60. fVSCHVolId=0;
  61. fVerboseLevel = 1;
  62. }
  63. // -------------------------------------------------------------------------
  64. // ----- Destructor ----------------------------------------------------
  65. MpdZdc::~MpdZdc() {
  66. if (fZdcCollection) {
  67. fZdcCollection->Delete();
  68. delete fZdcCollection;
  69. }
  70. }
  71. // -------------------------------------------------------------------------
  72. // ----- Public method Intialize ---------------------------------------
  73. void MpdZdc::Initialize() {
  74. // Init function
  75. FairDetector::Initialize();
  76. //FairRun* sim = FairRun::Instance();
  77. //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
  78. fVSCVolId = gMC->VolId("zdc01s");
  79. fVSCHVolId = gMC->VolId("ScH");
  80. }
  81. // -------------------------------------------------------------------------
  82. void MpdZdc::BeginEvent(){
  83. // Begin of the event
  84. }
  85. //_____________________________________________________________________________
  86. MpdZdcPoint* MpdZdc::GetHit(Int_t i) const
  87. {
  88. // Returns the hit for the specified layer.
  89. // ---
  90. //return (Ex03CalorHit*)fCalCollection->At(i);
  91. return (MpdZdcPoint*)fZdcCollection->At(i);
  92. }
  93. //_____________________________________________________________________________
  94. MpdZdcPoint* MpdZdc::GetHit(Int_t vsc, Int_t mod, Int_t zdc) const
  95. {
  96. // Returns the hit for the specified vsc and module0.
  97. // ---
  98. MpdZdcPoint *hit;
  99. Int_t nofHits = fZdcCollection->GetEntriesFast();
  100. //cout <<"NAVetoHit: nofHits " <<nofHits << " vsc " << vsc << " mod " << mod <<endl;
  101. //for (Int_t i=0; i<nofHits; i++) {
  102. // hit = GetHit(i);
  103. // cout <<"GetHit method-> " <<i <<" " <<hit->GetCopy() <<" " <<hit->GetCopyMother() <<" " <<hit->GetCopyZdc()<<" " <<hit->GetEnergyLoss() <<" " <<hit->GetZ() <<endl;
  104. //}
  105. for (Int_t i=0; i<nofHits; i++) {
  106. hit = GetHit(i);
  107. //cout <<"fEdep " <<i <<" " <<hit->GetEdep() <<endl;
  108. //int iVSCId = hit->GetVSCId();
  109. //int iMODId = hit->GetMODId();
  110. //if(iVSCId == vsc && iMODId == mod)
  111. if(hit->GetCopy() == vsc && hit->GetCopyMother() == mod && hit->GetCopyZdc() ==zdc)
  112. return hit;
  113. }
  114. return 0;
  115. }
  116. /*
  117. //_____________________________________________________________________________
  118. MpdZdcPoint* MpdZdc::GetHitH(Int_t vsch, Int_t modh) const
  119. {
  120. // Returns the hit for the specified vsc and module0.
  121. // ---
  122. MpdZdcPoint *hit;
  123. Int_t nofHits = fZdcCollection->GetEntriesFast();
  124. //cout <<"NAVetoHit: nofHits " <<nofHits << " vsc " << vsc << " mod " << mod <<endl;
  125. //for (Int_t i=0; i<nofHits; i++) {
  126. //hit = GetHit(i);
  127. //if(hit->GetCopyMotherH()<0 || hit->GetCopyH()<0) cout <<"GetHitH -59 -> " <<i <<" " <<hit->GetCopyMotherH() <<" " <<hit->GetEnergyLoss() <<" " <<hit->GetCopyH() <<endl;
  128. //}
  129. for (Int_t i=0; i<nofHits; i++) {
  130. hit = GetHit(i);
  131. //cout <<"fEdep " <<i <<" " <<hit->GetEdep() <<endl;
  132. //int iVSCId = hit->GetVSCId();
  133. //int iMODId = hit->GetMODId();
  134. //if(iVSCId == vsc && iMODId == mod)
  135. //cout <<"GetHitH " <<i <<" " <<hit->GetCopyMotherH() <<" " <<hit->GetEnergyLoss() <<" " <<hit->GetCopyH() <<endl;
  136. //if(hit->GetCopyH() == vsch && hit->GetCopyMotherH() == modh)
  137. if(hit->GetCopy() == vsch && hit->GetCopyMother() == modh)
  138. return hit;
  139. }
  140. return 0;
  141. }
  142. */
  143. // ----- Public method ProcessHits --------------------------------------
  144. Bool_t MpdZdc::ProcessHits(FairVolume* vol) {
  145. // if (TMath::Abs(gMC->TrackCharge()) <= 0) return kFALSE;
  146. //copyNoVTYVEC,copyNoVMOD,copyNoVZDC;
  147. //,copyNoVTYVECH,copyNoVMODH,copyNoVZDCH; copyNoVSCCom
  148. Int_t copyNoVSC{0}, copyNoVSCH{0}, copyNoVTYVECCom{0}, copyNoVMODCom{0}, copyNoVZDCCom{0};
  149. Int_t ivol;
  150. TLorentzVector tPos1, tMom1;
  151. TLorentzVector tPos, tMom;
  152. //Int_t copyNo;
  153. //Int_t ivol1;
  154. Int_t module0;//, module_h; //module0, module0 with hole
  155. Int_t slice;//, slice_h;
  156. Int_t zdc;//,zdc_h;
  157. Double_t time=0;
  158. Double_t length =0;
  159. //TParticle* part;
  160. //Double_t charge;
  161. Double_t QCF=1; //quenching for Birk
  162. Double_t BirkConst = 12.6; //0.126 mm/MeV for polystyrene
  163. //0.126 *(0.1/0.001) = 12.6 cm/GeV
  164. //(0.126 mm/MeV - from Wikipedia, 0.07943mm/MeV ô– in Geant4,G4EmSaturation.cc (push))
  165. //#define EDEBUG
  166. #ifdef EDEBUG
  167. static Int_t lEDEBUGcounter=0;
  168. if (lEDEBUGcounter<1)
  169. std::cout << "EDEBUG-- MpdZdc::ProcessHits: entered" << gMC->CurrentVolPath() << endl;
  170. #endif
  171. if (gMC->CurrentVolID(copyNoVSC) != fVSCVolId &&
  172. gMC->CurrentVolID(copyNoVSCH) != fVSCHVolId) {
  173. return kFALSE;
  174. }
  175. /*
  176. if (gMC->CurrentVolID(copyNoVSC) == fVSCVolId) {
  177. gMC->CurrentVolOffID(1, slice);
  178. gMC->CurrentVolOffID(2, module0);
  179. gMC->CurrentVolOffID(3, zdc);
  180. //copyNo = (copyNoVMOD-1)*60 + copyNoVTYVEC; //1-6480
  181. copyNoVTYVECCom = slice; copyNoVMODCom = module0; copyNoVZDCCom = zdc;
  182. }
  183. if (gMC->CurrentVolID(copyNoVSCH) == fVSCHVolId) {
  184. gMC->CurrentVolOffID(1, slice);
  185. gMC->CurrentVolOffID(2, module0);
  186. gMC->CurrentVolOffID(3, zdc);
  187. copyNoVTYVECCom = slice; copyNoVMODCom = module0; copyNoVZDCCom = zdc;
  188. }
  189. */
  190. ivol = vol->getMCid();
  191. if (gMC->CurrentVolID(copyNoVSC) == fVSCVolId || gMC->CurrentVolID(copyNoVSCH) == fVSCHVolId) {
  192. gMC->CurrentVolOffID(1, slice);
  193. gMC->CurrentVolOffID(2, module0);
  194. gMC->CurrentVolOffID(3, zdc);
  195. //copyNo = (copyNoVMOD-1)*60 + copyNoVTYVEC; //1-6480
  196. copyNoVTYVECCom = slice;
  197. copyNoVMODCom = module0;
  198. copyNoVZDCCom = zdc;
  199. }
  200. //cout <<"BEGIN copyNoVZDCCom " <<copyNoVZDCCom <<endl;
  201. if (gMC->IsTrackEntering()) {
  202. ResetParameters();
  203. fELoss = 0.;
  204. time = 0.;
  205. length = 0.;
  206. gMC->TrackPosition(tPos);
  207. gMC->TrackMomentum(tMom);
  208. #ifdef EDEBUG
  209. gMC->TrackPosition(tPos1);
  210. gMC->TrackMomentum(tMom1);
  211. #endif
  212. /*
  213. //Int_t copyNo;
  214. ivol1 = gMC->CurrentVolID(copyNo);
  215. gMC->CurrentVolOffID(2, module0);
  216. gMC->CurrentVolOffID(1, slice);
  217. */
  218. }//if (gMC->IsTrackEntering())
  219. if ( gMC->IsTrackInside()) {
  220. //if (gMC->Edep() <= 0. ) return kFALSE;
  221. //if (gMC->Edep() > 0. ) fELoss += gMC->Edep();
  222. ////if (gMC->Edep() <= 1.0e-20 ) return kFALSE;
  223. ////if (gMC->Edep() > 1.0e-20 ) fELoss += gMC->Edep();
  224. gMC->TrackPosition(tPos);
  225. gMC->TrackMomentum(tMom);
  226. length += gMC->TrackStep();
  227. //fELoss +=gMC->Edep();
  228. //Birk corrections
  229. QCF = 1.+(BirkConst/gMC->TrackStep())*gMC->Edep();
  230. fELoss +=(gMC->Edep())/QCF;
  231. time += gMC->TrackTime() * 1.0e09;
  232. if ( gMC->IsTrackStop() || gMC->IsTrackDisappeared() ) {
  233. // part = gMC->GetStack()->GetCurrentTrack();
  234. //charge = part->GetPDG()->Charge() / 3. ;
  235. // Create MpdZdcPoint
  236. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  237. //time = gMC->TrackTime() * 1.0e09;
  238. //length = gMC->TrackLength();
  239. //gMC->TrackPosition(tPos);
  240. //gMC->TrackMomentum(tMom);
  241. //gMC->CurrentVolOffID(2, module0);
  242. //gMC->CurrentVolOffID(1, slice);
  243. if(fELoss>0) {
  244. //std::cout << "INSIDE MpdZdc::ProcessHits: TrackID:" <<part->GetPdgCode() <<" " << fTrackID << " " <<fELoss <<" " << gMC->CurrentVolPath() << " " << tPos.Z() <<" " <<tMom.Pz() <<" " << ivol <<" " <<gMC->CurrentVolOffName(2) << " " <<gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<std::endl;
  245. //cout <<"CHECK INSIDE ivol,copyNo,iCell " <<part->GetPdgCode() <<" " << fTrackID <<" " <<ivol <<" " <<copyNoVZDCCom <<" " <<copyNoVMODCom <<" " <<copyNoVTYVECCom <<" " <<tPos.Z() <<" " <<tMom.Pz() <<" " <<fELoss <<endl;
  246. /*
  247. if(copyNoVTYVECCom==slice && copyNoVMODCom==module0 && copyNoVZDCCom==zdc) {//module0
  248. if ( !GetHit(slice,module0,zdc) ) {
  249. //cout <<"INSIDE NO GETHIT 1-> " <<slice <<" " <<module0 <<" " <<zdc <<endl;
  250. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, -59, -59, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  251. new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  252. //cout <<"INSIDE NO GETHIT 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  253. }
  254. else {
  255. //cout <<"INSIDE GETHIT 1-> " <<slice <<" " <<module0 <<" " <<zdc <<endl;
  256. GetHit(slice,module0,zdc)->AddVSC(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  257. //cout <<"INSIDE GETHIT 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  258. }
  259. }//if(copyNoVTYVECCom==copyNoVTYVEC
  260. */
  261. if(copyNoVTYVECCom==slice && copyNoVMODCom==module0 && copyNoVZDCCom==zdc) {//module0
  262. if ( !GetHit(slice,module0,zdc) ) {
  263. //cout <<"INSIDE NO GETHIT 1-> " <<slice <<" " <<module0 <<" " <<zdc <<endl;
  264. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, -59, -59, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  265. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  266. AddHit(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  267. //cout <<"INSIDE NO GETHIT 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  268. }
  269. else {
  270. //cout <<"INSIDE GETHIT 1-> " <<slice <<" " <<module0 <<" " <<zdc <<endl;
  271. GetHit(slice,module0,zdc)->AddVSC(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  272. //cout <<"INSIDE GETHIT 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  273. }
  274. }//if(copyNoVTYVECCom==copyNoVTYVEC
  275. /*
  276. if(copyNoVTYVECCom==slice_h && copyNoVMODCom==module_h && copyNoVZDCCom==zdc_h) {//module0 with hole
  277. if ( !GetHitH(slice_h,module_h) ) {
  278. cout <<"INSIDE NO GetHitH 1-> " <<slice_h <<" " <<module_h <<" " <<zdc_h <<endl;
  279. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, -59, -59, slice_h, module_h, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  280. new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice_h, module_h, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  281. cout <<"INSIDE NO GetHitH 2-> " <<slice_h <<" "<<module_h <<" " <<zdc_h <<endl;
  282. }
  283. else {
  284. cout <<"INSIDE GetHitH 1-> " <<slice_h <<" "<<module_h <<" " <<zdc_h <<endl;
  285. GetHitH(slice_h,module_h)->AddVSCH(fTrackID, ivol, slice_h, module_h, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  286. cout <<"INSIDE GetHitH 2-> " <<slice_h <<" "<<module_h <<" " <<zdc_h <<endl;
  287. }
  288. }//if(copyNoVTYVECCom==copyNoVTYVEC
  289. */
  290. /*
  291. std::cout << "INSIDE MpdZdc::ProcessHits: TrackID:" <<part->GetPdgCode() <<" " << fTrackID << " " <<fELoss <<" " << gMC->CurrentVolPath() << " " << tPos.Z() <<" " << ivol << "=="<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<" " <<gMC->CurrentVolOffName(2) << " " <<gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<std::endl;
  292. cout <<"CHECK INSIDE ivol,copyNo,iCell " <<part->GetPdgCode() <<" " << fTrackID <<" " <<ivol <<" " <<module0 <<" " <<slice <<" " <<tPos.Z() <<" " <<tMom.Pz() <<" " <<fELoss <<endl;
  293. AddHit(fTrackID, ivol, slice, module0, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  294. */
  295. /*
  296. std::cout << "INSIDE MpdZdc::ProcessHits: TrackID:" <<part->GetPdgCode() <<" " << fTrackID << " " <<fELoss <<" " << gMC->CurrentVolPath() << " " << tPos.Z() <<" " << ivol << "=="<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<" " <<gMC->CurrentVolOffName(2) << " " <<gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<std::endl;
  297. cout <<"CHECK INSIDE ivol,copyNo,iCell " <<part->GetPdgCode() <<" " << fTrackID <<" " <<ivol <<" " <<module0 <<" " <<slice <<" " <<tPos.Z() <<" " <<tMom.Pz() <<" " <<fELoss <<endl;
  298. AddHit(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  299. */
  300. }//if(fELoss>0)
  301. }//if ( gMC->IsTrackStop() || gMC->IsTrackDisappeared() )
  302. }//if ( gMC->IsTrackInside())
  303. if ( gMC->IsTrackExiting()) {
  304. //if (gMC->Edep() <= 0. ) return kFALSE;
  305. //if (gMC->Edep() > 0. ) fELoss += gMC->Edep();
  306. ////if (gMC->Edep() <= 1.0e-20 ) return kFALSE;
  307. ////if (gMC->Edep() > 1.0e-20 ) fELoss += gMC->Edep();
  308. //part = gMC->GetStack()->GetCurrentTrack();
  309. // Create MpdZdcPoint
  310. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  311. time += gMC->TrackTime() * 1.0e09;
  312. length += gMC->TrackLength();
  313. //fELoss +=gMC->Edep();
  314. //Birk corrections
  315. QCF = 1.+(BirkConst/gMC->TrackStep())*gMC->Edep();
  316. fELoss +=(gMC->Edep())/QCF;
  317. gMC->TrackPosition(tPos);
  318. gMC->TrackMomentum(tMom);
  319. //gMC->CurrentVolOffID(2, module0);
  320. //gMC->CurrentVolOffID(1, slice);
  321. if(fELoss>0) {
  322. //std::cout << "EXIT MpdZdc::ProcessHits: TrackID:" <<part->GetPdgCode() <<" " << fTrackID << " " <<fELoss <<" " << gMC->CurrentVolPath() << " " << tPos.Z() <<" " <<tMom.Pz() <<" " << ivol <<" " <<gMC->CurrentVolOffName(2) << " " <<gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<std::endl;
  323. //cout <<"CHECK EXIT ivol,copyNo,iCell " <<part->GetPdgCode() <<" " << fTrackID <<" " <<ivol <<" " <<copyNoVZDCCom <<" " <<copyNoVMODCom <<" " <<copyNoVTYVECCom <<" " <<slice <<" " <<module0 <<" " <<zdc <<" " <<tPos.Z() <<" " <<tMom.Pz() <<" " <<fELoss <<endl;
  324. /*
  325. if(copyNoVTYVECCom==slice && copyNoVMODCom==module0 && copyNoVZDCCom==zdc) {
  326. cout <<"EXIT slice module0 " <<slice <<" " <<module0 <<" " <<zdc <<endl;
  327. if ( !GetHit(slice,module0,zdc) ) {
  328. cout <<"EXIT NO GETHIT 1-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  329. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, -59, -59, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  330. new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  331. cout <<"EXIT NO GETHIT 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  332. }
  333. else {
  334. cout <<"EXIT GetHit 1-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  335. GetHit(slice,module0,zdc)->AddVSC(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  336. cout <<"EXIT GetHit 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  337. }
  338. }//if(copyNoVTYVECCom==copyNoVTYVEC
  339. */
  340. if(copyNoVTYVECCom==slice && copyNoVMODCom==module0 && copyNoVZDCCom==zdc) {
  341. //cout <<"EXIT slice module0 " <<slice <<" " <<module0 <<" " <<zdc <<endl;
  342. if ( !GetHit(slice,module0,zdc) ) {
  343. //cout <<"EXIT NO GETHIT 1-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  344. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, -59, -59, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  345. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  346. AddHit(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  347. //cout <<"EXIT NO GETHIT 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  348. }
  349. else {
  350. //cout <<"EXIT GetHit 1-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  351. GetHit(slice,module0,zdc)->AddVSC(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  352. //cout <<"EXIT GetHit 2-> " <<slice <<" "<<module0 <<" " <<zdc <<endl;
  353. }
  354. }//if(copyNoVTYVECCom==copyNoVTYVEC
  355. /*
  356. if(copyNoVTYVECCom==slice_h && copyNoVMODCom==module_h && copyNoVZDCCom==zdc_h) {
  357. cout <<"EXIT slice_h module_h " <<slice <<" " <<module0 <<" " <<slice_h <<" " <<module_h <<" " <<zdc <<" " <<zdc_h <<endl;
  358. if ( !GetHitH(slice_h,module_h) ) {
  359. cout <<"EXIT NO GetHitH cxx 1 -> " <<slice_h <<" " <<module_h <<" " <<zdc_h <<endl;
  360. //new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, -59, -59, slice_h, module_h, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  361. new((*fZdcCollection)[fHitNb++]) MpdZdcPoint(fTrackID, ivol, slice_h, module_h, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  362. cout <<"EXIT NO GetHitH cxx 2 -> " <<slice_h <<" " <<module_h <<" " <<zdc_h <<endl;
  363. }
  364. else {
  365. cout <<"EXIT GetHitH cxx 1 -> " <<slice_h <<" " <<module_h <<" " <<zdc_h <<endl;
  366. GetHitH(slice_h,module_h)->AddVSCH(fTrackID, ivol, slice_h, module_h, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  367. cout <<"EXIT GetHitH cxx 2 -> " <<slice_h <<" " <<module_h <<" " <<zdc_h <<endl;
  368. }
  369. }//if(copyNoVTYVECCom==copyNoVTYVEC
  370. */
  371. /*
  372. std::cout << "EXIT MpdZdc::ProcessHits: TrackID:" <<part->GetPdgCode() <<" " << fTrackID << " " <<fELoss <<" " << gMC->CurrentVolPath() << " " << tPos.Z() <<" " << ivol << "=="<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<" " <<gMC->CurrentVolOffName(2) << " " <<gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<std::endl;
  373. cout <<"CHECK EXIT ivol,copyNo,iCell " <<part->GetPdgCode() <<" " << fTrackID <<" " <<ivol <<" " <<module0 <<" " <<slice <<" " <<tPos.X() <<" " <<tPos.Y() <<" " <<tPos.Z() <<" " <<tMom.Pz() <<" " <<fELoss <<endl;
  374. AddHit(fTrackID, ivol, slice, module0, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  375. */
  376. /*
  377. std::cout << "EXIT MpdZdc::ProcessHits: TrackID:" <<part->GetPdgCode() <<" " << fTrackID << " " <<fELoss <<" " << gMC->CurrentVolPath() << " " << tPos.Z() <<" " << ivol << "=="<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<" " <<gMC->CurrentVolOffName(2) << " " <<gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<std::endl;
  378. cout <<"CHECK EXIT ivol,copyNo,iCell " <<part->GetPdgCode() <<" " << fTrackID <<" " <<ivol <<" " <<module0 <<" " <<slice <<" " <<tPos.X() <<" " <<tPos.Y() <<" " <<tPos.Z() <<" " <<tMom.Pz() <<" " <<fELoss <<endl;
  379. AddHit(fTrackID, ivol, slice, module0, zdc, TVector3(tPos.X(), tPos.Y(), tPos.Z()),TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),time, length, fELoss);
  380. */
  381. }//if(fELoss>0)
  382. }//if ( gMC->IsTrackExiting()) {
  383. Int_t points = gMC->GetStack()->GetCurrentTrack()->GetMother(1);
  384. // Int_t nZdcPoints = (points & (1<<30)) >> 30;
  385. // nZdcPoints ++;
  386. // if (nZdcPoints > 1) nZdcPoints = 1;
  387. // points = ( points & ( ~ (1<<30) ) ) | (nZdcPoints << 30);
  388. points = ( points & ( ~ (1<<30) ) ) | (1 << 30);
  389. gMC->GetStack()->GetCurrentTrack()->SetMother(1,points);
  390. ((MpdStack*)gMC->GetStack())->AddPoint(kZDC);
  391. //}
  392. //}
  393. // Int_t copyNo;
  394. // gMC->CurrentVolID(copyNo);
  395. // TString nam = gMC->GetMC()->GetName();
  396. // cout<<"name "<<gMC->GetMC()->GetName()<<endl;
  397. // ResetParameters();
  398. return kTRUE;
  399. }
  400. // ----------------------------------------------------------------------------
  401. // ----- Public method EndOfEvent -----------------------------------------
  402. void MpdZdc::EndOfEvent() {
  403. if (fVerboseLevel) Print();
  404. Reset();
  405. }
  406. // ----- Public method Register -------------------------------------------
  407. void MpdZdc::Register() {
  408. FairRootManager::Instance()->Register("ZdcPoint","Zdc", fZdcCollection, kTRUE);
  409. }
  410. // ----- Public method GetCollection --------------------------------------
  411. TClonesArray* MpdZdc::GetCollection(Int_t iColl) const {
  412. if (iColl == 0) return fZdcCollection;
  413. return NULL;
  414. }
  415. // ----- Public method Print ----------------------------------------------
  416. void MpdZdc::Print() const {
  417. Int_t nHits = fZdcCollection->GetEntriesFast();
  418. cout << "-I- MpdZdc: " << nHits << " points registered in this event."
  419. << endl;
  420. if (fVerboseLevel>1)
  421. for (Int_t i=0; i<nHits; i++) (*fZdcCollection)[i]->Print();
  422. }
  423. // ----- Public method Reset ----------------------------------------------
  424. void MpdZdc::Reset() {
  425. fZdcCollection->Delete();
  426. fPosIndex = 0;
  427. }
  428. // guarda in FairRootManager::CopyClones
  429. // ----- Public method CopyClones -----------------------------------------
  430. void MpdZdc::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset ) {
  431. Int_t nEntries = cl1->GetEntriesFast();
  432. //cout << "-I- MpdZdc: " << nEntries << " entries to add." << endl;
  433. TClonesArray& clref = *cl2;
  434. MpdZdcPoint* oldpoint = NULL;
  435. for (Int_t i=0; i<nEntries; i++) {
  436. oldpoint = (MpdZdcPoint*) cl1->At(i);
  437. Int_t index = oldpoint->GetTrackID() + offset;
  438. oldpoint->SetTrackID(index);
  439. new (clref[fPosIndex]) MpdZdcPoint(*oldpoint);
  440. fPosIndex++;
  441. }
  442. cout << " -I- MpdZdc: " << cl2->GetEntriesFast() << " merged entries."
  443. << endl;
  444. }
  445. // ----- Public method ConstructGeometry ----------------------------------
  446. void MpdZdc::ConstructGeometry() {
  447. TString fileName = GetGeometryFileName();
  448. if(fileName.EndsWith(".root"))
  449. {
  450. LOG(INFO) << "Constructing ZDC geometry from ROOT file " << fileName.Data();
  451. ConstructRootGeometry();
  452. }
  453. /*
  454. else if ( fileName.EndsWith(".geo") )
  455. {
  456. FairLogger::GetLogger()->Info(MESSAGE_ORIGIN, "Constructing ZDC geometry from ASCII file %s", fileName.Data());
  457. ConstructAsciiGeometry();
  458. }
  459. else FairLogger::GetLogger()->Fatal(MESSAGE_ORIGIN, "Geometry format of ZDC file %S not supported.", fileName.Data());
  460. */
  461. FairGeoLoader* geoLoad = FairGeoLoader::Instance();
  462. FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  463. MpdZdcGeo* zdcGeo = new MpdZdcGeo();
  464. zdcGeo->setGeomFile(GetGeometryFileName());
  465. geoFace->addGeoModule(zdcGeo);
  466. Bool_t rc = geoFace->readSet(zdcGeo);
  467. if (rc) zdcGeo->create(geoLoad->getGeoBuilder());
  468. TList* volList = zdcGeo->getListOfVolumes();
  469. // store geo parameter
  470. FairRun *fRun = FairRun::Instance();
  471. FairRuntimeDb *rtdb= FairRun::Instance()->GetRuntimeDb();
  472. MpdZdcGeoPar* par=(MpdZdcGeoPar*)(rtdb->getContainer("MpdZdcGeoPar"));
  473. TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
  474. TObjArray *fPassNodes = par->GetGeoPassiveNodes();
  475. TListIter iter(volList);
  476. FairGeoNode* node = NULL;
  477. FairGeoVolume *aVol=NULL;
  478. while( (node = (FairGeoNode*)iter.Next()) ) {
  479. aVol = dynamic_cast<FairGeoVolume*> ( node );
  480. if ( node->isSensitive() ) {
  481. fSensNodes->AddLast( aVol );
  482. }else{
  483. fPassNodes->AddLast( aVol );
  484. }
  485. }
  486. par->setChanged();
  487. par->setInputVersion(fRun->GetRunId(),1);
  488. ProcessNodes ( volList );
  489. }
  490. //Check if Sensitive-----------------------------------------------------------
  491. Bool_t MpdZdc::CheckIfSensitive(std::string name) {
  492. TString tsname = name;
  493. if (tsname.Contains("zdc01s") || tsname.Contains("ScH")) {
  494. return kTRUE;
  495. }
  496. return kFALSE;
  497. }
  498. //-----------------------------------------------------------------------------
  499. // ----- Private method AddHit --------------------------------------------
  500. MpdZdcPoint* MpdZdc::AddHit(Int_t trackID, Int_t detID, Int_t copyNo, Int_t copyNoMother, Int_t copyNoZdc,
  501. TVector3 pos, TVector3 mom, Double_t time,
  502. Double_t length, Double_t eLoss) {
  503. TClonesArray& clref = *fZdcCollection;
  504. Int_t size = clref.GetEntriesFast();
  505. //cout <<"AddHit " <<trackID <<" " <<detID <<" " <<copyNoMother <<" " <<copyNo <<" " <<pos.Z() <<" " <<eLoss <<endl;
  506. //cout <<"-------------------------- " <<endl;
  507. //return new(clref[size]) MpdZdcPoint(trackID, detID, copyNo, copyNoMother,copyNo_h, copyNoMother_h,pos, mom, time, length, eLoss);
  508. return new(clref[size]) MpdZdcPoint(trackID, detID, copyNo, copyNoMother, copyNoZdc, pos, mom, time, length, eLoss);
  509. }
  510. // ----
  511. ClassImp(MpdZdc)