1 | // ESAF : Euso Simulation and Analysis Framework |
---|
2 | // $Id: PhotonsDatabaseBuilder.cc 2789 2008-02-11 07:26:05Z naumov $ |
---|
3 | // A.Thea, J.Watts created Mar, 19 2004 |
---|
4 | |
---|
5 | #include "PhotonsDatabaseBuilder.hh" |
---|
6 | #include "Photon.hh" |
---|
7 | #include "OpticalSystem.hh" |
---|
8 | #include "OpticsFactory.hh" |
---|
9 | #include "TTimeStamp.h" |
---|
10 | #include <TMath.h> |
---|
11 | |
---|
12 | ClassImp(PhotonsDatabaseBuilder) |
---|
13 | |
---|
14 | //______________________________________________________________________________ |
---|
15 | PhotonsDatabaseBuilder::PhotonsDatabaseBuilder() { |
---|
16 | // ctor |
---|
17 | |
---|
18 | fSamples = (Int_t)Conf()->GetNum("PhotonsDatabaseBuilder.fSamples"); |
---|
19 | fOptics = OpticsFactory::Get()->MakeOpticalSystem(); |
---|
20 | fRootFileName = Conf()->GetStr("PhotonsDatabaseBuilder.fRootFileName"); |
---|
21 | } |
---|
22 | |
---|
23 | //______________________________________________________________________________ |
---|
24 | PhotonsDatabaseBuilder::~PhotonsDatabaseBuilder() { |
---|
25 | // dtor |
---|
26 | |
---|
27 | delete fOptics; |
---|
28 | |
---|
29 | } |
---|
30 | |
---|
31 | //______________________________________________________________________________ |
---|
32 | void PhotonsDatabaseBuilder::CloseRoot() { |
---|
33 | // close root file |
---|
34 | if ( fRootFile ) { |
---|
35 | fRootFile->Write(); |
---|
36 | fRootFile->Close(); |
---|
37 | cout << "Root file closed\n"; |
---|
38 | delete fRootFile; fRootFile = 0; |
---|
39 | } |
---|
40 | } |
---|
41 | |
---|
42 | //______________________________________________________________________________ |
---|
43 | void PhotonsDatabaseBuilder::OpenRoot() { |
---|
44 | // open rootfile with the correct name |
---|
45 | CloseRoot(); |
---|
46 | |
---|
47 | cout << "Opening root file " << fRootFileName << endl; |
---|
48 | fRootFile = new TFile((fRootFileName+".root").c_str(),"CREATE"); |
---|
49 | if (fRootFile->IsZombie()) { |
---|
50 | cout << endl; |
---|
51 | cout << "Error opening file " << fRootFileName << endl; |
---|
52 | cout << "Probably it already exists or there is no " << endl; |
---|
53 | cout << "write access permission to this directory" << endl; |
---|
54 | cout << endl; |
---|
55 | |
---|
56 | // try adding date and time to the name |
---|
57 | TTimeStamp time; |
---|
58 | TString now = time.AsString("s"); |
---|
59 | now[now.Last(' ')]=':'; |
---|
60 | TString newFileName = Form("%s-", fRootFileName.c_str())+now+".root"; |
---|
61 | fRootFile = new TFile(newFileName,"CREATE"); |
---|
62 | if (fRootFile->IsZombie()) { |
---|
63 | throw invalid_argument("Cannot open root file"); |
---|
64 | } |
---|
65 | cout << "Root file " << newFileName << " successfully opened\n"; |
---|
66 | } |
---|
67 | |
---|
68 | } |
---|
69 | //______________________________________________________________________________ |
---|
70 | void PhotonsDatabaseBuilder::InitDBTree() { |
---|
71 | // create the tree |
---|
72 | |
---|
73 | fTree = new OADBTree; |
---|
74 | fPhotons = new OADBPhotons(); |
---|
75 | |
---|
76 | cout << fPupil.GetNumPos() << " spots to scan" << endl; |
---|
77 | fPhotons->SetSize(fPupil.GetNumPos()); |
---|
78 | fTree->Branch("photons","OADBPhotons",&fPhotons); |
---|
79 | |
---|
80 | // fill the header |
---|
81 | OADBHeader &header = fTree->GetHeader(); |
---|
82 | |
---|
83 | header.SetXYStep( fPupil.GetXYStep() ); |
---|
84 | header.SetOpticsRadius( fOptics->Radius() ); |
---|
85 | header.SetEntranceDiscRadius( fPupil.GetEntranceDiscRadius() ); |
---|
86 | |
---|
87 | header.SetFirstLensBottom( fOptics->FirstLensBottom() ); |
---|
88 | header.SetFirstLensTop( fOptics->FirstLensTop() ); |
---|
89 | header.SetSecondLensBottom( fOptics->SecondLensBottom() ); |
---|
90 | header.SetSecondLensTop( fOptics->SecondLensTop() ); |
---|
91 | header.SetEntranceZ(fOptics->Bottom()); |
---|
92 | header.SetExitZ(fOptics->Top()); |
---|
93 | header.SetPosZ(fOptics->Position().Z()); |
---|
94 | |
---|
95 | vector<Double_t> dummyD; |
---|
96 | dummyD = fPupil.GetAngles(); |
---|
97 | header.SetNumAngles(dummyD.size()); |
---|
98 | for(size_t i(0); i<dummyD.size(); i++ ) |
---|
99 | header.SetAngle(i, dummyD[i]); |
---|
100 | cout << header.GetNumAngles() << " angles saved in the header" << endl; |
---|
101 | |
---|
102 | dummyD = fPupil.GetWavelengths(); |
---|
103 | header.SetNumWavelengths(dummyD.size()); |
---|
104 | for(size_t i(0); i<dummyD.size(); i++ ) |
---|
105 | header.SetWavelength(i, dummyD[i]); |
---|
106 | |
---|
107 | vector<Int_t> dummyI = fPupil.GetRows(); |
---|
108 | header.SetNumRows(dummyI.size()); |
---|
109 | for(size_t i(0); i<dummyI.size(); i++ ) |
---|
110 | header.SetNumCols(i, dummyI[i]); |
---|
111 | |
---|
112 | } |
---|
113 | |
---|
114 | //______________________________________________________________________________ |
---|
115 | void PhotonsDatabaseBuilder::FillDBTree() { |
---|
116 | // calls fTree->Fill() and then clears the photon container |
---|
117 | |
---|
118 | fTree->Fill(); |
---|
119 | fPhotons->Clean(); |
---|
120 | |
---|
121 | // check rootfile |
---|
122 | if ( fRootFile != fTree->GetCurrentFile() ) { |
---|
123 | fRootFile = fTree->GetCurrentFile(); |
---|
124 | cout << Form("RootFile capacity limit reached. Moving to %s", |
---|
125 | fRootFile->GetName()) << endl; |
---|
126 | } |
---|
127 | |
---|
128 | } |
---|
129 | |
---|
130 | //______________________________________________________________________________ |
---|
131 | void PhotonsDatabaseBuilder::CloseDBTree() { |
---|
132 | // write and close the tree |
---|
133 | |
---|
134 | cout << "Saving rootfile" << endl; |
---|
135 | fTree->Write(); |
---|
136 | |
---|
137 | delete fTree; fTree = 0; |
---|
138 | delete fPhotons; fPhotons = 0; |
---|
139 | |
---|
140 | } |
---|
141 | |
---|
142 | //______________________________________________________________________________ |
---|
143 | Bool_t PhotonsDatabaseBuilder::Sample( Photon *ph ) { |
---|
144 | // resample the same photon and save the most rated path in fLastPath |
---|
145 | |
---|
146 | vector<pair<Photon,Int_t> > buffer; |
---|
147 | |
---|
148 | // fire fSamples Photons |
---|
149 | for(Int_t k(0); k< fSamples; k++){ |
---|
150 | // reset the photon |
---|
151 | ph->parent->CreatePhoton(); |
---|
152 | |
---|
153 | // transport the ph |
---|
154 | Photon* newph = fOptics->Transport(ph); |
---|
155 | // if survives and it is out of the second lens add to the list |
---|
156 | // FIXME it does not take in account the phs that go out from the sides |
---|
157 | if ( newph && TMath::Abs(newph->pos[Z]-fOptics->SecondLensTop()) < TOLERANCE) { |
---|
158 | vector< pair<Photon, Int_t> >::iterator it; |
---|
159 | Bool_t isNew(kTRUE); |
---|
160 | |
---|
161 | // look for a similar ph in the buffer |
---|
162 | for(it=buffer.begin(); it!=buffer.end(); it++){ |
---|
163 | Float_t dirDiff = (it->first.dir.Unit()-newph->dir.Unit()).Mag(); |
---|
164 | Float_t posDiff = (it->first.pos-newph->pos).Mag(); |
---|
165 | if (dirDiff < TOLERANCE && posDiff < TOLERANCE ){ |
---|
166 | isNew = kFALSE; |
---|
167 | it->second++; |
---|
168 | break; |
---|
169 | } |
---|
170 | } |
---|
171 | // this ph is new, add an entry |
---|
172 | if (isNew) { |
---|
173 | pair<Photon,Int_t> dummy(*newph, 1); |
---|
174 | buffer.push_back(dummy); |
---|
175 | } |
---|
176 | } |
---|
177 | } |
---|
178 | vector< pair<Photon, Int_t> >::iterator it; |
---|
179 | fLastPath.second = 0; |
---|
180 | |
---|
181 | for (it = buffer.begin(); it != buffer.end(); it++) { |
---|
182 | if (it->second > fLastPath.second) { |
---|
183 | // calculate transmission and put it in the path |
---|
184 | fLastPath.first = it->first; |
---|
185 | fLastPath.second = it->second/(Float_t)fSamples; |
---|
186 | } |
---|
187 | } |
---|
188 | |
---|
189 | if ( fLastPath.second == 0 ) |
---|
190 | return kFALSE; |
---|
191 | else |
---|
192 | return kTRUE; |
---|
193 | |
---|
194 | } |
---|
195 | |
---|
196 | //______________________________________________________________________________ |
---|
197 | void PhotonsDatabaseBuilder::Go() { |
---|
198 | // build the database |
---|
199 | |
---|
200 | OpenRoot(); |
---|
201 | InitDBTree(); |
---|
202 | for(; !fPupil.IsWlEnd() ;fPupil.NextWl() ) { |
---|
203 | cout << "Switching to wavelength " << fPupil.GetCurWl() << endl; |
---|
204 | for(; !fPupil.IsAngleEnd() ;fPupil.NextAngle() ) { |
---|
205 | cout << "Switching to angle " << TMath::RadToDeg()*fPupil.GetCurAngle() << endl; |
---|
206 | fPhotons->SetAngle(fPupil.GetCurAngle(), fPupil.GetCurAngleId()); |
---|
207 | fPhotons->SetWavelength(fPupil.GetCurWl(), fPupil.GetCurWlId()); |
---|
208 | while ( Photon *ph = fPupil.GetPhoton() ) { |
---|
209 | // find the "top rated photon" |
---|
210 | if ( Sample( ph ) ) { |
---|
211 | // get the "top rated photon" and save it |
---|
212 | fPhotons->SavePhoton(fPupil.GetLastPosId(), &(fLastPath.first), fLastPath.second ); |
---|
213 | } else { |
---|
214 | //if there are no ph in the database, save a ghost |
---|
215 | fPhotons->SaveGhost(fPupil.GetLastPosId()); |
---|
216 | } |
---|
217 | } |
---|
218 | |
---|
219 | FillDBTree(); |
---|
220 | } |
---|
221 | } |
---|
222 | |
---|
223 | |
---|
224 | CloseDBTree(); |
---|
225 | CloseRoot(); |
---|
226 | } |
---|