Reference Manual

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