ROOTPWA
testIntegral.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2010
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 // $Rev:: $: revision of last commit
24 // $Author:: $: author of last commit
25 // $Date:: $: date of last commit
26 //
27 // Description:
28 // basic test program for amplitude integral matrix
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <algorithm>
39 
40 #include <boost/assign/list_of.hpp>
41 
42 #include "TROOT.h"
43 #include "TFile.h"
44 
45 #include "reportingUtils.hpp"
46 #include "ampIntegralMatrix.h"
47 
48 
49 using namespace std;
50 using namespace boost::assign;
51 using namespace rpwa;
52 
53 
54 int
55 main(int argc,
56  char** argv)
57 {
58  // switch on debug output
59  ampIntegralMatrix::setDebug(true);
60 
61  if (0) {
62  // test integral calculation
63  // get file list
64  // vector<string> binAmpFileNames = filesMatchingGlobPattern("/data/compass/hadronData/massBins/2004/Q3PiData/r481.trunk/1260.1300/PSPAMPS/SYM/*.amp", true);
65  // sort(binAmpFileNames.begin(), binAmpFileNames.end());
66  // make sure order of waves is the same as used by int
67  vector<string> binAmpFileNames = list_of
68  ("1-0-+0+f0980_00_pi-.amp" )
69  ("1-0-+0+rho770_11_pi-.amp" )
70  ("1-0-+0+sigma_00_pi-.amp" )
71  ("1-1++0+f21270_12_pi-.amp" )
72  ("1-1++0+rho770_01_pi-.amp" )
73  ("1-1-+0-rho770_11_pi-.amp" )
74  ("1-1++0+rho770_21_pi-.amp" )
75  ("1-1++0+sigma_10_pi-.amp" )
76  ("1-1++1+f21270_12_pi-.amp" )
77  ("1-1++1-rho770_01_pi-.amp" )
78  ("1-1++1+rho770_01_pi-.amp" )
79  ("1-1-+1-rho770_11_pi-.amp" )
80  ("1-1-+1+rho770_11_pi-.amp" )
81  ("1-1++1+rho770_21_pi-.amp" )
82  ("1-1++1+sigma_10_pi-.amp" )
83  ("1-2-+0+f21270_02_pi-.amp" )
84  ("1-2++0-f21270_12_pi-.amp" )
85  ("1-2-+0+f21270_22_pi-.amp" )
86  ("1-2-+0+rho770_11_pi-.amp" )
87  ("1-2++0-rho770_21_pi-.amp" )
88  ("1-2-+0+rho770_31_pi-.amp" )
89  ("1-2-+0+sigma_20_pi-.amp" )
90  ("1-2-+1-f21270_02_pi-.amp" )
91  ("1-2-+1+f21270_02_pi-.amp" )
92  ("1-2++1-f21270_12_pi-.amp" )
93  ("1-2++1+f21270_12_pi-.amp" )
94  ("1-2-+1+f21270_22_pi-.amp" )
95  ("1-2-+1+rho770_11_pi-.amp" )
96  ("1-2++1+rho770_21_pi-.amp" )
97  ("1-2-+1+rho770_31_pi-.amp" )
98  ("1-2-+1+sigma_20_pi-.amp" )
99  ("1-3++0+f21270_12_pi-.amp" )
100  ("1-3++0+rho31690_03_pi-.amp")
101  ("1-3++0+rho770_21_pi-.amp" )
102  ("1-3++1+f21270_12_pi-.amp" )
103  ("1-3++1+rho31690_03_pi-.amp")
104  ("1-3++1+rho770_21_pi-.amp" )
105  ("1-4-+0+rho770_31_pi-.amp" )
106  ("1-4++1+f21270_32_pi-.amp" )
107  ("1-4-+1+rho770_31_pi-.amp" )
108  ("1-4++1+rho770_41_pi-.amp" );
109  for (size_t i = 0; i < binAmpFileNames.size(); ++i)
110  binAmpFileNames[i] = "/data/compass/hadronData/massBins/2004/Q3PiData/r481.trunk/1260.1300/PSPAMPS/SYM/" + binAmpFileNames[i];
111  vector<string> rootAmpFileNames;
112  ampIntegralMatrix integral;
113  integral.integrate(binAmpFileNames, rootAmpFileNames);
114  integral.writeAscii("testIntegral2.int");
115  }
116 
117 
118  if (1) {
119  // test I/O and copying
120  ampIntegralMatrix integral;
121  // ascii I/O
122  integral.readAscii("testIntegral.int");
123  ampIntegralMatrix integral2(integral);
124  integral2.writeAscii("testIntegral2.int");
125 #ifdef USE_STD_COMPLEX_TREE_LEAFS
126  // root I/O
127  // force loading predefined std::complex dictionary
128  gROOT->ProcessLine("#include <complex>");
129  {
130  TFile* outFile = TFile::Open("testIntegral.root", "RECREATE");
131  printInfo << "writing integral to 'testIntegral.root'" << endl;
132  integral.Write("integral");
133  outFile->Close();
134  }
135  {
136  TFile* inFile = TFile::Open("testIntegral.root", "READ");
137  ampIntegralMatrix* integral3 = 0;
138  printInfo << "reading integral from 'testIntegral.root'" << endl;
139  inFile->GetObject("integral", integral3);
140  if (not integral3)
141  printErr << "cannot find integral 'integral'" << endl;
142  else {
143  if (*integral3 == integral)
144  printSucc << "ROOT file integral is identical" << endl;
145  else
146  printErr << "ROOT file intergral differs" << endl;
147  integral3->writeAscii("testIntegral3.int");
148  }
149  inFile->Close();
150  }
151 #endif // USE_STD_COMPLEX_TREE_LEAFS
152  }
153 
154 
155  if (0) {
156  // test arthmetic operations
157  ampIntegralMatrix integral, integrals[6];
158  const unsigned int nmbInt = sizeof(integrals) / sizeof(integrals[0]);
159  integral.readAscii("testIntegral.int");
160  for (unsigned int i = 0; i < nmbInt; ++i)
161  integrals[i].readAscii("testIntegral.int");
162  integrals[0] = integral + integral;
163  integrals[1] = integrals[0] - integral;
164  integrals[2] = 2 * integral - integral;
165  integrals[3] = integral * 0.1 + integral * 0.9;
166  integrals[4] = integral / 2 + integral / 2;
167  integrals[5] = (integral + 2 * integral) / 3;
168  integrals[5].writeAscii("testIntegral1.int");
169  for (unsigned int i = 1; i < nmbInt; ++i)
170  if (integrals[i] == integral)
171  printSucc << "integrals[" << i << "] is identical" << endl;
172  else
173  printErr << "integrals[" << i << "] differs" << endl;
174  }
175 
176 
177  if (0) {
178  // test renormalization
179  ampIntegralMatrix integral;
180  integral.readAscii("testIntegral.int");
181  const unsigned int nmbEvents = integral.nmbEvents();
182  integral.renormalize(10000);
183  integral.renormalize(nmbEvents); // scale back to original
184  integral.writeAscii("testIntegral2.int");
185  cout << integral;
186  //integral.print(cout, true);
187  }
188 
189 
190 }
Generated by   doxygen 1.8.1.2