123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532 |
- //------------------------------------------------------------------------------------------------------------------------
- // -------------------------------------------------------------------------
- // ----- MpdStrawECT source file -----
- // -------------------------------------------------------------------------
- #include <iostream>
- #include "TClonesArray.h"
- #include "TLorentzVector.h"
- #include "TParticle.h"
- #include "TVirtualMC.h"
- #include "TGeoManager.h"
- #include "TGeoMatrix.h"
- #include "FairGeoInterface.h"
- #include "FairGeoLoader.h"
- #include "FairGeoNode.h"
- #include "FairGeoRootBuilder.h"
- #include "MpdStack.h"
- #include "MpdStrawECTGeo.h"
- #include "FairRootManager.h"
- #include "MpdStrawECT.h"
- #include "MpdStrawECTPoint.h"
- #include "FairRuntimeDb.h"
- #include "MpdStrawECTGeoPar.h"
- #include "TObjArray.h"
- #include "FairRun.h"
- #include "FairVolume.h"
- #include "TMath.h"
- #include "TParticlePDG.h"
- #include "TGeoManager.h"
- #include "TGDMLParse.h"
- #include "FairGeoMedia.h"
- #include "TString.h"
- class FairVolume;
- //------------------------------------------------------------------------------------------------------------------------
- MpdStrawECT::MpdStrawECT()
- : FairDetector("strawECT", kTRUE)
- {
- fStrawECTCollection = new TClonesArray("MpdStrawECTPoint");
- fPosIndex = 0;
- fVerboseLevel = 1;
- }
- //------------------------------------------------------------------------------------------------------------------------
- MpdStrawECT::MpdStrawECT(const char* name, Bool_t active)
- : FairDetector(name, active)
- {
- fStrawECTCollection = new TClonesArray("MpdStrawECTPoint");
- fPosIndex = 0;
- fVerboseLevel = 1;
- }
- //------------------------------------------------------------------------------------------------------------------------
- MpdStrawECT::~MpdStrawECT()
- {
- if(fStrawECTCollection){ fStrawECTCollection->Delete(); delete fStrawECTCollection; }
- }
- //------------------------------------------------------------------------------------------------------------------------
- int MpdStrawECT::DistAndPoints(TVector3 p1, TVector3 p2, TVector3 p3, TVector3 p4, TVector3& pa, TVector3& pb) {
- TVector3 A = p2 - p1;
- TVector3 B = p4 - p3;
- TVector3 C = p1 - p3;
- if (p3 == p4) {
- pb = p4;
- TVector3 unit = A.Unit();
- Double_t dist = unit.Dot(-C);
- pa = p1 + dist * unit;
- return 0;
- }
- Double_t numer = C.Dot(B) * B.Dot(A) - C.Dot(A) * B.Dot(B);
- Double_t denom = A.Dot(A) * B.Dot(B) - B.Dot(A) * B.Dot(A);
- if (denom == 0) {
- // parallel lines
- pa.SetXYZ(-10000,-10000,-10000);
- pb.SetXYZ(-10000,-10000,-10000);
- return 1;
- }
- Double_t ma = numer / denom;
- Double_t mb = ( C.Dot(B) + B.Dot(A)*ma ) / B.Dot(B);
- pa = p1 + A*ma;
- pb = p3 + B*mb;
- return 0;
- }
- //------------------------------------------------------------------------------------------------------------------------
- TVector3 MpdStrawECT::GlobalToLocal(TVector3& global) {
- Double_t globPos[3];
- Double_t localPos[3];
- global.GetXYZ(globPos);
- gMC->Gmtod(globPos, localPos, 1);
- return TVector3(localPos);
- }
- //------------------------------------------------------------------------------------------------------------------------
- TVector3 MpdStrawECT::LocalToGlobal(TVector3& local) {
- Double_t globPos[3];
- Double_t localPos[3];
- local.GetXYZ(localPos);
- gMC->Gdtom(localPos, globPos, 1);
- return TVector3(globPos);
- }
- //----------------------------------------------------------------------------------------------------------------------
- Bool_t MpdStrawECT::ProcessHits(FairVolume* vol) {
- // Set parameters at entrance of volume. Reset ELoss.
- if(gMC->IsTrackEntering()) {
- ResetParameters();
- fELoss = 0.;
- fTime = gMC->TrackTime() * 1.0e09;
- fLength = gMC->TrackLength();
- fIsPrimary = 0;
- fCharge = -1;
- fPdgId = 0;
- TLorentzVector PosIn;
- gMC->TrackPosition(PosIn);
- fPosIn.SetXYZ(PosIn.X(), PosIn.Y(), PosIn.Z());
- gMC->TrackMomentum(fMom);
- TParticle* part = 0;
- part = gMC->GetStack()->GetCurrentTrack();
- if (part) {
- fIsPrimary = (Int_t)part->IsPrimary();
- fCharge = (Int_t)part->GetPDG()->Charge();
- fPdgId = (Int_t)part->GetPdgCode();
- }
- fVolumeID = vol->getMCid();
- fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
- }
- // Sum energy loss for all steps in the active volume
- fELoss += gMC->Edep();
- // Create MpdStrawECTPoint at EXIT of active volume;
- if ( (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) && fELoss > 0 ) {
- fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
- TLorentzVector PosOut;
- gMC->TrackPosition(PosOut);
- fPosOut.SetXYZ(PosOut.X(), PosOut.Y(), PosOut.Z());
- //-------- define id for volumes from vol. path -------------------------------
- TString vol_path = gMC->CurrentVolPath();
- Int_t module_num = -1;
- Int_t submodule_num = -1;
- Int_t layer_num = -1;
- TString layer_type = "";
- Int_t straw_num = -1;
- TString key_module = "module_";
- TString key_submodule = "sub_module_";
- TString key_layer_radial = "layer_radial_";
- TString key_layer_stereoR = "layer_stereoR_";
- TString key_layer_stereoL = "layer_stereoL_";
- TString key_straw_tube = "straw_tube_";
- TString key_delim = "/";
- Int_t current_pos = 0;
- Int_t pos_in = 0;
- Int_t pos_out = 0;
- //module number
- pos_in = vol_path.Index(key_module, current_pos);
- pos_in += key_module.Length();
- pos_out = vol_path.Index(key_delim, pos_in);
- current_pos = pos_out;
- module_num = ((TString)(vol_path(pos_in, pos_out-pos_in))).Atoi();
- //submodule number
- pos_in = vol_path.Index(key_submodule, current_pos);
- pos_in += key_submodule.Length();
- pos_out = vol_path.Index(key_delim, pos_in);
- current_pos = pos_out;
- submodule_num = ((TString)(vol_path(pos_in, pos_out-pos_in))).Atoi();
- //layer number
- pos_in = vol_path.Index(key_layer_radial, current_pos);
- if(pos_in != -1) {
- pos_in += key_layer_radial.Length();
- layer_type = "radial";
- }
- if(pos_in == -1) {
- pos_in = vol_path.Index(key_layer_stereoR, current_pos);
- if(pos_in != -1) {
- pos_in += key_layer_stereoR.Length();
- layer_type = "stereoR";
- }
- }
- if(pos_in == -1) {
- pos_in = vol_path.Index(key_layer_stereoL, current_pos);
- if(pos_in != -1) {
- pos_in += key_layer_stereoL.Length();
- layer_type = "stereoL";
- }
- }
- if(pos_in != -1) {
- pos_out = vol_path.Index(key_delim, pos_in);
- current_pos = pos_out;
- layer_num = ((TString)(vol_path(pos_in, pos_out-pos_in))).Atoi();
- }
- //straw number
- pos_in = vol_path.Index(key_straw_tube, current_pos);
- pos_in += key_straw_tube.Length();
- pos_out = vol_path.Index(key_delim, pos_in);
- current_pos = pos_out;
- straw_num = ((TString)(vol_path(pos_in, pos_out-pos_in))).Atoi();
- //cout << "vol_path = " << vol_path << "\n";
- //cout << "module_num = " << module_num << "\n";
- //cout << "submodule_num = " << submodule_num << "\n";
- //cout << "layer_type = " << layer_type << "\n";
- //cout << "layer_num = " << layer_num << "\n";
- //cout << "straw_num = " << straw_num << "\n";
- //------------------------------------------------------------------------------
- MpdStrawECTPoint *p =
- AddHit(fTrackID, fVolumeID, fPosIn, TVector3(fMom.Px(), fMom.Py(), fMom.Pz()),
- fTime, (fLength+gMC->TrackLength())/2, fELoss);
- p->SetModule(module_num);
- p->SetSubmodule(submodule_num);
- p->SetLayer(layer_num);
- p->SetLayerType(layer_type);
- p->SetStraw(straw_num);
- ((MpdStack*)gMC->GetStack())->AddPoint(kECT);
- }
- return kTRUE;
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdStrawECT::EndOfEvent()
- {
- if(fVerboseLevel) Print();
- fStrawECTCollection->Delete();
- fPosIndex = 0;
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdStrawECT::Register(){ FairRootManager::Instance()->Register("StrawECTPoint", "StrawECT", fStrawECTCollection, kTRUE); }
- //------------------------------------------------------------------------------------------------------------------------
- TClonesArray* MpdStrawECT::GetCollection(Int_t iColl) const {
- if(iColl == 0) return fStrawECTCollection;
- return NULL;
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdStrawECT::Print() const
- {
- Int_t nHits = fStrawECTCollection->GetEntriesFast();
- cout << "-I- MpdStrawECT: " << nHits << " points registered in this event." << endl;
- if(fVerboseLevel > 1)
- for(Int_t i=0; i<nHits; i++) (*fStrawECTCollection)[i]->Print();
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdStrawECT::Reset(){ fStrawECTCollection->Delete(); ResetParameters(); }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdStrawECT::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)
- {
- Int_t nEntries = cl1->GetEntriesFast();
- cout << "-I- MpdStrawECT: " << nEntries << " entries to add." << endl;
- TClonesArray& clref = *cl2;
- MpdStrawECTPoint* oldpoint = NULL;
- for(Int_t i=0; i<nEntries; i++) {
- oldpoint = (MpdStrawECTPoint*) cl1->At(i);
- Int_t index = oldpoint->GetTrackID() + offset;
- oldpoint->SetTrackID(index);
- new (clref[fPosIndex]) MpdStrawECTPoint(*oldpoint);
- fPosIndex++;
- }
- cout << "-I- MpdStrawECT: " << cl2->GetEntriesFast() << " merged entries." << endl;
- }
- //------------------------------------------------------------------------------------------------------------------------
- void MpdStrawECT::ConstructGeometry() {
- TString fileName = GetGeometryFileName();
- if ( fileName.EndsWith(".root") ) {
- gLogger->Info(MESSAGE_ORIGIN,
- "Constructing straw ECT geometry from ROOT file %s",
- fileName.Data());
- ConstructRootGeometry();
- }
- else if ( fileName.EndsWith(".geo") ) {
- gLogger->Info(MESSAGE_ORIGIN,
- "Constructing straw ECT geometry from ASCII file %s",
- fileName.Data());
- ConstructAsciiGeometry();
- }
- else if ( fileName.EndsWith(".gdml") ) {
- gLogger->Info(MESSAGE_ORIGIN,
- "Constructing straw ECT geometry from GDML file %s",
- fileName.Data());
- //ConstructGDMLGeometry();
- }
- else {
- gLogger->Fatal(MESSAGE_ORIGIN,
- "Geometry format of straw ECT file %s not supported.",
- fileName.Data());
- }
- }
- // ----- ConstructAsciiGeometry -------------------------------------------
- void MpdStrawECT::ConstructAsciiGeometry() {
- Int_t count=0;
- Int_t count_tot=0;
- FairGeoLoader* geoLoad = FairGeoLoader::Instance();
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- MpdStrawECTGeo* sttGeo = new MpdStrawECTGeo();
- sttGeo->setGeomFile(GetGeometryFileName());
- geoFace->addGeoModule(sttGeo);
- Bool_t rc = geoFace->readSet(sttGeo);
- if(rc) sttGeo->create(geoLoad->getGeoBuilder());
- TList* volList = sttGeo->getListOfVolumes();
- // store geo parameter
- FairRun *fRun = FairRun::Instance();
- FairRuntimeDb *rtdb = FairRun::Instance()->GetRuntimeDb();
- MpdStrawECTGeoPar* par =(MpdStrawECTGeoPar*)(rtdb->getContainer("MpdStrawECTGeoPar"));
- TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
- TObjArray *fPassNodes = par->GetGeoPassiveNodes();
- FairGeoNode *node = NULL;
- FairGeoVolume *aVol = NULL;
- TListIter iter(volList);
- while((node = (FairGeoNode*)iter.Next())) {
- aVol = dynamic_cast<FairGeoVolume*> (node);
- if(node->isSensitive()){ fSensNodes->AddLast(aVol); count++; }
- else fPassNodes->AddLast(aVol);
- count_tot++;
- }
- par->setChanged();
- par->setInputVersion(fRun->GetRunId(), 1);
- ProcessNodes(volList);
- }
- // ----------------------------------------------------------------------------
- // ----- ConstructGDMLGeometry -------------------------------------------
- void MpdStrawECT::ConstructGDMLGeometry()
- {
- TFile *old = gFile;
- TGDMLParse parser;
- TGeoVolume* gdmlTop;
- // Before importing GDML
- Int_t maxInd = gGeoManager->GetListOfMedia()->GetEntries() - 1;
- gdmlTop = parser.GDMLReadFile(GetGeometryFileName());
- // Cheating - reassigning media indices after GDML import (need to fix this in TGDMLParse class!!!)
- // for (Int_t i=0; i<gGeoManager->GetListOfMedia()->GetEntries(); i++)
- // gGeoManager->GetListOfMedia()->At(i)->Dump();
- // After importing GDML
- Int_t j = gGeoManager->GetListOfMedia()->GetEntries() - 1;
- Int_t curId;
- TGeoMedium* m;
- do {
- m = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(j);
- curId = m->GetId();
- m->SetId(curId+maxInd);
- j--;
- } while (curId > 1);
- // LOG(DEBUG) << "====================================================================" << FairLogger::endl;
- // for (Int_t i=0; i<gGeoManager->GetListOfMedia()->GetEntries(); i++)
- // gGeoManager->GetListOfMedia()->At(i)->Dump();
- Int_t newMaxInd = gGeoManager->GetListOfMedia()->GetEntries() - 1;
- gGeoManager->GetTopVolume()->AddNode(gdmlTop, 1, 0);
- ExpandNodeForGdml(gGeoManager->GetTopVolume()->GetNode(gGeoManager->GetTopVolume()->GetNdaughters()-1));
- for (Int_t k = maxInd+1; k < newMaxInd+1; k++) {
- TGeoMedium* medToDel = (TGeoMedium*)(gGeoManager->GetListOfMedia()->At(maxInd+1));
- LOG(DEBUG) << " removing media " << medToDel->GetName() << " with id " << medToDel->GetId() << " (k=" << k << ")" << FairLogger::endl;
- gGeoManager->GetListOfMedia()->Remove(medToDel);
- }
- gGeoManager->SetAllIndex();
- gFile = old;
- }
- void MpdStrawECT::ExpandNodeForGdml(TGeoNode* node)
- {
- LOG(DEBUG) << "----------------------------------------- ExpandNodeForGdml for node " << node->GetName() << FairLogger::endl;
- TGeoVolume* curVol = node->GetVolume();
- LOG(DEBUG) << " volume: " << curVol->GetName() << FairLogger::endl;
- if (curVol->IsAssembly()) {
- LOG(DEBUG) << " skipping volume-assembly" << FairLogger::endl;
- }
- else {
- TGeoMedium* curMed = curVol->GetMedium();
- TGeoMaterial* curMat = curVol->GetMaterial();
- TGeoMedium* curMedInGeoManager = gGeoManager->GetMedium(curMed->GetName());
- TGeoMaterial* curMatOfMedInGeoManager = curMedInGeoManager->GetMaterial();
- TGeoMaterial* curMatInGeoManager = gGeoManager->GetMaterial(curMat->GetName());
- // Current medium and material assigned to the volume from GDML
- LOG(DEBUG2) << " curMed\t\t\t\t" << curMed << "\t" << curMed->GetName() << "\t" << curMed->GetId() << FairLogger::endl;
- LOG(DEBUG2) << " curMat\t\t\t\t" << curMat << "\t" << curMat->GetName() << "\t" << curMat->GetIndex() << FairLogger::endl;
- // Medium and material found in the gGeoManager - either the pre-loaded one or one from GDML
- LOG(DEBUG2) << " curMedInGeoManager\t\t" << curMedInGeoManager
- << "\t" << curMedInGeoManager->GetName() << "\t" << curMedInGeoManager->GetId() << FairLogger::endl;
- LOG(DEBUG2) << " curMatOfMedInGeoManager\t\t" << curMatOfMedInGeoManager
- << "\t" << curMatOfMedInGeoManager->GetName() << "\t" << curMatOfMedInGeoManager->GetIndex() << FairLogger::endl;
- LOG(DEBUG2) << " curMatInGeoManager\t\t" << curMatInGeoManager
- << "\t" << curMatInGeoManager->GetName() << "\t" << curMatInGeoManager->GetIndex() << FairLogger::endl;
- TString matName = curMat->GetName();
- TString medName = curMed->GetName();
- if (curMed->GetId() != curMedInGeoManager->GetId()) {
- if (fFixedMedia.find(medName) == fFixedMedia.end()) {
- LOG(DEBUG) << " Medium needs to be fixed" << FairLogger::endl;
- fFixedMedia[medName] = curMedInGeoManager;
- Int_t ind = curMat->GetIndex();
- gGeoManager->RemoveMaterial(ind);
- LOG(DEBUG) << " removing material " << curMat->GetName()
- << " with index " << ind << FairLogger::endl;
- for (Int_t i=ind; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++) {
- TGeoMaterial* m = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(i);
- m->SetIndex(m->GetIndex()-1);
- }
- LOG(DEBUG) << " Medium fixed" << FairLogger::endl;
- }
- else {
- LOG(DEBUG) << " Already fixed medium found in the list " << FairLogger::endl;
- }
- }
- else {
- if (fFixedMedia.find(medName) == fFixedMedia.end()) {
- LOG(DEBUG) << " There is no correct medium in the memory yet" << FairLogger::endl;
- FairGeoLoader* geoLoad = FairGeoLoader::Instance();
- FairGeoInterface* geoFace = geoLoad->getGeoInterface();
- FairGeoMedia* geoMediaBase = geoFace->getMedia();
- FairGeoBuilder* geobuild = geoLoad->getGeoBuilder();
- FairGeoMedium* curMedInGeo = geoMediaBase->getMedium(medName);
- if (curMedInGeo == 0) {
- LOG(FATAL) << " Media not found in Geo file: " << medName << FairLogger::endl;
- //! This should not happen.
- //! This means that somebody uses material in GDML that is not in the media.geo file.
- //! Most probably this is the sign to the user to check materials' names in the CATIA model.
- }
- else {
- LOG(DEBUG) << " Found media in Geo file" << medName << FairLogger::endl;
- Int_t nmed = geobuild->createMedium(curMedInGeo);
- fFixedMedia[medName] = (TGeoMedium*)gGeoManager->GetListOfMedia()->Last();
- gGeoManager->RemoveMaterial(curMatOfMedInGeoManager->GetIndex());
- LOG(DEBUG) << " removing material " << curMatOfMedInGeoManager->GetName()
- << " with index " << curMatOfMedInGeoManager->GetIndex() << FairLogger::endl;
- for (Int_t i=curMatOfMedInGeoManager->GetIndex(); i<gGeoManager->GetListOfMaterials()->GetEntries(); i++) {
- TGeoMaterial* m = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(i);
- m->SetIndex(m->GetIndex()-1);
- }
- }
- if (curMedInGeo->getSensitivityFlag()) {
- LOG(DEBUG) << " Adding sensitive " << curVol->GetName() << FairLogger::endl;
- AddSensitiveVolume(curVol);
- }
- }
- else {
- LOG(DEBUG) << " Already fixed medium found in the list" << FairLogger::endl;
- LOG(DEBUG) << "!!! Sensitivity: " << fFixedMedia[medName]->GetParam(0) << FairLogger::endl;
- if (fFixedMedia[medName]->GetParam(0) == 1) {
- LOG(DEBUG) << " Adding sensitive " << curVol->GetName() << FairLogger::endl;
- AddSensitiveVolume(curVol);
- }
- }
- }
- curVol->SetMedium(fFixedMedia[medName]);
- gGeoManager->SetAllIndex();
- // gGeoManager->GetListOfMaterials()->Print();
- // gGeoManager->GetListOfMedia()->Print();
- }
- //! Recursevly go down the tree of nodes
- if (curVol->GetNdaughters() != 0) {
- TObjArray* NodeChildList = curVol->GetNodes();
- TGeoNode* curNodeChild;
- for (Int_t j=0; j<NodeChildList->GetEntriesFast(); j++) {
- curNodeChild = (TGeoNode*)NodeChildList->At(j);
- ExpandNodeForGdml(curNodeChild);
- }
- }
- }
- //-----------------------------------------------------------------------------
- //Check if Sensitive-----------------------------------------------------------
- Bool_t MpdStrawECT::CheckIfSensitive(std::string name) {
- TString tsname = name;
- if (tsname.Contains("straw_gas")) {
- return kTRUE;
- }
- return kFALSE;
- }
- //-----------------------------------------------------------------------------
- MpdStrawECTPoint* MpdStrawECT::AddHit(Int_t trackID, Int_t detID, TVector3 pos,
- TVector3 mom, Double_t time, Double_t length, Double_t eLoss) {
- TClonesArray& clref = *fStrawECTCollection;
- Int_t size = clref.GetEntriesFast();
- //std::cout << "ELoss: " << eLoss << "\n";
- return new(clref[size]) MpdStrawECTPoint(trackID, detID, pos, mom, time, length, eLoss);
- }
- //------------------------------------------------------------------------------------------------------------------------
- ClassImp(MpdStrawECT)
|