read_dst_example.C

Go to the documentation of this file.

// read_dst_example.C: loop over DST events and print basic information on vertices and global tracks
// works with simulated and experimental DSTs; MC information is read only if it is present in the file
void read_dst_example(TString dstFileName = "bmndst.root", Int_t nPrint = 3)
{
    // run-level information is stored once per file
    TFile* file = TFile::Open(dstFileName);
    DstRunHeader* runHeader = (DstRunHeader*)file->Get("DstRunHeader");
    cout << "Period " << runHeader->GetPeriodId() << ", run " << runHeader->GetRunId() << ": " << runHeader->GetBeam()
         << " + " << runHeader->GetTarget() << " at " << runHeader->GetBeamEnergy() << " AGeV, field scale "
         << runHeader->GetFieldScale() << endl;
 
    TTree* tree = (TTree*)file->Get("bmndata");
    tree->SetBranchStatus("*Cluster*", 0);   // skip the heavy cluster branches, they are not used here
 
    DstEventHeader* eventHeader = nullptr;
    CbmVertex* vertex = nullptr;
    TClonesArray *globalTracks = nullptr, *innerTracks = nullptr, *tof400Hits = nullptr;
    tree->SetBranchAddress("DstEventHeader.", &eventHeader);
    tree->SetBranchAddress("MpdVertex.", &vertex);
    tree->SetBranchAddress("BmnGlobalTrack", &globalTracks);
    tree->SetBranchAddress("StsVector", &innerTracks);
    tree->SetBranchAddress("BmnTof400Hit", &tof400Hits);
 
    // MC information exists only in simulated DSTs
    BmnMCInfoDst* mcInfo = nullptr;
    Bool_t isSimulation = (tree->GetBranch("BmnMCInfo.") != nullptr);
    if (isSimulation)
        tree->SetBranchAddress("BmnMCInfo.", &mcInfo);
 
    for (Long64_t iEvent = 0; iEvent < TMath::Min(tree->GetEntries(), (Long64_t)nPrint); iEvent++) {
        tree->GetEntry(iEvent);
        printf("Event %u: vertex (%.2f, %.2f, %.2f) cm from %d tracks, %d global tracks\n", eventHeader->GetEventId(),
               vertex->GetX(), vertex->GetY(), vertex->GetZ(), vertex->GetNTracks(), globalTracks->GetEntriesFast());
        if (vertex->GetNTracks() < 2)
            continue;
 
        for (Int_t iTrack = 0; iTrack < TMath::Min(globalTracks->GetEntriesFast(), 3); iTrack++) {
            BmnGlobalTrack* track = (BmnGlobalTrack*)globalTracks->UncheckedAt(iTrack);
            CbmStsTrack* inner = (CbmStsTrack*)innerTracks->UncheckedAt(track->GetGemTrackIndex());
 
            FairTrackParam* par = track->GetParamFirst();   // at the z of MpdVertex
            Double_t p = 1. / TMath::Abs(par->GetQp());     // momentum for |Z| = 1, rigidity p/|Z| otherwise
            Int_t charge = (par->GetQp() > 0) ? +1 : -1;
            Double_t pz = p / TMath::Sqrt(1. + par->GetTx() * par->GetTx() + par->GetTy() * par->GetTy());
            Double_t pt = pz * TMath::Sqrt(par->GetTx() * par->GetTx() + par->GetTy() * par->GetTy());
 
            printf("  track %d: q = %+d, p = %.3f, pt = %.3f GeV/c, hits = %d, in vertex = %d, DCA = %.2f cm", iTrack,
                   charge, p, pt, inner->GetNStsHits(), inner->FromVertex() ? 1 : 0, track->GetDCAInVertex());
            if (track->GetTof1HitIndex() >= 0) {
                BmnTofHit* hit = (BmnTofHit*)tof400Hits->UncheckedAt(track->GetTof1HitIndex());
                Double_t beta = hit->GetLength() / hit->GetTimeStamp() / 29.9792458;   // length in cm, time in ns
                printf(", TOF400: t = %.2f ns, beta = %.3f, m2 = %.3f", hit->GetTimeStamp(), beta, track->GetMass2(1));
            }
            // the matched MC track is stored in the inner track
            if (isSimulation && inner->GetTrkID() >= 0) {
                CbmMCTrack* mcTrack = (CbmMCTrack*)mcInfo->GetMCTracks()->UncheckedAt(inner->GetTrkID());
                printf(", MC: pdg = %d, p = %.3f", mcTrack->GetPdgCode(), mcTrack->GetP());
            }
            printf("\n");
        }
    }
}