From 0476dc7a29c8839b106d3989600d67d64ab87719 Mon Sep 17 00:00:00 2001 From: Jorge Lopez Date: Wed, 28 Nov 2018 22:52:33 +0000 Subject: [PATCH 1/6] WIP Improving hit creation Using what's on AliROOT as a template. --- Detectors/TRD/simulation/src/Detector.cxx | 60 +++++++++++++++++++++-- 1 file changed, 57 insertions(+), 3 deletions(-) diff --git a/Detectors/TRD/simulation/src/Detector.cxx b/Detectors/TRD/simulation/src/Detector.cxx index 5aa11561268a5..df1bc9196c8e5 100644 --- a/Detectors/TRD/simulation/src/Detector.cxx +++ b/Detectors/TRD/simulation/src/Detector.cxx @@ -59,6 +59,60 @@ bool Detector::ProcessHits(FairVolume* v) return false; } + // Inside sensitive volume ? + bool drRegion = false; + bool amRegion = false; + const TString cIdSensDr = "J"; + const TString cIdSensAm = "K"; + TString cIdCurrent = fMC->CurrentVolName(); + std::cout << "TRD::Detector::ProcessHits() \t cIdCurrent = " << cIdCurrent << std::endl; + if (cIdCurrent[1] == cIdSensDr) { + drRegion = true; + } + if (cIdCurrent[2] == cIdSensAm) { + amRegion = true; + } + + if (!drRegion && !amRegion) { + return false; + } + + TString cIdPath = gGeoManager->GetPath(); + std::cout << "TRD::Detector::ProcessHits() \t cIdPath = " << cIdPath << std::endl; + // Example of cIdPath to get IdSector + // /cave_1/B077_1/BSEGMO13_1/BTRD13_1/UTR3_1/UTS3_1/UTI3_1/UT29_1/UD29_1/UE29_1/UK29_1 + // ** + // cIdPath[21] = 1; + // cIdPath[22] = 3; // then + // cIdSector = 130; + + // The sector number (0 - 17), according to standard coordinate system + char cIdSector[3]; + cIdSector[0] = cIdPath[21]; + cIdSector[1] = cIdPath[22]; + cIdSector[2] = 0; + int sector = std::stoi(cIdSector); + + // The plane and chamber number + char cIdChamber[3]; + cIdChamber[0] = cIdCurrent[2]; + cIdChamber[1] = cIdCurrent[3]; + cIdChamber[2] = 0; + const int idChamber = mGeom->getDetectorSec(std::stoi(cIdSector)); //(std::stoi(cIdChamber) % (kNlayer*kNstack)); // getDetectorSec(int det) + const int det = mGeom->getDetector(mGeom->getLayer(idChamber), mGeom->getStack(idChamber), sector); + + // Special hits if track is entering + if (drRegion && fMC->IsTrackEntering()) { + // Create a track reference at the entrance of each + // chamber that contains the momentum components of the particle + + } + else if(amRegion && fMC->IsTrackExiting()) { + // Create a track reference at the exit of each + // chamber that contains the momentum components of the particle + + } + // just record position and basic quantities for the moment // TODO: needs to be interpreted properly float x, y, z; @@ -66,9 +120,9 @@ bool Detector::ProcessHits(FairVolume* v) float enDep = fMC->Edep(); float time = fMC->TrackTime() * 1.0e09; - auto stack = (o2::Data::Stack*)fMC->GetStack(); - auto trackID = stack->GetCurrentTrackNumber(); - auto sensID = v->getMCid(); + o2::Data::Stack* stack = (o2::Data::Stack*) fMC->GetStack(); + const int trackID = stack->GetCurrentTrackNumber(); + const int sensID = v->getMCid(); addHit(x, y, z, time, enDep, trackID, sensID); stack->addHit(GetDetId()); From f5396f49637924249c33d0db85a79ceebbea00f5 Mon Sep 17 00:00:00 2001 From: Jorge Lopez Date: Thu, 29 Nov 2018 13:23:45 +0000 Subject: [PATCH 2/6] [WIP] TRD hits update (O2-453) --- Detectors/TRD/simulation/src/Detector.cxx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Detectors/TRD/simulation/src/Detector.cxx b/Detectors/TRD/simulation/src/Detector.cxx index df1bc9196c8e5..730453721822c 100644 --- a/Detectors/TRD/simulation/src/Detector.cxx +++ b/Detectors/TRD/simulation/src/Detector.cxx @@ -65,7 +65,7 @@ bool Detector::ProcessHits(FairVolume* v) const TString cIdSensDr = "J"; const TString cIdSensAm = "K"; TString cIdCurrent = fMC->CurrentVolName(); - std::cout << "TRD::Detector::ProcessHits() \t cIdCurrent = " << cIdCurrent << std::endl; + // LOG(DEBUG) << "TRD::Detector::ProcessHits() \t cIdCurrent = " << cIdCurrent; if (cIdCurrent[1] == cIdSensDr) { drRegion = true; } @@ -78,7 +78,7 @@ bool Detector::ProcessHits(FairVolume* v) } TString cIdPath = gGeoManager->GetPath(); - std::cout << "TRD::Detector::ProcessHits() \t cIdPath = " << cIdPath << std::endl; + // LOG(DEBUG) << "TRD::Detector::ProcessHits() \t cIdPath = " << cIdPath; // Example of cIdPath to get IdSector // /cave_1/B077_1/BSEGMO13_1/BTRD13_1/UTR3_1/UTS3_1/UTI3_1/UT29_1/UD29_1/UE29_1/UK29_1 // ** @@ -91,14 +91,15 @@ bool Detector::ProcessHits(FairVolume* v) cIdSector[0] = cIdPath[21]; cIdSector[1] = cIdPath[22]; cIdSector[2] = 0; - int sector = std::stoi(cIdSector); - - // The plane and chamber number + // The plane and chamber number char cIdChamber[3]; cIdChamber[0] = cIdCurrent[2]; cIdChamber[1] = cIdCurrent[3]; cIdChamber[2] = 0; - const int idChamber = mGeom->getDetectorSec(std::stoi(cIdSector)); //(std::stoi(cIdChamber) % (kNlayer*kNstack)); // getDetectorSec(int det) + + int sector = std::stoi(cIdSector); + // const idChamber = (std::stoi(cIdChamber) % (kNlayer*kNstack)); + const int idChamber = mGeom->getDetectorSec(sector); const int det = mGeom->getDetector(mGeom->getLayer(idChamber), mGeom->getStack(idChamber), sector); // Special hits if track is entering From 75f4891f47d66174445ac3c88f07abb65c798c0e Mon Sep 17 00:00:00 2001 From: Jorge Lopez Date: Fri, 30 Nov 2018 16:24:46 +0000 Subject: [PATCH 3/6] [WIP] TRD hits update (O2-453) - Added a behaviour close to what is done in AliRoot - Still need to add TR. --- .../include/TRDSimulation/Detector.h | 3 + Detectors/TRD/simulation/src/Detector.cxx | 99 ++++++++++--------- 2 files changed, 58 insertions(+), 44 deletions(-) diff --git a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h index be402839c32f4..5b4441f6c544f 100644 --- a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h +++ b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h @@ -93,6 +93,9 @@ class Detector : public o2::Base::DetImpl float mGasNobleFraction; float mGasDensity; + bool mTRon; // Switch for TR simulation + int mWion; // Ionization potential + TRDGeometry* mGeom = nullptr; template diff --git a/Detectors/TRD/simulation/src/Detector.cxx b/Detectors/TRD/simulation/src/Detector.cxx index 730453721822c..0403d853d5874 100644 --- a/Detectors/TRD/simulation/src/Detector.cxx +++ b/Detectors/TRD/simulation/src/Detector.cxx @@ -34,8 +34,19 @@ Detector::Detector(const Detector& rhs) mFoilDensity(rhs.mFoilDensity), mGasNobleFraction(rhs.mGasNobleFraction), mGasDensity(rhs.mGasDensity), - mGeom(rhs.mGeom) + mGeom(rhs.mGeom), + mTRon(false) { + if (TRDCommonParam::Instance()->IsXenon()) { + mWion = 23.53; // Ionization energy XeCO2 (85/15) + } + else if (TRDCommonParam::Instance()->IsArgon()) { + mWion = 27.21; // Ionization energy ArCO2 (82/18) + } + else { + LOG(FATAL) << "Wrong gas mixture"; + // add hard exit here! + } } Detector::~Detector() @@ -52,7 +63,12 @@ void Detector::InitializeO2Detector() bool Detector::ProcessHits(FairVolume* v) { // very rudimentatary hit creation - // TODO: needs upgrade to the level of AliROOT + /* TODO: needs upgrade to the level of AliROOT + + - Add primary ionization (fluka, geant?) see AliRoot + - Add TR + + */ // If not charged track or already stopped or disappeared, just return. if ((!fMC->TrackCharge()) || fMC->IsTrackDisappeared()) { @@ -77,57 +93,52 @@ bool Detector::ProcessHits(FairVolume* v) return false; } - TString cIdPath = gGeoManager->GetPath(); - // LOG(DEBUG) << "TRD::Detector::ProcessHits() \t cIdPath = " << cIdPath; - // Example of cIdPath to get IdSector - // /cave_1/B077_1/BSEGMO13_1/BTRD13_1/UTR3_1/UTS3_1/UTI3_1/UT29_1/UD29_1/UE29_1/UK29_1 - // ** - // cIdPath[21] = 1; - // cIdPath[22] = 3; // then - // cIdSector = 130; - - // The sector number (0 - 17), according to standard coordinate system - char cIdSector[3]; - cIdSector[0] = cIdPath[21]; - cIdSector[1] = cIdPath[22]; - cIdSector[2] = 0; - // The plane and chamber number - char cIdChamber[3]; - cIdChamber[0] = cIdCurrent[2]; - cIdChamber[1] = cIdCurrent[3]; - cIdChamber[2] = 0; - - int sector = std::stoi(cIdSector); - // const idChamber = (std::stoi(cIdChamber) % (kNlayer*kNstack)); - const int idChamber = mGeom->getDetectorSec(sector); - const int det = mGeom->getDetector(mGeom->getLayer(idChamber), mGeom->getStack(idChamber), sector); - + // 0: InFlight 1: Entering 2: Exiting + int trkStat = 0; // Special hits if track is entering if (drRegion && fMC->IsTrackEntering()) { // Create a track reference at the entrance of each // chamber that contains the momentum components of the particle - + /* Do what AliRoot does here + fMC->TrackMomentum(mom); + AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD); + */ + // Create the hits from TR photons if electron/positron is entering the drift volume + const int ele = TMath::Abs(fMC->TrackPid()) == 11; // electron PDG code. + if (mTRon && ele) { + // Do TR simulation here + // Emulate what AliTRDsimTR does here + // CreateTRhit(); // See AliTRDv1.cxx + } + trkStat = 1; } else if(amRegion && fMC->IsTrackExiting()) { // Create a track reference at the exit of each // chamber that contains the momentum components of the particle - + /* Do what AliRoot does here + fMC->TrackMomentum(mom); + AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD); + trkStat = 2; + */ + trkStat = 2; } - - // just record position and basic quantities for the moment - // TODO: needs to be interpreted properly - float x, y, z; - fMC->TrackPosition(x, y, z); - - float enDep = fMC->Edep(); - float time = fMC->TrackTime() * 1.0e09; - o2::Data::Stack* stack = (o2::Data::Stack*) fMC->GetStack(); - const int trackID = stack->GetCurrentTrackNumber(); - const int sensID = v->getMCid(); - addHit(x, y, z, time, enDep, trackID, sensID); - stack->addHit(GetDetId()); - - return true; + // Calculate the charge according to GEANT Edep + // Create a new dEdx hit + const double enDep = fMC->Edep(); + // Store those hits with enDep bigger than the ionization potential of the gas mixture for in-flight tracks + // or store hits of tracks that are entering or exiting + if ((enDep > mWion) || trkStat) { + double x,y,z; + fMC->TrackPosition(x,y,z); + double time = fMC->TrackTime() * 1e9; + o2::Data::Stack* stack = (o2::Data::Stack*) fMC->GetStack(); + const int trackID = stack->GetCurrentTrackNumber(); + const int sensID = v->getMCid(); + addHit(x, y, z, time, enDep, trackID, sensID); + stack->addHit(GetDetId()); + return true; + } + return false; } void Detector::Register() From ae4976b1d55a7c9a0cce50038e1ac6b12c580392 Mon Sep 17 00:00:00 2001 From: Jorge Lopez Date: Fri, 30 Nov 2018 16:54:59 +0000 Subject: [PATCH 4/6] fix wrong type --- .../TRD/simulation/include/TRDSimulation/Detector.h | 2 +- Detectors/TRD/simulation/src/Detector.cxx | 13 +++++-------- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h index 5b4441f6c544f..c6d60a39bc496 100644 --- a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h +++ b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h @@ -94,7 +94,7 @@ class Detector : public o2::Base::DetImpl float mGasDensity; bool mTRon; // Switch for TR simulation - int mWion; // Ionization potential + float mWion; // Ionization potential TRDGeometry* mGeom = nullptr; diff --git a/Detectors/TRD/simulation/src/Detector.cxx b/Detectors/TRD/simulation/src/Detector.cxx index 0403d853d5874..1eb4279e5e38d 100644 --- a/Detectors/TRD/simulation/src/Detector.cxx +++ b/Detectors/TRD/simulation/src/Detector.cxx @@ -39,11 +39,9 @@ Detector::Detector(const Detector& rhs) { if (TRDCommonParam::Instance()->IsXenon()) { mWion = 23.53; // Ionization energy XeCO2 (85/15) - } - else if (TRDCommonParam::Instance()->IsArgon()) { + } else if (TRDCommonParam::Instance()->IsArgon()) { mWion = 27.21; // Ionization energy ArCO2 (82/18) - } - else { + } else { LOG(FATAL) << "Wrong gas mixture"; // add hard exit here! } @@ -78,8 +76,8 @@ bool Detector::ProcessHits(FairVolume* v) // Inside sensitive volume ? bool drRegion = false; bool amRegion = false; - const TString cIdSensDr = "J"; - const TString cIdSensAm = "K"; + const TString cIdSensDr = "J"; + const TString cIdSensAm = "K"; TString cIdCurrent = fMC->CurrentVolName(); // LOG(DEBUG) << "TRD::Detector::ProcessHits() \t cIdCurrent = " << cIdCurrent; if (cIdCurrent[1] == cIdSensDr) { @@ -111,8 +109,7 @@ bool Detector::ProcessHits(FairVolume* v) // CreateTRhit(); // See AliTRDv1.cxx } trkStat = 1; - } - else if(amRegion && fMC->IsTrackExiting()) { + } else if(amRegion && fMC->IsTrackExiting()) { // Create a track reference at the exit of each // chamber that contains the momentum components of the particle /* Do what AliRoot does here From 16f8d3eb1c4073869b927bae5db67533d6d00ee0 Mon Sep 17 00:00:00 2001 From: Jorge Lopez Date: Mon, 3 Dec 2018 10:22:31 +0000 Subject: [PATCH 5/6] [WIP] fix the formatting (O2-453) --- .../TRD/simulation/include/TRDSimulation/Detector.h | 6 +++--- Detectors/TRD/simulation/src/Detector.cxx | 12 ++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h index c6d60a39bc496..41a7e99a13a36 100644 --- a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h +++ b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h @@ -93,9 +93,9 @@ class Detector : public o2::Base::DetImpl float mGasNobleFraction; float mGasDensity; - bool mTRon; // Switch for TR simulation - float mWion; // Ionization potential - + bool mTRon; // Switch for TR simulation + float mWion; // Ionization potential + TRDGeometry* mGeom = nullptr; template diff --git a/Detectors/TRD/simulation/src/Detector.cxx b/Detectors/TRD/simulation/src/Detector.cxx index 1eb4279e5e38d..22c25acf66990 100644 --- a/Detectors/TRD/simulation/src/Detector.cxx +++ b/Detectors/TRD/simulation/src/Detector.cxx @@ -101,15 +101,15 @@ bool Detector::ProcessHits(FairVolume* v) fMC->TrackMomentum(mom); AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD); */ - // Create the hits from TR photons if electron/positron is entering the drift volume + // Create the hits from TR photons if electron/positron is entering the drift volume const int ele = TMath::Abs(fMC->TrackPid()) == 11; // electron PDG code. if (mTRon && ele) { // Do TR simulation here - // Emulate what AliTRDsimTR does here + // Emulate what AliTRDsimTR does here // CreateTRhit(); // See AliTRDv1.cxx } trkStat = 1; - } else if(amRegion && fMC->IsTrackExiting()) { + } else if (amRegion && fMC->IsTrackExiting()) { // Create a track reference at the exit of each // chamber that contains the momentum components of the particle /* Do what AliRoot does here @@ -125,10 +125,10 @@ bool Detector::ProcessHits(FairVolume* v) // Store those hits with enDep bigger than the ionization potential of the gas mixture for in-flight tracks // or store hits of tracks that are entering or exiting if ((enDep > mWion) || trkStat) { - double x,y,z; - fMC->TrackPosition(x,y,z); + double x, y, z; + fMC->TrackPosition(x, y, z); double time = fMC->TrackTime() * 1e9; - o2::Data::Stack* stack = (o2::Data::Stack*) fMC->GetStack(); + o2::Data::Stack* stack = (o2::Data::Stack*)fMC->GetStack(); const int trackID = stack->GetCurrentTrackNumber(); const int sensID = v->getMCid(); addHit(x, y, z, time, enDep, trackID, sensID); From d22183af6ef4df33c1e40ada13c784babbd93505 Mon Sep 17 00:00:00 2001 From: Jorge Lopez Date: Mon, 3 Dec 2018 11:42:31 +0000 Subject: [PATCH 6/6] [WIP] fix the formatting (O2-453) --- Detectors/TRD/simulation/include/TRDSimulation/Detector.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h index 41a7e99a13a36..d58e1fd64fe86 100644 --- a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h +++ b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h @@ -95,7 +95,7 @@ class Detector : public o2::Base::DetImpl bool mTRon; // Switch for TR simulation float mWion; // Ionization potential - + TRDGeometry* mGeom = nullptr; template