MpdNDet.cxx 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550
  1. /** MpdNDet.cxx
  2. *@author Mikhail Prokudin
  3. **
  4. ** Defines the active detector NDet with geometry coded here.
  5. ** Layers, holes, fibers,steel tapes implemented
  6. **/
  7. #include "MpdNDet.h"
  8. #include "MpdNDetPointLite.h"
  9. #include "MpdNDetPoint.h"
  10. #include "MpdNDetParam.h"
  11. #include "FairGeoInterface.h"
  12. #include "FairGeoLoader.h"
  13. #include "FairGeoNode.h"
  14. #include "FairGeoRootBuilder.h"
  15. #include "FairRuntimeDb.h"
  16. #include "FairRootManager.h"
  17. #include "FairRun.h"
  18. #include "FairRunAna.h"
  19. #include "FairMCTrack.h"
  20. #include "FairStack.h"
  21. #include "FairVolume.h"
  22. #include "FairGeoMedium.h"
  23. #include "FairGeoMedia.h"
  24. #include "TClonesArray.h"
  25. #include "TGeoMCGeometry.h"
  26. #include "TGeoManager.h"
  27. #include "TParticle.h"
  28. #include "TVirtualMC.h"
  29. #include "TGeoBBox.h"
  30. #include "TGeoPgon.h"
  31. #include "TGeoTube.h"
  32. #include "TGeoMatrix.h"
  33. #include "TList.h"
  34. #include <iostream>
  35. #include <stdlib.h>
  36. using namespace std;
  37. #define kN kNumberOfNdetSensitiveVolumes //sensitive air and scintillator
  38. // ----- Default constructor -------------------------------------------
  39. MpdNDet::MpdNDet() : FairDetector("NDET", kTRUE) {
  40. // MpdNDetPoint::Class() ->IgnoreTObjectStreamer();
  41. // MpdNDetPointLite::Class()->IgnoreTObjectStreamer();
  42. // FairMCTrack::Class() ->IgnoreTObjectStreamer();
  43. fNDetCollection = new TClonesArray("MpdNDetPoint");
  44. fLiteCollection = new TClonesArray("MpdNDetPointLite");
  45. fPosIndex = 0;
  46. fVerboseLevel = 1;
  47. Int_t i;
  48. for(i=kN-1;i>-1;i--)
  49. fVolArr[i]=-1111;
  50. }
  51. // -------------------------------------------------------------------------
  52. // ----- Standard constructor ------------------------------------------
  53. MpdNDet::MpdNDet(const char* name, Bool_t active, const char* fileGeo)
  54. : FairDetector(name, active)
  55. {
  56. /** MpdNDet constructor:
  57. ** reads geometry parameters from the ascii file <fileGeo>,
  58. ** creates the Ndet geometry container MpdNdetInf
  59. ** and initializes basic geometry parameters needed to construct
  60. ** TGeo geometry
  61. **/
  62. Int_t i;
  63. Int_t j;
  64. TString nm;
  65. Info("MpdNDet","Geometry is read from file %s.", fileGeo);
  66. for(i=kN-1;i>-1;i--)
  67. fVolArr[i]=-1111;
  68. fNDetCollection = new TClonesArray("MpdNDetPoint");
  69. fLiteCollection = new TClonesArray("MpdNDetPointLite");
  70. fPosIndex = 0;
  71. fVerboseLevel = 1;
  72. MpdNDetParam* conf=new MpdNDetParam("NDetParam",fileGeo);
  73. fZSize=conf->GetDouble("zsize");
  74. fAirInnerRad=conf->GetDouble("airinnerradius");
  75. fAirThickness=conf->GetDouble("airthickness");
  76. fScinThickness=conf->GetDouble("scinthickness");
  77. fHCut=conf->GetDouble("hcut");
  78. fECut=conf->GetDouble("ecut");
  79. fZDivisions=conf->GetInteger("zdivisions");
  80. fRDivisions=conf->GetInteger("rdivisions");// phi divisions
  81. }
  82. // -------------------------------------------------------------------------
  83. void MpdNDet::Initialize()
  84. {
  85. FairDetector::Initialize();
  86. /*
  87. FairRun* sim = FairRun::Instance();
  88. FairRuntimeDb* rtdb=sim->GetRuntimeDb();
  89. FairGeoEcalPar *par=new FairGeoEcalPar();
  90. fInf->FillGeoPar(par,0);
  91. rtdb->addContainer(par);
  92. */
  93. }
  94. // ----- Destructor ----------------------------------------------------
  95. MpdNDet::~MpdNDet()
  96. {
  97. if (fNDetCollection) {
  98. fNDetCollection->Delete();
  99. delete fNDetCollection;
  100. fNDetCollection=NULL;
  101. }
  102. if (fLiteCollection) {
  103. fLiteCollection->Delete();
  104. delete fLiteCollection;
  105. fLiteCollection=NULL;
  106. }
  107. }
  108. // -------------------------------------------------------------------------
  109. // ----- Private method SetNDetCuts ------------------------------------
  110. void MpdNDet::SetNdetCuts(Int_t medium)
  111. {
  112. /** Set GEANT3 tracking energy cuts for selected medium **/
  113. if (fECut > 0) {
  114. gMC->Gstpar(medium,"CUTGAM",fECut);
  115. gMC->Gstpar(medium,"CUTELE",fECut);
  116. gMC->Gstpar(medium,"BCUTE" ,fECut);
  117. gMC->Gstpar(medium,"BCUTM" ,fECut);
  118. }
  119. if (fHCut > 0) {
  120. gMC->Gstpar(medium,"CUTNEU",fHCut);
  121. gMC->Gstpar(medium,"CUTHAD",fHCut);
  122. gMC->Gstpar(medium,"CUTMUO",fHCut);
  123. gMC->Gstpar(medium,"PPCUTM",fHCut);
  124. }
  125. }
  126. // -------------------------------------------------------------------------
  127. void MpdNDet::FinishPrimary()
  128. {
  129. fFirstNumber=fLiteCollection->GetEntriesFast();
  130. }
  131. //_____________________________________________________________________________
  132. void MpdNDet::ChangeHit(MpdNDetPointLite* oldHit)
  133. {
  134. Double_t edep = gMC->Edep();
  135. Double_t el=oldHit->GetEnergyLoss();
  136. Double_t ttime=gMC->TrackTime()*1.0e9;
  137. oldHit->SetEnergyLoss(el+edep);
  138. if(ttime<oldHit->GetTime())
  139. oldHit->SetTime(ttime);
  140. }
  141. //_____________________________________________________________________________
  142. void MpdNDet::SetSpecialPhysicsCuts()
  143. {
  144. /** Change the special tracking cuts for
  145. ** two NDet media, Scintillator and Lead
  146. **/
  147. FairRun* fRun = FairRun::Instance();
  148. if (strcmp(fRun->GetName(),"TGeant3") == 0) {
  149. Int_t mediumID;
  150. mediumID = gGeoManager->GetMedium("SensAir")->GetId();
  151. SetNdetCuts(mediumID);
  152. mediumID = gGeoManager->GetMedium("NDetScin")->GetId();
  153. SetNdetCuts(mediumID);
  154. }
  155. }
  156. // ----- Public method ProcessHits --------------------------------------
  157. Bool_t MpdNDet::ProcessHits(FairVolume* vol)
  158. {
  159. /** Fill MC point for sensitive NDet volumes **/
  160. fELoss = gMC->Edep();
  161. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  162. fTime = gMC->TrackTime()*1.0e09;
  163. fLength = gMC->TrackLength();
  164. if (vol->getVolumeId()==fStructureId)
  165. if (gMC->IsTrackEntering())
  166. {
  167. FillWallPoint();
  168. TParticle* p=((FairStack*)gMC->GetStack())->GetParticle(fTrackID);
  169. Int_t points=p->GetMother(1);
  170. Int_t ep=(points&(15<<24))>>24;
  171. ep++; if (ep>15) ep=15;
  172. points=(points&(~(15<<24)))|(ep<<24);
  173. p->SetMother(1, points);
  174. ResetParameters();
  175. return kTRUE;
  176. }
  177. else
  178. return kFALSE;
  179. if (fELoss>0)
  180. {
  181. Int_t i;
  182. TParticle* p=gMC->GetStack()->GetCurrentTrack();
  183. Double_t x, y, z;
  184. Double_t px;
  185. Double_t py;
  186. Double_t dx;
  187. Int_t mphi;
  188. Int_t mz;
  189. gMC->TrackPosition(x, y, z);
  190. // cout << "Id: " << p->GetPdgCode() << " (" << x << ", " << y << ", ";
  191. // cout << z << "): ";
  192. // cout << endl;
  193. gMC->CurrentVolOffID(0, mphi); mphi--;
  194. gMC->CurrentVolOffID(1, mz); mz--;
  195. Int_t id=mphi+1+mz*1000;
  196. fVolumeID=id;
  197. FillLitePoint(0);
  198. // type=fInf->GetType(mx, my);
  199. // cx=cell%type;
  200. // cy=cell/type;
  201. // px=mx*fModuleSize-fEcalSize[0]/2.0+cx*fModuleSize/type+1.0;
  202. // py=my*fModuleSize-fEcalSize[1]/2.0+cy*fModuleSize/type+1.0;
  203. // cout << "(" << px << ", " << py << "|" << type << "): ";
  204. /* for(i=0;i<5;i++)
  205. {
  206. Int_t t;
  207. gMC->CurrentVolOffID(i, t);
  208. cout << i << ": " << gMC->CurrentVolOffName(i) << " " << t << "; ";
  209. }
  210. cout << endl;
  211. */
  212. }
  213. TParticle* p=((FairStack*)gMC->GetStack())->GetParticle(fTrackID);
  214. Int_t points=p->GetMother(1);
  215. Int_t ep=(points&(15<<24))>>24;
  216. ep++; if (ep>15) ep=15;
  217. points=(points&(~(15<<24)))|(ep<<24);
  218. p->SetMother(1, points);
  219. ResetParameters();
  220. return kTRUE;
  221. }
  222. /** returns type of volume **/
  223. Int_t MpdNDet::GetVolType(Int_t volnum)
  224. {
  225. Int_t i;
  226. for(i=kN-1;i>-1;i--) {
  227. if (fVolArr[i]==volnum) break;
  228. }
  229. return i;
  230. }
  231. //-----------------------------------------------------------------------------
  232. void MpdNDet::FillWallPoint()
  233. {
  234. /** Fill MC points on the NDet front wall **/
  235. gMC->TrackPosition(fPos);
  236. gMC->TrackMomentum(fMom);
  237. fVolumeID = -1;
  238. Double_t mass = gMC->TrackMass();
  239. // Calculate kinetic energy
  240. Double_t etot = TMath::Sqrt( fMom.Px()*fMom.Px() +
  241. fMom.Py()*fMom.Py() +
  242. fMom.Pz()*fMom.Pz() +
  243. mass * mass );// - mass;
  244. fELoss = etot;
  245. // Create MpdNDetPoint at the entrance of neutron detector
  246. // for particles with pz>0 coming through the front wall
  247. AddHit(fTrackID, fVolumeID, TVector3(fPos.X(), fPos.Y(), fPos.Z()),
  248. TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, fLength,
  249. fELoss);
  250. fTrackID=gMC->GetStack()->GetCurrentTrackNumber();
  251. ((FairStack*)gMC->GetStack())->AddPoint(kNDET);
  252. }
  253. //-----------------------------------------------------------------------------
  254. MpdNDetPointLite* MpdNDet::FindHit(Int_t VolId, Int_t TrackId)
  255. {
  256. for(Int_t i=fFirstNumber;i<fLiteCollection->GetEntriesFast();i++)
  257. {
  258. MpdNDetPointLite* point=(MpdNDetPointLite*)fLiteCollection->At(i);
  259. if (point->GetTrackID()==TrackId&&point->GetDetectorID()==VolId)
  260. return point;
  261. }
  262. return NULL;
  263. }
  264. //-----------------------------------------------------------------------------
  265. Bool_t MpdNDet::FillLitePoint(Int_t volnum)
  266. {
  267. /** Fill MC points inside the NDet for non-zero deposited energy **/
  268. //Search for input track
  269. static Double_t rmin=fAirInnerRad-0.001;
  270. static Double_t rmax=fAirInnerRad+fAirThickness+0.001;
  271. TParticle* part=gMC->GetStack()->GetCurrentTrack();
  272. fTrackID=gMC->GetStack()->GetCurrentTrackNumber();
  273. Double_t vx=part->Vx(); vx*=vx;
  274. Double_t vy=part->Vy(); vy*=vy;
  275. Double_t r=vx; r+=vy; r=TMath::Sqrt(r);
  276. while (part->GetFirstMother()>=0&&\
  277. r<=rmax&&r>=rmin&&\
  278. TMath::Abs(part->Vz())<=fZSize)
  279. {
  280. fTrackID=part->GetFirstMother();
  281. part =((FairStack*)gMC->GetStack())->GetParticle(fTrackID);
  282. vx=part->Vx(); vx*=vx;
  283. vy=part->Vy(); vy*=vy;
  284. r=vx; r+=vy; r=TMath::Sqrt(r);
  285. }
  286. #ifdef _DECAL
  287. if (fTrackID<0) cout<<"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!fTrackID="<<fTrackID<<endl;
  288. #endif
  289. MpdNDetPointLite* oldHit;
  290. MpdNDetPointLite* newHit;
  291. if ((oldHit=FindHit(fVolumeID,fTrackID))!=NULL)
  292. ChangeHit(oldHit);
  293. else {
  294. // Create MpdNDetPoint for scintillator volumes
  295. newHit = AddLiteHit(fTrackID, fVolumeID, fTime, fELoss);
  296. ((FairStack*)gMC->GetStack())->AddPoint(kNDET);
  297. }
  298. return kTRUE;
  299. }
  300. // ----- Public method EndOfEvent --------------------------------------
  301. void MpdNDet::EndOfEvent() {
  302. if (fVerboseLevel) Print();
  303. fNDetCollection->Clear();
  304. fLiteCollection->Clear();
  305. fPosIndex = 0;
  306. fFirstNumber=0;
  307. }
  308. // -------------------------------------------------------------------------
  309. // ----- Public method GetCollection -----------------------------------
  310. TClonesArray* MpdNDet::GetCollection(Int_t iColl) const
  311. {
  312. if (iColl == 0) return fNDetCollection;
  313. return NULL;
  314. // if (iColl == 1) return fLiteCollection;
  315. // return NULL;
  316. }
  317. // -------------------------------------------------------------------------
  318. // ----- Public method Reset -------------------------------------------
  319. void MpdNDet::Reset()
  320. {
  321. fNDetCollection->Clear();
  322. fLiteCollection->Clear();
  323. ResetParameters();
  324. fFirstNumber=0;
  325. }
  326. // -------------------------------------------------------------------------
  327. // ----- Public method Print -------------------------------------------
  328. void MpdNDet::Print() const
  329. {
  330. Int_t nHits = fNDetCollection->GetEntriesFast();
  331. Int_t nLiteHits;
  332. Int_t i;
  333. cout << "-I- MpdNDet: " << nHits << " points registered in this event.";
  334. cout << endl;
  335. nLiteHits = fLiteCollection->GetEntriesFast();
  336. cout << "-I- MpdNDet: " << nLiteHits << " lite points registered in this event.";
  337. cout << endl;
  338. if (fVerboseLevel>1) {
  339. for (i=0; i<nHits; i++) (*fNDetCollection)[i]->Print();
  340. for (i=0; i<nLiteHits; i++) (*fLiteCollection)[i]->Print();
  341. }
  342. }
  343. // -------------------------------------------------------------------------
  344. // ----- Public method CopyClones --------------------------------------
  345. void MpdNDet::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)
  346. {
  347. Int_t nEntries = cl1->GetEntriesFast();
  348. Int_t i;
  349. Int_t index;
  350. cout << "-I- MpdNDet: " << nEntries << " entries to add." << endl;
  351. TClonesArray& clref = *cl2;
  352. if (cl1->GetClass()==MpdNDetPoint::Class()) {
  353. MpdNDetPoint* oldpoint = NULL;
  354. for (i=0; i<nEntries; i++) {
  355. oldpoint = (MpdNDetPoint*) cl1->At(i);
  356. index = oldpoint->GetTrackID()+offset;
  357. oldpoint->SetTrackID(index);
  358. new (clref[fPosIndex]) MpdNDetPoint(*oldpoint);
  359. fPosIndex++;
  360. }
  361. cout << "-I- MpdNDet: " << cl2->GetEntriesFast() << " merged entries."
  362. << endl;
  363. }
  364. else if (cl1->GetClass()==MpdNDetPointLite::Class()) {
  365. MpdNDetPointLite* oldpoint = NULL;
  366. for (i=0; i<nEntries; i++) {
  367. oldpoint = (MpdNDetPointLite*) cl1->At(i);
  368. index = oldpoint->GetTrackID()+offset;
  369. oldpoint->SetTrackID(index);
  370. new (clref[fPosIndex]) MpdNDetPointLite(*oldpoint);
  371. fPosIndex++;
  372. }
  373. cout << "-I- MpdNDet: " << cl2->GetEntriesFast() << " merged entries."
  374. << endl;
  375. }
  376. }
  377. // -------------------------------------------------------------------------
  378. // ----- Public method Register ----------------------------------------
  379. void MpdNDet::Register()
  380. {
  381. FairRootManager::Instance()->Register("NDetPoint","NDet",fNDetCollection,kTRUE);
  382. FairRootManager::Instance()->Register("NDetPointLite","NDetLite",fLiteCollection,kTRUE);
  383. }
  384. // -------------------------------------------------------------------------
  385. // ----- Public method ConstructGeometry -------------------------------
  386. void MpdNDet::ConstructGeometry()
  387. {
  388. TGeoVolume *volume;
  389. TGeoVolume *vol1;
  390. TGeoVolume *vol2;
  391. Float_t *buf = 0;
  392. Int_t i;
  393. Double_t par[10];
  394. // create SensVacuum which is defined in the media file
  395. /** Initialize all media **/
  396. InitMedia();
  397. par[0]=fAirInnerRad;
  398. par[1]=fAirInnerRad+fAirThickness;
  399. par[2]=fZSize+0.001;
  400. volume=gGeoManager->Volume("NDet", "TUBE", gGeoManager->GetMedium("SensAir")->GetId(), par, 3);
  401. gGeoManager->Node("NDet", 1, "cave", 0.0, 0.0, 0.0, 0, kTRUE, buf, 0);
  402. // An ugly way!!!
  403. // Need to make a two volumes for each Ndet
  404. AddSensitiveVolume(volume);
  405. fStructureId=volume->GetNumber();
  406. par[0]=fAirInnerRad+0.001;
  407. par[1]=fAirInnerRad+fScinThickness+0.001;
  408. par[2]=fZSize;
  409. volume=gGeoManager->Volume("NDetScin", "TUBE", gGeoManager->GetMedium("NDetScin")->GetId(), par, 3);
  410. gGeoManager->Node("NDetScin", 1, "NDet", 0.0, 0.0, 0.0, 0, kTRUE, buf, 0);
  411. vol1=gGeoManager->Division("NDetColumn","NDetScin", 3, -1,-fZSize,fZSize*2/fZDivisions,0,"S");
  412. vol2=gGeoManager->Division("NDetRaw","NDetColumn" , 2,fRDivisions,0.0,-1,0,"N");
  413. AddSensitiveVolume(vol2);
  414. fScinId=vol2->GetNumber();
  415. //gGeoManager->Node("EcalStructure", 1, "Ecal", 0.0, 0.0, 0.0, 0, kTRUE, buf, 0);
  416. }
  417. // -------------------------------------------------------------------------
  418. // ----- Public method BeginEvent -----------------------------------------
  419. void MpdNDet::BeginEvent()
  420. {
  421. }
  422. // -------------------------------------------------------------------------
  423. // -------------------------------------------------------------------------
  424. // ----- Private method AddHit -----------------------------------------
  425. MpdNDetPoint* MpdNDet::AddHit(Int_t trackID, Int_t detID, TVector3 pos,
  426. TVector3 mom, Double_t time, Double_t length,
  427. Double_t eLoss)
  428. {
  429. TClonesArray& clref = *fNDetCollection;
  430. Int_t size = clref.GetEntriesFast();
  431. return new(clref[size]) MpdNDetPoint(trackID, detID, pos, mom,
  432. time, length, eLoss);
  433. }
  434. // -------------------------------------------------------------------------
  435. // ----- Private method AddHit -----------------------------------------
  436. MpdNDetPointLite* MpdNDet::AddLiteHit(Int_t trackID, Int_t detID, Double32_t time, Double32_t eLoss)
  437. {
  438. TClonesArray& clref = *fLiteCollection;
  439. Int_t size = clref.GetEntriesFast();
  440. return new(clref[size]) MpdNDetPointLite(trackID, detID, time, eLoss);
  441. }
  442. // -------------------------------------------------------------------------
  443. // ----- Private method InitMedium ---------------------------------------
  444. Int_t MpdNDet::InitMedium(const char* name)
  445. {
  446. static FairGeoLoader *geoLoad=FairGeoLoader::Instance();
  447. static FairGeoInterface *geoFace=geoLoad->getGeoInterface();
  448. static FairGeoMedia *media=geoFace->getMedia();
  449. static FairGeoBuilder *geoBuild=geoLoad->getGeoBuilder();
  450. FairGeoMedium *FairMedium=media->getMedium(name);
  451. if (!FairMedium)
  452. {
  453. Fatal("InitMedium","Material %s not defined in media file.", name);
  454. return -1111;
  455. }
  456. return geoBuild->createMedium(FairMedium);
  457. }
  458. // -------------------------------------------------------------------------
  459. // ----- Private method InitMedia ----------------------------------------
  460. void MpdNDet::InitMedia()
  461. {
  462. Info("InitMedia", "Initializing media.");
  463. InitMedium("SensAir");
  464. InitMedium("NDetScin");
  465. }
  466. // -------------------------------------------------------------------------
  467. // ----- Public method GetCellCoordInf ----------------------------------------
  468. // TODO: !!!
  469. Bool_t MpdNDet::GetCellCoordInf(Int_t fVolID, Float_t &x, Float_t &y, Int_t& tenergy)
  470. {
  471. return kFALSE;
  472. }
  473. // ------------------------------------------------------------------------------
  474. Bool_t MpdNDet::GetCellCoord(Int_t fVolID, Float_t &x, Float_t &y, Int_t& tenergy)
  475. {
  476. return GetCellCoordInf(fVolID, x, y, tenergy);
  477. }
  478. ClassImp(MpdNDet)