ROOTPWA
massDepFit.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2009 Sebastian Neubert
4 //
5 // This file is part of rootpwa
6 //
7 // rootpwa is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // rootpwa is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //-----------------------------------------------------------
22 // File and Version Information:
23 // $Id$
24 //
25 // Description:
26 // fitting program for massdependent fit rootpwa
27 // minimizes massDepFitLikeli function
28 //
29 // BE CAREFULL: Apart from the dynamic width the sqrts of the phase space factors are needed
30 // So at most places in this code ps acctually is the sqrt!!!
31 //
32 // Author List:
33 // Sebastian Neubert TUM (original author)
34 //
35 //
36 //-----------------------------------------------------------
37 
38 
39 #include <iostream>
40 #include <iomanip>
41 #include <vector>
42 #include <algorithm>
43 #include <string>
44 #include <complex>
45 #include <cassert>
46 #include <time.h>
47 
48 #include "TTree.h"
49 #include "TF1.h"
50 #include "TFile.h"
51 #include "TGraph.h"
52 #include "TGraph2D.h"
53 #include "TGraphErrors.h"
54 #include "TMultiGraph.h"
55 #include "TString.h"
56 #include "TComplex.h"
57 #include "TRandom3.h"
58 #include "Math/Minimizer.h"
59 #include "Math/Factory.h"
60 
61 #include "reportingUtils.hpp"
62 #include "fitResult.h"
63 #include "pwacomponent.h"
64 #include "massDepFitLikeli.h"
65 #include "TStopwatch.h"
66 #include "libconfig.h++"
67 
68 #define MASSSCALE 0.001
69 
70 
71 using namespace std;
72 using namespace libconfig;
73 using namespace ROOT::Math;
74 using namespace rpwa;
75 
76 
77 void
78 usage(const string& progName,
79  const int errCode = 0)
80 {
81  cerr << "usage:" << endl
82  << progName
83  << " -c configfile -i inputfile [-o outfile -l # -u #"
84  << " -M minimizer [-m algorithm] -t # -q -h -P -C] [-S fitResultFiles]" << endl
85  << " where:" << endl
86  << " -c file path to config File" << endl
87  << " -i file path to input file" << endl
88  << " -o file path to output file (default: 'mDep.result.root')" << endl
89  // << " -r # rank of spin density matrix (default: 1)" << endl
90  << " -l # -u # lower and upper mass range used for fit" << endl
91  << " -M name minimizer (default: Minuit2)" << endl
92  << " -m name minimization algorithm (optional, default: Migrad)" << endl
93  << " available minimizers: Minuit: Migrad, Simplex, Minimize, Migrad_imp" << endl
94  << " Minuit2: Migrad, Simplex, Combined, Scan, Fumili" << endl
95  << " GSLMultiMin: ConjugateFR, ConjugatePR, BFGS, BFGS2, SteepestDescent" << endl
96  << " GSLMultiFit: -" << endl
97  << " GSLSimAn: -" << endl
98  << " Linear: Robust" << endl
99  << " Fumili: -" << endl
100  << " -t # minimizer tolerance (default: 1e-10)" << endl
101  << " -q run quietly (default: false)" << endl
102  << " -h print help" << endl
103  << " -P plotting only - no fit" << endl
104  << " -C switch OFF covariances between real and imag part" << endl
105  << " -S files Systematic error plotting. give list of files" << endl
106 
107  << endl;
108  exit(errCode);
109 }
110 
111 // function loops through fitResults and puts phasespace values into a graph for interpolation
112 // THIS CONTAINS NOW THE RIGHT VALUE NOT!!! THE SQRT!!!
113 TGraph*
114 getPhaseSpace(TTree* tree, TF1* fsps,const std::string& wave){
115  unsigned int n=tree->GetEntries();
116  TGraph* graph=new TGraph(n);
117  fitResult* res=0;
118  tree->SetBranchAddress("fitResult_v2",&res);
119  for(unsigned int i=0; i<n; ++i){
120  tree->GetEntry(i);
121  double m=res->massBinCenter();
122  double ps=res->phaseSpaceIntegral(wave);
123  ps*=ps; // remember that phaseSpaceIntegral returns sqrt of integral!!!
124  //ps*=fsps->Eval(m);
125  graph->SetPoint(i,m,ps);
126  }
127  return graph;
128 }
129 
130 // changes status of variables (fixed/released)
131 // fixed values from config remain fixed
132 // parameters are taken from current status of fitter
133 // level
134 // 0 = release only couplings
135 // 1 = release couplings and masses
136 // 2 = release couplings, masses and widths
137 void releasePars(Minimizer* minimizer, const pwacompset& compset,
138  const vector<string>& anchorwave_reso,
139  const vector<string>& anchorwave_channel,
140  int level){
141  // copy state
142  unsigned int npar=minimizer->NDim();
143  double par[npar];
144  for(unsigned int i=0;i<npar;++i)par[i]=minimizer->X()[i];
145  minimizer->Clear();
146 
147  unsigned int parcount=0;
148  for(unsigned int ic=0;ic<compset.n();++ic){
149  const pwacomponent& comp=*compset[ic];
150  TString name(comp.name());
151  double mmin,mmax,gmin,gmax;
152  comp.getLimits(mmin,mmax,gmin,gmax);
153  if(comp.fixM() || level==0)minimizer->SetFixedVariable(parcount,
154  (name+"_M").Data() ,
155  par[parcount]);
156  else minimizer->SetLimitedVariable(parcount,
157  (name+"_M").Data(),
158  par[parcount],
159  5.0,
160  mmin,mmax);
161  if(level==0 && !comp.fixM()) printInfo << minimizer->VariableName(parcount)
162  << " fixed to " << par[parcount] << endl;
163  ++parcount;
164  if(comp.fixGamma() || level < 2)minimizer->SetFixedVariable(parcount,
165  (name+"_Gamma").Data() ,
166  par[parcount]);
167  else minimizer->SetLimitedVariable(parcount,
168  (name+"_Gamma").Data(),
169  par[parcount],
170  5.0,
171  gmin,gmax);
172  if(level<2 && !comp.fixGamma()) printInfo << minimizer->VariableName(parcount)
173  << " fixed to " << par[parcount] << endl;
174  ++parcount;
175 
176  std::map<std::string,pwachannel >::const_iterator it=comp.channels().begin();
177  while(it!=comp.channels().end()){
178  minimizer->SetVariable(parcount,(name + "_ReC" + it->first).Data() , par[parcount], 10.0);
179  ++parcount;
180  // fix one phase
181  if(find(anchorwave_reso.begin(),anchorwave_reso.end(),name)!=anchorwave_reso.end() && find(anchorwave_channel.begin(),anchorwave_channel.end(),it->first)!=anchorwave_channel.end()){
182  minimizer->SetFixedVariable(parcount,(name + "_ImC" + it->first).Data() , 0.0);
183  }
184  else {minimizer->SetVariable(parcount,(name + "_ImC" + it->first).Data() , par[parcount], 0.10);}
185  ++parcount;
186  ++it;
187  } // end loop over channels
188  }// end loop over components
189  // set phase space
190  unsigned int nfreePS=compset.nFreePSPar();
191  for(unsigned int ifreePS=0;ifreePS<nfreePS;++ifreePS){
192  double val,lower,upper;
193  val=par[parcount];
194  compset.getFreePSLimits(ifreePS,lower,upper);
195  TString name("PSP_"); name+=+ifreePS;
196  minimizer->SetLimitedVariable(parcount,
197  name.Data(),
198  val, 0.0001 ,lower,upper);
199  }
200 
201 
202 
203  const unsigned int nfree=minimizer->NFree();
204  printInfo << nfree << " Free Parameters in fit" << endl;
205 
206 
207 }
208 
209 
210 
211 int
212 main(int argc,
213  char** argv)
214 {
215  // --------------------------------------------------------------------------
216  // internal parameters
217  const string valTreeName = "pwa";
218  const string valBranchName = "fitResult_v2";
219  // double defaultStartValue = 0.01;
220 // bool useFixedStartValues = false;
221 // double startValStep = 0.0005;
222  const unsigned int maxNmbOfIterations = 20000;
223  const bool runHesse = true;
224  const bool runMinos = false;
225  bool onlyPlotting = false;
226  bool sysPlotting = false;
227  bool doCov = true;
228 
229  //unsigned int maxParNameLength = 20; // maximum length of parameter names
230 // int startValSeed = 1234567;
231  // parse command line options
232  const string progName = argv[0];
233  double massBinMin = 0; // [MeV/c^2]
234  double massBinMax = 5000; // [MeV/c^2]
235 
236  string inFileName = "fitresult.root"; // input filename
237  string outFileName = "mDep.result.root"; // output filename
238  //string normIntFileName = ""; // file with normalization integrals
239  string minimizerType[2] = {"Minuit2", "Migrad"}; // minimizer, minimization algorithm
240  double minimizerTolerance = 1e-10; // minimizer tolerance
241  bool quiet = false;
242 
243  string configFile; // configuration file
244 
245 extern char* optarg;
246  extern int optind;
247  // extern int optind;
248  int ca;
249  while ((ca = getopt(argc, argv, "c:i:o:u:l:M:m:t:qhPCS")) != -1)
250  switch (ca) {
251  case 'c':
252  configFile = optarg;
253  break;
254  case 'o':
255  outFileName = optarg;
256  break;
257  case 'i':
258  inFileName = optarg;
259  break;
260  case 'M':
261  minimizerType[0] = optarg;
262  break;
263  case 'm':
264  minimizerType[1] = optarg;
265  break;
266  case 't':
267  minimizerTolerance = atof(optarg);
268  break;
269  case 'l':
270  massBinMin = atof(optarg);
271  break;
272  case 'u':
273  massBinMax = atof(optarg);
274  break;
275  case 'q':
276  quiet = true;
277  break;
278  case 'h':
279  usage(progName);
280  break;
281  case 'P':
282  onlyPlotting=true;
283  break;
284  case 'C':
285  doCov=false;
286  break;
287  case 'S':
288  sysPlotting=true;
289  break;
290  }
291 
292 
293  // open input file and get results tree
294  TFile* infile=TFile::Open(inFileName.c_str());
295  if(infile==NULL){
296  cerr << "Input file " << inFileName <<" not found."<< endl;
297  return 1;
298  }
299  TTree* tree=(TTree*)infile->Get(valTreeName.c_str());
300  if(tree==NULL){
301  cerr << "Input tree " << valTreeName <<" not found."<< endl;
302  return 1;
303  }
304 
305 
306  vector<TTree*> sysTrees;
307  if(sysPlotting){
308  // add this fit
309  sysTrees.push_back(tree);
310  // open files with fits
311  for(int i=optind;i<argc;++i){
312  // open input file and get results tree
313  TFile* infile=TFile::Open(argv[i]);
314  if(infile==NULL){
315  cerr << "Systematics Input file " << inFileName <<" not found."<< endl;
316  return 1;
317  }
318  TTree* systree=(TTree*)infile->Get(valTreeName.c_str());
319  if(systree==NULL){
320  cerr << "Input tree " << valTreeName <<" not found."<< endl;
321  return 1;
322  }
323  sysTrees.push_back(systree);
324  }
325  printInfo << sysTrees.size() << " files for systematics found " << endl;
326  }// end if sysPlotting
327 
328  printInfo << "creating and setting up likelihood function" << endl;
329  printInfo << "doCovariances = " << doCov << endl;
330 
331 
332 
333  // Setup Component Set (Resonances + Background)
334  pwacompset compset;
335  Config Conf;
336  Conf.readFile(configFile.c_str());
337  const Setting& root = Conf.getRoot();
338  bool check=true;
339 
340  // overall Phasespace
341  double gps=0;
342  double pslower=0;
343  double psupper=0;
344  if(Conf.exists("ps")){
345  const Setting &pss = root["ps"];
346  check&=pss.lookupValue("formfactor",gps);
347  check&=pss.lookupValue("lower",pslower);
348  check&=pss.lookupValue("upper",psupper);
349 
350  }
351  // add overall phase space
352  TF1* fPS=new TF1("fps","[1] * ((x-[0])/1000.)^5*(1.+((x-[0])/1000.)*[2])*exp(-[3]* ((x-[0])/1000.)^2)",900,3000);
353  fPS->SetParameter(0,698); //5pi threshold
354  fPS->SetParLimits(0,698,698);
355  fPS->SetParameter(1,308.7E-6); // normalization
356  fPS->SetParLimits(1,308.7E-6,308.7E-6); //
357  fPS->SetParameter(2,4.859); // correction
358  fPS->SetParLimits(2,4.859,4.859);
359  fPS->SetParameter(3,gps); // Damping
360  fPS->SetParLimits(3,pslower,psupper);
361  //fPS->SetParameter(3,10); // Normalization
362  //TF1* fPS=new TF1("fps","1.",900,4000);
363 
364 
365  // Resonances
366  if(Conf.exists("components.resonances")){
367  const Setting &bws = root["components"]["resonances"];
368  // loop through breitwigners
369  int nbw=bws.getLength();
370  printInfo << "found " << nbw << " Resonances in config" << endl;
371  for(int ibw = 0; ibw < nbw; ++ibw) {
372  const Setting &bw = bws[ibw];
373  string jpc;
374  string name;
375  double mass=-1;double ml,mu;int mfix;
376  double width=-1;double wl,wu;int wfix;int wdyn;
377 
378  check&=bw.lookupValue("name", name);
379  check&=bw.lookupValue("jpc", jpc);
380  const Setting &massSet = bw["mass"];
381  check&=massSet.lookupValue("val", mass);
382  check&=massSet.lookupValue("lower", ml);
383  check&=massSet.lookupValue("upper", mu);
384  check&=massSet.lookupValue("fix", mfix);
385  const Setting &widthSet = bw["width"];
386  check&=widthSet.lookupValue("val", width);
387  check&=widthSet.lookupValue("lower", wl);
388  check&=widthSet.lookupValue("upper", wu);
389  check&=widthSet.lookupValue("fix", wfix);
390  bool checkdyn=widthSet.lookupValue("dyn", wdyn);
391  if(!checkdyn)wdyn=0;
392  cout << "---------------------------------------------------------------------" << endl;
393  cout << name << " JPC = " << jpc << endl;
394  cout << "mass(limits) = " << mass <<" ("<<ml<<","<<mu<<") MeV/c^2";
395  if(mfix==1)cout<<" -- FIXED";
396  cout<< endl;
397  cout << "width(limits) = " << width <<" ("<<wl<<","<<wu<<") MeV/c^2";
398  if(wfix==1)cout<<" -- FIXED";
399  if(wdyn!=0)cout<<" -- DYNAMIC WIDTH";
400  else cout<<" -- CONST WIDTH";
401  cout<< endl;
402  const Setting &channelSet = bw["decaychannels"];
403  unsigned int nCh=channelSet.getLength();
404  cout << "Decaychannels (coupling):" << endl;
405  std::map<std::string,pwachannel > channels;
406  for(unsigned int iCh=0;iCh<nCh;++iCh){
407  const Setting &ch = channelSet[iCh];
408  string amp;
409  double cRe=0;
410  double cIm=0;
411  check&=ch.lookupValue("amp",amp);
412  check&=ch.lookupValue("coupling_Re",cRe);
413  check&=ch.lookupValue("coupling_Im",cIm);
414  complex<double> C(cRe,cIm);
415  cout << " " << amp << " " << C << endl;
416  channels[amp]=pwachannel(C,getPhaseSpace(tree,fPS,amp));
417  }// end loop over channels
418  if(!check){
419  printErr << "Bad config value lookup! Check your config file!" << endl;
420  return 1;
421  }
422  pwacomponent* comp1=new pwacomponent(name,mass,width,channels);
423  cerr << "created component" << endl;
424  comp1->setLimits(ml,mu,wl,wu);
425  comp1->setFixed(mfix,wfix);
426  if(wdyn==0)comp1->setConstWidth();
427  compset.add(comp1);
428  cout << "CHECK val(m0)="<< comp1->val(mass) << endl;
429  }// end loop over resonances
430  }
431  cout << endl;
432  // Background components
433  if(Conf.exists("components.background")){
434  const Setting &bws = root["components"]["background"];
435  // loop through breitwigners
436  int nbw=bws.getLength();
437  printInfo << "found " << nbw << " Background components in config" << endl;
438  for(int ibw = 0; ibw < nbw; ++ibw) {
439  const Setting &bw = bws[ibw];
440  string name;
441  double mass=-1;double ml,mu;int mfix;
442  double width=-1;double wl,wu;int wfix;
443 
444  check&=bw.lookupValue("name", name);
445  const Setting &massSet = bw["m0"];
446  check&=massSet.lookupValue("val", mass);
447  check&=massSet.lookupValue("lower", ml);
448  check&=massSet.lookupValue("upper", mu);
449  check&=massSet.lookupValue("fix", mfix);
450  const Setting &widthSet = bw["g"];
451  check&=widthSet.lookupValue("val", width);
452  check&=widthSet.lookupValue("lower", wl);
453  check&=widthSet.lookupValue("upper", wu);
454  check&=widthSet.lookupValue("fix", wfix);
455  cout << "---------------------------------------------------------------------" << endl;
456  cout << name << endl;
457  cout << "mass-offset(limits) = " << mass <<" ("<<ml<<","<<mu<<") MeV/c^2";
458  if(mfix==1)cout<<" -- FIXED";
459  cout<< endl;
460  cout << "g(limits) = " << width <<" ("<<wl<<","<<wu<<") MeV/c^2";
461  if(wfix==1)cout<<" -- FIXED";
462  cout<< endl;
463  std::map<std::string,pwachannel > channels;
464  string amp;
465  double cRe=0;
466  double cIm=0;
467  double mIso1=0;
468  double mIso2=0;
469  check&=bw.lookupValue("amp",amp);
470  check&=bw.lookupValue("coupling_Re",cRe);
471  check&=bw.lookupValue("coupling_Im",cIm);
472  check&=bw.lookupValue("mIsobar1",mIso1);
473  check&=bw.lookupValue("mIsobar2",mIso2);
474  complex<double> C(cRe,cIm);
475  cout << "Decaychannel (coupling):" << endl;
476  cout << " " << amp << " " << C << endl;
477  cout << " Isobar masses: " << mIso1<<" "<< mIso2<< endl;
478  channels[amp]=pwachannel(C,getPhaseSpace(tree,fPS,amp));
479 
480  if(!check){
481  printErr << "Bad config value lookup! Check your config file!" << endl;
482  return 1;
483  }
484  pwabkg* bkg=new pwabkg(name,mass,width,channels);
485  bkg->setIsobars(mIso1,mIso2);
486  bkg->setLimits(ml,mu,wl,wu);
487  bkg->setFixed(mfix,wfix);
488  compset.add(bkg);
489  }// end loop over background
490  }// endif
491 
492 
493  cout << "---------------------------------------------------------------------" << endl << endl;
494 
495  compset.setPS(fPS);
496  compset.doMapping();
497 
498  cout << "---------------------------------------------------------------------" << endl << endl;
499 
500  // set anchorwave
501  vector<string> anchorwave_channel;
502  vector<string> anchorwave_reso;
503  if(Conf.exists("components.anchorwave")){
504  const Setting &anc = root["components"]["anchorwave"];
505  // loop through breitwigners
506  unsigned int nanc=anc.getLength();
507  for(unsigned int ianc=0;ianc<nanc;++ianc){
508  string ch,re;
509  const Setting &anco = anc[ianc];
510  anco.lookupValue("channel",ch);
511  anco.lookupValue("resonance",re);
512  cout << "Ancorwave: "<< endl;
513  cout << " " << re << endl;
514  cout << " " << ch << endl;
515  anchorwave_channel.push_back(ch);
516  anchorwave_reso.push_back(re);
517  }
518  }
519 
520 
521  cout << "---------------------------------------------------------------------" << endl << endl;
522 
523 
524 
525  massDepFitLikeli L;
526  L.init(tree,fPS,&compset,massBinMin,massBinMax,doCov);
527 
528  const unsigned int nmbPar = L.NDim();
529  // double par[nmbPar];
530  // for(unsigned int ip=0;ip<nmbPar;++ip)par[ip]=1.4;
531 
532 
533  // TStopwatch watch;
534  // L.DoEval(par);
535  // watch.Stop();
536 
537 
538  //printInfo << "TESTCALL TO LIKELIHOOD takes " << maxPrecisionAlign(watch.CpuTime()) << " s" << endl;
539 
540  printInfo << nmbPar << " Parameters in fit" << endl;
541 
542  // ---------------------------------------------------------------------------
543  // setup minimizer
544  printInfo << "creating and setting up minimizer " << minimizerType[0] << " using algorithm " << minimizerType[1] << endl;
545  Minimizer* minimizer = Factory::CreateMinimizer(minimizerType[0], minimizerType[1]);
546  if (!minimizer) {
547  printErr << "could not create minimizer! exiting!" << endl;
548  throw;
549  }
550  minimizer->SetFunction(L);
551  minimizer->SetPrintLevel((quiet) ? 0 : 3);
552 
553  // ---------------------------------------------------------------------------
554 
555  // Set startvalues
556  unsigned int parcount=0;
557  for(unsigned int ic=0;ic<compset.n();++ic){
558  const pwacomponent& comp=*compset[ic];
559  TString name(comp.name());
560  double mmin,mmax,gmin,gmax;
561  comp.getLimits(mmin,mmax,gmin,gmax);
562  if(comp.fixM())minimizer->SetFixedVariable(parcount++,
563  (name+"_M").Data() ,
564  comp.m0());
565  else minimizer->SetLimitedVariable(parcount++,
566  (name+"_M").Data(),
567  comp.m0(),
568  0.10,
569  mmin,mmax);
570  if(comp.fixGamma())minimizer->SetFixedVariable(parcount++,
571  (name+"_Gamma").Data() ,
572  comp.gamma());
573  else minimizer->SetLimitedVariable(parcount++,
574  (name+"_Gamma").Data(),
575  comp.gamma(),
576  0.01,
577  gmin,gmax);
578  std::map<std::string,pwachannel >::const_iterator it=comp.channels().begin();
579  while(it!=comp.channels().end()){
580  minimizer->SetVariable(parcount++,(name + "_ReC" + it->first).Data() , it->second.C().real(), 0.10);
581 
582  // fix one phase
583  if(find(anchorwave_reso.begin(),anchorwave_reso.end(),name)!=anchorwave_reso.end() && find(anchorwave_channel.begin(),anchorwave_channel.end(),it->first)!=anchorwave_channel.end()){
584  minimizer->SetFixedVariable(parcount++,(name + "_ImC" + it->first).Data() , 0.0);
585  }
586  else {minimizer->SetVariable(parcount++,(name + "_ImC" + it->first).Data() , it->second.C().imag(), 0.10);}
587 
588  ++it;
589  } // end loop over channels
590  }// end loop over components
591  // set phase space
592  unsigned int nfreePS=compset.nFreePSPar();
593  for(unsigned int ifreePS=0;ifreePS<nfreePS;++ifreePS){
594  double val,lower,upper;
595  val=compset.getFreePSPar(ifreePS);
596  compset.getFreePSLimits(ifreePS,lower,upper);
597  TString name("PSP_"); name+=+ifreePS;
598  minimizer->SetLimitedVariable(parcount++,
599  name.Data(),
600  val, 0.0001 ,lower,upper);
601  }
602 
603 
604 
605  const unsigned int nfree=minimizer->NFree();
606  printInfo << nfree << " Free Parameters in fit" << endl;
607 
608 
609  // find minimum of likelihood function
610  double chi2=0;
611  if(onlyPlotting) printInfo << "Plotting mode, skipping minimzation!" << endl;
612  else {
613  printInfo << "performing minimization. MASSES AND WIDTHS FIXED" << endl;
614 
615  minimizer->SetMaxIterations(maxNmbOfIterations);
616  minimizer->SetMaxFunctionCalls(maxNmbOfIterations*5);
617  minimizer->SetTolerance (minimizerTolerance);
618  // only do couplings
619  TStopwatch fitW;
620  // releasePars(minimizer,compset,anchorwave_reso,anchorwave_channel,0);
621  bool success = minimizer->Minimize();
622  if(!success)printWarn << "minimization failed." << endl;
623  else printInfo << "minimization successful." << endl;
624  printInfo << "Minimization took " << maxPrecisionAlign(fitW.CpuTime()) << " s" << endl;
625  //release masses
626  releasePars(minimizer,compset,anchorwave_reso,anchorwave_channel,1);
627  printInfo << "performing minimization. MASSES RELEASED" << endl;
628  fitW.Start();
629  success &= minimizer->Minimize();
630  if(!success)printWarn << "minimization failed." << endl;
631  else printInfo << "minimization successful." << endl;
632  printInfo << "Minimization took " << maxPrecisionAlign(fitW.CpuTime()) << " s" << endl;
633  //release widths
634  releasePars(minimizer,compset,anchorwave_reso,anchorwave_channel,2);
635  printInfo << "performing minimization. ALL RELEASED" << endl;
636  fitW.Start();
637  success &= minimizer->Minimize();
638  printInfo << "Minimization took " << maxPrecisionAlign(fitW.CpuTime()) << " s" << endl;
639 
640  const double* par=minimizer->X();
641  compset.setPar(par);
642  cerr << compset << endl;
643  if (success){
644  printInfo << "minimization finished successfully." << endl;
645  chi2=minimizer->MinValue();
646  }
647  else
648  printWarn << "minimization failed." << endl;
649  if (runHesse) {
650  printInfo << "calculating Hessian matrix." << endl;
651  success = minimizer->Hesse(); // comes only with ROOT 5.24+
652  if (!success)
653  printWarn << "calculation of Hessian matrix failed." << endl;
654  }
655  }
656  printInfo << "minimization stopped after " << minimizer->NCalls() << " function calls. minimizer status summary:" << endl
657  << " total number of parameters .......................... " << minimizer->NDim() << endl
658  << " number of free parameters ........................... " << minimizer->NFree() << endl
659  << " maximum allowed number of iterations ................ " << minimizer->MaxIterations() << endl
660  << " maximum allowed number of function calls ............ " << minimizer->MaxFunctionCalls() << endl
661  << " minimizer status .................................... " << minimizer->Status() << endl
662  << " minimizer provides error and error matrix ........... " << minimizer->ProvidesError() << endl
663  << " minimizer has performed detailed error validation ... " << minimizer->IsValidError() << endl
664  << " estimated distance to minimum ....................... " << minimizer->Edm() << endl
665  << " statistical scale used for error calculation ........ " << minimizer->ErrorDef() << endl
666  << " minimizer strategy .................................. " << minimizer->Strategy() << endl
667  << " absolute tolerance .................................. " << minimizer->Tolerance() << endl;
668 
669 
670  // ---------------------------------------------------------------------------
671  // print results
672  //map<TString, double> errormap;
673  printInfo << "minimization result:" << endl;
674  for (unsigned int i = 0; i< nmbPar; ++i) {
675  cout << " parameter [" << setw(3) << i << "] ";
676  cout << minimizer->VariableName(i) << " " ;
677  // << setw(maxParNameLength); //<< L.parName(i) << " = ";
678  //if (parIsFixed[i])
679  // cout << minimizer->X()[i] << " (fixed)" << endl;
680  //else {
681  cout << setw(12) << maxPrecisionAlign(minimizer->X()[i]) << " +- "
682  << setw(12) << maxPrecisionAlign(minimizer->Errors()[i]);
683  //errormap[minimizer]=minimizer->Errors()[i];
684 
685 
686  if (runMinos && (i == 156)) { // does not work for all parameters
687  double minosErrLow = 0;
688  double minosErrUp = 0;
689  const bool success = minimizer->GetMinosError(i, minosErrLow, minosErrUp);
690  if (success)
691  cout << " Minos: " << "[" << minosErrLow << ", +" << minosErrUp << "]" << endl;
692  } else
693  cout << endl;
694  }
695 
696  cout << "---------------------------------------------------------------------" << endl;
697  // Reduced chi2
698 
699  printInfo << chi2 << " chi2" << endl;
700  unsigned int numdata=L.NDataPoints();
701  // numDOF
702  unsigned int numDOF=numdata-nfree;
703  printInfo << numDOF << " degrees of freedom" << endl;
704  double redChi2 = chi2/(double)numDOF;
705  printInfo << redChi2 << " chi2/nDF" << endl;
706  cout << "---------------------------------------------------------------------" << endl;
707 
708 
709  // write out results
710  // Likelihood and such
711  const Setting& fitqualS= root["fitquality"];
712  Setting& chi2S=fitqualS["chi2"];
713  chi2S=chi2;
714  Setting& ndfS=fitqualS["ndf"];
715  ndfS=(int)numDOF;
716  Setting& redchi2S=fitqualS["redchi2"];
717  redchi2S=redChi2;
718  const Setting& psS= root["ps"];
719  Setting& ffS=psS["formfactor"];
720  ffS=(double)fPS->GetParameter(3);
721 
722  if(nfreePS>0){
723  Setting& ffeS=psS["error"];
724  ffeS=minimizer->Errors()[minimizer->VariableIndex("PSP_0")];
725  }
726 
727  // Setup Component Set (Resonances + Background)
728  const Setting& bws= root["components"]["resonances"];
729  const Setting& bkgs= root["components"]["background"];
730  unsigned int nbws=bws.getLength();
731  unsigned int nbkgs=bkgs.getLength();
732  // loop over components
733  unsigned int nc=compset.n();
734  for(unsigned int ic=0;ic<nc;++ic){
735  const pwacomponent* comp=compset[ic];
736  string name=comp->name();
737  // search corresponding setting
738 
739  string sname;
740  bool found=false;
741  for(unsigned int is=0;is<nbws;++is){
742  const Setting& bw = bws[is];
743  bw.lookupValue("name", sname);
744  if(sname==name){
745  found=true;
746  // set values to this setting
747  Setting& sm = bw["mass"];
748  Setting& smval = sm["val"];
749  smval = comp->m0();
750  Setting& smerr = sm["error"];
751  TString merrname=name+"_M";
752  smerr=minimizer->Errors()[minimizer->VariableIndex(merrname.Data())];
753 
754  Setting& sw = bw["width"];
755  Setting& swval = sw["val"];
756  swval = comp->gamma();
757 
758  Setting& swerr = sw["error"];
759  TString werrname=name+"_Gamma";
760  swerr=minimizer->Errors()[minimizer->VariableIndex(werrname.Data())];
761 
762  cout << name
763  << " mass="<<double(smval)<<" +- "<<double(smerr)
764  << " width="<<double(swval)<<" +- "<<double(swerr)<< endl;
765 
766  // loop through channel and fix couplings
767  const Setting& sChn=bw["decaychannels"];
768  unsigned int nCh=sChn.getLength();
769  const std::map<std::string,pwachannel >& ch=comp->channels();
770  std::map<std::string,pwachannel>::const_iterator it=ch.begin();
771  for(;it!=ch.end();++it){
772  std::complex<double> c= it->second.C();
773  string ampname=it->first;
774  // loop through channels in setting
775  for(unsigned int isc=0;isc<nCh;++isc){
776  Setting& sCh=sChn[isc];
777  string amp; sCh.lookupValue("amp",amp);
778  if(amp==ampname){
779  Setting& sRe=sCh["coupling_Re"];
780  sRe=c.real();
781  Setting& sIm=sCh["coupling_Im"];
782  sIm=c.imag();
783  break;
784  } // endif
785  } // end loop through cannels in setting
786 
787  } // end loop through channels of component
788 
789  break;
790  }
791  }
792 
793  // loop over background settings
794  if(!found){
795  for(unsigned int is=0;is<nbkgs;++is){
796  const Setting& bw = bkgs[is];
797  bw.lookupValue("name", sname);
798  if(sname==name){
799  Setting& sm = bw["m0"];
800  Setting& smval = sm["val"];
801  smval = comp->m0();
802  Setting& sw = bw["g"];
803  Setting& swval = sw["val"];
804  swval = comp->gamma();
805 
806  const pwachannel& ch=comp->channels().begin()->second;
807  std::complex<double> c=ch.C();
808  Setting& sRe=bw["coupling_Re"];
809  sRe=c.real();
810  Setting& sIm=bw["coupling_Im"];
811  sIm=c.imag();
812  break;
813  }
814  }
815  }
816 
817 
818 
819 
820  }
821 
822 
823 
824 
825  // bool check=true;
826 // // Resonances
827 // if(Conf.exists("components.resonances")){
828 // const Setting &bws = root["components"]["resonances"];
829 // // loop through breitwigners
830 // int nbw=bws.getLength();
831 // printInfo << "found " << nbw << " Resonances in config" << endl;
832 // for(int ibw = 0; ibw < nbw; ++ibw) {
833 // const Setting &bw = bws[ibw];
834 // string jpc;
835 // string name;
836 // double mass=-1;double ml,mu;int mfix;
837 // double width=-1;double wl,wu;int wfix;
838 
839 // check&=bw.lookupValue("name", name);
840 // check&=bw.lookupValue("jpc", jpc);
841 // const Setting &massSet = bw["mass"];
842 // check&=massSet.lookupValue("val", mass);
843 // check&=massSet.lookupValue("lower", ml);
844 // check&=massSet.lookupValue("upper", mu);
845 // check&=massSet.lookupValue("fix", mfix);
846 // const Setting &widthSet = bw["width"];
847 // check&=widthSet.lookupValue("val", width);
848 // check&=widthSet.lookupValue("lower", wl);
849 // check&=widthSet.lookupValue("upper", wu);
850 // check&=widthSet.lookupValue("fix", wfix);
851 // cout << "---------------------------------------------------------------------" << endl;
852 // cout << name << " JPC = " << jpc << endl;
853 // cout << "mass(limits) = " << mass <<" ("<<ml<<","<<mu<<") MeV/c^2";
854 // if(mfix==1)cout<<" -- FIXED";
855 // cout<< endl;
856 // cout << "width(limits) = " << width <<" ("<<wl<<","<<wu<<") MeV/c^2";
857 // if(wfix==1)cout<<" -- FIXED";
858 // cout<< endl;
859 // const Setting &channelSet = bw["decaychannels"];
860 // unsigned int nCh=channelSet.getLength();
861 // cout << "Decaychannels (coupling):" << endl;
862 // std::map<std::string,pwachannel > channels;
863 // for(unsigned int iCh=0;iCh<nCh;++iCh){
864 // const Setting &ch = channelSet[iCh];
865 // string amp;
866 // double cRe=0;
867 // double cIm=0;
868 // check&=ch.lookupValue("amp",amp);
869 // check&=ch.lookupValue("coupling_Re",cRe);
870 // check&=ch.lookupValue("coupling_Im",cIm);
871 // complex<double> C(cRe,cIm);
872 // cout << " " << amp << " " << C << endl;
873 // channels[amp]=pwachannel(C,getPhaseSpace(tree,amp));
874 // }// end loop over channels
875 // if(!check){
876 // printErr << "Bad config value lookup! Check your config file!" << endl;
877 // return 1;
878 // }
879 // pwacomponent* comp1=new pwacomponent(name,mass,width,channels);
880 // comp1->setLimits(ml,mu,wl,wu);
881 // comp1->setFixed(mfix,wfix);
882 // compset.add(comp1);
883 // }// end loop over resonances
884 // }
885 // cout << endl;
886 // // Background components
887 // if(Conf.exists("components.background")){
888 // const Setting &bws = root["components"]["background"];
889 // // loop through breitwigners
890 // int nbw=bws.getLength();
891 // printInfo << "found " << nbw << " Background components in config" << endl;
892 // for(int ibw = 0; ibw < nbw; ++ibw) {
893 // const Setting &bw = bws[ibw];
894 // string name;
895 // double mass=-1;double ml,mu;int mfix;
896 // double width=-1;double wl,wu;int wfix;
897 
898 // check&=bw.lookupValue("name", name);
899 // const Setting &massSet = bw["m0"];
900 // check&=massSet.lookupValue("val", mass);
901 // check&=massSet.lookupValue("lower", ml);
902 // check&=massSet.lookupValue("upper", mu);
903 // check&=massSet.lookupValue("fix", mfix);
904 // const Setting &widthSet = bw["g"];
905 // check&=widthSet.lookupValue("val", width);
906 // check&=widthSet.lookupValue("lower", wl);
907 // check&=widthSet.lookupValue("upper", wu);
908 // check&=widthSet.lookupValue("fix", wfix);
909 // cout << "---------------------------------------------------------------------" << endl;
910 // cout << name << endl;
911 // cout << "mass-offset(limits) = " << mass <<" ("<<ml<<","<<mu<<") MeV/c^2";
912 // if(mfix==1)cout<<" -- FIXED";
913 // cout<< endl;
914 // cout << "g(limits) = " << width <<" ("<<wl<<","<<wu<<") MeV/c^2";
915 // if(wfix==1)cout<<" -- FIXED";
916 // cout<< endl;
917 // std::map<std::string,pwachannel > channels;
918 // string amp;
919 // double cRe=0;
920 // double cIm=0;
921 // double mIso1=0;
922 // double mIso2=0;
923 // check&=bw.lookupValue("amp",amp);
924 // check&=bw.lookupValue("coupling_Re",cRe);
925 // check&=bw.lookupValue("coupling_Im",cIm);
926 // check&=bw.lookupValue("mIsobar1",mIso1);
927 // check&=bw.lookupValue("mIsobar2",mIso2);
928 // complex<double> C(cRe,cIm);
929 // cout << "Decaychannel (coupling):" << endl;
930 // cout << " " << amp << " " << C << endl;
931 // cout << " Isobar masses: " << mIso1<<" "<< mIso2<< endl;
932 // channels[amp]=pwachannel(C,getPhaseSpace(tree,amp));
933 
934 // if(!check){
935 // printErr << "Bad config value lookup! Check your config file!" << endl;
936 // return 1;
937 // }
938 // pwabkg* bkg=new pwabkg(name,mass,width,channels);
939 // bkg->setIsobars(mIso1,mIso2);
940 // bkg->setLimits(ml,mu,wl,wu);
941 // bkg->setFixed(mfix,wfix);
942 // compset.add(bkg);
943 // }// end loop over background
944 // }// endif
945 
946  string outconfig(outFileName);
947  outconfig.append(".conf");
948  Conf.writeFile(outconfig.c_str());
949 
950  cerr << "Fitting finished... Start building graphs ... " << endl;
951 
952  int syscolor=kAzure-9;
953 
954  std::vector<std::string> wl=compset.wavelist();
955  std::map<std::string, unsigned int> wmap;
956  unsigned int ndatabins=tree->GetEntries();
957 
958  std::vector<TGraphErrors*> datagraphs;
959  std::vector<TGraphErrors*> intenssysgraphs;
960 
961  std::vector<TMultiGraph*> graphs;
962 
963  for(unsigned int iw=0; iw<wl.size();++iw){
964  wmap[wl[iw]]=iw;
965  graphs.push_back(new TMultiGraph);
966 
967  intenssysgraphs.push_back(new TGraphErrors(ndatabins));
968  string name=("sys_");name.append(wl[iw]);
969  intenssysgraphs[iw]->SetName(name.c_str());
970  intenssysgraphs[iw]->SetTitle(name.c_str());
971  intenssysgraphs[iw]->SetLineColor(syscolor);
972  intenssysgraphs[iw]->SetFillColor(syscolor);
973  intenssysgraphs[iw]->SetDrawOption("2");
974  graphs[iw]->Add(intenssysgraphs[iw],"2");
975 
976 
977 
978  graphs[iw]->SetName(wl[iw].c_str());
979  graphs[iw]->SetTitle(wl[iw].c_str());
980  graphs[iw]->SetDrawOption("AP");
981  datagraphs.push_back(new TGraphErrors(ndatabins));
982  name="data_";name.append(wl[iw]);
983  datagraphs[iw]->SetName(name.c_str());
984  datagraphs[iw]->SetTitle(name.c_str());
985  datagraphs[iw]->SetDrawOption("AP");
986  //datagraphs[iw]->SetLineColor(kRed);
987  //datagraphs[iw]->SetMarkerColor(kRed);
988 
989  graphs[iw]->Add(datagraphs[iw],"P");
990  }
991 
992 
993 
994 
995  // build fitgraphs
996  unsigned int nbins=ndatabins;//200;
997  //double mmin=1200.;
998  //double md=10.;
999  std::vector<TGraph*> fitgraphs;
1000  std::vector<TGraph*> absphasegraphs;
1001 
1002 
1003  for(unsigned int iw=0; iw<wl.size();++iw){
1004  fitgraphs.push_back(new TGraph(nbins));
1005  string name("fit_");name.append(wl[iw]);
1006  fitgraphs[iw]->SetName(name.c_str());
1007  fitgraphs[iw]->SetTitle(name.c_str());
1008  fitgraphs[iw]->SetLineColor(kRed);
1009  fitgraphs[iw]->SetLineWidth(2);
1010  fitgraphs[iw]->SetMarkerColor(kRed);
1011  fitgraphs[iw]->SetDrawOption("AP");
1012  //fitgraphs[iw]->SetMarkerStyle(22);
1013  graphs[iw]->Add(fitgraphs[iw],"cp");
1014  graphs[iw]->Add(getPhaseSpace(tree,fPS,wl[iw]));
1015 
1016  absphasegraphs.push_back(new TGraph(nbins));
1017  name="absphase_";name.append(wl[iw]);
1018  absphasegraphs[iw]->SetName(name.c_str());
1019  absphasegraphs[iw]->SetTitle(name.c_str());
1020  absphasegraphs[iw]->SetLineColor(kRed);
1021  absphasegraphs[iw]->SetLineWidth(2);
1022  absphasegraphs[iw]->SetMarkerColor(kRed);
1023  absphasegraphs[iw]->SetDrawOption("AP");
1024 
1025 
1026 
1027  }
1028 
1029  std::vector<TGraph*> compgraphs; // individual components
1030  // loop over components and build graphs
1031  for(unsigned int ic=0;ic<compset.n();++ic){
1032  const pwacomponent* c=compset[ic];
1033  std::map<std::string,pwachannel >::const_iterator it=c->channels().begin();
1034  while(it!=c->channels().end()){
1035  string name=c->name();name.append("__");
1036  name.append(it->first);
1037  TGraph* gcomp=new TGraph(nbins);
1038  gcomp->SetName(name.c_str());
1039  gcomp->SetTitle(name.c_str());
1040  unsigned int color=kBlue;
1041  if(dynamic_cast<const pwabkg*>(c)!=NULL)color=kMagenta;
1042  gcomp->SetLineColor(color);
1043  gcomp->SetMarkerColor(color);
1044 
1045  compgraphs.push_back(gcomp);
1046  graphs[wmap[it->first]]->Add(gcomp,"cp");
1047  ++it;
1048  }// end loop over channels
1049 
1050  }// end loop over components
1051 
1052 
1053  std::vector<TGraphErrors*> phasedatagraphs;
1054  std::vector<TGraphErrors*> phasesysgraphs;
1055 
1056 
1057  std::vector<TGraphErrors*> realdatagraphs;
1058  std::vector<TGraphErrors*> realsysgraphs;
1059  std::vector<TGraphErrors*> imagdatagraphs;
1060  std::vector<TGraphErrors*> imagsysgraphs;
1061 
1062  std::vector<TGraph*> realfitgraphs;
1063  std::vector<TGraph*> imagfitgraphs;
1064 
1065  std::vector<TMultiGraph*> phasegraphs;
1066  std::vector<TMultiGraph*> overlapRegraphs;
1067  std::vector<TMultiGraph*> overlapImgraphs;
1068 
1069  std::vector<TGraph*> phasefitgraphs;
1070  unsigned int c=0;
1071 
1072 
1073  for(unsigned int iw=0; iw<wl.size();++iw){
1074  for(unsigned int iw2=iw+1; iw2<wl.size();++iw2){
1075 
1076 
1077 
1078  phasegraphs.push_back(new TMultiGraph);
1079  overlapImgraphs.push_back(new TMultiGraph);
1080  overlapRegraphs.push_back(new TMultiGraph);
1081  string name("dPhi_");name.append(wl[iw]);
1082  name.append("---");name.append(wl[iw2]);
1083  phasegraphs[c]->SetName(name.c_str());
1084  phasegraphs[c]->SetTitle(name.c_str());
1085  name="Re_";name.append(wl[iw]);
1086  name.append("---");name.append(wl[iw2]);
1087  overlapRegraphs[c]->SetName(name.c_str());
1088  overlapRegraphs[c]->SetTitle(name.c_str());
1089  name="Im_";name.append(wl[iw]);
1090  name.append("---");name.append(wl[iw2]);
1091  overlapImgraphs[c]->SetName(name.c_str());
1092  overlapImgraphs[c]->SetTitle(name.c_str());
1093 
1094  phasesysgraphs.push_back(new TGraphErrors(nbins));
1095  name=("dPhi_sys_");name.append(wl[iw]);
1096  name.append("---");name.append(wl[iw2]);
1097  phasesysgraphs[c]->SetName(name.c_str());
1098  phasesysgraphs[c]->SetTitle(name.c_str());
1099  phasesysgraphs[c]->SetLineColor(syscolor);
1100  phasesysgraphs[c]->SetFillColor(syscolor);
1101  phasesysgraphs[c]->SetDrawOption("2");
1102  //phasesysgraphs[c]->SetLineWidth(1);
1103  //phasesysgraphs[c]->SetFillStyle(1001);
1104 
1105  phasegraphs[c]->Add(phasesysgraphs[c],"2");
1106 
1107 
1108  phasedatagraphs.push_back(new TGraphErrors(nbins));
1109  name=("dPhi_data_");name.append(wl[iw]);
1110  name.append("---");name.append(wl[iw2]);
1111  phasedatagraphs[c]->SetName(name.c_str());
1112  phasedatagraphs[c]->SetTitle(name.c_str());
1113  phasegraphs[c]->Add(phasedatagraphs[c],"p");
1114 
1115 
1116  realsysgraphs.push_back(new TGraphErrors(nbins));
1117  name=("RE_sys_");name.append(wl[iw]);
1118  name.append("---");name.append(wl[iw2]);
1119  realsysgraphs[c]->SetName(name.c_str());
1120  realsysgraphs[c]->SetTitle(name.c_str());
1121  realsysgraphs[c]->SetLineColor(syscolor);
1122  realsysgraphs[c]->SetFillColor(syscolor);
1123  realsysgraphs[c]->SetDrawOption("2");
1124  overlapRegraphs[c]->Add(realsysgraphs[c],"2");
1125 
1126  imagsysgraphs.push_back(new TGraphErrors(nbins));
1127  name=("IM_sys_");name.append(wl[iw]);
1128  name.append("---");name.append(wl[iw2]);
1129  imagsysgraphs[c]->SetName(name.c_str());
1130  imagsysgraphs[c]->SetTitle(name.c_str());
1131  imagsysgraphs[c]->SetLineColor(syscolor);
1132  imagsysgraphs[c]->SetFillColor(syscolor);
1133  imagsysgraphs[c]->SetDrawOption("2");
1134  overlapImgraphs[c]->Add(imagsysgraphs[c],"2");
1135 
1136 
1137  realdatagraphs.push_back(new TGraphErrors(nbins));
1138  name=("RE_data_");name.append(wl[iw]);
1139  name.append("---");name.append(wl[iw2]);
1140  realdatagraphs[c]->SetName(name.c_str());
1141  realdatagraphs[c]->SetTitle(name.c_str());
1142  overlapRegraphs[c]->Add(realdatagraphs[c],"p");
1143 
1144  imagdatagraphs.push_back(new TGraphErrors(nbins));
1145  name=("IM_data_");name.append(wl[iw]);
1146  name.append("---");name.append(wl[iw2]);
1147  imagdatagraphs[c]->SetName(name.c_str());
1148  imagdatagraphs[c]->SetTitle(name.c_str());
1149  //imagdatagraphs[c]->SetLineStyle(2);
1150  overlapImgraphs[c]->Add(imagdatagraphs[c],"p");
1151 
1152  ++c;
1153  }
1154  }
1155 
1156  c=0;
1157  for(unsigned int iw=0; iw<wl.size();++iw){
1158  for(unsigned int iw2=iw+1; iw2<wl.size();++iw2){
1159  phasefitgraphs.push_back(new TGraph(nbins));
1160  string name("dPhi_fit_");name.append(wl[iw]);
1161  name.append("---");name.append(wl[iw2]);
1162  phasefitgraphs[c]->SetName(name.c_str());
1163  phasefitgraphs[c]->SetTitle(name.c_str());
1164  phasefitgraphs[c]->SetLineColor(kRed);
1165  phasefitgraphs[c]->SetMarkerColor(kRed);
1166  phasefitgraphs[c]->SetDrawOption("P");
1167  phasefitgraphs[c]->SetMarkerStyle(24);
1168  phasefitgraphs[c]->SetMarkerSize(0.2);
1169  phasegraphs[c]->Add(phasefitgraphs[c],"cp");
1170 
1171  realfitgraphs.push_back(new TGraph(nbins));
1172  name=("Re_fit_");name.append(wl[iw]);
1173  name.append("---");name.append(wl[iw2]);
1174  realfitgraphs[c]->SetName(name.c_str());
1175  realfitgraphs[c]->SetTitle(name.c_str());
1176  realfitgraphs[c]->SetLineColor(kRed);
1177  realfitgraphs[c]->SetMarkerColor(kRed);
1178  realfitgraphs[c]->SetDrawOption("AP");
1179  realfitgraphs[c]->SetMarkerStyle(24);
1180  realfitgraphs[c]->SetMarkerSize(0.2);
1181  overlapRegraphs[c]->Add(realfitgraphs[c],"cp");
1182 
1183  imagfitgraphs.push_back(new TGraph(nbins));
1184  name=("Im_fit_");name.append(wl[iw]);
1185  name.append("---");name.append(wl[iw2]);
1186  imagfitgraphs[c]->SetName(name.c_str());
1187  imagfitgraphs[c]->SetTitle(name.c_str());
1188  imagfitgraphs[c]->SetLineColor(kRed);
1189  //imagfitgraphs[c]->SetLineStyle(2);
1190  imagfitgraphs[c]->SetMarkerColor(kRed);
1191  imagfitgraphs[c]->SetDrawOption("AP");
1192  imagfitgraphs[c]->SetMarkerStyle(24);
1193  imagfitgraphs[c]->SetMarkerSize(0.2);
1194  overlapImgraphs[c]->Add(imagfitgraphs[c],"cp");
1195 
1196  ++c;
1197  }
1198  }
1199 
1200  std::vector<TGraph2D*> phase2d;
1201  c=0;
1202  for(unsigned int iw=0; iw<wl.size();++iw){
1203  for(unsigned int iw2=iw+1; iw2<wl.size();++iw2){
1204  phase2d.push_back(new TGraph2D(nbins));
1205  string name("dPhi_2d_");name.append(wl[iw]);
1206  name.append("---");name.append(wl[iw2]);
1207  phase2d[c]->SetName(name.c_str());
1208  phase2d[c]->SetTitle(name.c_str());
1209  phase2d[c]->SetLineColor(kRed);
1210  phase2d[c]->SetMarkerColor(kRed);
1211  //phasegraphs[c]->Add(phasefitgraphs[c],"cp");
1212  ++c;
1213  }
1214  }
1215 
1216 
1217 
1218 
1219 
1220 
1221  // get data
1222  fitResult* rho=0;
1223  tree->SetBranchAddress(valBranchName.c_str(),&rho);
1224  vector<double> prevps(wl.size());
1225  //double mprev=0;
1226  vector<double> prevphase(wl.size());
1227  double binwidth=MASSSCALE*30; // half binwidth
1228  //double w=2*30/10;
1229  for(unsigned int i=0;i<ndatabins;++i){
1230  tree->GetEntry(i);
1231  double m=rho->massBinCenter();
1232  double mScaled=m*MASSSCALE;
1233  //cout << "MASS: "<<m << endl;
1234  double fsps=sqrt(fPS->Eval(rho->massBinCenter()));
1235  unsigned int c=0;
1236  for(unsigned int iw=0; iw<wl.size();++iw){
1237 
1238  //cout << wl[iw] << endl;
1239 
1240  double ps=rho->phaseSpaceIntegral(wl[iw].c_str())*fsps;
1241 
1242  datagraphs[iw]->SetPoint(i,mScaled,rho->intensity(wl[iw].c_str()));
1243  datagraphs[iw]->SetPointError(i,binwidth,rho->intensityErr(wl[iw].c_str()));
1244  fitgraphs[iw]->SetPoint(i,mScaled,compset.intensity(wl[iw],m));
1245  double absphase=compset.phase(wl[iw],m)*TMath::RadToDeg();
1246  if(i>0){
1247  double absp=absphase+360;
1248  double absm=absphase-360;
1249  if(fabs(absphase-prevphase[iw])>fabs(absp-prevphase[iw])){
1250  absphase=absp;
1251  }
1252  else if(fabs(absphase-prevphase[iw])>fabs(absm-prevphase[iw])){
1253  absphase=absm;
1254  }
1255  }
1256  prevphase[iw]=absphase;
1257  absphasegraphs[iw]->SetPoint(i,mScaled,absphase);
1258  if(sysPlotting){
1259  double maxIntens=-10000000;
1260  double minIntens=10000000;
1261  for(unsigned int iSys=0;iSys<sysTrees.size();++iSys){
1262  // get data
1263  fitResult* rhoSys=0;
1264  if(iSys==0)rhoSys=rho;
1265  else {
1266  sysTrees[iSys]->SetBranchAddress(valBranchName.c_str(),&rhoSys);
1267  sysTrees[iSys]->GetEntry(i);
1268  }
1269 // rename one wave
1270  string mywave1=wl[iw];
1271 
1272 
1273  if(mywave1=="1-1++0+pi-_01_rho1700=pi-+_10_pi1300=pi+-_00_sigma.amp" && iSys>0)mywave1="1-1++0+pi-_01_eta11600=pi-+_10_pi1300=pi+-_00_sigma.amp";
1274 
1275  // check if waves are in fit
1276  if(rhoSys->waveIndex(mywave1)==-1){
1277  delete rhoSys;
1278  continue;
1279  }
1280  double myI=rhoSys->intensity(mywave1.c_str());
1281  if(maxIntens<myI)maxIntens=myI;
1282  if(minIntens>myI)minIntens=myI;
1283  if(iSys>0)delete rhoSys;
1284  } // end loop over systematic trees
1285 
1286  intenssysgraphs[iw]->SetPoint(i,mScaled,(maxIntens+minIntens)*0.5);
1287  intenssysgraphs[iw]->SetPointError(i,binwidth,(maxIntens-minIntens)*0.5);
1288  }
1289 
1290  //cout << "getting phases" << endl;
1291 
1292  // second loop to get phase differences
1293  unsigned int wi1=rho->waveIndex(wl[iw].c_str());
1294 
1295  for(unsigned int iw2=iw+1; iw2<wl.size();++iw2){
1296  //double ps2=rho->phaseSpaceIntegral(wl[iw2].c_str())*fsps;
1297 
1298  unsigned int wi2=rho->waveIndex(wl[iw2].c_str());
1299  complex<double> r=rho->spinDensityMatrixElem(wi1,wi2);
1300  TMatrixT<double> rCov=rho->spinDensityMatrixElemCov(wi1,wi2);
1301 
1302  realdatagraphs[c]->SetPoint(i,
1303  mScaled,
1304  r.real());
1305  realdatagraphs[c]->SetPointError(i,
1306  binwidth,
1307  sqrt(rCov[0][0]));
1308  imagdatagraphs[c]->SetPoint(i,
1309  mScaled,
1310  r.imag());
1311 
1312  imagdatagraphs[c]->SetPointError(i,
1313  binwidth,
1314  sqrt(rCov[1][1]));
1315 
1316 
1317  double dataphi=rho->phase(wl[iw].c_str(),wl[iw2].c_str());
1318 
1319  phasedatagraphs[c]->SetPoint(i,mScaled,dataphi);
1320 
1321  TVector2 v;v.SetMagPhi(1,rho->phase(wl[iw].c_str(),
1322  wl[iw2].c_str())/TMath::RadToDeg());
1323  phase2d[c]->SetPoint(i,v.X(),v.Y(),mScaled);
1324 
1325  phasedatagraphs[c]->SetPointError(i,binwidth,
1326  rho->phaseErr(wl[iw].c_str(),
1327  wl[iw2].c_str()));
1328  double fitphase=compset.phase(wl[iw],wl[iw2],m)*TMath::RadToDeg();
1329 
1330  if(sysPlotting){
1331  //cerr << "start sysplotting" << endl;
1332  // loop over systematics files
1333  double maxPhase=-10000;
1334  double minPhase=10000;
1335  double maxRe=-10000000;
1336  double minRe=10000000;
1337  double maxIm=-10000000;
1338  double minIm=10000000;
1339 
1340  for(unsigned int iSys=0;iSys<sysTrees.size();++iSys){
1341  //cerr << iSys;
1342  // get data
1343  fitResult* rhoSys=0;
1344  if(iSys==0)rhoSys=rho;
1345  else {
1346  sysTrees[iSys]->SetBranchAddress(valBranchName.c_str(),&rhoSys);
1347  if(i<sysTrees[iSys]->GetEntries()) sysTrees[iSys]->GetEntry(i);
1348  else{
1349  delete rhoSys;
1350  continue;
1351  }
1352  }
1353  // rename one wave
1354  string mywave1=wl[iw];
1355  string mywave2=wl[iw2];
1356 
1357 if(mywave1=="1-1++0+pi-_01_rho1700=pi-+_10_pi1300=pi+-_00_sigma.amp" && iSys>0)mywave1="1-1++0+pi-_01_eta11600=pi-+_10_pi1300=pi+-_00_sigma.amp";
1358 if(mywave2=="1-1++0+pi-_01_rho1700=pi-+_10_pi1300=pi+-_00_sigma.amp" && iSys>0)mywave2="1-1++0+pi-_01_eta11600=pi-+_10_pi1300=pi+-_00_sigma.amp";
1359 
1360  // check if waves are in fit
1361  if(rhoSys==NULL || rhoSys->waveIndex(mywave1.c_str())==-1 || rhoSys->waveIndex(mywave2.c_str()) ==-1){
1362  delete rhoSys;
1363  continue;
1364  }
1365  // get correct wave indices!!!
1366  unsigned int wi1Sys=rhoSys->waveIndex(mywave1.c_str());
1367  unsigned int wi2Sys=rhoSys->waveIndex(mywave2.c_str());
1368 
1369  double myphi=rhoSys->phase(mywave1.c_str(),mywave2.c_str());
1370  double myphiplus=myphi+360;
1371  double myphiminus=myphi-360;
1372  // translate by 2pi to get closest solution to fit
1373  if(fabs(myphiplus-dataphi)<fabs(myphi-dataphi)){
1374  myphi=myphiplus;
1375  //cout << "myphiminus" << endl;
1376  }
1377  if(fabs(myphiminus-dataphi)<fabs(myphi-dataphi)){
1378  myphi=myphiminus;
1379  //cout << "myphiplus" << endl;
1380  }
1381 
1382  if(myphi>maxPhase)maxPhase=myphi;
1383  if(myphi<minPhase)minPhase=myphi;
1384 
1385  // real and imag part:
1386  complex<double> r=rhoSys->spinDensityMatrixElem(wi1Sys,wi2Sys);
1387  if(maxRe<r.real())maxRe=r.real();
1388  if(minRe>r.real())minRe=r.real();
1389  if(maxIm<r.imag())maxIm=r.imag();
1390  if(minIm>r.imag())minIm=r.imag();
1391  if(iSys>0)delete rhoSys;
1392  rhoSys=0;
1393  }// end loop over sys trees
1394  // cerr << "loop over systrees finished" << endl;
1395 
1396  phasesysgraphs[c]->SetPoint(i,mScaled,(maxPhase+minPhase)*0.5);
1397  phasesysgraphs[c]->SetPointError(i,binwidth,(maxPhase-minPhase)*0.5);
1398 
1399  realsysgraphs[c]->SetPoint(i,mScaled,(maxRe+minRe)*0.5);
1400  realsysgraphs[c]->SetPointError(i,binwidth,(maxRe-minRe)*0.5);
1401  imagsysgraphs[c]->SetPoint(i,mScaled,(maxIm+minIm)*0.5);
1402  imagsysgraphs[c]->SetPointError(i,binwidth,(maxIm-minIm)*0.5);
1403 
1404 
1405  }// end if sysplotting
1406  //cerr << "sysplotting finished" << endl;
1407 
1408  phasefitgraphs[c]->SetPoint(i,mScaled,fitphase);
1409 
1410  complex<double> fitval=compset.overlap(wl[iw],wl[iw2],m);
1411 
1412  realfitgraphs[c]->SetPoint(i,mScaled,fitval.real());
1413  imagfitgraphs[c]->SetPoint(i,mScaled,fitval.imag());
1414 
1415 
1416  c++;
1417  }// end inner loop over waves
1418 
1419  prevps[iw]=ps;
1420 
1421 
1422  } // end loop over waves
1423  //cerr << "outer loop over waves finished" << endl;
1424  // loop over components to fill individual graphs
1425  unsigned int compcount=0;
1426  for(unsigned int ic=0;ic<compset.n();++ic){
1427  const pwacomponent* c=compset[ic];
1428  std::map<std::string,pwachannel >::const_iterator it=c->channels().begin();
1429  while(it!=c->channels().end()){
1430  double I=norm(c->val(m)*it->second.C())*it->second.ps(m)*compset.ps(m);
1431  compgraphs[compcount]->SetPoint(i,mScaled,I);
1432  ++compcount;
1433  ++it;
1434  } // end loop over channels
1435  }// end loop over components
1436 
1437  cerr << "Finished plotting mass-bin " << m << endl;
1438  //mprev=m;
1439  }
1440  cerr << "Finished Loop Over DataBins" << endl;
1441 
1442 
1443 
1444 // for(unsigned int im=0;im<nbins;++im){ // fine loop in masses -> fits
1445 
1446 // double m=mmin+im*md;
1447 // for(unsigned int iw=0; iw<wl.size();++iw){
1448 // fitgraphs[iw]->SetPoint(im,m,compset.intensity(wl[iw],m));
1449 // } // end loop over waves
1450 
1451 
1452 // }// end loop over mass bins
1453 
1454 
1455 
1456 
1457 
1458  TFile* outfile=TFile::Open(outFileName.c_str(),"RECREATE");
1459  for(unsigned int iw=0; iw<wl.size();++iw){
1460  TGraph* g=(TGraph*)graphs[iw]->GetListOfGraphs()->At(3);
1461  for(unsigned int ib=0;ib<nbins;++ib){
1462  double m,ps;g->GetPoint(ib,m,ps);
1463  g->SetPoint(ib,m,ps*500.);
1464  }
1465 
1466  graphs[iw]->Write();
1467  //absphasegraphs[iw]->Write();
1468  }
1469 
1470 
1471  for(unsigned int iw=0; iw<phasegraphs.size();++iw){
1472 
1474 
1475  unsigned int refbin=6;
1476  double m;
1477  double predph;
1478  phasedatagraphs[iw]->GetPoint(refbin,m,predph);
1479  double prefph;
1480  phasefitgraphs[iw]->GetPoint(refbin,m,prefph);
1481  for(unsigned int ib=refbin+1;ib<nbins;++ib){
1482  double dph; phasedatagraphs[iw]->GetPoint(ib,m,dph);
1483  double fph; phasefitgraphs[iw]->GetPoint(ib,m,fph);
1484  double dp,dm;dp=dph+360;dm=dph-360;
1485  double fp,fm;fp=fph+360;fm=fph-360;
1486  if(1){
1487  if(fabs(fp-prefph)<fabs(fph-prefph) && fabs(fp-prefph)<fabs(fm-prefph))
1488  phasefitgraphs[iw]->SetPoint(ib,m,fp);
1489  else if(fabs(fm-prefph)<fabs(fph-prefph) && fabs(fm-prefph)<fabs(fp-prefph))
1490  phasefitgraphs[iw]->SetPoint(ib,m,fm);
1491  phasefitgraphs[iw]->GetPoint(ib,m,prefph);
1492 
1493  if(fabs(dp-prefph)<fabs(dph-prefph) && fabs(dp-prefph)<fabs(dm-prefph))
1494  phasedatagraphs[iw]->SetPoint(ib,m,dp);
1495  else if(fabs(dm-prefph)<fabs(dph-prefph) && fabs(dm-prefph)<fabs(dp-prefph))
1496  phasedatagraphs[iw]->SetPoint(ib,m,dm);
1497 
1498 
1499 
1500  phasedatagraphs[iw]->GetPoint(ib,m,predph);
1501 
1502 
1503 
1504  // put systematic error closest to fit/data
1505  double sph; phasesysgraphs[iw]->GetPoint(ib,m,sph);
1506  double sp,sm;sp=sph+360;sm=sph-360;
1507  if(fabs(sp-prefph)<fabs(sph-prefph) && fabs(sp-prefph)<fabs(sm-prefph))
1508  phasesysgraphs[iw]->SetPoint(ib,m,sp);
1509  else if(fabs(sm-prefph)<fabs(sph-prefph) && fabs(sm-prefph)<fabs(sp-prefph))
1510  phasesysgraphs[iw]->SetPoint(ib,m,sm);
1511 
1512 
1513  }
1514  }
1515  // backward:
1516  phasedatagraphs[iw]->GetPoint(refbin,m,predph);
1517  phasefitgraphs[iw]->GetPoint(refbin,m,prefph);
1518  for(unsigned int i=0;i<refbin;++i){
1519  unsigned int ib=refbin-i-1;
1520  double dph; phasedatagraphs[iw]->GetPoint(ib,m,dph);
1521  double fph; phasefitgraphs[iw]->GetPoint(ib,m,fph);
1522  double dp,dm;dp=dph+360;dm=dph-360;
1523  double fp,fm;fp=fph+360;fm=fph-360;
1524  if(1){
1525 
1526  if(fabs(fp-prefph)<fabs(fph-prefph) && fabs(fp-prefph)<fabs(fm-prefph)){
1527  phasefitgraphs[iw]->SetPoint(ib,m,fp);
1528 
1529  }
1530  else if(fabs(fm-prefph)<fabs(fph-prefph) && fabs(fm-prefph)<fabs(fp-prefph)){
1531  phasefitgraphs[iw]->SetPoint(ib,m,fm);
1532 
1533  }
1534 
1535  phasefitgraphs[iw]->GetPoint(ib,m,prefph);
1536 
1537 
1538 
1539  if(fabs(dp-prefph)<fabs(dph-prefph) && fabs(dp-prefph)<fabs(dm-prefph)){
1540  phasedatagraphs[iw]->SetPoint(ib,m,dp);
1541  }
1542 
1543  else if(fabs(dm-prefph)<fabs(dph-prefph) && fabs(dm-prefph)<fabs(dp-prefph)){
1544  phasedatagraphs[iw]->SetPoint(ib,m,dm);
1545  }
1546 
1547  phasedatagraphs[iw]->GetPoint(ib,m,predph);
1548 
1549 
1550  // put systematic error closest to fit
1551  double sph; phasesysgraphs[iw]->GetPoint(ib,m,sph);
1552  double sp,sm;sp=sph+360;sm=sph-360;
1553  if(fabs(sp-predph)<fabs(sph-predph) && fabs(sp-predph)<fabs(sm-predph))
1554  phasesysgraphs[iw]->SetPoint(ib,m,sp);
1555  else if(fabs(sm-predph)<fabs(sph-predph) && fabs(sm-predph)<fabs(sp-predph))
1556  phasesysgraphs[iw]->SetPoint(ib,m,sm);
1557 
1558 
1559  }
1560  }
1561 
1562 
1563  phasegraphs[iw]->Write();
1564  phasesysgraphs[iw]->Write();
1565  overlapRegraphs[iw]->Write();
1566  overlapImgraphs[iw]->Write();
1567  //phase2d[iw]->Write();
1568  } // end loop over waves
1569 
1570  fPS->Write("fPS");
1571 
1572 outfile->Close();
1573 
1574  return 0;
1575 
1576 }
Generated by   doxygen 1.8.1.2