CMSSW Reference Manual

00001 //   COCOA class implementation file
00002 //Id:  NtupleManager.cc
00003 //CAT: Analysis
00004 //
00005 //   History: v1.0 
00006 //   Luca Scodellaro
00007 #include <map>
00008 #include <fstream>
00009  
00010 #include "Alignment/CocoaAnalysis/interface/NtupleManager.h"
00011 #include "Alignment/CocoaAnalysis/interface/FittedEntry.h"
00012 // #include <iostream>
00013 // #include <iomanip>
00014 #include "Alignment/CocoaModel/interface/Measurement.h"
00015 #include "Alignment/CocoaModel/interface/OpticalObject.h"
00016 #include "Alignment/CocoaModel/interface/Entry.h"
00017 #include "Alignment/CocoaModel/interface/Model.h"
00018 // #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
00019 // #include "Alignment/CocoaUtilities/interface/GlobalOptionMgr.h"
00020 #include "TFile.h" 
00021 #include "TTree.h"
00022 #include "TClonesArray.h"
00023 #include "TMath.h"
00024 
00025 NtupleManager* NtupleManager::instance = 0;
00026 
00027 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00028 //@@  Gets the only instance of Model
00029 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00030 NtupleManager* NtupleManager::getInstance()
00031 {
00032   if(!instance) {
00033     instance = new NtupleManager;
00034   }
00035   return instance;
00036 }
00037 
00038 
00039 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00040 //@@ Book ntuple
00041 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00042 void NtupleManager::BookNtuple()
00043 {
00044   theRootFile = new TFile("report.root","RECREATE","Simple ROOT Ntuple");
00045 
00046   CocoaTree = new TTree("CocoaTree","CocoaTree"); 
00047 
00048   CocoaTree->Branch("Chi2Measurements",&Chi2Measurements,"Chi2Measurements/D");
00049   CocoaTree->Branch("Chi2CalibratedParameters",&Chi2CalibratedParameters,"Chi2CalibratedParameters/D");
00050   CocoaTree->Branch("NDegreesOfFreedom",&NDegreesOfFreedom,"NDegreesOfFreedom/I");
00051 
00052   CloneFitParam = new TClonesArray("FitParam");
00053   CocoaTree->Branch("FitParameters",&CloneFitParam,32000,2);
00054   CocoaTree->Branch("NFitParameters",&NFitParameters,"NFitParameters/I");
00055 
00056   CloneOptObject = new TClonesArray("OptObject");
00057   CocoaTree->Branch("OptObjects",&CloneOptObject,32000,2);
00058   CocoaTree->Branch("NOptObjects",&NOptObjects,"NOptObjects/I");
00059 
00060   CloneSensor2DMeas = new TClonesArray("Sensor2DMeas");
00061   CocoaTree->Branch("Sensor2DMeasurements",&CloneSensor2DMeas,32000,2);
00062   CocoaTree->Branch("NSensor2D",&NSensor2D,"NSensor2D/I");
00063   
00064   CloneDistancemeterMeas = new TClonesArray("DistancemeterMeas");
00065   CocoaTree->Branch("DistancemeterMeasurements",&CloneDistancemeterMeas,32000,2);
00066   CocoaTree->Branch("NDistancemeter",&NDistancemeter,"NDistancemeter/I");
00067   
00068   CloneDistancemeter1DimMeas = new TClonesArray("Distancemeter1DimMeas");
00069   CocoaTree->Branch("Distancemeter1DimMeasurements",&CloneDistancemeter1DimMeas,32000,2);
00070   CocoaTree->Branch("NDistancemeter1Dim",&NDistancemeter1Dim,"NDistancemeter1Dim/I");
00071 
00072   CloneTiltmeterMeas = new TClonesArray("TiltmeterMeas");
00073   CocoaTree->Branch("TiltmeterMeasurements",&CloneTiltmeterMeas,32000,2);
00074   CocoaTree->Branch("NTiltmeter",&NTiltmeter,"NTiltmeter/I");
00075 
00076   CloneCopsMeas = new TClonesArray("CopsMeas");
00077   CocoaTree->Branch("CopsMeasurements",&CloneCopsMeas,32000,2);
00078   CocoaTree->Branch("NCops",&NCops,"NCops/I");
00079 
00080   theRootFile->Add(CocoaTree);
00081 
00082 //   FitParametersTree = new TTree("FitParametersTree","FitParametersTree");
00083 //   FitParametersTree->Branch("NFitParameters",&NFitParameters,"NFitParameters/I");
00084 //   BookFitParameters = false;
00085 //   theRootFile->Add(FitParametersTree);
00086 
00087 //   MeasurementsTree = new TTree("MeasurementsTree","MeasurementsTree");
00088 //   MeasurementsTree->Branch("NMeasurements",&NMeasurements,"NMeasurements/I");
00089 //   BookMeasurements = false;
00090 //   theRootFile->Add(MeasurementsTree);
00091  
00092 
00093 }
00094 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00095 //@@ Init ntuple
00096 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00097 void NtupleManager::InitNtuple()
00098 {
00099   CloneFitParam->Clear();
00100 
00101   Chi2Measurements = 0.;
00102   Chi2CalibratedParameters = 0.;
00103   NDegreesOfFreedom = 0;
00104   NFitParameters = 0; 
00105   NOptObjects = 0; 
00106   NSensor2D = 0; 
00107   NDistancemeter = 0; 
00108   NDistancemeter1Dim = 0; 
00109   NTiltmeter = 0;
00110   NCops = 0;
00111 }
00112 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00113 //@@ Fill ntuple tree
00114 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00115 void NtupleManager::FillNtupleTree()
00116 {
00117   CocoaTree->Fill();
00118 }
00119 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00120 //@@ Close ntuple
00121 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00122 void NtupleManager::WriteNtuple()
00123 {
00124   theRootFile->Write();
00125   theRootFile->Close();
00126 }
00127 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00128 //@@ Close ntuple
00129 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00130 void NtupleManager::FillChi2()
00131 {
00132   double chi2meas = 0; 
00133   double chi2cal = 0;
00134   ALIint nMeas = 0, nUnk = 0;
00135 
00136   //----- Calculate the chi2 of measurements
00137   std::vector< Measurement* >::const_iterator vmcite;
00138   for ( vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); vmcite++) {
00139     for ( ALIuint ii = 0; ii < ALIuint((*vmcite)->dim()); ii++ ){
00140       nMeas++;
00141       double c2 = ( (*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii) ) / (*vmcite)->sigma(ii);
00142       chi2meas += c2*c2;
00143     }
00144   }
00145 
00146   //----- Calculate the chi2 of calibrated parameters
00147   std::vector< Entry* >::iterator veite;
00148   for ( veite = Model::EntryList().begin();
00149         veite != Model::EntryList().end(); veite++ ) {
00150     if ( (*veite)->quality() == 2 ) nUnk++;
00151     if ( (*veite)->quality() == 1 ) {
00152 //       std::cout << " " << (*veite)->valueDisplacementByFitting() << " " 
00153 //              << (*veite)->value << " " << (*veite)->sigma() << std::endl;
00154       double c2 = (*veite)->valueDisplacementByFitting() / (*veite)->sigma();
00155       chi2cal += c2*c2;
00156     }
00157   }
00158 
00159   Chi2Measurements = chi2meas;
00160   Chi2CalibratedParameters = chi2cal;
00161   NDegreesOfFreedom = nMeas - nUnk;
00162 
00163 }
00164 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00165 //@@ Fill ntuple with fitted parameters
00166 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00167 void NtupleManager::FillFitParameters(MatrixMeschach* AtWAMatrix)
00168 {
00169 
00170 //   double ParValue[1000], ParError[1000];
00171   int theMinEntryQuality = 1; 
00172   int ii = 0;
00173   std::vector<Entry*>::const_iterator vecite; 
00174   for ( vecite = Model::EntryList().begin();
00175     vecite != Model::EntryList().end(); vecite++ ) {
00176 
00177     //--- Only for good quality parameters (='unk')
00178     if ( (*vecite)->quality() >= theMinEntryQuality ) {
00179 
00180       ALIint ipos = (*vecite)->fitPos();
00181       FittedEntry* fe = new FittedEntry( (*vecite), ipos, sqrt(AtWAMatrix->Mat()->me[ipos][ipos]));
00182 //       if (!BookFitParameters) {
00183 //      CocoaTree->Branch("NFitParameters",&NFitParameters,"NFitParameters/I:");
00184 //      ALIstring partype = fe->getName() + "/D";
00185 //      FitParametersTree->Branch(fe->getName().c_str(), &ParValue[ii], partype.c_str());
00186 //      ALIstring parerrname = fe->getName() + "_err";
00187 //      ALIstring parerrtype = parerrname + "/D";
00188 //      FitParametersTree->Branch(parerrname.c_str(), &ParError[ii], parerrtype.c_str());
00189 //       }
00190 //       ParValue[ii] = fe->getValue();
00191 //       ParError[ii] = fe->getSigma();
00192      std::cout << "EEE " << (*vecite)->ValueDimensionFactor() << " " << (*vecite)->SigmaDimensionFactor() << " " << fe->getOptOName() << " " << fe->getEntryName() << " " << fe->getName() << " " << fe->getOrder() << " " << fe->getQuality() << " " << (*vecite)->type() << " " << std::endl;
00193       FitParamA = new( (*CloneFitParam)[ii] ) FitParam();
00194       FitParamA->Name = fe->getName();
00195       if (fe->getQuality()==1) FitParamA->Quality = "Calibrated";
00196       else if (fe->getQuality()==2) FitParamA->Quality = "Unknown";
00197       for (int no = 0; no<NOptObjects; no++) {
00198         OptObject* optobj = (OptObject*) CloneOptObject->At(no);
00199         if (optobj->Name==fe->getOptOName()) FitParamA->OptObjectIndex = no;
00200       }
00201       float DF = 1.;
00202       if ((*vecite)->type()=="centre" || (*vecite)->type()=="length") DF = 1000.;
00203       FitParamA->InitialValue = DF*fe->getOrigValue()*(*vecite)->ValueDimensionFactor();
00204       FitParamA->InitialSigma = DF*fe->getOrigSigma()*(*vecite)->SigmaDimensionFactor();
00205       FitParamA->FittedValue  = DF*fe->getValue()*(*vecite)->ValueDimensionFactor();
00206       FitParamA->FittedSigma  = DF*fe->getSigma()*(*vecite)->SigmaDimensionFactor();
00207       ii++;
00208 
00209     }
00210 
00211   }
00212 //   BookFitParameters = true;
00213   NFitParameters = ii;
00214 //   FitParametersTree->Fill();
00215 
00216   /*  
00217   //---------- Loop sets of entries
00218   std::vector< FittedEntriesSet* > theFittedEntriesSets;
00219   std::vector< FittedEntriesSet* >::const_iterator vfescite;
00220   std::vector< FittedEntry* >::const_iterator vfecite;
00221   ALIint jj = 1;
00222   for( vfescite = theFittedEntriesSets.begin(); vfescite != theFittedEntriesSets.end(); vfescite++) {
00223     //---------- Loop entries
00224     if( vfescite == theFittedEntriesSets.begin() ) {
00225     //----- dump entries names if first set 
00226       ALIint ii = 0;
00227       for( vfecite = ((*vfescite)->FittedEntries()).begin(); vfecite != ((*vfescite)->FittedEntries()).end(); vfecite++) {
00228         ALIstring partype = (*vfecite)->getName() + "/D";
00229         FitParametersTree->Branch((*vfecite)->getName().c_str(), &ParValue[ii], partype.c_str());
00230         ALIstring parerrname = (*vfecite)->getName() + "_err";
00231         ALIstring parerrtype = parerrname + "/D";
00232         FitParametersTree->Branch(parerrname.c_str(), &ParError[ii], parerrtype.c_str());
00233         ii++;
00234       }
00235     }
00236     ALIint ii = 0;
00237     for( vfecite = ((*vfescite)->FittedEntries()).begin(); vfecite != ((*vfescite)->FittedEntries()).end(); vfecite++) {
00238       ParValue[ii] = (*vfecite)->getValue();
00239       ParError[ii] = (*vfecite)->getSigma();
00240       ii++;
00241     }
00242     NFitParameters = ii;
00243     FitParametersTree->Fill();
00244     jj++;
00245   }
00246   */
00247 
00248 }
00249 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00250 //@@ Fill ntuple with optical object positions and orientations
00251 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00252 void NtupleManager::FillOptObjects(MatrixMeschach* AtWAMatrix)
00253 {
00254 
00255   int ii = 0;
00256   std::vector< OpticalObject* >::const_iterator vecobj;
00257   for ( vecobj = Model::OptOList().begin();
00258         vecobj != Model::OptOList().end(); vecobj++ ) {
00259     OptObjectA = new( (*CloneOptObject)[ii] ) OptObject();
00260 
00261     OptObjectA->Name =  (*vecobj)->name();
00262     OptObjectA->Type = (*vecobj)->type();
00263 
00264     if (!(*vecobj)->parent()) {
00265       OptObjectA->Parent = ii;
00266       ii++;
00267       continue;
00268     }
00269 
00270     int pp = 0;
00271     std::vector< OpticalObject* >::const_iterator vecobj2;
00272     for ( vecobj2 = Model::OptOList().begin();
00273           vecobj2 != Model::OptOList().end(); vecobj2++ ) {
00274       if ((*vecobj2)->name()==(*vecobj)->parent()->name()) {
00275         OptObjectA->Parent = pp;
00276         continue;
00277       }
00278       pp++;
00279     }
00280     
00281     OptObjectA->CentreGlobal[0] = 1000.*(*vecobj)->centreGlobal().x();
00282     OptObjectA->CentreGlobal[1] = 1000.*(*vecobj)->centreGlobal().y();
00283     OptObjectA->CentreGlobal[2] = 1000.*(*vecobj)->centreGlobal().z();
00284 
00285     OptObjectA->CentreLocal[0] = 1000.*(*vecobj)->centreLocal().x();
00286     OptObjectA->CentreLocal[1] = 1000.*(*vecobj)->centreLocal().y();
00287     OptObjectA->CentreLocal[2] = 1000.*(*vecobj)->centreLocal().z();
00288 
00289     OptObjectA->AnglesLocal[0] = (*vecobj)->getEntryRMangle(XCoor);
00290     OptObjectA->AnglesLocal[1] = (*vecobj)->getEntryRMangle(YCoor);
00291     OptObjectA->AnglesLocal[2] = (*vecobj)->getEntryRMangle(ZCoor);
00292 
00293     double theta[3];
00294     GetGlobalAngles((*vecobj)->rmGlob(), theta);
00295     for (int i = 0; i<3; i++) OptObjectA->AnglesGlobal[i] = theta[i];
00296 
00297     ii++;
00298 
00299   }
00300 
00301   NOptObjects = ii;
00302 
00303 
00304 }
00305 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00306 //@@ Fill ntuple with measurements
00307 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00308 void NtupleManager::FillMeasurements()
00309 {
00310   //---------- Loop Measurements
00311   int ss = 0, dd = 0, d1 = 0, tt = 0, cc = 0;
00312   std::vector< Measurement* >::const_iterator vmcite;
00313   for ( vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); vmcite++) {
00314     std::vector<ALIstring> optonamelist = (*vmcite)->OptONameList(); 
00315     int last = optonamelist.size() - 1; ALIstring LastOptOName = optonamelist[last];
00316     int optoind = -999;
00317     for (int no = 0; no<NOptObjects; no++) {
00318       OptObject* optobj = (OptObject*) CloneOptObject->At(no);
00319       if (optobj->Name==LastOptOName) optoind = no;
00320     }
00321     //std::cout << "DimSens " << (*vmcite)->type() << " " << (*vmcite)->sigma(0) << " " << LastOptOName << " " << optoind << std::endl;
00322     if ((*vmcite)->type()=="SENSOR2D") {
00323       Sensor2DMeasA = new( (*CloneSensor2DMeas)[ss] ) Sensor2DMeas();
00324       Sensor2DMeasA->Name = (*vmcite)->name();
00325       Sensor2DMeasA->OptObjectIndex = optoind;
00326       for (uint i = 0; i<(*vmcite)->dim(); i++) {
00327         Sensor2DMeasA->Position[i] = 1000.*(*vmcite)->value()[i];
00328         Sensor2DMeasA->PosError[i] = 1000.*(*vmcite)->sigma()[i];
00329         Sensor2DMeasA->SimulatedPosition[i] = 1000.*(*vmcite)->valueSimulated(i);
00330       }
00331       ss++;
00332     }
00333     if ((*vmcite)->type()=="DISTANCEMETER") {
00334       DistancemeterMeasA = new( (*CloneDistancemeterMeas)[dd] ) DistancemeterMeas();
00335       DistancemeterMeasA->Name = (*vmcite)->name();
00336       DistancemeterMeasA->OptObjectIndex = optoind;
00337       DistancemeterMeasA->Distance = 1000.*(*vmcite)->value()[0];
00338       DistancemeterMeasA->DisError = 1000.*(*vmcite)->sigma()[0];
00339       DistancemeterMeasA->SimulatedDistance = 1000.*(*vmcite)->valueSimulated(0);
00340       dd++;
00341     }
00342     if ((*vmcite)->type()=="DISTANCEMETER1DIM") {
00343       Distancemeter1DimMeasA = new( (*CloneDistancemeter1DimMeas)[d1] ) Distancemeter1DimMeas();
00344       Distancemeter1DimMeasA->Name = (*vmcite)->name();
00345       Distancemeter1DimMeasA->OptObjectIndex = optoind;
00346       Distancemeter1DimMeasA->Distance = 1000.*(*vmcite)->value()[0];
00347       Distancemeter1DimMeasA->DisError = 1000.*(*vmcite)->sigma()[0];
00348       Distancemeter1DimMeasA->SimulatedDistance = 1000.*(*vmcite)->valueSimulated(0);
00349       d1++;
00350 
00351     }
00352     if ((*vmcite)->type()=="TILTMETER") {
00353       TiltmeterMeasA = new( (*CloneTiltmeterMeas)[tt] ) TiltmeterMeas();
00354       TiltmeterMeasA->Name = (*vmcite)->name();
00355       TiltmeterMeasA->OptObjectIndex = optoind;
00356       TiltmeterMeasA->Angle    = (*vmcite)->value()[0];
00357       TiltmeterMeasA->AngError = (*vmcite)->sigma()[0];
00358       TiltmeterMeasA->SimulatedAngle = (*vmcite)->valueSimulated(0);
00359       tt++;
00360     }
00361     if ((*vmcite)->type()=="COPS") {
00362       CopsMeasA = new( (*CloneCopsMeas)[cc] ) CopsMeas();
00363       CopsMeasA->Name = (*vmcite)->name();
00364       CopsMeasA->OptObjectIndex = optoind;
00365       for (uint i = 0; i<(*vmcite)->dim(); i++) {
00366         CopsMeasA->Position[i] = 1000.*(*vmcite)->value()[i];
00367         CopsMeasA->PosError[i] = 1000.*(*vmcite)->sigma()[i];
00368         CopsMeasA->SimulatedPosition[i] = 1000.*(*vmcite)->valueSimulated(i);
00369       }
00370       cc++;
00371     }
00372   }
00373   NSensor2D = ss; 
00374   NDistancemeter = dd; 
00375   NDistancemeter1Dim = d1; 
00376   NTiltmeter = tt; 
00377   NCops = cc; 
00378   //   MeasurementsTree->Fill();
00379 }
00380 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00381 //@@ Get global angles from global matrix rotation
00382 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00383 void NtupleManager::GetGlobalAngles(const HepRotation& rmGlob, double *theta)
00384 {
00385 
00386   double xx = rmGlob.xx(); if (fabs(xx)<1.e-08) xx = 0.;
00387   double xy = rmGlob.xy(); if (fabs(xy)<1.e-08) xy = 0.;
00388   double xz = rmGlob.xz(); if (fabs(xz)<1.e-08) xz = 0.;
00389   double yx = rmGlob.yx(); if (fabs(yx)<1.e-08) yx = 0.;
00390   double yy = rmGlob.yy(); if (fabs(yy)<1.e-08) yy = 0.;
00391   double yz = rmGlob.yz(); if (fabs(yz)<1.e-08) yz = 0.;
00392   double zx = rmGlob.zx(); if (fabs(zx)<1.e-08) zx = 0.;
00393   double zy = rmGlob.zy(); if (fabs(zy)<1.e-08) zy = 0.;
00394   double zz = rmGlob.zz(); if (fabs(zz)<1.e-08) zz = 0.;
00395 
00396   double beta = asin(-zx);
00397 
00398   double alpha, gamma;
00399   if (fabs(zx)!=1.) {
00400   
00401     double sinalpha = zy/cos(beta);
00402     double cosalpha = zz/cos(beta);
00403     if (cosalpha>=0) alpha = asin(sinalpha);
00404     else alpha = TMath::Pi() - asin(sinalpha);
00405     if (alpha>TMath::Pi()) alpha -= 2*TMath::Pi();
00406   
00407     double singamma = yx/cos(beta);
00408     double cosgamma = xx/cos(beta);
00409     if (cosgamma>=0) gamma = asin(singamma);
00410     else gamma = TMath::Pi() - asin(singamma);
00411     if (gamma>TMath::Pi()) gamma -= 2*TMath::Pi();
00412     
00413   } else {
00414 
00415     alpha = 0.;
00416     
00417     double singamma = yz/sin(beta);
00418     double cosgamma = yy;
00419     if (cosgamma>=0) gamma = asin(singamma);
00420     else gamma = TMath::Pi() - asin(singamma);
00421     if (gamma>TMath::Pi()) gamma -= 2*TMath::Pi();
00422 
00423   }
00424 
00425   int GotGlobalAngles = 0;
00426   if (fabs(xy-(sin(alpha)*sin(beta)*cos(gamma)-sin(gamma)*cos(alpha)))>1.e-08)
00427     GotGlobalAngles += 1;
00428   if (fabs(xz-(cos(alpha)*sin(beta)*cos(gamma)+sin(gamma)*sin(alpha)))>1.e-08) 
00429     GotGlobalAngles += 10;
00430   if (fabs(yy-(sin(alpha)*sin(beta)*sin(gamma)+cos(gamma)*cos(alpha)))>1.e-08)
00431     GotGlobalAngles += 100;
00432   if (fabs(yz-(cos(alpha)*sin(beta)*sin(gamma)-cos(gamma)*sin(alpha)))>1.e-08) 
00433     GotGlobalAngles += 1000;
00434   if (GotGlobalAngles>0) 
00435     std::cout << "NtupleManager Warning: cannot get global rotation: " 
00436               << GotGlobalAngles << std::endl;
00437 
00438   theta[0] = alpha;
00439   theta[1] = beta;
00440   theta[2] = gamma;
00441 
00442 }
00443