Using files simulated with edep-sim

This wiki will explain how to run GArSoft on files that where simulated with edep-sim

Type of files output by edep-sim

edep-sim is using its own data model. Different classes handles the object access and the streaming IO into ROOT. The classes are available here:
GArSoft data model is different from edep-sim and therefore in order to run the digitization and reconstruction on the simulated files, one needs to convert the edep-sim data model into GArSoft simulation data products.

edep-sim converter module

A converter module, available here, has been written in order to create simulation data product in GArSoft. It can be configured using this fhicl file for single events files or this fhicl file for overlayed spill events.

Different parameters needs to be changed

  • physics.producers.edepconvert.EDepSimFile: "" $\
ightarrow$ needs the path to the edep-sim output file
  • physics.producers.edepconvert.GhepFile: "" ["$\\rightarrow$"] if events generated with genie, needs the path to the ghep genie file
  • physics.producers.edepconvert.OverlayFile: false ["$\\rightarrow$"] flag if the file is containing overlayed events
  • services.Geometry.GDML: "nd_hall_mpd_only/nd_hall_mpd_only_ECalOctagon_60l_UniformMagnet.gdml" ["$\\rightarrow$"] path to the geometry gdml file

The module will create the MCTruth and MCParticle data products in each event, as well it is look for the sensitive detector names and create data products for the TPC and the ECAL. Optionally the muon detector is present.
Here is an example of how it is done. This is a bit hardcoded and would be great if a solution could be found to avoid this.

        for (auto d = fEvent->SegmentDetectors.begin(); d != fEvent->SegmentDetectors.end(); ++d)
            if( d->first == "TPC_Drift1" || d->first == "TPC_Drift2" )
                //GAr deposits
                for (std::vector<TG4HitSegment>::const_iterator h = d->second.begin(); h != d->second.end(); ++h)
                    const TG4HitSegment *hit = &(*h);
                    int trackID = hit->GetPrimaryId();
                    double edep = hit->GetEnergyDeposit() * CLHEP::MeV / CLHEP::GeV;
                    double time = (hit->GetStart().T() + hit->GetStop().T())/2 / CLHEP::s;
                    double x = (hit->GetStart().X() + hit->GetStop().X())/2 /CLHEP::cm;
                    double y = (hit->GetStart().Y() + hit->GetStop().Y())/2 /CLHEP::cm;
                    double z = (hit->GetStart().Z() + hit->GetStop().Z())/2 /CLHEP::cm;
                    double stepLength = hit->GetTrackLength() /CLHEP::cm;
                    if(edep <= 0 || edep < fEnergyCut || stepLength <= 0)
                    TGeoNode *node = fGeo->FindNode(x, y, z);//Node in cm...
                    std::string VolumeName  = node->GetVolume()->GetName();
                    std::string volmaterial = node->GetMedium()->GetMaterial()->GetName();
                    if ( ! std::regex_match(volmaterial, std::regex(fTPCMaterial)) ) continue;
                    fGArDeposits.emplace_back(trackID, time, edep, x, y, z, stepLength, (trackID > 0));

To run the module, one needs the command

art -c conversionedepjob.fcl -n ${NEVENTS} -o ${filename_out}

The number of events is important or it will only run on 1 event. A couple of files are attached to this wiki page to provide an example of a bash script to run the conversion.