MpdTgem.cxx 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357
  1. //------------------------------------------------------------------------------------------------------------------------
  2. // -------------------------------------------------------------------------
  3. // ----- MpdTgem source file -----
  4. // -------------------------------------------------------------------------
  5. #include <iostream>
  6. #include "TClonesArray.h"
  7. #include "TLorentzVector.h"
  8. #include "TParticle.h"
  9. #include "TVirtualMC.h"
  10. #include "TGeoManager.h"
  11. #include "TGeoMatrix.h"
  12. #include "FairGeoInterface.h"
  13. #include "FairGeoLoader.h"
  14. #include "FairGeoNode.h"
  15. #include "FairGeoRootBuilder.h"
  16. #include "MpdTgemGeo.h"
  17. #include "FairRootManager.h"
  18. #include "MpdTgem.h"
  19. #include "MpdTgemPoint.h"
  20. #include "FairRuntimeDb.h"
  21. #include "MpdTgemGeoPar.h"
  22. #include "TObjArray.h"
  23. #include "FairRun.h"
  24. #include "FairVolume.h"
  25. #include "TMath.h"
  26. #include "TParticlePDG.h"
  27. class FairVolume;
  28. //------------------------------------------------------------------------------------------------------------------------
  29. MpdTgem::MpdTgem()
  30. : FairDetector("Tgem", kTRUE)
  31. {
  32. fTgemCollection = new TClonesArray("MpdTgemPoint");
  33. fPosIndex = 0;
  34. fVerboseLevel = 1;
  35. ResetParameters();
  36. }
  37. //------------------------------------------------------------------------------------------------------------------------
  38. MpdTgem::MpdTgem(const char* name, Bool_t active)
  39. : FairDetector(name, active)
  40. {
  41. fTgemCollection = new TClonesArray("MpdTgemPoint");
  42. fPosIndex = 0;
  43. fVerboseLevel = 1;
  44. ResetParameters();
  45. }
  46. //------------------------------------------------------------------------------------------------------------------------
  47. MpdTgem::~MpdTgem()
  48. {
  49. if(fTgemCollection){fTgemCollection->Delete(); delete fTgemCollection; }
  50. }
  51. //------------------------------------------------------------------------------------------------------------------------
  52. int MpdTgem::DistAndPoints(TVector3 p3, TVector3 p4, TVector3& pa, TVector3& pb) {
  53. pa=(p3+p4)*0.5;
  54. pb=pa;
  55. //pa=p3; //del
  56. //pb=pa; //del
  57. return 0;
  58. }
  59. //------------------------------------------------------------------------------------------------------------------------
  60. TVector3 MpdTgem::GlobalToLocal(TVector3& global) {
  61. Double_t globPos[3];
  62. Double_t localPos[3];
  63. global.GetXYZ(globPos);
  64. gMC->Gmtod(globPos, localPos, 1);
  65. return TVector3(localPos);
  66. }
  67. //------------------------------------------------------------------------------------------------------------------------
  68. TVector3 MpdTgem::LocalToGlobal(TVector3& local) {
  69. Double_t globPos[3];
  70. Double_t localPos[3];
  71. local.GetXYZ(localPos);
  72. gMC->Gdtom(localPos, globPos, 1);
  73. return TVector3(globPos);
  74. }
  75. //----------------------------------------------------------------------------------------------------------------------
  76. Bool_t MpdTgem::ProcessHits(FairVolume* vol)
  77. {
  78. //AZ
  79. static Int_t first = 1, nLayMod = 30; //!!! geometry dependent
  80. static Double_t z0[4] = {156.5, 186.5, 0, 0};
  81. first=0; //=
  82. if (first) {
  83. first = 0;
  84. TGeoHMatrix matrix;
  85. if (gGeoManager) {
  86. if (gGeoManager->GetCurrentNavigator()->CheckPath("/cave_1/stt01layerradial_1")) { //EL
  87. if (gMC->GetTransformation("/cave_1/stt01layerradial_1",matrix))
  88. z0[0] = TMath::Abs (matrix.GetTranslation()[2]);
  89. }
  90. //matrix.Print();
  91. //cout << gGeoManager->GetPath() << " " << matrix.GetTranslation()[2] << endl;
  92. if (gGeoManager->GetCurrentNavigator()->CheckPath("/cave_1/stt02layerradial_1")) { //EL
  93. if (gMC->GetTransformation("/cave_1/stt02layerradial_1",matrix))
  94. z0[1] = TMath::Abs (matrix.GetTranslation()[2]);
  95. }
  96. if (gGeoManager->GetCurrentNavigator()->CheckPath("/cave_1/stt03layerradial_1")) {
  97. if (gMC->GetTransformation("/cave_1/stt03layerradial_1",matrix))
  98. z0[2] = TMath::Abs (matrix.GetTranslation()[2]);
  99. }
  100. if (gGeoManager->GetCurrentNavigator()->CheckPath("/cave_1/stt04layerradial_1")) {
  101. if (gMC->GetTransformation("/cave_1/stt04layerradial_1",matrix))
  102. z0[3] = TMath::Abs (matrix.GetTranslation()[2]);
  103. }
  104. }
  105. }
  106. //AZ
  107. Int_t gap, cell, module, region, wheel, gasgap;
  108. static Int_t sameVol = 0;
  109. TString Volname = vol->getRealName(); // EL
  110. sscanf(&(Volname[11]),"%d",&gasgap);
  111. //==++cout <<gasgap<<"===============-----------------Volname="<< Volname << endl;
  112. sscanf(&(gMC->CurrentVolPath()[15]),"%d",&wheel);
  113. //=cout<<"************--"<<wheel<<" "<<gMC->CurrentVolPath()<<endl;
  114. // Set parameters at entrance of volume. Reset ELoss.
  115. // Check if this is not due to passing through the anode wire.
  116. // If anode wire, do not reset anything.
  117. if (gMC->IsTrackEntering()) {
  118. sameVol = 0; //== 1 ==
  119. gMC->CurrentVolOffID(2, module);
  120. sscanf(&(gMC->CurrentVolPath()[15]),"%d",&wheel);
  121. //cout<<"************--"<<wheel<<" "<<gMC->CurrentVolPath()<<endl;
  122. // cout<<"-----!!!!!!!!!!!!!! "<<gMC->GetStack()->GetCurrentTrackNumber()<<" "<<fTrackID<<" "<<gMC->CurrentVolID(gap)+1000*module<<" "<<fVolumeID<<endl;
  123. if ((gMC->GetStack()->GetCurrentTrackNumber() != fTrackID ||
  124. wheel != fwheel)) {
  125. //if(fTrackIDTmp != -1) {ResetParameters();sameVol = 0; cout<<" !!!!!!!!==========ERROR==========!!!!!!!!! "<<fPdgId<<endl; return kTRUE;}
  126. // gMC->CurrentVolID(gap)+1000*module != fVolumeID) {
  127. ResetParameters();
  128. sameVol = 0;
  129. fELoss = 0.;
  130. fTime = gMC->TrackTime() * 1.0e09;
  131. fLength = gMC->TrackLength();
  132. fIsPrimary = 0;
  133. fCharge = -1;
  134. fPdgId = 0;
  135. if(gasgap == 1){
  136. TLorentzVector PosIn;
  137. gMC->TrackPosition(PosIn);
  138. fPosIn.SetXYZ(PosIn.X(), PosIn.Y(), PosIn.Z());
  139. fPosInTmp.SetXYZ(PosIn.X(), PosIn.Y(), PosIn.Z());
  140. //==++cout <<"===============-----------------"<< Volname <<" PosIn.Z()="<<PosIn.Z()<< endl;
  141. gMC->TrackMomentum(fMom);
  142. TParticle* part = 0;
  143. part = gMC->GetStack()->GetCurrentTrack();
  144. if (part) {
  145. fIsPrimary = (Int_t)part->IsPrimary();
  146. fCharge = (Int_t)part->GetPDG()->Charge();
  147. fPdgId = (Int_t)part->GetPdgCode();
  148. }
  149. fVolumeID = gMC->CurrentVolID(gap) + 1000 * module;
  150. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  151. fTrackIDTmp = fTrackID;
  152. fwheel = wheel; // wheel number
  153. fwheelTmp = fwheel;
  154. }
  155. }
  156. }
  157. // Sum energy loss for all steps in the active volume
  158. fELoss += gMC->Edep();
  159. // Create MpdTgemPoint at EXIT of active volume;
  160. //if (gMC->IsTrackExiting() && fELoss > 0) {
  161. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  162. if ((gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) && fELoss > 0 && gasgap == 2 &&
  163. fTrackID == fTrackIDTmp && fwheelTmp == fwheel) {
  164. //if(fTrackID != fTrackIDTmp || fwheelTmp != fwheel) { return 0;}
  165. //==++ cout<<"************* gMC->GetStack()->GetCurrentTrackNumber()="<<fTrackID<<" "<<"fPdgId="<<fPdgId<<endl;
  166. //= Volname = vol->getName();
  167. //= cout <<"===============-----------------"<< Volname << endl;
  168. //AZ region = Volname[5] - '0'; //?????????????????????????
  169. //gMC->CurrentVolID(gap);
  170. //gMC->CurrentVolOffID(1, cell);
  171. //gMC->CurrentVolOffID(2, module);
  172. TLorentzVector PosOut;
  173. gMC->TrackPosition(PosOut);
  174. fPosOut.SetXYZ(PosOut.X(), PosOut.Y(), PosOut.Z());
  175. //AZ Get layer number
  176. //== sscanf(&(Volname[7]),"%d",&wheel); // wheel number
  177. //==++cout <<"===============-----------------"<< Volname <<" PosOut.Z()="<<PosOut.Z()<< endl;
  178. //= cout <<"===============-----------------"<< wheel << endl;
  179. Double_t dz = TMath::Abs(PosOut.Z()) - z0[wheel-1];
  180. //Int_t lay = TMath::Nint(dz/1.) + 1;
  181. //Int_t lay = TMath::Nint(dz/1.) + 2; // geometry with triplets
  182. Int_t lay = TMath::Nint(dz/0.85) + 1; // geometry with triplets
  183. if (wheel == 2 || wheel == 4) lay += nLayMod;
  184. //= region = lay * 1000 + cell; // lay*1000+tubeID
  185. cell=1; //=
  186. //== cout<<"__________________________wheel="<<wheel<<endl;
  187. region = wheel * 10000 + cell; //= fwheel*10000+segmentID
  188. /*
  189. Double_t posG[3], posL[3];
  190. fPosOut.GetXYZ(posG);
  191. gMC->Gmtod(posG, posL, 1);
  192. cout << lay << " " << cell << " " << fTrackID << " " << posL[0] << " "
  193. << posL[1] << " " << posL[2] << " " << gap << " " << module << endl; */
  194. //AZ
  195. // Straw line defined in local coordinates
  196. TVector3 p1(0,0, -10); // -10,10 - arbitrary number in Z axis...
  197. TVector3 p2(0,0, 10);
  198. // Conversion to global coordinates
  199. p1 = LocalToGlobal(p1);
  200. p2 = LocalToGlobal(p2);
  201. Double_t phi = TMath::ATan2 (p2.Y()-p1.Y(),p2.X()-p1.X()); //AZ
  202. // "will-be-filled-out-soon" Points of closest approach
  203. TVector3 trackPosition(0,0,0); // trackPosition => point on track, fPos => point on straw
  204. // calculate points of closest approach between track and straw
  205. int result = DistAndPoints(fPosInTmp, fPosOut, fPos, trackPosition);
  206. //TVector3 dist = trackPosition - fPos;
  207. //std::cout << "Dist: " << dist.Mag() << std::endl;
  208. //TVector3 trackPosition(0,0,0);
  209. //if (result == 0) { // track and straw should not be parallel!
  210. // convert point coordinates to back to global system
  211. // TVector3 dist = pa - pb;
  212. // std::cout << "distance: " << dist.Mag() << std::endl;
  213. //}
  214. //AZ fVolumeID = ((region-1)<<24);
  215. // Remove point if it is from the same track and in the same tube (to account for the cases
  216. // when the track goes through the anode wire - probably will not work if secondaries are
  217. // produced inside the anode wire (they are transported before the parent track?).
  218. if (sameVol) fTgemCollection->RemoveAt(fTgemCollection->GetEntriesFast()-1);
  219. //== cout<<"*******************************Charge="<<fCharge<<endl;
  220. MpdTgemPoint *p = AddHit(fTrackID, region, fPos, fPos.Perp(),
  221. TVector3(fMom.Px(), fMom.Py(), fMom.Pz()),
  222. fTime, (fLength+gMC->TrackLength())/2, fELoss,
  223. fIsPrimary, fCharge, fPdgId, trackPosition);
  224. //= p->SetPhi(phi); //AZ
  225. //==((FairStack*)gMC->GetStack())->AddPoint(kTgem); //==
  226. //ResetParameters();
  227. fTrackIDTmp = -1;
  228. fwheelTmp = -1;
  229. }
  230. return kTRUE;
  231. }
  232. //------------------------------------------------------------------------------------------------------------------------
  233. void MpdTgem::EndOfEvent()
  234. {
  235. if(fVerboseLevel) Print();
  236. fTgemCollection->Clear();
  237. fPosIndex = 0;
  238. }
  239. //------------------------------------------------------------------------------------------------------------------------
  240. void MpdTgem::Register(){ FairRootManager::Instance()->Register("TgemPoint", "Tgem", fTgemCollection, kTRUE); }
  241. //------------------------------------------------------------------------------------------------------------------------
  242. TClonesArray* MpdTgem::GetCollection(Int_t iColl) const
  243. {
  244. if(iColl == 0) return fTgemCollection;
  245. return NULL;
  246. }
  247. //------------------------------------------------------------------------------------------------------------------------
  248. void MpdTgem::Print() const
  249. {
  250. Int_t nHits = fTgemCollection->GetEntriesFast();
  251. cout << "-I- MpdTgem: " << nHits << " points registered in this event." << endl;
  252. if(fVerboseLevel > 1)
  253. for(Int_t i=0; i<nHits; i++) (*fTgemCollection)[i]->Print();
  254. }
  255. //------------------------------------------------------------------------------------------------------------------------
  256. void MpdTgem::Reset(){ fTgemCollection->Clear(); ResetParameters(); }
  257. //------------------------------------------------------------------------------------------------------------------------
  258. void MpdTgem::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)
  259. {
  260. Int_t nEntries = cl1->GetEntriesFast();
  261. cout << "-I- MpdTgem: " << nEntries << " entries to add." << endl;
  262. TClonesArray& clref = *cl2;
  263. MpdTgemPoint* oldpoint = NULL;
  264. for(Int_t i=0; i<nEntries; i++)
  265. {
  266. oldpoint = (MpdTgemPoint*) cl1->At(i);
  267. Int_t index = oldpoint->GetTrackID() + offset;
  268. oldpoint->SetTrackID(index);
  269. new (clref[fPosIndex]) MpdTgemPoint(*oldpoint);
  270. fPosIndex++;
  271. }
  272. cout << "-I- MpdTgem: " << cl2->GetEntriesFast() << " merged entries." << endl;
  273. }
  274. //------------------------------------------------------------------------------------------------------------------------
  275. void MpdTgem::ConstructGeometry()
  276. {
  277. Int_t count=0;
  278. Int_t count_tot=0;
  279. FairGeoLoader* geoLoad = FairGeoLoader::Instance();
  280. FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  281. MpdTgemGeo* tgemGeo = new MpdTgemGeo();
  282. tgemGeo->setGeomFile(GetGeometryFileName());
  283. geoFace->addGeoModule(tgemGeo);
  284. Bool_t rc = geoFace->readSet(tgemGeo);
  285. if(rc) tgemGeo->create(geoLoad->getGeoBuilder());
  286. TList* volList = tgemGeo->getListOfVolumes();
  287. volList->Print();
  288. // store geo parameter
  289. FairRun *fRun = FairRun::Instance();
  290. FairRuntimeDb *rtdb = FairRun::Instance()->GetRuntimeDb();
  291. MpdTgemGeoPar* par =(MpdTgemGeoPar*)(rtdb->getContainer("MpdTgemGeoPar"));
  292. TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
  293. TObjArray *fPassNodes = par->GetGeoPassiveNodes();
  294. FairGeoNode *node = NULL;
  295. FairGeoVolume *aVol = NULL;
  296. TListIter iter(volList);
  297. while((node = (FairGeoNode*)iter.Next()))
  298. {
  299. aVol = dynamic_cast<FairGeoVolume*> (node);
  300. if(node->isSensitive()){ fSensNodes->AddLast(aVol); count++; }
  301. else fPassNodes->AddLast(aVol);
  302. count_tot++;
  303. }
  304. par->setChanged();
  305. par->setInputVersion(fRun->GetRunId(), 1);
  306. ProcessNodes(volList);
  307. }
  308. //------------------------------------------------------------------------------------------------------------------------
  309. MpdTgemPoint* MpdTgem::AddHit(Int_t trackID, Int_t detID, TVector3 pos, Double_t radius,
  310. TVector3 mom, Double_t time, Double_t length,
  311. Double_t eLoss, Int_t isPrimary, Double_t charge, Int_t pdgId, TVector3 trackPos)
  312. {
  313. TClonesArray& clref = *fTgemCollection;
  314. Int_t size = clref.GetEntriesFast();
  315. //std::cout << "ELoss: " << eLoss << "\n";
  316. return new(clref[size]) MpdTgemPoint(trackID, detID, pos, radius, mom, time, length, eLoss, isPrimary, charge, pdgId, trackPos);
  317. }
  318. //------------------------------------------------------------------------------------------------------------------------
  319. ClassImp(MpdTgem)