1 | // A simple macro used to estimate the impact of possible EUSO rotations around X and Y axis within |
---|
2 | // the cone defined by Beta aperture angle on the direction, Hmax, and energy. |
---|
3 | // Author: Dmitry V.Naumov |
---|
4 | // Date: 15 November, 2003 |
---|
5 | // Use: root[].L TildedEUSO.C |
---|
6 | // root[] FillTree(Nevents) // where Nevents is the number of with different Theta, Phi, Beta |
---|
7 | // root[] ReadTree() |
---|
8 | void FillTree(int evMax) |
---|
9 | { |
---|
10 | double Theta, Phi,ErrorTheta,ErrorPhi,Alpha, Beta; |
---|
11 | double ThetaDeg, PhiDeg,ErrorThetaDeg,ErrorPhiDeg, AlphaDeg, BetaDeg, Dir, ErrorH,ErrorEnergy; |
---|
12 | double deg=TMath::RadToDeg(), BetaMax=5*TMath::DegToRad(); |
---|
13 | double AlphaX,AlphaY; |
---|
14 | int Events(0); |
---|
15 | TRotation r; |
---|
16 | TVector3 n(0,0,1),tilted,Omega,OmegaPrime; |
---|
17 | TMath math; |
---|
18 | // Creating the tree |
---|
19 | TTree tree("tree","Tildet EUSO"); |
---|
20 | TFile f("output/tilded_euso.root","recreate"); |
---|
21 | tree.Branch("Alpha",&AlphaDeg,"AlphaDeg/D"); |
---|
22 | tree.Branch("Beta",&BetaDeg,"BetaDeg/D"); |
---|
23 | tree.Branch("Theta",&ThetaDeg,"ThetaDeg/D"); |
---|
24 | tree.Branch("Phi",&PhiDeg,"PhiDeg/D"); |
---|
25 | tree.Branch("ErrorTheta",&ErrorThetaDeg,"ThetaDeg/D"); |
---|
26 | tree.Branch("ErrorPhi",&ErrorPhiDeg,"PhiDeg/D"); |
---|
27 | tree.Branch("Dir",&Dir,"Dir/D"); |
---|
28 | tree.Branch("ErrorH",&ErrorH,"ErrorH/D"); |
---|
29 | tree.Branch("ErrorEnergy",&ErrorEnergy,"ErrorEnergy/D"); |
---|
30 | |
---|
31 | while(Events<evMax){ |
---|
32 | Theta = math.PiOver2()*gRandom->Rndm(); |
---|
33 | Phi = math.TwoPi()*gRandom->Rndm(); |
---|
34 | Omega.SetXYZ(sin(Theta)*cos(Phi),sin(Theta)*sin(Phi),cos(Theta)); |
---|
35 | Beta = BetaMax*gRandom->Rndm(); |
---|
36 | Events++; |
---|
37 | for(int i=0;i<1000;i++) { |
---|
38 | Alpha = math.TwoPi()*gRandom->Rndm(); |
---|
39 | AlphaX = acos(sqrt(pow(sin(Beta),2)*pow(sin(Alpha),2)+pow(cos(Beta),2))); |
---|
40 | AlphaY = asin(sin(Beta)*cos(Alpha)); |
---|
41 | r.RotateY(AlphaY); |
---|
42 | r.RotateX(-AlphaX); |
---|
43 | OmegaPrime = r*Omega; |
---|
44 | Dir = Omega.Angle(OmegaPrime)*deg; |
---|
45 | ErrorTheta = OmegaPrime.Theta()-Theta; |
---|
46 | ErrorPhi = OmegaPrime.Phi()-Phi; |
---|
47 | AlphaDeg = Alpha*deg; |
---|
48 | BetaDeg = Beta*deg; |
---|
49 | ThetaDeg = Theta*deg; |
---|
50 | PhiDeg = Phi*deg; |
---|
51 | ErrorThetaDeg = ErrorTheta*deg; |
---|
52 | ErrorPhiDeg = ErrorPhi*deg; |
---|
53 | ErrorH = OmegaPrime.CosTheta()/cos(Theta)-1; |
---|
54 | ErrorEnergy = pow(cos(math.Pi()/6+Beta)/cos(TMath::Pi()/6),3)-1; // Maximum Error for 30 degrees |
---|
55 | tree.Fill(); |
---|
56 | r.SetToIdentity(); |
---|
57 | } |
---|
58 | } |
---|
59 | tree.Write(); |
---|
60 | f.Close(); |
---|
61 | } |
---|
62 | void ReadTree() |
---|
63 | { |
---|
64 | TFile *f = new TFile("output/tilded_euso.root",""); |
---|
65 | TTree *tree = (TTree*)f->Get("tree"); |
---|
66 | double Alpha, Beta, Theta, Phi,ErrorTheta,ErrorPhi,Dir, ErrorH,ErrorEnergy; |
---|
67 | tree->SetBranchAddress("Alpha",&Alpha); |
---|
68 | tree->SetBranchAddress("Beta",&Beta); |
---|
69 | tree->SetBranchAddress("Theta",&Theta); |
---|
70 | tree->SetBranchAddress("Phi",&Phi); |
---|
71 | tree->SetBranchAddress("ErrorTheta",&ErrorTheta); |
---|
72 | tree->SetBranchAddress("ErrorPhi",&ErrorPhi); |
---|
73 | tree->SetBranchAddress("Dir",&Dir); |
---|
74 | tree->SetBranchAddress("ErrorH",&ErrorH); |
---|
75 | tree->SetBranchAddress("ErrorEnergy",&ErrorEnergy); |
---|
76 | tree->StartViewer(); |
---|
77 | TH2F *t1 = new TH2F("t1","Altitude Relative Error",5,0,5,100,-1,1); |
---|
78 | TH2F *t2 = new TH2F("t2","Theta Error (deg)",5,0,5,100,-7,7); |
---|
79 | TH2F *t3 = new TH2F("t3","Phi Error (deg)",5,0,5,100,-20,20); |
---|
80 | TH2F *t4 = new TH2F("t4","Direction Error (deg)",5,0,5,100,0,10); |
---|
81 | tree->Draw("ErrorH:Beta>>t1","",""); |
---|
82 | t1->Draw(); |
---|
83 | c1->Print("output/ErrorH.gif"); |
---|
84 | t1->FitSlicesY(); |
---|
85 | t1_2->Draw(); |
---|
86 | c1->Print("output/ErrorHSigma.gif"); |
---|
87 | tree->Draw("ErrorEnergy:Beta"); c1->Print("output/ErrorEnergy.gif"); |
---|
88 | tree->Draw("ErrorTheta:Beta>>t2","",""); |
---|
89 | t2->Draw(); |
---|
90 | c1->Print("output/ErrorTheta.gif"); |
---|
91 | t2->FitSlicesY(); |
---|
92 | t2_2->Draw(); |
---|
93 | c1->Print("output/ErrorThetaSigma.gif"); |
---|
94 | tree->Draw("ErrorPhi:Beta>>t3","",""); |
---|
95 | t3->Draw(); |
---|
96 | c1->Print("output/ErrorPhi.gif"); |
---|
97 | t3->FitSlicesY(); |
---|
98 | t3_2->Draw(); |
---|
99 | c1->Print("output/ErrorPhiSigma.gif"); |
---|
100 | tree->Draw("Dir:Beta>>t4","",""); |
---|
101 | t4->Draw(); |
---|
102 | c1->Print("output/Dir.gif"); |
---|
103 | t4->FitSlicesY(); |
---|
104 | t4_2->Draw(); |
---|
105 | c1->Print("output/DirSigma.gif"); |
---|
106 | } |
---|
107 | |
---|
108 | void TildedEuso(int evMax) |
---|
109 | { |
---|
110 | FillTree(evMax); |
---|
111 | ReadTree(); |
---|
112 | } |
---|