read_digi_example.C

Go to the documentation of this file.

// read_digi_example.C: loop over the events of a DIGI file and print the event header, the GEM digits per station,
// the FHCal signal and, for simulated files, the MC particle of the first GEM digit
// works with simulated and experimental DIGI files: the branch names are chosen by the type of the file
void read_digi_example(TString digiFileName = "bmnsim.root", Int_t nPrint = 3)
{
    TFile* file = TFile::Open(digiFileName);
    TTree* tree = (TTree*)file->Get("bmndata");
    Bool_t isSimulation = (tree->GetBranch("MCTrack") != nullptr);
 
    // run information is stored only in experimental files
    DigiRunHeader* runHeader = (DigiRunHeader*)file->Get("DigiRunHeader");
    if (runHeader)
        cout << "Period " << runHeader->GetPeriodId() << ", run " << runHeader->GetRunId() << endl;
 
    // the same detector data are stored in branches with different names in simulated and experimental files
    TString gemName = isSimulation ? "BmnGemStripDigit" : "GEM";
    TString fhcalName = isSimulation ? "FHCalDigit" : "FHCalDigi";
    vector<TString> branches = {gemName, fhcalName};
    if (isSimulation)
        branches.insert(branches.end(), {"MCEventHeader.", "BmnGemStripDigitMatch", "StsPoint", "MCTrack"});
    else
        branches.push_back("BmnEventHeader.");
    tree->SetBranchStatus("*", 0);   // read only the needed branches
    for (auto& name : branches)
        tree->SetBranchStatus(name + "*", 1);
 
    TClonesArray *gemDigits = nullptr, *fhcalDigits = nullptr;
    tree->SetBranchAddress(gemName, &gemDigits);
    tree->SetBranchAddress(fhcalName, &fhcalDigits);
    BmnEventHeader* eventHeader = nullptr;   // experimental files
    FairMCEventHeader* mcHeader = nullptr;   // simulated files
    TClonesArray *gemMatches = nullptr, *gemPoints = nullptr, *mcTracks = nullptr;
    if (isSimulation) {
        tree->SetBranchAddress("MCEventHeader.", &mcHeader);
        tree->SetBranchAddress("BmnGemStripDigitMatch", &gemMatches);
        tree->SetBranchAddress("StsPoint", &gemPoints);   // MC points of the GEM detectors
        tree->SetBranchAddress("MCTrack", &mcTracks);
    } else
        tree->SetBranchAddress("BmnEventHeader.", &eventHeader);
 
    for (Long64_t iEvent = 0; iEvent < TMath::Min(tree->GetEntries(), (Long64_t)nPrint); iEvent++) {
        tree->GetEntry(iEvent);
        if (isSimulation)
            printf("Event %u (b = %.1f fm):", mcHeader->GetEventID(), mcHeader->GetB());
        else
            printf("Event %u (%s):", eventHeader->GetEventId(), eventHeader->GetEventTimeTS().AsString("s"));
 
        // GEM: number of strip digits per station
        map<Int_t, Int_t> nGemDigits;
        for (Int_t i = 0; i < gemDigits->GetEntriesFast(); i++)
            nGemDigits[((BmnStripDigit*)gemDigits->UncheckedAt(i))->GetStation()]++;
        printf(" GEM digits per station:");
        for (auto& [station, n] : nGemDigits)
            printf(" %d", n);
 
        // FHCal: the digit classes of both file types derive from BmnFHCalDigit
        Double_t fhcalSignal = 0;
        for (Int_t i = 0; i < fhcalDigits->GetEntriesFast(); i++)
            fhcalSignal += ((BmnFHCalDigit*)fhcalDigits->UncheckedAt(i))->GetSignal();
        printf("; FHCal: %d digits, total signal %.1f", fhcalDigits->GetEntriesFast(), fhcalSignal);
 
        // simulated files: MC particle with the largest contribution to the first GEM digit
        if (isSimulation && gemMatches->GetEntriesFast() > 0) {
            BmnMatch* match = (BmnMatch*)gemMatches->UncheckedAt(0);
            if (match->GetNofLinks() > 0) {
                FairMCPoint* point = (FairMCPoint*)gemPoints->UncheckedAt(match->GetMatchedLink().GetIndex());
                CbmMCTrack* track = (CbmMCTrack*)mcTracks->UncheckedAt(point->GetTrackID());
                printf("; first GEM digit: PDG %d", track->GetPdgCode());
            }
        }
        printf("\n");
    }
}