CocoaFit/src/NtupleManager.cc

0001 //   COCOA class implementation file
0002 //Id:  NtupleManager.cc
0003 //CAT: Analysis
0004 //
0005 //   History: v1.0
0006 //   Luca Scodellaro
0007 #include "TROOT.h"
0008 #include <cstdlib>
0009 #include <fstream>
0010
0011 #include "Alignment/CocoaModel/interface/Model.h"
0012 #include "Alignment/CocoaFit/interface/NtupleManager.h"
0013 #include "Alignment/CocoaFit/interface/FittedEntry.h"
0014 #include "Alignment/CocoaModel/interface/Measurement.h"
0015 #include "Alignment/CocoaModel/interface/OpticalObject.h"
0016 #include "Alignment/CocoaModel/interface/Entry.h"
0017 // #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
0018 // #include "Alignment/CocoaUtilities/interface/GlobalOptionMgr.h"
0019 #include "TFile.h"
0020 #include "TTree.h"
0021 #include "TClonesArray.h"
0022 //#include "TMath.h"
0023
0024 NtupleManager* NtupleManager::instance = nullptr;
0025
0026 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0027 //@@  Gets the only instance of Model
0028 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0029 NtupleManager* NtupleManager::getInstance() {
0030   if (!instance) {
0031     instance = new NtupleManager;
0032   }
0033   return instance;
0034 }
0035
0036 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0037 //@@ Book ntuple
0038 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0039 void NtupleManager::BookNtuple() {
0040   theRootFile = new TFile("report.root", "RECREATE", "Simple ROOT Ntuple");
0041
0042   CocoaTree = new TTree("CocoaTree", "CocoaTree");
0043
0044   CocoaTree->Branch("Chi2Measurements", &Chi2Measurements, "Chi2Measurements/D");
0045   CocoaTree->Branch("Chi2CalibratedParameters", &Chi2CalibratedParameters, "Chi2CalibratedParameters/D");
0046   CocoaTree->Branch("NDegreesOfFreedom", &NDegreesOfFreedom, "NDegreesOfFreedom/I");
0047
0048   CloneFitParam = new TClonesArray("FitParam");
0049   CocoaTree->Branch("FitParameters", &CloneFitParam, 32000, 2);
0050   CocoaTree->Branch("NFitParameters", &NFitParameters, "NFitParameters/I");
0051
0052   CloneOptObject = new TClonesArray("OptObject");
0053   CocoaTree->Branch("OptObjects", &CloneOptObject, 32000, 2);
0054   CocoaTree->Branch("NOptObjects", &NOptObjects, "NOptObjects/I");
0055
0056   CloneSensor2DMeas = new TClonesArray("Sensor2DMeas");
0057   CocoaTree->Branch("Sensor2DMeasurements", &CloneSensor2DMeas, 32000, 2);
0058   CocoaTree->Branch("NSensor2D", &NSensor2D, "NSensor2D/I");
0059
0060   CloneDistancemeterMeas = new TClonesArray("DistancemeterMeas");
0061   CocoaTree->Branch("DistancemeterMeasurements", &CloneDistancemeterMeas, 32000, 2);
0062   CocoaTree->Branch("NDistancemeter", &NDistancemeter, "NDistancemeter/I");
0063
0064   CloneDistancemeter1DimMeas = new TClonesArray("Distancemeter1DimMeas");
0065   CocoaTree->Branch("Distancemeter1DimMeasurements", &CloneDistancemeter1DimMeas, 32000, 2);
0066   CocoaTree->Branch("NDistancemeter1Dim", &NDistancemeter1Dim, "NDistancemeter1Dim/I");
0067
0068   CloneTiltmeterMeas = new TClonesArray("TiltmeterMeas");
0069   CocoaTree->Branch("TiltmeterMeasurements", &CloneTiltmeterMeas, 32000, 2);
0070   CocoaTree->Branch("NTiltmeter", &NTiltmeter, "NTiltmeter/I");
0071
0072   CloneCopsMeas = new TClonesArray("CopsMeas");
0073   CocoaTree->Branch("CopsMeasurements", &CloneCopsMeas, 32000, 2);
0074   CocoaTree->Branch("NCops", &NCops, "NCops/I");
0075
0076   theRootFile->Add(CocoaTree);
0077
0078   //   FitParametersTree = new TTree("FitParametersTree","FitParametersTree");
0079   //   FitParametersTree->Branch("NFitParameters",&NFitParameters,"NFitParameters/I");
0080   //   BookFitParameters = false;
0081   //   theRootFile->Add(FitParametersTree);
0082
0083   //   MeasurementsTree = new TTree("MeasurementsTree","MeasurementsTree");
0084   //   MeasurementsTree->Branch("NMeasurements",&NMeasurements,"NMeasurements/I");
0085   //   BookMeasurements = false;
0086   //   theRootFile->Add(MeasurementsTree);
0087 }
0088 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0089 //@@ Init ntuple
0090 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0091 void NtupleManager::InitNtuple() {
0092   CloneFitParam->Clear();
0093
0094   Chi2Measurements = 0.;
0095   Chi2CalibratedParameters = 0.;
0096   NDegreesOfFreedom = 0;
0097   NFitParameters = 0;
0098   NOptObjects = 0;
0099   NSensor2D = 0;
0100   NDistancemeter = 0;
0101   NDistancemeter1Dim = 0;
0102   NTiltmeter = 0;
0103   NCops = 0;
0104 }
0105 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0106 //@@ Fill ntuple tree
0107 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0108 void NtupleManager::FillNtupleTree() { CocoaTree->Fill(); }
0109 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0110 //@@ Close ntuple
0111 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0112 void NtupleManager::WriteNtuple() {
0113   theRootFile->Write();
0114   theRootFile->Close();
0115 }
0116 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0117 //@@ Close ntuple
0118 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0119 void NtupleManager::FillChi2() {
0120   double chi2meas = 0;
0121   double chi2cal = 0;
0122   ALIint nMeas = 0, nUnk = 0;
0123
0124   //----- Calculate the chi2 of measurements
0125   std::vector<Measurement*>::const_iterator vmcite;
0126   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
0127     for (ALIuint ii = 0; ii < ALIuint((*vmcite)->dim()); ii++) {
0128       nMeas++;
0129       double c2 = ((*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii)) / (*vmcite)->sigma(ii);
0130       chi2meas += c2 * c2;
0131     }
0132   }
0133
0134   //----- Calculate the chi2 of calibrated parameters
0135   std::vector<Entry*>::iterator veite;
0136   for (veite = Model::EntryList().begin(); veite != Model::EntryList().end(); ++veite) {
0137     if ((*veite)->quality() == 2)
0138       nUnk++;
0139     if ((*veite)->quality() == 1) {
0140       //       std::cout << " " << (*veite)->valueDisplacementByFitting() << " "
0141       //        << (*veite)->value << " " << (*veite)->sigma() << std::endl;
0142       double c2 = (*veite)->valueDisplacementByFitting() / (*veite)->sigma();
0143       chi2cal += c2 * c2;
0144     }
0145   }
0146
0147   Chi2Measurements = chi2meas;
0148   Chi2CalibratedParameters = chi2cal;
0149   NDegreesOfFreedom = nMeas - nUnk;
0150 }
0151 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0152 //@@ Fill ntuple with fitted parameters
0153 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0154 void NtupleManager::FillFitParameters(MatrixMeschach* AtWAMatrix) {
0155   //   double ParValue[1000], ParError[1000];
0156   int theMinEntryQuality = 1;
0157   int ii = 0;
0158   std::vector<Entry*>::const_iterator vecite;
0159   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
0160     //--- Only for good quality parameters (='unk')
0161     if ((*vecite)->quality() >= theMinEntryQuality) {
0162       ALIint ipos = (*vecite)->fitPos();
0163       FittedEntry* fe = new FittedEntry((*vecite), ipos, sqrt(AtWAMatrix->Mat()->me[ipos][ipos]));
0164       //       if (!BookFitParameters) {
0165       //    CocoaTree->Branch("NFitParameters",&NFitParameters,"NFitParameters/I:");
0166       //    ALIstring partype = fe->getName() + "/D";
0167       //    FitParametersTree->Branch(fe->getName().c_str(), &ParValue[ii], partype.c_str());
0168       //    ALIstring parerrname = fe->getName() + "_err";
0169       //    ALIstring parerrtype = parerrname + "/D";
0170       //    FitParametersTree->Branch(parerrname.c_str(), &ParError[ii], parerrtype.c_str());
0171       //       }
0172       //       ParValue[ii] = fe->getValue();
0173       //       ParError[ii] = fe->getSigma();
0174       std::cout << "EEE " << (*vecite)->ValueDimensionFactor() << " " << (*vecite)->SigmaDimensionFactor() << " "
0175                 << fe->getOptOName() << " " << fe->getEntryName() << " " << fe->getName() << " " << fe->getOrder()
0176                 << " " << fe->getQuality() << " " << (*vecite)->type() << " " << std::endl;
0177       FitParamA = new ((*CloneFitParam)[ii]) FitParam();
0178       FitParamA->Name = fe->getName();
0179       if (fe->getQuality() == 1)
0180         FitParamA->Quality = "Calibrated";
0181       else if (fe->getQuality() == 2)
0182         FitParamA->Quality = "Unknown";
0183       for (int no = 0; no < NOptObjects; no++) {
0184         OptObject* optobj = (OptObject*)CloneOptObject->At(no);
0185         if (optobj->Name == fe->getOptOName())
0186           FitParamA->OptObjectIndex = no;
0187       }
0188       float DF = 1.;
0189       if ((*vecite)->type() == "centre" || (*vecite)->type() == "length")
0190         DF = 1000.;
0191       FitParamA->InitialValue = DF * fe->getOrigValue() * (*vecite)->ValueDimensionFactor();
0192       FitParamA->InitialSigma = DF * fe->getOrigSigma() * (*vecite)->SigmaDimensionFactor();
0193       FitParamA->FittedValue = DF * fe->getValue() * (*vecite)->ValueDimensionFactor();
0194       FitParamA->FittedSigma = DF * fe->getSigma() * (*vecite)->SigmaDimensionFactor();
0195       ii++;
0196     }
0197   }
0198   //   BookFitParameters = true;
0199   NFitParameters = ii;
0200   //   FitParametersTree->Fill();
0201
0202   /*
0203   //---------- Loop sets of entries
0204   std::vector< FittedEntriesSet* > theFittedEntriesSets;
0205   std::vector< FittedEntriesSet* >::const_iterator vfescite;
0206   std::vector< FittedEntry* >::const_iterator vfecite;
0207   ALIint jj = 1;
0208   for( vfescite = theFittedEntriesSets.begin(); vfescite != theFittedEntriesSets.end(); vfescite++) {
0209     //---------- Loop entries
0210     if( vfescite == theFittedEntriesSets.begin() ) {
0211     //----- dump entries names if first set
0212       ALIint ii = 0;
0213       for( vfecite = ((*vfescite)->FittedEntries()).begin(); vfecite != ((*vfescite)->FittedEntries()).end(); vfecite++) {
0214     ALIstring partype = (*vfecite)->getName() + "/D";
0215     FitParametersTree->Branch((*vfecite)->getName().c_str(), &ParValue[ii], partype.c_str());
0216     ALIstring parerrname = (*vfecite)->getName() + "_err";
0217     ALIstring parerrtype = parerrname + "/D";
0218     FitParametersTree->Branch(parerrname.c_str(), &ParError[ii], parerrtype.c_str());
0219         ii++;
0220       }
0221     }
0222     ALIint ii = 0;
0223     for( vfecite = ((*vfescite)->FittedEntries()).begin(); vfecite != ((*vfescite)->FittedEntries()).end(); vfecite++) {
0224       ParValue[ii] = (*vfecite)->getValue();
0225       ParError[ii] = (*vfecite)->getSigma();
0226       ii++;
0227     }
0228     NFitParameters = ii;
0229     FitParametersTree->Fill();
0230     jj++;
0231   }
0232   */
0233 }
0234 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0235 //@@ Fill ntuple with optical object positions and orientations
0236 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0237 void NtupleManager::FillOptObjects(MatrixMeschach* AtWAMatrix) {
0238   int ii = 0;
0239   std::vector<OpticalObject*>::const_iterator vecobj;
0240   for (vecobj = Model::OptOList().begin(); vecobj != Model::OptOList().end(); ++vecobj) {
0241     OptObjectA = new ((*CloneOptObject)[ii]) OptObject();
0242
0243     OptObjectA->Name = (*vecobj)->name();
0244     OptObjectA->Type = (*vecobj)->type();
0245
0246     if (!(*vecobj)->parent()) {
0247       OptObjectA->Parent = ii;
0248       ii++;
0249       continue;
0250     }
0251
0252     int pp = 0;
0253     std::vector<OpticalObject*>::const_iterator vecobj2;
0254     for (vecobj2 = Model::OptOList().begin(); vecobj2 != Model::OptOList().end(); ++vecobj2) {
0255       if ((*vecobj2)->name() == (*vecobj)->parent()->name()) {
0256         OptObjectA->Parent = pp;
0257         continue;
0258       }
0259       pp++;
0260     }
0261
0262     OptObjectA->CentreGlobal[0] = 1000. * (*vecobj)->centreGlobal().x();
0263     OptObjectA->CentreGlobal[1] = 1000. * (*vecobj)->centreGlobal().y();
0264     OptObjectA->CentreGlobal[2] = 1000. * (*vecobj)->centreGlobal().z();
0265
0266     OptObjectA->CentreLocal[0] = 1000. * (*vecobj)->centreLocal().x();
0267     OptObjectA->CentreLocal[1] = 1000. * (*vecobj)->centreLocal().y();
0268     OptObjectA->CentreLocal[2] = 1000. * (*vecobj)->centreLocal().z();
0269
0270     OptObjectA->AnglesLocal[0] = (*vecobj)->getEntryRMangle(XCoor);
0271     OptObjectA->AnglesLocal[1] = (*vecobj)->getEntryRMangle(YCoor);
0272     OptObjectA->AnglesLocal[2] = (*vecobj)->getEntryRMangle(ZCoor);
0273
0274     double theta[3];
0275     GetGlobalAngles((*vecobj)->rmGlob(), theta);
0276     for (int i = 0; i < 3; i++)
0277       OptObjectA->AnglesGlobal[i] = theta[i];
0278
0279     ii++;
0280   }
0281
0282   NOptObjects = ii;
0283 }
0284 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0285 //@@ Fill ntuple with measurements
0286 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0287 void NtupleManager::FillMeasurements() {
0288   //---------- Loop Measurements
0289   int ss = 0, dd = 0, d1 = 0, tt = 0, cc = 0;
0290   std::vector<Measurement*>::const_iterator vmcite;
0291   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
0292     std::vector<ALIstring> optonamelist = (*vmcite)->OptONameList();
0293     int last = optonamelist.size() - 1;
0294     ALIstring LastOptOName = optonamelist[last];
0295     int optoind = -999;
0296     for (int no = 0; no < NOptObjects; no++) {
0297       OptObject* optobj = (OptObject*)CloneOptObject->At(no);
0298       if (optobj->Name == LastOptOName)
0299         optoind = no;
0300     }
0301     //std::cout << "DimSens " << (*vmcite)->type() << " " << (*vmcite)->sigma(0) << " " << LastOptOName << " " << optoind << std::endl;
0302     if ((*vmcite)->type() == "SENSOR2D") {
0303       Sensor2DMeasA = new ((*CloneSensor2DMeas)[ss]) Sensor2DMeas();
0304       Sensor2DMeasA->Name = (*vmcite)->name();
0305       Sensor2DMeasA->OptObjectIndex = optoind;
0306       for (ALIuint i = 0; i < (*vmcite)->dim(); i++) {
0307         Sensor2DMeasA->Position[i] = 1000. * (*vmcite)->value()[i];
0308         Sensor2DMeasA->PosError[i] = 1000. * (*vmcite)->sigma()[i];
0309         Sensor2DMeasA->SimulatedPosition[i] = 1000. * (*vmcite)->valueSimulated(i);
0310       }
0311       ss++;
0312     }
0313     if ((*vmcite)->type() == "DISTANCEMETER") {
0314       DistancemeterMeasA = new ((*CloneDistancemeterMeas)[dd]) DistancemeterMeas();
0315       DistancemeterMeasA->Name = (*vmcite)->name();
0316       DistancemeterMeasA->OptObjectIndex = optoind;
0317       DistancemeterMeasA->Distance = 1000. * (*vmcite)->value()[0];
0318       DistancemeterMeasA->DisError = 1000. * (*vmcite)->sigma()[0];
0319       DistancemeterMeasA->SimulatedDistance = 1000. * (*vmcite)->valueSimulated(0);
0320       dd++;
0321     }
0322     if ((*vmcite)->type() == "DISTANCEMETER1DIM") {
0323       Distancemeter1DimMeasA = new ((*CloneDistancemeter1DimMeas)[d1]) Distancemeter1DimMeas();
0324       Distancemeter1DimMeasA->Name = (*vmcite)->name();
0325       Distancemeter1DimMeasA->OptObjectIndex = optoind;
0326       Distancemeter1DimMeasA->Distance = 1000. * (*vmcite)->value()[0];
0327       Distancemeter1DimMeasA->DisError = 1000. * (*vmcite)->sigma()[0];
0328       Distancemeter1DimMeasA->SimulatedDistance = 1000. * (*vmcite)->valueSimulated(0);
0329       d1++;
0330     }
0331     if ((*vmcite)->type() == "TILTMETER") {
0332       TiltmeterMeasA = new ((*CloneTiltmeterMeas)[tt]) TiltmeterMeas();
0333       TiltmeterMeasA->Name = (*vmcite)->name();
0334       TiltmeterMeasA->OptObjectIndex = optoind;
0335       TiltmeterMeasA->Angle = (*vmcite)->value()[0];
0336       TiltmeterMeasA->AngError = (*vmcite)->sigma()[0];
0337       TiltmeterMeasA->SimulatedAngle = (*vmcite)->valueSimulated(0);
0338       tt++;
0339     }
0340     if ((*vmcite)->type() == "COPS") {
0341       CopsMeasA = new ((*CloneCopsMeas)[cc]) CopsMeas();
0342       CopsMeasA->Name = (*vmcite)->name();
0343       CopsMeasA->OptObjectIndex = optoind;
0344       for (ALIuint i = 0; i < (*vmcite)->dim(); i++) {
0345         CopsMeasA->Position[i] = 1000. * (*vmcite)->value()[i];
0346         CopsMeasA->PosError[i] = 1000. * (*vmcite)->sigma()[i];
0347         CopsMeasA->SimulatedPosition[i] = 1000. * (*vmcite)->valueSimulated(i);
0348       }
0349       cc++;
0350     }
0351   }
0352   NSensor2D = ss;
0353   NDistancemeter = dd;
0354   NDistancemeter1Dim = d1;
0355   NTiltmeter = tt;
0356   NCops = cc;
0357   //   MeasurementsTree->Fill();
0358 }
0359 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0360 //@@ Get global angles from global matrix rotation
0361 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0362 void NtupleManager::GetGlobalAngles(const CLHEP::HepRotation& rmGlob, double* theta) {
0363   double xx = rmGlob.xx();
0364   if (fabs(xx) < 1.e-08)
0365     xx = 0.;
0366   double xy = rmGlob.xy();
0367   if (fabs(xy) < 1.e-08)
0368     xy = 0.;
0369   double xz = rmGlob.xz();
0370   if (fabs(xz) < 1.e-08)
0371     xz = 0.;
0372   double yx = rmGlob.yx();
0373   if (fabs(yx) < 1.e-08)
0374     yx = 0.;
0375   double yy = rmGlob.yy();
0376   if (fabs(yy) < 1.e-08)
0377     yy = 0.;
0378   double yz = rmGlob.yz();
0379   if (fabs(yz) < 1.e-08)
0380     yz = 0.;
0381   double zx = rmGlob.zx();
0382   if (fabs(zx) < 1.e-08)
0383     zx = 0.;
0384   double zy = rmGlob.zy();
0385   if (fabs(zy) < 1.e-08)
0386     zy = 0.;
0387   double zz = rmGlob.zz();
0388   if (fabs(zz) < 1.e-08)
0389     zz = 0.;
0390
0391   double beta = asin(-zx);
0392
0393   double alpha, gamma;
0394   if (fabs(zx) != 1.) {
0395     double sinalpha = zy / cos(beta);
0396     double cosalpha = zz / cos(beta);
0397     if (cosalpha >= 0)
0398       alpha = asin(sinalpha);
0399     else
0400       alpha = M_PI - asin(sinalpha);
0401     if (alpha > M_PI)
0402       alpha -= 2 * M_PI;
0403
0404     double singamma = yx / cos(beta);
0405     double cosgamma = xx / cos(beta);
0406     if (cosgamma >= 0)
0407       gamma = asin(singamma);
0408     else
0409       gamma = M_PI - asin(singamma);
0410     if (gamma > M_PI)
0411       gamma -= 2 * M_PI;
0412
0413   } else {
0414     alpha = 0.;
0415
0416     double singamma = yz / sin(beta);
0417     double cosgamma = yy;
0418     if (cosgamma >= 0)
0419       gamma = asin(singamma);
0420     else
0421       gamma = M_PI - asin(singamma);
0422     if (gamma > M_PI)
0423       gamma -= 2 * M_PI;
0424   }
0425
0426   int GotGlobalAngles = 0;
0427   if (fabs(xy - (sin(alpha) * sin(beta) * cos(gamma) - sin(gamma) * cos(alpha))) > 1.e-08)
0428     GotGlobalAngles += 1;
0429   if (fabs(xz - (cos(alpha) * sin(beta) * cos(gamma) + sin(gamma) * sin(alpha))) > 1.e-08)
0430     GotGlobalAngles += 10;
0431   if (fabs(yy - (sin(alpha) * sin(beta) * sin(gamma) + cos(gamma) * cos(alpha))) > 1.e-08)
0432     GotGlobalAngles += 100;
0433   if (fabs(yz - (cos(alpha) * sin(beta) * sin(gamma) - cos(gamma) * sin(alpha))) > 1.e-08)
0434     GotGlobalAngles += 1000;
0435   if (GotGlobalAngles > 0)
0436     std::cout << "NtupleManager Warning: cannot get global rotation: " << GotGlobalAngles << std::endl;
0437
0438   theta[0] = alpha;
0439   theta[1] = beta;
0440   theta[2] = gamma;
0441 }