test2_p.C

макрос для получения импульса - Bahytbek Mauyey, 11/10/2017 03:36 PM

Download (3.98 KB)

 
1
// Global declarations
2
TGraph* gr, igr;
3

    
4
// C++ defined function
5
Double_t f(Double_t *x, Double_t *par)
6
{
7
        return igr->Eval(x[0]);
8
}
9

    
10
ofstream out("CMom_B11_N15.txt");
11

    
12
void test2_p()
13
{
14
        // Read the text file into the tree
15
        // TFile file("myfile.root","RECREATE");
16
        TTree* MyTree = new TTree("MyTree", "MyTree");
17
        Double_t tetta, crossSection;
18
        Double_t CMmom, shorty,phi, px,py,pz,pz_B,px_B,py_B;
19
        MyTree->ReadFile("cos_tetta_cross.txt", "tetta/D:crossSection/D");
20
 TBranch *bpt = MyTree->Branch("px",&px,"px/F");
21
        
22
        MyTree->SetBranchAddress("tetta", &tetta);
23
        MyTree->SetBranchAddress("crossSection", &crossSection);
24
        /*
25
        MyTree->SetBranchAddress("pz", &pz);
26
        MyTree->SetBranchAddress("px", &px);
27
        MyTree->SetBranchAddress("py", &py);
28
        MyTree->SetBranchAddress("pz_B", &pz_B);
29
        MyTree->SetBranchAddress("px_B", &px_B);
30
        MyTree->SetBranchAddress("py_B", &py_B);
31
        */
32
        // Create two vectors to be filled from the tree
33
        const Int_t nEntries = MyTree->GetEntries(); // number of entries in the input file
34
        TVectorD tet(nEntries);
35
        TVectorD sigma(nEntries);
36

    
37
        // Fill two vectors from the tree
38
        for (int i=0; i<nEntries; i++) {
39
                MyTree->GetEntry(i);
40
                tet(i) = tetta;
41
                sigma(i) = crossSection;
42
        }
43

    
44
        // Create the graph
45
        gr = new TGraph(tet, sigma);
46
        igr = new TGraph(sigma,tet);
47
        // Draw the graph
48
        TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);
49
        c1->cd();
50
        igr->Draw();
51
        igr->SetLineColor(kRed);
52
        igr->SetLineWidth(3);
53
        gPad->SetGrid(1, 1);
54

    
55
        // ROOT-TF1 object with the reference to the C++ defined function
56
        TF1* func = new TF1("func", f, 0, 1, 0);
57

    
58
        // Draw the func
59
        TCanvas *c2 = new TCanvas("c2", "c2", 800, 600);
60
        c2->cd();
61
        func->Draw();
62
        func->SetLineColor(kGreen);
63
        func->SetLineWidth(3);
64
        gPad->SetGrid(1, 1);
65

    
66
        // Uniform random number generator
67
        TRandom3* rndGen = new TRandom3();
68

    
69
        TH1D* histo1 = new TH1D("histo1", "histo1", 1800, 0., 1.);
70
        TH1D* histo2 = new TH1D("histo2", "histo2", 1800, 0., 180.);
71
        TH1D* histo3 = new TH1D("histo3", "histo3", 1800, 0., 180.);
72
        TH1D* histo4 = new TH1D("histo4", "histo4", 1800, 0., 180.);
73
        TH1D* histo5 = new TH1D("histo5", "histo5", 1800, 0., 180.);
74

    
75
        UInt_t nIterations = 500000;
76
        // Loop over all iterations
77
        for (UInt_t iIter=0; iIter<nIterations; iIter++)
78
        {
79
                if (iIter%100000 == 0) cout << "Iteration " << iIter << endl;
80

    
81
                Double_t x = rndGen->Rndm();        // from 0.0 to 1.0
82
                Double_t y = func->Eval(x);
83
   //___________fill px py pz -----------------
84
        
85
Double_t         invariant;
86
                Double_t AMU=931.5 ;
87
                Double_t m11B=AMU*11.;
88
                Double_t m15N=AMU*15;
89
                Double_t T_N15=43;
90
Double_t CMmom, shorty,phi, px,py,pz;
91
                invariant = (m11B+m15N)*(m11B+m15N)+2*m11B*T_N15;
92

    
93
Double_t shorty=(invariant-m11B*m11B-m15N*m15N)*(invariant-m11B*m11B-m15N*m15N);
94
CMmom = sqrt((shorty-4*m15N*m15N*m11B*m11B)/(4*invariant));
95
                //phi= CLHEP::RandFlat::shoot(0.0,2*double(CLHEP::pi));
96
                Double_t PI = 3.14159265358;
97
                phi = rand()% 2*PI;
98
                pz=CMmom*cos(y);
99
                px=CMmom*sin(y)*cos(phi*180/PI) ;//cos(phi);
100
                py=CMmom*sin(y)*sin(phi*180/PI); //sin(phi);
101
                pz_B=-pz;
102
                px_B=-px;
103
                py_B=-py;
104
                
105
                out<<pz<<"\t"<<px<<"\t"<< py<<"\t"<<-pz<<"\t"<<-px<<"\t"<< -py<<endl;
106
                
107
                // Histograms filling
108
                histo1->Fill(x);
109
                histo2->Fill(y);
110
                histo3->Fill(pz);
111
                histo4->Fill(px);
112
                histo5->Fill(py);
113
                //file.Write();
114
                
115
        }
116

    
117
        // Draw
118
        TCanvas* canv1 = new TCanvas("canv1", "canv1");
119
        canv1->cd();
120
        histo1->Draw();
121
        gStyle->SetOptStat(111111);
122
        gPad->SetGrid(1, 1);
123

    
124
        // Draw
125
        TCanvas* canv2 = new TCanvas("canv2", "canv2");
126
        canv2->cd();
127
        histo2->Draw();
128
        gStyle->SetOptStat(111111);
129
        gPad->SetGrid(1, 1);
130
        
131
                TCanvas* canv3 = new TCanvas("canv3", "canv3");
132
        canv3->cd();
133
        histo3->Draw();
134
        gStyle->SetOptStat(111111);
135
        gPad->SetGrid(1, 1);
136
        
137
                TCanvas* canv4 = new TCanvas("canv4", "canv4");
138
        canv3->cd();
139
        histo4->Draw();
140
        gStyle->SetOptStat(111111);
141
        gPad->SetGrid(1, 1);
142
        
143
                TCanvas* canv5 = new TCanvas("canv5", "canv5");
144
        canv3->cd();
145
        histo5->Draw();
146
        gStyle->SetOptStat(111111);
147
        gPad->SetGrid(1, 1);
148
        
149
}