From 905081dd3f4f6b060040ff133905a72bbc17b535 Mon Sep 17 00:00:00 2001 From: lialavezzi Date: Wed, 21 Dec 2022 11:20:29 +0100 Subject: [PATCH] in converter/convertHits.cc: added the conversion of MC primary tracks to EDM4Hep. In install_standalone.sh: changed the wget file location and set the ROME revision to download, to fix installation problems. In converter/convertTracks.cc fixed bug related to skipped vector, added the exception catching --- converter/convertHits.cc | 57 ++++++++++++++++++++++++++++++++++---- converter/convertTracks.cc | 12 +++++++- install_standalone.sh | 6 ++-- 3 files changed, 66 insertions(+), 9 deletions(-) diff --git a/converter/convertHits.cc b/converter/convertHits.cc index 9a69f8b..c5f80ae 100755 --- a/converter/convertHits.cc +++ b/converter/convertHits.cc @@ -28,6 +28,7 @@ #include "podio/EventStore.h" #include "podio/ROOTWriter.h" #include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/MCParticleCollection.h" // using namespace std; @@ -54,13 +55,14 @@ int main(int argc,char** argv) std::vector *hitspshw = new std::vector(); std::vector *hitsphcv = new std::vector(); std::vector *hitsphcvrd = new std::vector(); - std::vector *tracks = new std::vector(); + std::vector * primtracks = new std::vector(); bool hitChIsPresent=false; bool hitPxIsPresent=false; bool hitSVXIsPresent=false; bool hitPSHWIsPresent=false; bool hitPHCVIsPresent=false; bool hitPHCVRadIsPresent=false; + bool primtrackIsPresent=false; fo.GetListOfKeys()->Print(); @@ -109,12 +111,17 @@ int main(int argc,char** argv) std::cout << "Number of events: " << hit_tree->GetEntries() << std::endl; } - + if (br2.CompareTo(key->GetName()) == 0) { fo.GetObject(key->GetName(), track_tree); - track_tree->SetBranchAddress("Tracks",&tracks); + if (track_tree->FindBranch("Tracks")!=0x0) { + primtrackIsPresent=true; + track_tree->SetBranchAddress("Tracks",&primtracks); + std::cout<<"Found MCParticle tracks"<GetName() << std::endl; std::cout << "Number of events: " << track_tree->GetEntries() << std::endl; + } } @@ -124,6 +131,7 @@ int main(int argc,char** argv) std::cout << "PSHW hit " << hitPSHWIsPresent << std::endl; std::cout << "PHCV hit " << hitPHCVIsPresent << std::endl; std::cout << "PHCV hit " << hitPHCVRadIsPresent << std::endl; + std::cout << "MC primary track " << primtrackIsPresent << std::endl; @@ -163,14 +171,18 @@ int main(int argc,char** argv) edm4hep::SimTrackerHitCollection *s_pshwtrackerHits = new edm4hep::SimTrackerHitCollection(); l_evtstore->registerCollection("S_PSHWtrackerHits",s_pshwtrackerHits); l_writer->registerForWrite("S_PSHWtrackerHits"); + // MC + edm4hep::MCParticleCollection *s_mcparticleprimTracks = new edm4hep::MCParticleCollection(); + l_evtstore->registerCollection("S_MCParticlePrimTracks",s_mcparticleprimTracks); + l_writer->registerForWrite("S_MCParticlePrimTracks"); + - // event loop ------------------- int nevt = hit_tree->GetEntries(); std::cout << "nof events " << nevt << std::endl; for(int ievt=0; ievtGetEntry(ievt); + track_tree->GetEntry(ievt); hit_tree->GetEntry(ievt); // DCH --------------------------------- @@ -359,6 +371,41 @@ int main(int argc,char** argv) } } + // MC tracks --------------------------------------- + if (primtrackIsPresent) { + int ntracks = primtracks->size(); + std::cout << "event " << ievt << " has nof MC tracks " << ntracks << std::endl; + + // loop on MC tracks + for (int itrack=0; itrack < ntracks; itrack++) { + G4int s_parentid = primtracks->at(itrack)->GetParentID(); + if(s_parentid!=0) continue; // CHECK only primaries! + + // Get methods + G4int s_trackid = primtracks->at(itrack)->GetTrackID(); + G4String s_name = primtracks->at(itrack)->GetParticleName(); + G4double s_charge = primtracks->at(itrack)->GetParticleCharge(); + G4int s_pdg = primtracks->at(itrack)->GetPDGCode(); + G4ThreeVector s_vertex = primtracks->at(itrack)->GetPosStart(); + G4ThreeVector s_end = primtracks->at(itrack)->GetPosEnd(); + G4ThreeVector s_momentum = primtracks->at(itrack)->GetMomentum(); + + // convert to EDM DCH hit ............................ !! CHECK all the variables and units !! + auto l_track = s_mcparticleprimTracks->create(); + + l_track.setPDG(s_pdg); //PDG code of the particle + l_track.setCharge(s_charge); //particle charge + //production vertex of the particle in [mm] + edm4hep::Vector3d vertex(s_vertex.x(), s_vertex.y(), s_vertex.z()); + l_track.setVertex(vertex); + //endpoint of the particle in [mm] + edm4hep::Vector3d endpoint(s_end.x(), s_end.y(), s_end.z()); + l_track.setEndpoint(endpoint); + //particle 3-momentum at the production vertex in [GeV] + edm4hep::Vector3f momentum(s_momentum.x(), s_momentum.y(), s_momentum.z()); + l_track.setMomentum(momentum); + } + } // for each event write output if (l_writer != NULL) l_writer->writeEvent(); diff --git a/converter/convertTracks.cc b/converter/convertTracks.cc index 33e6ddf..ab82f22 100755 --- a/converter/convertTracks.cc +++ b/converter/convertTracks.cc @@ -277,7 +277,17 @@ int main(int argc,char** argv) // radius @ innermost state double tmpradius = 10000; // initialized to extremely big value by choice [mm] for(int ihit = 0; ihit < statevector->GetEntries(); ihit++) { - if(skipped.at(ihit)) continue; + + // CHECK fixed bug 2022-12-02 + // if(skipped.at(ihit)) continue; + try { + bool isthere = skipped.at(ihit); + if(isthere==true) continue; + } + catch(std::out_of_range) { + std::cout << "event " << ievt << "out of range error because skipped size is " << skipped.size() << std::endl; + } + TVector3 *state = (TVector3*) statevector->At(ihit); if(state==NULL) continue; // std::cout << state << " " << counter << " " << "hit " << hitindex.at(ihit) << " detector id " << detid.at(ihit) << " is skipped " << skipped.at(ihit) diff --git a/install_standalone.sh b/install_standalone.sh index 48b1cb2..54e93d5 100644 --- a/install_standalone.sh +++ b/install_standalone.sh @@ -58,9 +58,9 @@ cd $STANDALONE_INSTALL_DIR/LOCAL/ROME echo "download the ROME release" git clone https://bitbucket.org/muegamma/rome3.git rome cd $STANDALONE_INSTALL_DIR/LOCAL/ROME/rome -git checkout master +#git checkout master #git tag -#git checkout v3.2.15.1 +git checkout v3.2.15.1 export ROMESYS=$STANDALONE_INSTALL_DIR/LOCAL/ROME/rome export PATH=$ROMESYS/bin:$PATH cd $ROMESYS @@ -98,7 +98,7 @@ $ROMESYS/bin/romebuilder.exe -i GMC.xml #### Getting the gdml file from the ideadr box -wget https://cernbox.cern.ch/index.php/s/KxGYRFnkcob09z1/download -O g4-IDEA_reco.gdml +wget https://cernbox.cern.ch/remote.php/dav/public-files/KxGYRFnkcob09z1/g4-IDEA_reco.gdml # set the right PATH in the xml files string1="path_to_simulation"