00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "DQM/Physics/src/EwkDQM.h"
00010 #include "DataFormats/Candidate/interface/Candidate.h"
00011 #include "FWCore/Framework/interface/Event.h"
00012 #include "FWCore/Framework/interface/EventSetup.h"
00013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00014
00015 #include "DQMServices/Core/interface/DQMStore.h"
00016 #include "DQMServices/Core/interface/MonitorElement.h"
00017 #include "FWCore/ServiceRegistry/interface/Service.h"
00018
00019 #include "DataFormats/Common/interface/Handle.h"
00020
00021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00022 #include "FWCore/Common/interface/TriggerNames.h"
00023
00024
00025 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00026 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00027 #include "DataFormats/MuonReco/interface/Muon.h"
00028 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00029 #include "DataFormats/JetReco/interface/Jet.h"
00030 #include "DataFormats/METReco/interface/MET.h"
00031 #include "DataFormats/TrackReco/interface/Track.h"
00032 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00033 #include "DataFormats/VertexReco/interface/Vertex.h"
00034 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00035
00036
00037 #include "DataFormats/Math/interface/LorentzVector.h"
00038 #include "TLorentzVector.h"
00039
00040 #include <vector>
00041
00042 #include <string>
00043 #include <cmath>
00044 using namespace std;
00045 using namespace edm;
00046 using namespace reco;
00047
00048
00049
00050
00051 EwkDQM::EwkDQM(const ParameterSet& parameters) {
00052 eJetMin_ = parameters.getUntrackedParameter<double>("EJetMin", 999999.);
00053
00054
00055
00056 theElecTriggerPathToPass_ = parameters.getParameter<string>("elecTriggerPathToPass");
00057 theMuonTriggerPathToPass_ = parameters.getParameter<string>("muonTriggerPathToPass");
00058
00059
00060 theTriggerResultsCollection_ = parameters.getParameter<InputTag>("triggerResultsCollection");
00061 theMuonCollectionLabel_ = parameters.getParameter<InputTag>("muonCollection");
00062 theElectronCollectionLabel_ = parameters.getParameter<InputTag>("electronCollection");
00063
00064 thePFJetCollectionLabel_ = parameters.getParameter<InputTag>("PFJetCollection");
00065 theCaloMETCollectionLabel_ = parameters.getParameter<InputTag>("caloMETCollection");
00066
00067
00068 isValidHltConfig_ = false;
00069
00070
00071 h_vertex_number = 0;
00072 h_vertex_chi2 = 0;
00073 h_vertex_numTrks = 0;
00074 h_vertex_sumTrks = 0;
00075 h_vertex_d0 = 0;
00076
00077 h_jet_count = 0;
00078 h_jet_et = 0;
00079 h_jet_pt = 0;
00080 h_jet_eta = 0;
00081 h_jet_phi = 0;
00082 h_jet2_et = 0;
00083
00084 h_jet2_eta = 0;
00085 h_jet2_phi = 0;
00086
00087 h_e1_et = 0;
00088 h_e2_et = 0;
00089 h_e1_eta = 0;
00090 h_e2_eta = 0;
00091 h_e1_phi = 0;
00092 h_e2_phi = 0;
00093
00094 h_m1_pt = 0;
00095 h_m2_pt = 0;
00096 h_m1_eta = 0;
00097 h_m2_eta = 0;
00098 h_m1_phi = 0;
00099 h_m2_phi = 0;
00100
00101
00102
00103
00104
00105 h_met = 0;
00106 h_met_phi = 0;
00107
00108 h_e_invWMass = 0;
00109 h_m_invWMass = 0;
00110 h_mumu_invMass = 0;
00111 h_ee_invMass = 0;
00112
00113 theDbe = Service<DQMStore>().operator->();
00114
00115 }
00116
00117 EwkDQM::~EwkDQM() {
00118 }
00119
00120
00121 void EwkDQM::beginJob() {
00122
00123 char chtitle[256] = "";
00124
00125 logTraceName = "EwkAnalyzer";
00126
00127 LogTrace(logTraceName)<<"Parameters initialization";
00128 theDbe->setCurrentFolder("Physics/EwkDQM");
00129
00130 const float pi = 4*atan(1);
00131
00132
00133 h_vertex_number = theDbe->book1D("vertex_number", "Number of event vertices in collection", 10,-0.5, 9.5 );
00134 h_vertex_chi2 = theDbe->book1D("vertex_chi2" , "Event Vertex #chi^{2}/n.d.o.f." , 20, 0.0, 2.0 );
00135 h_vertex_numTrks = theDbe->book1D("vertex_numTrks", "Event Vertex, number of tracks" , 20, -0.5, 59.5 );
00136 h_vertex_sumTrks = theDbe->book1D("vertex_sumTrks", "Event Vertex, sum of track pt" , 20, 0.0, 100.0 );
00137 h_vertex_d0 = theDbe->book1D("vertex_d0" , "Event Vertex d0" , 20, 0.0, 0.05);
00138
00139 snprintf(chtitle, 255, "Number of %s (E_{T} > 15 GeV);Number of Jets",thePFJetCollectionLabel_.label().data());
00140 h_jet_count = theDbe->book1D("jet_count", chtitle, 8, -0.5, 7.5);
00141
00142 snprintf(chtitle, 255, "Leading jet E_{T} (from %s);E_{T}(1^{st} jet) (GeV)",
00143 thePFJetCollectionLabel_.label().data());
00144 h_jet_et = theDbe->book1D("jet_et", chtitle, 20, 0., 200.0);
00145
00146 snprintf(chtitle, 255, "Leading jet p_{T} (from %s);p_{T}(1^{st} jet) (GeV/c)", thePFJetCollectionLabel_.label().data());
00147 h_jet_pt = theDbe->book1D("jet_pt", chtitle, 20, 0., 200.0);
00148
00149 snprintf(chtitle, 255, "Leading jet #eta (from %s); #eta (1^{st} jet)",thePFJetCollectionLabel_.label().data());
00150 h_jet_eta = theDbe->book1D("jet_eta", chtitle, 20, -10., 10.0);
00151 snprintf(chtitle, 255, "Leading jet #phi (from %s); #phi(1^{st} jet)", thePFJetCollectionLabel_.label().data());
00152 h_jet_phi = theDbe->book1D("jet_phi", chtitle, 22, -1.1*pi, 1.1*pi);
00153
00154 snprintf(chtitle, 255, "2^{nd} leading jet E_{T} (from %s);E_{T}(2^{nd} jet) (GeV)",thePFJetCollectionLabel_.label().data());
00155 h_jet2_et = theDbe->book1D("jet2_et", chtitle, 20, 0., 200.0);
00156
00157
00158 snprintf(chtitle, 255, "2^{nd} leading jet #eta (from %s); #eta (2^{nd} jet)",thePFJetCollectionLabel_.label().data());
00159 h_jet2_eta = theDbe->book1D("jet2_eta", chtitle, 20, -10., 10.0);
00160 snprintf(chtitle, 255, "2^{nd} leading jet #phi (from %s); #phi(2^{nd} jet)", thePFJetCollectionLabel_.label().data());
00161 h_jet2_phi = theDbe->book1D("jet2_phi", chtitle, 22, -1.1*pi, 1.1*pi);
00162
00163 h_e1_et = theDbe->book1D("e1_et", "E_{T} of Leading Electron;E_{T} (GeV)" , 20, 0.0 , 100.0);
00164 h_e2_et = theDbe->book1D("e2_et", "E_{T} of Second Electron;E_{T} (GeV)" , 20, 0.0 , 100.0);
00165 h_e1_eta = theDbe->book1D("e1_eta", "#eta of Leading Electron;#eta" , 20, -4.0 , 4.0);
00166 h_e2_eta = theDbe->book1D("e2_eta", "#eta of Second Electron;#eta" , 20, -4.0 , 4.0);
00167 h_e1_phi = theDbe->book1D("e1_phi", "#phi of Leading Electron;#phi" , 22, -1.1*pi, 1.1*pi );
00168 h_e2_phi = theDbe->book1D("e2_phi", "#phi of Second Electron;#phi" , 22, -1.1*pi, 1.1*pi );
00169 h_m1_pt = theDbe->book1D("m1_pt", "p_{T} of Leading Muon;p_{T}(1^{st} #mu) (GeV)", 20, 0.0 , 100.0);
00170 h_m2_pt = theDbe->book1D("m2_pt", "p_{T} of Second Muon;p_{T}(2^{nd} #mu) (GeV)" , 20, 0.0 , 100.0);
00171 h_m1_eta = theDbe->book1D("m1_eta", "#eta of Leading Muon;#eta(1^{st} #mu)" , 20, -4.0 , 4.0);
00172 h_m2_eta = theDbe->book1D("m2_eta", "#eta of Second Muon;#eta(2^{nd} #mu)" , 20, -4.0 , 4.0);
00173 h_m1_phi = theDbe->book1D("m1_phi", "#phi of Leading Muon;#phi(1^{st} #mu)" , 20, (-1.-1./10.)*pi, (1.+1./10.)*pi);
00174 h_m2_phi = theDbe->book1D("m2_phi", "#phi of Second Muon;#phi(2^{nd} #mu)" , 20, (-1.-1./10.)*pi, (1.+1./10.)*pi);
00175
00176
00177
00178 snprintf(chtitle, 255, "Missing E_{T} (%s); GeV", theCaloMETCollectionLabel_.label().data());
00179 h_met = theDbe->book1D("met", chtitle, 20, 0.0 , 100);
00180 h_met_phi = theDbe->book1D("met_phi", "Missing E_{T} #phi;#phi(MET)" , 22, (-1.-1./10.)*pi, (1.+1./10.)*pi );
00181
00182 h_e_invWMass = theDbe->book1D("we_invWMass", "W-> e #nu Transverse Mass;M_{T} (GeV)" , 20, 0.0, 140.0);
00183 h_m_invWMass = theDbe->book1D("wm_invWMass", "W-> #mu #nu Transverse Mass;M_{T} (GeV)" , 20, 0.0, 140.0);
00184
00185 h_mumu_invMass = theDbe->book1D("z_mm_invMass", "#mu#mu Invariant Mass;InvMass (GeV)" , 20, 40.0, 140.0 );
00186 h_ee_invMass = theDbe->book1D("z_ee_invMass", "ee Invariant Mass;InvMass (Gev)" , 20, 40.0, 140.0 );
00187 }
00188
00189
00190
00194 void EwkDQM::beginRun( const edm::Run& theRun, const edm::EventSetup& theSetup ) {
00195
00196
00197 bool isConfigChanged = false;
00198
00199
00200 const std::string hltProcessName( theTriggerResultsCollection_.process() );
00201 isValidHltConfig_ = hltConfigProvider_.init( theRun, theSetup, hltProcessName, isConfigChanged );
00202
00203 }
00204
00205
00206 void EwkDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
00207
00208
00209 if( ! isValidHltConfig_ ) return;
00210
00211
00212 LogTrace(logTraceName)<<"Analysis of event # ";
00213
00214 Handle<TriggerResults> HLTresults;
00215 iEvent.getByLabel(theTriggerResultsCollection_, HLTresults);
00216 if ( !HLTresults.isValid() ) return;
00217
00218 const edm::TriggerNames & trigNames = iEvent.triggerNames(*HLTresults);
00219
00220
00221 std::vector<std::string> eleTrigPathNames;
00222 std::vector<std::string> muTrigPathNames;
00223 eleTrigPathNames.push_back(theElecTriggerPathToPass_);
00224 muTrigPathNames.push_back(theMuonTriggerPathToPass_);
00225
00226
00227 bool passed_electron_HLT = false;
00228 bool passed_muon_HLT = false;
00229 for (unsigned int i=0; i<HLTresults->size(); i++) {
00230 const std::string trigName = trigNames.triggerName(i);
00231
00232 for(unsigned int index=0; index<eleTrigPathNames.size() && !passed_electron_HLT; index++) {
00233 size_t trigPath = trigName.find(eleTrigPathNames[index]);
00234 if (trigPath==0) {
00235 passed_electron_HLT = HLTresults->accept(i);
00236 }
00237 }
00238
00239 for(unsigned int index=0; index<muTrigPathNames.size() && !passed_muon_HLT; index++) {
00240 size_t trigPath = trigName.find(muTrigPathNames[index]);
00241 if (trigPath==0) {
00242 passed_muon_HLT = HLTresults->accept(i);
00243 }
00244 }
00245 }
00246
00247
00248 if ( !(passed_electron_HLT || passed_muon_HLT) ) return;
00249
00251
00252 Handle<VertexCollection> vertexHandle;
00253 iEvent.getByLabel("offlinePrimaryVertices", vertexHandle);
00254 if ( !vertexHandle.isValid() ) return;
00255 VertexCollection vertexCollection = *(vertexHandle.product());
00256 int vertex_number = vertexCollection.size();
00257 VertexCollection::const_iterator v = vertexCollection.begin();
00258 double vertex_chi2 = v->normalizedChi2();
00259 double vertex_d0 = sqrt(v->x()*v->x()+v->y()*v->y());
00260
00261
00262 double vertex_numTrks = v->tracksSize();
00263 double vertex_sumTrks = 0.0;
00264 for (Vertex::trackRef_iterator vertex_curTrack = v->tracks_begin(); vertex_curTrack!=v->tracks_end(); vertex_curTrack++) {
00265 vertex_sumTrks += (*vertex_curTrack)->pt();
00266 }
00267
00269
00270 Handle< View<MET> > caloMETCollection;
00271 iEvent.getByLabel(theCaloMETCollectionLabel_, caloMETCollection);
00272 if ( !caloMETCollection.isValid() ) return;
00273 float missing_et = caloMETCollection->begin()->et();
00274 float met_phi = caloMETCollection->begin()->phi();
00275
00276
00278
00279 Handle<GsfElectronCollection> electronCollection;
00280 iEvent.getByLabel(theElectronCollectionLabel_, electronCollection);
00281 if ( !electronCollection.isValid() ) return;
00282
00283
00284 float electron_et = -8.0;
00285 float electron_eta = -8.0;
00286 float electron_phi = -8.0;
00287 float electron2_et = -9.0;
00288 float electron2_eta = -9.0;
00289 float electron2_phi = -9.0;
00290 float ee_invMass = -9.0;
00291 TLorentzVector e1, e2;
00292
00293
00294 if( passed_electron_HLT ) {
00295
00296 for (reco::GsfElectronCollection::const_iterator recoElectron=electronCollection->begin(); recoElectron!=electronCollection->end(); recoElectron++){
00297
00298
00299 if ( recoElectron->et() < 20 || fabs(recoElectron->eta())>2.5 ) continue;
00300
00301
00302 if ( recoElectron->deltaPhiSuperClusterTrackAtVtx() > 0.58 ||
00303 recoElectron->deltaEtaSuperClusterTrackAtVtx() > 0.01 ||
00304 recoElectron->sigmaIetaIeta() > 0.027 ) continue;
00305
00306 if (recoElectron->et() > electron_et){
00307 electron2_et = electron_et;
00308 electron2_eta = electron_eta;
00309 electron2_phi = electron_phi;
00310 electron_et = recoElectron->et();
00311 electron_eta = recoElectron->eta();
00312 electron_phi = recoElectron->phi();
00313 e1 = TLorentzVector(recoElectron->momentum().x(),recoElectron->momentum().y(),recoElectron->momentum().z(),recoElectron->p());
00314 } else if (recoElectron->et() > electron2_et) {
00315 electron2_et = recoElectron->et();
00316 electron2_eta = recoElectron->eta();
00317 electron2_phi = recoElectron->phi();
00318 e2 = TLorentzVector(recoElectron->momentum().x(),recoElectron->momentum().y(),recoElectron->momentum().z(),recoElectron->p());
00319 }
00320 }
00321 if (electron2_et>0.0) {
00322 TLorentzVector pair=e1+e2;
00323 ee_invMass = pair.M();
00324 }
00325 }
00327
00328
00329
00331
00332 Handle<MuonCollection> muonCollection;
00333 iEvent.getByLabel(theMuonCollectionLabel_,muonCollection);
00334 if ( !muonCollection.isValid() ) return;
00335
00336
00337 float mm_invMass = -9.0;
00338 float muon_pt = -9.0;
00339 float muon_eta = -9.0;
00340 float muon_phi = -9.0;
00341 float muon2_pt = -9.0;
00342 float muon2_eta = -9.0;
00343 float muon2_phi = -9.0;
00344 TLorentzVector m1, m2;
00345
00346 if( passed_muon_HLT ) {
00347 for (reco::MuonCollection::const_iterator recoMuon=muonCollection->begin(); recoMuon!=muonCollection->end(); recoMuon++){
00348
00349
00350 if ( recoMuon->pt() < 20 || !recoMuon->isGlobalMuon() ) continue;
00351
00352 if ( recoMuon->globalTrack()->normalizedChi2() > 10 ) continue;
00353
00354 if (recoMuon->pt() > muon_pt){
00355 muon2_pt = muon_pt;
00356 muon2_eta = muon_eta;
00357 muon2_phi = muon_phi;
00358 muon_pt = recoMuon->pt();
00359 muon_eta = recoMuon->eta();
00360 muon_phi = recoMuon->phi();
00361 m1 = TLorentzVector(recoMuon->momentum().x(),recoMuon->momentum().y(),recoMuon->momentum().z(),recoMuon->p());
00362 } else if (recoMuon->pt() > muon2_pt) {
00363 muon2_pt = recoMuon->pt();
00364 muon2_eta = recoMuon->eta();
00365 muon2_phi = recoMuon->phi();
00366 m2 = TLorentzVector(recoMuon->momentum().x(),recoMuon->momentum().y(),recoMuon->momentum().z(),recoMuon->p());
00367 }
00368 }
00369 }
00370 if (muon2_pt>0.0) {
00371 TLorentzVector pair=m1+m2;
00372 mm_invMass = pair.M();
00373 }
00375
00376
00378
00379
00380
00381 Handle<View<Jet> > PFJetCollection;
00382
00383 iEvent.getByLabel (thePFJetCollectionLabel_,PFJetCollection);
00384
00385 if ( !PFJetCollection.isValid() ) return;
00386
00387 unsigned int muonCollectionSize = muonCollection->size();
00388
00389 unsigned int PFJetCollectionSize = PFJetCollection->size();
00390 int jet_count = 0;
00391
00392
00393
00394 float jet_et = -80.0;
00395 float jet_pt = -80.0;
00396 float jet_eta = -80.0;
00397 float jet_phi = -80.0;
00398 float jet2_et = -90.0;
00399 float jet2_eta = -90.0;
00400 float jet2_phi = -90.0;
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418 for (unsigned int i=0; i<PFJetCollectionSize; i++) {
00419 const Jet& jet = PFJetCollection->at(i);
00420
00421 double minDistance=99999;
00422 for (unsigned int j=0; j<muonCollectionSize; j++) {
00423 const Muon& mu = muonCollection->at(j);
00424 double distance = sqrt( (mu.eta()-jet.eta())*(mu.eta()-jet.eta()) +(mu.phi()-jet.phi())*(mu.phi()-jet.phi()) );
00425 if (minDistance>distance) minDistance=distance;
00426 }
00427 if (minDistance<0.3) continue;
00428
00429
00430
00431 if ( electron_et>0.0 && fabs(jet.eta()-electron_eta ) < 0.2 && calcDeltaPhi(jet.phi(), electron_phi ) < 0.2) continue;
00432 if ( electron2_et>0.0&& fabs(jet.eta()-electron2_eta) < 0.2 && calcDeltaPhi(jet.phi(), electron2_phi) < 0.2) continue;
00433
00434
00435
00436
00437
00438
00439
00440
00441 if (jet.et() < eJetMin_) continue;
00442 jet_count ++;
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455 if (jet.et() > jet_et) {
00456 jet2_et = jet_et;
00457 jet2_eta = jet_eta;
00458 jet2_phi = jet_phi;
00459
00460 jet_et = jet.et();
00461
00462 jet_pt = jet.pt();
00463 jet_eta = jet.eta();
00464 jet_phi = jet.phi()*(Geom::pi()/180.);
00465 }
00466 else if (jet.et() > jet2_et) {
00467
00468 jet2_et = jet.et();
00469
00470
00471 jet2_eta = jet.eta();
00472 jet2_phi = jet.phi();
00473 }
00474
00475
00476 }
00478
00479
00480
00482
00484
00485 bool fill_e1 = false;
00486 bool fill_e2 = false;
00487 bool fill_m1 = false;
00488 bool fill_m2 = false;
00489 bool fill_met = false;
00490
00491
00492 if (ee_invMass>0.0) {
00493 h_ee_invMass ->Fill(ee_invMass);
00494 fill_e1 = true;
00495 fill_e2 = true;
00496 }
00497
00498
00499 if (mm_invMass > 0.0) {
00500 h_mumu_invMass->Fill(mm_invMass);
00501 fill_m1 = true;
00502 fill_m2 = true; h_jet2_et ->Fill(jet2_et);
00503 }
00504
00505
00506 if (electron_et>0.0&&missing_et>20.0) {
00507 float dphiW = fabs(met_phi-electron_phi);
00508 float W_mt_e = sqrt(2*missing_et*electron_et*(1-cos(dphiW)));
00509 h_e_invWMass ->Fill(W_mt_e);
00510 fill_e1 = true;
00511 fill_met = true;
00512 }
00513
00514
00515 if (muon_pt>0.0&&missing_et>20.0) {
00516 float dphiW = fabs(met_phi-muon_phi);
00517 float W_mt_m = sqrt(2*missing_et*muon_pt*(1-cos(dphiW)));
00518 h_m_invWMass ->Fill(W_mt_m);
00519 fill_m1 = true;
00520 fill_met = true;
00521 }
00522
00523
00524
00525 if (jet_et>-10.0) {
00526 h_jet_et ->Fill(jet_et);
00527 h_jet_count->Fill(jet_count);
00528 }
00529
00530 if (jet_pt>0.) {
00531 h_jet_pt ->Fill(jet_pt);
00532 }
00533
00534 if (jet_eta>-50.) {
00535 h_jet_eta ->Fill(jet_eta);
00536 }
00537
00538 if (jet_phi>-10.) {
00539 h_jet_phi ->Fill(jet_phi);
00540 }
00541
00542 if (jet2_et>-10.0) {
00543 h_jet2_et ->Fill(jet2_et);
00544 }
00545
00546
00547
00548
00549
00550 if (jet2_eta>-50.) {
00551 h_jet2_eta ->Fill(jet2_eta);
00552 }
00553
00554 if (jet2_phi>-10.) {
00555 h_jet2_phi ->Fill(jet2_phi);
00556 }
00557
00558
00559
00560 if (fill_e1 || fill_m1) {
00561 h_vertex_number->Fill(vertex_number);
00562 h_vertex_chi2->Fill(vertex_chi2);
00563 h_vertex_d0 ->Fill(vertex_d0);
00564 h_vertex_numTrks->Fill(vertex_numTrks);
00565 h_vertex_sumTrks->Fill(vertex_sumTrks);
00566 }
00567
00568 if (fill_e1) {
00569 h_e1_et ->Fill(electron_et);
00570 h_e1_eta ->Fill(electron_eta);
00571 h_e1_phi ->Fill(electron_phi);
00572 }
00573 if (fill_e2) {
00574 h_e2_et ->Fill(electron2_et);
00575 h_e2_eta ->Fill(electron2_eta);
00576 h_e2_phi ->Fill(electron2_phi);
00577 }
00578 if (fill_m1) {
00579 h_m1_pt ->Fill(muon_pt);
00580 h_m1_eta ->Fill(muon_eta);
00581 h_m1_phi ->Fill(muon_phi);
00582 }
00583 if (fill_m2) {
00584 h_m2_pt ->Fill(muon2_pt);
00585 h_m2_eta ->Fill(muon2_eta);
00586 h_m2_phi ->Fill(muon2_phi);
00587 }
00588 if (fill_met) {
00589 h_met ->Fill(missing_et);
00590 h_met_phi ->Fill(met_phi);
00591 }
00593 }
00594
00595
00596 void EwkDQM::endJob(void) {}
00597
00598
00599
00600 double EwkDQM::calcDeltaPhi(double phi1, double phi2) {
00601
00602 double deltaPhi = phi1 - phi2;
00603
00604 if (deltaPhi < 0) deltaPhi = -deltaPhi;
00605
00606 if (deltaPhi > 3.1415926) {
00607 deltaPhi = 2 * 3.1415926 - deltaPhi;
00608 }
00609
00610 return deltaPhi;
00611 }