@@ -35,17 +35,30 @@ ElectronScattering_DurhamData::ElectronScattering_DurhamData(nuiskey samplekey)
3535 fSettings .SetAllowedTypes (" FIX,FREE,SHAPE/DIAG/NORM/MASK" " FIX/DIAG"
3636 fSettings .SetXTitle (" q0"
3737 fSettings .SetYTitle (" #sigma"
38+ // fIsRawEvents = true;
3839
3940 FinaliseSampleSettings ();
4041
42+ // Plot Setup -------------------------------------------------------
43+ SetDataFromName (fSettings .GetS (" name"
44+ SetCovarFromDiagonal ();
45+
4146 // Scaling Setup ---------------------------------------------------
4247 // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
4348 // fScaleFactor = ((GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / TotalIntegratedFlux());
44- fScaleFactor = 1.0 ;
49+ EnuMin = fZLowLim ;
50+ EnuMax = fZHighLim ;
4551
46- // Plot Setup -------------------------------------------------------
47- SetDataFromName (fSettings .GetS (" name"
48- SetCovarFromDiagonal ();
52+ double sigscale = GetEventHistogram ()->Integral () * 1E-38 / double (fNEvents ) / TotalIntegratedFlux ();
53+ double dangle = 2 * M_PI * fabs ((1 . - cos (fYLowLim * M_PI / 180 .)) - (1 . - cos (fYHighLim * M_PI / 180 .)));
54+ fScaleFactor = sigscale / dangle / fZCenter ;
55+
56+
57+ std::cout << " Event Integral = " GetEventHistogram ()->Integral (" width"
58+ std::cout << " Flux Integral = " TotalIntegratedFlux (" width" " " fZLowLim << " " fZHighLim << std::endl;
59+ std::cout << " DAngle = "
60+ std::cout << " sigscale = "
61+ std::cout << " fZCenter = " fZCenter << std::endl;
4962
5063 // Finish up
5164 FinaliseMeasurement ();
@@ -80,7 +93,13 @@ void ElectronScattering_DurhamData::SetDataFromName(std::string name) {
8093 thetabinedges.push_back ( fYCenter + thetawidth * (double (i)) );
8194 }
8295 for (int i = -nebins; i <= nebins; i++) {
83- ebinedges.push_back ( fZCenter + ewidth * (double (i)) );
96+ double newval = fZCenter + ewidth * (double (i));
97+
98+ if (newval < 0.0 ) newval = 0.0 ;
99+ if (newval < GetEventHistogram ()->GetXaxis ()->GetXmin ()) newval = GetEventHistogram ()->GetXaxis ()->GetXmin ();
100+ if (newval > GetEventHistogram ()->GetXaxis ()->GetXmax ()) newval = GetEventHistogram ()->GetXaxis ()->GetXmax ();
101+ if (std::find (ebinedges.begin (), ebinedges.end (), newval) != ebinedges.end ()) continue ;
102+ ebinedges.push_back (newval);
84103 }
85104
86105 // Determine target
@@ -102,6 +121,7 @@ void ElectronScattering_DurhamData::SetDataFromName(std::string name) {
102121 std::vector<double > errorx;
103122 std::vector<double > pointy;
104123 std::vector<double > errory;
124+ double scalef = 1 .E -38 * 1 .E5 ;
105125
106126 while (std::getline (mask >> std::ws, line, ' \n '
107127 // std::cout << "Line = " << line << std::endl;
@@ -116,23 +136,71 @@ void ElectronScattering_DurhamData::SetDataFromName(std::string name) {
116136 if (sstring.compare (lineentries[7 ])) continue ;
117137
118138 std::cout << " Registering data point : "
119-
120139 std::cout << " Adding Graph Point : " GeneralUtils::StrToDbl (lineentries[4 ]) << " " GeneralUtils::StrToDbl (lineentries[5 ]) << std::endl;
121- pointx.push_back (GeneralUtils::StrToDbl (lineentries[4 ]));
122- errorx.push_back (0.0 );
123- pointy.push_back (GeneralUtils::StrToDbl (lineentries[5 ]));
124- errory.push_back (GeneralUtils::StrToDbl (lineentries[6 ]));
140+
141+ // Loop through x and y points and find a place to insert
142+ if (pointx.empty ()) {
143+ pointx.push_back (GeneralUtils::StrToDbl (lineentries[4 ]));
144+ errorx.push_back (0.0 );
145+ pointy.push_back (GeneralUtils::StrToDbl (lineentries[5 ]) * scalef);
146+ errory.push_back (GeneralUtils::StrToDbl (lineentries[6 ]) * scalef);
147+
148+ } else {
149+ for (size_t j = 0 ; j < pointx.size (); j++) {
150+
151+ if (GeneralUtils::StrToDbl (lineentries[4 ]) < pointx[j] && j == 0 ) {
152+ std::cout << " Inserting at start point iterator "
153+ pointx.insert (pointx.begin () + j, GeneralUtils::StrToDbl (lineentries[4 ]));
154+ errorx.insert (errorx.begin () + j, 0.0 );
155+ pointy.insert (pointy.begin () + j, GeneralUtils::StrToDbl (lineentries[5 ]) * scalef);
156+ errory.insert (errory.begin () + j, GeneralUtils::StrToDbl (lineentries[6 ]) * scalef);
157+ break ;
158+
159+ } else if (GeneralUtils::StrToDbl (lineentries[4 ]) > pointx[j] && j == pointx.size () - 1 ) {
160+ std::cout << " Pushing back data point "
161+ pointx.push_back (GeneralUtils::StrToDbl (lineentries[4 ]));
162+ errorx.push_back (0.0 );
163+ pointy.push_back (GeneralUtils::StrToDbl (lineentries[5 ]) * scalef);
164+ errory.push_back (GeneralUtils::StrToDbl (lineentries[6 ]) * scalef);
165+ break ;
166+
167+ } else if (GeneralUtils::StrToDbl (lineentries[4 ]) > pointx[j - 1 ] && GeneralUtils::StrToDbl (lineentries[4 ]) < pointx[j]) {
168+ std::cout << " Inserting at point iterator = "
169+
170+ pointx.insert (pointx.begin () + j, GeneralUtils::StrToDbl (lineentries[4 ]));
171+ errorx.insert (errorx.begin () + j, 0.0 );
172+ pointy.insert (pointy.begin () + j, GeneralUtils::StrToDbl (lineentries[5 ]) * scalef);
173+ errory.insert (errory.begin () + j, GeneralUtils::StrToDbl (lineentries[6 ]) * scalef);
174+ break ;
175+ }
176+ }
177+
178+
179+
180+ }
181+
182+ // pointx.push_back(GeneralUtils::StrToDbl(lineentries[4]));
183+ // errorx.push_back(0.0);
184+ // pointy.push_back(GeneralUtils::StrToDbl(lineentries[5]));
185+ // errory.push_back(GeneralUtils::StrToDbl(lineentries[6]));
125186
126187 i++;
127188 }
128189
190+ for (int i = 0 ; i < pointx.size (); i++) {
191+ std::cout << " Q0 Point " " = "
192+ }
193+
194+
129195 fDataGraph = new TGraphErrors (pointx.size (), &pointx[0 ], &pointy[0 ], &errorx[0 ], &errory[0 ]);
196+ fDataGraph ->SetNameTitle ((fName + " _data_GRAPH" c_str (), (fName + " _data_GRAPH" c_str ());
197+
130198
131- // Now form an effective data and mc histogram
199+ // Now form an effective data and mc histogram
132200 std::vector<double > q0binedges;
133201 const double * x = fDataGraph ->GetX ();
134202
135- // Loop over graph and get mid way point between each data point.
203+ // Loop over graph and get mid way point between each data point.
136204 for (int i = 0 ; i < fDataGraph ->GetN (); i++) {
137205 std::cout << " X Point = "
138206
@@ -147,9 +215,13 @@ void ElectronScattering_DurhamData::SetDataFromName(std::string name) {
147215 // Set half distance to point below
148216 q0binedges.push_back (x[i] - ((x[i] - x[i - 1 ]) / 2.0 ));
149217 }
150-
151218 }
152219
220+ // Bubble Sort
221+
222+
223+
224+
153225 for (int i = 0 ; i < q0binedges.size (); i++) {
154226 std::cout << " Q0 Edge " " = "
155227 }
@@ -161,8 +233,8 @@ void ElectronScattering_DurhamData::SetDataFromName(std::string name) {
161233 std::cout << " theta Edge " " = "
162234 }
163235
164- // Form the data hist, mchist, etc
165- fDataHist = new TH1D (" electron_data " , " electron_data "
236+ // Form the data hist, mchist, etc
237+ fDataHist = new TH1D (( fName + " _data " ). c_str (), ( fName + " _data " ). c_str () ,
166238 q0binedges.size () - 1 , &q0binedges[0 ]);
167239 fMCHist = (TH1D*) fDataHist ->Clone (" MC"
168240
@@ -179,25 +251,34 @@ void ElectronScattering_DurhamData::SetDataFromName(std::string name) {
179251
180252
181253
182- fMCScan_Q0vsThetavsE = new TH3D (" mc_q0vsthetavse " " mc_q0vsthetavse "
254+ fMCScan_Q0vsThetavsE = new TH3D (( fName + " _MC_q0vsthetavse " ). c_str () , " MC_q0vsthetavse "
183255 q0binedges.size () - 1 , &q0binedges[0 ],
184256 thetabinedges.size () - 1 , &thetabinedges[0 ],
185257 ebinedges.size () - 1 , &ebinedges[0 ]);
186258 fMCScan_Q0vsThetavsE ->Reset ();
187- fMCScan_Q0vsTheta = new TH2D (" mc_q0vstheta " " mc_q0vstheta "
259+ fMCScan_Q0vsTheta = new TH2D (( fName + " _MC_q0vstheta " ). c_str () , " MC_q0vstheta "
188260 q0binedges.size () - 1 , &q0binedges[0 ],
189261 thetabinedges.size () - 1 , &thetabinedges[0 ]);
190262 fMCScan_Q0vsTheta ->Reset ();
191263
192- fMCScan_Q0vsE = new TH2D (" mc_q0vse " " mc_q0vse "
264+ fMCScan_Q0vsE = new TH2D (( fName + " _MC_q0vse " ). c_str () , " MC_q0vse "
193265 q0binedges.size () - 1 , &q0binedges[0 ],
194266 ebinedges.size () - 1 , &ebinedges[0 ]);
195267 fMCScan_Q0vsE ->Reset ();
268+
269+ fXLowLim = fMCScan_Q0vsThetavsE ->GetXaxis ()->GetBinLowEdge (1 );
270+ fXHighLim = fMCScan_Q0vsThetavsE ->GetXaxis ()->GetBinLowEdge (fMCScan_Q0vsThetavsE ->GetNbinsX () + 2 );
271+
272+ fYLowLim = fMCScan_Q0vsThetavsE ->GetYaxis ()->GetBinLowEdge (1 );
273+ fYHighLim = fMCScan_Q0vsThetavsE ->GetYaxis ()->GetBinLowEdge (fMCScan_Q0vsThetavsE ->GetNbinsY () + 2 );
274+
275+ fZLowLim = fMCScan_Q0vsThetavsE ->GetZaxis ()->GetBinLowEdge (1 );
276+ fZHighLim = fMCScan_Q0vsThetavsE ->GetZaxis ()->GetBinLowEdge (fMCScan_Q0vsThetavsE ->GetNbinsZ () + 2 );
196277}
197278
198279
199280// ********************************************************************
200- void ElectronScattering_DurhamData::FillEventVariables (FitEvent *event) {
281+ void ElectronScattering_DurhamData::FillEventVariables (FitEvent * event) {
201282// ********************************************************************
202283
203284 if (event->NumFSParticle (11 ) == 0 )
@@ -208,28 +289,31 @@ void ElectronScattering_DurhamData::FillEventVariables(FitEvent *event) {
208289
209290 double q0 = fabs (ein->fP .E () - eout->fP .E ()) / 1000.0 ;
210291 double E = ein->fP .E () / 1000.0 ;
211- double theta = ein->fP .Vect ().Angle (eout->fP .Vect ()) * 180 ;
292+ double theta = ein->fP .Vect ().Angle (eout->fP .Vect ()) * 180 . / M_PI ;
212293
213294 fXVar = q0;
214295 fYVar = theta;
215296 fZVar = E;
216297
217- std::cout << " Got Event " " " " "
298+ // std::cout << "Got Event " << q0 << " " << theta << " " << E << std::endl;
218299
219300 return ;
220301};
221302
222303// ********************************************************************
223- bool ElectronScattering_DurhamData::isSignal (FitEvent *event) {
304+ bool ElectronScattering_DurhamData::isSignal (FitEvent * event) {
224305// ********************************************************************
225306
226307 if (event->NumFSParticle (11 ) == 0 )
227308 return false ;
228309
310+ // std::cout << "fXVar = " << fXVar << " " << fXLowLim << " " << fXHighLim << std::endl;
311+ // std::cout << "fYVar = " << fYVar << " " << fYLowLim << " " << fYHighLim << std::endl;
312+ // std::cout << "fZVar = " << fZVar << " " << fZLowLim << " " << fZHighLim << std::endl;
229313
230- // if (fXVar < fXLowLim or fXVar > fXHighLim) return false;
231- // if (fYVar < fYLowLim or fYVar > fYHighLim) return false;
232- // if (fZVar < fZLowLim or fZVar > fZHighLim) return false;
314+ if (fXVar < fXLowLim or fXVar > fXHighLim ) return false ;
315+ if (fYVar < fYLowLim or fYVar > fYHighLim ) return false ;
316+ if (fZVar < fZLowLim or fZVar > fZHighLim ) return false ;
233317
234318 return true ;
235319};
@@ -238,30 +322,44 @@ bool ElectronScattering_DurhamData::isSignal(FitEvent *event) {
238322void ElectronScattering_DurhamData::FillHistograms () {
239323// ********************************************************************
240324
325+ Measurement1D::FillHistograms ();
326+
241327 if (Signal) {
242328 fMCScan_Q0vsThetavsE ->Fill (fXVar , fYVar , fZVar );
243-
244- fMCHist ->Fill (fXVar );
245- fMCScan_Q0vsTheta ->Fill (fXVar , fYVar );
246- fMCScan_Q0vsE ->Fill (fXVar , fZVar );
247- }
329+ fMCScan_Q0vsTheta ->Fill (fXVar , fYVar );
330+ fMCScan_Q0vsE ->Fill (fXVar , fZVar );
331+ }
248332}
333+ // Weight = 1.0;
334+ // if (Signal) {
335+ // fMCHist->Fill(fXVar, Weight);
336+ // fMCFine->Fill(fXVar, Weight);
337+ // fMCStat->Fill(fXVar, 1.0);
338+
339+ // if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight);
340+ // }
341+ // }
249342
250343void ElectronScattering_DurhamData::ResetAll () {
344+ Measurement1D::ResetAll ();
251345 fMCScan_Q0vsThetavsE ->Reset ();
346+ fMCScan_Q0vsTheta ->Reset ();
347+ fMCScan_Q0vsE ->Reset ();
252348}
253349
254350void ElectronScattering_DurhamData::ApplyNormScale (double norm) {
351+ Measurement1D::ApplyNormScale (norm);
255352 fMCScan_Q0vsThetavsE ->Scale (1.0 / norm);
256353 fMCScan_Q0vsTheta ->Scale (1.0 / norm);
257354 fMCScan_Q0vsE ->Scale (1.0 / norm);
258- fMCHist ->Scale (1.0 / norm);
259355}
260356
261357// ********************************************************************
262358void ElectronScattering_DurhamData::ScaleEvents () {
263359// ********************************************************************
360+ Measurement1D::ScaleEvents ();
264361
362+ /*
265363 fMCScan_Q0vsThetavsE->Scale(fScaleFactor, "width");
266364
267365 // Project into fMCScan_Q0vsTheta
@@ -300,7 +398,7 @@ void ElectronScattering_DurhamData::ScaleEvents() {
300398 }
301399
302400 fMCHist->Scale(fDataHist->Integral() / fMCHist->Integral());
303-
401+ */
304402}
305403
306404// ********************************************************************
@@ -312,13 +410,13 @@ int ElectronScattering_DurhamData::GetNDOF() {
312410
313411void ElectronScattering_DurhamData::Write (std::string drawOpts) {
314412
413+ Measurement1D::Write (drawOpts);
414+
315415 fMCScan_Q0vsThetavsE ->Write ();
316416 fMCScan_Q0vsTheta ->Write ();
317417 fMCScan_Q0vsE ->Write ();
318418 fDataGraph ->Write ();
319419
320- Measurement1D::Write (drawOpts);
321-
322420}
323421
324422double ElectronScattering_DurhamData::GetLikelihood () {
0 commit comments