K600-ALBA/plotK600.C at master · KevinCWLi/K600-ALBA · 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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
int main()
{

    gROOT->Reset();

    TFile f("Analyse.root","recreate");

    double nrofThreads = 24;

    //char rootFileDirectory[] = "K600-buildMT";
    char rootFileDirectory[] = "K600-buildMT2";
    //char rootFileDirectory[] = "K600-buildMT-1332";
    //char rootFileDirectory[] = "/Users/kevinchingweili/Physics/GEANT4/K600/master/K600-build";

    char name[512];
    char condition[512];
    char command[512];


    TChain *DataTreeChain = new TChain("DataTreeSim","DataTreeSim");
    TChain *GeometryAnalysisTree = new TChain("GeometryAnalysisTree","GeometryAnalysisTree");
    TChain *InputVariableTree = new TChain("InputVariableTree","InputVariableTree");

    for(Int_t i=0; i<nrofThreads; i++)
    {
        sprintf(name,"K600Output_t%d.root", i);
        //sprintf(name,"%s/K600Output.root", rootFileDirectory, i);

        DataTreeChain->Add(name);
        GeometryAnalysisTree->Add(name);
        InputVariableTree->Add(name);
    }


    TCanvas * c1 = new TCanvas("c1", "c1", 1000, 700);


    ////================================================////

    std::vector<std::shared_ptr<TH1F>> hEnergy_energy;
    std::vector<double> integrals;
    std::vector<double> efficiencies;

    for(int i=0; i<9; i++)
    {
        double initialParticleKineticEnergy = 0.0;

        if(i==0)
        {
            initialParticleKineticEnergy = 0.25;
        }
        else if(i==1)
        {
            initialParticleKineticEnergy = 0.5;
        }
        else if(i==2)
        {
            initialParticleKineticEnergy = 1;
        }
        else if(i==3)
        {
            initialParticleKineticEnergy = 1.332;
        }
        else if(i==4)
        {
            initialParticleKineticEnergy = 2;
        }
        else if(i==5)
        {
            initialParticleKineticEnergy = 5;
        }
        else if(i==6)
        {
            initialParticleKineticEnergy = 10;
        }
        else if(i==7)
        {
            initialParticleKineticEnergy = 15;
        }
        else if(i==8)
        {
            initialParticleKineticEnergy = 20;
        }

        //initialParticleKineticEnergy *= 1.0e+03;

        sprintf(name,"hEnergy_energy%d", i);
        hEnergy_energy.push_back(std::shared_ptr<TH1F>(new TH1F(name, "", 22000, 0.0, 22000.0)));

        double lowerE, higherE;
        lowerE = initialParticleKineticEnergy-0.010;
        higherE = initialParticleKineticEnergy+0.010;

        int lowerEBin, higherEBin;
        lowerEBin = (int) 1000*lowerE;
        higherEBin = (int) 1000*higherE;

        sprintf(condition,"(%f)<InitialParticleKineticEnergy && InitialParticleKineticEnergy<(%f)", lowerE, higherE);
        sprintf(command,"CLOVER_Energy>>%s", name);

        DataTreeChain->Draw(command, condition, "");
        //std::cout << "condition: " << condition << std::endl;
        //std::cout << "lowerEBin: " << lowerEBin << std::endl;
        //std::cout << "higherEBin: " << higherEBin << std::endl;

        DataTreeChain->Draw(command, condition, "");

        integrals.push_back(hEnergy_energy.back()->Integral(lowerEBin, higherEBin));
        efficiencies.push_back(100.0*integrals.back()/1000000.0);

        std::cout << "Energy: " << initialParticleKineticEnergy << " [MeV], Efficiency: " << efficiencies.back() << std::endl;
    }


    f.Write();
    c1->Close();

    return 0;

}