-
Notifications
You must be signed in to change notification settings - Fork 0
/
example_dipoleHiRes.cc
172 lines (137 loc) · 5.67 KB
/
example_dipoleHiRes.cc
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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
#include <iostream>
#include <cmath>
#include "dipoleHiRes.h"
#include "simuevents.h"
#include "common.h"
#include "maptools.h"
#include "projmap.h"
#include "TROOT.h"
#include "TRint.h"
#include "TCanvas.h"
#include "TStyle.h"
#ifdef gcc323
char* operator+( std::streampos&, char* );
#endif
using namespace std;
void Usage(string myName)
{
cout << endl;
cout << " Synopsis : " << endl;
cout << myName << " <events file> <simu flag>" << endl << endl;
cout << " Description :" << endl;
cout << myName << " is based upon comparison between data and a large quantity of events generated "
<< "isotropically by Monte Carlo simulation to determine the presence of a dipole. A binning technisque "
<< "that considers the event counts for the full range of opening angles from the center of each proposed "
<< "dipole distribution is used. The <simu flag> must be set to 0 if one wants to use the events stored "
<< "in <events file> and set to 1 if one wants to use simulated events according to a dipole pattern. In "
<< "the latter case <events file> will be ignored. The <events file> must contain the following fields :"
<< endl;
DumpFields();
cout << endl;
exit(0);
}
int main(int argc, char* argv[])
{
////////////////////////////////////////////////////////////////////////////
// //
// To start (initialization) //
// //
////////////////////////////////////////////////////////////////////////////
// Command line
if(argc != 3) Usage(argv[0]);
string eventFile = argv[1];
int isSimu = atoi(argv[2]);
if( !CheckFile(eventFile) ) {cerr << "File: " << eventFile << " not found" << endl; exit(0);}
// ROOT
int fargc = 1;
string extension;
TRint* rint = new TRint("Fit Dipole HiRes", &fargc, argv);
extension = ".png";
gROOT->SetStyle("Plain");
gStyle->SetTitleFont(30,"TITLE");
////////////////////////////////////////////////////////////////////////////
// //
// dipole (simulation) & events (file or simulation) //
// //
////////////////////////////////////////////////////////////////////////////
unsigned int nSide = 128;
double latSite = kConstantsTK::AugerSouthLatitude;
double lonSite = kConstantsTK::AugerSouthLongitude;
double thetaMax = 60.;
unsigned int nVal = 10000;
vector<double> thVal(nVal);
vector<double> pthVal(nVal);
for(unsigned int i=0; i<nVal; i++)
{
thVal[i] = i*180./(nVal-1);
pthVal[i] = sin(thVal[i]*M_PI/180)*cos(thVal[i]*M_PI/180);
if (thVal[i] > thetaMax) pthVal[i] = 0;
}
double lDipole, bDipole, ampDipole;
ampDipole = 0.05;
lDipole = 0.;
bDipole = 45.;
vector<TEvent> events;
if(isSimu == 1)
{
// Simulation of my dipole
unsigned int nSide = 128;
THealpixMap mapSimu(nSide);
vector<long> iPix(mapSimu.NPix());
vector<double> lPix; // Deg
vector<double> bPix; // Deg
for(unsigned int i = 0; i < mapSimu.NPix(); i++) iPix[i] = i;
mapSimu.GiveLB(iPix,lPix,bPix);
vector<double> uvDipole = ll2uv(lDipole, bDipole);
vector<double> dipoleMap(mapSimu.NPix());
for(unsigned int i=0; i<mapSimu.NPix(); i++)
{
vector<double> uvPix = ll2uv(lPix[i], bPix[i]);
dipoleMap[i] = ampDipole*(uvPix[0]*uvDipole[0]+uvPix[1]*uvDipole[1]+uvPix[2]*uvDipole[2]);
}
mapSimu = 1.+dipoleMap;
unsigned int nb = 50000;
events = SimulateEvents(mapSimu, nb, thVal, pthVal, latSite, lonSite);
cout << "Number of simulated events: " << events.size() << endl;
}
else
{
cout << "Reading events file " << eventFile << endl;
events = GetEvents(eventFile);
long nEvents = events.size();
if(nEvents == 0) {cout << "Program Failed: No events read. Exiting." << endl; exit(0);}
}
////////////////////////////////////////////////////////////////////////////
// //
// Analysis //
// //
////////////////////////////////////////////////////////////////////////////
double dCosTheta = 0.05;
TDipoleHiRes dipole(events, lDipole, bDipole, dCosTheta);
// Events Opening Angle
dipole.ComputeEventsOpeningAngle();
// Background Opening Angle
unsigned int nSimu = 100;
THealpixMap map(nSide);
map = map+1;
dipole.ComputeCoverageOpeningAngle(lonSite, latSite, map, nSimu, thVal, pthVal);
dipole.ComputeTrueOpeningAngle();
dipole.fitDipole();
dipole.DrawEvents();
dipole.DrawCov();
dipole.DrawBoth();
dipole.DrawTrueFit();
vector<double> fitParameters = dipole.GetFitParameters();
double yIntercept = fitParameters[0];
double slope = fitParameters[1];
vector<double> fitParametersErrors = dipole.GetFitParametersErrors();
double yInterceptError = fitParametersErrors[0];
double slopeError = fitParametersErrors[1];
double amplitude, varAmplitude, dAmplitude;
amplitude = slope/yIntercept;
varAmplitude = (slope*yInterceptError/pow(yIntercept,2))*(slope*yInterceptError/pow(yIntercept,2))+(slopeError/yIntercept)*(slopeError/yIntercept);
dAmplitude = sqrt(varAmplitude);
cout << "amplitude = " << amplitude << " +/- " << dAmplitude << endl;
cout << "Program Finished Normally" << endl;
rint->Run(kTRUE);
}