O2DPG/MC/config/PWGDQ/ini/tests/Generator_InjectedInclusiveJpsiMidy_Pythia8_TriggerGap.C at 46715b61e890de10fe08071e0ab86dc62bf1ca32 · shreyasiacharya/O2DPG · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
int External()
{
    int checkPdgSignal[] = {443, 100443};
    int checkPdgDecay = 11;
    double rapiditymin = -4.3;
    double rapiditymax = -2.3;
    std::string path{"o2sim_Kine.root"};
    std::cout << "Check for\nsignal PDG " << checkPdgSignal << "\ndecay PDG " << checkPdgDecay << "\n";
    TFile file(path.c_str(), "READ");
    if (file.IsZombie())
    {
        std::cerr << "Cannot open ROOT file " << path << "\n";
        return 1;
    }

    auto tree = (TTree *)file.Get("o2sim");
    std::vector<o2::MCTrack> *tracks{};
    tree->SetBranchAddress("MCTrack", &tracks);

    int nLeptons{};
    int nAntileptons{};
    int nLeptonPairs{};
    int nLeptonPairsToBeDone{};
    int nSignalJpsi{};
    int nSignalPsi2S{};
    int nSignalJpsiWithinAcc{};
    int nSignalPsi2SWithinAcc{};
    auto nEvents = tree->GetEntries();
    o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
    Bool_t isInjected = kFALSE;

    for (int i = 0; i < nEvents; i++)
    {
        tree->GetEntry(i);
        for (auto &track : *tracks)
        {
            auto pdg = track.GetPdgCode();
            auto rapidity = track.GetRapidity();
            auto idMoth = track.getMotherTrackId();
            if (pdg == checkPdgDecay)
            {
                // count leptons
                nLeptons++;
            }
            else if (pdg == -checkPdgDecay)
            {
                // count anti-leptons
                nAntileptons++;
            }
            else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1])
            {
                if (idMoth < 0)
                {
                    // count signal PDG
                    pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++;
                    // count signal PDG within acceptance
                    if (rapidity > rapiditymin && rapidity < rapiditymax)
                    {
                        pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;
                    }
                }
                auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
                auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
                if (child0 != nullptr && child1 != nullptr)
                {
                    // check for parent-child relations
                    auto pdg0 = child0->GetPdgCode();
                    auto pdg1 = child1->GetPdgCode();
                    std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
                    if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1)
                    {
                        nLeptonPairs++;
                        if (child0->getToBeDone() && child1->getToBeDone())
                        {
                            nLeptonPairsToBeDone++;
                        }
                    }
                }
            }
        }
    }
    std::cout << "#events: " << nEvents << "\n"
              << "#leptons: " << nLeptons << "\n"
              << "#antileptons: " << nAntileptons << "\n"
              << "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n"
              << "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n"
              << "#lepton pairs: " << nLeptonPairs << "\n"
              << "#lepton pairs to be done: " << nLeptonPairs << "\n";

    if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0)
    {
        std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n";
        return 1;
    }
    if (nLeptonPairs != nLeptonPairsToBeDone)
    {
        std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n";
        return 1;
    }

    return 0;
}