StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Calib_SC_GL.C
1 //
3 // Calib_SC_GL.C
4 // Macro for calibrating the SpaceCharge and GridLeak
5 // distortion corrections.
6 //
7 // Author: G. Van Buren, BNL
8 // Jan 6, 2006
9 //
10 // Synopsis:
11 // The general idea is to
12 // 1) Determine the luminosity scaler which gives the smallest variance
13 // in a linear fit to measured SpaceCharge
14 // 2) Fit GridLeak using the selected luminosity scaler
15 // 3) Use the two distinct constraints on SC & GL from these two fits
16 // to arrive at initial values
17 // 4) Fit SpaceCharge & GridLeak simultaneously
18 // 5) Plot the fits
19 //
20 // Parameters:
21 // 1) Input specification with info on data to be calibrated
22 // 2) Any additional cuts to help calibrate (i.e. "run!=6082047")
23 // 3) A specific scaler identifier (if you know which scaler you want to use)
24 // Possible values:
25 // -1 : automatically determine which of the predefined scalers is best
26 // 0..n : force the calibration to use specific scaler with this id (see
27 // 'dets' in code below, or StDetectorDbMaker/St_spaceChargeCorC.cxx)
28 // -n : perform principal components analysis with up to 'n' scalers
29 // "det" : force the calibration to use a specific scaler, e.g. "zdcx"
30 // "det1:det2:..." : perform principal components analysis using a specific
31 // set of scalers, e.g. "zdcx:zdcx*zdcx" allows a function of
32 // sc = A*(zdcx-offset)+B*(zdcx^2)
33 // 4) Debug mode
34 // 0 : default, limited output
35 // 1 : outlier removal info and chisquare contours
36 // 2 : detailed fit results
37 // 3 : more detailed fit results
38 // 5) If specified, the cuts for including data in the plots versus luminosity;
39 // same as parameter 2 if omitted. A looser cut here than in parameter 2
40 // allows to see the predictive power of whether the fits describe data
41 // outside what was used in the calibration
42 //
43 // Details:
44 //
45 // The macro now only runs compiled, e.g. execute as Calib_SC_GL.C+(...)
46 //
47 // PCA does not determine errors on individual terms, so errors are reported
48 // on the overall scale, e.g. (1.0 +/- error)*(...)
49 //
50 // Plots show:
51 // - small brown points : single file (ntuple) results from data
52 // - violet diamonds : averages for each dataset
53 // - grey lines : independent fits for SpaceCharge & GridLeak
54 // - black lines : simultaneous SpaceCharge & GridLeak fits
55 // (sometimes obscures the grey lines)
56 // - red crosses & ellipses : final values and 1-sigma ellipse
57 // - blue lines : constraint on SCxGL and 1-sigma bands
58 // - pink lines : outlier removal bands
59 //
60 // There are three input specification methods.
61 // 1) 1 An empty (or 0) input argument will use data files already opened
62 // in ROOT, and the used calibrations are read from each data file
63 // and used to determine the different input parameter sets:
64 // root hists*/*
65 // .x Calib_SC_GL.C+
66 // 2) A '@' at the first character of the input argument will open all files
67 // matching a simple wildcard pattern (does not work well with multiple
68 // wildcards), and the used calibrations are read from each data file
69 // and used to determine the different input parameter sets:
70 // root Calib_SC_GL.C+("@hists/*")
71 // 3) The input specification can be the name of an input file, which should
72 // contain lines with the path/name (simple wildcards allowed) of dataset
73 // files. Each line beginning with a '@' character will be assumed to be
74 // a distinct dataset with common used calibrations:
75 // @histsSetA/*
76 // @histsSetB/*
77 // @histsSetC/*
78 // A # at the beginning of a line will cause that dataset to be skipped:
79 // #@histsSetC/*
80 // If this file is named input.dat, it could be analyzed with:
81 // root Calib_SC_GL.C+("input.dat")
82 // If you know you want to use only bbce+bbcw:
83 // root Calib_SC_GL.C+("input.dat","",4)
84 //
85 // Note: Older input file formats are supported for backward compatibility
86 // with hist files that do not have the used corrections stored in them.
87 // For these, each line of the file must contain the following info:
88 // - dataset file specification
89 // - scaler detector used in that dataset (see dets strings below)
90 // - SpaceCharge rate used in that dataset
91 // - GridLeak multiplier used in that dataset
92 // - (optional) SpaceCharge offset used in that dataset (set "!" below)
93 // Example 3 dataset file which used bbce+bbcw would look like this:
94 // histsSetA/* 4 1.7e-8 9.0
95 // histsSetB/* 4 1.7e-8 12.0
96 // histsSetC/* 4 1.7e-8 15.0
97 // A ! at the beginning of a file spec will cause the optional SpaceCharge
98 // offset to be read in on that line:
99 // !histsSetC/* 4 1.7e-8 15.0 550
100 // A # can be used to skip these old format lines as well
101 //
103 
104 #if !defined(__CINT__) || defined(__MAKECINT__)
105 #include "TFile.h"
106 #include "TKey.h"
107 #include "TCanvas.h"
108 #include "TChain.h"
109 #include "TChainElement.h"
110 #include "TGraphErrors.h"
111 #include "TF3.h"
112 #include "TH2.h"
113 #include "TProfile.h"
114 #include "TArrayD.h"
115 #include "TCut.h"
116 #include "TMarker.h"
117 #include "TPolyMarker.h"
118 #include "TEllipse.h"
119 #include "TFitter.h"
120 #include "TMinuit.h"
121 #include "TMath.h"
122 #include "TStyle.h"
123 #include "TString.h"
124 #include "TPrincipal.h"
125 #include "TROOT.h"
126 #include <Stiostream.h>
127 #endif
128 
129 // Switches for users? Needs better arrangement
130 Int_t EWmode = 0; // -1 east, 0 both, 1 west
131 Int_t EW_ASYMMETRY = kFALSE;
132 
133 // Switches and constants for experts
134 Bool_t USE_OLD_STDDEV = kTRUE; // method for determining STDDEV
135 Bool_t FORCE_GLO_SO = kTRUE; // force GLO = SO
136 Bool_t FORCE_g3_same = kTRUE;
137 Bool_t FORCE_g4_same = kTRUE;
138 Bool_t FORCE_g5_same = kTRUE; // if g4_same, then g3_same and g5_same are equivalent
139 Bool_t NO_PLOTS = kFALSE;
140 Bool_t BY_SECTOR = kFALSE;
141 Bool_t FIX_GL = kFALSE;
142 double GUESS_GL = 9.0; // auto-sets to 0.5 for runs>20000000 (iTPC) unless FIX_GL is true
143 double GUESS_g5 = 15.0;
144 double MAX_DEV = 4.0;
145 double REF_CONST1 = 0.318;
146 double RUN_CORRELATE = 0.70; // Assume that RUN_CORRELATE fraction of variance
147  // among results from same run/file is correlated
148  // (lower values may help with failed fits)
149 
150 // Main routine:
151 int Calib_SC_GL(const char* input=0, const char* cuts=0, int scaler=-1, int debug=0, const char* gcuts=0);
152 int Calib_SC_GL(const char* input, const char* cuts, const char* scalerstr, int debug=0, const char* gcuts=0);
153 
154 // Fitting functions:
155 Double_t funcGapf(Double_t* x, Double_t* pars);
156 Double_t funcSC(Double_t* x, Double_t* pars);
157 Double_t funcSC2(Double_t* x, Double_t* pars);
158 void fnchGapf(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
159 void fnchSC(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
160 void fnchSCGapf(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
161 void fnchSCGapfSec(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
162 
163 // Helper functions:
164 int Init(const char* input);
165 int Waiting();
166 void SetMinMax(int n, Double_t* ar, double& min, double& max, double inclusion, double margin=0.075);
167 int SetMinuitPars(const int n, TString* names, Double_t* starts, Double_t* steps, int debug, Bool_t* fixing=0);
168 void GetMinuitPar(const int n);
169 void GetMinuitPars(const int n, Bool_t* fixing=0, Double_t* fixedValues=0);
170 void Log2Lin(int i);
171 int FitWithOutlierRemoval(int npar, int debug);
172 void DrawErrorContours(int npar, Bool_t* fix=0);
173 TString PCA(int Nmax=32, int debug=0);
174 void PrintResult(double scp, double escp, double sop, double esop,
175  double glp, double eglp, double ewp, double eewp, const char* det);
176 void PrintResult(double scp, double escp, double sop, double esop,
177  double* glp, double* eglp, const char* det);
178 
179 // Global parameters:
180 const int nLimit=150;
181 const int nfip=128;
182 const int nsca=32;
183 const int npos=nfip*nsca;
184 const int nMeasuresMax=nfip*256;
185 const int nMeasuresHalf=nMeasuresMax/2;
186 int nMeasures = 0;
187 int nMeasuresI = 0;
188 int nfi = 0;
189 
190 // Global inputs
191 TCut cut,gcut;
192 TChain* SCi[nfip];
193 TChain* SCall = 0;
194 Bool_t allZeros = kTRUE;
195 TString scastr;
196 
197 // Graphics
198 TCanvas* cSummary = 0;
199 TCanvas* conSC[nsca];
200 TCanvas* conGL = 0;
201 TCanvas* conSCGL = 0;
202 TCanvas* cSC = 0;
203 TCanvas* cGL = 0;
204 TH1* histo = 0;
205 
206 // Colors
207 Int_t colorData = kOrange+2;
208 Int_t colorMiFit = kGray;
209 Int_t colorFit = kBlack;
210 Int_t colorConstraint = kBlue;
211 Int_t colorGrid = kSpring+8;
212 Int_t colorAggregate = kViolet+1;
213 Int_t colorOutlier = kMagenta-9;
214 Int_t colorFinal = kRed;
215 
216 // Fitting
217 TF1* scgl_fit = 0;
218 TVirtualFitter* minuit = 0;
219 TMinuit* minuit2 = 0;
220 double arglist[10];
221 Double_t* fitPars = 0;
222 Double_t* fitParErrs = 0;
223 char** parName = 0;
224 Double_t fitParsSC[3];
225 Double_t fitParErrsSC[3];
226 Double_t fitParsGL[3];
227 Double_t fitParErrsGL[3];
228 Double_t fitParsSCGL[9];
229 Double_t fitParErrsSCGL[9];
230 Double_t fitParsSCGLsec[19];
231 Double_t fitParErrsSCGLsec[19];
232 char* parNameSC[3];
233 char* parNameGL[3];
234 char* parNameSCGL[9];
235 char* parNameSCGLsec[19];
236 Bool_t unfittableSec[12];
237 bool atLimit[32];
238 
239 // Global fit data and parameters
240 Double_t m_sc[nMeasuresMax];
241 Double_t m_scS[12][nMeasuresMax];
242 Double_t m_sc2[nMeasuresMax];
243 Double_t m_sc2S[12][nMeasuresMax];
244 Double_t m_usc[nMeasuresMax];
245 Double_t m_ugl[nMeasuresMax];
246 Double_t m_gapf[nMeasuresMax];
247 Double_t m_gapfS[12][nMeasuresMax];
248 Double_t m_runs[nMeasuresMax];
249 Double_t m_L[nMeasuresMax];
250 Double_t m_c1[nMeasuresMax];
251 Int_t m_set[nMeasuresMax];
252 Int_t m_runIdx[nMeasuresMax];
253 Double_t devsSC[nMeasuresMax];
254 Double_t devsGL[nMeasuresMax];
255 Double_t maxvarSC[nMeasuresMax];
256 Double_t maxvarGL[nMeasuresMax];
257 Double_t devsSCsec[12][nMeasuresMax];
258 Double_t devsGLsec[12][nMeasuresMax];
259 Double_t devs_set[nfip];
260 Double_t* devs;
261 Bool_t outliers0[nMeasuresMax];
262 Bool_t outliersSC[nMeasuresMax];
263 Bool_t outliersGL[nMeasuresMax];
264 Bool_t* outliers;
265 Double_t nno[nfip]; // number of non-outliers per set
266 
267 double STDDEV = 1.0;
268 double CHISQ = 1.0;
269 double VAR = 1.0;
270 double VAR_all = 1.0;
271 double STDDEV_SC = 0.001;
272 double STDDEV_GL = 0.035;
273 
274 // Data for each chain
275 int glmode[nfip]; // mode of obtaining used GL
276 double ugl[nfip]; // GL used
277 
278 
279 // Plot data:
280 double xgl[nMeasuresMax];
281 double ygl[nMeasuresMax];
282 double xgo[nMeasuresMax];
283 double ygo[nMeasuresMax];
284 double xsc[nMeasuresMax];
285 double ysc[nMeasuresMax];
286 double yso[nMeasuresMax];
287 
288 // PCA globals
289 TString pcadets[nsca];
290 double pcacoef[nsca];
291 double pcaoffset;
292 int pcaN = -1;
293 const int maxp = 1024;
294 TPrincipal* ppl[maxp];
295 Bool_t DO_PCA = kFALSE;
296 Bool_t ITER0 = kFALSE;
297 
298 
299 
301 // Main routine
303 
304 int Calib_SC_GL(const char* input, const char* cuts, const char* scalerstr, int debug, const char* gcuts) {
305  scastr = scalerstr;
306  return (scastr.Contains(":") ?
307  Calib_SC_GL(input,cuts, -999,debug,gcuts) :
308  Calib_SC_GL(input,cuts,nsca-1,debug,gcuts) );
309 }
310 
311 int Calib_SC_GL(const char* input, const char* cuts, int scaler, int debug, const char* gcuts) {
312 
313  // Define useful luminosity scaler detectors
314  // copied from StRoot/StDetectorDbMaker/St_spaceChargeCorC.cxx
315  TString dets[nsca];
316  dets[0] = "vpdx" ;
317  dets[1] = "bbcx" ;
318  dets[2] = "zdcx" ;
319  dets[3] = "zdce+zdcw";
320  dets[4] = "bbce+bbcw";
321  dets[5] = "zdce" ;
322  dets[6] = "zdcw" ;
323  dets[7] = "bbce" ;
324  dets[8] = "bbcw" ;
325  dets[9] = "bbcyb" ;
326  dets[10]= "bbcbb" ;
327  dets[11]= "vpde" ;
328  dets[12]= "vpdw" ;
329  dets[13]= "zdcxnk" ;
330  dets[14]= "zdcenk" ;
331  dets[15]= "zdcwnk" ;
332  //dets[21] = "zdcc" ;
333  //dets[22] = "bbcc" ;
334  //dets[23] = "bbce+bbcw-4*(bbcyb+bbcbb)"; // Just for fun.
335  // reserve dets[nsca-1] for special uses (PCA, manually provided scalers)
336 
337  if (BY_SECTOR && EWmode==0) {
338  printf("Error: must do east or west only when processing by sector!");
339  return 1;
340  }
341  if (EW_ASYMMETRY && EWmode!=0) {
342  printf("Error: must do both east and west when processing full asymmetry!");
343  return 1;
344  }
345  if (EW_ASYMMETRY && BY_SECTOR) {
346  printf("Error: by sector not yet implemented when processing full asymmetry!");
347  return 1;
348  }
349 
350  int status = Init(input);
351  if (status) return status;
352 
353  cut = ((cuts) ? cuts : "");
354  gcut = ((gcuts) ? gcuts : cut.GetTitle());
355 
356  DO_PCA = (scaler < -1);
357 
358  if (DO_PCA) { // Doing a PCA analysis
359  if (!ITER0) {
360  Bool_t no_plots = NO_PLOTS;
361  NO_PLOTS = kTRUE; // Don't do plots for iteration 0
362  ITER0 = kTRUE;
363  printf("\n*** Running PCA iteration 0 ***\n\n");
364  // Iteration 0 uses a single input dataset to define PCA
365  status = Calib_SC_GL(input,cuts,scaler,debug,gcuts);
366  if (status) return status;
367  // Iteration 1 uses all input datasets to define PCA
368  NO_PLOTS = no_plots;
369  ITER0 = kFALSE;
370  printf("\n*** Running PCA iteration 1 ***\n\n");
371  delete minuit; // Need to delete the minuit instance that was used in iteration 0
372  } else {
373  // This line is for iteration 0 only
374  fitPars = fitParsSCGL;
375  }
376  dets[nsca-1] = PCA(1-scaler,debug); // dets[nsca-1] reserved for such use
377  scaler = nsca-1;
378  } else {
379  if (scastr.Length()) {
380  dets[nsca-1] = scastr; // dets[nsca-1] reserved for such use
381  scaler = nsca - 1;
382  }
383  }
384 
385  int i,j,k;
386  double temp1,vsc_min = 1e10;
387  int jmin=-1;
388 
389  // Asymmetrical calibration mode variables
390  TString EWstr = (EWmode < 0 ? " EAST" : (EWmode > 0 ? " WEST" : ""));
391  TString scvarstr = (EWmode < 0 ? "sce" : (EWmode > 0 ? "scw" : (EW_ASYMMETRY ? "scw" : "sc")));
392  TString uscvarstr = (EWmode < 0 ? "usce" : "usc"); // "usc" is same for west
393  const char* scvar = scvarstr.Data();
394  const char* uscvar = uscvarstr.Data();
395  if (EWmode != 0) printf("\nAsymmetrical calibration mode:%s\n\n",EWstr.Data());
396 
397 
398  // Data arrays
399  double asg[nfip];
400  double glk[nfip];
401  double ago[nfip];
402  double glo[nfip];
403  double spc[nfip];
404  double spo[nfip];
405  double sce[nsca];
406  double sof[nsca];
407  double ggl[nsca];
408  double sceE[nsca];
409  double sofE[nsca];
410  double gglE[nsca];
411  double vsc[nsca];
412  double ssc[nsca];
413  double vgl;
414 
415  // Define useful graphics tools
416  TMarker mark;
417  mark.SetMarkerColor(colorFinal);
418  mark.SetMarkerStyle(28);
419  TPolyMarker mark2;
420  TEllipse ellip;
421  ellip.SetLineColor(colorFinal);
422  ellip.SetLineStyle(2);
423  ellip.SetFillStyle(0);
424  double zero = 0;
425 
426 
427  // Prepare for SC fit
428  devs = devsSC;
429  outliers = outliers0;
430  fitPars = fitParsSC;
431  fitParErrs = fitParErrsSC;
432  parName = parNameSC;
433 
434 
435  // Set some dimensional variables
436  Int_t fittableCnt[12];
437  memset(fittableCnt,0,12*sizeof(Int_t));
438  int startSec = (EWmode < 0 ? 13 : 1);
439  nMeasures = 0;
440  for (i=0;i<nfi;i++) { // parameter sets
441  TString SCvarstr = (EW_ASYMMETRY ? "usc:usce:scw:sce" : Form("%s:%s:gapf",uscvar,scvar));
442  nMeasuresI = SCi[i]->Draw(SCvarstr.Data(),cut,"goff");
443  if (EW_ASYMMETRY) {
444  // east values offset by nMeasuresHalf
445  memcpy(&(m_usc [nMeasures ]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
446  memcpy(&(m_usc [nMeasures+nMeasuresHalf]),SCi[i]->GetV2(),nMeasuresI*sizeof(Double_t));
447  memcpy(&(m_sc [nMeasures ]),SCi[i]->GetV3(),nMeasuresI*sizeof(Double_t));
448  memcpy(&(m_sc [nMeasures+nMeasuresHalf]),SCi[i]->GetV4(),nMeasuresI*sizeof(Double_t));
449  Bool_t useGapfew = (SCi[i]->Draw("1",cut&&"gapfw!=0||gapfe!=0","goff") > 0);
450  SCi[i]->Draw("gapfw:gapfe:gapf",cut,"goff");
451  if (useGapfew) {
452  // gapfe & gapfw are filled (non-zero)
453  memcpy(&(m_gapf[nMeasures ]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
454  memcpy(&(m_gapf[nMeasures+nMeasuresHalf]),SCi[i]->GetV2(),nMeasuresI*sizeof(Double_t));
455  } else {
456  printf("Set %3d: no e/w gapf, using global gapf\n",i);
457  memcpy(&(m_gapf[nMeasures ]),SCi[i]->GetV3(),nMeasuresI*sizeof(Double_t));
458  memcpy(&(m_gapf[nMeasures+nMeasuresHalf]),SCi[i]->GetV3(),nMeasuresI*sizeof(Double_t));
459  }
460  } else {
461  memcpy(&(m_usc [nMeasures]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
462  memcpy(&(m_sc [nMeasures]),SCi[i]->GetV2(),nMeasuresI*sizeof(Double_t));
463  memcpy(&(m_gapf[nMeasures]),SCi[i]->GetV3(),nMeasuresI*sizeof(Double_t));
464  }
465  if (BY_SECTOR) {
466  for (j=0;j<12;j++) {
467  k = startSec + j;
468  SCvarstr = Form("sc%d:gapf%d",k,k);
469  TCut testFittable = Form("sc%d!=0&&gapf%d!=0",k,k);
470  Int_t fittableCntS = SCi[i]->Draw(SCvarstr.Data(),cut&&testFittable,"goff");
471  if (fittableCntS > 0) {
472  fittableCnt[j] += fittableCntS;
473  SCi[i]->Draw(SCvarstr.Data(),cut,"goff");
474  memcpy(&(m_scS [j][nMeasures]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
475  memcpy(&(m_gapfS[j][nMeasures]),SCi[i]->GetV2(),nMeasuresI*sizeof(Double_t));
476  } else {
477  if (debug>0) printf("Set %3d : Unfittable sector: %d\n",i,k);
478  memset(&(m_scS [j][nMeasures]),0,nMeasuresI*sizeof(Double_t));
479  memset(&(m_gapfS[j][nMeasures]),0,nMeasuresI*sizeof(Double_t));
480  }
481  }
482  }
483  SCi[i]->Draw("const1:run",cut,"goff");
484  memcpy(&(m_c1 [nMeasures]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
485  memcpy(&(m_runs[nMeasures]),SCi[i]->GetV2(),nMeasuresI*sizeof(Double_t));
486  // possibly improve to (event+run*1e-6) as a more unique ID
487  // for the case of two files from the same run
488  if (glmode[i] == 2) {
489  SCi[i]->Draw("ugl",cut,"goff");
490  memcpy(&(m_ugl[nMeasures]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
491  ugl[i] = m_ugl[nMeasures];
492  } else for (k=0; k<nMeasuresI; k++) m_ugl[k+nMeasures] = ugl[i]; // if not in the ntuple
493  if (!FIX_GL && (m_runs[0] > 20000000)) GUESS_GL = 0.5;
494  for (k=0; k<nMeasuresI; k++) {
495  m_set[k+nMeasures] = i;
496  m_runIdx[k+nMeasures] = k+nMeasures;
497  for (j=0; j<k+nMeasures; j++) { // find identical input runs/files
498  if (m_runs[k+nMeasures] == m_runs[j]) { m_runIdx[k+nMeasures] = j; break; }
499  }
500  }
501  nMeasures += nMeasuresI;
502  }
503  for (i=0;BY_SECTOR && i<12;i++) {
504  unfittableSec[i] = (fittableCnt[i] < 8);
505  if (unfittableSec[i]) printf("Unfittable sector: %d\n",startSec + i);
506  }
507 
508  // Loop over available scaler detectors
509  for (j=0;j<nsca;j++) {
510  if (scaler>=0 && j!=scaler) continue;
511  if (dets[j].Length()<1) continue;
512  const char* dt = dets[j].Data();
513 
514 
516  // 2D fit for SpaceCharge
517 
518  // Set dimensional variables
519  nMeasures = 0;
520  for (i=0;i<nfi;i++) { // parameter sets
521  SCi[i]->Draw(dt,cut,"goff");
522  nMeasuresI = SCi[i]->GetSelectedRows();
523  memcpy(&(m_L[nMeasures]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
524  nMeasures += nMeasuresI;
525  }
526  STDDEV = STDDEV_SC;
527 
528  // Set starting values and step sizes for parameters
529  double sce_init = GUESS_g5 + GUESS_GL; // approximate guess on (g5 + GL)
530  if (DO_PCA) {
531  if (!ITER0) {
532  // Do a sanity check on g5 from iteration 0 before using it
533  if (TMath::Abs(TMath::Log(fitParsSCGL[1]/GUESS_g5)) > 0.9) {
534  if (debug>0) printf("\nWARNING: g5 from iteration 0 at or near limit - not sane\n");
535  sce_init = GUESS_g5 + fitParsSCGL[4];
536  } else {
537  sce_init = fitParsSCGL[1] + fitParsSCGL[4];
538  }
539  }
540  } else {
541  // Small denominators can cause large values that throw off the mean
542  // where small is in the context of the RMS of the distribution instead
543  // of the range, which is even more sensitive to outliers than the RMS
544  SCall->Draw(Form("%s",dt),cut,"goff");
545  histo = SCall->GetHistogram();
546  double dtrms = histo->GetRMS();
547  TCut avoidSmallDenominators = Form("abs(%s)>%g",dt,dtrms * 0.4); // 0.4 found by trial & error
548  SCall->Draw(Form("%s/(%s)",scvar,dt),avoidSmallDenominators&&cut,"goff");
549  histo = SCall->GetHistogram();
550  sce_init *= (histo->GetMean());
551  }
552  if (debug) printf("\nUsing initial value for SCe of %g\n",sce_init);
553 
554  TString sname[3] = {"SO","log(SCe)","log(g5)"};
555  Double_t sstart[3] = {100., TMath::Log(sce_init), TMath::Log(GUESS_g5)};
556  Double_t sstep[3] = {1000., TMath::Abs(0.01*TMath::Log(sce_init)), 0.01};
557  if (DO_PCA) { sstart[0] = 1e-3; sstep[0] = 1e-5; }
558  SetMinuitPars(3,sname,sstart,sstep,debug);
559  minuit->SetFCN(fnchSC); // Keep g3r=g4r=1 to stabilize initial fit
560 
561  // Perform the fit
562  printf("\nSpaceCharge fit results {scaler: %s}:\n",dt);
563  status = FitWithOutlierRemoval(3,debug);
564  if (status==4) {
565  if (debug) printf("Fit failed with initial g5=15. Trying 12...\n");
566  GUESS_g5 = 12.0;
567  sstart[2] = TMath::Log(GUESS_g5);
568  SetMinuitPars(3,sname,sstart,sstep,debug);
569  status = FitWithOutlierRemoval(3,debug);
570  }
571  if (status==4) {
572  if (debug) printf("Fit failed with initial g5=12. Trying 20...\n");
573  GUESS_g5 = 20.0;
574  sstart[2] = TMath::Log(GUESS_g5);
575  SetMinuitPars(3,sname,sstart,sstep,debug);
576  status = FitWithOutlierRemoval(3,debug);
577  }
578  if (status) {
579  printf("Fit failed for sc, err = %d\nTrying next scaler (if any)...\n\n",status);
580  Waiting();
581  continue;
582  }
583  if (!NO_PLOTS && debug>0) {
584  if (conSC[j] ==0) {
585  conSC[j] = new TCanvas(Form("conSC_%d",j),Form("SC fit contours for %s",dt),600,400);
586  conSC[j]->Divide(2,2);
587  } else conSC[j]->cd();
588  DrawErrorContours(3);
589  }
590  Log2Lin(1);
591  Log2Lin(2);
592  sof[j] = fitPars[0];
593  sce[j] = fitPars[1];
594  ggl[j] = fitPars[2];
595  sofE[j] = fitParErrs[0];
596  sceE[j] = fitParErrs[1];
597  gglE[j] = fitParErrs[2];
598  vsc[j] = VAR; // Excludes outliers
599  ssc[j] = STDDEV;
600 
601  if (debug>0) printf("VAR_all = %6.4g for %s\n",VAR_all,dt);
602  if (VAR_all < vsc_min) { // Select best scaler based on variance without excluding outliers
603  vsc_min = VAR_all;
604  jmin = j;
605  memcpy(outliersSC,outliers0,nMeasures*sizeof(Double_t));
606  }
607 
608  }
609 
611  // Set up for using the best scaler
612 
613  if (jmin<0) { printf("ERROR - no good scaler found\n"); return 1; }
614  const char* detbest = dets[jmin].Data();
615  if (scaler<0) { // were looking for the best scaler
616  printf("*** Best scaler = %s [ID = %d]\n\n",detbest,jmin);
617  }
618  fitParsSC[0] = sof[jmin];
619  fitParsSC[1] = sce[jmin];
620  fitParsSC[2] = ggl[jmin];
621  fitParErrsSC[0] = sofE[jmin];
622  fitParErrsSC[1] = sceE[jmin];
623  fitParErrsSC[2] = gglE[jmin];
624  STDDEV_SC = ssc[jmin];
625 
626 
628  // 3D fit for GridLeak
629 
630  // Do this fit even if fixing GL to get an estimate for g2
631 
632  // Prepare for GL fit
633  devs = devsGL;
634  outliers = outliersGL;
635  fitPars = fitParsGL;
636  fitParErrs = fitParErrsGL;
637  parName = parNameGL;
638  STDDEV = STDDEV_GL;
639 
640  // Set dimensional variables
641  nMeasures = 0;
642  for (i=0;i<nfi;i++) {
643  SCi[i]->Draw(detbest,cut,"goff");
644  nMeasuresI = SCi[i]->GetSelectedRows();
645  memcpy(&(m_L[nMeasures]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
646  nMeasures += nMeasuresI;
647  }
648 
649  // Set starting values and step sizes for parameters
650  double scale_init = GUESS_GL / (ggl[jmin] + GUESS_GL); // SC*GL = SC*(g5+GL) * GL/(g5+GL)
651  double scXgl = sce[jmin]*scale_init;
652  double escXgl = sceE[jmin]*scale_init;
653  double GLO = sof[jmin];
654  double eGLO = sofE[jmin];
655  TString gname[3] = {"g2","g1/g2 = SCxGL","GLO"};
656  Double_t gstart[3] = {1.0, scXgl, GLO};
657  Double_t gstep[3] = {0.01, escXgl, eGLO};
658  if (DO_PCA) {
659  gstart[1] = GUESS_GL;
660  gstep[1] = 0.1*GUESS_GL;
661  }
662  SetMinuitPars(3,gname,gstart,gstep,debug);
663  minuit->SetFCN(fnchGapf);
664 
665  // Perform the fit
666  printf("\nGridLeak fit results {scaler: %s}:\n",detbest);
667  status = FitWithOutlierRemoval(3,debug);
668  if (status) {
669  printf("Fit failed for gapf, err = %d\n",status);
670  if (!FIX_GL) return 1;
671  printf("Continuing with a fixed GL.");
672  for (k=0;k<3;k++) { fitPars[k] = gstart[k]; fitParErrs[k] = gstep[k]; }
673  }
674  if (!NO_PLOTS && debug>0) {
675  if (conGL ==0) {
676  conGL = new TCanvas("conGL","GL fit contours",600,400);
677  conGL->Divide(2,2);
678  } else conGL->cd();
679  DrawErrorContours(3);
680  }
681 
682  vgl = VAR;
683  STDDEV_GL = STDDEV;
684 
685  if (FIX_GL) {
686  // use scXgl and GLO estimates from before the 3D fit for GridLeak
687  escXgl = 0;
688  } else {
689  scXgl = fitPars[1];
690  escXgl = fitParErrs[1];
691  GLO = fitPars[2];
692  eGLO = fitParErrs[2];
693  }
694 
695 
697  // Combining SC & GL fits...we have our first guesses of solutions!!!
698 
699  double scp = (sce[jmin]-scXgl)/ggl[jmin];
700  double escp = TMath::Sqrt(sceE[jmin]*sceE[jmin] +
701  escXgl*escXgl +
702  scp*scp*gglE[jmin]*gglE[jmin]
703  ) / ggl[jmin];
704 
705  double sop = sof[jmin];
706  double esop = sofE[jmin];
707  double ewp = 1.0;
708  double eewp = 0.0;
709 
710  double scale_final = scXgl/(sce[jmin] - scXgl);
711  double glp = ggl[jmin]*scale_final;
712  double eglp = 0;
713  if (!FIX_GL) {
714  scale_final *= glp;
715  eglp = TMath::Sqrt(TMath::Power(glp*gglE[jmin]/ggl[jmin],2) +
716  TMath::Power(scale_final*escXgl*sce[jmin]/(scXgl*scXgl),2) +
717  TMath::Power(scale_final*sceE[jmin]*scXgl,2));
718  }
719  if (debug>0) {
720  printf("\n*** FIRST PASS CALIBRATION VALUES: ***\n");
721  PrintResult(scp, escp, sop, esop, glp, eglp, ewp, eewp, detbest);
722  printf("USING STDDEV sc => %f :: gapf => %f\n",STDDEV_SC,STDDEV_GL);
723  }
724 
726  // Try again with one unified full 3D fit
727  // No further outlier removal (use already-determined outliers)
728 
729  // Prepare for SC+GL fit
730  fitPars = fitParsSCGL;
731  fitParErrs = fitParErrsSCGL;
732  parName = parNameSCGL;
733 
734  // Set starting values and step sizes for parameters
735  int npar = 9;
736  TString fname[9] = {"g2","log(g5)","log(SC)","SO","log(GL)","GLO","log(g5r)","log(g4r)","log(ewratio)"};
737  Double_t fstart[9] = {fitParsGL[0], TMath::Log(fitParsSC[2]), TMath::Log(scp),
738  sop, TMath::Log(GUESS_GL), GLO, 0, 0, TMath::Log(ewp)};
739  Double_t fstep[9] = {fitParErrsGL[0], fitParErrsSC[2]/fitParsSC[2], escp/scp,
740  esop, 0.1, eGLO, 0.001, 0.001, 0.001};
741  Bool_t ffix[9] = {kFALSE, kFALSE, kFALSE, kFALSE, FIX_GL,
742  FORCE_GLO_SO, (FORCE_g3_same || FORCE_g5_same), FORCE_g4_same, !EW_ASYMMETRY};
743  SetMinuitPars(9,fname,fstart,fstep,debug,ffix);
744  minuit->SetFCN(fnchSCGapf);
745 
746  // Perform the fit
747  printf("\nSpaceCharge & GridLeak fit results {scaler: %s}:\n",detbest);
748  status = minuit->ExecuteCommand("MINIMIZE", arglist, 1);
749  if (status) {
750  printf("Fit failed for sc+gl, err = %d\n",status);
751  return 1;
752  }
753  // covarAdjust is a fudge factor for how the covariance (which is negative)
754  // between pars 2 and 4 modifies the variance on the product of their exponentials
755  Double_t covarAdjust = 1;
756  if (!(ffix[2] || ffix[4]))
757  covarAdjust += (minuit->GetCovarianceMatrixElement(2,4) /
758  (minuit->GetCovarianceMatrixElement(2,2) +
759  minuit->GetCovarianceMatrixElement(4,4)));
760  GetMinuitPars(npar,ffix,fstart);
761  status = minuit->ExecuteCommand("SET LIM", 0, 0); // Remove limits before MINOS
762  if (status) {
763  printf("SetLimit failed for sc+gl, err = %d\n",status);
764  return 1;
765  }
766  double eplus,eminus,eparab,gcc;
767  for (k=0;k<npar;k++) {
768  if (ffix[k]) { fitPars[k] = fstart[k]; fitParErrs[k] = 0; continue; }
769  if (atLimit[k]) {
770  printf("%s\t:\t%g\t+/- %g\t(AT LIMIT!)\n",parName[k],fitPars[k],fitParErrs[k]);
771  continue;
772  }
773  for (i=0;i<npar;i++) if (i!=k && !ffix[i]) minuit->FixParameter(i);
774  status = minuit->ExecuteCommand("MINIMIZE", arglist, 1);
775  if (status) {
776  printf("Fit failed for sc+gl, err = %d\n",status);
777  return 1;
778  }
779  arglist[1] = (double) k+1;
780  status = minuit->ExecuteCommand("MINOS", arglist, 2);
781  if (status) {
782  printf("Fit errors failed for sc+gl, err = %d\n",status);
783  return 1;
784  }
785  GetMinuitPar(k);
786  minuit2->mnerrs(k,eplus,eminus,eparab,gcc);
787  if (eplus!=0.0 && eminus!=0.0) fitParErrs[k] = 0.5*(eplus-eminus);
788  else if (eplus!=0.0 || eminus!=0.0) fitParErrs[k] = eparab;
789  if (debug>0) printf("%s\t:\t%g\t+%g/%g or +/- %g\n",parName[k],fitPars[k],eplus,eminus,eparab);
790  printf("%s\t:\t%g\t+/- %g\n",parName[k],fitPars[k],fitParErrs[k]);
791  for (i=0;i<npar;i++) if (i!=k && !ffix[i]) minuit->ReleaseParameter(i);
792  }
793  if (!NO_PLOTS && debug>0) {
794  if (conSCGL ==0) {
795  conSCGL = new TCanvas("conSCGL","SCGL fit contours",600,900);
796  conSCGL->Divide(3,5);
797  } else conSCGL->cd();
798  DrawErrorContours(npar,ffix);
799  }
800  Double_t bstart[19] = {fitPars[0], fitPars[1], fitPars[2], fitPars[3],
801  fitPars[4], fitPars[4], fitPars[4], fitPars[4], fitPars[4], fitPars[4],
802  fitPars[4], fitPars[4], fitPars[4], fitPars[4], fitPars[4], fitPars[4],
803  fitPars[5], fitPars[6], fitPars[7]};
804  Double_t bstep[19] = {fitParErrs[0], fitParErrs[1], fitParErrs[2],
805  fitParErrs[3], fitParErrs[4], fitParErrs[4], fitParErrs[4],
806  fitParErrs[4], fitParErrs[4], fitParErrs[4], fitParErrs[4],
807  fitParErrs[4], fitParErrs[4], fitParErrs[4], fitParErrs[4],
808  fitParErrs[4], fitParErrs[5], fitParErrs[6], fitParErrs[7]};
809  Log2Lin(1);
810  Log2Lin(2);
811  Log2Lin(4);
812  Log2Lin(6);
813  Log2Lin(7);
814  Log2Lin(8);
815  scp = fitPars[2];
816  sop = fitPars[3];
817  glp = fitPars[4];
818  escp = fitParErrs[2];
819  esop = fitParErrs[3];
820  eglp = fitParErrs[4];
821  if (FORCE_GLO_SO) {
822  GLO = sop;
823  eGLO = esop;
824  fitPars[5] = GLO;
825  fitParErrs[5] = eGLO;
826  } else {
827  GLO = fitPars[5];
828  eGLO = fitParErrs[5];
829  }
830  if (FORCE_g4_same) {
831  fitPars[7] = 1.0;
832  fitParErrs[7] = 0.0;
833  }
834  if (FORCE_g5_same) {
835  fitPars[6] = 1.0;
836  fitParErrs[6] = 0.0;
837  } else if (FORCE_g3_same) {
838  fitPars[6] = 1.0/fitPars[7];
839  fitParErrs[6] = fitPars[6]*(fitParErrs[7]/fitPars[7]);
840  }
841  if (EW_ASYMMETRY) {
842  ewp = fitPars[8];
843  eewp = fitParErrs[8];
844  } else {
845  ewp = 1.0;
846  eewp = 0.0;
847  }
848  double scXgl_final = scp*glp;
849  double escXgl_final = scXgl_final*TMath::Sqrt(covarAdjust*
850  (TMath::Power(fitParErrs[4]/fitPars[4],2)+
851  TMath::Power(fitParErrs[2]/fitPars[2],2)));
852  printf(Form("\n***%s FINAL CALIBRATION VALUES: ***\n",EWstr.Data()));
853  PrintResult(scp, escp, sop, esop, glp, eglp, ewp, eewp, detbest);
854 
855 
856 
858  // FIT GL BY SECTOR
860 
861  double glpS[12];
862  double eglpS[12];
863  double scpS, sopS, escpS, esopS, GLOS, eGLOS;
864  scpS = 0;
865  sopS = 0;
866  if (BY_SECTOR) {
867 
868  // Try again with one unified full 3D fit
869  // No further outlier removal (use already-determined outliers)
870 
871  // Prepare for SC+GL fit
872  fitPars = fitParsSCGLsec;
873  fitParErrs = fitParErrsSCGLsec;
874  parName = parNameSCGLsec;
875 
876  // Set starting values and step sizes for parameters
877  npar = 19;
878  TString bname[19] = {"g2","log(g5)","log(SC)","SO","log(GLa)","log(GLb)",
879  "log(GLc)", "log(GLd)","log(GLe)","log(GLf)","log(GLg)","log(GLh)",
880  "log(GLi)","log(GLj)","log(GLk)","log(GLl)","GLO","log(g5r)","log(g4r)"};
881  Bool_t bfix[19] = {kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE,
882  kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE, kFALSE,
883  FORCE_GLO_SO, (FORCE_g3_same || FORCE_g5_same), FORCE_g4_same};
884  for (k=0;k<12;k++) bfix[k+4] = unfittableSec[k];
885  SetMinuitPars(19,bname,bstart,bstep,debug,bfix);
886  minuit->SetFCN(fnchSCGapfSec);
887 
888  // Perform the fit
889  printf("\nSpaceCharge & GridLeak fit by sector results {scaler: %s}:\n",detbest);
890  status = minuit->ExecuteCommand("MINIMIZE", arglist, 1);
891  if (status) {
892  printf("Fit failed for sc+gl by sector, err = %d\n",status);
893  return 1;
894  }
895  GetMinuitPars(npar,bfix,bstart);
896  status = minuit->ExecuteCommand("SET LIM", 0, 0); // Remove limits before MINOS
897  if (status) {
898  printf("SetLimit failed for sc+gl by sector, err = %d\n",status);
899  return 1;
900  }
901  for (k=0;k<npar;k++) {
902  if (bfix[k]) { fitPars[k] = bstart[k]; fitParErrs[k] = 0; continue; }
903  if (atLimit[k]) {
904  printf("%s\t:\t%g\t+/- %g\t(AT LIMIT!)\n",parName[k],fitPars[k],fitParErrs[k]);
905  continue;
906  }
907  for (i=0;i<npar;i++) if (i!=k && !bfix[i]) minuit->FixParameter(i);
908  status = minuit->ExecuteCommand("MINIMIZE", arglist, 1);
909  if (status) {
910  printf("Fit failed for sc+gl by sector, err = %d\n",status);
911  return 1;
912  }
913  arglist[1] = (double) k+1;
914  status = minuit->ExecuteCommand("MINOS", arglist, 2);
915  if (status) {
916  printf("Fit errors failed for sc+gl by sector, err = %d\n",status);
917  return 1;
918  }
919  GetMinuitPar(k);
920  minuit2->mnerrs(k,eplus,eminus,eparab,gcc);
921  if (eplus!=0.0 && eminus!=0.0) fitParErrs[k] = 0.5*(eplus-eminus);
922  else if (eplus!=0.0 || eminus!=0.0) fitParErrs[k] = eparab;
923  if (debug>0) printf("%s\t:\t%g\t+%g/%g or +/- %g\n",parName[k],fitPars[k],eplus,eminus,eparab);
924  printf("%s\t:\t%g\t+/- %g\n",parName[k],fitPars[k],fitParErrs[k]);
925  for (i=0;i<npar;i++) if (i!=k && !bfix[i]) minuit->ReleaseParameter(i);
926  }
927  Log2Lin(1);
928  Log2Lin(2);
929  Log2Lin(4);
930  Log2Lin(5);
931  Log2Lin(6);
932  Log2Lin(7);
933  Log2Lin(8);
934  Log2Lin(9);
935  Log2Lin(10);
936  Log2Lin(11);
937  Log2Lin(12);
938  Log2Lin(13);
939  Log2Lin(14);
940  Log2Lin(15);
941  Log2Lin(17);
942  Log2Lin(18);
943 
944  scpS = fitPars[2];
945  sopS = fitPars[3];
946  for (i=0;i<12;i++) {
947  glpS[i]=fitPars[4+i];
948  eglpS[i]=fitParErrs[4+i];
949  }
950  escpS = fitParErrs[2];
951  esopS = fitParErrs[3];
952  if (FORCE_GLO_SO) {
953  GLOS = sop;
954  eGLOS = esop;
955  fitPars[16] = GLOS;
956  fitParErrs[16] = eGLOS;
957  } else {
958  GLOS = fitPars[16];
959  eGLOS = fitParErrs[16];
960  }
961  if (FORCE_g4_same) {
962  fitPars[18] = 1.0;
963  fitParErrs[18] = 0.0;
964  }
965  if (FORCE_g5_same) {
966  fitPars[17] = 1.0;
967  fitParErrs[17] = 0.0;
968  } else if (FORCE_g3_same) {
969  fitPars[17] = 1.0/fitPars[18];
970  fitParErrs[17] = fitPars[17]*(fitParErrs[18]/fitPars[18]);
971  }
972  printf(Form("\n***%s FINAL BY-SECTOR CALIBRATION VALUES: ***\n",EWstr.Data()));
973  PrintResult(scpS, escpS, sopS, esopS, glpS, eglpS, detbest);
974 
975 
976 
977  }
978 
979 
981 
982 
983 
984 
985 
986 
987 
988  if (NO_PLOTS) return 0;
989 
990 
992  // Individual datasets compared to their fits
993  // Magenta lines show outlier exclusion zones
994  // gapf is presented as if no correction was applied,
995  // because it would otherwise depend on both the
996  // used luminosity scaler and the estimated best one
997 
998  if (fitParsSCGL[6] == 1.0 && fitParsSCGL[7] == 1.0) {
999  memcpy(m_sc2,m_sc,nMeasures*sizeof(Double_t));
1000  for (k=0; k<12; k++)
1001  memcpy(m_sc2S[k],m_scS[k],nMeasures*sizeof(Double_t));
1002  } else {
1003  // Need to shift data points for display purposes
1004  // (cannot shift the fit curves point-by-point)
1005  for (i=0; i<nMeasures; i++) {
1006  m_sc2[i] = m_sc[i] - m_usc[i] *
1007  (1.0 - ( ( fitParsSCGL[1] + m_ugl[i])/
1008  (fitParsSCGL[7]*(fitParsSCGL[6]*fitParsSCGL[1] + m_ugl[i]))));
1009  if (!BY_SECTOR) continue;
1010  for (k=0; k<12; k++)
1011  m_sc2S[k][i] = m_scS[k][i] - m_usc[i] *
1012  (1.0 - ( ( fitParsSCGLsec[1] + m_ugl[i])/
1013  (fitParsSCGLsec[18]*(fitParsSCGLsec[17]*fitParsSCGLsec[1] + m_ugl[i]))));
1014  }
1015  }
1016 
1017  double sc_min[12],sc_max[12],sc_offset[12],gapf_min[12],gapf_max[12],gapf_offset[12],Lmin,Lmax;
1018  SetMinMax(nMeasures,m_L,Lmin,Lmax,m_L[0]);
1019  for (k=0; k<(BY_SECTOR?12:1); k++) {
1020  SetMinMax(nMeasures,(BY_SECTOR?m_sc2S[k]:m_sc2),sc_min[k],sc_max[k],0);
1021  // Offset is minimum of either 0.2*(max-min) or 5*MAX_DEV*VAR,
1022  sc_offset[k] = 0.002*TMath::Max(1,
1023  TMath::Max(TMath::Nint((sc_max[k]-sc_min[k])*0.2/0.002),
1024  TMath::Nint(MAX_DEV*vsc[jmin]*5./0.002)));
1025  sc_max[k] += (nfi-1)*sc_offset[k];
1026  SetMinMax(nMeasures,(BY_SECTOR?m_gapfS[k]:m_gapf),gapf_min[k],gapf_max[k],0,0.25);
1027  gapf_max[k] += (BY_SECTOR ?
1028  fitParsSCGLsec[0]*scpS*glpS[k]*(Lmax-sopS) :
1029  fitParsSCGL[0]*scp*glp*(Lmax-sop) );
1030  gapf_offset[k] = 0.05*TMath::Max(1,
1031  TMath::Max(TMath::Nint((gapf_max[k]-gapf_min[k])*0.2/0.05),
1032  TMath::Nint(MAX_DEV*vgl*5./0.05)));
1033  gapf_max[k] += (nfi-1)*gapf_offset[k];
1034  }
1035 
1036  // Same as funcGapf and funcSC except non-luminosity dimensions become
1037  // extra parameters, and we're adjusting gapf back to uncorrected values
1038  TF1* fnGapf = 0;
1039  TF1* fnSC = 0;
1040  if (BY_SECTOR) {
1041  fnGapf = new TF1("fnGapf","[0] * ( ([2]*[4]) * (x - [16]) ) + [10]",Lmin,Lmax);
1042  fnSC = new TF1("fnSC","([2]*([1] + [4])) * (x - [3]) / ([18]*([17]*[1] + [8])) + [10]",Lmin,Lmax);
1043  fnGapf->SetParameters(fitParsSCGLsec);
1044  fnSC->SetParameters(fitParsSCGLsec);
1045  } else {
1046  fnGapf = new TF1("fnGapf","[0] * ( ([2]*[4]) * (x - [5]) ) + [10]",Lmin,Lmax);
1047  fnSC = new TF1("fnSC","([2]*([1] + [4])) * (x - [3]) / ([7]*([6]*[1] + [8])) + [10]",Lmin,Lmax);
1048  fnGapf->SetParameters(fitParsSCGL);
1049  fnSC->SetParameters(fitParsSCGL);
1050  }
1051  fnGapf->SetLineWidth(1);
1052  fnSC->SetLineWidth(1);
1053  TF1* fiGapf = new TF1("fiGapf","[0] * ( ([1]) * (x - [2]) ) + [10]",Lmin,Lmax);
1054  TF1* fiSC = new TF1("fiSC","[1] * (x - [0]) / ([2] + [8]) + [10]",Lmin,Lmax);
1055  fiGapf->SetParameters(fitParsGL);
1056  fiSC->SetParameters(fitParsSC);
1057  fiGapf->SetLineWidth(1);
1058  fiSC->SetLineWidth(1);
1059  fiGapf->SetLineColor(colorMiFit);
1060  fiSC->SetLineColor(colorMiFit);
1061  if (BY_SECTOR) {
1062  if (!cGL) {
1063  cGL = new TCanvas("cGL","GridLeak Fits by Sector",30,30,1000,750);
1064  cGL->Divide(4,3);
1065  }
1066  if (!cSC) {
1067  cSC = new TCanvas("cSC","SpaceCharge Fits",60,60,1000,750);
1068  cSC->Divide(4,3);
1069  }
1070  TH2D* htGLsec[12];
1071  TH2D* htSCsec[12];
1072  for (k=0; k<12; k++) {
1073  if (unfittableSec[k]) continue;
1074  htGLsec[k] = new TH2D(Form("htGL%d",k),
1075  Form("Sector %d: adjusted #font[32]{gapf} vs. %s for all sets, offset by %4.2f",
1076  k+startSec,detbest,gapf_offset[k]),
1077  1,Lmin,Lmax,1,gapf_min[k],gapf_max[k]);
1078  cGL->cd(k+1);
1079  htGLsec[k]->Draw();
1080  htSCsec[k] = new TH2D(Form("htSC%d",k),
1081  Form("Sector %d: #font[32]{sc} vs. %s for all sets, offset by %5.3f",
1082  k+startSec,detbest,sc_offset[k]),
1083  1,Lmin,Lmax,1,sc_min[k],sc_max[k]);
1084  cSC->cd(k+1);
1085  htSCsec[k]->Draw();
1086 
1087  fnGapf->SetParameter(4,fitParsSCGLsec[4+k]);
1088  fnSC->SetParameter(4,fitParsSCGLsec[4+k]);
1089  for (i=0; i<nfi; i++) {
1090  cGL->cd(k+1);
1091  SCi[i]->Draw(Form("gapf%d+%s*(%g)*%s+%g:%s",
1092  k+startSec,
1093  (glmode[i] == 2 ? "ugl" : Form("(%g)",ugl[i])),
1094  fitParsSCGLsec[0],uscvar,i*gapf_offset[k],detbest),
1095  gcut,"same");
1096  // should use fitParsSCGL[0] to match fnGapf and outlier bands,
1097  // but fitParsGL[0] to match fiGapf
1098  fiGapf->SetParameter(10,0+i*gapf_offset[k]);
1099  fiGapf->DrawCopy("same");
1100  fnGapf->SetParameter(10,0+i*gapf_offset[k]);
1101  fnGapf->SetLineColor(colorFit);
1102  fnGapf->DrawCopy("same");
1103  fnGapf->SetParameter(10,MAX_DEV*vgl+i*gapf_offset[k]);
1104  fnGapf->SetLineColor(colorOutlier);
1105  fnGapf->DrawCopy("same");
1106  fnGapf->SetParameter(10,-MAX_DEV*vgl+i*gapf_offset[k]);
1107  fnGapf->DrawCopy("same");
1108 
1109  cSC->cd(k+1);
1110  if (fitParsSCGL[6] == 1.0 && fitParsSCGL[7] == 1.0) {
1111  SCi[i]->Draw(Form("%s+%g:%s",scvar,i*sc_offset[k],detbest),gcut,"same");
1112  } else {
1113  SCi[i]->Draw(Form("%s-%s*(1-(%g+ugl)/(%g*(%g+ugl)))+%g:%s",scvar,uscvar,
1114  fitParsSCGLsec[1],fitParsSCGLsec[18],fitParsSCGLsec[17]*fitParsSCGLsec[1],
1115  i*sc_offset[k],detbest),gcut,"same");
1116  }
1117  fiSC->SetParameter(8,ugl[i]);
1118  fiSC->SetParameter(10,0+i*sc_offset[k]);
1119  fiSC->DrawCopy("same");
1120  fnSC->SetParameter(8,ugl[i]);
1121  fnSC->SetParameter(10,0+i*sc_offset[k]);
1122  fnSC->SetLineColor(colorFit);
1123  fnSC->DrawCopy("same");
1124  fnSC->SetParameter(10,MAX_DEV*vsc[jmin]+i*sc_offset[k]);
1125  fnSC->SetLineColor(colorOutlier);
1126  fnSC->DrawCopy("same");
1127  fnSC->SetParameter(10,-MAX_DEV*vsc[jmin]+i*sc_offset[k]);
1128  fnSC->DrawCopy("same");
1129  }
1130  }
1131  } else {
1132  if (!cGL) cGL = new TCanvas("cGL","GridLeak Fits",30,30,500,500);
1133  TH2D* htGL = new TH2D("htGL",
1134  Form("adjusted #font[32]{gapf} vs. %s for all sets, offset by %4.2f",
1135  detbest,gapf_offset[0]),
1136  1,Lmin,Lmax,1,gapf_min[0],gapf_max[0]);
1137  cGL->cd();
1138  htGL->Draw();
1139  if (!cSC) cSC = new TCanvas("cSC","SpaceCharge Fits",60,60,500,500);
1140  TH2D* htSC = new TH2D("htSC",
1141  Form("#font[32]{sc} vs. %s for all sets, offset by %5.3f",
1142  detbest,sc_offset[0]),
1143  1,Lmin,Lmax,1,sc_min[0],sc_max[0]);
1144  cSC->cd();
1145  htSC->Draw();
1146  for (i=0; i<nfi; i++) {
1147  cGL->cd();
1148  SCi[i]->Draw(Form("gapf+%s*(%g)*%s+%g:%s",
1149  (glmode[i] == 2 ? "ugl" : Form("(%g)",ugl[i])),
1150  fitParsSCGL[0],uscvar,i*gapf_offset[0],detbest),
1151  gcut,"same");
1152  // should use fitParsSCGL[0] to match fnGapf and outlier bands,
1153  // but fitParsGL[0] to match fiGapf
1154  fiGapf->SetParameter(10,0+i*gapf_offset[0]);
1155  fiGapf->DrawCopy("same");
1156  fnGapf->SetParameter(10,0+i*gapf_offset[0]);
1157  fnGapf->SetLineColor(colorFit);
1158  fnGapf->DrawCopy("same");
1159  fnGapf->SetParameter(10,MAX_DEV*vgl+i*gapf_offset[0]);
1160  fnGapf->SetLineColor(colorOutlier);
1161  fnGapf->DrawCopy("same");
1162  fnGapf->SetParameter(10,-MAX_DEV*vgl+i*gapf_offset[0]);
1163  fnGapf->DrawCopy("same");
1164 
1165  cSC->cd();
1166  if (fitParsSCGL[6] == 1.0 && fitParsSCGL[7] == 1.0) {
1167  SCi[i]->Draw(Form("%s+%g:%s",scvar,i*sc_offset[0],detbest),gcut,"same");
1168  } else {
1169  SCi[i]->Draw(Form("%s-%s*(1-(%g+ugl)/(%g*(%g+ugl)))+%g:%s",scvar,uscvar,
1170  fitParsSCGL[1],fitParsSCGL[7],fitParsSCGL[6]*fitParsSCGL[1],
1171  i*sc_offset[0],detbest),gcut,"same");
1172  }
1173  fiSC->SetParameter(8,ugl[i]);
1174  fiSC->SetParameter(10,0+i*sc_offset[0]);
1175  fiSC->DrawCopy("same");
1176  fnSC->SetParameter(8,ugl[i]);
1177  fnSC->SetParameter(10,0+i*sc_offset[0]);
1178  fnSC->SetLineColor(colorFit);
1179  fnSC->DrawCopy("same");
1180  fnSC->SetParameter(10,MAX_DEV*vsc[jmin]+i*sc_offset[0]);
1181  fnSC->SetLineColor(colorOutlier);
1182  fnSC->DrawCopy("same");
1183  fnSC->SetParameter(10,-MAX_DEV*vsc[jmin]+i*sc_offset[0]);
1184  fnSC->DrawCopy("same");
1185  }
1186  }
1187 
1188 
1189 
1191  // For zero space charge run, do not do the full fits,
1192  // just report back the SpaceCharge quantities to try
1193  // for each of the GridLeaks used:
1194 
1195  if (allZeros) {
1196  printf("\n\n*** The following calibration values may not be trusted at this time... ***");
1197  printf("\n\n*** Try the following calibration values: ***\n");
1198  for (i=0; i<nfi; i++) {
1199  printf("sc = %6.4g * ((%s) - (%6.4g))",sce[jmin],detbest,sof[jmin]);
1200  printf(" with GL = %5.2f\n\n",ugl[i]); // glmode[i]==2 ?
1201  }
1202  return 0;
1203  }
1204 
1205 
1206  if (cSummary==0) cSummary = new TCanvas("cSummary","Calib SC and GL");
1207  cSummary->Divide(2,2,0.01,0.025);
1208 
1209  Double_t min,max,ymin,ymax,ymin2,ymax2;
1210  Double_t xx[2];
1211 
1213  // Plot for gapf/(L-GLO) vs. usc*GLu/(L-GLO) (~SC*GL) => g1/g2 @ zero crossing
1214 
1215  TF1* miGL = new TF1("igapfL_of_SCGL","[0]*([1]-x)",0,5e-5);
1216  miGL->SetParameters(fitParsGL);
1217  miGL->SetLineWidth(1);
1218  miGL->SetLineColor(colorMiFit);
1219  TF1* myGL = new TF1("gapfL_of_SCGL","[0]*([2]*[4]-x)",0,5e-5);
1220  myGL->SetParameters(fitParsSCGL);
1221  myGL->SetLineWidth(1);
1222  myGL->SetLineColor(colorFit);
1223  // average used SC*GL converted from used lum scaler to best
1224  memset(asg,0,nfip*sizeof(Double_t));
1225  memset(glk,0,nfip*sizeof(Double_t));
1226  memset(nno,0,nfip*sizeof(Double_t));
1227  k = 0;
1228  for (i=0;i<nMeasures;i++) { // parameter sets
1229  if (outliersGL[i]) continue;
1230  int ifi = m_set[i];
1231  // gapf = g5*(L-GLO) - g2*usc*GLu = g2*(SC*GL*(L-GLO) - usc*GLu)
1232  // SC*GL = usc*GLu/(L-GLO) @ gapf = 0
1233  xgl[k] = m_usc[i]*m_ugl[i]/(m_L[i]-fitParsSCGL[5]);
1234  ygl[k] = m_gapf[i]/(m_L[i]-fitParsSCGL[5]);
1235  asg[ifi] += xgl[k];
1236  glk[ifi] += ygl[k];
1237  nno[ifi]++;
1238  k++;
1239  }
1240  for (i=0;i<nfi;i++) {
1241  asg[i] /= nno[i];
1242  glk[i] /= nno[i];
1243  }
1244  SetMinMax(k,xgl,min,max,scXgl_final);
1245  SetMinMax(k,ygl,ymin,ymax,zero);
1246  miGL->SetRange(min,max);
1247  myGL->SetRange(min,max);
1248 
1249  cSummary->cd(1);
1250  myGL->Draw();
1251  histo = myGL->GetHistogram();
1252  histo->SetTitle("#font[32]{gapf}/#font[32]{L} vs. SC#timesGL");
1253  histo->SetYTitle("#font[32]{gapf} / (#font[32]{L}_{best} - GLO)");
1254  histo->SetXTitle("#font[32]{sc}_{used} GL_{used} / (#font[32]{L}_{best} - GLO)");
1255  histo->GetYaxis()->SetTitleOffset(1.25);
1256  histo->SetMinimum(ymin);
1257  histo->SetMaximum(ymax);
1258  miGL->Draw("same");
1259  myGL->Draw("same");
1260  mark2.SetMarkerColor(colorData);
1261  mark2.SetMarkerStyle(1);
1262  mark2.DrawPolyMarker(k,xgl,ygl);
1263  mark2.SetMarkerColor(colorAggregate);
1264  mark2.SetMarkerStyle(27);
1265  mark2.DrawPolyMarker(nfi,asg,glk);
1266 
1267  ellip.DrawEllipse(scXgl_final,zero,escXgl_final,0.004*(ymax-ymin),0,360,0);
1268  mark.DrawMarker(scXgl_final,zero);
1269 
1270 
1272  // Plot for gapf vs. L - usc*GL*g2/g1 => GLO @ zero crossing
1273 
1274  TF1* miGO = new TF1("igapf_of_GLO","[0]*[1]*(x-[2])",0,5e-6);
1275  miGO->SetParameters(fitParsGL);
1276  miGO->SetLineWidth(1);
1277  miGO->SetLineColor(colorMiFit);
1278  TF1* myGO = new TF1("gapf_of_GLO","[0]*[2]*[4]*(x-[5])",0,5e-6);
1279  myGO->SetParameters(fitParsSCGL);
1280  myGO->SetLineWidth(1);
1281  myGO->SetLineColor(colorFit);
1282  // average used GLO converted from used lum scaler to best
1283  memset(ago,0,nfip*sizeof(Double_t));
1284  memset(glo,0,nfip*sizeof(Double_t));
1285  k = 0;
1286  for (i=0;i<nMeasures;i++) { // parameter sets
1287  if (outliersGL[i]) continue;
1288  int ifi = m_set[i];
1289  // GLO = L - (usc*GLu)/(SC*GL) @ gapf = 0
1290  xgo[k] = m_L[i] - (m_usc[i]*m_ugl[i])/(fitParsSCGL[2]*fitParsSCGL[4]);
1291  ygo[k] = m_gapf[i];
1292  ago[ifi] += xgo[k];
1293  glo[ifi] += ygo[k];
1294  k++;
1295  }
1296  for (i=0;i<nfi;i++) {
1297  ago[i] /= nno[i];
1298  glo[i] /= nno[i];
1299  }
1300  SetMinMax(nfi,ago,min,max,GLO);
1301  SetMinMax(k,ygo,ymin,ymax,zero);
1302  miGO->SetRange(min,max);
1303  myGO->SetRange(min,max);
1304 
1305  cSummary->cd(3);
1306  myGO->Draw();
1307  histo = myGO->GetHistogram();
1308  histo->SetTitle(Form("#font[32]{gapf} vs. GLO%s",(FORCE_GLO_SO ? " #equiv SO" : "")));
1309  histo->SetYTitle("#font[32]{gapf}");
1310  histo->SetXTitle("#font[32]{L}_{best} - (#font[32]{sc}_{used} GL_{used} g_{2} / g_{GL})");
1311  histo->GetYaxis()->SetTitleOffset(1.25);
1312  histo->SetMinimum(ymin);
1313  histo->SetMaximum(ymax);
1314  miGO->Draw("same");
1315  myGO->Draw("same");
1316  mark2.SetMarkerColor(colorData);
1317  mark2.SetMarkerStyle(1);
1318  mark2.DrawPolyMarker(nMeasures,xgo,ygo);
1319  mark2.SetMarkerColor(colorAggregate);
1320  mark2.SetMarkerStyle(27);
1321  mark2.DrawPolyMarker(nfi,ago,glo);
1322 
1323  ellip.DrawEllipse(GLO,zero,eGLO,0.004*(ymax-ymin),0,360,0);
1324  mark.DrawMarker(GLO,zero);
1325 
1326 
1328  // Plots for SC and SO vs. GL
1329 
1330  TF1* miSO = new TF1("iSO_of_GL","[0]",0,100);
1331  TF1* miSC = new TF1("iSC_of_GL","[1]/([2]+[0])",0,100);
1332  miSO->SetParameters(fitParsSC);
1333  miSC->SetParameters(fitParsSC);
1334  miSC->SetParameter(0,fitParsSCGL[4]);
1335  miSO->SetLineWidth(1);
1336  miSC->SetLineWidth(1);
1337  miSO->SetLineColor(colorMiFit);
1338  miSC->SetLineColor(colorMiFit);
1339  TF1* mySO = new TF1("SO_of_GL","[3]",0,100);
1340  TF1* mySC = new TF1("SC_of_GL","[2]",0,100);
1341  mySO->SetParameters(fitParsSCGL);
1342  mySC->SetParameters(fitParsSCGL);
1343  mySO->SetLineWidth(1);
1344  mySC->SetLineWidth(1);
1345  mySO->SetLineColor(colorFit);
1346  mySC->SetLineColor(colorFit);
1347  memset(nno,0,nfip*sizeof(Double_t));
1348  memset(spo,0,nfip*sizeof(Double_t));
1349  memset(spc,0,nfip*sizeof(Double_t));
1350  k = 0;
1351  for (i=0;i<nMeasures;i++) { // parameter sets
1352  if (outliersSC[i]) continue;
1353  int ifi = m_set[i];
1354  xx[0] = m_L[i];
1355  xx[1] = m_ugl[i];
1356  xsc[k] = m_ugl[i];
1357  // Basically, sc = SC*(L-SO), so...
1358  // SO = L-sc/SC
1359  // SC = sc/(L-SO)
1360  // More explicitly,
1361  // sc_real = [(sc_obs - sc_used) * g4r * (g5r * g5 + ugl) + sc_used * (g5 + ugl)] / (g5 + GL)
1362  temp1 = ((m_sc[i] - m_usc[i]) * fitParsSCGL[7] * (fitParsSCGL[6] * fitParsSCGL[1] + m_ugl[i]) +
1363  m_usc[i] * ( fitParsSCGL[1] + m_ugl[i])) /
1364  ( fitParsSCGL[1] + fitParsSCGL[4]);
1365  yso[k] = m_L[i] - (temp1 / fitParsSCGL[2]);
1366  ysc[k] = temp1 / (m_L[i] - fitParsSCGL[3]);
1367  spo[ifi] += yso[k];
1368  spc[ifi] += ysc[k];
1369  nno[ifi]++;
1370  k++;
1371  }
1372  j = 0;
1373  for (i=0;i<nfi;i++) {
1374  spo[i] /= nno[i];
1375  spc[i] /= nno[i];
1376  if (spc[i] - miSC->Eval(ugl[i]) > 0) j++;
1377  }
1378  i = (nfi>5 ? 1 : 0);
1379  if (j<=i || j>=nfi-i) {
1380  // too many spc[] points above or below the curve
1381  printf("WARNING! Symptoms of non-linearity in chosen luminosity scaler!\n");
1382  printf(" Suggestion: use ");
1383  if (DO_PCA) printf("more PCA variables\n\n");
1384  else printf("PCA method\n\n");
1385  }
1386  SetMinMax(nfi,ugl,min,max,glp);
1387  SetMinMax(k,yso,ymin,ymax,sop);
1388  SetMinMax(k,ysc,ymin2,ymax2,scp);
1389  miSO->SetRange(min,max);
1390  miSC->SetRange(min,max);
1391  mySO->SetRange(min,max);
1392  mySC->SetRange(min,max);
1393 
1394  cSummary->cd(2);
1395  mySC->Draw();
1396  histo = mySC->GetHistogram();
1397  histo->SetTitle("SC vs. GL");
1398  histo->SetYTitle("[#font[32]{sc} / (#font[32]{L}_{best} - SO)] (g_{5} + GL_{used}) / (g_{5} + GL)");
1399  histo->SetXTitle("GL_{used}");
1400  histo->GetYaxis()->SetTitleOffset(1.25);
1401  histo->SetMinimum(ymin2);
1402  histo->SetMaximum(ymax2);
1403  miSC->Draw("same");
1404  mySC->Draw("same");
1405  mark2.SetMarkerColor(colorData);
1406  mark2.SetMarkerStyle(1);
1407  mark2.DrawPolyMarker(k,xsc,ysc);
1408  mark2.SetMarkerColor(colorAggregate);
1409  mark2.SetMarkerStyle(27);
1410  mark2.DrawPolyMarker(nfi,ugl,spc);
1411 
1412 
1413  if (!scgl_fit) scgl_fit = new TF1("scgl_fit","[0]/x",-5.,100.);
1414  scgl_fit->SetParameter(0,scXgl_final);
1415  scgl_fit->SetLineColor(colorConstraint);
1416  scgl_fit->SetLineWidth(1);
1417  scgl_fit->DrawCopy("same");
1418  scgl_fit->SetLineStyle(7);
1419  scgl_fit->SetParameter(0,scXgl_final+escXgl_final);
1420  scgl_fit->DrawCopy("same");
1421  scgl_fit->SetParameter(0,scXgl_final-escXgl_final);
1422  scgl_fit->DrawCopy("same");
1423 
1424  ellip.DrawEllipse(glp,scp,eglp,escp,0,360,0);
1425  mark.DrawMarker(glp,scp);
1426 
1427  cSummary->cd(4);
1428  mySO->Draw();
1429  histo = mySO->GetHistogram();
1430  histo->SetTitle("SO vs. GL");
1431  histo->SetYTitle("#font[32]{L}_{best} - [#font[32]{sc} (g_{5} + GL_{used}) / (g_{5} + GL)]");
1432  histo->SetXTitle("GL_{used}");
1433  histo->GetYaxis()->SetTitleOffset(1.25);
1434  histo->SetMinimum(ymin);
1435  histo->SetMaximum(ymax);
1436  miSO->Draw("same");
1437  mySO->Draw("same");
1438  mark2.SetMarkerColor(colorData);
1439  mark2.SetMarkerStyle(1);
1440  mark2.DrawPolyMarker(k,xsc,yso);
1441  mark2.SetMarkerColor(colorAggregate);
1442  mark2.SetMarkerStyle(27);
1443  mark2.DrawPolyMarker(nfi,ugl,spo);
1444 
1445  ellip.DrawEllipse(glp,sop,eglp,esop,0,360,0);
1446  mark.DrawMarker(glp,sop);
1447 
1449  // Done
1450 
1451 
1452  // Set colors/markers for any further fun
1453  for (i=0;i<nfi;i++) {
1454  int color = 60 + (int) (40.0*((float) i)/((float) (nfi-1)));
1455  SCi[i]->SetMarkerColor(color);
1456  SCi[i]->SetLineColor(color);
1457  SCi[i]->SetMarkerStyle(22+i);
1458  }
1459 
1460 
1461  return 0;
1462 }
1463 
1465 // Fitting functions
1467 
1468 Double_t funcGapf(Double_t* x, Double_t* pars) {
1469  // x[0] = Luminosity to use
1470  // x[1] = usc*ugl
1471  // x[2] = const1
1472  // pars[0] = g2
1473  // pars[1] = g1 / g2 = SC * GL
1474  // pars[2] = GLO
1475  return (x[2] / REF_CONST1) * pars[0] * ( pars[1] * (x[0] - pars[2]) - x[1] );
1476 }
1477 
1478 Double_t funcSC(Double_t* x, Double_t* pars) {
1479  // x[0] = Luminosity
1480  // x[1] = ugl
1481  // pars[0] = SO
1482  // pars[1] = SCe = SC * (g5 + GL)
1483  // pars[2] = g5
1484  return pars[1] * (x[0] - pars[0]) / (x[1] + pars[2]);
1485 }
1486 
1487 Double_t funcSC2(Double_t* x, Double_t* pars) {
1488  // x[0] = Luminosity
1489  // x[1] = ugl
1490  // x[2] = usc
1491  // pars[0] = SO
1492  // pars[1] = SCe = SC * (g5 + GL)
1493  // pars[2] = g5
1494  // pars[3] = g5r = g5' / g5
1495  // pars[4] = g4r = g4' / g4
1496  return x[2] +
1497  (pars[1] * (x[0] - pars[0]) - (pars[2] + x[1]) * x[2]) /
1498  (pars[4] * (pars[3] * pars[2] + x[1]));
1499 }
1500 
1501 void fnchGapf(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
1502  //calculate chisquare
1503  gin = 0; iflag = 0;
1504  double chisq = 0;
1505  Double_t x[3];
1506  nMeasuresI = 0;
1507  for (int i=0;i<nMeasures;i++) {
1508  if (outliers[i]) continue;
1509  if (m_gapf[i] == 0) { outliers[i]=kTRUE; devs[i] = 0; continue; }
1510  x[0] = m_L[i];
1511  x[1] = m_usc[i] * m_ugl[i];
1512  x[2] = m_c1[i];
1513  devs[i] = m_gapf[i] - funcGapf(x,par);
1514  maxvarGL[i] = devs[i]*devs[i];
1515  int k = m_runIdx[i];
1516  if (k!=i) {
1517  if (maxvarGL[i]>maxvarGL[k]) {
1518  for (int j=k;j<i;j++) {
1519  if (m_runIdx[j]==k) maxvarGL[j] = maxvarGL[i];
1520  }
1521  } else maxvarGL[i] = maxvarGL[k];
1522  }
1523  }
1524  for (int i=0;i<nMeasures;i++) {
1525  if (outliers[i]) continue;
1526  chisq += RUN_CORRELATE*maxvarGL[i] + (1.0-RUN_CORRELATE)*devs[i]*devs[i];
1527  nMeasuresI++;
1528  }
1529  CHISQ = (chisq/(STDDEV*STDDEV))/(nMeasuresI-npar); // chisq/dof
1530  f = CHISQ;
1531 }
1532 
1533 void fnchSC(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
1534  //calculate chisquare
1535  gin = 0; iflag = 0;
1536  double chisq = 0;
1537  Double_t x[2];
1538  Double_t pars[3];
1539  pars[0] = par[0];
1540  pars[1] = TMath::Exp(par[1]);
1541  pars[2] = TMath::Exp(par[2]);
1542  nMeasuresI = 0;
1543  for (int i=0;i<nMeasures;i++) {
1544  if (outliers[i]) continue;
1545  if (m_sc[i] == 0) { outliers[i]=kTRUE; devs[i] = 0; continue; }
1546  x[0] = m_L[i];
1547  x[1] = m_ugl[i];
1548  devs[i] = m_sc[i] - funcSC(x,pars);
1549  maxvarSC[i] = devs[i]*devs[i];
1550  int k = m_runIdx[i];
1551  if (k!=i) {
1552  if (maxvarSC[i]>maxvarSC[k]) {
1553  for (int j=k;j<i;j++) {
1554  if (m_runIdx[j]==k) maxvarSC[j] = maxvarSC[i];
1555  }
1556  } else maxvarSC[i] = maxvarSC[k];
1557  }
1558  }
1559  for (int i=0;i<nMeasures;i++) {
1560  if (outliers[i]) continue;
1561  chisq += RUN_CORRELATE*maxvarSC[i] + (1.0-RUN_CORRELATE)*devs[i]*devs[i];
1562  nMeasuresI++;
1563  }
1564  CHISQ = (chisq/(STDDEV*STDDEV))/(nMeasuresI-npar); // chisq/dof
1565  f = CHISQ;
1566 }
1567 
1568 void fnchSCGapf(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
1569  //calculate chisquare
1570  gin = 0; iflag = 0;
1571  double chisqGL = 0;
1572  double chisqSC = 0;
1573  Double_t x[3];
1574  nMeasuresI = 0;
1575  Double_t parGL[3];
1576  Double_t parSC[5];
1577  Double_t SC = TMath::Exp(par[2]);
1578  Double_t GL = TMath::Exp(par[4]);
1579  Double_t g5 = TMath::Exp(par[1]);
1580  Double_t g5r= TMath::Exp(par[6]); // g5r = g5' / g5
1581  Double_t g4r= TMath::Exp(par[7]); // g4r = g4' / g4
1582  Double_t ewr= TMath::Exp(par[8]);
1583  parGL[0] = par[0]; // g2
1584  parGL[1] = SC*GL; // g1 / g2 = SC * GL
1585  parGL[2] = (FORCE_GLO_SO ? par[3] : par[5]); // GLO
1586  parSC[0] = par[3]; // SO
1587  parSC[1] = SC*(g5 + GL); // SCe = SC * (g5 + GL)
1588  parSC[2] = g5; // g5
1589  parSC[3] = (FORCE_g5_same ? 1. : (FORCE_g3_same ? 1./g4r : g5r));
1590  parSC[4] = (FORCE_g4_same ? 1. : g4r);
1591  for (int l=0;l<(EW_ASYMMETRY ? 2 : 1);l++) {
1592  if (l) { // east-side
1593  parGL[1] *= ewr;
1594  parSC[1] *= ewr;
1595  }
1596 
1597  for (int i=0;i<nMeasures;i++) {
1598  int m = i + l*nMeasuresHalf;
1599  x[0] = m_L[i];
1600  x[1] = m_usc[m] * m_ugl[i];
1601  x[2] = m_c1[i];
1602  if (!outliersGL[i]) {
1603  devsGL[i] = m_gapf[m] - funcGapf(x,parGL);
1604  maxvarGL[i] = devsGL[i]*devsGL[i];
1605  int k = m_runIdx[i];
1606  if (k!=i) {
1607  if (maxvarGL[i]>maxvarGL[k]) {
1608  for (int j=k;j<i;j++) {
1609  if (m_runIdx[j]==k) maxvarGL[j] = maxvarGL[i];
1610  }
1611  } else maxvarGL[i] = maxvarGL[k];
1612  }
1613  }
1614 
1615  x[1] = m_ugl[i];
1616  x[2] = m_usc[m];
1617  if (!outliersSC[i]) {
1618  devsSC[i] = m_sc[m] - funcSC2(x,parSC);
1619  maxvarSC[i] = devsSC[i]*devsSC[i];
1620  int k = m_runIdx[i];
1621  if (k!=i) {
1622  if (maxvarSC[i]>maxvarSC[k]) {
1623  for (int j=k;j<i;j++) {
1624  if (m_runIdx[j]==k) maxvarSC[j] = maxvarSC[i];
1625  }
1626  } else maxvarSC[i] = maxvarSC[k];
1627  }
1628  }
1629  }
1630  for (int i=0;i<nMeasures;i++) {
1631  if (!outliersGL[i]) {
1632  chisqGL += RUN_CORRELATE*maxvarGL[i] + (1.0-RUN_CORRELATE)*devsGL[i]*devsGL[i];
1633  nMeasuresI++;
1634  }
1635  if (!outliersSC[i]) {
1636  chisqSC += RUN_CORRELATE*maxvarSC[i] + (1.0-RUN_CORRELATE)*devsSC[i]*devsSC[i];
1637  nMeasuresI++;
1638  }
1639  }
1640  }
1641  CHISQ = ((chisqGL/(STDDEV_GL*STDDEV_GL))+(chisqSC/(STDDEV_SC*STDDEV_SC)))/(nMeasuresI-npar); // chisq/dof
1642  f = CHISQ;
1643 }
1644 
1645 void fnchSCGapfSec(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
1646  //calculate chisquare
1647  gin = 0; iflag = 0;
1648  double chisqGL = 0;
1649  double chisqSC = 0;
1650  Double_t x[3];
1651  nMeasuresI = 0;
1652  Double_t parGL[3];
1653  Double_t parSC[5];
1654  Double_t SC = TMath::Exp(par[2]);
1655  Double_t g5 = TMath::Exp(par[1]);
1656  Double_t g5r= TMath::Exp(par[17]); // g5r = g5' / g5
1657  Double_t g4r= TMath::Exp(par[18]); // g4r = g4' / g4
1658  parGL[0] = par[0]; // g2
1659  parGL[2] = (FORCE_GLO_SO ? par[3] : par[16]); // GLO
1660  parSC[0] = par[3]; // SO
1661  parSC[2] = g5; // g5
1662  parSC[3] = (FORCE_g5_same ? 1. : (FORCE_g3_same ? 1./g4r : g5r));
1663  parSC[4] = (FORCE_g4_same ? 1. : g4r);
1664  for (int sec = 0; sec<12; sec++) {
1665  if (unfittableSec[sec]) continue;
1666 
1667  Double_t GL = TMath::Exp(par[4+sec]);
1668  parGL[1] = SC*GL; // g1 / g2 = SC * GL
1669  parSC[1] = SC*(g5 + GL); // SCe = SC * (g5 + GL)
1670 
1671  for (int i=0;i<nMeasures;i++) {
1672  x[0] = m_L[i];
1673  x[1] = m_usc[i] * m_ugl[i];
1674  x[2] = m_c1[i];
1675  if (!outliersGL[i] && m_gapfS[sec][i]!=0.0) {
1676  devsGLsec[sec][i] = m_gapfS[sec][i] - funcGapf(x,parGL);
1677  maxvarGL[i] = devsGLsec[sec][i]*devsGLsec[sec][i];
1678  int k = m_runIdx[i];
1679  if (k!=i) {
1680  if (maxvarGL[i]>maxvarGL[k]) {
1681  for (int j=k;j<i;j++) {
1682  if (m_runIdx[j]==k) maxvarGL[j] = maxvarGL[i];
1683  }
1684  } else maxvarGL[i] = maxvarGL[k];
1685  }
1686  }
1687 
1688  x[1] = m_ugl[i];
1689  x[2] = m_usc[i];
1690  if (!outliersSC[i] && m_scS[sec][i]!=0.0) {
1691  devsSCsec[sec][i] = m_scS[sec][i] - funcSC2(x,parSC);
1692  maxvarSC[i] = devsSCsec[sec][i]*devsSCsec[sec][i];
1693  int k = m_runIdx[i];
1694  if (k!=i) {
1695  if (maxvarSC[i]>maxvarSC[k]) {
1696  for (int j=k;j<i;j++) {
1697  if (m_runIdx[j]==k) maxvarSC[j] = maxvarSC[i];
1698  }
1699  } else maxvarSC[i] = maxvarSC[k];
1700  }
1701  }
1702  }
1703  for (int i=0;i<nMeasures;i++) {
1704  if (!outliersGL[i] && m_gapfS[sec][i]!=0.0) {
1705  chisqGL += RUN_CORRELATE*maxvarGL[i] + (1.0-RUN_CORRELATE)*devsGLsec[sec][i]*devsGLsec[sec][i];
1706  nMeasuresI++;
1707  }
1708  if (!outliersSC[i] && m_scS[sec][i]!=0.0) {
1709  chisqSC += RUN_CORRELATE*maxvarSC[i] + (1.0-RUN_CORRELATE)*devsSCsec[sec][i]*devsSCsec[sec][i];
1710  nMeasuresI++;
1711  }
1712  }
1713  }
1714  CHISQ = ((chisqGL/(STDDEV_GL*STDDEV_GL))+(chisqSC/(STDDEV_SC*STDDEV_SC)))/(nMeasuresI-npar); // chisq/dof
1715  f = CHISQ;
1716 }
1717 
1719 // Helper functions
1721 
1722 
1723 int Init(const char* input) {
1724  if (SCall) return 0;
1725 
1726  gStyle->SetGridColor(colorGrid);
1727  gStyle->SetOptStat(0);
1728  gStyle->SetPalette(1);
1729 
1730  int i,k;
1731  char fname[256];
1732  double temp;
1733  TString fis[nfip];// file names
1734  SCall = new TChain("SC","ChainAll");
1735 
1736  if (!input || input[0] == 0 || input[0] == '@') {
1737 
1738  nfi = 0;
1739  int nfii[nfip];
1740  memset(nfii,0,nfip*sizeof(int));
1741  TSeqCollection* listOfFiles = 0;
1742  TObjArray SCGLcombos;
1743  TFile* fileI = 0;
1744  if (!input) { // Use already opened files
1745  listOfFiles = gROOT->GetListOfFiles();
1746  } else { // Unified file specification
1747  fis[0] = input;
1748  fis[0].Remove(0,1);
1749  SCall->Add(fis[0].Data());
1750  listOfFiles = SCall->GetListOfFiles();
1751  }
1752  Int_t nfiles = (listOfFiles ? listOfFiles->GetEntries() : 0);
1753  if (nfiles == 0) {
1754  printf("Error: no input specified!\n");
1755  return 1;
1756  }
1757  for (i=nfiles-1;i>=0;i--) {
1758  if (!input) {
1759  fileI = (TFile*) (listOfFiles->At(i));
1760  SCall->Add(fileI->GetName());
1761  } else {
1762  TChainElement* elem = (TChainElement*) (listOfFiles->At(i));
1763  if (elem->GetEntries() == 0) continue;
1764  fileI = new TFile(elem->GetTitle());
1765  }
1766  TKey* corrKey = fileI->GetKey("SCcorrection");
1767  if (!corrKey) {
1768  printf("File is missing its corrections (excluding from grouping): %s\n",fileI->GetName());
1769  continue;
1770  } else {
1771  TString corrStr = corrKey->ReadObj()->GetTitle();
1772  corrKey = fileI->GetFile()->GetKey("GLcorrection");
1773  (corrStr += ":") += corrKey->ReadObj()->GetTitle();
1774  corrKey = fileI->GetFile()->GetKey("SCEWRatio");
1775  if (corrKey)
1776  (corrStr += ":") += corrKey->ReadObj()->GetTitle();
1777  TObject* existing = SCGLcombos.FindObject(corrStr.Data());
1778  if (existing) {
1779  k = SCGLcombos.IndexOf(existing);
1780  int kent = SCi[k]->GetEntries();
1781  SCi[k]->Add(fileI->GetName());
1782  if (SCi[k]->GetEntries() == kent)
1783  printf("Skipping file with 0 ntuple entries: %s\n",fileI->GetName());
1784  else nfii[k]++;
1785  } else {
1786  SCi[nfi] = new TChain("SC",Form("Chain%d : %s",nfi,corrStr.Data()));
1787  SCi[nfi]->Add(fileI->GetName());
1788  if (SCi[nfi]->GetEntries()>0) {
1789  SCGLcombos.AddLast(new TNamed(corrStr.Data(),corrStr.Data()));
1790  nfii[nfi] = 1;
1791  nfi++;
1792  if (nfi>nfip) {
1793  printf("Warning! Exceeding maximum number of sets! Excluding the rest from grouping.\n");
1794  break;
1795  }
1796  } else {
1797  delete SCi[nfi];
1798  printf("Skipping file with 0 ntuple entries: %s\n",fileI->GetName());
1799  }
1800  }
1801  }
1802  fileI->Close();
1803  }
1804  for (i=0;i<nfi;i++) {
1805  printf("Set %3d: %d files added with...\n",i,nfii[i]);
1806  TString corrStr = SCGLcombos.At(i)->GetName();
1807  corrStr.Remove(corrStr.First(':'));
1808  if (corrStr.Length()) allZeros = kFALSE;
1809  printf(" used sc = %s\n",corrStr.Data());
1810  corrStr = SCGLcombos.At(i)->GetName();
1811  corrStr.Remove(0,corrStr.First(':')+1);
1812  if (corrStr.CountChar(':')) {
1813  corrStr.Remove(0,corrStr.First(':')+1);
1814  printf(" used ewratio = %s\n",corrStr.Data());
1815  corrStr = SCGLcombos.At(i)->GetName();
1816  corrStr.Remove(0,corrStr.First(':')+1);
1817  corrStr.Remove(corrStr.First(':'));
1818  }
1819  ugl[i] = corrStr.Atof();
1820  glmode[i] = 1;
1821  printf(" used GL = %g\n",ugl[i]);
1822  }
1823  } else {
1824  // Read input file
1825  nfi = nfip;
1826  ifstream dats(input);
1827  if (!(dats.good())) {
1828  printf("Error: problems opening input file specified!\n");
1829  return 2;
1830  }
1831  for (i=0;i<nfi && dats.good();i++) {
1832  dats >> fname;
1833  if (dats.eof()) { nfi=i; continue; }
1834  fis[i] = fname;
1835  if (fis[i].Contains('@')) {
1836  fis[i].Remove(fis[i].First('@'),1);
1837  } else { // depracated; supported only to read in ugl
1838  dats >> temp >> temp >> ugl[i];
1839  if (ugl[i]) allZeros = kFALSE;
1840  if (fis[i].Contains('!')) { // Defunct
1841  dats >> temp;
1842  fis[i].Remove(fis[i].First('!'),1);
1843  }
1844  }
1845 
1846  // Skip lines beginning with #
1847  if (fname[0]=='#') { i--; continue; }
1848 
1849  // Build TChains from histogram files
1850  SCi[i] = new TChain("SC",Form("Chain%d : %s",i,fis[i].Data()));
1851  int added_files = 0;
1852  if ((added_files = SCi[i]->Add(fis[i].Data())) < 1) {
1853  printf("Warning: no files added from %s (Set %d)\n",fis[i].Data(),i);
1854  } else {
1855  printf("Set %3d: %d files added from %s\n",i,added_files,fis[i].Data());
1856  TKey* corrKey = SCi[i]->GetFile()->GetKey("SCcorrection");
1857  if (corrKey) {
1858  TString corrStr = corrKey->ReadObj()->GetTitle();
1859  if (corrStr.Length()) allZeros = kFALSE;
1860  printf(" used sc = %s\n",corrStr.Data());
1861  corrKey = SCi[i]->GetFile()->GetKey("SCEWRatio");
1862  if (corrKey) {
1863  corrStr = corrKey->ReadObj()->GetTitle();
1864  printf(" used ewratio = %s\n",corrStr.Data());
1865  }
1866  corrKey = SCi[i]->GetFile()->GetKey("GLcorrection");
1867  corrStr = corrKey->ReadObj()->GetTitle();
1868  ugl[i] = corrStr.Atof();
1869  }
1870  printf(" used GL = %g\n",ugl[i]);
1871  // glmode: from ntuple (2), from file (1), from input.dat (0)
1872  glmode[i] = (SCi[i]->GetBranch("ugl") ? 2 :
1873  (corrKey ? 1 : 0));
1874  SCall->Add(fis[i].Data());
1875  }
1876  }
1877  }
1878  printf("Found %d dataset specifications.\n",nfi);
1879 
1880  for (i=0;i<nfi;i++) {
1881  if (BY_SECTOR && !(SCi[i]->GetBranch("sc1"))) {
1882  printf("ERROR: BY_SECTOR specified, but not possible with Set %d!\n",i);
1883  return 3;
1884  }
1885  SCi[i]->SetMarkerStyle(7);
1886  SCi[i]->SetMarkerColor(colorData);
1887  }
1888  SCall->SetMarkerStyle(7);
1889 
1890  memset(conSC,0,nsca*sizeof(TCanvas*));
1891 
1892  // Minimization items
1893  for (k=0;k<3;k++) parNameSC[k] = new char[16];
1894  for (k=0;k<3;k++) parNameGL[k] = new char[16];
1895  for (k=0;k<9;k++) parNameSCGL[k] = new char[16];
1896  for (k=0;k<19;k++) parNameSCGLsec[k] = new char[16];
1897  arglist[0] = 5000;
1898 
1899  return 0;
1900 }
1901 
1902 int Waiting() {
1903  int tempi;
1904  printf("Waiting...\n");
1905  cin >> tempi;
1906  return tempi;
1907 }
1908 
1909 void SetMinMax(int n, Double_t* ar, double& min, double& max,
1910  double inclusion, double margin) {
1911  min = 1e10;
1912  max = -1e10;
1913  for (int i=0; i<n; i++) {
1914  if (min > ar[i]) min = ar[i];
1915  if (max < ar[i]) max = ar[i];
1916  }
1917  if (min > inclusion) min = inclusion;
1918  else if (max < inclusion) max = inclusion;
1919  double margins = margin*(max-min);
1920  min -= margins;
1921  max += margins;
1922 }
1923 
1924 int SetMinuitPars(const int n, TString* names, Double_t* starts, Double_t* steps, int debug, Bool_t* fixing) {
1925  TVirtualFitter::SetDefaultFitter("Minuit");
1926  TVirtualFitter::SetMaxIterations(500000000);
1927  minuit = TVirtualFitter::Fitter(0,26);
1928  minuit2 = ((TFitter*) minuit)->GetMinuit();
1929  minuit2->SetPrintLevel(debug - 2);
1930 
1931  int i, maxlen = 0;
1932  for (i=0; i<n; i++) if (names[i].Length() > maxlen) maxlen = names[i].Length();
1933  for (i=0; i<n; i++) {
1934  names[i].Append(' ',maxlen-names[i].Length());
1935  if (fixing && fixing[i]) {
1936  minuit->SetParameter(i, names[i].Data(), starts[i], 0, 0, 0);
1937  minuit->FixParameter(i);
1938  } else {
1939  // Rules for limits (default is no limits):
1940  double lower = 0.;
1941  double upper = 0.;
1942  if (names[i].BeginsWith("log(g5")) {
1943  lower = starts[i] - 1.0;
1944  upper = starts[i] + 1.0;
1945  } else if (names[i].BeginsWith("log")) {
1946  lower = starts[i] - 3.0;
1947  upper = starts[i] + 3.0;
1948  } else if (names[i].Contains("g2")) {
1949  lower = starts[i] * 1e-1;
1950  upper = starts[i] * 1e1;
1951  }
1952  minuit->SetParameter(i, names[i].Data(), starts[i], steps[i], lower, upper);
1953  }
1954  }
1955  return maxlen;
1956 }
1957 
1958 void GetMinuitPar(const int n) {
1959  // Get a single parameter values from Minuit and ignore the errors and limits
1960  double temp1,temp2,temp3;
1961  minuit->GetParameter(n,parName[n],fitPars[n],temp1,temp2,temp3);
1962 }
1963 
1964 void GetMinuitPars(const int n, Bool_t* fixing, Double_t* fixedValues) {
1965  // Get all parameter values and errors from Minuit, unless they are a fixed value,
1966  // and record if parameters are at or near their limits in atLimit array
1967  double lowerLimit,upperLimit;
1968  int i;
1969  for (i=0; i<n; i++) {
1970  if (fixing && fixing[i]) {
1971  fitPars[i] = fixedValues[i];
1972  fitParErrs[i] = 0;
1973  atLimit[i] = false;
1974  } else {
1975  minuit->GetParameter(i,parName[i],fitPars[i],fitParErrs[i],lowerLimit,upperLimit);
1976  if (lowerLimit==0 && upperLimit==0) { atLimit[i] = false; continue; }
1977  double limitsCenter = 0.5*(lowerLimit+upperLimit);
1978  double limitsWindow = 0.49 * (upperLimit-lowerLimit);
1979  atLimit[i] = (TMath::Abs(fitPars[i]-limitsCenter) > limitsWindow);
1980  }
1981  }
1982 }
1983 
1984 void Log2Lin(int i) {
1985  fitPars[i] = TMath::Exp(fitPars[i]);
1986  fitParErrs[i] *= fitPars[i];
1987 }
1988 
1989 int FitWithOutlierRemoval(int npar, int debug) {
1990  // Iterates the fit, until a stable set of
1991  // outliers at > MAX_DEV*VAR are removed
1992  int i,k,status,nOutliers;
1993  double temp1,temp2;
1994  Double_t VAR_sets,VAR_single;
1995  Double_t nSet[nfip];
1996  Double_t VAR_singles[nfip];
1997  int iter = 0;
1998  int nBytes = nMeasuresMax*sizeof(Bool_t);
1999  int nBytes2 = nfip*sizeof(Double_t);
2000  memset(outliers,0,nBytes);
2001  Bool_t outliers1[nMeasuresMax];
2002  Bool_t outliers2[nMeasuresMax];
2003  memset(devs,0,nMeasuresMax*sizeof(Double_t));
2004 
2005  while (1) {
2006  if (iter>1) memcpy(outliers2,outliers1,nBytes);
2007  if (iter>0) memcpy(outliers1,outliers,nBytes);
2008  status = minuit->ExecuteCommand("MINIMIZE", arglist, 1);
2009  if (status) return status;
2010 
2011  // One more call to the fnch to get the deviations correct
2012  GetMinuitPars(npar);
2013  (*(minuit->GetFCN()))(k,&temp1,temp2,fitPars,status);
2014 
2015  // Calculate apprpriate errors
2016  VAR = 0;
2017  VAR_sets = 0;
2018  VAR_single = 0;
2019  memset(devs_set,0,nBytes2);
2020  memset(nSet,0,nBytes2);
2021  memset(VAR_singles,0,nBytes2);
2022  for (i=0; i<nMeasures; i++) {
2023  if (outliers[i]) continue;
2024  VAR += devs[i]*devs[i];
2025  int ifi = m_set[i];
2026  nSet[ifi]++;
2027  devs_set[ifi] += devs[i];
2028  }
2029  for (i=0; i<nfi; i++) {
2030  devs_set[i] /= nSet[i];
2031  VAR_sets += devs_set[i]*devs_set[i];
2032  }
2033  for (i=0; i<nMeasures; i++) {
2034  if (outliers[i]) continue;
2035  int ifi = m_set[i];
2036  VAR_singles[ifi] += (devs[i]-devs_set[ifi])*(devs[i]-devs_set[ifi]);
2037  }
2038  for (i=0; i<nfi; i++) {
2039  VAR_singles[i] = TMath::Sqrt(VAR_singles[i] / nSet[i]);
2040  VAR_single += VAR_singles[i];
2041  }
2042  // VAR_single is the average variance of the single sets:
2043  VAR_single /= ((Double_t) nfi);
2044  // VAR_set is the vaiance of one set to another
2045  VAR_sets = TMath::Sqrt(VAR_sets / ((Double_t) nfi-1));
2046  // VAR is the combined observed variance of all data
2047  VAR = TMath::Sqrt(VAR/((Double_t) nMeasuresI-1));
2048 
2049  // Use an equivalent variance^2 = (nfi-2) * VAR_single^2 + (<n_per_set> -2) * VAR_sets^2
2050  double N_sets = nfi - 2.0;
2051  double N_sing = (((Double_t) nMeasuresI)/((Double_t) nfi)) - 2.0;
2052  if (N_sets <= 0 || N_sing <=0) {
2053  printf("WARNING: could not estimate point-to-point errors due to limited samples\n");
2054  } else {
2055  STDDEV = TMath::Sqrt(N_sets* VAR_single*VAR_single + N_sing * VAR_sets*VAR_sets);
2056  }
2057  if (USE_OLD_STDDEV) STDDEV = VAR;
2058 
2059  // Determine outliers from this pass and compare to previous passes
2060  nOutliers = 0;
2061  for (i=0; i<nMeasures; i++) {
2062  if (devs[i] == 0 || TMath::Abs(devs[i]) > MAX_DEV*VAR) { outliers[i] = kTRUE; nOutliers++; }
2063  else outliers[i] = kFALSE;
2064  }
2065  if (iter>0 && memcmp(outliers,outliers1,nBytes) == 0) break;
2066  if (iter>1 && memcmp(outliers,outliers2,nBytes) == 0) break;
2067  if (iter > nMeasures/10) {
2068  printf("\nWARNING: Outlier removal at 10 percent of data points...stopping here.\n");
2069  break;
2070  }
2071  iter++;
2072  }
2073  if (debug>0) printf("Finished outlier removal: %d (of %d) data points removed (%d iterations)\n",
2074  nOutliers,nMeasures,iter+1);
2075 
2076  // Constrain errors by fixing other parameters
2077  GetMinuitPars(npar);
2078  status = minuit->ExecuteCommand("SET LIM", 0, 0);
2079  if (status) return status;
2080  double eplus,eminus,eparab,gcc;
2081  for (k=0;k<npar;k++) {
2082  if (atLimit[k]) {
2083  printf("%s\t:\t%g\t+/- %g\t(AT LIMIT!)\n",parName[k],fitPars[k],fitParErrs[k]);
2084  continue;
2085  }
2086  for (i=0;i<npar;i++) if (i!=k) minuit->FixParameter(i);
2087  status = minuit->ExecuteCommand("MINIMIZE", arglist, 1);
2088  if (status) return status;
2089  arglist[1] = (double) k+1;
2090  status = minuit->ExecuteCommand("MINOS", arglist, 2);
2091  if (status) return status;
2092  GetMinuitPar(k);
2093  minuit2->mnerrs(k,eplus,eminus,eparab,gcc);
2094  if (debug>0) printf("%s\t:\t%g\t+%g/%g or +/- %g\n",parName[k],fitPars[k],eplus,eminus,eparab);
2095  if (eplus!=0.0 && eminus!=0.0) fitParErrs[k] = 0.5*(eplus-eminus);
2096  else if (eplus!=0.0 || eminus!=0.0) fitParErrs[k] = eparab;
2097  printf("%s\t:\t%g\t+/- %g\n",parName[k],fitPars[k],fitParErrs[k]);
2098  for (i=0;i<npar;i++) if (i!=k) minuit->ReleaseParameter(i);
2099  }
2100 
2101  // One more determination of variance using final fit,
2102  // with and without outlier removal
2103  memcpy(outliers1,outliers,nBytes);
2104  memset(outliers,0,nBytes);
2105  (*(minuit->GetFCN()))(k,&temp1,temp2,fitPars,status);
2106  memcpy(outliers,outliers1,nBytes);
2107  VAR = 0;
2108  VAR_all = 0;
2109  for (i=0; i<nMeasures; i++) {
2110  if (!outliers[i]) VAR += devs[i]*devs[i];
2111  VAR_all += devs[i]*devs[i];
2112  }
2113  VAR = TMath::Sqrt(VAR/((Double_t) nMeasuresI-1));
2114  VAR_all = TMath::Sqrt(VAR_all/((Double_t) nMeasures-1));
2115 
2116  return status;
2117 }
2118 
2119 void DrawErrorContours(int npar, Bool_t* fix) {
2120  int i,j,k,l,status = 0;
2121  Double_t* gin = 0;
2122  Double_t x0,y0,chi;
2123  Int_t indx = 0;
2124  Double_t lows[128];
2125  Double_t highs[128];
2126  TString pName[32];
2127  for (k=0;k<npar;k++) {
2128  if (fix && fix[k]) continue;
2129  lows [k] = fitPars[k]-2.5*fitParErrs[k];
2130  highs[k] = fitPars[k]+2.5*fitParErrs[k];
2131  pName[k] = parName[k];
2132  pName[k].Remove(TString::kTrailing,' ');
2133  }
2134  TVirtualPad* origPad = gPad;
2135  double plmin = (USE_OLD_STDDEV ? 0.75 : 1.0/nfi);
2136  double plmax = plmin*100.;
2137  for (i=0;i<npar;i++) {
2138  if (fix && fix[i]) continue;
2139  x0 = fitPars[i];
2140  for (j=0;j<i;j++) {
2141  if (fix && fix[j]) continue;
2142  if (indx==15) {
2143  printf("Stopping error contour drawing at limit (5 parameters)!\n");
2144  return;
2145  }
2146  y0 = fitPars[j];
2147  TH2D* plot = new TH2D(Form("h%s_%d_%d",origPad->GetName(),i,j),
2148  Form("%s vs. %s",pName[j].Data(),pName[i].Data()),
2149  40,lows[i],highs[i],40,lows[j],highs[j]);
2150  for (k=0;k<40;k++) {
2151  fitPars[i] = lows[i]+(((double) k)+0.5)*(highs[i]-lows[i])/40.;
2152  for (l=0;l<40;l++) {
2153  fitPars[j] = lows[j]+(((double) l)+0.5)*(highs[j]-lows[j])/40.;
2154  (*(minuit->GetFCN()))(npar, gin, chi, fitPars, status);
2155  plot->SetBinContent(k+1,l+1,chi);
2156  }
2157  }
2158  fitPars[i] = x0;
2159  fitPars[j] = y0;
2160  indx++;
2161  origPad->cd(indx)->SetLogz();
2162  plot->SetMinimum(plmin);
2163  plot->SetMaximum(plmax);
2164  plot->Draw("zcol");
2165  }
2166  }
2167 }
2168 
2169 TString PCA(int Nmax, int debug) {
2170  printf("Running principal components analysis...\n");
2171 
2172  TString scl[maxp];
2173  TString scb[maxp];
2174  Double_t rmsl[maxp];
2175  long long rmsi[maxp];
2176  TPrincipal* pp = 0;
2177 
2178  int map[nsca];
2179  Double_t coef[nsca];
2180  Double_t x[nsca];
2181  Double_t* xvals[nsca];
2182  Double_t p[nsca];
2183  TString delim=":";
2184  TString scas = (scastr.Length() ? scastr :
2185  "zdcx:zdcw:zdce:bbcx:bbcw:bbce:bbcyb:bbcbb:vpdx:vpde:vpdw:zdcxnk:zdcwnk:zdcenk");
2186  scas += ":sc";
2187  TObjArray* scaA = scas.Tokenize(delim);
2188  int i,n,N1 = scaA->GetEntries();
2189  TH1F rmsp("rmsp","rmsp",500,-0.01,0.01);
2190  int iter0fi = nfi-1;
2191 
2192  if (pcaN >= 0) for (i=0;i<maxp;i++) if (ppl[i]) delete ppl[i];
2193  memset(ppl,0,maxp*sizeof(TPrincipal*));
2194 
2195  // Read the file once and store the values in a dummy TPrincipal
2196  //ppl[0] = SCpca->Principal(scas.Data(),cut);
2197  ppl[0] = new TPrincipal(N1,"D");
2198  ppl[0]->SetName("principal0");
2199  for (n=0; n<N1; n++) {
2200  xvals[n] = new Double_t[nMeasuresMax];
2201  if (ITER0) {
2202  SCi[iter0fi]->Draw(scaA->At(n)->GetName(),cut,"goff");
2203  nMeasures = SCi[iter0fi]->GetSelectedRows();
2204  memcpy(xvals[n],SCi[iter0fi]->GetV1(),nMeasures*sizeof(Double_t));
2205  } else {
2206  nMeasures = 0;
2207  for (i=0; i<nfi; i++) {
2208  SCi[i]->Draw(Form("%s*%g",scaA->At(n)->GetName(),
2209  (n < N1-1 ? 1.0 : (fitPars[1]+ugl[i])/(fitPars[1]+fitPars[4]))),
2210  cut,"goff");
2211  /* Using fitPars[6] and fitPars[7] can be dangerous if first iteration fit wasn't very good
2212  if (n < N1-1)
2213  SCi[i]->Draw(Form("%s",scaA->At(n)->GetName()),cut,"goff");
2214  else if (fitPars[6] == 1.0 && fitPars[7] == 1.0)
2215  SCi[i]->Draw(Form("sc*%g",(fitPars[1]+ugl[i])/(fitPars[1]+fitPars[4])),cut,"goff");
2216  else
2217  SCi[i]->Draw(Form("(sc-usc)*(%g)+usc*(%g)",
2218  fitPars[7]*(fitPars[6]*fitPars[1]+ugl[i])/(fitPars[1]+fitPars[4]),
2219  (fitPars[1]+ugl[i])/(fitPars[1]+fitPars[4])),
2220  cut,"goff");
2221  */
2222  nMeasuresI = SCi[i]->GetSelectedRows();
2223  memcpy(&(xvals[n][nMeasures]),SCi[i]->GetV1(),nMeasuresI*sizeof(Double_t));
2224  nMeasures += nMeasuresI;
2225  }
2226  }
2227  }
2228  int count,Count = nMeasures;
2229  for (count=0; count<Count;count++) {
2230  for (n=0; n<N1; n++) x[n] = xvals[n][count];
2231  ppl[0]->AddRow(x);
2232  }
2233  unsigned int combo = (unsigned int) (TMath::Power(2,N1-1));
2234  if (debug>0) printf("Npoints = %d, Nsca = %d, Possible combinations = %d\n",Count,N1-1,combo-1);
2235 
2236  // For bug: https://sft.its.cern.ch/jira/browse/ROOT-8238
2237  if (debug>1 && Count<100 && gROOT->GetVersionInt()<60900)
2238  printf("Please ignore any warnings about nbins from TH1::TH1 ...\n");
2239 
2240  // Try all combinations
2241  unsigned int l,m,one = 1;
2242  l = 0;
2243  rmsl[l] = 1e28;
2244  double rmsmin = rmsl[l];
2245  for (m=1;m<combo;m++) {
2246 
2247  int N = 0;
2248  for (n=0;n<N1;n++) if ((n==N1-1) || ((m>>n) & one)) { map[N]=n; N++; }
2249  if (N>Nmax) continue;
2250  l++;
2251  if (l == (unsigned int) maxp) {
2252  printf("Error! Maximum number of combinations %d exceeded!\n",maxp);
2253  break;
2254  }
2255 
2256  pp = new TPrincipal(N,"D");
2257  pp->SetName(Form("principal%d",l));
2258  ppl[l] = pp;
2259  for (count=0;count<Count;count++) {
2260  for (n=0; n<N; n++) x[n] = xvals[map[n]][count];
2261  pp->AddRow(x);
2262  }
2263  pp->MakePrincipals();
2264  pp->MakeHistograms(Form("PCA_SC%d",l),(debug>1 ? "xp" : "p"));
2265  TH1F* pca_sc = (TH1F*) (gROOT->FindObject(Form("PCA_SC%d_p%03d",l,N-1)));
2266 
2267  const Double_t* mean = pp->GetMeanValues()->GetMatrixArray();
2268  const Double_t* evals = pp->GetEigenValues()->GetMatrixArray();
2269  const Double_t* evecs = pp->GetEigenVectors()->GetMatrixArray();
2270  const Double_t* cova = pp->GetCovarianceMatrix()->GetMatrixArray();
2271  for (n=0; n<N-1; n++) coef[n] = -evecs[(n+1)*N - 1]/evecs[N*N-1];
2272 
2273  double trace = 0;
2274  for (n=0; n<N; n++) trace += cova[n*(N+1)];
2275 
2276  Double_t val_at_means = 0;
2277  for (n=0;n<N-1;n++) val_at_means += mean[n]*coef[n];
2278  Double_t zero_offset = val_at_means - mean[N-1];
2279  Double_t converted = zero_offset / coef[0];
2280 
2281  rmsl[l] = pca_sc->GetRMS();
2282 
2283  if (rmsl[l]<1e-8) {
2284  printf("TRYING TO FIND MY OWN RMS...%g == ",rmsl[l]);
2285  rmsp.Reset();
2286  for (int row=0;row<count;row++) {
2287  pp->X2P(pp->GetRow(row),p);
2288  rmsp.Fill(p[N-1]);
2289  }
2290  rmsl[l] = rmsp.GetRMS();
2291  printf("%g\n",rmsl[l]);
2292  }
2293 
2294  scl[l] = Form("RMS = %g :: ",rmsl[l]);
2295  scl[l] += Form("EV = %g :: ",evals[N-1]);
2296  scl[l] += Form("\n %s = ",scaA->At(N1-1)->GetName());
2297  scb[l] = Form("(%g*(%s",coef[0],scaA->At(map[0])->GetName());
2298  if (zero_offset!=0) {
2299  scb[l] += (converted>0 ? Form("-%g",converted) : Form("+%g",-converted));
2300  }
2301  scb[l] += "))";
2302  for (n=1;n<N-1;n++) scb[l] += Form("+(%g*(%s))",coef[n],scaA->At(map[n])->GetName());
2303  scl[l] += scb[l];
2304 
2305  if (rmsl[l] > 0 && rmsl[l] < rmsmin) {
2306  rmsmin = rmsl[l];
2307  pcaoffset = zero_offset;
2308  pcaN = N-1;
2309  for (n=0;n<N-1;n++) {
2310  pcadets[n] = scaA->At(map[n])->GetName();
2311  pcacoef[n] = coef[n];
2312  }
2313  }
2314 
2315 
2316  } // m loop
2317  if (debug>0) printf("Max limit on Nsca = %d, Explored combinations = %d\n",
2318  TMath::Min(Nmax,N1)-1,l);
2319  l++;
2320 
2321  TMath::Sort((long long) l,rmsl,rmsi);
2322 
2323  if (debug>1) {
2324  for (m=1;m<l;m++) {
2325  //printf("%d :: %d :: %s\n",l,rmsi[l],scl[rmsi[l]].Data());
2326  printf("%s\n",scl[rmsi[m]].Data());
2327  }
2328  }
2329 
2330  if (debug<1) {
2331  for (m=0;m<l;m++) {
2332  delete ppl[m]; // cleans up histograms too
2333  ppl[m]=0;
2334  }
2335  }
2336 
2337  for (n=0; n<N1; n++) delete xvals[n];
2338 
2339 
2340  return scb[rmsi[l-1]];
2341 }
2342 
2343 void PrintResult(double scp, double escp, double sop, double esop,
2344  double glp, double eglp, double ewp, double eewp, const char* det) {
2345  printf("sc = (");
2346  if (!DO_PCA && det) {
2347  printf("%6.4g +/- %6.4g) * ((%s) - (%6.4g +/- %6.4g)",scp,escp,det,sop,esop);
2348  } else {
2349  double lsop = (sop + pcaoffset)/pcacoef[0];
2350  double lesop = esop/pcacoef[0];
2351  printf("1.0 +/- %6.4g)*(%6.4g*(%s-(%6.4g +/- %6.4g))",escp/scp,pcacoef[0]*scp,
2352  pcadets[0].Data(),lsop,lesop);
2353  for (int n=1;n<pcaN;n++) printf("+(%6.4g*(%s))",pcacoef[n]*scp,pcadets[n].Data());
2354  }
2355  printf(")\n");
2356  if (EW_ASYMMETRY) printf(" with EWratio = %5.3f +/-%5.3f\n",ewp,eewp);
2357  printf(" with GL = %5.2f +/- %5.2f\n\n",glp,eglp);
2358 }
2359 
2360 void PrintResult(double scp, double escp, double sop, double esop,
2361  double* glp, double* eglp, const char* det) {
2362  printf("sc = (");
2363  if (!DO_PCA && det) {
2364  printf("%6.4g +/- %6.4g) * ((%s) - (%6.4g +/- %6.4g)",scp,escp,det,sop,esop);
2365  } else {
2366  double lsop = (sop + pcaoffset)/pcacoef[0];
2367  double lesop = esop/pcacoef[0];
2368  printf("1.0 +/- %6.4g)*(%6.4g*(%s-(%6.4g +/- %6.4g))",escp/scp,pcacoef[0]*scp,
2369  pcadets[0].Data(),lsop,lesop);
2370  for (int n=1;n<pcaN;n++) printf("+(%6.4g*(%s))",pcacoef[n]*scp,pcadets[n].Data());
2371  }
2372  printf(")\n with...\n");
2373  for (int n=0;n<12; n++) {
2374  if (!unfittableSec[n])
2375  printf("GL (%2d) = %5.2f +/- %5.2f\n",n+(EWmode<0 ? 13 : 1),glp[n],eglp[n]);
2376  }
2377  printf("\n");
2378 }
2379 
2381 // Revision 2.14 2021/09/23 genevb
2382 // A few improvements:
2383 // - Do not execute MINOS error calculation for parameters who are
2384 // at their limits from the standard fit as this was allowing the
2385 // limits to be exceeded
2386 // - PCA : check on sanity of g5 from iteration 0 before using in iteration 1
2387 // - FIX_GL : 3D fit for GridLeak now used only to provide g2 in this case
2388 //
2389 // Revision 2.13 2021/08/13 genevb
2390 // First git version, which doesn't autoincrement revision numbering
2391 // Several small improvements:
2392 // - sce_init calculation more robust by avoiding small denominators
2393 // - apply all initial guesses for gapf params if fit fails for FIX_GL
2394 // - VAR must be re-calculated at the end of FitWithOutlierRemoval()
2395 // because the individual parameter fits can change the results
2396 // (apparently for the better from what I can tell)
2397 // - function lines on plots vs. luminosity need to go to negative
2398 // luminosities in PCA fits for data with observed negative SpaceCharge
2399 //
2400 // $Id: Calib_SC_GL.C,v 2.12 2021/04/21 21:51:52 genevb Exp $
2401 // $Log: Calib_SC_GL.C,v $
2402 // Revision 2.12 2021/04/21 21:51:52 genevb
2403 // Better initial values for gapf fit, and don't do second PCA iteration if first fails (needed return values)
2404 //
2405 // Revision 2.11 2021/04/16 16:11:57 genevb
2406 // Introduce FIX_GL switch for fixed GL fits
2407 //
2408 // Revision 2.10 2019/10/21 15:20:57 genevb
2409 // Avoid Minuit complaint when asking about a fixed parameter
2410 //
2411 // Revision 2.9 2019/10/16 15:41:01 genevb
2412 // Fix issues with constant (fixed) parameter values
2413 //
2414 // Revision 2.8 2019/09/06 15:29:12 genevb
2415 // Include no-killer ZDC scalers in PCA
2416 //
2417 // Revision 2.7 2016/06/22 20:51:55 genevb
2418 // PCA: delete classes for debug=0, create x hists for debug>1 only
2419 //
2420 // Revision 2.6 2014/11/19 22:11:53 genevb
2421 // Print used ewratio from input files
2422 //
2423 // Revision 2.5 2014/11/19 19:19:56 genevb
2424 // Introduce EW asymmetry, and GL by sector
2425 //
2426 // Revision 2.4 2014/08/12 18:59:02 genevb
2427 // Introduce correlations between input datasets, MIGRAD=>MINIMIZE
2428 //
2429 // Revision 2.3 2014/06/10 19:16:16 genevb
2430 // Fix previous commit log: Introduce basic asymmetry, note of caution on 'guesses'
2431 //
2432 // Revision 2.2 2014/06/10 19:10:24 genevb
2433 // Introduce basic asymmetry, note of caution on 'guesses'
2434 //
2435 // Revision 2.1 2013/03/29 23:54:22 genevb
2436 // Modified fits with additional flexibility, use self-documented files for input, more scalers, bug fixes
2437 //
2438 // Revision 2.0 2012/06/13 04:45:02 genevb
2439 // Switch to universal fits, add PCA and flexible mixtures of scalers, chisq contour plots
2440 //
2441 // Revision 1.5 2010/01/06 17:56:43 genevb
2442 // Improvements to the fitter
2443 //
2444 // Revision 1.4 2008/04/30 14:54:38 genevb
2445 // Use gapf; get initial guess on linear fit
2446 //
2447 // Revision 1.3 2008/04/07 19:41:33 genevb
2448 // Minor updates to documentation, graphing
2449 //
2450 // Revision 1.2 2008/04/01 22:17:52 genevb
2451 // Improvements in documentation, fitting, graphing
2452 //
2453 // Revision 1.1 2006/05/23 16:14:42 genevb
2454 // Moved macro to calib directory
2455 //
2456 // Revision 3.2 2006/01/11 19:02:35 genevb
2457 // Better documentation
2458 //
2459 // Revision 3.1 2006/01/07 00:35:14 genevb
2460 // Introduce macro
2461 //
TDataSet::First
virtual TDataSet * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TDataSet.cxx:403
elem
Definition: dbStruct.hh:78
Generated by   doxygen 1.8.5