test/UnitTests/TestCrossSectionHandler.cc

0001 #include <cppunit/TestFixture.h>
0002 #include <cppunit/extensions/HelperMacros.h>
0003 #include <cppunit/TestResult.h>
0004 #include <cppunit/TestRunner.h>
0005 #include <cppunit/ui/text/TestRunner.h>
0006 #include <cppunit/TestResultCollector.h>
0007 #include <cppunit/TextTestProgressListener.h>
0008 #include <cppunit/CompilerOutputter.h>
0009
0010 #include <algorithm>
0011 #include <iterator>
0012
0013 #include "MuonAnalysis/MomentumScaleCalibration/interface/CrossSectionHandler.h"
0014
0015 #ifndef TestCrossSectionHandler_cc
0016 #define TestCrossSectionHandler_cc
0017
0018 class TestCrossSectionHandler : public CppUnit::TestFixture {
0019 public:
0020   TestCrossSectionHandler() {}
0021   void setUp() {
0022     crossSection.push_back(1.233);
0023     crossSection.push_back(2.07);
0024     crossSection.push_back(6.33);
0025     crossSection.push_back(13.9);
0026     crossSection.push_back(2.169);
0027     crossSection.push_back(127.2);
0028
0029     // First case: fit all the Upsilons
0030     fill_n(back_inserter(resfind_1), 6, 0);
0031     resfind_1[1] = 1;
0032     resfind_1[2] = 1;
0033     resfind_1[3] = 1;
0034     crossSectionHandler_1 = new CrossSectionHandler(crossSection, resfind_1);
0035
0036     // Second case: fit only the Upsilon(3S)
0037     fill_n(back_inserter(resfind_2), 6, 0);
0038     resfind_2[1] = 1;
0039     crossSectionHandler_2 = new CrossSectionHandler(crossSection, resfind_2);
0040
0041     // Third case: fit the Upsilon(3S) and the J/Psi and Psi(2S)
0042     fill_n(back_inserter(resfind_3), 6, 0);
0043     resfind_3[1] = 1;
0044     resfind_3[4] = 1;
0045     resfind_3[5] = 1;
0046     crossSectionHandler_3 = new CrossSectionHandler(crossSection, resfind_3);
0047
0048     // Fourth case: fit the Upsilon(3S) and (2S) and the J/Psi and Psi(2S)
0049     fill_n(back_inserter(resfind_4), 6, 0);
0050     resfind_4[1] = 1;
0051     resfind_4[2] = 1;
0052     resfind_4[4] = 1;
0053     resfind_4[5] = 1;
0054     crossSectionHandler_4 = new CrossSectionHandler(crossSection, resfind_4);
0055
0056     // Fifth case: fit nothing
0057     fill_n(back_inserter(resfind_5), 6, 0);
0058     crossSectionHandler_5 = new CrossSectionHandler(crossSection, resfind_5);
0059
0060     // Sixth case: fit everything
0061     fill_n(back_inserter(resfind_6), 6, 1);
0062     crossSectionHandler_6 = new CrossSectionHandler(crossSection, resfind_6);
0063   }
0064
0065   void tearDown() {
0066     delete crossSectionHandler_1;
0067     delete crossSectionHandler_2;
0068     delete crossSectionHandler_3;
0069     delete crossSectionHandler_4;
0070     delete crossSectionHandler_5;
0071     delete crossSectionHandler_6;
0072   }
0073
0074   void testConstructor() {
0075     CPPUNIT_ASSERT(crossSectionHandler_1->parNum_ == 2);
0076     CPPUNIT_ASSERT(crossSectionHandler_2->parNum_ == 0);
0077     CPPUNIT_ASSERT(crossSectionHandler_3->parNum_ == 2);
0078     CPPUNIT_ASSERT(crossSectionHandler_4->parNum_ == 3);
0079     CPPUNIT_ASSERT(crossSectionHandler_5->parNum_ == 0);
0080     CPPUNIT_ASSERT(crossSectionHandler_6->parNum_ == 5);
0081     CPPUNIT_ASSERT(crossSectionHandler_1->vars_.size() == crossSectionHandler_1->parNum_);
0082     CPPUNIT_ASSERT(crossSectionHandler_2->vars_.size() == crossSectionHandler_2->parNum_);
0083     CPPUNIT_ASSERT(crossSectionHandler_3->vars_.size() == crossSectionHandler_3->parNum_);
0084     CPPUNIT_ASSERT(crossSectionHandler_4->vars_.size() == crossSectionHandler_4->parNum_);
0085     CPPUNIT_ASSERT(crossSectionHandler_5->vars_.size() == crossSectionHandler_5->parNum_);
0086     CPPUNIT_ASSERT(crossSectionHandler_6->vars_.size() == crossSectionHandler_6->parNum_);
0087   }
0088
0089   void testComputeRelativeCrossSections() {
0090     crossSectionHandler_1->computeRelativeCrossSections(crossSection, resfind_1);
0091     crossSectionHandler_2->computeRelativeCrossSections(crossSection, resfind_2);
0092     crossSectionHandler_3->computeRelativeCrossSections(crossSection, resfind_3);
0093     crossSectionHandler_4->computeRelativeCrossSections(crossSection, resfind_4);
0094     crossSectionHandler_5->computeRelativeCrossSections(crossSection, resfind_5);
0095     crossSectionHandler_6->computeRelativeCrossSections(crossSection, resfind_6);
0096
0097     std::vector<std::vector<double> > relativeCrossSections;
0098     relativeCrossSections.push_back(expandCrossSectionVec(crossSectionHandler_1->relativeCrossSectionVec_, resfind_1));
0099     relativeCrossSections.push_back(expandCrossSectionVec(crossSectionHandler_2->relativeCrossSectionVec_, resfind_2));
0100     relativeCrossSections.push_back(expandCrossSectionVec(crossSectionHandler_3->relativeCrossSectionVec_, resfind_3));
0101     relativeCrossSections.push_back(expandCrossSectionVec(crossSectionHandler_4->relativeCrossSectionVec_, resfind_4));
0102     relativeCrossSections.push_back(expandCrossSectionVec(crossSectionHandler_5->relativeCrossSectionVec_, resfind_5));
0103     relativeCrossSections.push_back(expandCrossSectionVec(crossSectionHandler_6->relativeCrossSectionVec_, resfind_6));
0104
0105     checkRelativeCrossSections(relativeCrossSections);
0106   }
0107
0108   void checkRelativeCrossSections(const std::vector<std::vector<double> >& relativeCrossSections) {
0109     // First case: fit all the Upsilons
0110     double norm = crossSection[1] + crossSection[2] + crossSection[3];
0111     CPPUNIT_ASSERT(relativeCrossSections[0].size() == 6);
0112     CPPUNIT_ASSERT(float(relativeCrossSections[0][0]) == float(0.));
0113     CPPUNIT_ASSERT(float(relativeCrossSections[0][1]) == float(crossSection[1] / norm));
0114     CPPUNIT_ASSERT(float(relativeCrossSections[0][2]) == float(crossSection[2] / norm));
0115     CPPUNIT_ASSERT(float(relativeCrossSections[0][3]) == float(crossSection[3] / norm));
0116     CPPUNIT_ASSERT(float(relativeCrossSections[0][4]) == float(0.));
0117     CPPUNIT_ASSERT(float(relativeCrossSections[0][5]) == float(0.));
0118
0119     // Second case: fit only the Upsilon(3S)
0120     CPPUNIT_ASSERT(relativeCrossSections[1].size() == 6);
0121     CPPUNIT_ASSERT(float(relativeCrossSections[1][0]) == float(0.));
0122     CPPUNIT_ASSERT(float(relativeCrossSections[1][1]) == float(1.));
0123     CPPUNIT_ASSERT(float(relativeCrossSections[1][2]) == float(0.));
0124     CPPUNIT_ASSERT(float(relativeCrossSections[1][3]) == float(0.));
0125     CPPUNIT_ASSERT(float(relativeCrossSections[1][4]) == float(0.));
0126     CPPUNIT_ASSERT(float(relativeCrossSections[1][5]) == float(0.));
0127
0128     // Third case: fit the Upsilon(3S) and the J/Psi and Psi(2S)
0129     CPPUNIT_ASSERT(relativeCrossSections[2].size() == 6);
0130     norm = crossSection[1] + crossSection[4] + crossSection[5];
0131     CPPUNIT_ASSERT(float(relativeCrossSections[2][0]) == float(0.));
0132     CPPUNIT_ASSERT(float(relativeCrossSections[2][1]) == float(crossSection[1] / norm));
0133     CPPUNIT_ASSERT(float(relativeCrossSections[2][2]) == float(0.));
0134     CPPUNIT_ASSERT(float(relativeCrossSections[2][3]) == float(0.));
0135     CPPUNIT_ASSERT(float(relativeCrossSections[2][4]) == float(crossSection[4] / norm));
0136     CPPUNIT_ASSERT(float(relativeCrossSections[2][5]) == float(crossSection[5] / norm));
0137
0138     // Fourth case: fit the Upsilon(3S) and (2S) and the J/Psi and Psi(2S)
0139     CPPUNIT_ASSERT(relativeCrossSections[3].size() == 6);
0140     norm = crossSection[1] + crossSection[2] + crossSection[4] + crossSection[5];
0141     CPPUNIT_ASSERT(float(relativeCrossSections[3][0]) == float(0.));
0142     CPPUNIT_ASSERT(float(relativeCrossSections[3][1]) == float(crossSection[1] / norm));
0143     CPPUNIT_ASSERT(float(relativeCrossSections[3][2]) == float(crossSection[2] / norm));
0144     CPPUNIT_ASSERT(float(relativeCrossSections[3][3]) == float(0.));
0145     CPPUNIT_ASSERT(float(relativeCrossSections[3][4]) == float(crossSection[4] / norm));
0146     CPPUNIT_ASSERT(float(relativeCrossSections[3][5]) == float(crossSection[5] / norm));
0147
0148     // Fifth case: fit nothing
0149     CPPUNIT_ASSERT(relativeCrossSections[4].size() == 6);
0150     CPPUNIT_ASSERT(float(relativeCrossSections[4][0]) == float(0.));
0151     CPPUNIT_ASSERT(float(relativeCrossSections[4][1]) == float(0.));
0152     CPPUNIT_ASSERT(float(relativeCrossSections[4][2]) == float(0.));
0153     CPPUNIT_ASSERT(float(relativeCrossSections[4][3]) == float(0.));
0154     CPPUNIT_ASSERT(float(relativeCrossSections[4][4]) == float(0.));
0155     CPPUNIT_ASSERT(float(relativeCrossSections[4][5]) == float(0.));
0156
0157     // Sixth case: fit everything
0158     CPPUNIT_ASSERT(relativeCrossSections[5].size() == 6);
0159     norm = crossSection[0] + crossSection[1] + crossSection[2] + crossSection[3] + crossSection[4] + crossSection[5];
0160     CPPUNIT_ASSERT(float(relativeCrossSections[5][0]) == float(crossSection[0] / norm));
0161     CPPUNIT_ASSERT(float(relativeCrossSections[5][1]) == float(crossSection[1] / norm));
0162     CPPUNIT_ASSERT(float(relativeCrossSections[5][2]) == float(crossSection[2] / norm));
0163     CPPUNIT_ASSERT(float(relativeCrossSections[5][3]) == float(crossSection[3] / norm));
0164     CPPUNIT_ASSERT(float(relativeCrossSections[5][4]) == float(crossSection[4] / norm));
0165     CPPUNIT_ASSERT(float(relativeCrossSections[5][5]) == float(crossSection[5] / norm));
0166   }
0167
0168   void testImposeConstraint() {
0169     // First case: fit all the Upsilons
0170     crossSectionHandler_1->computeRelativeCrossSections(crossSection, resfind_1);
0171     crossSectionHandler_1->imposeConstraint();
0172     CPPUNIT_ASSERT(float(crossSectionHandler_1->vars_[0]) == float(crossSection[2] / crossSection[1]));
0173     CPPUNIT_ASSERT(float(crossSectionHandler_1->vars_[1]) == float(crossSection[3] / crossSection[2]));
0174
0175     // Second case: fit only the Upsilon(3S). Nothing to test. The vars_ vector is empty (tested in the constructor test)
0176
0177     // Third case: fit the Upsilon(3S) and the J/Psi and Psi(2S)
0178     crossSectionHandler_3->computeRelativeCrossSections(crossSection, resfind_3);
0179     crossSectionHandler_3->imposeConstraint();
0180     CPPUNIT_ASSERT(float(crossSectionHandler_3->vars_[0]) == float(crossSection[4] / crossSection[1]));
0181     CPPUNIT_ASSERT(float(crossSectionHandler_3->vars_[1]) == float(crossSection[5] / crossSection[4]));
0182
0183     // Fourth case: fit the Upsilon(3S) and (2S) and the J/Psi and Psi(2S)
0184     crossSectionHandler_4->computeRelativeCrossSections(crossSection, resfind_4);
0185     crossSectionHandler_4->imposeConstraint();
0186     CPPUNIT_ASSERT(float(crossSectionHandler_4->vars_[0]) == float(crossSection[2] / crossSection[1]));
0187     CPPUNIT_ASSERT(float(crossSectionHandler_4->vars_[1]) == float(crossSection[4] / crossSection[2]));
0188     CPPUNIT_ASSERT(float(crossSectionHandler_4->vars_[2]) == float(crossSection[5] / crossSection[4]));
0189
0190     // Fifth case: fit nothing. vars_.size() = 0 tested in the constructor test.
0191
0192     // Sixth case: fit everything
0193     crossSectionHandler_6->computeRelativeCrossSections(crossSection, resfind_6);
0194     crossSectionHandler_6->imposeConstraint();
0195     CPPUNIT_ASSERT(float(crossSectionHandler_6->vars_[0]) == float(crossSection[1] / crossSection[0]));
0196     CPPUNIT_ASSERT(float(crossSectionHandler_6->vars_[1]) == float(crossSection[2] / crossSection[1]));
0197     CPPUNIT_ASSERT(float(crossSectionHandler_6->vars_[2]) == float(crossSection[3] / crossSection[2]));
0198     CPPUNIT_ASSERT(float(crossSectionHandler_6->vars_[3]) == float(crossSection[4] / crossSection[3]));
0199     CPPUNIT_ASSERT(float(crossSectionHandler_6->vars_[4]) == float(crossSection[5] / crossSection[4]));
0200   }
0201
0202   void testRelativeCrossSections() {
0203     std::vector<std::vector<double> > relativeCrossSections;
0204     relativeCrossSections.push_back(getRelativeCrossSections(crossSectionHandler_1, resfind_1));
0205     relativeCrossSections.push_back(getRelativeCrossSections(crossSectionHandler_2, resfind_2));
0206     relativeCrossSections.push_back(getRelativeCrossSections(crossSectionHandler_3, resfind_3));
0207     relativeCrossSections.push_back(getRelativeCrossSections(crossSectionHandler_4, resfind_4));
0208     relativeCrossSections.push_back(getRelativeCrossSections(crossSectionHandler_5, resfind_5));
0209     relativeCrossSections.push_back(getRelativeCrossSections(crossSectionHandler_6, resfind_6));
0210
0211     checkRelativeCrossSections(relativeCrossSections);
0212   }
0213
0214   std::vector<double> getRelativeCrossSections(CrossSectionHandler* crossSectionHandler,
0215                                                const std::vector<int> resfind) {
0216     crossSectionHandler->computeRelativeCrossSections(crossSection, resfind);
0217     crossSectionHandler->imposeConstraint();
0218     std::vector<double>::const_iterator it = crossSectionHandler->vars_.begin();
0219     double* variables = new double[crossSectionHandler->vars_.size()];
0220     unsigned int i = 0;
0221     for (; it != crossSectionHandler->vars_.end(); ++it, ++i) {
0222       variables[i] = *it;
0223     }
0224     std::vector<double> vars(crossSectionHandler->relativeCrossSections(variables, resfind));
0225
0226     delete[] variables;
0227     return vars;
0228   }
0229
0230   void testSetParameters() {
0231     TMinuit rmin(crossSectionHandler_1->parNum_);
0232     std::vector<int> crossSectionOrder(6, 0);
0233     crossSectionOrder[1] = 1;
0234
0235     crossSectionHandler_1->setParameters(
0236         Start, Step, Mini, Maxi, ind, parname, crossSection, crossSectionOrder, resfind_1);
0237     CPPUNIT_ASSERT(Start[0] == crossSectionHandler_1->vars_[0]);
0238     CPPUNIT_ASSERT(Start[1] == crossSectionHandler_1->vars_[1]);
0239     CPPUNIT_ASSERT(parname[0] == "cross section var 1");
0240     CPPUNIT_ASSERT(parname[1] == "cross section var 2");
0241     CPPUNIT_ASSERT(ind[0] == 0);
0242     CPPUNIT_ASSERT(ind[1] == 1);
0243   }
0244
0245   void testReleaseParameters() {
0246     TMinuit rmin(crossSectionHandler_1->parNum_);
0247     std::vector<int> crossSectionOrder(6, 0);
0248     crossSectionOrder[1] = 1;
0249     crossSectionHandler_1->setParameters(
0250         Start, Step, Mini, Maxi, ind, parname, crossSection, crossSectionOrder, resfind_1);
0251     int ierror;
0252     for (unsigned int ipar = 0; ipar < crossSectionHandler_1->parNum(); ++ipar) {
0253       rmin.mnparm(ipar, parname[ipar], Start[ipar], Step[ipar], Mini[ipar], Maxi[ipar], ierror);
0254       rmin.FixParameter(ipar);
0255     }
0256
0257     std::vector<int> parfix(6, 0);
0258     // Test by setting two parameters to order 1, but one is for the Z which is not fitted and will be ignored.
0259     int* ind = new int[6];
0260     ind[0] = 0;
0261     ind[1] = 1;
0262     ind[2] = 0;
0263     ind[3] = 0;
0264     ind[4] = 0;
0265     ind[5] = 1;
0266     crossSectionHandler_1->releaseParameters(rmin, resfind_1, parfix, ind, 0, 0);
0267     CPPUNIT_ASSERT(rmin.GetNumFixedPars() == 1);
0268     crossSectionHandler_1->releaseParameters(rmin, resfind_1, parfix, ind, 1, 0);
0269     CPPUNIT_ASSERT(rmin.GetNumFixedPars() == 0);
0270   }
0271
0272   std::vector<double> expandCrossSectionVec(const std::vector<double>& relativeCrossSectionVec,
0273                                             const std::vector<int>& resfind) {
0274     std::vector<double> relCrossSec;
0275     unsigned int smallerVectorIndex = 0;
0276     std::vector<int>::const_iterator it = resfind.begin();
0277     for (; it != resfind.end(); ++it) {
0278       if (*it == 0) {
0279         relCrossSec.push_back(0.);
0280       } else {
0281         relCrossSec.push_back(relativeCrossSectionVec[smallerVectorIndex]);
0282         ++smallerVectorIndex;
0283       }
0284     }
0285     return relCrossSec;
0286   }
0287
0288   // Data members
0289   CrossSectionHandler* crossSectionHandler_1;
0290   CrossSectionHandler* crossSectionHandler_2;
0291   CrossSectionHandler* crossSectionHandler_3;
0292   CrossSectionHandler* crossSectionHandler_4;
0293   CrossSectionHandler* crossSectionHandler_5;
0294   CrossSectionHandler* crossSectionHandler_6;
0295   std::vector<int> resfind_1;
0296   std::vector<int> resfind_2;
0297   std::vector<int> resfind_3;
0298   std::vector<int> resfind_4;
0299   std::vector<int> resfind_5;
0300   std::vector<int> resfind_6;
0301   std::vector<double> crossSection;
0302
0303   double Start[100];
0304   double Step[100];
0305   double Mini[100];
0306   double Maxi[100];
0307   int ind[100];
0308   TString parname[100];
0309
0310   // Declare and build the test suite
0311   CPPUNIT_TEST_SUITE(TestCrossSectionHandler);
0312   CPPUNIT_TEST(testConstructor);
0313   CPPUNIT_TEST(testComputeRelativeCrossSections);
0314   CPPUNIT_TEST(testImposeConstraint);
0315   CPPUNIT_TEST(testRelativeCrossSections);
0316   CPPUNIT_TEST(testSetParameters);
0317   CPPUNIT_TEST(testReleaseParameters);
0318   CPPUNIT_TEST_SUITE_END();
0319 };
0320
0321 // Register the test suite in the registry.
0322 // This way we will have to only pass the registry to the runner
0323 // and it will contain all the registered test suites.
0324 CPPUNIT_TEST_SUITE_REGISTRATION(TestCrossSectionHandler);
0325
0326 #endif