CbmStt.cxx 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561
  1. // -------------------------------------------------------------------------
  2. // ----- CbmTof source file -----
  3. // ----- Created 28/07/04 by V. Friese -----
  4. // -------------------------------------------------------------------------
  5. #include <iostream>
  6. #include "TClonesArray.h"
  7. #include "TLorentzVector.h"
  8. #include "TParticle.h"
  9. #include "TVirtualMC.h"
  10. //#include "TGeant3.h"
  11. #include "FairGeoInterface.h"
  12. #include "FairGeoLoader.h"
  13. #include "FairGeoNode.h"
  14. #include "FairGeoRootBuilder.h"
  15. #include "CbmGeoStt.h"
  16. #include "FairRootManager.h"
  17. #include "CbmStt.h"
  18. #include "CbmSttPoint.h"
  19. #include "FairRuntimeDb.h"
  20. #include "CbmGeoSttPar.h"
  21. #include "TObjArray.h"
  22. #include "FairRun.h"
  23. #include "FairVolume.h"
  24. #include "FairStack.h"
  25. class FairVolume;
  26. // TODO: read this from geant initialization
  27. #define innerStrawDiameter 1.
  28. //#define redefineLambdaChargedDecay 0
  29. // ----- Default constructor -------------------------------------------
  30. CbmStt::CbmStt()
  31. {
  32. fSttCollection = new TClonesArray("CbmSttPoint");
  33. fPosIndex = 0;
  34. fVerboseLevel = 1;
  35. fIsInitialized = kFALSE;
  36. fSensNodes = NULL;
  37. fPassNodes = NULL;
  38. fParameters = NULL;
  39. // fRun = NULL;
  40. }
  41. // -------------------------------------------------------------------------
  42. // ----- Standard constructor ------------------------------------------
  43. CbmStt::CbmStt(const char* name, Bool_t active)
  44. : FairDetector(name, active)
  45. {
  46. fSttCollection = new TClonesArray("CbmSttPoint");
  47. fPosIndex = 0;
  48. fVerboseLevel = 1;
  49. fIsInitialized = kFALSE;
  50. fSensNodes = NULL;
  51. fPassNodes = NULL;
  52. fParameters = NULL;
  53. // fRun = NULL;
  54. }
  55. // -------------------------------------------------------------------------
  56. // ----- Destructor ----------------------------------------------------
  57. CbmStt::~CbmStt()
  58. {
  59. if (fSttCollection)
  60. {
  61. fSttCollection->Delete();
  62. delete fSttCollection;
  63. }
  64. }
  65. // -------------------------------------------------------------------------
  66. // ----- Private method GetSquaredDistanceFromWire -----------------------
  67. float CbmStt::GetSquaredDistanceFromWire()
  68. {
  69. TLorentzVector
  70. entryPosition;
  71. float
  72. positionInMother[3],
  73. positionInStraw[3];
  74. gMC->TrackPosition(entryPosition);
  75. positionInMother[0] = entryPosition.X();
  76. positionInMother[1] = entryPosition.Y();
  77. positionInMother[2] = entryPosition.Z();
  78. gMC->Gmtod(positionInMother, positionInStraw, 1);
  79. return positionInStraw[0] * positionInStraw[0] +
  80. positionInStraw[1] * positionInStraw[1];
  81. }
  82. // -------------------------------------------------------------------------
  83. // ----- Public method ProcessHits --------------------------------------
  84. Bool_t CbmStt::InitProcessHits()
  85. {
  86. //cout << "in process hits" << endl;
  87. if (!fIsInitialized)
  88. {
  89. fIsInitialized = kTRUE;
  90. //FairRun *fRun = FairRun::Instance();
  91. FairRuntimeDb
  92. *rtdb= FairRun::Instance()->GetRuntimeDb();
  93. if (!rtdb)
  94. {
  95. cout << "-I- CbmStt: No runtime database found." << endl;
  96. return kFALSE;
  97. }
  98. fParameters = (CbmGeoSttPar*)(rtdb->getContainer("CbmGeoSttPar"));
  99. if (!fParameters)
  100. {
  101. cout << "-I- CbmStt: No geometry container CbmGeoSttPar found in runtime database." << endl;
  102. return kFALSE;
  103. }
  104. fSensNodes = fParameters->GetGeoSensitiveNodes();
  105. if (!fSensNodes)
  106. {
  107. cout << "-I- CbmStt: No sensitive nodes found in geometry container." << endl;
  108. return kFALSE;
  109. }
  110. fPassNodes = fParameters->GetGeoPassiveNodes();
  111. if (!fPassNodes)
  112. {
  113. cout << "-I- CbmStt: No passive nodes found in geometry container." << endl;
  114. return kFALSE;
  115. }
  116. }
  117. return kTRUE;
  118. }
  119. bool CbmStt::Split(string &aDest, string &aSrc, char aDelim)
  120. {
  121. if(aSrc.empty())
  122. return false;
  123. string::size_type
  124. pos = aSrc.find(aDelim);
  125. aDest = aSrc.substr(0, pos);
  126. if(pos != string::npos)
  127. aSrc = aSrc.substr(pos + 1);
  128. else
  129. aSrc = "";
  130. return true;
  131. }
  132. string CbmStt::GetStringPart(string &aSrc, Int_t part, char aDelim)
  133. {
  134. string
  135. retval = "",
  136. sub;
  137. int
  138. counter = 0;
  139. while(Split(sub, aSrc, aDelim))
  140. {
  141. if (counter == part)
  142. {
  143. retval = sub;
  144. break;
  145. }
  146. counter++;
  147. }
  148. return retval;
  149. }
  150. // ----- Public method ProcessHits --------------------------------------
  151. Bool_t CbmStt::ProcessHits(FairVolume* vol)
  152. {
  153. if (!InitProcessHits())
  154. {
  155. return kFALSE;
  156. }
  157. if (gMC->TrackCharge() != 0.)
  158. {
  159. if ( gMC->IsTrackEntering() )
  160. {
  161. if (sqrt(GetSquaredDistanceFromWire()) > (innerStrawDiameter / 4.))
  162. {
  163. // Set parameters at entrance of volume. Reset ELoss.
  164. fELoss = 0.;
  165. fTime = gMC->TrackTime() * 1.0e09;
  166. fLength = gMC->TrackLength();
  167. gMC->TrackPosition(fPos);
  168. gMC->TrackMomentum(fMomIn);
  169. gMC->TrackPosition(fpostotin);// da cancellare
  170. Double_t globalPos[3] = {0., 0., 0.}; // stt1 modified
  171. Double_t localPos[3] = {0., 0., 0.}; // stt1 modified
  172. globalPos[0] = fPos.X();
  173. globalPos[1] = fPos.Y();
  174. globalPos[2] = fPos.Z();
  175. gMC->Gmtod(globalPos, localPos, 1);
  176. fPosInLocal.SetXYZM(localPos[0], localPos[1], localPos[2], 0.0);
  177. }
  178. }
  179. // Sum energy loss for all steps in the active volume
  180. fELoss += gMC->Edep();
  181. // Create CbmSttPoint at exit of active volume -- but not into the wire
  182. if (gMC->IsTrackExiting() && (sqrt(GetSquaredDistanceFromWire()) > (innerStrawDiameter / 4.)))
  183. {
  184. fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
  185. fVolumeID = vol->getMCid();
  186. fMass = gMC->TrackMass(); // mass (GeV)
  187. gMC->TrackPosition(fPosOut);
  188. gMC->TrackMomentum(fMomOut);
  189. gMC->TrackPosition(fpostotout);// da cancellare
  190. Double_t globalPos[3] = {0., 0., 0.}; // stt1 modified
  191. Double_t localPos[3] = {0., 0., 0.}; // stt1 modified
  192. gMC->Gdtom(localPos, globalPos, 1);
  193. fPos.SetXYZM(globalPos[0], globalPos[1], globalPos[2], 0.0);
  194. globalPos[0] = fPosOut.X();
  195. globalPos[1] = fPosOut.Y();
  196. globalPos[2] = fPosOut.Z();
  197. gMC->Gmtod(globalPos, localPos, 1);
  198. fPosOutLocal.SetXYZM(localPos[0], localPos[1], localPos[2], 0.0);
  199. // string basename("stt1tube");
  200. string basename;
  201. TString volumename;
  202. string
  203. hashmark("#"),
  204. volName,
  205. fullName,
  206. number,
  207. specialname,
  208. volPath(gMC->CurrentVolPath()),
  209. volPath2(gMC->CurrentVolPath());
  210. volumename = volPath;
  211. volName = GetStringPart(volPath, 2, '/');
  212. number = GetStringPart(volName, 1, '_');
  213. Int_t start =volPath2.find("stt",0);
  214. specialname = volPath2.substr(start,volPath2.find("_",start)-start);
  215. if(volumename.Contains("stt01")) basename = "stt01tube";
  216. if(volumename.Contains("stt02")) basename = "stt02tube";
  217. if(volumename.Contains("stt03")) basename = "stt03tube";
  218. if(volumename.Contains("stt04")) basename = "stt04tube";
  219. if(volumename.Contains("stt05")) basename = "stt05tube";
  220. if(volumename.Contains("stt06")) basename = "stt06tube";
  221. if(volumename.Contains("stt07")) basename = "stt07tube";
  222. if(volumename.Contains("stt08")) basename = "stt08tube";
  223. if(volumename.Contains("stt09")) basename = "stt09tube";
  224. if(volumename.Contains("stt10")) basename = "stt10tube";
  225. if(volumename.Contains("stt11")) basename = "stt11tube";
  226. if(volumename.Contains("stt12")) basename = "stt12tube";
  227. if(volumename.Contains("stt13")) basename = "stt13tube";
  228. if(volumename.Contains("stt14")) basename = "stt14tube";
  229. if(volumename.Contains("stt15")) basename = "stt15tube";
  230. fullName = basename + hashmark + number;
  231. // cout << "volname: " << volName << " " << volumename << endl;
  232. // cout << "number: " << number << endl;
  233. // cout << gMC->CurrentVolPath() << endl;
  234. // cout << "fullname: " << fullName << endl;
  235. FairGeoNode *vol_obj = dynamic_cast<FairGeoNode*> (fPassNodes->FindObject(fullName.c_str()));
  236. if(number=="0") {
  237. vol_obj = dynamic_cast<FairGeoNode*> (fPassNodes->FindObject(specialname.c_str()));
  238. //cout<<">>>>"<<endl;
  239. //cout<<"special "<<specialname.c_str()<<endl;
  240. }
  241. if (!vol_obj)
  242. {
  243. cout << "-I- CbmStt: No volume " << fullName.c_str() << " found in geometry container." << endl;
  244. return kFALSE;
  245. }
  246. FairGeoRotation rotation = vol_obj->getLabTransform()->getRotMatrix();
  247. FairGeoVector
  248. originalVector(0., 0., 1.),
  249. rotatedVector = rotation * originalVector;
  250. // cout << "positionc: " << fPos.X() << " " << fPos.Y() << " " << fPos.Z() << endl; // da cancellare
  251. // cout << "position: " << fpostot.X() << " " << fpostot.Y() << " " << fpostot.Z() << endl;
  252. // if(sqrt(fPosInLocal.X()*fPosInLocal.X() + fPosInLocal.Y()*fPosInLocal.Y()) < 0.45) {
  253. // cout << "position in : " << sqrt(fPosInLocal.X()*fPosInLocal.X() + fPosInLocal.Y()*fPosInLocal.Y()) << endl;
  254. // }
  255. // if(sqrt(fPosOutLocal.X()*fPosOutLocal.X() + fPosOutLocal.Y()*fPosOutLocal.Y()) < 0.45) {
  256. // cout << "position out: " << sqrt(fPosOutLocal.X()*fPosOutLocal.X() + fPosOutLocal.Y()*fPosOutLocal.Y()) << endl;
  257. // }
  258. fpostot.SetXYZM((fpostotin.X() + fpostotout.X())/2., (fpostotin.Y() + fpostotout.Y())/2., (fpostotin.Z() + fpostotout.Z())/2.,0.0); // da cancellare
  259. // cout << "in : " << fpostotin.X() << " " << fpostotin.Y() << " " << fpostotin.Z() << endl;
  260. // cout << "out: " << fpostotout.X() << " " << fpostotout.Y() << " " << fpostotout.Z() << endl;
  261. // cout << "tot: " << fpostot.X() << " " << fpostot.Y() << " " << fpostot.Z() << endl;
  262. // cout << (fpostotin.X() + fpostotout.X())/2. << " " << (fpostotin.Y() + fpostotout.Y())/2. << " " << (fpostotin.Z() + fpostotout.Z())/2. << endl;
  263. AddHit(fTrackID, fVolumeID,
  264. TVector3(fPos.X(), fPos.Y(), fPos.Z()),
  265. TVector3(fPosInLocal.X(), fPosInLocal.Y(), fPosInLocal.Z()),
  266. TVector3(fPosOutLocal.X(), fPosOutLocal.Y(), fPosOutLocal.Z()),
  267. TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()),
  268. TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()),
  269. TVector3(rotatedVector.getX(), rotatedVector.getY(), rotatedVector.getZ()),
  270. fTime, fLength, fELoss, fMass, TVector3(fpostot.X(), fpostot.Y(), fpostot.Z())); // da cancellare fpostot
  271. // Increment number of stt points for TParticle
  272. Int_t
  273. points = gMC->GetStack()->GetCurrentTrack()->GetMother(1);
  274. Int_t
  275. nSttPoints = (points & (15<<20)) >> 20;
  276. nSttPoints ++;
  277. if (nSttPoints > 15)
  278. nSttPoints = 15;
  279. points = ( points & ( ~ (15<<20) ) ) | (nSttPoints << 20);
  280. gMC->GetStack()->GetCurrentTrack()->SetMother(1,points);
  281. ((FairStack*)gMC->GetStack())->AddPoint(kECT);
  282. ResetParameters();
  283. }
  284. }
  285. return kTRUE;
  286. }
  287. // -------------------------------------------------------------------------
  288. // ----- Public method EndOfEvent --------------------------------------
  289. void CbmStt::EndOfEvent()
  290. {
  291. if (fVerboseLevel)
  292. Print();
  293. fSttCollection->Clear();
  294. fPosIndex = 0;
  295. }
  296. // -------------------------------------------------------------------------
  297. // ----- Public method Register ----------------------------------------
  298. void CbmStt::Register()
  299. {
  300. FairRootManager::Instance()->Register("STTPoint", "Stt", fSttCollection, kTRUE);
  301. }
  302. // -------------------------------------------------------------------------
  303. // ----- Public method GetCollection -----------------------------------
  304. TClonesArray* CbmStt::GetCollection(Int_t iColl) const
  305. {
  306. if (iColl == 0)
  307. return fSttCollection;
  308. else
  309. return NULL;
  310. }
  311. // -------------------------------------------------------------------------
  312. // ----- Public method Print -------------------------------------------
  313. void CbmStt::Print() const
  314. {
  315. Int_t
  316. nHits = fSttCollection->GetEntriesFast();
  317. cout << "-I- CbmStt: " << nHits << " points registered in this event." << endl;
  318. if (fVerboseLevel>1)
  319. for (Int_t i=0; i<nHits; i++)
  320. (*fSttCollection)[i]->Print();
  321. }
  322. // -------------------------------------------------------------------------
  323. // ----- Public method Reset -------------------------------------------
  324. void CbmStt::Reset()
  325. {
  326. fSttCollection->Clear();
  327. ResetParameters();
  328. }
  329. // -------------------------------------------------------------------------
  330. // ----- Public method CopyClones --------------------------------------
  331. void CbmStt::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)
  332. {
  333. Int_t
  334. nEntries = cl1->GetEntriesFast();
  335. cout << "-I- CbmStt: " << nEntries << " entries to add." << endl;
  336. TClonesArray& clref = *cl2;
  337. CbmSttPoint
  338. *oldpoint = NULL;
  339. for (Int_t i=0; i<nEntries; i++)
  340. {
  341. oldpoint = (CbmSttPoint*) cl1->At(i);
  342. Int_t
  343. index = oldpoint->GetTrackID() + offset;
  344. oldpoint->SetTrackID(index);
  345. new (clref[fPosIndex]) CbmSttPoint(*oldpoint);
  346. fPosIndex++;
  347. }
  348. cout << "-I- CbmStt: " << cl2->GetEntriesFast() << " merged entries."
  349. << endl;
  350. }
  351. // -------------------------------------------------------------------------
  352. // ----- Public method ConstructGeometry -------------------------------
  353. void CbmStt::ConstructGeometry()
  354. {
  355. if (!InitProcessHits())
  356. {
  357. return;
  358. }
  359. // get pointer to the instantons which interface
  360. // to monte carlo
  361. FairGeoLoader
  362. *geoLoad = FairGeoLoader::Instance();
  363. FairGeoInterface
  364. *geoFace = geoLoad->getGeoInterface();
  365. CbmGeoStt
  366. *sttGeo = new CbmGeoStt();
  367. sttGeo->setGeomFile(GetGeometryFileName());
  368. geoFace->addGeoModule(sttGeo);
  369. Bool_t
  370. rc = geoFace->readSet(sttGeo);
  371. if (rc)
  372. sttGeo->create(geoLoad->getGeoBuilder());
  373. TList*
  374. volList = sttGeo->getListOfVolumes();
  375. FairRun *fRun = FairRun::Instance();
  376. // fRun = FairRun::Instance();
  377. TListIter
  378. iter(volList);
  379. FairGeoNode
  380. *node = NULL;
  381. FairGeoVolume
  382. *aVol = NULL;
  383. while( (node = (FairGeoNode*)iter.Next()) )
  384. {
  385. aVol = dynamic_cast<FairGeoVolume*> ( node );
  386. if ( node->isSensitive() )
  387. {
  388. if (fSensNodes)
  389. {
  390. fSensNodes->AddLast( aVol );
  391. }
  392. }
  393. else
  394. {
  395. if (fPassNodes)
  396. {
  397. fPassNodes->AddLast( aVol );
  398. }
  399. }
  400. }
  401. fParameters->setChanged();
  402. fParameters->setInputVersion(fRun->GetRunId(),1);
  403. ProcessNodes ( volList );
  404. }
  405. // -------------------------------------------------------------------------
  406. // ----- Private method AddHit -----------------------------------------
  407. CbmSttPoint* CbmStt::AddHit(Int_t trackID, Int_t detID, TVector3 pos,
  408. TVector3 posInLocal, TVector3 posOutLocal,
  409. TVector3 momIn, TVector3 momOut, TVector3 wireDir,
  410. Double_t time, Double_t length, Double_t eLoss, Double_t mass, TVector3 postot) // da cancellare postot
  411. {
  412. TClonesArray&
  413. clref = *fSttCollection;
  414. Int_t
  415. size = clref.GetEntriesFast();
  416. return new(clref[size]) CbmSttPoint(trackID, detID, pos, posInLocal, posOutLocal,
  417. momIn, momOut, wireDir, time, length, eLoss, mass, postot);
  418. }
  419. // -------------------------------------------------------------------------
  420. ClassImp(CbmStt)