source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/macros/App/TildedEuso.C @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 4.0 KB
Line 
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()
8void 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}
62void 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
108void TildedEuso(int evMax)
109{
110  FillTree(evMax);
111  ReadTree();
112}
Note: See TracBrowser for help on using the repository browser.